/[PAMELA software]/quicklook/SatelliteInclination/src/geopack.f
ViewVC logotype

Contents of /quicklook/SatelliteInclination/src/geopack.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (show annotations) (download) (vendor branch)
Thu Feb 8 00:49:33 2007 UTC (17 years, 10 months ago) by cafagna
Branch: MAIN, first
CVS Tags: v1r0, HEAD
Changes since 1.1: +0 -0 lines
Firse release of the Satellite inclination quicklook

1 c<pre>
2 c
3 c ##########################################################################
4 c # #
5 c # GEOPACK-2005 #
6 c # (MAIN SET OF FORTRAN CODES) #
7 c # #
8 c ##########################################################################
9 C
10 c
11 c This collection of subroutines is a result of several upgrades of the original package
12 c written by N. A. Tsyganenko in 1978-1979. This version is dated May 04, 2005. On that
13 c date, the IGRF coefficients were updated according to the recently published table of
14 c IGRF-10 coefficients, so that the main field model now extends through 2010 (a linear
15 c extrapolation is used for 2005 - 2010, based on the table of secular velocities). For
16 c more details, see http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html (revision of 03/22/2005).
17 c
18 c
19 c Prefatory notes to the version of April 22, 2003:
20 c
21 c This package represents an in-depth revision of the previous version, with significant
22 c changes in the format of calling statements. Users should familiarize themselves with
23 c the new formats and rules, and accordingly adjust their source codes, as specified
24 c below. Please consult the documentation file geopack-2005.doc (also available from this
25 c site) for detailed descriptions of individual subroutines.
26 c
27 c The following changes were made to the previous release of GEOPACK (of Jan 5, 2001).
28 c
29 c (1) Subroutine IGRF, calculating the Earth's main field:
30
31 c (a) Two versions of this subroutine are provided here. In the first one (IGRF_GSM)
32 c both input (position) and output (field components) are in the Geocentric Solar-
33 c Magnetospheric Cartesian coordinates, while the second one (IGRF_GEO) uses sphe-
34 c rical geographical (geocentric) coordinates, as in the older releases.
35
36 c (b) updating of all expansion coefficients is now made separately in the s/r RECALC,
37 c which also takes into account the secular change of the coefficients within
38 c a given year (at the Earth's surface, the rate of the change can reach 7 nT/month).
39
40 c (c) the optimal length of spherical harmonic expansions is now automatically set
41 c inside the code, based on the radial distance, so that the deviation from the
42 c full-length approximation does not exceed 0.01 nT. (In the previous versions,
43 c the upper limit NM of the order of harmonics had to be specified by users),
44 c
45 c (2) Subroutine DIP, calculating the Earth's field in the dipole approximation:
46
47 c (a) no longer accepts the tilt angle via the list of formal parameters. Instead,
48 c the sine SPS and cosine CPS of that angle are now forwarded into DIP via the
49 c first common block /GEOPACK1/. Accordingly, there are two options: (i) to
50 c calculate SPS and CPS by calling RECALC before calling DIP, or (ii) to specify
51 c them explicitly. In the last case, SPS and CPS should be specified AFTER the
52 c invocation of RECALC (otherwise they will be overridden by those returned by
53 c RECALC).
54
55 c (b) the Earth's dipole moment is now calculated by RECALC, based on the table of
56 c the IGRF coefficients and their secular variation rates, for a given year and
57 c the day of the year, and the obtained value of the moment is forwarded into DIP
58 c via the second common block /GEOPACK2/. (In the previous versions, only a single
59 c fixed value was provided for the geodipole moment, corresponding to the most
60 c recent epoch).
61 c
62 c (3) Subroutine RECALC now consolidates in one module all calculations needed to
63 c initialize and update the values of coefficients and quantities that vary in
64 c time, either due to secular changes of the main geomagnetic field or as a result
65 c of Earth's diurnal rotation and orbital motion around Sun. That allowed us to
66 c simplify the codes and make them more compiler-independent.
67 c
68 c (4) Subroutine GEOMAG is now identical in its structure to other coordinate trans-
69 c formation subroutines. It no longer invokes RECALC from within GEOMAG, but uses
70 c precalculated values of the rotation matrix elements, obtained by a separate
71 c external invocation of RECALC. This eliminates possible interference of the
72 c two subroutines in the old version of the package.
73 c
74 c (5) Subroutine TRACE (and the subsidiary modules STEP and RHAND):
75 c
76 c (a) no longer needs to specify the highest order of spherical harmonics in the
77 c main geomagnetic field expansion - it is now calculated automatically inside the
78 c IGRF_GSM (or IGRF_GEO) subroutine.
79 c
80 c (b) the internal field model can now be explicitly chosen by specifying the para-
81 c meter INNAME (either IGRF_GSM or DIP).
82 c
83 c (6) A new subroutine BCARSP was added, providing a conversion of Cartesian field
84 c components into spherical ones (operation, inverse to that performed by the sub-
85 c routine BSPCAR).
86 c
87 c (7) Two new subroutines were added, SHUETAL_MGNP and T96_MGNP, providing the position
88 c of the magnetopause, according to the model of Shue et al. [1998] and the one
89 c used in the T96 magnetospheric magnetic field model.
90 c
91 c
92 c----------------------------------------------------------------------------------
93 c
94 SUBROUTINE IGRF_GSM (XGSM,YGSM,ZGSM,HXGSM,HYGSM,HZGSM)
95 c
96 C CALCULATES COMPONENTS OF THE MAIN (INTERNAL) GEOMAGNETIC FIELD IN THE GEOCENTRIC SOLAR
97 C MAGNETOSPHERIC COORDINATE SYSTEM, USING IAGA INTERNATIONAL GEOMAGNETIC REFERENCE MODEL
98 C COEFFICIENTS (e.g., http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html Revised: 22 March, 2005)
99 c
100 C
101 C BEFORE THE FIRST CALL OF THIS SUBROUTINE, OR IF THE DATE/TIME (IYEAR,IDAY,IHOUR,MIN,ISEC)
102 C WAS CHANGED, THE MODEL COEFFICIENTS AND GEO-GSM ROTATION MATRIX ELEMENTS SHOULD BE UPDATED
103 c BY CALLING THE SUBROUTINE RECALC
104 C
105 C-----INPUT PARAMETERS:
106 C
107 C XGSM,YGSM,ZGSM - CARTESIAN GSM COORDINATES (IN UNITS RE=6371.2 KM)
108 C
109 C-----OUTPUT PARAMETERS:
110 C
111 C HXGSM,HYGSM,HZGSM - CARTESIAN GSM COMPONENTS OF THE MAIN GEOMAGNETIC FIELD IN NANOTESLA
112 C
113 C LAST MODIFICATION: MAY 4, 2005.
114 C THIS VERSION OF THE CODE ACCEPTS DATES FROM 1965 THROUGH 2010.
115 c
116 C AUTHOR: N. A. TSYGANENKO
117 C
118 C
119 COMMON /GEOPACK2/ G(105),H(105),REC(105)
120
121 DIMENSION A(14),B(14)
122
123 CALL GEOGSM (XGEO,YGEO,ZGEO,XGSM,YGSM,ZGSM,-1)
124 RHO2=XGEO**2+YGEO**2
125 R=SQRT(RHO2+ZGEO**2)
126 C=ZGEO/R
127 RHO=SQRT(RHO2)
128 S=RHO/R
129 IF (S.LT.1.E-5) THEN
130 CF=1.
131 SF=0.
132 ELSE
133 CF=XGEO/RHO
134 SF=YGEO/RHO
135 ENDIF
136 C
137 PP=1./R
138 P=PP
139 C
140 C IN THIS NEW VERSION, THE OPTIMAL VALUE OF THE PARAMETER NM (MAXIMAL ORDER OF THE SPHERICAL
141 C HARMONIC EXPANSION) IS NOT USER-PRESCRIBED, BUT CALCULATED INSIDE THE SUBROUTINE, BASED
142 C ON THE VALUE OF THE RADIAL DISTANCE R:
143 C
144 IRP3=R+2
145 NM=3+30/IRP3
146 IF (NM.GT.13) NM=13
147
148 K=NM+1
149 DO 150 N=1,K
150 P=P*PP
151 A(N)=P
152 150 B(N)=P*N
153
154 P=1.
155 D=0.
156 BBR=0.
157 BBT=0.
158 BBF=0.
159
160 DO 200 M=1,K
161 IF(M.EQ.1) GOTO 160
162 MM=M-1
163 W=X
164 X=W*CF+Y*SF
165 Y=Y*CF-W*SF
166 GOTO 170
167 160 X=0.
168 Y=1.
169 170 Q=P
170 Z=D
171 BI=0.
172 P2=0.
173 D2=0.
174 DO 190 N=M,K
175 AN=A(N)
176 MN=N*(N-1)/2+M
177 E=G(MN)
178 HH=H(MN)
179 W=E*Y+HH*X
180 BBR=BBR+B(N)*W*Q
181 BBT=BBT-AN*W*Z
182 IF(M.EQ.1) GOTO 180
183 QQ=Q
184 IF(S.LT.1.E-5) QQ=Z
185 BI=BI+AN*(E*X-HH*Y)*QQ
186 180 XK=REC(MN)
187 DP=C*Z-S*Q-XK*D2
188 PM=C*Q-XK*P2
189 D2=Z
190 P2=Q
191 Z=DP
192 190 Q=PM
193 D=S*D+C*P
194 P=S*P
195 IF(M.EQ.1) GOTO 200
196 BI=BI*MM
197 BBF=BBF+BI
198 200 CONTINUE
199 C
200 BR=BBR
201 BT=BBT
202 IF(S.LT.1.E-5) GOTO 210
203 BF=BBF/S
204 GOTO 211
205 210 IF(C.LT.0.) BBF=-BBF
206 BF=BBF
207
208 211 HE=BR*S+BT*C
209 HXGEO=HE*CF-BF*SF
210 HYGEO=HE*SF+BF*CF
211 HZGEO=BR*C-BT*S
212
213 CALL GEOGSM (HXGEO,HYGEO,HZGEO,HXGSM,HYGSM,HZGSM,1)
214
215 RETURN
216 END
217 C
218 c==========================================================================================
219 C
220 c
221 SUBROUTINE IGRF_GEO (R,THETA,PHI,BR,BTHETA,BPHI)
222 c
223 C CALCULATES COMPONENTS OF THE MAIN (INTERNAL) GEOMAGNETIC FIELD IN THE SPHERICAL GEOGRAPHIC
224 C (GEOCENTRIC) COORDINATE SYSTEM, USING IAGA INTERNATIONAL GEOMAGNETIC REFERENCE MODEL
225 C COEFFICIENTS (e.g., http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html, revised: 22 March, 2005)
226 C
227 C BEFORE THE FIRST CALL OF THIS SUBROUTINE, OR IF THE DATE (IYEAR AND IDAY) WAS CHANGED,
228 C THE MODEL COEFFICIENTS SHOULD BE UPDATED BY CALLING THE SUBROUTINE RECALC
229 C
230 C-----INPUT PARAMETERS:
231 C
232 C R, THETA, PHI - SPHERICAL GEOGRAPHIC (GEOCENTRIC) COORDINATES:
233 C RADIAL DISTANCE R IN UNITS RE=6371.2 KM, COLATITUDE THETA AND LONGITUDE PHI IN RADIANS
234 C
235 C-----OUTPUT PARAMETERS:
236 C
237 C BR, BTHETA, BPHI - SPHERICAL COMPONENTS OF THE MAIN GEOMAGNETIC FIELD IN NANOTESLA
238 C (POSITIVE BR OUTWARD, BTHETA SOUTHWARD, BPHI EASTWARD)
239 C
240 C LAST MODIFICATION: MAY 4, 2005.
241 C THIS VERSION OF THE CODE ACCEPTS DATES FROM 1965 THROUGH 2010.
242 c
243 C AUTHOR: N. A. TSYGANENKO
244 C
245 C
246 COMMON /GEOPACK2/ G(105),H(105),REC(105)
247
248 DIMENSION A(14),B(14)
249
250 C=COS(THETA)
251 S=SIN(THETA)
252 CF=COS(PHI)
253 SF=SIN(PHI)
254 C
255 PP=1./R
256 P=PP
257 C
258 C IN THIS NEW VERSION, THE OPTIMAL VALUE OF THE PARAMETER NM (MAXIMAL ORDER OF THE SPHERICAL
259 C HARMONIC EXPANSION) IS NOT USER-PRESCRIBED, BUT CALCULATED INSIDE THE SUBROUTINE, BASED
260 C ON THE VALUE OF THE RADIAL DISTANCE R:
261 C
262 IRP3=R+2
263 NM=3+30/IRP3
264 IF (NM.GT.13) NM=13
265
266 K=NM+1
267 DO 150 N=1,K
268 P=P*PP
269 A(N)=P
270 150 B(N)=P*N
271
272 P=1.
273 D=0.
274 BBR=0.
275 BBT=0.
276 BBF=0.
277
278 DO 200 M=1,K
279 IF(M.EQ.1) GOTO 160
280 MM=M-1
281 W=X
282 X=W*CF+Y*SF
283 Y=Y*CF-W*SF
284 GOTO 170
285 160 X=0.
286 Y=1.
287 170 Q=P
288 Z=D
289 BI=0.
290 P2=0.
291 D2=0.
292 DO 190 N=M,K
293 AN=A(N)
294 MN=N*(N-1)/2+M
295 E=G(MN)
296 HH=H(MN)
297 W=E*Y+HH*X
298 BBR=BBR+B(N)*W*Q
299 BBT=BBT-AN*W*Z
300 IF(M.EQ.1) GOTO 180
301 QQ=Q
302 IF(S.LT.1.E-5) QQ=Z
303 BI=BI+AN*(E*X-HH*Y)*QQ
304 180 XK=REC(MN)
305 DP=C*Z-S*Q-XK*D2
306 PM=C*Q-XK*P2
307 D2=Z
308 P2=Q
309 Z=DP
310 190 Q=PM
311 D=S*D+C*P
312 P=S*P
313 IF(M.EQ.1) GOTO 200
314 BI=BI*MM
315 BBF=BBF+BI
316 200 CONTINUE
317 C
318 BR=BBR
319 BTHETA=BBT
320 IF(S.LT.1.E-5) GOTO 210
321 BPHI=BBF/S
322 RETURN
323 210 IF(C.LT.0.) BBF=-BBF
324 BPHI=BBF
325
326 RETURN
327 END
328 C
329 c==========================================================================================
330 c
331 SUBROUTINE DIP (XGSM,YGSM,ZGSM,BXGSM,BYGSM,BZGSM)
332 C
333 C CALCULATES GSM COMPONENTS OF A GEODIPOLE FIELD WITH THE DIPOLE MOMENT
334 C CORRESPONDING TO THE EPOCH, SPECIFIED BY CALLING SUBROUTINE RECALC (SHOULD BE
335 C INVOKED BEFORE THE FIRST USE OF THIS ONE AND IN CASE THE DATE/TIME WAS CHANGED).
336 C
337 C--INPUT PARAMETERS: XGSM,YGSM,ZGSM - GSM COORDINATES IN RE (1 RE = 6371.2 km)
338 C
339 C--OUTPUT PARAMETERS: BXGSM,BYGSM,BZGSM - FIELD COMPONENTS IN GSM SYSTEM, IN NANOTESLA.
340 C
341 C LAST MODIFICATION: MAY 4, 2005
342 C
343 C AUTHOR: N. A. TSYGANENKO
344 C
345 COMMON /GEOPACK1/ AAA(10),SPS,CPS,BBB(23)
346 COMMON /GEOPACK2/ G(105),H(105),REC(105)
347
348 DIPMOM=SQRT(G(2)**2+G(3)**2+H(3)**2)
349
350 P=XGSM**2
351 U=ZGSM**2
352 V=3.*ZGSM*XGSM
353 T=YGSM**2
354 Q=DIPMOM/SQRT(P+T+U)**5
355 BXGSM=Q*((T+U-2.*P)*SPS-V*CPS)
356 BYGSM=-3.*YGSM*Q*(XGSM*SPS+ZGSM*CPS)
357 BZGSM=Q*((P+T-2.*U)*CPS-V*SPS)
358 RETURN
359 END
360
361 C*******************************************************************
362 c
363 SUBROUTINE SUN (IYEAR,IDAY,IHOUR,MIN,ISEC,GST,SLONG,SRASN,SDEC)
364 C
365 C CALCULATES FOUR QUANTITIES NECESSARY FOR COORDINATE TRANSFORMATIONS
366 C WHICH DEPEND ON SUN POSITION (AND, HENCE, ON UNIVERSAL TIME AND SEASON)
367 C
368 C------- INPUT PARAMETERS:
369 C IYR,IDAY,IHOUR,MIN,ISEC - YEAR, DAY, AND UNIVERSAL TIME IN HOURS, MINUTES,
370 C AND SECONDS (IDAY=1 CORRESPONDS TO JANUARY 1).
371 C
372 C------- OUTPUT PARAMETERS:
373 C GST - GREENWICH MEAN SIDEREAL TIME, SLONG - LONGITUDE ALONG ECLIPTIC
374 C SRASN - RIGHT ASCENSION, SDEC - DECLINATION OF THE SUN (RADIANS)
375 C ORIGINAL VERSION OF THIS SUBROUTINE HAS BEEN COMPILED FROM:
376 C RUSSELL, C.T., COSMIC ELECTRODYNAMICS, 1971, V.2, PP.184-196.
377 C
378 C LAST MODIFICATION: MARCH 31, 2003 (ONLY SOME NOTATION CHANGES)
379 C
380 C ORIGINAL VERSION WRITTEN BY: Gilbert D. Mead
381 C
382 DOUBLE PRECISION DJ,FDAY
383 DATA RAD/57.295779513/
384 C
385 IF(IYEAR.LT.1901.OR.IYEAR.GT.2099) RETURN
386 FDAY=DFLOAT(IHOUR*3600+MIN*60+ISEC)/86400.D0
387 DJ=365*(IYEAR-1900)+(IYEAR-1901)/4+IDAY-0.5D0+FDAY
388 T=DJ/36525.
389 VL=DMOD(279.696678+0.9856473354*DJ,360.D0)
390 GST=DMOD(279.690983+.9856473354*DJ+360.*FDAY+180.,360.D0)/RAD
391 G=DMOD(358.475845+0.985600267*DJ,360.D0)/RAD
392 SLONG=(VL+(1.91946-0.004789*T)*SIN(G)+0.020094*SIN(2.*G))/RAD
393 IF(SLONG.GT.6.2831853) SLONG=SLONG-6.2831853
394 IF (SLONG.LT.0.) SLONG=SLONG+6.2831853
395 OBLIQ=(23.45229-0.0130125*T)/RAD
396 SOB=SIN(OBLIQ)
397 SLP=SLONG-9.924E-5
398 C
399 C THE LAST CONSTANT IS A CORRECTION FOR THE ANGULAR ABERRATION DUE TO
400 C THE ORBITAL MOTION OF THE EARTH
401 C
402 SIND=SOB*SIN(SLP)
403 COSD=SQRT(1.-SIND**2)
404 SC=SIND/COSD
405 SDEC=ATAN(SC)
406 SRASN=3.141592654-ATAN2(COS(OBLIQ)/SOB*SC,-COS(SLP)/COSD)
407 RETURN
408 END
409 C
410 C================================================================================
411 c
412 SUBROUTINE SPHCAR (R,THETA,PHI,X,Y,Z,J)
413 C
414 C CONVERTS SPHERICAL COORDS INTO CARTESIAN ONES AND VICA VERSA
415 C (THETA AND PHI IN RADIANS).
416 C
417 C J>0 J<0
418 C-----INPUT: J,R,THETA,PHI J,X,Y,Z
419 C----OUTPUT: X,Y,Z R,THETA,PHI
420 C
421 C NOTE: AT THE POLES (X=0 AND Y=0) WE ASSUME PHI=0 (WHEN CONVERTING
422 C FROM CARTESIAN TO SPHERICAL COORDS, I.E., FOR J<0)
423 C
424 C LAST MOFIFICATION: APRIL 1, 2003 (ONLY SOME NOTATION CHANGES AND MORE
425 C COMMENTS ADDED)
426 C
427 C AUTHOR: N. A. TSYGANENKO
428 C
429 IF(J.GT.0) GOTO 3
430 SQ=X**2+Y**2
431 R=SQRT(SQ+Z**2)
432 IF (SQ.NE.0.) GOTO 2
433 PHI=0.
434 IF (Z.LT.0.) GOTO 1
435 THETA=0.
436 RETURN
437 1 THETA=3.141592654
438 RETURN
439 2 SQ=SQRT(SQ)
440 PHI=ATAN2(Y,X)
441 THETA=ATAN2(SQ,Z)
442 IF (PHI.LT.0.) PHI=PHI+6.28318531
443 RETURN
444 3 SQ=R*SIN(THETA)
445 X=SQ*COS(PHI)
446 Y=SQ*SIN(PHI)
447 Z=R*COS(THETA)
448 RETURN
449 END
450 C
451 C===========================================================================
452 c
453 SUBROUTINE BSPCAR (THETA,PHI,BR,BTHETA,BPHI,BX,BY,BZ)
454 C
455 C CALCULATES CARTESIAN FIELD COMPONENTS FROM SPHERICAL ONES
456 C-----INPUT: THETA,PHI - SPHERICAL ANGLES OF THE POINT IN RADIANS
457 C BR,BTHETA,BPHI - SPHERICAL COMPONENTS OF THE FIELD
458 C-----OUTPUT: BX,BY,BZ - CARTESIAN COMPONENTS OF THE FIELD
459 C
460 C LAST MOFIFICATION: APRIL 1, 2003 (ONLY SOME NOTATION CHANGES)
461 C
462 C WRITTEN BY: N. A. TSYGANENKO
463 C
464 S=SIN(THETA)
465 C=COS(THETA)
466 SF=SIN(PHI)
467 CF=COS(PHI)
468 BE=BR*S+BTHETA*C
469 BX=BE*CF-BPHI*SF
470 BY=BE*SF+BPHI*CF
471 BZ=BR*C-BTHETA*S
472 RETURN
473 END
474 c
475 C==============================================================================
476 C
477 SUBROUTINE BCARSP (X,Y,Z,BX,BY,BZ,BR,BTHETA,BPHI)
478 C
479 CALCULATES SPHERICAL FIELD COMPONENTS FROM THOSE IN CARTESIAN SYSTEM
480 C
481 C-----INPUT: X,Y,Z - CARTESIAN COMPONENTS OF THE POSITION VECTOR
482 C BX,BY,BZ - CARTESIAN COMPONENTS OF THE FIELD VECTOR
483 C-----OUTPUT: BR,BTHETA,BPHI - SPHERICAL COMPONENTS OF THE FIELD VECTOR
484 C
485 C NOTE: AT THE POLES (THETA=0 OR THETA=PI) WE ASSUME PHI=0,
486 C AND HENCE BTHETA=BX, BPHI=BY
487 C
488 C WRITTEN AND ADDED TO THIS PACKAGE: APRIL 1, 2003,
489 C AUTHOR: N. A. TSYGANENKO
490 C
491 RHO2=X**2+Y**2
492 R=SQRT(RHO2+Z**2)
493 RHO=SQRT(RHO2)
494
495 IF (RHO.NE.0.) THEN
496 CPHI=X/RHO
497 SPHI=Y/RHO
498 ELSE
499 CPHI=1.
500 SPHI=0.
501 ENDIF
502
503 CT=Z/R
504 ST=RHO/R
505
506 BR=(X*BX+Y*BY+Z*BZ)/R
507 BTHETA=(BX*CPHI+BY*SPHI)*CT-BZ*ST
508 BPHI=BY*CPHI-BX*SPHI
509
510 RETURN
511 END
512 C
513 c=====================================================================================
514 C
515 SUBROUTINE RECALC (IYEAR,IDAY,IHOUR,MIN,ISEC)
516 C
517 C 1. PREPARES ELEMENTS OF ROTATION MATRICES FOR TRANSFORMATIONS OF VECTORS BETWEEN
518 C SEVERAL COORDINATE SYSTEMS, MOST FREQUENTLY USED IN SPACE PHYSICS.
519 C
520 C 2. PREPARES COEFFICIENTS USED IN THE CALCULATION OF THE MAIN GEOMAGNETIC FIELD
521 C (IGRF MODEL)
522 C
523 C THIS SUBROUTINE SHOULD BE INVOKED BEFORE USING THE FOLLOWING SUBROUTINES:
524 C IGRF_GEO, IGRF_GSM, DIP, GEOMAG, GEOGSM, MAGSM, SMGSM, GSMGSE, GEIGEO.
525 C
526 C THERE IS NO NEED TO REPEATEDLY INVOKE RECALC, IF MULTIPLE CALCULATIONS ARE MADE
527 C FOR THE SAME DATE AND TIME.
528 C
529 C-----INPUT PARAMETERS:
530 C
531 C IYEAR - YEAR NUMBER (FOUR DIGITS)
532 C IDAY - DAY OF YEAR (DAY 1 = JAN 1)
533 C IHOUR - HOUR OF DAY (00 TO 23)
534 C MIN - MINUTE OF HOUR (00 TO 59)
535 C ISEC - SECONDS OF MINUTE (00 TO 59)
536 C
537 C-----OUTPUT PARAMETERS: NONE (ALL OUTPUT QUANTITIES ARE PLACED
538 C INTO THE COMMON BLOCKS /GEOPACK1/ AND /GEOPACK2/)
539 C
540 C OTHER SUBROUTINES CALLED BY THIS ONE: SUN
541 C
542 C AUTHOR: N.A. TSYGANENKO
543 C DATE: DEC.1, 1991
544 C
545 C REVISION OF MAY 3, 2005:
546 C The table of IGRF coefficients was extended to include those for the epoch 2005
547 c the maximal order of spherical harmonics was also increased up to 13
548 c (for details, see http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html)
549 c
550 C REVISION OF APRIL 3, 2003:
551 c The code now includes preparation of the model coefficients for the subroutines
552 c IGRF and GEOMAG. This eliminates the need for the SAVE statements, used in the
553 c old versions, making the codes easier and more compiler-independent.
554 C
555 COMMON /GEOPACK1/ ST0,CT0,SL0,CL0,CTCL,STCL,CTSL,STSL,SFI,CFI,SPS,
556 * CPS,SHI,CHI,HI,PSI,XMUT,A11,A21,A31,A12,A22,A32,A13,A23,A33,DS3,
557 * CGST,SGST,BA(6)
558 C
559 C THE COMMON BLOCK /GEOPACK1/ CONTAINS ELEMENTS OF THE ROTATION MATRICES AND OTHER
560 C PARAMETERS RELATED TO THE COORDINATE TRANSFORMATIONS PERFORMED BY THIS PACKAGE
561 C
562 COMMON /GEOPACK2/ G(105),H(105),REC(105)
563 C
564 C THE COMMON BLOCK /GEOPACK2/ CONTAINS COEFFICIENTS OF THE IGRF FIELD MODEL, CALCULATED
565 C FOR A GIVEN YEAR AND DAY FROM THEIR STANDARD EPOCH VALUES. THE ARRAY REC CONTAINS
566 C COEFFICIENTS USED IN THE RECURSION RELATIONS FOR LEGENDRE ASSOCIATE POLYNOMIALS.
567 C
568 DIMENSION G65(105),H65(105),G70(105),H70(105),G75(105),H75(105),
569 + G80(105),H80(105),G85(105),H85(105),G90(105),H90(105),G95(105),
570 + H95(105),G00(105),H00(105),G05(105),H05(105),DG05(45),DH05(45)
571 c
572 DATA G65/0.,-30334.,-2119.,-1662.,2997.,1594.,1297.,-2038.,1292.,
573 *856.,957.,804.,479.,-390.,252.,-219.,358.,254.,-31.,-157.,-62.,
574 *45.,61.,8.,-228.,4.,1.,-111.,75.,-57.,4.,13.,-26.,-6.,13.,1.,13.,
575 *5.,-4.,-14.,0.,8.,-1.,11.,4.,8.,10.,2.,-13.,10.,-1.,-1.,5.,1.,-2.,
576 *-2.,-3.,2.,-5.,-2.,4.,4.,0.,2.,2.,0.,39*0./
577 DATA H65/0.,0.,5776.,0.,-2016.,114.,0.,-404.,240.,-165.,0.,148.,
578 *-269.,13.,-269.,0.,19.,128.,-126.,-97.,81.,0.,-11.,100.,68.,-32.,
579 *-8.,-7.,0.,-61.,-27.,-2.,6.,26.,-23.,-12.,0.,7.,-12.,9.,-16.,4.,
580 *24.,-3.,-17.,0.,-22.,15.,7.,-4.,-5.,10.,10.,-4.,1.,0.,2.,1.,2.,
581 *6.,-4.,0.,-2.,3.,0.,-6.,39*0./
582 c
583 DATA G70/0.,-30220.,-2068.,-1781.,3000.,1611.,1287.,-2091.,1278.,
584 *838.,952.,800.,461.,-395.,234.,-216.,359.,262.,-42.,-160.,-56.,
585 *43.,64.,15.,-212.,2.,3.,-112.,72.,-57.,1.,14.,-22.,-2.,13.,-2.,
586 *14.,6.,-2.,-13.,-3.,5.,0.,11.,3.,8.,10.,2.,-12.,10.,-1.,0.,3.,
587 *1.,-1.,-3.,-3.,2.,-5.,-1.,6.,4.,1.,0.,3.,-1.,39*0./
588 DATA H70/0.,0.,5737.,0.,-2047.,25.,0.,-366.,251.,-196.,0.,167.,
589 *-266.,26.,-279.,0.,26.,139.,-139.,-91.,83.,0.,-12.,100.,72.,-37.,
590 *-6.,1.,0.,-70.,-27.,-4.,8.,23.,-23.,-11.,0.,7.,-15.,6.,-17.,6.,
591 *21.,-6.,-16.,0.,-21.,16.,6.,-4.,-5.,10.,11.,-2.,1.,0.,1.,1.,3.,
592 *4.,-4.,0.,-1.,3.,1.,-4.,39*0./
593 c
594 DATA G75/0.,-30100.,-2013.,-1902.,3010.,1632.,1276.,-2144.,1260.,
595 *830.,946.,791.,438.,-405.,216.,-218.,356.,264.,-59.,-159.,-49.,
596 *45.,66.,28.,-198.,1.,6.,-111.,71.,-56.,1.,16.,-14.,0.,12.,-5.,
597 *14.,6.,-1.,-12.,-8.,4.,0.,10.,1.,7.,10.,2.,-12.,10.,-1.,-1.,4.,
598 *1.,-2.,-3.,-3.,2.,-5.,-2.,5.,4.,1.,0.,3.,-1.,39*0./
599 DATA H75/0.,0.,5675.,0.,-2067.,-68.,0.,-333.,262.,-223.,0.,191.,
600 *-265.,39.,-288.,0.,31.,148.,-152.,-83.,88.,0.,-13.,99.,75.,-41.,
601 *-4.,11.,0.,-77.,-26.,-5.,10.,22.,-23.,-12.,0.,6.,-16.,4.,-19.,6.,
602 *18.,-10.,-17.,0.,-21.,16.,7.,-4.,-5.,10.,11.,-3.,1.,0.,1.,1.,3.,
603 *4.,-4.,-1.,-1.,3.,1.,-5.,39*0./
604 c
605 DATA G80/0.,-29992.,-1956.,-1997.,3027.,1663.,1281.,-2180.,1251.,
606 *833.,938.,782.,398.,-419.,199.,-218.,357.,261.,-74.,-162.,-48.,
607 *48.,66.,42.,-192.,4.,14.,-108.,72.,-59.,2.,21.,-12.,1.,11.,-2.,
608 *18.,6.,0.,-11.,-7.,4.,3.,6.,-1.,5.,10.,1.,-12.,9.,-3.,-1.,7.,2.,
609 *-5.,-4.,-4.,2.,-5.,-2.,5.,3.,1.,2.,3.,0.,39*0./
610 DATA H80/0.,0.,5604.,0.,-2129.,-200.,0.,-336.,271.,-252.,0.,212.,
611 *-257.,53.,-297.,0.,46.,150.,-151.,-78.,92.,0.,-15.,93.,71.,-43.,
612 *-2.,17.,0.,-82.,-27.,-5.,16.,18.,-23.,-10.,0.,7.,-18.,4.,-22.,9.,
613 *16.,-13.,-15.,0.,-21.,16.,9.,-5.,-6.,9.,10.,-6.,2.,0.,1.,0.,3.,
614 *6.,-4.,0.,-1.,4.,0.,-6.,39*0./
615 c
616 DATA G85/0.,-29873.,-1905.,-2072.,3044.,1687.,1296.,-2208.,1247.,
617 *829.,936.,780.,361.,-424.,170.,-214.,355.,253.,-93.,-164.,-46.,
618 *53.,65.,51.,-185.,4.,16.,-102.,74.,-62.,3.,24.,-6.,4.,10.,0.,21.,
619 *6.,0.,-11.,-9.,4.,4.,4.,-4.,5.,10.,1.,-12.,9.,-3.,-1.,7.,1.,-5.,
620 *-4.,-4.,3.,-5.,-2.,5.,3.,1.,2.,3.,0.,39*0./
621 DATA H85/0.,0.,5500.,0.,-2197.,-306.,0.,-310.,284.,-297.,0.,232.,
622 *-249.,69.,-297.,0.,47.,150.,-154.,-75.,95.,0.,-16.,88.,69.,-48.,
623 *-1.,21.,0.,-83.,-27.,-2.,20.,17.,-23.,-7.,0.,8.,-19.,5.,-23.,11.,
624 *14.,-15.,-11.,0.,-21.,15.,9.,-6.,-6.,9.,9.,-7.,2.,0.,1.,0.,3.,
625 *6.,-4.,0.,-1.,4.,0.,-6.,39*0./
626 c
627 DATA G90/0., -29775., -1848., -2131., 3059., 1686., 1314.,
628 * -2239., 1248., 802., 939., 780., 325., -423.,
629 * 141., -214., 353., 245., -109., -165., -36.,
630 * 61., 65., 59., -178., 3., 18., -96.,
631 * 77., -64., 2., 26., -1., 5., 9.,
632 * 0., 23., 5., -1., -10., -12., 3.,
633 * 4., 2., -6., 4., 9., 1., -12.,
634 * 9., -4., -2., 7., 1., -6., -3.,
635 * -4., 2., -5., -2., 4., 3., 1.,
636 * 3., 3., 0., 39*0./
637
638 DATA H90/0., 0., 5406., 0., -2279., -373., 0.,
639 * -284., 293., -352., 0., 247., -240., 84.,
640 * -299., 0., 46., 154., -153., -69., 97.,
641 * 0., -16., 82., 69., -52., 1., 24.,
642 * 0., -80., -26., 0., 21., 17., -23.,
643 * -4., 0., 10., -19., 6., -22., 12.,
644 * 12., -16., -10., 0., -20., 15., 11.,
645 * -7., -7., 9., 8., -7., 2., 0.,
646 * 2., 1., 3., 6., -4., 0., -2.,
647 * 3., -1., -6., 39*0./
648
649 DATA G95/0., -29692., -1784., -2200., 3070., 1681., 1335.,
650 * -2267., 1249., 759., 940., 780., 290., -418.,
651 * 122., -214., 352., 235., -118., -166., -17.,
652 * 68., 67., 68., -170., -1., 19., -93.,
653 * 77., -72., 1., 28., 5., 4., 8.,
654 * -2., 25., 6., -6., -9., -14., 9.,
655 * 6., -5., -7., 4., 9., 3., -10.,
656 * 8., -8., -1., 10., -2., -8., -3.,
657 * -6., 2., -4., -1., 4., 2., 2.,
658 * 5., 1., 0., 39*0./
659
660 DATA H95/0., 0., 5306., 0., -2366., -413., 0.,
661 * -262., 302., -427., 0., 262., -236., 97.,
662 * -306., 0., 46., 165., -143., -55., 107.,
663 * 0., -17., 72., 67., -58., 1., 36.,
664 * 0., -69., -25., 4., 24., 17., -24.,
665 * -6., 0., 11., -21., 8., -23., 15.,
666 * 11., -16., -4., 0., -20., 15., 12.,
667 * -6., -8., 8., 5., -8., 3., 0.,
668 * 1., 0., 4., 5., -5., -1., -2.,
669 * 1., -2., -7., 39*0./
670
671 DATA G00/0.,-29619.4, -1728.2, -2267.7, 3068.4, 1670.9, 1339.6,
672 * -2288., 1252.1, 714.5, 932.3, 786.8, 250., -403.,
673 * 111.3, -218.8, 351.4, 222.3, -130.4, -168.6, -12.9,
674 * 72.3, 68.2, 74.2, -160.9, -5.9, 16.9, -90.4,
675 * 79.0, -74.0, 0., 33.3, 9.1, 6.9, 7.3,
676 * -1.2, 24.4, 6.6, -9.2, -7.9, -16.6, 9.1,
677 * 7.0, -7.9, -7., 5., 9.4, 3., - 8.4,
678 * 6.3, -8.9, -1.5, 9.3, -4.3, -8.2, -2.6,
679 * -6., 1.7, -3.1, -0.5, 3.7, 1., 2.,
680 * 4.2, 0.3, -1.1, 2.7, -1.7, -1.9, 1.5,
681 * -0.1, 0.1, -0.7, 0.7, 1.7, 0.1, 1.2,
682 * 4.0, -2.2, -0.3, 0.2, 0.9, -0.2, 0.9,
683 * -0.5, 0.3, -0.3, -0.4, -0.1, -0.2, -0.4,
684 * -0.2, -0.9, 0.3, 0.1, -0.4, 1.3, -0.4,
685 * 0.7, -0.4, 0.3, -0.1, 0.4, 0., 0.1/
686
687
688 DATA H00/0., 0., 5186.1, 0., -2481.6, -458.0, 0.,
689 * -227.6, 293.4, -491.1, 0., 272.6, -231.9, 119.8,
690 * -303.8, 0., 43.8, 171.9, -133.1, -39.3, 106.3,
691 * 0., -17.4, 63.7, 65.1, -61.2, 0.7, 43.8,
692 * 0., -64.6, -24.2, 6.2, 24., 14.8, -25.4,
693 * -5.8, 0.0, 11.9, -21.5, 8.5, -21.5, 15.5,
694 * 8.9, -14.9, -2.1, 0.0, -19.7, 13.4, 12.5,
695 * -6.2, -8.4, 8.4, 3.8, -8.2, 4.8, 0.0,
696 * 1.7, 0.0, 4.0, 4.9, -5.9, -1.2, -2.9,
697 * 0.2, -2.2, -7.4, 0.0, 0.1, 1.3, -0.9,
698 * -2.6, 0.9, -0.7, -2.8, -0.9, -1.2, -1.9,
699 * -0.9, 0.0, -0.4, 0.3, 2.5, -2.6, 0.7,
700 * 0.3, 0.0, 0.0, 0.3, -0.9, -0.4, 0.8,
701 * 0.0, -0.9, 0.2, 1.8, -0.4, -1.0, -0.1,
702 * 0.7, 0.3, 0.6, 0.3, -0.2, -0.5, -0.9/
703 *
704
705 DATA G05/0.,-29556.8, -1671.8, -2340.5, 3047., 1656.9, 1335.7,
706 * -2305.3, 1246.8, 674.4, 919.8, 798.2, 211.5, -379.5,
707 * 100.2, -227.6, 354.4, 208.8, -136.6, -168.3, -14.1,
708 * 72.9, 69.6, 76.6, -151.1, -15.0, 14.7, -86.4,
709 * 79.8, -74.4, -1.4, 38.6, 12.3, 9.4, 5.5,
710 * 2.0, 24.8, 7.7, -11.4, -6.8, -18.0, 10.0,
711 * 9.4, -11.4, -5.0, 5.6, 9.8, 3.6, -7.0,
712 * 5.0, -10.8, -1.3, 8.7, -6.7, -9.2, -2.2,
713 * -6.3, 1.6, -2.5, -0.1, 3.0, 0.3, 2.1,
714 * 3.9, -0.1, -2.2, 2.9, -1.6, -1.7, 1.5,
715 * -0.2, 0.2, -0.7, 0.5, 1.8, 0.1, 1.0,
716 * 4.1, -2.2, -0.3, 0.3, 0.9, -0.4, 1.0,
717 * -0.4, 0.5, -0.3, -0.4, 0.0, -0.4, 0.0,
718 * -0.2, -0.9, 0.3, 0.3, -0.4, 1.2, -0.4,
719 * 0.7, -0.3, 0.4, -0.1, 0.4, -0.1, -0.3/
720
721 DATA H05/0., 0.0, 5080.0, 0.0, -2594.9, -516.7, 0.0,
722 * -200.4, 269.3, -524.5, 0.0, 281.4, -225.8, 145.7,
723 * -304.7, 0.0, 42.7, 179.8, -123.0, -19.5, 103.6,
724 * 0.0, -20.2, 54.7, 63.7, -63.4, 0.0, 50.3,
725 * 0.0, -61.4, -22.5, 6.9, 25.4, 10.9, -26.4,
726 * -4.8, 0.0, 11.2, -21.0, 9.7, -19.8, 16.1,
727 * 7.7, -12.8, -0.1, 0.0, -20.1, 12.9, 12.7,
728 * -6.7, -8.1, 8.1, 2.9, -7.9, 5.9, 0.0,
729 * 2.4, 0.2, 4.4, 4.7, -6.5, -1.0, -3.4,
730 * -0.9, -2.3, -8.0, 0.0, 0.3, 1.4, -0.7,
731 * -2.4, 0.9, -0.6, -2.7, -1.0, -1.5, -2.0,
732 * -1.4, 0.0, -0.5, 0.3, 2.3, -2.7, 0.6,
733 * 0.4, 0.0, 0.0, 0.3, -0.8, -0.4, 1.0,
734 * 0.0, -0.7, 0.3, 1.7, -0.5, -1.0, 0.0,
735 * 0.7, 0.2, 0.6, 0.4, -0.2, -0.5, -1.0/
736
737 DATA DG05/0.0, 8.8, 10.8, -15.0, -6.9, -1.0, -0.3,
738 * -3.1, -0.9, -6.8, -2.5, 2.8, -7.1, 5.9,
739 * -3.2, -2.6, 0.4, -3.0, -1.2, 0.2, -0.6,
740 * -0.8, 0.2, -0.2, 2.1, -2.1, -0.4, 1.3,
741 * -0.4, 0.0, -0.2, 1.1, 0.6, 0.4, -0.5,
742 * 0.9, -0.2, 0.2, -0.2, 0.2, -0.2, 0.2,
743 * 0.5, -0.7, 0.5/
744
745 DATA DH05/0.0, 0.0, -21.3, 0.0, -23.3, -14.0, 0.0,
746 * 5.4, -6.5, -2.0, 0.0, 2.0, 1.8, 5.6,
747 * 0.0, 0.0, 0.1, 1.8, 2.0, 4.5, -1.0,
748 * 0.0, -0.4, -1.9, -0.4, -0.4, -0.2, 0.9,
749 * 0.0, 0.8, 0.4, 0.1, 0.2, -0.9, -0.3,
750 * 0.3, 0.0, -0.2, 0.2, 0.2, 0.4, 0.2,
751 * -0.3, 0.5, 0.4/
752 C
753 C
754 IY=IYEAR
755 C
756 C WE ARE RESTRICTED BY THE INTERVAL 1965-2010, FOR WHICH THE IGRF COEFFICIENTS
757 c ARE KNOWN; IF IYEAR IS OUTSIDE THIS INTERVAL, THEN THE SUBROUTINE USES THE
758 C NEAREST LIMITING VALUE AND PRINTS A WARNING:
759 C
760 IF(IY.LT.1965) THEN
761 IY=1965
762 WRITE (*,10) IYEAR,IY
763 ENDIF
764
765 IF(IY.GT.2010) THEN
766 IY=2010
767 WRITE (*,10) IYEAR,IY
768 ENDIF
769
770 C
771 C CALCULATE THE ARRAY REC, CONTAINING COEFFICIENTS FOR THE RECURSION RELATIONS,
772 C USED IN THE IGRF SUBROUTINE FOR CALCULATING THE ASSOCIATE LEGENDRE POLYNOMIALS
773 C AND THEIR DERIVATIVES:
774 c
775 DO 20 N=1,14
776 N2=2*N-1
777 N2=N2*(N2-2)
778 DO 20 M=1,N
779 MN=N*(N-1)/2+M
780 20 REC(MN)=FLOAT((N-M)*(N+M-2))/FLOAT(N2)
781 C
782 IF (IY.LT.1970) GOTO 50 !INTERPOLATE BETWEEN 1965 - 1970
783 IF (IY.LT.1975) GOTO 60 !INTERPOLATE BETWEEN 1970 - 1975
784 IF (IY.LT.1980) GOTO 70 !INTERPOLATE BETWEEN 1975 - 1980
785 IF (IY.LT.1985) GOTO 80 !INTERPOLATE BETWEEN 1980 - 1985
786 IF (IY.LT.1990) GOTO 90 !INTERPOLATE BETWEEN 1985 - 1990
787 IF (IY.LT.1995) GOTO 100 !INTERPOLATE BETWEEN 1990 - 1995
788 IF (IY.LT.2000) GOTO 110 !INTERPOLATE BETWEEN 1995 - 2000
789 IF (IY.LT.2005) GOTO 120 !INTERPOLATE BETWEEN 2000 - 2005
790 C
791 C EXTRAPOLATE BEYOND 2005:
792 C
793
794 DT=FLOAT(IY)+FLOAT(IDAY-1)/365.25-2005.
795 DO 40 N=1,105
796 G(N)=G05(N)
797 H(N)=H05(N)
798 IF (N.GT.45) GOTO 40
799 G(N)=G(N)+DG05(N)*DT
800 H(N)=H(N)+DH05(N)*DT
801 40 CONTINUE
802 GOTO 300
803 C
804 C INTERPOLATE BETWEEEN 1965 - 1970:
805 C
806 50 F2=(FLOAT(IY)+FLOAT(IDAY-1)/365.25-1965)/5.
807 F1=1.-F2
808 DO 55 N=1,66
809 G(N)=G65(N)*F1+G70(N)*F2
810 55 H(N)=H65(N)*F1+H70(N)*F2
811 GOTO 300
812 C
813 C INTERPOLATE BETWEEN 1970 - 1975:
814 C
815 60 F2=(FLOAT(IY)+FLOAT(IDAY-1)/365.25-1970)/5.
816 F1=1.-F2
817 DO 65 N=1,66
818 G(N)=G70(N)*F1+G75(N)*F2
819 65 H(N)=H70(N)*F1+H75(N)*F2
820 GOTO 300
821 C
822 C INTERPOLATE BETWEEN 1975 - 1980:
823 C
824 70 F2=(FLOAT(IY)+FLOAT(IDAY-1)/365.25-1975)/5.
825 F1=1.-F2
826 DO 75 N=1,66
827 G(N)=G75(N)*F1+G80(N)*F2
828 75 H(N)=H75(N)*F1+H80(N)*F2
829 GOTO 300
830 C
831 C INTERPOLATE BETWEEN 1980 - 1985:
832 C
833 80 F2=(FLOAT(IY)+FLOAT(IDAY-1)/365.25-1980)/5.
834 F1=1.-F2
835 DO 85 N=1,66
836 G(N)=G80(N)*F1+G85(N)*F2
837 85 H(N)=H80(N)*F1+H85(N)*F2
838 GOTO 300
839 C
840 C INTERPOLATE BETWEEN 1985 - 1990:
841 C
842 90 F2=(FLOAT(IY)+FLOAT(IDAY-1)/365.25-1985)/5.
843 F1=1.-F2
844 DO 95 N=1,66
845 G(N)=G85(N)*F1+G90(N)*F2
846 95 H(N)=H85(N)*F1+H90(N)*F2
847 GOTO 300
848 C
849 C INTERPOLATE BETWEEN 1990 - 1995:
850 C
851 100 F2=(FLOAT(IY)+FLOAT(IDAY-1)/365.25-1990)/5.
852 F1=1.-F2
853 DO 105 N=1,66
854 G(N)=G90(N)*F1+G95(N)*F2
855 105 H(N)=H90(N)*F1+H95(N)*F2
856 GOTO 300
857 C
858 C INTERPOLATE BETWEEN 1995 - 2000:
859 C
860 110 F2=(FLOAT(IY)+FLOAT(IDAY-1)/365.25-1995)/5.
861 F1=1.-F2
862 DO 115 N=1,105 ! THE 2000 COEFFICIENTS (G00) GO THROUGH THE ORDER 13, NOT 10
863 G(N)=G95(N)*F1+G00(N)*F2
864 115 H(N)=H95(N)*F1+H00(N)*F2
865 GOTO 300
866 C
867 C INTERPOLATE BETWEEN 2000 - 2005:
868 C
869 120 F2=(FLOAT(IY)+FLOAT(IDAY-1)/365.25-2000)/5.
870 F1=1.-F2
871 DO 125 N=1,105
872 G(N)=G00(N)*F1+G05(N)*F2
873 125 H(N)=H00(N)*F1+H05(N)*F2
874 GOTO 300
875 C
876 C COEFFICIENTS FOR A GIVEN YEAR HAVE BEEN CALCULATED; NOW MULTIPLY
877 C THEM BY SCHMIDT NORMALIZATION FACTORS:
878 C
879 300 S=1.
880 DO 130 N=2,14
881 MN=N*(N-1)/2+1
882 S=S*FLOAT(2*N-3)/FLOAT(N-1)
883 G(MN)=G(MN)*S
884 H(MN)=H(MN)*S
885 P=S
886 DO 130 M=2,N
887 AA=1.
888 IF (M.EQ.2) AA=2.
889 P=P*SQRT(AA*FLOAT(N-M+1)/FLOAT(N+M-2))
890 MNN=MN+M-1
891 G(MNN)=G(MNN)*P
892 130 H(MNN)=H(MNN)*P
893
894 G10=-G(2)
895 G11= G(3)
896 H11= H(3)
897 C
898 C NOW CALCULATE THE COMPONENTS OF THE UNIT VECTOR EzMAG IN GEO COORD.SYSTEM:
899 C SIN(TETA0)*COS(LAMBDA0), SIN(TETA0)*SIN(LAMBDA0), AND COS(TETA0)
900 C ST0 * CL0 ST0 * SL0 CT0
901 C
902 SQ=G11**2+H11**2
903 SQQ=SQRT(SQ)
904 SQR=SQRT(G10**2+SQ)
905 SL0=-H11/SQQ
906 CL0=-G11/SQQ
907 ST0=SQQ/SQR
908 CT0=G10/SQR
909 STCL=ST0*CL0
910 STSL=ST0*SL0
911 CTSL=CT0*SL0
912 CTCL=CT0*CL0
913 C
914 CALL SUN (IY,IDAY,IHOUR,MIN,ISEC,GST,SLONG,SRASN,SDEC)
915 C
916 C S1,S2, AND S3 ARE THE COMPONENTS OF THE UNIT VECTOR EXGSM=EXGSE IN THE
917 C SYSTEM GEI POINTING FROM THE EARTH'S CENTER TO THE SUN:
918 C
919 S1=COS(SRASN)*COS(SDEC)
920 S2=SIN(SRASN)*COS(SDEC)
921 S3=SIN(SDEC)
922 CGST=COS(GST)
923 SGST=SIN(GST)
924 C
925 C DIP1, DIP2, AND DIP3 ARE THE COMPONENTS OF THE UNIT VECTOR EZSM=EZMAG
926 C IN THE SYSTEM GEI:
927 C
928 DIP1=STCL*CGST-STSL*SGST
929 DIP2=STCL*SGST+STSL*CGST
930 DIP3=CT0
931 C
932 C NOW CALCULATE THE COMPONENTS OF THE UNIT VECTOR EYGSM IN THE SYSTEM GEI
933 C BY TAKING THE VECTOR PRODUCT D x S AND NORMALIZING IT TO UNIT LENGTH:
934 C
935 Y1=DIP2*S3-DIP3*S2
936 Y2=DIP3*S1-DIP1*S3
937 Y3=DIP1*S2-DIP2*S1
938 Y=SQRT(Y1*Y1+Y2*Y2+Y3*Y3)
939 Y1=Y1/Y
940 Y2=Y2/Y
941 Y3=Y3/Y
942 C
943 C THEN IN THE GEI SYSTEM THE UNIT VECTOR Z = EZGSM = EXGSM x EYGSM = S x Y
944 C HAS THE COMPONENTS:
945 C
946 Z1=S2*Y3-S3*Y2
947 Z2=S3*Y1-S1*Y3
948 Z3=S1*Y2-S2*Y1
949 C
950 C THE VECTOR EZGSE (HERE DZ) IN GEI HAS THE COMPONENTS (0,-SIN(DELTA),
951 C COS(DELTA)) = (0.,-0.397823,0.917462); HERE DELTA = 23.44214 DEG FOR
952 C THE EPOCH 1978 (SEE THE BOOK BY GUREVICH OR OTHER ASTRONOMICAL HANDBOOKS).
953 C HERE THE MOST ACCURATE TIME-DEPENDENT FORMULA IS USED:
954 C
955 DJ=FLOAT(365*(IY-1900)+(IY-1901)/4 +IDAY)
956 * -0.5+FLOAT(IHOUR*3600+MIN*60+ISEC)/86400.
957 T=DJ/36525.
958 OBLIQ=(23.45229-0.0130125*T)/57.2957795
959 DZ1=0.
960 DZ2=-SIN(OBLIQ)
961 DZ3=COS(OBLIQ)
962 C
963 C THEN THE UNIT VECTOR EYGSE IN GEI SYSTEM IS THE VECTOR PRODUCT DZ x S :
964 C
965 DY1=DZ2*S3-DZ3*S2
966 DY2=DZ3*S1-DZ1*S3
967 DY3=DZ1*S2-DZ2*S1
968 C
969 C THE ELEMENTS OF THE MATRIX GSE TO GSM ARE THE SCALAR PRODUCTS:
970 C CHI=EM22=(EYGSM,EYGSE), SHI=EM23=(EYGSM,EZGSE), EM32=(EZGSM,EYGSE)=-EM23,
971 C AND EM33=(EZGSM,EZGSE)=EM22
972 C
973 CHI=Y1*DY1+Y2*DY2+Y3*DY3
974 SHI=Y1*DZ1+Y2*DZ2+Y3*DZ3
975 HI=ASIN(SHI)
976 C
977 C TILT ANGLE: PSI=ARCSIN(DIP,EXGSM)
978 C
979 SPS=DIP1*S1+DIP2*S2+DIP3*S3
980 CPS=SQRT(1.-SPS**2)
981 PSI=ASIN(SPS)
982 C
983 C THE ELEMENTS OF THE MATRIX MAG TO SM ARE THE SCALAR PRODUCTS:
984 C CFI=GM22=(EYSM,EYMAG), SFI=GM23=(EYSM,EXMAG); THEY CAN BE DERIVED AS FOLLOWS:
985 C
986 C IN GEO THE VECTORS EXMAG AND EYMAG HAVE THE COMPONENTS (CT0*CL0,CT0*SL0,-ST0)
987 C AND (-SL0,CL0,0), RESPECTIVELY. HENCE, IN GEI THE COMPONENTS ARE:
988 C EXMAG: CT0*CL0*COS(GST)-CT0*SL0*SIN(GST)
989 C CT0*CL0*SIN(GST)+CT0*SL0*COS(GST)
990 C -ST0
991 C EYMAG: -SL0*COS(GST)-CL0*SIN(GST)
992 C -SL0*SIN(GST)+CL0*COS(GST)
993 C 0
994 C THE COMPONENTS OF EYSM IN GEI WERE FOUND ABOVE AS Y1, Y2, AND Y3;
995 C NOW WE ONLY HAVE TO COMBINE THE QUANTITIES INTO SCALAR PRODUCTS:
996 C
997 EXMAGX=CT0*(CL0*CGST-SL0*SGST)
998 EXMAGY=CT0*(CL0*SGST+SL0*CGST)
999 EXMAGZ=-ST0
1000 EYMAGX=-(SL0*CGST+CL0*SGST)
1001 EYMAGY=-(SL0*SGST-CL0*CGST)
1002 CFI=Y1*EYMAGX+Y2*EYMAGY
1003 SFI=Y1*EXMAGX+Y2*EXMAGY+Y3*EXMAGZ
1004 C
1005 XMUT=(ATAN2(SFI,CFI)+3.1415926536)*3.8197186342
1006 C
1007 C THE ELEMENTS OF THE MATRIX GEO TO GSM ARE THE SCALAR PRODUCTS:
1008 C
1009 C A11=(EXGEO,EXGSM), A12=(EYGEO,EXGSM), A13=(EZGEO,EXGSM),
1010 C A21=(EXGEO,EYGSM), A22=(EYGEO,EYGSM), A23=(EZGEO,EYGSM),
1011 C A31=(EXGEO,EZGSM), A32=(EYGEO,EZGSM), A33=(EZGEO,EZGSM),
1012 C
1013 C ALL THE UNIT VECTORS IN BRACKETS ARE ALREADY DEFINED IN GEI:
1014 C
1015 C EXGEO=(CGST,SGST,0), EYGEO=(-SGST,CGST,0), EZGEO=(0,0,1)
1016 C EXGSM=(S1,S2,S3), EYGSM=(Y1,Y2,Y3), EZGSM=(Z1,Z2,Z3)
1017 C AND THEREFORE:
1018 C
1019 A11=S1*CGST+S2*SGST
1020 A12=-S1*SGST+S2*CGST
1021 A13=S3
1022 A21=Y1*CGST+Y2*SGST
1023 A22=-Y1*SGST+Y2*CGST
1024 A23=Y3
1025 A31=Z1*CGST+Z2*SGST
1026 A32=-Z1*SGST+Z2*CGST
1027 A33=Z3
1028 C
1029 10 FORMAT(//1X,
1030 *'**** RECALC WARNS: YEAR IS OUT OF INTERVAL 1965-2010: IYEAR=',I4,
1031 * /,6X,'CALCULATIONS WILL BE DONE FOR IYEAR=',I4,/)
1032 RETURN
1033 END
1034 c
1035 c====================================================================
1036 C
1037 SUBROUTINE GEOMAG (XGEO,YGEO,ZGEO,XMAG,YMAG,ZMAG,J)
1038 C
1039 C CONVERTS GEOGRAPHIC (GEO) TO DIPOLE (MAG) COORDINATES OR VICA VERSA.
1040 C
1041 C J>0 J<0
1042 C-----INPUT: J,XGEO,YGEO,ZGEO J,XMAG,YMAG,ZMAG
1043 C-----OUTPUT: XMAG,YMAG,ZMAG XGEO,YGEO,ZGEO
1044 C
1045 C
1046 C ATTENTION: SUBROUTINE RECALC MUST BE INVOKED BEFORE GEOMAG IN TWO CASES:
1047 C /A/ BEFORE THE FIRST TRANSFORMATION OF COORDINATES
1048 C /B/ IF THE VALUES OF IYEAR AND/OR IDAY HAVE BEEN CHANGED
1049 C
1050 C
1051 C LAST MOFIFICATION: MARCH 30, 2003 (INVOCATION OF RECALC INSIDE THIS S/R WAS REMOVED)
1052 C
1053 C AUTHOR: N. A. TSYGANENKO
1054 C
1055 COMMON /GEOPACK1/ ST0,CT0,SL0,CL0,CTCL,STCL,CTSL,STSL,AB(19),BB(8)
1056
1057 IF(J.GT.0) THEN
1058 XMAG=XGEO*CTCL+YGEO*CTSL-ZGEO*ST0
1059 YMAG=YGEO*CL0-XGEO*SL0
1060 ZMAG=XGEO*STCL+YGEO*STSL+ZGEO*CT0
1061 ELSE
1062 XGEO=XMAG*CTCL-YMAG*SL0+ZMAG*STCL
1063 YGEO=XMAG*CTSL+YMAG*CL0+ZMAG*STSL
1064 ZGEO=ZMAG*CT0-XMAG*ST0
1065 ENDIF
1066
1067 RETURN
1068 END
1069 c
1070 c=========================================================================================
1071 c
1072 SUBROUTINE GEIGEO (XGEI,YGEI,ZGEI,XGEO,YGEO,ZGEO,J)
1073 C
1074 C CONVERTS EQUATORIAL INERTIAL (GEI) TO GEOGRAPHICAL (GEO) COORDS
1075 C OR VICA VERSA.
1076 C J>0 J<0
1077 C----INPUT: J,XGEI,YGEI,ZGEI J,XGEO,YGEO,ZGEO
1078 C----OUTPUT: XGEO,YGEO,ZGEO XGEI,YGEI,ZGEI
1079 C
1080 C ATTENTION: SUBROUTINE RECALC MUST BE INVOKED BEFORE GEIGEO IN TWO CASES:
1081 C /A/ BEFORE THE FIRST TRANSFORMATION OF COORDINATES
1082 C /B/ IF THE CURRENT VALUES OF IYEAR,IDAY,IHOUR,MIN,ISEC HAVE BEEN CHANGED
1083 C
1084 C LAST MODIFICATION: MARCH 31, 2003
1085
1086 C AUTHOR: N. A. TSYGANENKO
1087 C
1088 COMMON /GEOPACK1/ A(27),CGST,SGST,B(6)
1089 C
1090 IF(J.GT.0) THEN
1091 XGEO=XGEI*CGST+YGEI*SGST
1092 YGEO=YGEI*CGST-XGEI*SGST
1093 ZGEO=ZGEI
1094 ELSE
1095 XGEI=XGEO*CGST-YGEO*SGST
1096 YGEI=YGEO*CGST+XGEO*SGST
1097 ZGEI=ZGEO
1098 ENDIF
1099
1100 RETURN
1101 END
1102 C
1103 C=======================================================================================
1104 C
1105 SUBROUTINE MAGSM (XMAG,YMAG,ZMAG,XSM,YSM,ZSM,J)
1106 C
1107 C CONVERTS DIPOLE (MAG) TO SOLAR MAGNETIC (SM) COORDINATES OR VICA VERSA
1108 C
1109 C J>0 J<0
1110 C----INPUT: J,XMAG,YMAG,ZMAG J,XSM,YSM,ZSM
1111 C----OUTPUT: XSM,YSM,ZSM XMAG,YMAG,ZMAG
1112 C
1113 C ATTENTION: SUBROUTINE RECALC MUST BE INVOKED BEFORE MAGSM IN TWO CASES:
1114 C /A/ BEFORE THE FIRST TRANSFORMATION OF COORDINATES
1115 C /B/ IF THE VALUES OF IYEAR,IDAY,IHOUR,MIN,ISEC HAVE BEEN CHANGED
1116 C
1117 C LAST MODIFICATION: MARCH 31, 2003
1118 C
1119 C AUTHOR: N. A. TSYGANENKO
1120 C
1121 COMMON /GEOPACK1/ A(8),SFI,CFI,B(7),AB(10),BA(8)
1122 C
1123 IF (J.GT.0) THEN
1124 XSM=XMAG*CFI-YMAG*SFI
1125 YSM=XMAG*SFI+YMAG*CFI
1126 ZSM=ZMAG
1127 ELSE
1128 XMAG=XSM*CFI+YSM*SFI
1129 YMAG=YSM*CFI-XSM*SFI
1130 ZMAG=ZSM
1131 ENDIF
1132
1133 RETURN
1134 END
1135 C
1136 C=======================================================================================
1137 C
1138 SUBROUTINE GSMGSE (XGSM,YGSM,ZGSM,XGSE,YGSE,ZGSE,J)
1139 C
1140 C CONVERTS GEOCENTRIC SOLAR MAGNETOSPHERIC (GSM) COORDS TO SOLAR ECLIPTIC (GSE) ONES
1141 C OR VICA VERSA.
1142 C J>0 J<0
1143 C-----INPUT: J,XGSM,YGSM,ZGSM J,XGSE,YGSE,ZGSE
1144 C----OUTPUT: XGSE,YGSE,ZGSE XGSM,YGSM,ZGSM
1145 C
1146 C ATTENTION: SUBROUTINE RECALC MUST BE INVOKED BEFORE GSMGSE IN TWO CASES:
1147 C /A/ BEFORE THE FIRST TRANSFORMATION OF COORDINATES
1148 C /B/ IF THE VALUES OF IYEAR,IDAY,IHOUR,MIN,ISEC HAVE BEEN CHANGED
1149 C
1150 C LAST MODIFICATION: MARCH 31, 2003
1151 C
1152 C AUTHOR: N. A. TSYGANENKO
1153 C
1154 COMMON /GEOPACK1/ A(12),SHI,CHI,AB(13),BA(8)
1155 C
1156 IF (J.GT.0) THEN
1157 XGSE=XGSM
1158 YGSE=YGSM*CHI-ZGSM*SHI
1159 ZGSE=YGSM*SHI+ZGSM*CHI
1160 ELSE
1161 XGSM=XGSE
1162 YGSM=YGSE*CHI+ZGSE*SHI
1163 ZGSM=ZGSE*CHI-YGSE*SHI
1164 ENDIF
1165
1166 RETURN
1167 END
1168 C
1169 C=====================================================================================
1170 C
1171 SUBROUTINE SMGSM (XSM,YSM,ZSM,XGSM,YGSM,ZGSM,J)
1172 C
1173 C CONVERTS SOLAR MAGNETIC (SM) TO GEOCENTRIC SOLAR MAGNETOSPHERIC
1174 C (GSM) COORDINATES OR VICA VERSA.
1175 C J>0 J<0
1176 C-----INPUT: J,XSM,YSM,ZSM J,XGSM,YGSM,ZGSM
1177 C----OUTPUT: XGSM,YGSM,ZGSM XSM,YSM,ZSM
1178 C
1179 C ATTENTION: SUBROUTINE RECALC MUST BE INVOKED BEFORE SMGSM IN TWO CASES:
1180 C /A/ BEFORE THE FIRST TRANSFORMATION OF COORDINATES
1181 C /B/ IF THE VALUES OF IYEAR,IDAY,IHOUR,MIN,ISEC HAVE BEEN CHANGED
1182 C
1183 C LAST MODIFICATION: MARCH 31, 2003
1184 C
1185 C AUTHOR: N. A. TSYGANENKO
1186 C
1187 COMMON /GEOPACK1/ A(10),SPS,CPS,B(15),AB(8)
1188
1189 IF (J.GT.0) THEN
1190 XGSM=XSM*CPS+ZSM*SPS
1191 YGSM=YSM
1192 ZGSM=ZSM*CPS-XSM*SPS
1193 ELSE
1194 XSM=XGSM*CPS-ZGSM*SPS
1195 YSM=YGSM
1196 ZSM=XGSM*SPS+ZGSM*CPS
1197 ENDIF
1198
1199 RETURN
1200 END
1201 C
1202 C==========================================================================================
1203 C
1204 SUBROUTINE GEOGSM (XGEO,YGEO,ZGEO,XGSM,YGSM,ZGSM,J)
1205 C
1206 C CONVERTS GEOGRAPHIC (GEO) TO GEOCENTRIC SOLAR MAGNETOSPHERIC (GSM) COORDINATES
1207 C OR VICA VERSA.
1208 C
1209 C J>0 J<0
1210 C----- INPUT: J,XGEO,YGEO,ZGEO J,XGSM,YGSM,ZGSM
1211 C---- OUTPUT: XGSM,YGSM,ZGSM XGEO,YGEO,ZGEO
1212 C
1213 C ATTENTION: SUBROUTINE RECALC MUST BE INVOKED BEFORE GEOGSM IN TWO CASES:
1214 C /A/ BEFORE THE FIRST TRANSFORMATION OF COORDINATES
1215 C /B/ IF THE VALUES OF IYEAR,IDAY,IHOUR,MIN,ISEC HAVE BEEN CHANGED
1216 C
1217 C LAST MODIFICATION: MARCH 31, 2003
1218 C
1219 C AUTHOR: N. A. TSYGANENKO
1220 C
1221 COMMON /GEOPACK1/AA(17),A11,A21,A31,A12,A22,A32,A13,A23,A33,D,B(8)
1222 C
1223 IF (J.GT.0) THEN
1224 XGSM=A11*XGEO+A12*YGEO+A13*ZGEO
1225 YGSM=A21*XGEO+A22*YGEO+A23*ZGEO
1226 ZGSM=A31*XGEO+A32*YGEO+A33*ZGEO
1227 ELSE
1228 XGEO=A11*XGSM+A21*YGSM+A31*ZGSM
1229 YGEO=A12*XGSM+A22*YGSM+A32*ZGSM
1230 ZGEO=A13*XGSM+A23*YGSM+A33*ZGSM
1231 ENDIF
1232
1233 RETURN
1234 END
1235 C
1236 C=====================================================================================
1237 C
1238 SUBROUTINE RHAND (X,Y,Z,R1,R2,R3,IOPT,PARMOD,EXNAME,INNAME)
1239 C
1240 C CALCULATES THE COMPONENTS OF THE RIGHT HAND SIDE VECTOR IN THE GEOMAGNETIC FIELD
1241 C LINE EQUATION (a subsidiary subroutine for the subroutine STEP)
1242 C
1243 C LAST MODIFICATION: MARCH 31, 2003
1244 C
1245 C AUTHOR: N. A. TSYGANENKO
1246 C
1247 DIMENSION PARMOD(10)
1248 C
1249 C EXNAME AND INNAME ARE NAMES OF SUBROUTINES FOR THE EXTERNAL AND INTERNAL
1250 C PARTS OF THE TOTAL FIELD
1251 C
1252 COMMON /GEOPACK1/ A(15),PSI,AA(10),DS3,BB(8)
1253
1254 CALL EXNAME (IOPT,PARMOD,PSI,X,Y,Z,BXGSM,BYGSM,BZGSM)
1255 CALL INNAME (X,Y,Z,HXGSM,HYGSM,HZGSM)
1256
1257 BX=BXGSM+HXGSM
1258 BY=BYGSM+HYGSM
1259 BZ=BZGSM+HZGSM
1260 B=DS3/SQRT(BX**2+BY**2+BZ**2)
1261 R1=BX*B
1262 R2=BY*B
1263 R3=BZ*B
1264 RETURN
1265 END
1266 C
1267 C===================================================================================
1268 C
1269 SUBROUTINE STEP (X,Y,Z,DS,ERRIN,IOPT,PARMOD,EXNAME,INNAME)
1270 C
1271 C RE-CALCULATES {X,Y,Z}, MAKING A STEP ALONG A FIELD LINE.
1272 C DS IS THE STEP SIZE, ERRIN IS PERMISSIBLE ERROR VALUE, IOPT SPECIFIES THE EXTERNAL
1273 C MODEL VERSION, THE ARRAY PARMOD CONTAINS INPUT PARAMETERS FOR THAT MODEL
1274 C EXNAME IS THE NAME OF THE EXTERNAL FIELD SUBROUTINE
1275 C INNAME IS THE NAME OF THE INTERNAL FIELD SUBROUTINE (EITHER DIP OR IGRF)
1276 C
1277 C ALL THE PARAMETERS ARE INPUT ONES; OUTPUT IS THE RENEWED TRIPLET X,Y,Z
1278 C
1279 C LAST MODIFICATION: MARCH 31, 2003
1280 C
1281 C AUTHOR: N. A. TSYGANENKO
1282 C
1283 DIMENSION PARMOD(10)
1284 COMMON /GEOPACK1/ A(26),DS3,B(8)
1285 EXTERNAL EXNAME,INNAME
1286
1287 1 DS3=-DS/3.
1288 c (nico) I add this to avoid a seg fault
1289 WRITE(100,*)
1290 CALL RHAND (X,Y,Z,R11,R12,R13,IOPT,PARMOD,EXNAME,INNAME)
1291 CALL RHAND (X+R11,Y+R12,Z+R13,R21,R22,R23,IOPT,PARMOD,EXNAME,
1292 * INNAME)
1293 CALL RHAND (X+.5*(R11+R21),Y+.5*(R12+R22),Z+.5*
1294 *(R13+R23),R31,R32,R33,IOPT,PARMOD,EXNAME,INNAME)
1295 CALL RHAND (X+.375*(R11+3.*R31),Y+.375*(R12+3.*R32
1296 *),Z+.375*(R13+3.*R33),R41,R42,R43,IOPT,PARMOD,EXNAME,INNAME)
1297 CALL RHAND (X+1.5*(R11-3.*R31+4.*R41),Y+1.5*(R12-
1298 *3.*R32+4.*R42),Z+1.5*(R13-3.*R33+4.*R43),
1299 *R51,R52,R53,IOPT,PARMOD,EXNAME,INNAME)
1300 ERRCUR=ABS(R11-4.5*R31+4.*R41-.5*R51)+ABS(R12-4.5*R32+4.*R42-.5*
1301 *R52)+ABS(R13-4.5*R33+4.*R43-.5*R53)
1302 IF (ERRCUR.LT.ERRIN) GOTO 2
1303 DS=DS*.5
1304 GOTO 1
1305 2 X=X+.5*(R11+4.*R41+R51)
1306 Y=Y+.5*(R12+4.*R42+R52)
1307 Z=Z+.5*(R13+4.*R43+R53)
1308 IF(ERRCUR.LT.ERRIN*.04.AND.ABS(DS).LT.1.33) DS=DS*1.5
1309 RETURN
1310 END
1311 C
1312 C==============================================================================
1313 C
1314 SUBROUTINE TRACE (XI,YI,ZI,DIR,RLIM,R0,IOPT,PARMOD,EXNAME,INNAME,
1315 *XF,YF,ZF,XX,YY,ZZ,L)
1316 C
1317 C TRACES A FIELD LINE FROM AN ARBITRARY POINT OF SPACE TO THE EARTH'S
1318 C SURFACE OR TO A MODEL LIMITING BOUNDARY.
1319 C
1320 C THE HIGHEST ORDER OF SPHERICAL HARMONICS IN THE MAIN FIELD EXPANSION USED
1321 C IN THE MAPPING IS CALCULATED AUTOMATICALLY. IF INNAME=IGRF_GSM, THEN AN IGRF MODEL
1322 C FIELD WILL BE USED, AND IF INNAME=DIP, A PURE DIPOLE FIELD WILL BE USED.
1323
1324 C IN ANY CASE, BEFORE CALLING TRACE, ONE SHOULD INVOKE RECALC, TO CALCULATE CORRECT
1325 C VALUES OF THE IGRF COEFFICIENTS AND ALL QUANTITIES NEEDED FOR TRANSFORMATIONS
1326 C BETWEEN COORDINATE SYSTEMS INVOLVED IN THIS CALCULATIONS.
1327 C
1328 C ALTERNATIVELY, THE SUBROUTINE RECALC CAN BE INVOKED WITH THE DESIRED VALUES OF
1329 C IYEAR AND IDAY (TO SPECIFY THE DIPOLE MOMENT), WHILE THE VALUES OF THE DIPOLE
1330 C TILT ANGLE PSI (IN RADIANS) AND ITS SINE (SPS) AND COSINE (CPS) CAN BE EXPLICITLY
1331 C SPECIFIED AND FORWARDED TO THE COMMON BLOCK GEOPACK1 (11th, 12th, AND 16th ELEMENTS, RESP.)
1332 C
1333 C------------- INPUT PARAMETERS:
1334 C
1335 C XI,YI,ZI - GSM COORDS OF INITIAL POINT (IN EARTH RADII, 1 RE = 6371.2 km),
1336 C
1337 C DIR - SIGN OF THE TRACING DIRECTION: IF DIR=1.0 THEN WE MOVE ANTIPARALLEL TO THE
1338 C FIELD VECTOR (E.G. FROM NORTHERN TO SOUTHERN CONJUGATE POINT),
1339 C AND IF DIR=-1.0 THEN THE TRACING GOES IN THE OPPOSITE DIRECTION.
1340 C
1341 C R0 - RADIUS OF A SPHERE (IN RE) FOR WHICH THE FIELD LINE ENDPOINT COORDINATES
1342 C XF,YF,ZF SHOULD BE CALCULATED.
1343 C
1344 C RLIM - UPPER LIMIT OF THE GEOCENTRIC DISTANCE, WHERE THE TRACING IS TERMINATED.
1345 C
1346 C IOPT - A MODEL INDEX; CAN BE USED FOR SPECIFYING AN OPTION OF THE EXTERNAL FIELD
1347 C MODEL (E.G., INTERVAL OF THE KP-INDEX). ALTERNATIVELY, ONE CAN USE THE ARRAY
1348 C PARMOD FOR THAT PURPOSE (SEE BELOW); IN THAT CASE IOPT IS JUST A DUMMY PARAMETER.
1349 C
1350 C PARMOD - A 10-ELEMENT ARRAY CONTAINING MODEL PARAMETERS, NEEDED FOR A UNIQUE
1351 C SPECIFICATION OF THE EXTERNAL FIELD. THE CONCRETE MEANING OF THE COMPONENTS
1352 C OF PARMOD DEPENDS ON A SPECIFIC VERSION OF THE EXTERNAL FIELD MODEL.
1353 C
1354 C EXNAME - NAME OF A SUBROUTINE PROVIDING COMPONENTS OF THE EXTERNAL MAGNETIC FIELD
1355 C (E.G., T96_01).
1356 C INNAME - NAME OF A SUBROUTINE PROVIDING COMPONENTS OF THE INTERNAL MAGNETIC FIELD
1357 C (EITHER DIP OR IGRF_GSM).
1358 C
1359 C-------------- OUTPUT PARAMETERS:
1360 C
1361 C XF,YF,ZF - GSM COORDS OF THE LAST CALCULATED POINT OF A FIELD LINE
1362 C XX,YY,ZZ - ARRAYS, CONTAINING COORDS OF FIELD LINE POINTS. HERE THEIR MAXIMAL LENGTH WAS
1363 C ASSUMED EQUAL TO 999.
1364 C L - ACTUAL NUMBER OF THE CALCULATED FIELD LINE POINTS. IF L EXCEEDS 999, TRACING
1365 C TERMINATES, AND A WARNING IS DISPLAYED.
1366 C
1367 C
1368 C LAST MODIFICATION: MARCH 31, 2003.
1369 C
1370 C AUTHOR: N. A. TSYGANENKO
1371 C
1372 DIMENSION XX(1000),YY(1000),ZZ(1000), PARMOD(10)
1373 COMMON /GEOPACK1/ AA(26),DD,BB(8)
1374 EXTERNAL EXNAME,INNAME
1375 C
1376 ERR=0.0001
1377 L=0
1378 DS=0.5*DIR
1379 X=XI
1380 Y=YI
1381 Z=ZI
1382 DD=DIR
1383 AL=0.
1384 c
1385 c here we call RHAND just to find out the sign of the radial component of the field
1386 c vector, and to determine the initial direction of the tracing (i.e., either away
1387 c or towards Earth):
1388 c
1389 WRITE(100,*)
1390 CALL RHAND (X,Y,Z,R1,R2,R3,IOPT,PARMOD,EXNAME,INNAME)
1391 AD=0.01
1392 IF (X*R1+Y*R2+Z*R3.LT.0.) AD=-0.01
1393 C
1394 c |AD|=0.01 and its sign follows the rule:
1395 c (1) if DIR=1 (tracing antiparallel to B vector) then the sign of AD is the same as of Br
1396 c (2) if DIR=-1 (tracing parallel to B vector) then the sign of AD is opposite to that of Br
1397 c AD is defined in order to initialize the value of RR (radial distance at previous step):
1398
1399 RR=SQRT(X**2+Y**2+Z**2)+AD
1400 1 L=L+1
1401 IF(L.GT.999) GOTO 7
1402 XX(L)=X
1403 YY(L)=Y
1404 ZZ(L)=Z
1405 RYZ=Y**2+Z**2
1406 R2=X**2+RYZ
1407 R=SQRT(R2)
1408
1409 c check if the line hit the outer tracing boundary; if yes, then terminate
1410 c the tracing (label 8):
1411
1412 IF (R.GT.RLIM.OR.RYZ.GT.1600.D0.OR.X.GT.20.D0) GOTO 8
1413 c
1414 c check whether or not the inner tracing boundary was crossed from outside,
1415 c if yes, then calculate the footpoint position by interpolation (go to label 6):
1416 c
1417 IF (R.LT.R0.AND.RR.GT.R) GOTO 6
1418
1419 c check if (i) we are moving outward, or (ii) we are still sufficiently
1420 c far from Earth (beyond R=5Re); if yes, proceed further:
1421 c
1422 IF (R.GE.RR.OR.R.GT.5.) GOTO 5
1423
1424 c now we moved closer inward (between R=3 and R=5); go to 3 and begin logging
1425 c previous values of X,Y,Z, to be used in the interpolation (after having
1426 c crossed the inner tracing boundary):
1427
1428 IF (R.GE.3.) GOTO 3
1429 c
1430 c we entered inside the sphere R=3: to avoid too large steps (and hence inaccurate
1431 c interpolated position of the footpoint), enforce the progressively smaller
1432 c stepsize values as we approach the inner boundary R=R0:
1433 c
1434 FC=0.2
1435 IF(R-R0.LT.0.05) FC=0.05
1436 AL=FC*(R-R0+0.2)
1437 DS=DIR*AL
1438 GOTO 4
1439 3 DS=DIR
1440 4 XR=X
1441 YR=Y
1442 ZR=Z
1443 5 RR=R
1444 CALL STEP (X,Y,Z,DS,ERR,IOPT,PARMOD,EXNAME,INNAME)
1445 GOTO 1
1446 c
1447 c find the footpoint position by interpolating between the current and previous
1448 c field line points:
1449 c
1450 6 R1=(R0-R)/(RR-R)
1451 X=X-(X-XR)*R1
1452 Y=Y-(Y-YR)*R1
1453 Z=Z-(Z-ZR)*R1
1454 GOTO 8
1455 7 WRITE (*,10)
1456 L=999
1457 8 XF=X
1458 YF=Y
1459 ZF=Z
1460 RETURN
1461 10 FORMAT(//,1X,'**** COMPUTATIONS IN THE SUBROUTINE TRACE ARE',
1462 *' TERMINATED: THE CURRENT NUMBER OF POINTS EXCEEDED 1000 ****'//)
1463 END
1464 c
1465 C====================================================================================
1466 C
1467 SUBROUTINE SHUETAL_MGNP(XN_PD,VEL,BZIMF,XGSM,YGSM,ZGSM,
1468 * XMGNP,YMGNP,ZMGNP,DIST,ID)
1469 C
1470 C FOR ANY POINT OF SPACE WITH COORDINATES (XGSM,YGSM,ZGSM) AND SPECIFIED CONDITIONS
1471 C IN THE INCOMING SOLAR WIND, THIS SUBROUTINE:
1472 C
1473 C (1) DETERMINES IF THE POINT (XGSM,YGSM,ZGSM) LIES INSIDE OR OUTSIDE THE
1474 C MODEL MAGNETOPAUSE OF SHUE ET AL. (JGR-A, V.103, P. 17691, 1998).
1475 C
1476 C (2) CALCULATES THE GSM POSITION OF A POINT {XMGNP,YMGNP,ZMGNP}, LYING AT THE MODEL
1477 C MAGNETOPAUSE AND ASYMPTOTICALLY TENDING TO THE NEAREST BOUNDARY POINT WITH
1478 C RESPECT TO THE OBSERVATION POINT {XGSM,YGSM,ZGSM}, AS IT APPROACHES THE MAGNETO-
1479 C PAUSE.
1480 C
1481 C INPUT: XN_PD - EITHER SOLAR WIND PROTON NUMBER DENSITY (PER C.C.) (IF VEL>0)
1482 C OR THE SOLAR WIND RAM PRESSURE IN NANOPASCALS (IF VEL<0)
1483 C BZIMF - IMF BZ IN NANOTESLAS
1484 C
1485 C VEL - EITHER SOLAR WIND VELOCITY (KM/SEC)
1486 C OR ANY NEGATIVE NUMBER, WHICH INDICATES THAT XN_PD STANDS
1487 C FOR THE SOLAR WIND PRESSURE, RATHER THAN FOR THE DENSITY
1488 C
1489 C XGSM,YGSM,ZGSM - GSM POSITION OF THE OBSERVATION POINT IN EARTH RADII
1490 C
1491 C OUTPUT: XMGNP,YMGNP,ZMGNP - GSM POSITION OF THE BOUNDARY POINT
1492 C DIST - DISTANCE (IN RE) BETWEEN THE OBSERVATION POINT (XGSM,YGSM,ZGSM)
1493 C AND THE MODEL NAGNETOPAUSE
1494 C ID - POSITION FLAG: ID=+1 (-1) MEANS THAT THE OBSERVATION POINT
1495 C LIES INSIDE (OUTSIDE) OF THE MODEL MAGNETOPAUSE, RESPECTIVELY.
1496 C
1497 C OTHER SUBROUTINES USED: T96_MGNP
1498 C
1499 c AUTHOR: N.A. TSYGANENKO,
1500 C DATE: APRIL 4, 2003.
1501 C
1502 IF (VEL.LT.0.) THEN
1503 PD=XN_PD
1504 ELSE
1505 PD=1.94E-6*XN_PD*VEL**2 ! PD IS THE SOLAR WIND DYNAMIC PRESSURE (IN nPa)
1506 ENDIF
1507
1508 c
1509 c DEFINE THE ANGLE PHI, MEASURED DUSKWARD FROM THE NOON-MIDNIGHT MERIDIAN PLANE;
1510 C IF THE OBSERVATION POINT LIES ON THE X AXIS, THE ANGLE PHI CANNOT BE UNIQUELY
1511 C DEFINED, AND WE SET IT AT ZERO:
1512 c
1513 IF (YGSM.NE.0..OR.ZGSM.NE.0.) THEN
1514 PHI=ATAN2(YGSM,ZGSM)
1515 ELSE
1516 PHI=0.
1517 ENDIF
1518 C
1519 C FIRST, FIND OUT IF THE OBSERVATION POINT LIES INSIDE THE SHUE ET AL BDRY
1520 C AND SET THE VALUE OF THE ID FLAG:
1521 C
1522 ID=-1
1523 R0=(10.22+1.29*TANH(0.184*(BZIMF+8.14)))*PD**(-.15151515)
1524 ALPHA=(0.58-0.007*BZIMF)*(1.+0.024*ALOG(PD))
1525 R=SQRT(XGSM**2+YGSM**2+ZGSM**2)
1526 RM=R0*(2./(1.+XGSM/R))**ALPHA
1527 IF (R.LE.RM) ID=+1
1528 C
1529 C NOW, FIND THE CORRESPONDING T96 MAGNETOPAUSE POSITION, TO BE USED AS
1530 C A STARTING APPROXIMATION IN THE SEARCH OF A CORRESPONDING SHUE ET AL.
1531 C BOUNDARY POINT:
1532 C
1533 CALL T96_MGNP (PD,-1.,XGSM,YGSM,ZGSM,XMT96,YMT96,ZMT96,DIST,ID96)
1534 C
1535 RHO2=YMT96**2+ZMT96**2
1536 R=SQRT(RHO2+XMT96**2)
1537 ST=SQRT(RHO2)/R
1538 CT=XMT96/R
1539 C
1540 C NOW, USE NEWTON'S ITERATIVE METHOD TO FIND THE NEAREST POINT AT THE
1541 C SHUE ET AL.'S BOUNDARY:
1542 C
1543 NIT=0
1544
1545 1 T=ATAN2(ST,CT)
1546 RM=R0*(2./(1.+CT))**ALPHA
1547
1548 F=R-RM
1549 GRADF_R=1.
1550 GRADF_T=-ALPHA/R*RM*ST/(1.+CT)
1551 GRADF=SQRT(GRADF_R**2+GRADF_T**2)
1552
1553 DR=-F/GRADF**2
1554 DT= DR/R*GRADF_T
1555
1556 R=R+DR
1557 T=T+DT
1558 ST=SIN(T)
1559 CT=COS(T)
1560
1561 DS=SQRT(DR**2+(R*DT)**2)
1562
1563 NIT=NIT+1
1564
1565 IF (NIT.GT.1000) THEN
1566 PRINT *,
1567 *' BOUNDARY POINT COULD NOT BE FOUND; ITERATIONS DO NOT CONVERGE'
1568 ENDIF
1569
1570 IF (DS.GT.1.E-4) GOTO 1
1571
1572 XMGNP=R*COS(T)
1573 RHO= R*SIN(T)
1574
1575 YMGNP=RHO*SIN(PHI)
1576 ZMGNP=RHO*COS(PHI)
1577
1578 DIST=SQRT((XGSM-XMGNP)**2+(YGSM-YMGNP)**2+(ZGSM-ZMGNP)**2)
1579
1580 RETURN
1581 END
1582 C
1583 C=======================================================================================
1584 C
1585 SUBROUTINE T96_MGNP (XN_PD,VEL,XGSM,YGSM,ZGSM,XMGNP,YMGNP,ZMGNP,
1586 * DIST,ID)
1587 C
1588 C FOR ANY POINT OF SPACE WITH GIVEN COORDINATES (XGSM,YGSM,ZGSM), THIS SUBROUTINE DEFINES
1589 C THE POSITION OF A POINT (XMGNP,YMGNP,ZMGNP) AT THE T96 MODEL MAGNETOPAUSE, HAVING THE
1590 C SAME VALUE OF THE ELLIPSOIDAL TAU-COORDINATE, AND THE DISTANCE BETWEEN THEM. THIS IS
1591 C NOT THE SHORTEST DISTANCE D_MIN TO THE BOUNDARY, BUT DIST ASYMPTOTICALLY TENDS TO D_MIN,
1592 C AS THE OBSERVATION POINT GETS CLOSER TO THE MAGNETOPAUSE.
1593 C
1594 C INPUT: XN_PD - EITHER SOLAR WIND PROTON NUMBER DENSITY (PER C.C.) (IF VEL>0)
1595 C OR THE SOLAR WIND RAM PRESSURE IN NANOPASCALS (IF VEL<0)
1596 C VEL - EITHER SOLAR WIND VELOCITY (KM/SEC)
1597 C OR ANY NEGATIVE NUMBER, WHICH INDICATES THAT XN_PD STANDS
1598 C FOR THE SOLAR WIND PRESSURE, RATHER THAN FOR THE DENSITY
1599 C
1600 C XGSM,YGSM,ZGSM - COORDINATES OF THE OBSERVATION POINT IN EARTH RADII
1601 C
1602 C OUTPUT: XMGNP,YMGNP,ZMGNP - GSM POSITION OF THE BOUNDARY POINT, HAVING THE SAME
1603 C VALUE OF TAU-COORDINATE AS THE OBSERVATION POINT (XGSM,YGSM,ZGSM)
1604 C DIST - THE DISTANCE BETWEEN THE TWO POINTS, IN RE,
1605 C ID - POSITION FLAG; ID=+1 (-1) MEANS THAT THE POINT (XGSM,YGSM,ZGSM)
1606 C LIES INSIDE (OUTSIDE) THE MODEL MAGNETOPAUSE, RESPECTIVELY.
1607 C
1608 C THE PRESSURE-DEPENDENT MAGNETOPAUSE IS THAT USED IN THE T96_01 MODEL
1609 C (TSYGANENKO, JGR, V.100, P.5599, 1995; ESA SP-389, P.181, OCT. 1996)
1610 C
1611 c AUTHOR: N.A. TSYGANENKO
1612 C DATE: AUG.1, 1995, REVISED APRIL 3, 2003.
1613 C
1614 C
1615 C DEFINE SOLAR WIND DYNAMIC PRESSURE (NANOPASCALS, ASSUMING 4% OF ALPHA-PARTICLES),
1616 C IF NOT EXPLICITLY SPECIFIED IN THE INPUT:
1617
1618 IF (VEL.LT.0.) THEN
1619 PD=XN_PD
1620 ELSE
1621 PD=1.94E-6*XN_PD*VEL**2
1622 C
1623 ENDIF
1624 C
1625 C RATIO OF PD TO THE AVERAGE PRESSURE, ASSUMED EQUAL TO 2 nPa:
1626
1627 RAT=PD/2.0
1628 RAT16=RAT**0.14
1629
1630 C (THE POWER INDEX 0.14 IN THE SCALING FACTOR IS THE BEST-FIT VALUE OBTAINED FROM DATA
1631 C AND USED IN THE T96_01 VERSION)
1632 C
1633 C VALUES OF THE MAGNETOPAUSE PARAMETERS FOR PD = 2 nPa:
1634 C
1635 A0=70.
1636 S00=1.08
1637 X00=5.48
1638 C
1639 C VALUES OF THE MAGNETOPAUSE PARAMETERS, SCALED BY THE ACTUAL PRESSURE:
1640 C
1641 A=A0/RAT16
1642 S0=S00
1643 X0=X00/RAT16
1644 XM=X0-A
1645 C
1646 C (XM IS THE X-COORDINATE OF THE "SEAM" BETWEEN THE ELLIPSOID AND THE CYLINDER)
1647 C
1648 C (FOR DETAILS ON THE ELLIPSOIDAL COORDINATES, SEE THE PAPER:
1649 C N.A.TSYGANENKO, SOLUTION OF CHAPMAN-FERRARO PROBLEM FOR AN
1650 C ELLIPSOIDAL MAGNETOPAUSE, PLANET.SPACE SCI., V.37, P.1037, 1989).
1651 C
1652 IF (YGSM.NE.0..OR.ZGSM.NE.0.) THEN
1653 PHI=ATAN2(YGSM,ZGSM)
1654 ELSE
1655 PHI=0.
1656 ENDIF
1657 C
1658 RHO=SQRT(YGSM**2+ZGSM**2)
1659 C
1660 IF (XGSM.LT.XM) THEN
1661 XMGNP=XGSM
1662 RHOMGNP=A*SQRT(S0**2-1)
1663 YMGNP=RHOMGNP*SIN(PHI)
1664 ZMGNP=RHOMGNP*COS(PHI)
1665 DIST=SQRT((XGSM-XMGNP)**2+(YGSM-YMGNP)**2+(ZGSM-ZMGNP)**2)
1666 IF (RHOMGNP.GT.RHO) ID=+1
1667 IF (RHOMGNP.LE.RHO) ID=-1
1668 RETURN
1669 ENDIF
1670 C
1671 XKSI=(XGSM-X0)/A+1.
1672 XDZT=RHO/A
1673 SQ1=SQRT((1.+XKSI)**2+XDZT**2)
1674 SQ2=SQRT((1.-XKSI)**2+XDZT**2)
1675 SIGMA=0.5*(SQ1+SQ2)
1676 TAU=0.5*(SQ1-SQ2)
1677 C
1678 C NOW CALCULATE (X,Y,Z) FOR THE CLOSEST POINT AT THE MAGNETOPAUSE
1679 C
1680 XMGNP=X0-A*(1.-S0*TAU)
1681 ARG=(S0**2-1.)*(1.-TAU**2)
1682 IF (ARG.LT.0.) ARG=0.
1683 RHOMGNP=A*SQRT(ARG)
1684 YMGNP=RHOMGNP*SIN(PHI)
1685 ZMGNP=RHOMGNP*COS(PHI)
1686 C
1687 C NOW CALCULATE THE DISTANCE BETWEEN THE POINTS {XGSM,YGSM,ZGSM} AND {XMGNP,YMGNP,ZMGNP}:
1688 C (IN GENERAL, THIS IS NOT THE SHORTEST DISTANCE D_MIN, BUT DIST ASYMPTOTICALLY TENDS
1689 C TO D_MIN, AS WE ARE GETTING CLOSER TO THE MAGNETOPAUSE):
1690 C
1691 DIST=SQRT((XGSM-XMGNP)**2+(YGSM-YMGNP)**2+(ZGSM-ZMGNP)**2)
1692 C
1693 IF (SIGMA.GT.S0) ID=-1 ! ID=-1 MEANS THAT THE POINT LIES OUTSIDE
1694 IF (SIGMA.LE.S0) ID=+1 ! ID=+1 MEANS THAT THE POINT LIES INSIDE
1695 C THE MAGNETOSPHERE
1696 RETURN
1697 END
1698 C
1699 C -----------------------------------------------------------------------------
1700 C Subroutine SETPSI added by PJ 15.12.2003, copied from geopack2003.f 20.5.2005
1701 C -----------------------------------------------------------------------------
1702 SUBROUTINE SETPSI(PS)
1703 REAL PS
1704 COMMON /GEOPACK1/ AAA(10),SPS,CPS,BBB(23)
1705 REAL AAA,SPS,CPS,BBB
1706 SPS=SIN(PS)
1707 CPS=COS(PS)
1708 RETURN
1709 END
1710 C ----------------------------------------
1711 C
1712 C===================================================================================
1713 C
1714 c</pre>

  ViewVC Help
Powered by ViewVC 1.1.23