Commit e365d374 authored by thomas's avatar thomas
Browse files

dg: GBR testcase update

git-svn-id: https://geuz.org/svn/gmsh/trunk@13704 59c8cabf-10d6-4a28-b43f-3936b85956c4
parent ee4bfff0
......@@ -8,6 +8,9 @@ from math import sqrt
#----------------
mesh = "400K_ext"
#mesh = "100K"
#SET RADIUS: 6.37101e6 for 'sma_zoomed', else 6.371e6
R = 6.371e6
#================
#================
......@@ -151,15 +154,26 @@ evalPointsLatLon = [
[-21.8200, 151.1400, 'Bell', 8] ]
[-18.241667, 147.366667, 'Myrmidon Reef', 9], #NB: WAS -18.2452, 147.4100
tidalPointsAbbotLatLon = [
[-19.8187007904, 148.013000488, 'Deepwater West', 0],
[-19.9050006866, 148.005004883, 'Coastal West', 1],
[-19.8815994263, 148.074996948, 'Coastal Footprint', 2],
[-19.8484992981, 148.04699707, 'Inshore Footprint', 3],
[-19.8791999817, 147.970993042, 'Inshore West', 4],
[-19.831199646, 148.115005493, 'Deepwater East', 5],
[-19.9340991974, 148.130004883, 'Coastal East', 6],
[-19.8999004364, 148.158996582, 'Inshore East', 7] ]
#tidalPointsAbbotLatLon = [
#[-19.8187007904, 148.013000488, 'Deepwater West', 0],
#[-19.9050006866, 148.005004883, 'Coastal West', 1],
#[-19.8815994263, 148.074996948, 'Coastal Footprint', 2],
#[-19.8484992981, 148.04699707, 'Inshore Footprint', 3],
#[-19.8791999817, 147.970993042, 'Inshore West', 4],
#[-19.831199646, 148.115005493, 'Deepwater East', 5],
#[-19.9340991974, 148.130004883, 'Coastal East', 6],
#[-19.8999004364, 148.158996582, 'Inshore East', 7] ]
evalPointsGBROOSLatLon = [
[-14.33948333, 145.3436667, 'Lizard Island Slope', 0],
[-14.69953333, 145.6436, 'Lizard Island Shelf', 1],
[-18.31035, 147.1632667, 'Palm Passage Shelf', 2],
[-18.21963333, 147.3436833, 'Myrmidon Slope', 3],
[-21.03088333, 152.8793, 'Elusive Reef', 4],
[-22.40815, 151.9933667, 'Capricorn Channel', 5],
[-23.38343333,151.9868167, 'Heron Island North', 6],
[-23.5133, 151.9549167, 'Heron Island South', 7],
[-23.48338333, 152.1726333, 'One Tree East', 8] ]
tidalPointsSeafarerLatLon = [
[-18.466667, 146.866667, 'Rib Reef', 0],
......
......@@ -22,7 +22,7 @@ print ''
if (Msg.GetCommRank()==0):
print ''
print(colored('Compiling Libraries ....', "red"))
print(colored('Compiling Libraries ...', "red"))
functionC.buildLibraryFromFile ("Libraries/GBR.cc", "Libraries/lib_gbr.so");
# functionC.buildLibraryFromFile ("Libraries/REEF.cc", "Libraries/lib_reef.so");
......@@ -39,7 +39,7 @@ if ( len(sys.argv)>1 ):
simOutputDir = '/scratch/cthomas/Hydro_Realistic/gbr500K_sma_1svIn_bDf10_allTides_windOn_28112007_5wks/'
#simOutputDir = 'Output10K_4wks_windsOn_allTides_16112000'
LP_OutputDir = '/scratch/cthomas/temp'#simOutputDirName +'/LPT_Output/SettProg_Mort007'
LP_OutputDir = '/scratch/cthomas/temp'#simOutputDir +'/LPT_Output/SettProg_Mort007'
#LP_OutputDir = 'temp'
......@@ -102,8 +102,8 @@ if(not os.path.exists('./LPT_Functions')):
simSolutions = [simExportGap*solutionNumber for solutionNumber in range(0,simTotIter/simExportGap+1)]
print(colored('Loading particle-tracking parameters ...',"red"))
LP_dt = 120. #(double) in seconds
LP_TotTime = 3*7*24.0 #(double) in hours
LP_dt = 90. #(double) in seconds
LP_TotTime = 4*7*24.0 #(double) in hours
LP_Export = int( (12.0*3600.0)/LP_dt ) #(int) export particle positions every ... steps
LP_Export_ConnectivityMatrix = int( (24.0*3600.0)/LP_dt ) #(int) export connectivity matrix every ... steps
LP_TotIter = int(LP_TotTime*3600.0/LP_dt)
......@@ -139,7 +139,7 @@ settleParams = fullMatrixDouble(8,8)
#[0]: 0: settling off, 1: const. settling, 2: triangle settling, 3: progressive competence acquisition/loss
settleParams.set(0,0, 3.0)
#[1]: t_preComp (seconds) == for type-3: t_timeLag
# A. Gemmifera: 3.47
# G. Retiformis: 0.0, A. Gemmifera: 3.47
settleParams.set(0,1, 3.47*24.0*3600.0)
#[2]: t_postComp (for ALL settling types)
settleParams.set(0,2, 100.0*24.0*3600.0)
......@@ -148,10 +148,10 @@ settleParams.set(0,3, -1.0)
#[4]: (max if [0]=2) settling rate (PER dt) - **! must convert to .dt(-1) !**
settleParams.set(0,4, 1.0)
#[5]: (for type-3 only) rate of acquisition of competence (d-1)
# A. Gemmifera: 0.39
# G. Retiformis: 0.58, A. Gemmifera: 0.39
settleParams.set(0,5, 0.39)
#[6]: (for type-3 only) rate of loss of competence (d-1)
# A. Gemmifera: 0.145
# G. Retiformis: 0.096, A. Gemmifera: 0.145
settleParams.set(0,6, 0.145)
#[7]: settle over shallow, deep or both reef types
settleParams.set(0,7, settleOverWhichReefs)
......@@ -187,7 +187,7 @@ XYZ_DC.exportMsh("XYZ_DC",0,0)
bathDC = dgDofContainer(groups, 1);
if(os.path.exists("./Bath/"+simMesh+"_bath_smooth/"+simMesh+"_bath_smooth.idx")):
print(colored('Importing bathymetry.',"red"))
print(colored('Importing bathymetry ...',"red"))
bathDC.readMsh("./Bath/"+simMesh+"_bath_smooth/"+simMesh+"_bath_smooth.idx")
else:
print(colored('**** ERROR: Bathymetry file does not exist! ****',"red"))
......@@ -287,7 +287,7 @@ nbOfReefs = nbOfShallowReefs + nbOfDeepReefs
print "Maximum possible number of reefs: ",nbOfReefs, "(shallow:", nbOfShallowReefs, ", deep:", nbOfDeepReefs, ")"
# Record simulation info to file
LP_SimInfo = [str(nbOfDeepReefs)+'\n'+str(LP_dt)+'\n'+str((LP_StartTime-simStartTime)/3600)+'\n' +str(LP_TotTime)+'\n'+str(LP_Export)+'\n'+ +str(simItersPerLPIter)+'\n'+str(mortParams(0,0))+'\n' +str(mortParams(0,1)/(24.0*3600.0))+'\n'+str(settleParams(0,0))+'\n','Pre-comptetency period:'+str(settleParams(0,1)/(24.0*3600.0))+' (days)\n', 'Post-competent at:'+str(settleParams(0,2)/(24.0*3600.0))+' days)\n', '(max) settling rate:'+str(settleParams(0,4))+' per dt\nConcentration:'+str(concentration)+'/kmsq\n', 'Seeding over: '+str(seedOverWhichReefs)+'\n', 'Settling over: '+str(settleOverWhichReefs)+'\n', 'Highest shallow reef Id: '+str(nbOfShallowReefs)+'\n']
LP_SimInfo = [str(nbOfDeepReefs),'\n'+str(LP_dt),'\n',str((LP_StartTime-simStartTime)/3600),'\n', str(LP_TotTime),'\n',str(LP_Export),'\n', str(simItersPerLPIter),'\n',str(mortParams(0,0)),'\n' ,str(mortParams(0,1)/(24.0*3600.0)),'\n',str(settleParams(0,0)),'\n','Pre-comptetency period:'+str(settleParams(0,1)/(24.0*3600.0))+' (days)\n', 'Post-competent at:'+str(settleParams(0,2)/(24.0*3600.0))+' days)\n', '(max) settling rate:'+str(settleParams(0,4))+' per dt\nConcentration:'+str(concentration)+'/kmsq\n', 'Seeding over: '+str(seedOverWhichReefs)+'\n', 'Settling over: '+str(settleOverWhichReefs)+'\n', 'Highest shallow reef Id: '+str(nbOfShallowReefs)+'\n']
LP_SimInfoFile = open(LP_OutputDir+'/LP_runInfo.txt','w')
LP_SimInfoFile.writelines(LP_SimInfo)
LP_SimInfoFile.close()
......@@ -339,7 +339,7 @@ else:
nbOfShallowReefsFConst = functionConstant(nbOfShallowReefs)
reefsDualDC = dgDofContainer(groups, 1)
reefsDualFC = functionC(glib, "overwrite2", 1, [reefsShallowDC.getFunction(), reefsDeepDC.getFunction(), nbOfShallowReefsFConst])
reefsDualFC = functionC(glib, "overwrite2", 1, [bathDC.getFunction(), reefsShallowDC.getFunction(), reefsDeepDC.getFunction(), nbOfShallowReefsFConst])
reefsDualDC.interpolate(reefsDualFC)
reefsDualDC.exportMsh("LPT_Functions/REEFS_RAW_LPT_DUAL_"+simMesh,0,0)
......@@ -356,7 +356,7 @@ else:
# ***SEED AT A POINT*** (NB: Look at method used to print results if changing seeding method)
#cartCoords = [0]*3
#geoCoordTo3DCartesian(-19.9400,149.1100,cartCoords)
#geoCoordTo3DCartesian(-19.9400,149.1100,cartCoords,6.371e6)
#projectPoint(cartCoords)
#particles = particleArray(groups, 0, cartCoords[0], cartCoords[1], 0, 0)
#for location in range(0,10):
......
......@@ -258,11 +258,11 @@ void merge3(dataCacheMap *,fullMatrix<double> &sol,fullMatrix<double> &eta, full
}
//Function to merge shallow (hard) and deep (soft) DCs
void overwrite2(dataCacheMap *,fullMatrix<double> &sol, fullMatrix<double> &shallowDC, fullMatrix<double> &deepDC, fullMatrix<double> &nbOfShallowReefs){
void overwrite2(dataCacheMap *,fullMatrix<double> &sol, fullMatrix<double> &bath, fullMatrix<double> &shallowDC, fullMatrix<double> &deepDC, fullMatrix<double> &nbOfShallowReefs){
for (size_t i = 0; i < sol.size1(); i++) {
if (shallowDC(i,0) > 0.1) sol.set(i, 0, shallowDC(i,0));
else {
if (deepDC(i,0) > 0.1) sol.set(i, 0, deepDC(i,0)+nbOfShallowReefs(0,0));
if ((deepDC(i,0) > 0.1) && (bath(i,0) > 6.0)) sol.set(i, 0, deepDC(i,0)+nbOfShallowReefs(0,0));
else sol.set(i, 0, 0.0);
}
}
......@@ -280,7 +280,7 @@ void wind (dataCacheMap *, fullMatrix<double> &sol, fullMatrix<double> &xyz) {
void coriolis (dataCacheMap *, fullMatrix<double> &sol, fullMatrix<double> &xyz) {
double theta = -17.147569/180*M_PI;
for (size_t i = 0; i< sol.size1(); i++) {
//sol.set(i,0,2*7.292e-5*xyz(i,2)/6.37101e6);
//sol.set(i,0,2*7.292e-5*xyz(i,2)/6.371e6);
sol.set(i,0,2*7.292e-5*sin(theta));
}
}
......
#"""
#Simulations on the GBR with multirate explicit Runge-Kutta time-stepping.
#All the parameters can be found and modified in the file "InpuGBR.py"
# Simulations on the GBR with multirate explicit Runge-Kutta time-stepping.
# All the parameters can be found and modified in the file "InpuGBR.py"
#"""
#"""
#*** How to launch the testcase ***
#Simulation on 1 processor : rundgpy MultirateGBR.py
#Simulation on n processors : mpirun -np n rundgpy MultirateGBR.py ( !!! You need MPI !!! )
#**********************************
#"""
#"""
# TO START SIMULATION FROM A PRE-EXISTING SOLUTION:
# 1. Set the following in Input.py:
# A) HOW TO LAUNCH A SIMULATION
# Simulation on 1 processor : rundgpy MultirateGBR.py
# Simulation on n processors : mpirun -np n rundgpy MultirateGBR.py ( !!! You need MPI !!! )
#
#
# B) TO START SIMULATION FROM A PRE-EXISTING SOLUTION:
# 1. Set the following in Input.py:
# -> set flag startFromExistingSolution = 1
# -> set solutionFileName
# -> set stepNbOfInitialSolution
# -> set old_dt
# -> then, Ti = old Ti + stepNbOfInitialSolution*old_dt
# 2. Set the following in MultirateGBR.py:
# 2. Set the following in MultirateGBR.py:
# -> Nothing
#
# ** Must point to a 'merged' .idx file (use generateIdxFiles.py, and make sure you adjust nb_parts)
# ** Make sure idx file is complete (ie open it and check it looks right) - especially if previous job crashed out!
# Notes: Must point to a 'merged' .idx file (use generateIdxFiles.py, and make sure you adjust nb_parts)
# Make sure idx file is complete (ie open it and check it looks right) - especially if previous job crashed out!
#"""
# Import for python
......@@ -282,6 +277,7 @@ elif(gbrResidual == 1):
claw.addBoundaryCondition('Open Sea South', claw.newOutsideValueBoundary("Surface", CurrentsAndTideSouth))
claw.addBoundaryCondition('Coast',claw.newBoundaryWall())
claw.addBoundaryCondition('Island',claw.newBoundaryWall())
claw.addBoundaryCondition('Islands',claw.newBoundaryWall())
# **********************************************
......@@ -337,7 +333,7 @@ if (Msg.GetCommRank() == 0):
print ' - Final Time (Tf):', printDate(Tf)
print ' - Simulation Length (days):', (Tf-Ti)/3600.0/24.0
print ' - Reference Time Step:', printTime(dt), ', which is', str(dt), 's'
print ' - Number of Iterartions:', nbSteps
print ' - Number of Iterations:', nbSteps
print ' - Theoretical Speedup:', rk.getSpeedUp()
print ' - Norm of Solution at Ti:', norm
print ' - Mass of Solution at Ti:', totalMassRef(0, 0)
......@@ -355,7 +351,7 @@ if (Msg.GetCommRank() == 0):
# **********************************************
if (Msg.GetCommRank() == 0):
print 'Loading evaluator points.'
print 'Loading evaluator points ...'
# ********* Current Evaluator Points **********
#if (Msg.GetCommRank() == 0):
# print 'Current evaluator points.'
......@@ -364,7 +360,7 @@ for i in range(0,len(evalPointsLatLon)):
cartCoords = [0]*3
evalPointsCart[i][3] = evalPointsLatLon[i][2] #Site name
evalPointsCart[i][4] = evalPointsLatLon[i][3] #Site Id#
geoCoordTo3DCartesian(evalPointsLatLon[i][0],evalPointsLatLon[i][1],cartCoords)
geoCoordTo3DCartesian(evalPointsLatLon[i][0],evalPointsLatLon[i][1],cartCoords, R)
projectPoint(cartCoords)
for j in range(0,2):
evalPointsCart[i][j] = cartCoords[j]
......@@ -376,27 +372,29 @@ if (Msg.GetCommRank() == 0):
outEval = open(outputDir+'/Evaluator/evalCurrent_%02d'%(siteId)+'.dat',"w")
outEval.close()
# ********* Elevation Evaluator Points **********
# *** Abbot points ***
# ********* Elevation and current Evaluator Points **********
# *** GBROOS points ***
#if (Msg.GetCommRank() == 0):
# print ' "Abbot" Elevation evaluator points:'
tidalPointsAbbotCart = [ [0]*5 for i in range(len(tidalPointsAbbotLatLon)) ]
for i in range(0,len(tidalPointsAbbotLatLon)):
# print ' "GBROOS" Elevation & current evaluator points:'
evalPointsGBROOSCart = [ [0]*5 for i in range(len(evalPointsGBROOSLatLon)) ]
for i in range(0,len(evalPointsGBROOSLatLon)):
cartCoords = [0]*3
tidalPointsAbbotCart[i][3] = tidalPointsAbbotLatLon[i][2] #Site name
tidalPointsAbbotCart[i][4] = tidalPointsAbbotLatLon[i][3] #Site Id#
geoCoordTo3DCartesian(tidalPointsAbbotLatLon[i][0],tidalPointsAbbotLatLon[i][1],cartCoords)
evalPointsGBROOSCart[i][3] = evalPointsGBROOSLatLon[i][2] #Site name
evalPointsGBROOSCart[i][4] = evalPointsGBROOSLatLon[i][3] #Site Id#
geoCoordTo3DCartesian(evalPointsGBROOSLatLon[i][0],evalPointsGBROOSLatLon[i][1],cartCoords, R)
projectPoint(cartCoords)
for j in range(0,2):
tidalPointsAbbotCart[i][j] = cartCoords[j]
evalPointsGBROOSCart[i][j] = cartCoords[j]
# if (Msg.GetCommRank() == 0):
# print str(tidalPointsAbbotCart[i][3]),': ', tidalPointsAbbotLatLon[i][0], '(lat) |', tidalPointsAbbotLatLon[i][1], '(lon)'
# print 'x(flat): ', tidalPointsAbbotCart[i][0]/1000, '(km) | y(flat): ', tidalPointsAbbotCart[i][1]/1000,'(km) | z(flat): ', tidalPointsAbbotCart[i][2]/1000, '(km)'
# print str(evalPointsGBROOSCart[i][3]),': ', evalPointsGBROOSLatLon[i][0], '(lat) |', evalPointsGBROOSLatLon[i][1], '(lon)'
# print 'x(flat): ', evalPointsGBROOSCart[i][0]/1000, '(km) | y(flat): ', evalPointsGBROOSCart[i][1]/1000,'(km) | z(flat): ', evalPointsGBROOSCart[i][2]/1000, '(km)'
if (Msg.GetCommRank() == 0):
for siteId in range(0,len(tidalPointsAbbotCart)):
outElevEval = open(outputDir+'/Evaluator/evalElevAbbot_%02d'%(siteId)+'.dat',"w")
for siteId in range(0,len(evalPointsGBROOSCart)):
outElevEval = open(outputDir+'/Evaluator/evalElevGBROOS_%02d'%(siteId)+'.dat',"w")
outCurrentEval = open(outputDir+'/Evaluator/evalCurrentGBROOS_%02d'%(siteId)+'.dat',"w")
outElevEval.close()
# ********* Elevation Evaluator Points **********
# *** Seafarer points ***
#if (Msg.GetCommRank() == 0):
# print ' "Seafarer" Elevation evaluator points:'
......@@ -405,7 +403,7 @@ for i in range(0,len(tidalPointsSeafarerLatLon)):
cartCoords = [0]*3
tidalPointsSeafarerCart[i][3] = tidalPointsSeafarerLatLon[i][2] #Site name
tidalPointsSeafarerCart[i][4] = tidalPointsSeafarerLatLon[i][3] #Site Id#
geoCoordTo3DCartesian(tidalPointsSeafarerLatLon[i][0],tidalPointsSeafarerLatLon[i][1],cartCoords)
geoCoordTo3DCartesian(tidalPointsSeafarerLatLon[i][0],tidalPointsSeafarerLatLon[i][1],cartCoords, R)
projectPoint(cartCoords)
for j in range(0,2):
tidalPointsSeafarerCart[i][j] = cartCoords[j]
......@@ -425,7 +423,7 @@ for i in range(0,len(reefPointsLatLon)):
cartCoords = [0]*3
reefPointsCart[i][3] = reefPointsLatLon[i][2] #Site name
reefPointsCart[i][4] = reefPointsLatLon[i][3] #Site Id#
geoCoordTo3DCartesian(reefPointsLatLon[i][0],reefPointsLatLon[i][1],cartCoords)
geoCoordTo3DCartesian(reefPointsLatLon[i][0],reefPointsLatLon[i][1],cartCoords, R)
projectPoint(cartCoords)
for j in range(0,2):
reefPointsCart[i][j] = cartCoords[j]
......@@ -444,12 +442,18 @@ ElevationEval = dgFunctionEvaluator(groups, solution.getFunction())
res = fullMatrixDouble()
#Bathymetry Evaluator
outBathEval = open(outputDir+'/Evaluator/evalBathymetry.dat',"w")
outBathEval = open(outputDir+'/Evaluator/evalBathymetry_StdLocations.dat',"w")
outBathGBROOSEval = open(outputDir+'/Evaluator/evalBathymetry_GBROOSLocations.dat',"w")
for point in range(0,len(evalPointsLatLon)):
#Only write to file if we are on correct cpu:
BathEval.compute(evalPointsCart[point][0],evalPointsCart[point][1],evalPointsCart[point][2], res)
outBathEval.write(str(res(0,0))+'\n')
outBathEval.close()
for point in range(0,len(evalPointsGBROOSLatLon)):
#Only write to file if we are on correct cpu:
BathEval.compute(evalPointsGBROOSCart[point][0],evalPointsGBROOSCart[point][1],evalPointsGBROOSCart[point][2], res)
outBathGBROOSEval.write(str(res(0,0))+'\n')
outBathGBROOSEval.close()
# **********************************************
exporterSol.exportIdx(0, t)
......@@ -510,17 +514,22 @@ for i in range(firstStep, lastStep):
outEval = open(outputDir+'/Evaluator/evalCurrent_%02d'%(point)+'.dat',"a")
outEval.write(str(t/3600-Ti/3600)+'\t'+str(res.get(0,0))+'\t'+str(res.get(0,1))+'\n')
outEval.close()
#*** Elevation evaluator - Abbot: ***
for point in range(0,len(tidalPointsAbbotLatLon)):
#*** Elevation evaluator - GBROOS: ***
for point in range(0,len(evalPointsGBROOSLatLon)):
#Only write to file if we are on correct cpu:
evalPointSP = SPoint3(tidalPointsAbbotCart[point][0],tidalPointsAbbotCart[point][1],tidalPointsAbbotCart[point][2])
evalPointSP = SPoint3(evalPointsGBROOSCart[point][0],evalPointsGBROOSCart[point][1],evalPointsGBROOSCart[point][2])
evalPointElem = model.getMeshElementByCoord(evalPointSP)
whichPartition = evalPointElem.getPartition()
ElevationEval.compute(tidalPointsAbbotCart[point][0],tidalPointsAbbotCart[point][1],tidalPointsAbbotCart[point][2], res)
ElevationEval.compute(evalPointsGBROOSCart[point][0],evalPointsGBROOSCart[point][1],evalPointsGBROOSCart[point][2], res)
res2 = fullMatrixDouble()
FunctionEval.compute(evalPointsGBROOSCart[point][0],evalPointsGBROOSCart[point][1],evalPointsGBROOSCart[point][2], res2)
if (Msg.GetCommRank()+1 == whichPartition):
outElevEval = open(outputDir+'/Evaluator/evalElevAbbot_%02d'%(point)+'.dat',"a")
outElevEval = open(outputDir+'/Evaluator/evalElevGBROOS_%02d'%(point)+'.dat',"a")
outElevEval.write(str(t/3600-Ti/3600)+'\t'+str(res.get(0,0))+'\n')
outElevEval.close()
outEval = open(outputDir+'/Evaluator/evalCurrentGBROOS_%02d'%(point)+'.dat',"a")
outEval.write(str(t/3600-Ti/3600)+'\t'+str(res2.get(0,0))+'\t'+str(res2.get(0,1))+'\n')
outEval.close()
#*** Elevation evaluator - Seafarer: ***
for point in range(0,len(tidalPointsSeafarerLatLon)):
#Only write to file if we are on correct cpu:
......
......@@ -3,10 +3,11 @@ from dgpy import *
from gmshpy import *
from math import *
# !!! WARNING: OLD MESHES HAVE R=6.37101e6, NEW MESHES HAVE R=6.371e6 !!! #
# Projection Parameters
phi = 146.5892/180*pi
theta = -17.147569/180*pi
Radius = 6.37101e6;
Ox = cos(theta)*cos(phi)
Oy = cos(theta)*sin(phi)
Oz = sin(theta)
......@@ -56,12 +57,12 @@ def projectMesh (model):
project(f)
#takes geo lat/lon and outputs x,y,z(sphere)
def geoCoordTo3DCartesian(geoLat, geoLon, pointCart):
def geoCoordTo3DCartesian(geoLat, geoLon, pointCart, R):
#calculate radius at given latitude
a = 6378137.0 #Earth SemiMajor
b = 6356752.3 #Earth SemiMinor
#R = sqrt( ((a*a*(cos(geoLat*pi/180.0)))**2 + (b*b*(sin(geoLat*pi/180.0)))**2) / ((a*(cos(geoLat*pi/180.0)))**2 + (b*(sin(geoLat*pi/180.0)))**2) )
R = 6.371e6;
#R = 6.371e6;
#geo coords --> S.P.C.
pointTheta = 90 - geoLat
......@@ -74,8 +75,8 @@ def geoCoordTo3DCartesian(geoLat, geoLon, pointCart):
#Takes x,y,z(sphere) as input and overwrites with alpha,beta,0(stereo)
def threeDCartesianToStereo(point):
R = 6.371e6
def threeDCartesianToStereo(point, R):
#R = 6.371e6
alpha = (2.0 * R * point[0]) / (R + point[2])
beta = (2.0 * R * point[1]) / (R + point[2])
point[0] = alpha
......@@ -174,9 +175,9 @@ def stereoVectorToLonLat(uvLonLat, alpha, beta, uv_stereo):
uvLonLat.set( 0, 1, sol_v );
uvLonLat.set( 0, 2, 0.0 );
def stereoToLonLatPoint(lonLat, xyz):
def stereoToLonLatPoint(lonLat, xyz, R):
pi = acos(-1.0)
R = 6.371e6
#R = 6.371e6
for i in range(0,lonLat.size1()):
xi = xyz.get(i,0)
eta = xyz.get(i,1)
......@@ -189,10 +190,10 @@ def stereoToLonLatPoint(lonLat, xyz):
def lonLat_fn(lonLat, xyz):
for i in range(0,lonLat.size1()):
x = ePhiX * xyz(i,0) + eThetaX * xyz(i,1) + Ox*Radius;
y = ePhiY * xyz(i,0) + eThetaY * xyz(i,1) + Oy*Radius;
z = ePhiZ * xyz(i,0) + eThetaZ * xyz(i,1) + Oz*Radius;
r = sqrt(x*x+y*y+z*z);
x = ePhiX * xyz(i,0) + eThetaX * xyz(i,1) + Ox*r
y = ePhiY * xyz(i,0) + eThetaY * xyz(i,1) + Oy*r
z = ePhiZ * xyz(i,0) + eThetaZ * xyz(i,1) + Oz*r
lonLat.set(i, 0, atan2(y,x));
lonLat.set(i, 1, asin(z/r));
lonLat.set(i, 2, 0);
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment