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