##// END OF EJS Templates
MKFact No Hardcoded
rflores -
r1569:c96fc34a400b
parent child
Show More
@@ -1,668 +1,668
1 subroutine geobfield(tm,r,theta,phi,br,bt,bp,b)
1 subroutine geobfield(tm,r,theta,phi,br,bt,bp,b)
2 c*****evaluate bfield at geocentric coords r,theta,phi
2 c*****evaluate bfield at geocentric coords r,theta,phi
3 c and return components of field in r,theta and phi directions
3 c and return components of field in r,theta and phi directions
4 c and magnitude
4 c and magnitude
5 c*****theta and phi are in radians, r is in km, b is in gauss
5 c*****theta and phi are in radians, r is in km, b is in gauss
6 c ===============================================================
6 c ===============================================================
7 c
7 c
8 c essentially, this is an adaptation of igrfdemo. it interpolates
8 c essentially, this is an adaptation of igrfdemo. it interpolates
9 c smoothly between the coefficients from 1945 to 1985.
9 c smoothly between the coefficients from 1945 to 1985.
10 c
10 c
11 c subroutines:
11 c subroutines:
12 c getshc, interpshc, extrapshc, shval3
12 c getshc, interpshc, extrapshc, shval3
13 c
13 c
14 c igrfdemo is due to:
14 c igrfdemo is due to:
15 c a. zunde
15 c a. zunde
16 c usgs, ms 964, box 25046 federal center, denver, co 80225
16 c usgs, ms 964, box 25046 federal center, denver, co 80225
17 c
17 c
18 c ===============================================================
18 c ===============================================================
19 character*8 filmod(17)
19 character*8 filmod(17)
20 character*31 fqual
20 character*9 fqual
21 dimension gh1(220), gh2(220), gha(224), ext(3), dtemod(17)
21 dimension gh1(220), gh2(220), gha(224), ext(3), dtemod(17)
22 data ext /3*0./
22 data ext /3*0./
23 data filmod / 'dgrf45', 'dgrf50',
23 data filmod / 'dgrf45', 'dgrf50',
24 1 'dgrf55', 'dgrf60', 'dgrf65',
24 1 'dgrf55', 'dgrf60', 'dgrf65',
25 2 'dgrf70', 'dgrf75', 'dgrf80',
25 2 'dgrf70', 'dgrf75', 'dgrf80',
26 3 'dgrf85', 'dgrf90', 'dgrf95',
26 3 'dgrf85', 'dgrf90', 'dgrf95',
27 4 'dgrf00', 'dgrf05', 'dgrf10',
27 4 'dgrf00', 'dgrf05', 'dgrf10',
28 5 'dgrf15', 'igrf20', 'igrf20s'/
28 5 'dgrf15', 'igrf20', 'igrf20s'/
29 data fqual/"/usr/local/lib/faraday/bfmodel/"/
29 data fqual/"/bfmodel/"/
30 data dtemod / 1945., 1950., 1955., 1960.,
30 data dtemod / 1945., 1950., 1955., 1960.,
31 1 1965., 1970., 1975., 1980., 1985., 1990., 1995., 2000.,
31 1 1965., 1970., 1975., 1980., 1985., 1990., 1995., 2000.,
32 2 2005.,2010.,2015.,2020.,2025./
32 2 2005.,2010.,2015.,2020.,2025./
33 data iu/98/,ndt/17/
33 data iu/98/,ndt/17/
34 c data a2/40680925./, b2/40408588./
34 c data a2/40680925./, b2/40408588./
35 data a2/40680631.6/, b2/40408296.0/ ! updated 2010
35 data a2/40680631.6/, b2/40408296.0/ ! updated 2010
36 c
36 c
37 data rtd/57.29577951/
37 data rtd/57.29577951/
38 data tmp/0./,lp/0/
38 data tmp/0./,lp/0/
39 flat=90. - theta*rtd
39 flat=90. - theta*rtd
40 flon=rtd*phi
40 flon=rtd*phi
41 c*****if previous time is not equal to current time
41 c*****if previous time is not equal to current time
42 if(tm .ne. tmp)then
42 if(tm .ne. tmp)then
43 l=0
43 l=0
44 C for i=1,ndt
44 C for i=1,ndt
45 DO 1010 i=1,ndt
45 DO 1010 i=1,ndt
46 C exit for if(tm .lt. dtemod(i))
46 C exit for if(tm .lt. dtemod(i))
47 if(tm .lt. dtemod(i)) GO TO 1011
47 if(tm .lt. dtemod(i)) GO TO 1011
48 l=i
48 l=i
49 C end for
49 C end for
50 1010 CONTINUE
50 1010 CONTINUE
51 1011 CONTINUE
51 1011 CONTINUE
52 c write(*,-)tm,l
52 c write(*,-)tm,l
53 if(l .eq. 0)then
53 if(l .eq. 0)then
54 write(*,fmt='(" geobfield: time is before earliest model.")')
54 write(*,fmt='(" geobfield: time is before earliest model.")')
55 stop
55 stop
56 end if
56 end if
57 if(l .eq. ndt)then
57 if(l .eq. ndt)then
58 l=l-1
58 l=l-1
59 write(*,fmt='(" geobfield: warning - extrapolating beyond ",
59 write(*,fmt='(" geobfield: warning - extrapolating beyond ",
60 1 "last set of coefficients.")')
60 1 "last set of coefficients.")')
61 end if
61 end if
62 if(l .ne. lp)then
62 if(l .ne. lp)then
63 c*********if previous epoch not the same, read in new coefs
63 c*********if previous epoch not the same, read in new coefs
64 c write(*,fmt='(" read coefs"))')
64 c write(*,fmt='(" read coefs"))')
65 c write(*,*) fqual//filmod(l)
65 c write(*,*) fqual//filmod(l)
66 call getshc2 (iu, fqual//filmod(l), nmax1, erad, gh1, ier)
66 call getshc2 (iu, fqual//filmod(l), nmax1, erad, gh1, ier)
67 if(ier .ne. 0)then
67 if(ier .ne. 0)then
68 write(*,fmt='(" geobfield: read error=",i2," on ",a)')
68 write(*,fmt='(" geobfield: read error=",i2," on ",a)')
69 1 ier,filmod(l)
69 1 ier,filmod(l)
70 stop
70 stop
71 end if
71 end if
72
72
73 call getshc2 (iu, fqual//filmod(l+1), nmax2, erad, gh2, ier)
73 call getshc2 (iu, fqual//filmod(l+1), nmax2, erad, gh2, ier)
74 if(ier .ne. 0)then
74 if(ier .ne. 0)then
75 write(*,fmt='(" geobfield: read error=",i2," on ",a)')
75 write(*,fmt='(" geobfield: read error=",i2," on ",a)')
76 1 ier,filmod(l+1)
76 1 ier,filmod(l+1)
77 stop
77 stop
78 end if
78 end if
79 end if
79 end if
80 if (l .lt. ndt-1) then
80 if (l .lt. ndt-1) then
81 call interpshc (tm, dtemod(l), nmax1, gh1, dtemod(l+1),
81 call interpshc (tm, dtemod(l), nmax1, gh1, dtemod(l+1),
82 1 nmax2, gh2, nmax, gha)
82 1 nmax2, gh2, nmax, gha)
83 c write(*,fmt='(" interpolate"))')
83 c write(*,fmt='(" interpolate"))')
84 else
84 else
85 call extrapshc (tm, dtemod(l), nmax1, gh1, nmax2,
85 call extrapshc (tm, dtemod(l), nmax1, gh1, nmax2,
86 1 gh2, nmax, gha)
86 1 gh2, nmax, gha)
87 c write(*,fmt='(" extrapolate"))')
87 c write(*,fmt='(" extrapolate"))')
88 end if
88 end if
89 end if
89 end if
90 c tmp=tm
90 c tmp=tm
91 c lp=l
91 c lp=l
92 call shval3(2,flat,flon,r,erad,a2,b2,nmax,gha,0,ext,x,y,z)
92 call shval3(2,flat,flon,r,erad,a2,b2,nmax,gha,0,ext,x,y,z)
93 br=-z/100000.
93 br=-z/100000.
94 bp=y/100000.
94 bp=y/100000.
95 bt=-x/100000.
95 bt=-x/100000.
96 b=sqrt(br**2 + bp**2 + bt**2)
96 b=sqrt(br**2 + bp**2 + bt**2)
97
97
98 C write(*,*)flat,flon,r,erad,a2,b2,nmax,gha,ext
98 C write(*,*)flat,flon,r,erad,a2,b2,nmax,gha,ext
99 C write(*,*)tm,r,theta,phi,br,bt,bp,b
99 C write(*,*)tm,r,theta,phi,br,bt,bp,b
100 return
100 return
101 end
101 end
102 c
102 c
103 subroutine convrt(i, gdlat, gdalt, gclat, rkm)
103 subroutine convrt(i, gdlat, gdalt, gclat, rkm)
104
104
105 c convrt converts between geodetic and geocentric coordinates. the
105 c convrt converts between geodetic and geocentric coordinates. the
106 c reference geoid is that adopted by the iau in 1964. a=6378.16,
106 c reference geoid is that adopted by the iau in 1964. a=6378.16,
107 c b=6356.7746, f=1/298.25. the equations for conversion from
107 c b=6356.7746, f=1/298.25. the equations for conversion from
108 c geocentric to geodetic are from astron. j., vol 66, 1961, p. 15.
108 c geocentric to geodetic are from astron. j., vol 66, 1961, p. 15.
109
109
110 c i=1 geodetic to geocentric
110 c i=1 geodetic to geocentric
111 c i=2 geocentric to geodetic
111 c i=2 geocentric to geodetic
112 c gdlat geodetic latitude (degrees)
112 c gdlat geodetic latitude (degrees)
113 c gdalt altitude above geoid (km)
113 c gdalt altitude above geoid (km)
114 c gclat geocentric latitude (degrees)
114 c gclat geocentric latitude (degrees)
115 c rkm geocentric radial distance (km)
115 c rkm geocentric radial distance (km)
116 c
116 c
117 c data a/6378.16/, ab2/1.0067397/, ep2/.0067397/
117 c data a/6378.16/, ab2/1.0067397/, ep2/.0067397/
118 c update 2010
118 c update 2010
119 data a/6378.137/, ab2/1.0067396/, ep2/.0067396/
119 data a/6378.137/, ab2/1.0067396/, ep2/.0067396/
120 data dtr/.0174532925199/
120 data dtr/.0174532925199/
121
121
122 if (i .eq. 1) then
122 if (i .eq. 1) then
123
123
124 c .....geodetic to geocentric.....
124 c .....geodetic to geocentric.....
125
125
126 gdl = dtr*gdlat
126 gdl = dtr*gdlat
127 sinlat = sin(gdl)
127 sinlat = sin(gdl)
128 coslat = cos(gdl)
128 coslat = cos(gdl)
129 sl2 = sinlat*sinlat
129 sl2 = sinlat*sinlat
130 cl2 = ab2*coslat
130 cl2 = ab2*coslat
131 cl2 = cl2*cl2
131 cl2 = cl2*cl2
132 sinbet = sinlat/sqrt(sl2+cl2)
132 sinbet = sinlat/sqrt(sl2+cl2)
133 sb2 = amin1(sinbet*sinbet, 1.)
133 sb2 = amin1(sinbet*sinbet, 1.)
134 cosbet = sqrt(1. - sb2)
134 cosbet = sqrt(1. - sb2)
135 rgeoid = a/sqrt(1. + ep2*sb2)
135 rgeoid = a/sqrt(1. + ep2*sb2)
136 x = rgeoid*cosbet + gdalt*coslat
136 x = rgeoid*cosbet + gdalt*coslat
137 y = rgeoid*sinbet + gdalt*sinlat
137 y = rgeoid*sinbet + gdalt*sinlat
138 rkm = sqrt(x*x + y*y)
138 rkm = sqrt(x*x + y*y)
139 gclat = atan2(y,x)/dtr
139 gclat = atan2(y,x)/dtr
140 return
140 return
141
141
142 else if (i .eq. 2) then
142 else if (i .eq. 2) then
143
143
144 c .....geocentric to geodetic.....
144 c .....geocentric to geodetic.....
145
145
146 rer = rkm/a
146 rer = rkm/a
147 a2 = ((-1.4127348e-8/rer + .94339131e-8)/rer +
147 a2 = ((-1.4127348e-8/rer + .94339131e-8)/rer +
148 1 .33523288e-2)/rer
148 1 .33523288e-2)/rer
149 a4 = (((-1.2545063e-10/rer + .11760996e-9)/rer +
149 a4 = (((-1.2545063e-10/rer + .11760996e-9)/rer +
150 1 .11238084e-4)/rer - .2814244e-5)/rer
150 1 .11238084e-4)/rer - .2814244e-5)/rer
151 a6 = ((54.939685e-9/rer - 28.301730e-9)/rer +
151 a6 = ((54.939685e-9/rer - 28.301730e-9)/rer +
152 1 3.5435979e-9)/rer
152 1 3.5435979e-9)/rer
153 a8 = (((320./rer - 252.)/rer + 64.)/rer - 5.)
153 a8 = (((320./rer - 252.)/rer + 64.)/rer - 5.)
154 1 /rer*0.98008304e-12
154 1 /rer*0.98008304e-12
155 gcl = dtr*gclat
155 gcl = dtr*gclat
156 ccl = cos(gcl)
156 ccl = cos(gcl)
157 scl = sin(gcl)
157 scl = sin(gcl)
158 s2cl = 2.*scl*ccl
158 s2cl = 2.*scl*ccl
159 c2cl = 2.*ccl*ccl - 1.
159 c2cl = 2.*ccl*ccl - 1.
160 s4cl = 2.*s2cl*c2cl
160 s4cl = 2.*s2cl*c2cl
161 c4cl = 2.*c2cl*c2cl - 1.
161 c4cl = 2.*c2cl*c2cl - 1.
162 s8cl = 2.*s4cl*c4cl
162 s8cl = 2.*s4cl*c4cl
163 s6cl = s2cl*c4cl + c2cl*s4cl
163 s6cl = s2cl*c4cl + c2cl*s4cl
164 dltcl = s2cl*a2 + s4cl*a4 + s6cl*a6 + s8cl*a8
164 dltcl = s2cl*a2 + s4cl*a4 + s6cl*a6 + s8cl*a8
165 gdlat = gclat + dltcl/dtr
165 gdlat = gclat + dltcl/dtr
166 gdalt = rkm - a/sqrt(1.+ep2*scl*scl)
166 gdalt = rkm - a/sqrt(1.+ep2*scl*scl)
167 return
167 return
168
168
169 end if
169 end if
170 end
170 end
171 c
171 c
172 subroutine getshc2 (iu, fspec, nmax, erad, gh, ier)
172 subroutine getshc2 (iu, fspec, nmax, erad, gh, ier)
173
173
174 c ===============================================================
174 c ===============================================================
175 c
175 c
176 c version 1.01
176 c version 1.01
177 c
177 c
178 c reads spherical harmonic coefficients from the specified
178 c reads spherical harmonic coefficients from the specified
179 c file into an array.
179 c file into an array.
180 c
180 c
181 c input:
181 c input:
182 c iu - logical unit number
182 c iu - logical unit number
183 c fspec - file specification
183 c fspec - file specification
184 c
184 c
185 c output:
185 c output:
186 c nmax - maximum degree and order of model
186 c nmax - maximum degree and order of model
187 c erad - earth's radius associated with the spherical
187 c erad - earth's radius associated with the spherical
188 c harmonic coefficients, in the same units as
188 c harmonic coefficients, in the same units as
189 c elevation
189 c elevation
190 c gh - schmidt quasi-normal internal spherical
190 c gh - schmidt quasi-normal internal spherical
191 c harmonic coefficients
191 c harmonic coefficients
192 c ier - error number: = 0, no error
192 c ier - error number: = 0, no error
193 c = -2, records out of order
193 c = -2, records out of order
194 c = fortran run-time error number
194 c = fortran run-time error number
195 c
195 c
196 c a. zunde
196 c a. zunde
197 c usgs, ms 964, box 25046 federal center, denver, co 80225
197 c usgs, ms 964, box 25046 federal center, denver, co 80225
198 c
198 c
199 c ===============================================================
199 c ===============================================================
200
200
201 character fspec*(*)
201 character fspec*(*)
202 dimension gh(*)
202 dimension gh(*)
203
203
204 c ---------------------------------------------------------------
204 c ---------------------------------------------------------------
205 c open coefficient file. read past first header record.
205 c open coefficient file. read past first header record.
206 c read degree and order of model and earth's radius.
206 c read degree and order of model and earth's radius.
207 c ---------------------------------------------------------------
207 c ---------------------------------------------------------------
208
208
209
209
210 open (iu, file=fspec, status='old', iostat=ier, err=999 )
210 open (iu, file=fspec, status='old', iostat=ier, err=999 )
211
211
212
212
213 read (iu, *, iostat=ier, err=999)
213 read (iu, *, iostat=ier, err=999)
214 read (iu, *, iostat=ier, err=999) nmax, erad
214 read (iu, *, iostat=ier, err=999) nmax, erad
215
215
216
216
217 c ---------------------------------------------------------------
217 c ---------------------------------------------------------------
218 c read the coefficient file, arranged as follows:
218 c read the coefficient file, arranged as follows:
219 c
219 c
220 c n m g h
220 c n m g h
221 c ----------------------
221 c ----------------------
222 c / 1 0 gh(1) -
222 c / 1 0 gh(1) -
223 c / 1 1 gh(2) gh(3)
223 c / 1 1 gh(2) gh(3)
224 c / 2 0 gh(4) -
224 c / 2 0 gh(4) -
225 c / 2 1 gh(5) gh(6)
225 c / 2 1 gh(5) gh(6)
226 c nmax*(nmax+3)/2 / 2 2 gh(7) gh(8)
226 c nmax*(nmax+3)/2 / 2 2 gh(7) gh(8)
227 c records \ 3 0 gh(9) -
227 c records \ 3 0 gh(9) -
228 c \ . . . .
228 c \ . . . .
229 c \ . . . .
229 c \ . . . .
230 c nmax*(nmax+2) \ . . . .
230 c nmax*(nmax+2) \ . . . .
231 c elements in gh \ nmax nmax . .
231 c elements in gh \ nmax nmax . .
232 c
232 c
233 c n and m are, respectively, the degree and order of the
233 c n and m are, respectively, the degree and order of the
234 c coefficient.
234 c coefficient.
235 c ---------------------------------------------------------------
235 c ---------------------------------------------------------------
236
236
237 i = 0
237 i = 0
238 C for nn = 1, nmax
238 C for nn = 1, nmax
239 DO 1010 nn = 1, nmax
239 DO 1010 nn = 1, nmax
240 C for mm = 0, nn
240 C for mm = 0, nn
241 DO 1020 mm = 0, nn
241 DO 1020 mm = 0, nn
242 read (iu, *, iostat=ier, err=999) n, m, g, h
242 read (iu, *, iostat=ier, err=999) n, m, g, h
243 c write(*,*) n,m,g,h,nn,nmax,mm,ier
243 c write(*,*) n,m,g,h,nn,nmax,mm,ier
244 if (nn .ne. n .or. mm .ne. m) then
244 if (nn .ne. n .or. mm .ne. m) then
245 ier = -2
245 ier = -2
246 goto 999
246 goto 999
247 endif
247 endif
248 i = i + 1
248 i = i + 1
249 gh(i) = g
249 gh(i) = g
250 if (m .ne. 0) then
250 if (m .ne. 0) then
251 i = i + 1
251 i = i + 1
252 gh(i) = h
252 gh(i) = h
253 endif
253 endif
254 C end for
254 C end for
255 1020 CONTINUE
255 1020 CONTINUE
256 1021 CONTINUE
256 1021 CONTINUE
257 C end for
257 C end for
258 1010 CONTINUE
258 1010 CONTINUE
259 1011 CONTINUE
259 1011 CONTINUE
260
260
261 999 close (iu)
261 999 close (iu)
262
262
263 return
263 return
264 end
264 end
265 c
265 c
266 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ bfieldsr ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
266 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ bfieldsr ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
267 c
267 c
268 subroutine interpshc (date, dte1, nmax1, gh1, dte2,
268 subroutine interpshc (date, dte1, nmax1, gh1, dte2,
269 1 nmax2, gh2, nmax, gh)
269 1 nmax2, gh2, nmax, gh)
270
270
271 c ===============================================================
271 c ===============================================================
272 c
272 c
273 c version 1.01
273 c version 1.01
274 c
274 c
275 c interpolates linearly, in time, between two spherical
275 c interpolates linearly, in time, between two spherical
276 c harmonic models.
276 c harmonic models.
277 c
277 c
278 c input:
278 c input:
279 c date - date of resulting model (in decimal year)
279 c date - date of resulting model (in decimal year)
280 c dte1 - date of earlier model
280 c dte1 - date of earlier model
281 c nmax1 - maximum degree and order of earlier model
281 c nmax1 - maximum degree and order of earlier model
282 c gh1 - schmidt quasi-normal internal spherical
282 c gh1 - schmidt quasi-normal internal spherical
283 c harmonic coefficients of earlier model
283 c harmonic coefficients of earlier model
284 c dte2 - date of later model
284 c dte2 - date of later model
285 c nmax2 - maximum degree and order of later model
285 c nmax2 - maximum degree and order of later model
286 c gh2 - schmidt quasi-normal internal spherical
286 c gh2 - schmidt quasi-normal internal spherical
287 c harmonic coefficients of later model
287 c harmonic coefficients of later model
288 c
288 c
289 c output:
289 c output:
290 c gh - coefficients of resulting model
290 c gh - coefficients of resulting model
291 c nmax - maximum degree and order of resulting model
291 c nmax - maximum degree and order of resulting model
292 c
292 c
293 c a. zunde
293 c a. zunde
294 c usgs, ms 964, box 25046 federal center, denver, co 80225
294 c usgs, ms 964, box 25046 federal center, denver, co 80225
295 c
295 c
296 c ===============================================================
296 c ===============================================================
297
297
298 dimension gh1(*), gh2(*), gh(*)
298 dimension gh1(*), gh2(*), gh(*)
299
299
300 c ---------------------------------------------------------------
300 c ---------------------------------------------------------------
301 c the coefficients (gh) of the resulting model, at date
301 c the coefficients (gh) of the resulting model, at date
302 c date, are computed by linearly interpolating between the
302 c date, are computed by linearly interpolating between the
303 c coefficients of the earlier model (gh1), at date dte1,
303 c coefficients of the earlier model (gh1), at date dte1,
304 c and those of the later model (gh2), at date dte2. if one
304 c and those of the later model (gh2), at date dte2. if one
305 c model is smaller than the other, the interpolation is
305 c model is smaller than the other, the interpolation is
306 c performed with the missing coefficients assumed to be 0.
306 c performed with the missing coefficients assumed to be 0.
307 c ---------------------------------------------------------------
307 c ---------------------------------------------------------------
308
308
309 factor = (date - dte1) / (dte2 - dte1)
309 factor = (date - dte1) / (dte2 - dte1)
310
310
311 if (nmax1 .eq. nmax2) then
311 if (nmax1 .eq. nmax2) then
312 k = nmax1 * (nmax1 + 2)
312 k = nmax1 * (nmax1 + 2)
313 nmax = nmax1
313 nmax = nmax1
314 else if (nmax1 .gt. nmax2) then
314 else if (nmax1 .gt. nmax2) then
315 k = nmax2 * (nmax2 + 2)
315 k = nmax2 * (nmax2 + 2)
316 l = nmax1 * (nmax1 + 2)
316 l = nmax1 * (nmax1 + 2)
317 C for i = k + 1, l
317 C for i = k + 1, l
318 DO 1010 i = k + 1, l
318 DO 1010 i = k + 1, l
319 gh(i) = gh1(i) + factor * (-gh1(i))
319 gh(i) = gh1(i) + factor * (-gh1(i))
320 C end for
320 C end for
321 1010 CONTINUE
321 1010 CONTINUE
322 1011 CONTINUE
322 1011 CONTINUE
323 nmax = nmax1
323 nmax = nmax1
324 else
324 else
325 k = nmax1 * (nmax1 + 2)
325 k = nmax1 * (nmax1 + 2)
326 l = nmax2 * (nmax2 + 2)
326 l = nmax2 * (nmax2 + 2)
327 C for i = k + 1, l
327 C for i = k + 1, l
328 DO 1020 i = k + 1, l
328 DO 1020 i = k + 1, l
329 gh(i) = factor * gh2(i)
329 gh(i) = factor * gh2(i)
330 C end for
330 C end for
331 1020 CONTINUE
331 1020 CONTINUE
332 1021 CONTINUE
332 1021 CONTINUE
333 nmax = nmax2
333 nmax = nmax2
334 endif
334 endif
335
335
336 C for i = 1, k
336 C for i = 1, k
337 DO 1030 i = 1, k
337 DO 1030 i = 1, k
338 gh(i) = gh1(i) + factor * (gh2(i) - gh1(i))
338 gh(i) = gh1(i) + factor * (gh2(i) - gh1(i))
339 C end for
339 C end for
340 1030 CONTINUE
340 1030 CONTINUE
341 1031 CONTINUE
341 1031 CONTINUE
342
342
343 return
343 return
344 end
344 end
345 c
345 c
346 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ bfieldsr ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
346 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ bfieldsr ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
347 c
347 c
348 subroutine extrapshc (date, dte1, nmax1, gh1, nmax2,
348 subroutine extrapshc (date, dte1, nmax1, gh1, nmax2,
349 1 gh2, nmax, gh)
349 1 gh2, nmax, gh)
350
350
351 c ===============================================================
351 c ===============================================================
352 c
352 c
353 c version 1.01
353 c version 1.01
354 c
354 c
355 c extrapolates linearly a spherical harmonic model with a
355 c extrapolates linearly a spherical harmonic model with a
356 c rate-of-change model.
356 c rate-of-change model.
357 c
357 c
358 c input:
358 c input:
359 c date - date of resulting model (in decimal year)
359 c date - date of resulting model (in decimal year)
360 c dte1 - date of base model
360 c dte1 - date of base model
361 c nmax1 - maximum degree and order of base model
361 c nmax1 - maximum degree and order of base model
362 c gh1 - schmidt quasi-normal internal spherical
362 c gh1 - schmidt quasi-normal internal spherical
363 c harmonic coefficients of base model
363 c harmonic coefficients of base model
364 c nmax2 - maximum degree and order of rate-of-change
364 c nmax2 - maximum degree and order of rate-of-change
365 c model
365 c model
366 c gh2 - schmidt quasi-normal internal spherical
366 c gh2 - schmidt quasi-normal internal spherical
367 c harmonic coefficients of rate-of-change model
367 c harmonic coefficients of rate-of-change model
368 c
368 c
369 c output:
369 c output:
370 c gh - coefficients of resulting model
370 c gh - coefficients of resulting model
371 c nmax - maximum degree and order of resulting model
371 c nmax - maximum degree and order of resulting model
372 c
372 c
373 c a. zunde
373 c a. zunde
374 c usgs, ms 964, box 25046 federal center, denver, co 80225
374 c usgs, ms 964, box 25046 federal center, denver, co 80225
375 c
375 c
376 c ===============================================================
376 c ===============================================================
377
377
378 dimension gh1(*), gh2(*), gh(*)
378 dimension gh1(*), gh2(*), gh(*)
379
379
380 c ---------------------------------------------------------------
380 c ---------------------------------------------------------------
381 c the coefficients (gh) of the resulting model, at date
381 c the coefficients (gh) of the resulting model, at date
382 c date, are computed by linearly extrapolating the coef-
382 c date, are computed by linearly extrapolating the coef-
383 c ficients of the base model (gh1), at date dte1, using
383 c ficients of the base model (gh1), at date dte1, using
384 c those of the rate-of-change model (gh2), at date dte2. if
384 c those of the rate-of-change model (gh2), at date dte2. if
385 c one model is smaller than the other, the extrapolation is
385 c one model is smaller than the other, the extrapolation is
386 c performed with the missing coefficients assumed to be 0.
386 c performed with the missing coefficients assumed to be 0.
387 c ---------------------------------------------------------------
387 c ---------------------------------------------------------------
388
388
389 factor = (date - dte1)
389 factor = (date - dte1)
390
390
391 if (nmax1 .eq. nmax2) then
391 if (nmax1 .eq. nmax2) then
392 k = nmax1 * (nmax1 + 2)
392 k = nmax1 * (nmax1 + 2)
393 nmax = nmax1
393 nmax = nmax1
394 else if (nmax1 .gt. nmax2) then
394 else if (nmax1 .gt. nmax2) then
395 k = nmax2 * (nmax2 + 2)
395 k = nmax2 * (nmax2 + 2)
396 l = nmax1 * (nmax1 + 2)
396 l = nmax1 * (nmax1 + 2)
397 C for i = k + 1, l
397 C for i = k + 1, l
398 DO 1010 i = k + 1, l
398 DO 1010 i = k + 1, l
399 gh(i) = gh1(i)
399 gh(i) = gh1(i)
400 C end for
400 C end for
401 1010 CONTINUE
401 1010 CONTINUE
402 1011 CONTINUE
402 1011 CONTINUE
403 nmax = nmax1
403 nmax = nmax1
404 else
404 else
405 k = nmax1 * (nmax1 + 2)
405 k = nmax1 * (nmax1 + 2)
406 l = nmax2 * (nmax2 + 2)
406 l = nmax2 * (nmax2 + 2)
407 C for i = k + 1, l
407 C for i = k + 1, l
408 DO 1020 i = k + 1, l
408 DO 1020 i = k + 1, l
409 gh(i) = factor * gh2(i)
409 gh(i) = factor * gh2(i)
410 C end for
410 C end for
411 1020 CONTINUE
411 1020 CONTINUE
412 1021 CONTINUE
412 1021 CONTINUE
413 nmax = nmax2
413 nmax = nmax2
414 endif
414 endif
415
415
416 C for i = 1, k
416 C for i = 1, k
417 DO 1030 i = 1, k
417 DO 1030 i = 1, k
418 gh(i) = gh1(i) + factor * gh2(i)
418 gh(i) = gh1(i) + factor * gh2(i)
419 C end for
419 C end for
420 1030 CONTINUE
420 1030 CONTINUE
421 1031 CONTINUE
421 1031 CONTINUE
422
422
423 return
423 return
424 end
424 end
425 c
425 c
426 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ bfieldsr ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
426 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ bfieldsr ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
427 c
427 c
428 subroutine shval3 (igdgc, flat, flon, elev, erad, a2, b2,
428 subroutine shval3 (igdgc, flat, flon, elev, erad, a2, b2,
429 1 nmax, gh, iext, ext, x, y, z)
429 1 nmax, gh, iext, ext, x, y, z)
430
430
431 c ================================================================
431 c ================================================================
432 c
432 c
433 c version 1.01
433 c version 1.01
434 c
434 c
435 c calculates field components from spherical harmonic (sh)
435 c calculates field components from spherical harmonic (sh)
436 c models.
436 c models.
437 c
437 c
438 c input:
438 c input:
439 c igdgc - indicates coordinate system used; set equal to
439 c igdgc - indicates coordinate system used; set equal to
440 c 1 if geodetic, 2 if geocentric
440 c 1 if geodetic, 2 if geocentric
441 c flat - north latitude, in degrees
441 c flat - north latitude, in degrees
442 c flon - east longitude, in degrees
442 c flon - east longitude, in degrees
443 c elev - elevation above mean sea level (igdgc=1), or
443 c elev - elevation above mean sea level (igdgc=1), or
444 c radial distance from earth's center (igdgc=2)
444 c radial distance from earth's center (igdgc=2)
445 c erad - value of earth's radius associated with the sh
445 c erad - value of earth's radius associated with the sh
446 c coefficients, in same units as elev
446 c coefficients, in same units as elev
447 c a2,b2 - squares of semi-major and semi-minor axes of
447 c a2,b2 - squares of semi-major and semi-minor axes of
448 c the reference spheroid used for transforming
448 c the reference spheroid used for transforming
449 c between geodetic and geocentric coordinates or
449 c between geodetic and geocentric coordinates or
450 c components
450 c components
451 c nmax - maximum degree and order of coefficients
451 c nmax - maximum degree and order of coefficients
452 c gh - schmidt quasi-normal internal spherical
452 c gh - schmidt quasi-normal internal spherical
453 c harmonic coefficients
453 c harmonic coefficients
454 c iext - external coefficients flag (= 0 if none)
454 c iext - external coefficients flag (= 0 if none)
455 c ext - the three 1st-degree external coefficients
455 c ext - the three 1st-degree external coefficients
456 c (not used if iext = 0)
456 c (not used if iext = 0)
457 c
457 c
458 c output:
458 c output:
459 c x - northward component
459 c x - northward component
460 c y - eastward component
460 c y - eastward component
461 c z - vertically-downward component
461 c z - vertically-downward component
462 c
462 c
463 c based on subroutine "igrf" by d. r. barraclough and
463 c based on subroutine "igrf" by d. r. barraclough and
464 c s. r. c. malin, report no. 71/1, institute of geological
464 c s. r. c. malin, report no. 71/1, institute of geological
465 c sciences, u.k.
465 c sciences, u.k.
466 c
466 c
467 c norman w. peddie, u.s. geological survey, mail stop 964,
467 c norman w. peddie, u.s. geological survey, mail stop 964,
468 c federal center, box 25046, denver, colorado 80225
468 c federal center, box 25046, denver, colorado 80225
469 c
469 c
470 c ================================================================
470 c ================================================================
471
471
472 c the required sizes of the arrays used in this subroutine
472 c the required sizes of the arrays used in this subroutine
473 c depend on the value of nmax. the minimum dimensions
473 c depend on the value of nmax. the minimum dimensions
474 c needed are indicated in the table below. (note that this
474 c needed are indicated in the table below. (note that this
475 c version is dimensioned for nmax of 14 or less).
475 c version is dimensioned for nmax of 14 or less).
476 c
476 c
477 c minimum dimension
477 c minimum dimension
478 c --------------------------
478 c --------------------------
479 c nmax
479 c nmax
480 dimension sl(14), cl(14)
480 dimension sl(14), cl(14)
481 c (nmax * (nmax + 3)) / 2
481 c (nmax * (nmax + 3)) / 2
482 dimension p(119), q(119)
482 dimension p(119), q(119)
483 c nmax * (nmax + 2)
483 c nmax * (nmax + 2)
484 dimension gh(224)
484 dimension gh(224)
485 c 3
485 c 3
486 dimension ext(3)
486 dimension ext(3)
487
487
488 c ================================================================
488 c ================================================================
489
489
490 parameter (dtr = .01745329)
490 parameter (dtr = .01745329)
491
491
492 r = elev
492 r = elev
493 slat = sin (flat * dtr)
493 slat = sin (flat * dtr)
494 if (90. - flat .lt. .001) then
494 if (90. - flat .lt. .001) then
495 c 300 ft from n. pole
495 c 300 ft from n. pole
496 aa = 89.999
496 aa = 89.999
497 else if (90. + flat .lt. .001) then
497 else if (90. + flat .lt. .001) then
498 c 300 ft from s. pole
498 c 300 ft from s. pole
499 aa = -89.999
499 aa = -89.999
500 else
500 else
501 aa = flat
501 aa = flat
502 endif
502 endif
503 clat = cos (aa * dtr)
503 clat = cos (aa * dtr)
504 sl(1) = sin (flon * dtr)
504 sl(1) = sin (flon * dtr)
505 cl(1) = cos (flon * dtr)
505 cl(1) = cos (flon * dtr)
506 x = 0.
506 x = 0.
507 y = 0.
507 y = 0.
508 z = 0.
508 z = 0.
509 sd = 0.
509 sd = 0.
510 cd = 1.
510 cd = 1.
511 n = 0
511 n = 0
512 l = 1
512 l = 1
513 m = 1
513 m = 1
514 npq = (nmax * (nmax + 3)) / 2
514 npq = (nmax * (nmax + 3)) / 2
515
515
516 if (igdgc .eq. 1) then
516 if (igdgc .eq. 1) then
517 aa = a2 * clat * clat
517 aa = a2 * clat * clat
518 bb = b2 * slat * slat
518 bb = b2 * slat * slat
519 cc = aa + bb
519 cc = aa + bb
520 dd = sqrt (cc)
520 dd = sqrt (cc)
521 r = sqrt (elev * (elev + 2. * dd)
521 r = sqrt (elev * (elev + 2. * dd)
522 1 + (a2 * aa + b2 * bb) / cc)
522 1 + (a2 * aa + b2 * bb) / cc)
523 cd = (elev + dd) / r
523 cd = (elev + dd) / r
524 sd = (a2 - b2) / dd * slat * clat / r
524 sd = (a2 - b2) / dd * slat * clat / r
525
525
526 aa = slat
526 aa = slat
527 slat = slat * cd - clat * sd
527 slat = slat * cd - clat * sd
528 clat = clat * cd + aa * sd
528 clat = clat * cd + aa * sd
529 endif
529 endif
530
530
531 ratio = erad / r
531 ratio = erad / r
532
532
533 aa = sqrt (3.)
533 aa = sqrt (3.)
534 p(1) = 2. * slat
534 p(1) = 2. * slat
535 p(2) = 2. * clat
535 p(2) = 2. * clat
536 p(3) = 4.5 * slat * slat - 1.5
536 p(3) = 4.5 * slat * slat - 1.5
537 p(4) = 3. * aa * clat * slat
537 p(4) = 3. * aa * clat * slat
538 q(1) = -clat
538 q(1) = -clat
539 q(2) = slat
539 q(2) = slat
540 q(3) = -3. * clat * slat
540 q(3) = -3. * clat * slat
541 q(4) = aa * (slat * slat - clat * clat)
541 q(4) = aa * (slat * slat - clat * clat)
542
542
543 C for k = 1, npq
543 C for k = 1, npq
544 DO 1010 k = 1, npq
544 DO 1010 k = 1, npq
545 if (n .lt. m) then
545 if (n .lt. m) then
546 m = 0
546 m = 0
547 n = n + 1
547 n = n + 1
548 rr = ratio**(n + 2)
548 rr = ratio**(n + 2)
549 fn = n
549 fn = n
550 endif
550 endif
551 fm = m
551 fm = m
552 if (k .ge. 5) then
552 if (k .ge. 5) then
553 if (m .eq. n) then
553 if (m .eq. n) then
554 aa = sqrt (1. - .5 / fm)
554 aa = sqrt (1. - .5 / fm)
555 j = k - n - 1
555 j = k - n - 1
556 p(k) = (1. + 1. / fm) * aa * clat * p(j)
556 p(k) = (1. + 1. / fm) * aa * clat * p(j)
557 q(k) = aa * (clat * q(j) + slat / fm * p(j))
557 q(k) = aa * (clat * q(j) + slat / fm * p(j))
558 sl(m) = sl(m-1) * cl(1) + cl(m-1) * sl(1)
558 sl(m) = sl(m-1) * cl(1) + cl(m-1) * sl(1)
559 cl(m) = cl(m-1) * cl(1) - sl(m-1) * sl(1)
559 cl(m) = cl(m-1) * cl(1) - sl(m-1) * sl(1)
560 else
560 else
561 aa = sqrt (fn * fn - fm * fm)
561 aa = sqrt (fn * fn - fm * fm)
562 bb = sqrt ((fn - 1.)**2 - fm * fm) / aa
562 bb = sqrt ((fn - 1.)**2 - fm * fm) / aa
563 cc = (2. * fn - 1.) / aa
563 cc = (2. * fn - 1.) / aa
564 i = k - n
564 i = k - n
565 j = k - 2 * n + 1
565 j = k - 2 * n + 1
566 p(k) = (fn + 1.) * (cc * slat / fn * p(i) - bb
566 p(k) = (fn + 1.) * (cc * slat / fn * p(i) - bb
567 1 / (fn - 1.) * p(j))
567 1 / (fn - 1.) * p(j))
568 q(k) = cc * (slat * q(i) - clat / fn * p(i))
568 q(k) = cc * (slat * q(i) - clat / fn * p(i))
569 1 - bb * q(j)
569 1 - bb * q(j)
570 endif
570 endif
571 endif
571 endif
572
572
573 aa = rr * gh(l)
573 aa = rr * gh(l)
574
574
575 if (m .eq. 0) then
575 if (m .eq. 0) then
576 x = x + aa * q(k)
576 x = x + aa * q(k)
577 z = z - aa * p(k)
577 z = z - aa * p(k)
578 l = l + 1
578 l = l + 1
579 else
579 else
580 bb = rr * gh(l+1)
580 bb = rr * gh(l+1)
581 cc = aa * cl(m) + bb * sl(m)
581 cc = aa * cl(m) + bb * sl(m)
582 x = x + cc * q(k)
582 x = x + cc * q(k)
583 z = z - cc * p(k)
583 z = z - cc * p(k)
584 if (clat .gt. 0.) then
584 if (clat .gt. 0.) then
585 y = y + (aa * sl(m) - bb * cl(m)) * fm * p(k)
585 y = y + (aa * sl(m) - bb * cl(m)) * fm * p(k)
586 1 / ((fn + 1.) * clat)
586 1 / ((fn + 1.) * clat)
587 else
587 else
588 y = y + (aa * sl(m) - bb * cl(m)) * q(k)
588 y = y + (aa * sl(m) - bb * cl(m)) * q(k)
589 1 * slat
589 1 * slat
590 endif
590 endif
591 l = l + 2
591 l = l + 2
592 endif
592 endif
593 m = m + 1
593 m = m + 1
594 C end for
594 C end for
595 1010 CONTINUE
595 1010 CONTINUE
596 1011 CONTINUE
596 1011 CONTINUE
597
597
598 if (iext .ne. 0) then
598 if (iext .ne. 0) then
599 aa = ext(2) * cl(1) + ext(3) * sl(1)
599 aa = ext(2) * cl(1) + ext(3) * sl(1)
600 x = x - ext(1) * clat + aa * slat
600 x = x - ext(1) * clat + aa * slat
601 y = y + ext(2) * sl(1) - ext(3) * cl(1)
601 y = y + ext(2) * sl(1) - ext(3) * cl(1)
602 z = z + ext(1) * slat + aa * clat
602 z = z + ext(1) * slat + aa * clat
603 endif
603 endif
604
604
605 aa = x
605 aa = x
606 x = x * cd + z * sd
606 x = x * cd + z * sd
607 z = z * cd - aa * sd
607 z = z * cd - aa * sd
608
608
609
609
610 return
610 return
611 end
611 end
612 c
612 c
613 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ bfieldsr ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
613 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ bfieldsr ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
614 c
614 c
615 subroutine dihf (x, y, z, d, i, h, f)
615 subroutine dihf (x, y, z, d, i, h, f)
616
616
617 c ===============================================================
617 c ===============================================================
618 c
618 c
619 c version 1.01
619 c version 1.01
620 c
620 c
621 c computes the geomagnetic elements d, i, h, and f from
621 c computes the geomagnetic elements d, i, h, and f from
622 c x, y, and z.
622 c x, y, and z.
623 c
623 c
624 c input:
624 c input:
625 c x - northward component
625 c x - northward component
626 c y - eastward component
626 c y - eastward component
627 c z - vertically-downward component
627 c z - vertically-downward component
628 c
628 c
629 c output:
629 c output:
630 c d - declination
630 c d - declination
631 c i - inclination
631 c i - inclination
632 c h - horizontal intensity
632 c h - horizontal intensity
633 c f - total intensity
633 c f - total intensity
634 c
634 c
635 c a. zunde
635 c a. zunde
636 c usgs, ms 964, box 25046 federal center, denver, co 80225
636 c usgs, ms 964, box 25046 federal center, denver, co 80225
637 c
637 c
638 c ===============================================================
638 c ===============================================================
639
639
640 real i
640 real i
641 parameter ( sn = 0.0001 )
641 parameter ( sn = 0.0001 )
642
642
643 c ---------------------------------------------------------------
643 c ---------------------------------------------------------------
644 c if d and i cannot be determined, set equal to 999.0.
644 c if d and i cannot be determined, set equal to 999.0.
645 c ---------------------------------------------------------------
645 c ---------------------------------------------------------------
646
646
647 h2 = x*x + y*y
647 h2 = x*x + y*y
648 h = sqrt (h2)
648 h = sqrt (h2)
649 f = sqrt (h2 + z*z)
649 f = sqrt (h2 + z*z)
650 if (f .lt. sn) then
650 if (f .lt. sn) then
651 d = 999.
651 d = 999.
652 i = 999.
652 i = 999.
653 else
653 else
654 i = atan2 (z, h)
654 i = atan2 (z, h)
655 if (h .lt. sn) then
655 if (h .lt. sn) then
656 d = 999.
656 d = 999.
657 else
657 else
658 hpx = h + x
658 hpx = h + x
659 if (hpx .lt. sn) then
659 if (hpx .lt. sn) then
660 d = 180.
660 d = 180.
661 else
661 else
662 d = 2. * atan2 (y, hpx)
662 d = 2. * atan2 (y, hpx)
663 endif
663 endif
664 endif
664 endif
665 endif
665 endif
666
666
667 return
667 return
668 end
668 end
@@ -1,208 +1,222
1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 # All rights reserved.
2 # All rights reserved.
3 #
3 #
4 # Distributed under the terms of the BSD 3-clause license.
4 # Distributed under the terms of the BSD 3-clause license.
5 """schainpy is an open source library to read, write and process radar data
5 """schainpy is an open source library to read, write and process radar data
6
6
7 Signal Chain is a radar data processing library wich includes modules to read,
7 Signal Chain is a radar data processing library wich includes modules to read,
8 and write different files formats, besides modules to process and visualize the
8 and write different files formats, besides modules to process and visualize the
9 data.
9 data.
10 """
10 """
11
11
12 import os
12 import os
13 from setuptools import setup, Extension
13 from setuptools import setup, Extension
14 from setuptools.command.build_ext import build_ext as _build_ext
14 from setuptools.command.build_ext import build_ext as _build_ext
15 from schainpy import __version__
15 from schainpy import __version__
16
16
17 DOCLINES = __doc__.split("\n")
17 DOCLINES = __doc__.split("\n")
18
18
19 class build_ext(_build_ext):
19 class build_ext(_build_ext):
20 def finalize_options(self):
20 def finalize_options(self):
21 _build_ext.finalize_options(self)
21 _build_ext.finalize_options(self)
22 # Prevent numpy from thinking it is still in its setup process:
22 # Prevent numpy from thinking it is still in its setup process:
23 __builtins__.__NUMPY_SETUP__ = False
23 __builtins__.__NUMPY_SETUP__ = False
24 import numpy
24 import numpy
25 self.include_dirs.append(numpy.get_include())
25 self.include_dirs.append(numpy.get_include())
26
26
27 setup(
27 setup(
28 name = "schainpy",
28 name = "schainpy",
29 version = __version__,
29 version = __version__,
30 description = DOCLINES[0],
30 description = DOCLINES[0],
31 long_description = "\n".join(DOCLINES[2:]),
31 long_description = "\n".join(DOCLINES[2:]),
32 url = "https://github.com/JRO-Peru/schainpy",
32 url = "https://github.com/JRO-Peru/schainpy",
33 author = "Jicamarca Radio Observatory",
33 author = "Jicamarca Radio Observatory",
34 author_email = "jro-developers@jro.igp.gob.pe",
34 author_email = "jro-developers@jro.igp.gob.pe",
35 license="BSD-3-Clause",
35 license="BSD-3-Clause",
36 classifiers=[
36 classifiers=[
37 "Development Status :: 4 - Beta",
37 "Development Status :: 4 - Beta",
38 "Environment :: Console",
38 "Environment :: Console",
39 "Intended Audience :: Science/Research",
39 "Intended Audience :: Science/Research",
40 "License :: OSI Approved :: BSD License",
40 "License :: OSI Approved :: BSD License",
41 "Operating System :: MacOS :: MacOS X",
41 "Operating System :: MacOS :: MacOS X",
42 "Operating System :: POSIX :: Linux",
42 "Operating System :: POSIX :: Linux",
43 "Programming Language :: Python :: 2",
43 "Programming Language :: Python :: 2",
44 "Programming Language :: Python :: 2.7",
44 "Programming Language :: Python :: 2.7",
45 "Programming Language :: Python :: 3",
45 "Programming Language :: Python :: 3",
46 "Programming Language :: Python :: 3.5",
46 "Programming Language :: Python :: 3.5",
47 "Programming Language :: Python :: 3.6",
47 "Programming Language :: Python :: 3.6",
48 "Programming Language :: Python :: 3.7",
48 "Programming Language :: Python :: 3.7",
49 "Programming Language :: Python :: 3.8",
49 "Programming Language :: Python :: 3.8",
50 "Topic :: Scientific/Engineering",
50 "Topic :: Scientific/Engineering",
51 ],
51 ],
52 packages = {
52 packages = {
53 'schainpy',
53 'schainpy',
54 'schainpy.model',
54 'schainpy.model',
55 'schainpy.model.data',
55 'schainpy.model.data',
56 'schainpy.model.graphics',
56 'schainpy.model.graphics',
57 'schainpy.model.io',
57 'schainpy.model.io',
58 'schainpy.model.proc',
58 'schainpy.model.proc',
59 'schainpy.model.utils',
59 'schainpy.model.utils',
60 'schainpy.utils',
60 'schainpy.utils',
61 'schainpy.gui',
61 'schainpy.gui',
62 'schainpy.cli',
62 'schainpy.cli',
63 },
63 },
64 package_data = {'': ['schain.conf.template'],
64 package_data = {'': ['schain.conf.template'],
65 'schainpy.files': ['*.oga']
65 'schainpy.files': ['*.oga']
66 },
66 },
67 include_package_data = False,
67 include_package_data = False,
68 scripts = ['schainpy/gui/schainGUI'],
68 scripts = ['schainpy/gui/schainGUI'],
69 entry_points = {
69 entry_points = {
70 'console_scripts': [
70 'console_scripts': [
71 'schain = schainpy.cli.cli:main',
71 'schain = schainpy.cli.cli:main',
72 ],
72 ],
73 },
73 },
74 cmdclass = {'build_ext': build_ext},
74 cmdclass = {'build_ext': build_ext},
75 ext_modules=[
75 ext_modules=[
76 Extension("schainpy.model.data._noise", ["schainc/_noise.c"]),
76 Extension("schainpy.model.data._noise", ["schainc/_noise.c"]),
77 Extension("schainpy.model.data._HS_algorithm", ["schainc/_HS_algorithm.c"]),
77 Extension("schainpy.model.data._HS_algorithm", ["schainc/_HS_algorithm.c"]),
78 ],
78 ],
79 setup_requires = ["numpy"],
79 setup_requires = ["numpy"],
80 install_requires = [
80 install_requires = [
81 "scipy",
81 "scipy",
82 "h5py",
82 "h5py",
83 "matplotlib",
83 "matplotlib",
84 "pyzmq",
84 "pyzmq",
85 "fuzzywuzzy",
85 "fuzzywuzzy",
86 "click",
86 "click",
87 ],
87 ],
88 )
88 )
89
89
90 from numpy.distutils.core import Extension, setup
90 from numpy.distutils.core import Extension, setup
91
91
92 setup(name='schainpy',
92 setup(name='schainpy',
93 ext_modules = [
93 ext_modules = [
94 Extension("schainpy.model.proc.mkfact_short_2020_2",
94 Extension("schainpy.model.proc.mkfact_short_2020_2",
95 sources=[
95 sources=[
96 "schainf/mkfact/mkfact_short_2020_2.pyf",
96 "schainf/mkfact/mkfact_short_2020_2.pyf",
97 "schainf/mkfact/lmdif1.f",
97 "schainf/mkfact/lmdif1.f",
98 "schainf/mkfact/mkfact.f",
98 "schainf/mkfact/mkfact.f",
99 "schainf/mkfact/r1mach.f",
99 "schainf/mkfact/r1mach.f",
100 "schainf/mkfact/bfield2.f"])
101 ]
102 )
103
104 '''
105 setup(name='schainpy',
106 ext_modules = [
107 Extension("schainpy.model.proc.mkfact_short_2020_2",
108 sources=[
109 "schainf/mkfact/mkfact_short_2020_2.pyf",
110 "schainf/mkfact/lmdif1.f",
111 "schainf/mkfact/mkfact.f",
112 "schainf/mkfact/r1mach.f",
100 "schainf/mkfact/bfield2.f"]),
113 "schainf/mkfact/bfield2.f"]),
101 Extension("schainpy.model.proc.full_profile_profile",
114 Extension("schainpy.model.proc.full_profile_profile",
102 sources=[
115 sources=[
103 "schainf/full_profile/full_profile_profile.pyf",
116 "schainf/full_profile/full_profile_profile.pyf",
104 "schainf/full_profile/full_profile_profile.f",
117 "schainf/full_profile/full_profile_profile.f",
105 "schainf/full_profile/fitacf.f",
118 "schainf/full_profile/fitacf.f",
106 "schainf/full_profile/r1mach.f",
119 "schainf/full_profile/r1mach.f",
107 "schainf/full_profile/lmdif1.f",
120 "schainf/full_profile/lmdif1.f",
108 "schainf/full_profile/lagp.f",
121 "schainf/full_profile/lagp.f",
109 "schainf/full_profile/reader.c",
122 "schainf/full_profile/reader.c",
110 "schainf/full_profile/cbesi.f",
123 "schainf/full_profile/cbesi.f",
111 "schainf/full_profile/i1mach.f",
124 "schainf/full_profile/i1mach.f",
112 "schainf/full_profile/zeta.f",
125 "schainf/full_profile/zeta.f",
113 "schainf/full_profile/qc25f.f",
126 "schainf/full_profile/qc25f.f",
114 "schainf/full_profile/qwgtf.f",
127 "schainf/full_profile/qwgtf.f",
115 "schainf/full_profile/qcheb.f",
128 "schainf/full_profile/qcheb.f",
116 "schainf/full_profile/sgtsl.f",
129 "schainf/full_profile/sgtsl.f",
117 "schainf/full_profile/qk15w.f",
130 "schainf/full_profile/qk15w.f",
118 "schainf/full_profile/complex.c",
131 "schainf/full_profile/complex.c",
119 "schainf/full_profile/cbinu.f",
132 "schainf/full_profile/cbinu.f",
120 "schainf/full_profile/cseri.f",
133 "schainf/full_profile/cseri.f",
121 "schainf/full_profile/cwrsk.f",
134 "schainf/full_profile/cwrsk.f",
122 "schainf/full_profile/crati.f",
135 "schainf/full_profile/crati.f",
123 "schainf/full_profile/casyi.f",
136 "schainf/full_profile/casyi.f",
124 "schainf/full_profile/cbuni.f",
137 "schainf/full_profile/cbuni.f",
125 "schainf/full_profile/cuni2.f",
138 "schainf/full_profile/cuni2.f",
126 "schainf/full_profile/gamln.f",
139 "schainf/full_profile/gamln.f",
127 "schainf/full_profile/cuchk.f",
140 "schainf/full_profile/cuchk.f",
128 "schainf/full_profile/cbknu.f",
141 "schainf/full_profile/cbknu.f",
129 "schainf/full_profile/cshch.f",
142 "schainf/full_profile/cshch.f",
130 "schainf/full_profile/ckscl.f",
143 "schainf/full_profile/ckscl.f",
131 "schainf/full_profile/cuoik.f",
144 "schainf/full_profile/cuoik.f",
132 "schainf/full_profile/cunik.f",
145 "schainf/full_profile/cunik.f",
133 "schainf/full_profile/cuni1.f",
146 "schainf/full_profile/cuni1.f",
134 "schainf/full_profile/cairy.f",
147 "schainf/full_profile/cairy.f",
135 "schainf/full_profile/cmlri.f",
148 "schainf/full_profile/cmlri.f",
136 "schainf/full_profile/cunhj.f",
149 "schainf/full_profile/cunhj.f",
137 "schainf/full_profile/cacai.f",
150 "schainf/full_profile/cacai.f",
138 "schainf/full_profile/csisl.f",
151 "schainf/full_profile/csisl.f",
139 "schainf/full_profile/caxpy.f",
152 "schainf/full_profile/caxpy.f",
140 "schainf/full_profile/cs1s2.f",
153 "schainf/full_profile/cs1s2.f",
141 "schainf/full_profile/scabs1.f",
154 "schainf/full_profile/scabs1.f",
142 "schainf/full_profile/cdotu.f",
155 "schainf/full_profile/cdotu.f",
143 "schainf/full_profile/rs.f",
156 "schainf/full_profile/rs.f",
144 "schainf/full_profile/sppfa.f",
157 "schainf/full_profile/sppfa.f",
145 "schainf/full_profile/sdot.f",
158 "schainf/full_profile/sdot.f",
146 "schainf/full_profile/tred2.f",
159 "schainf/full_profile/tred2.f",
147 "schainf/full_profile/tql2.f",
160 "schainf/full_profile/tql2.f",
148 "schainf/full_profile/sppdi.f",
161 "schainf/full_profile/sppdi.f",
149 "schainf/full_profile/saxpy.f",
162 "schainf/full_profile/saxpy.f",
150 "schainf/full_profile/sscal.f",
163 "schainf/full_profile/sscal.f",
151 "schainf/full_profile/pythag.f",
164 "schainf/full_profile/pythag.f",
152 "schainf/full_profile/tql1.f",
165 "schainf/full_profile/tql1.f",
153 "schainf/full_profile/tred1.f"]),
166 "schainf/full_profile/tred1.f"]),
154 Extension("schainpy.model.proc.fitacf_acf2",
167 Extension("schainpy.model.proc.fitacf_acf2",
155 sources = [
168 sources = [
156 "schainf/acf2/fitacf_acf2.pyf",
169 "schainf/acf2/fitacf_acf2.pyf",
157 "schainf/acf2/full_profile_profile.f",
170 "schainf/acf2/full_profile_profile.f",
158 "schainf/acf2/fitacf.f",
171 "schainf/acf2/fitacf.f",
159 "schainf/acf2/r1mach.f",
172 "schainf/acf2/r1mach.f",
160 "schainf/acf2/lmdif1.f",
173 "schainf/acf2/lmdif1.f",
161 "schainf/acf2/lagp.f",
174 "schainf/acf2/lagp.f",
162 "schainf/acf2/reader.c",
175 "schainf/acf2/reader.c",
163 "schainf/acf2/cbesi.f",
176 "schainf/acf2/cbesi.f",
164 "schainf/acf2/i1mach.f",
177 "schainf/acf2/i1mach.f",
165 "schainf/acf2/zeta.f",
178 "schainf/acf2/zeta.f",
166 "schainf/acf2/qc25f.f",
179 "schainf/acf2/qc25f.f",
167 "schainf/acf2/qwgtf.f",
180 "schainf/acf2/qwgtf.f",
168 "schainf/acf2/qcheb.f",
181 "schainf/acf2/qcheb.f",
169 "schainf/acf2/sgtsl.f",
182 "schainf/acf2/sgtsl.f",
170 "schainf/acf2/qk15w.f",
183 "schainf/acf2/qk15w.f",
171 "schainf/acf2/complex.c",
184 "schainf/acf2/complex.c",
172 "schainf/acf2/cbinu.f",
185 "schainf/acf2/cbinu.f",
173 "schainf/acf2/cseri.f",
186 "schainf/acf2/cseri.f",
174 "schainf/acf2/cwrsk.f",
187 "schainf/acf2/cwrsk.f",
175 "schainf/acf2/crati.f",
188 "schainf/acf2/crati.f",
176 "schainf/acf2/casyi.f",
189 "schainf/acf2/casyi.f",
177 "schainf/acf2/cbuni.f",
190 "schainf/acf2/cbuni.f",
178 "schainf/acf2/cuni2.f",
191 "schainf/acf2/cuni2.f",
179 "schainf/acf2/gamln.f",
192 "schainf/acf2/gamln.f",
180 "schainf/acf2/cuchk.f",
193 "schainf/acf2/cuchk.f",
181 "schainf/acf2/cbknu.f",
194 "schainf/acf2/cbknu.f",
182 "schainf/acf2/cshch.f",
195 "schainf/acf2/cshch.f",
183 "schainf/acf2/ckscl.f",
196 "schainf/acf2/ckscl.f",
184 "schainf/acf2/cuoik.f",
197 "schainf/acf2/cuoik.f",
185 "schainf/acf2/cunik.f",
198 "schainf/acf2/cunik.f",
186 "schainf/acf2/cuni1.f",
199 "schainf/acf2/cuni1.f",
187 "schainf/acf2/cairy.f",
200 "schainf/acf2/cairy.f",
188 "schainf/acf2/cmlri.f",
201 "schainf/acf2/cmlri.f",
189 "schainf/acf2/cunhj.f",
202 "schainf/acf2/cunhj.f",
190 "schainf/acf2/cacai.f",
203 "schainf/acf2/cacai.f",
191 "schainf/acf2/csisl.f",
204 "schainf/acf2/csisl.f",
192 "schainf/acf2/caxpy.f",
205 "schainf/acf2/caxpy.f",
193 "schainf/acf2/cs1s2.f",
206 "schainf/acf2/cs1s2.f",
194 "schainf/acf2/scabs1.f",
207 "schainf/acf2/scabs1.f",
195 "schainf/acf2/cdotu.f",
208 "schainf/acf2/cdotu.f",
196 "schainf/acf2/rs.f",
209 "schainf/acf2/rs.f",
197 "schainf/acf2/sppfa.f",
210 "schainf/acf2/sppfa.f",
198 "schainf/acf2/sdot.f",
211 "schainf/acf2/sdot.f",
199 "schainf/acf2/tred2.f",
212 "schainf/acf2/tred2.f",
200 "schainf/acf2/tql2.f",
213 "schainf/acf2/tql2.f",
201 "schainf/acf2/sppdi.f",
214 "schainf/acf2/sppdi.f",
202 "schainf/acf2/saxpy.f",
215 "schainf/acf2/saxpy.f",
203 "schainf/acf2/sscal.f",
216 "schainf/acf2/sscal.f",
204 "schainf/acf2/pythag.f",
217 "schainf/acf2/pythag.f",
205 "schainf/acf2/tql1.f",
218 "schainf/acf2/tql1.f",
206 "schainf/acf2/tred1.f"])
219 "schainf/acf2/tred1.f"])
207 ]
220 ]
208 )
221 )
222 '''
General Comments 0
You need to be logged in to leave comments. Login now