-
Notifications
You must be signed in to change notification settings - Fork 0
/
gravmodelgen.py
148 lines (131 loc) · 6.01 KB
/
gravmodelgen.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
#21 NOV 2019
#gravmodelgen.py
#This is a simple script to parse a coefficient file into a c++ header file
#To run, use command "python gravmodelgen.py", add -h flag for help
#it uses fully normalized coefficient in the zero tide tide system
# The coefficients are stored by:
# index=((((NMAX<<1) - m + 1)*m)>>1) + n
import math
def main(infilename, headerfilename, maxdegree,earth_radius,earth_gravity_constant):
"""parse infilename into headerfilename c++ header file
The coefficents are fully-normalized.
indexing is by ((2*maxdegree-m+1)*m)/2+n
Args:
infilename(filename): the .COF file that contains the
spherical harmonic coefficents, download this from http://www2.csr.utexas.edu/grace/gravity/
headerfilename(string ending in .hpp): the c++ header file.
max_degree(positive int >1): max degree and order of the model
earth_radius(positive float): the radius of the earth for the model (m), ex 0.6378136300E+07
earth_gravity_constant(positive float): the GM of the earth for the model (m^3/s^2), ex 0.3986004415E+15
Stored in a Coeff, index=((2*maxdegree-m+1)*m)/2+n
template<int NMAX>//NMAX maximum degree and order
struct Coeff{
real_t earth_radius;
real_t earth_gravity_constant;
real_t _Cnm[Csize(NMAX, NMAX)];
real_t _Snm[Ssize(NMAX, NMAX)];"""
data= parseescof(infilename,maxdegree)
cofs= data
c_cofs=[0]*(((maxdegree+1)*(maxdegree+2))//2)
s_cofs=[0]*(((maxdegree+1)*(maxdegree+2))//2)
outstr = '#pragma once\n#include "geograv.hpp"\n//file generated from %s, by gravmodelgen.py\n'%(infilename)
for cof in cofs:
n= cof[0]
m= cof[1]
c= cof[2]
s= cof[3]
c_cofs[((2*maxdegree-m+1)*m)//2+n]= c
s_cofs[((2*maxdegree-m+1)*m)//2+n]= s
c_cofs[0]= 0;#manually set the earth mass term to zero
j2=c_cofs[2]
c_cofs[2]= 0;#manually set the J2 term to zero
outstr = outstr+header_file_model_code(headerfilename[:-4],maxdegree,earth_radius,earth_gravity_constant,j2,c_cofs,s_cofs)
with open(headerfilename,'w') as f:
f.write(outstr)
def header_file_model_code(modelname,max_degree, earth_radius,earth_gravity_constant,j2,c_cofs,s_cofs):
"""return the code defining the spherical coefficents and model
Stored in a Coeff, index=((2*maxdegree-m+1)*m)/2+n
template<int NMAX>//NMAX maximum degree and order
struct Coeff{
real_t earth_radius;
real_t earth_gravity_constant;
real_t _Cnm[Csize(NMAX, NMAX)];
real_t _Snm[Ssize(NMAX, NMAX)];
Args:
modelname(str, a valid C++ name): name of the model
max_degree(positive int >1): max degree and order of the model
earth_radius(positive float): the radius of the earth for the model (m), ex 0.6378136300E+07
earth_gravity_constant(positive float): the GM of the earth for the model (m^3/s^2), ex 0.3986004415E+15
j2(positive float): the j2 or c_cofs[2] term ex, -0.000484165371736
c_cofs,s_cofs(list of floats): coefficents of the model"""
head="constexpr geograv::Coeff<%d> %s={%sL,%sL,%sL,\n"%(max_degree,modelname,repr(earth_radius),repr(earth_gravity_constant),repr(j2))
modeltail= '};\n\n'
cs='{'
ss='{'
for i in range(len(c_cofs)):
cs+= repr(c_cofs[i])+'L,'
ss+= repr(s_cofs[i])+'L,'
cs= cs[:-1]+'}'
ss= ss[:-1]+'}'
sqrttable='{'
for i in range(max(2*max_degree+ 5, 15) + 1):
sqrttable+= repr(math.sqrt(i))+'L,'
sqrttable= sqrttable[:-1]+'}'
return head+cs+',\n'+ss+',\n'+sqrttable+modeltail
def parseescof(infilename, maxdegree):
"""return a list of lists from the infilename cof data file
n,m,C,S,
...
Args:
infilename(string ending in .COF): the .COF file that contains the
spherical harmonic coefficents, download this from https://www.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml
maxdegree(positive integer): maximum degree"""
with open(infilename,'r') as f:
filestr= f.read()
filelines= filestr.split('\n')
l=[]
C=0
S=0
n=0
m=0
i=0
while(n<=maxdegree and m<=maxdegree and i<len(filelines)):
l.append([n,m,C,S])
rowa= filelines[i].split()
try:
n=int(rowa[0])
m=int(rowa[1])
C=float(rowa[2].replace('D','E'))
S=float(rowa[3].replace('D','E'))
except:
try:
n=int(rowa[1])
m=int(rowa[2])
C=float(rowa[3].replace('D','E'))
S=float(rowa[4].replace('D','E'))
except:
pass
i+=1
return l
if __name__ == '__main__':
import argparse
parser = argparse.ArgumentParser(
formatter_class=argparse.RawDescriptionHelpFormatter,
description="""
#Nathan Zimmerberg
#21 NOV 2019
#gravmodelgen.py
#This is a simple script to parse a coefficient file into a c++ header file
#To run, use command "python gravmodelgen.py", add -h flag for help
#it uses fully normalized coefficient in the zero tide tide system
# The coefficients are stored by:
# index=((((NMAX<<1) - m + 1)*m)>>1) + n
""")
parser.add_argument('-f',type=str,default='GGM05S.GEO',help="""the file that contains the
sherical harmonic coefficents""")
parser.add_argument('-o',type=str,default='GGM05S.hpp',help='the c++ header filename to write the coefficents, and the name of the model')
parser.add_argument('-n',type=int,default=20,help='maximum number of degrees to use')
parser.add_argument('-r',type=float,default=0.6378136300E+07,help='radius of earth for model(m)')
parser.add_argument('-u',type=float,default=0.3986004415E+15,help='earth_gravity_constant for model(m^3/s^2)')
arg=parser.parse_args()
main(arg.f,arg.o,arg.n,arg.r,arg.u)