| 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> |
Parent Directory
| 
Revision Log