CDPP / AMDA_COTS

Commit 55c4dfc0ef263687d108e47e36d46a3ba70b0c2e

Authored by Benjamin Renard
1 parent e7cf8b1f
Exists in master and in 2 other branches add_gcc_archives, test_include_boost

Add COTS for Tsyganenko96 model magnetic field calculation

Showing 6 changed files with 4687 additions and 0 deletions   Show diff stats
Inline Side-by-side
1 netcdf/tmp/ 1 netcdf/tmp/
2 cspice/cspice/ 2 cspice/cspice/
3 build/ 3 build/
  4 +ceflib/CESR_CEFLIB/
  5 +gcovr/gcovr-3.2/
geopack/CMakeLists.txt 0 → 100644
  Diff comments   View file @55c4dfc
@@ -0,0 +1,22 @@ @@ -0,0 +1,22 @@
  1 +cmake_minimum_required(VERSION 2.6)
  2 +
  3 +PROJECT(geopack CXX Fortran)
  4 +
  5 +add_definitions( -DLINUX )
  6 +
  7 +set(CMAKE_BUILD_TYPE Release)
  8 +
  9 +set(CMAKE_CXX_FLAGS "-O3 -fPIC -fno-align-commons -fbounds-check -fbacktrace -ffixed-form")
  10 +
  11 +file(
  12 + GLOB
  13 + source_files
  14 + src/*.for
  15 +)
  16 +
  17 +add_library(geopack SHARED ${source_files})
  18 +
  19 +SET_TARGET_PROPERTIES(geopack PROPERTIES LINKER_LANGUAGE CXX)
  20 +
  21 +install (TARGETS geopack DESTINATION lib)
  22 +install (FILES src/include/geopack.hh DESTINATION include)
0 \ No newline at end of file 23 \ No newline at end of file
geopack/src/Geopack-2008.for 0 → 100644
  Diff comments   View file @55c4dfc
@@ -0,0 +1,2053 @@ @@ -0,0 +1,2053 @@
  1 +c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
  2 +c
  3 +c
  4 +c ##########################################################################
  5 +c # #
  6 +c # GEOPACK-2008 #
  7 +c # (MAIN SET OF FORTRAN CODES) #
  8 +c # #
  9 +c ##########################################################################
  10 +C
  11 +c
  12 +c This collection of subroutines is a result of several upgrades of the original package
  13 +c written by N. A. Tsyganenko in 1978-1979.
  14 +c
  15 +c PREFATORY NOTE TO THE VERSION OF FEBRUARY 8, 2008:
  16 +c
  17 +c To avoid inappropriate use of obsolete subroutines from earlier versions, a suffix 08 was
  18 +c added to the name of each subroutine in this release.
  19 +c
  20 +c A possibility has been added in this version to calculate vector components in the
  21 +c "Geocentric Solar Wind" (GSW) coordinate system, which, to our knowledge, was first
  22 +c introduced by Hones et al., Planet. Space Sci., v.34, p.889, 1986 (aka GSWM, see Appendix,
  23 +c Tsyganenko et al., JGRA, v.103(A4), p.6827, 1998). The GSW system is analogous to the
  24 +c standard GSM, except that its X-axis is antiparallel to the currently observed solar wind
  25 +c flow vector, rather than aligned with the Earth-Sun line. The orientation of axes in the
  26 +c GSW system can be uniquely defined by specifying three components (VGSEX,VGSEY,VGSEZ) of
  27 +c the solar wind velocity, and in the case of a strictly radial anti-sunward flow (VGSEY=
  28 +c VGSEZ=0) the GSW system becomes identical to the standard GSM, which fact was used here
  29 +c to minimize the number of subroutines in the package. To that end, instead of the special
  30 +c case of the GSM coordinates, this version uses a more general GSW system, and three more
  31 +c input parameters are added in the subroutine RECALC_08, the observed values (VGSEX,VGSEY,
  32 +c VGSEZ) of the solar wind velocity. Invoking RECALC_08 with VGSEY=VGSEZ=0 restores the
  33 +c standard (sunward) orientation of the X axis, which allows one to easily convert vectors
  34 +c between GSW and GSM, as well as to/from other standard and commonly used systems. For more
  35 +c details, see the documentation file GEOPACK-2008.DOC.
  36 +c
  37 +c Another modification allows users to have more control over the procedure of field line
  38 +c mapping using the subroutine TRACE_08. To that end, three new input parameters were added
  39 +c in that subroutine, allowing one to set (i) an upper limit, DSMAX, on the automatically
  40 +c adjusted step size, (ii) a permissible step error, ERR, and (iii) maximal length, LMAX,
  41 +c of arrays where field line point coordinates are stored. Minor changes were also made in
  42 +c the tracing subroutine, to make it more compact and easier for understanding, and to
  43 +c prevent the algorithm from making uncontrollable large number of multiple loops in some
  44 +c cases with plasmoid-like field structures.
  45 +c
  46 +C One more subroutine, named GEODGEO_08, was added to the package, allowing one to convert
  47 +c geodetic coordinates of a point in space (altitude above the Earth's WGS84 ellipsoid and
  48 +c geodetic latitude) to geocentric radial distance and colatitude, and vice versa.
  49 +c
  50 +C For a complete list of modifications made earlier in previous versions, see the
  51 +c documentation file GEOPACK-2008.DOC.
  52 +c
  53 +c----------------------------------------------------------------------------------
  54 +c
  55 + SUBROUTINE IGRF_GSW_08 (XGSW,YGSW,ZGSW,HXGSW,HYGSW,HZGSW)
  56 +c
  57 +C CALCULATES COMPONENTS OF THE MAIN (INTERNAL) GEOMAGNETIC FIELD IN THE GEOCENTRIC SOLAR-WIND
  58 +C (GSW) COORDINATE SYSTEM, USING IAGA INTERNATIONAL GEOMAGNETIC REFERENCE MODEL COEFFICIENTS
  59 +C (e.g., http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html, revised 22 March, 2005)
  60 +c
  61 +C THE GSW SYSTEM IS ESSENTIALLY SIMILAR TO THE STANDARD GSM (THE TWO SYSTEMS BECOME IDENTICAL
  62 +C TO EACH OTHER IN THE CASE OF STRICTLY ANTI-SUNWARD SOLAR WIND FLOW). FOR A DETAILED
  63 +C DEFINITION, SEE INTRODUCTORY COMMENTS FOR THE SUBROUTINE GSWGSE_08 .
  64 +C
  65 +C BEFORE THE FIRST CALL OF THIS SUBROUTINE, OR, IF THE DATE/TIME (IYEAR,IDAY,IHOUR,MIN,ISEC),
  66 +C OR THE SOLAR WIND VELOCITY COMPONENTS (VGSEX,VGSEY,VGSEZ) HAVE CHANGED, THE MODEL COEFFICIENTS
  67 +c AND GEO-GSW ROTATION MATRIX ELEMENTS SHOULD BE UPDATED BY CALLING THE SUBROUTINE RECALC_08.
  68 +C
  69 +C-----INPUT PARAMETERS:
  70 +C
  71 +C XGSW,YGSW,ZGSW - CARTESIAN GEOCENTRIC SOLAR-WIND COORDINATES (IN UNITS RE=6371.2 KM)
  72 +C
  73 +C-----OUTPUT PARAMETERS:
  74 +C
  75 +C HXGSW,HYGSW,HZGSW - CARTESIAN GEOCENTRIC SOLAR-WIND COMPONENTS OF THE MAIN GEOMAGNETIC
  76 +C FIELD IN NANOTESLA
  77 +C
  78 +C LAST MODIFICATION: FEB 07, 2008.
  79 +C THIS VERSION OF THE CODE ACCEPTS DATES FROM 1965 THROUGH 2015.
  80 +c
  81 +C AUTHOR: N. A. TSYGANENKO
  82 +C
  83 +C
  84 + COMMON /GEOPACK2/ G(105),H(105),REC(105)
  85 +
  86 + DIMENSION A(14),B(14)
  87 +
  88 + CALL GEOGSW_08 (XGEO,YGEO,ZGEO,XGSW,YGSW,ZGSW,-1)
  89 + RHO2=XGEO**2+YGEO**2
  90 + R=SQRT(RHO2+ZGEO**2)
  91 + C=ZGEO/R
  92 + RHO=SQRT(RHO2)
  93 + S=RHO/R
  94 + IF (S.LT.1.E-5) THEN
  95 + CF=1.
  96 + SF=0.
  97 + ELSE
  98 + CF=XGEO/RHO
  99 + SF=YGEO/RHO
  100 + ENDIF
  101 +C
  102 + PP=1./R
  103 + P=PP
  104 +C
  105 +C IN THIS VERSION, THE OPTIMAL VALUE OF THE PARAMETER NM (MAXIMAL ORDER OF THE SPHERICAL
  106 +C HARMONIC EXPANSION) IS NOT USER-PRESCRIBED, BUT CALCULATED INSIDE THE SUBROUTINE, BASED
  107 +C ON THE VALUE OF THE RADIAL DISTANCE R:
  108 +C
  109 + IRP3=R+2
  110 + NM=3+30/IRP3
  111 + IF (NM.GT.13) NM=13
  112 +
  113 + K=NM+1
  114 + DO 150 N=1,K
  115 + P=P*PP
  116 + A(N)=P
  117 +150 B(N)=P*N
  118 +
  119 + P=1.
  120 + D=0.
  121 + BBR=0.
  122 + BBT=0.
  123 + BBF=0.
  124 +
  125 + DO 200 M=1,K
  126 + IF(M.EQ.1) GOTO 160
  127 + MM=M-1
  128 + W=X
  129 + X=W*CF+Y*SF
  130 + Y=Y*CF-W*SF
  131 + GOTO 170
  132 +160 X=0.
  133 + Y=1.
  134 +170 Q=P
  135 + Z=D
  136 + BI=0.
  137 + P2=0.
  138 + D2=0.
  139 + DO 190 N=M,K
  140 + AN=A(N)
  141 + MN=N*(N-1)/2+M
  142 + E=G(MN)
  143 + HH=H(MN)
  144 + W=E*Y+HH*X
  145 + BBR=BBR+B(N)*W*Q
  146 + BBT=BBT-AN*W*Z
  147 + IF(M.EQ.1) GOTO 180
  148 + QQ=Q
  149 + IF(S.LT.1.E-5) QQ=Z
  150 + BI=BI+AN*(E*X-HH*Y)*QQ
  151 +180 XK=REC(MN)
  152 + DP=C*Z-S*Q-XK*D2
  153 + PM=C*Q-XK*P2
  154 + D2=Z
  155 + P2=Q
  156 + Z=DP
  157 +190 Q=PM
  158 + D=S*D+C*P
  159 + P=S*P
  160 + IF(M.EQ.1) GOTO 200
  161 + BI=BI*MM
  162 + BBF=BBF+BI
  163 +200 CONTINUE
  164 +C
  165 + BR=BBR
  166 + BT=BBT
  167 + IF(S.LT.1.E-5) GOTO 210
  168 + BF=BBF/S
  169 + GOTO 211
  170 +210 IF(C.LT.0.) BBF=-BBF
  171 + BF=BBF
  172 +
  173 +211 HE=BR*S+BT*C
  174 + HXGEO=HE*CF-BF*SF
  175 + HYGEO=HE*SF+BF*CF
  176 + HZGEO=BR*C-BT*S
  177 +C
  178 + CALL GEOGSW_08 (HXGEO,HYGEO,HZGEO,HXGSW,HYGSW,HZGSW,1)
  179 +C
  180 + RETURN
  181 + END
  182 +C
  183 +c==========================================================================================
  184 +C
  185 +c
  186 + SUBROUTINE IGRF_GEO_08 (R,THETA,PHI,BR,BTHETA,BPHI)
  187 +c
  188 +C CALCULATES COMPONENTS OF THE MAIN (INTERNAL) GEOMAGNETIC FIELD IN THE SPHERICAL GEOGRAPHIC
  189 +C (GEOCENTRIC) COORDINATE SYSTEM, USING IAGA INTERNATIONAL GEOMAGNETIC REFERENCE MODEL
  190 +C COEFFICIENTS (e.g., http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html, revised: 22 March, 2005)
  191 +C
  192 +C BEFORE THE FIRST CALL OF THIS SUBROUTINE, OR IF THE DATE (IYEAR AND IDAY) WAS CHANGED,
  193 +C THE MODEL COEFFICIENTS SHOULD BE UPDATED BY CALLING THE SUBROUTINE RECALC_08
  194 +C
  195 +C-----INPUT PARAMETERS:
  196 +C
  197 +C R, THETA, PHI - SPHERICAL GEOGRAPHIC (GEOCENTRIC) COORDINATES:
  198 +C RADIAL DISTANCE R IN UNITS RE=6371.2 KM, COLATITUDE THETA AND LONGITUDE PHI IN RADIANS
  199 +C
  200 +C-----OUTPUT PARAMETERS:
  201 +C
  202 +C BR, BTHETA, BPHI - SPHERICAL COMPONENTS OF THE MAIN GEOMAGNETIC FIELD IN NANOTESLA
  203 +C (POSITIVE BR OUTWARD, BTHETA SOUTHWARD, BPHI EASTWARD)
  204 +C
  205 +C LAST MODIFICATION: MAY 4, 2005.
  206 +C THIS VERSION OF THE CODE ACCEPTS DATES FROM 1965 THROUGH 2015.
  207 +c
  208 +C AUTHOR: N. A. TSYGANENKO
  209 +C
  210 +C
  211 + COMMON /GEOPACK2/ G(105),H(105),REC(105)
  212 +
  213 + DIMENSION A(14),B(14)
  214 +
  215 + C=COS(THETA)
  216 + S=SIN(THETA)
  217 + CF=COS(PHI)
  218 + SF=SIN(PHI)
  219 +C
  220 + PP=1./R
  221 + P=PP
  222 +C
  223 +C IN THIS NEW VERSION, THE OPTIMAL VALUE OF THE PARAMETER NM (MAXIMAL ORDER OF THE SPHERICAL
  224 +C HARMONIC EXPANSION) IS NOT USER-PRESCRIBED, BUT CALCULATED INSIDE THE SUBROUTINE, BASED
  225 +C ON THE VALUE OF THE RADIAL DISTANCE R:
  226 +C
  227 + IRP3=R+2
  228 + NM=3+30/IRP3
  229 + IF (NM.GT.13) NM=13
  230 +
  231 + K=NM+1
  232 + DO 150 N=1,K
  233 + P=P*PP
  234 + A(N)=P
  235 +150 B(N)=P*N
  236 +
  237 + P=1.
  238 + D=0.
  239 + BBR=0.
  240 + BBT=0.
  241 + BBF=0.
  242 +
  243 + DO 200 M=1,K
  244 + IF(M.EQ.1) GOTO 160
  245 + MM=M-1
  246 + W=X
  247 + X=W*CF+Y*SF
  248 + Y=Y*CF-W*SF
  249 + GOTO 170
  250 +160 X=0.
  251 + Y=1.
  252 +170 Q=P
  253 + Z=D
  254 + BI=0.
  255 + P2=0.
  256 + D2=0.
  257 + DO 190 N=M,K
  258 + AN=A(N)
  259 + MN=N*(N-1)/2+M
  260 + E=G(MN)
  261 + HH=H(MN)
  262 + W=E*Y+HH*X
  263 + BBR=BBR+B(N)*W*Q
  264 + BBT=BBT-AN*W*Z
  265 + IF(M.EQ.1) GOTO 180
  266 + QQ=Q
  267 + IF(S.LT.1.E-5) QQ=Z
  268 + BI=BI+AN*(E*X-HH*Y)*QQ
  269 +180 XK=REC(MN)
  270 + DP=C*Z-S*Q-XK*D2
  271 + PM=C*Q-XK*P2
  272 + D2=Z
  273 + P2=Q
  274 + Z=DP
  275 +190 Q=PM
  276 + D=S*D+C*P
  277 + P=S*P
  278 + IF(M.EQ.1) GOTO 200
  279 + BI=BI*MM
  280 + BBF=BBF+BI
  281 +200 CONTINUE
  282 +C
  283 + BR=BBR
  284 + BTHETA=BBT
  285 + IF(S.LT.1.E-5) GOTO 210
  286 + BPHI=BBF/S
  287 + RETURN
  288 +210 IF(C.LT.0.) BBF=-BBF
  289 + BPHI=BBF
  290 +
  291 + RETURN
  292 + END
  293 +C
  294 +c==========================================================================================
  295 +c
  296 + SUBROUTINE DIP_08 (XGSW,YGSW,ZGSW,BXGSW,BYGSW,BZGSW)
  297 +C
  298 +C CALCULATES GSW (GEOCENTRIC SOLAR-WIND) COMPONENTS OF GEODIPOLE FIELD WITH THE DIPOLE MOMENT
  299 +C CORRESPONDING TO THE EPOCH, SPECIFIED BY CALLING SUBROUTINE RECALC_08 (SHOULD BE
  300 +C INVOKED BEFORE THE FIRST USE OF THIS ONE, OR IF THE DATE/TIME, AND/OR THE OBSERVED
  301 +C SOLAR WIND DIRECTION, HAVE CHANGED.
  302 +C
  303 +C THE GSW COORDINATE SYSTEM IS ESSENTIALLY SIMILAR TO THE STANDARD GSM (THE TWO SYSTEMS BECOME
  304 +C IDENTICAL TO EACH OTHER IN THE CASE OF STRICTLY RADIAL ANTI-SUNWARD SOLAR WIND FLOW). ITS
  305 +C DETAILED DEFINITION IS GIVEN IN INTRODUCTORY COMMENTS FOR THE SUBROUTINE GSWGSE_08 .
  306 +
  307 +C--INPUT PARAMETERS: XGSW,YGSW,ZGSW - GSW COORDINATES IN RE (1 RE = 6371.2 km)
  308 +C
  309 +C--OUTPUT PARAMETERS: BXGSW,BYGSW,BZGSW - FIELD COMPONENTS IN GSW SYSTEM, IN NANOTESLA.
  310 +C
  311 +C LAST MODIFICATION: JAN 28, 2008.
  312 +C
  313 +C AUTHOR: N. A. TSYGANENKO
  314 +C
  315 + COMMON /GEOPACK1/ AA(10),SPS,CPS,BB(22)
  316 + COMMON /GEOPACK2/ G(105),H(105),REC(105)
  317 +C
  318 + DIPMOM=SQRT(G(2)**2+G(3)**2+H(3)**2)
  319 +C
  320 + P=XGSW**2
  321 + U=ZGSW**2
  322 + V=3.*ZGSW*XGSW
  323 + T=YGSW**2
  324 + Q=DIPMOM/SQRT(P+T+U)**5
  325 + BXGSW=Q*((T+U-2.*P)*SPS-V*CPS)
  326 + BYGSW=-3.*YGSW*Q*(XGSW*SPS+ZGSW*CPS)
  327 + BZGSW=Q*((P+T-2.*U)*CPS-V*SPS)
  328 + RETURN
  329 + END
  330 +
  331 +C*******************************************************************
  332 +c
  333 + SUBROUTINE SUN_08 (IYEAR,IDAY,IHOUR,MIN,ISEC,GST,SLONG,SRASN,SDEC)
  334 +C
  335 +C CALCULATES FOUR QUANTITIES NECESSARY FOR COORDINATE TRANSFORMATIONS
  336 +C WHICH DEPEND ON SUN POSITION (AND, HENCE, ON UNIVERSAL TIME AND SEASON)
  337 +C
  338 +C------- INPUT PARAMETERS:
  339 +C IYR,IDAY,IHOUR,MIN,ISEC - YEAR, DAY, AND UNIVERSAL TIME IN HOURS, MINUTES,
  340 +C AND SECONDS (IDAY=1 CORRESPONDS TO JANUARY 1).
  341 +C
  342 +C------- OUTPUT PARAMETERS:
  343 +C GST - GREENWICH MEAN SIDEREAL TIME, SLONG - LONGITUDE ALONG ECLIPTIC
  344 +C SRASN - RIGHT ASCENSION, SDEC - DECLINATION OF THE SUN (RADIANS)
  345 +C ORIGINAL VERSION OF THIS SUBROUTINE HAS BEEN COMPILED FROM:
  346 +C RUSSELL, C.T., COSMIC ELECTRODYNAMICS, 1971, V.2, PP.184-196.
  347 +C
  348 +C LAST MODIFICATION: MARCH 31, 2003 (ONLY SOME NOTATION CHANGES)
  349 +C
  350 +C ORIGINAL VERSION WRITTEN BY: Gilbert D. Mead
  351 +C
  352 + DOUBLE PRECISION DJ,FDAY
  353 + DATA RAD/57.295779513/
  354 +C
  355 + IF(IYEAR.LT.1901.OR.IYEAR.GT.2099) RETURN
  356 + FDAY=DFLOAT(IHOUR*3600+MIN*60+ISEC)/86400.D0
  357 + DJ=365*(IYEAR-1900)+(IYEAR-1901)/4+IDAY-0.5D0+FDAY
  358 + T=DJ/36525.
  359 + VL=DMOD(279.696678+0.9856473354*DJ,360.D0)
  360 + GST=DMOD(279.690983+.9856473354*DJ+360.*FDAY+180.,360.D0)/RAD
  361 + G=DMOD(358.475845+0.985600267*DJ,360.D0)/RAD
  362 + SLONG=(VL+(1.91946-0.004789*T)*SIN(G)+0.020094*SIN(2.*G))/RAD
  363 + IF(SLONG.GT.6.2831853) SLONG=SLONG-6.2831853
  364 + IF (SLONG.LT.0.) SLONG=SLONG+6.2831853
  365 + OBLIQ=(23.45229-0.0130125*T)/RAD
  366 + SOB=SIN(OBLIQ)
  367 + SLP=SLONG-9.924E-5
  368 +C
  369 +C THE LAST CONSTANT IS A CORRECTION FOR THE ANGULAR ABERRATION DUE TO
  370 +C EARTH'S ORBITAL MOTION
  371 +C
  372 + SIND=SOB*SIN(SLP)
  373 + COSD=SQRT(1.-SIND**2)
  374 + SC=SIND/COSD
  375 + SDEC=ATAN(SC)
  376 + SRASN=3.141592654-ATAN2(COS(OBLIQ)/SOB*SC,-COS(SLP)/COSD)
  377 + RETURN
  378 + END
  379 +C
  380 +C================================================================================
  381 +c
  382 + SUBROUTINE SPHCAR_08 (R,THETA,PHI,X,Y,Z,J)
  383 +C
  384 +C CONVERTS SPHERICAL COORDS INTO CARTESIAN ONES AND VICE VERSA
  385 +C (THETA AND PHI IN RADIANS).
  386 +C
  387 +C J>0 J<0
  388 +C-----INPUT: J,R,THETA,PHI J,X,Y,Z
  389 +C----OUTPUT: X,Y,Z R,THETA,PHI
  390 +C
  391 +C NOTE: AT THE POLES (X=0 AND Y=0) WE ASSUME PHI=0 WHEN CONVERTING
  392 +C FROM CARTESIAN TO SPHERICAL COORDS (I.E., FOR J<0)
  393 +C
  394 +C LAST MOFIFICATION: APRIL 1, 2003 (ONLY SOME NOTATION CHANGES AND MORE
  395 +C COMMENTS ADDED)
  396 +C
  397 +C AUTHOR: N. A. TSYGANENKO
  398 +C
  399 + IF(J.GT.0) GOTO 3
  400 + SQ=X**2+Y**2
  401 + R=SQRT(SQ+Z**2)
  402 + IF (SQ.NE.0.) GOTO 2
  403 + PHI=0.
  404 + IF (Z.LT.0.) GOTO 1
  405 + THETA=0.
  406 + RETURN
  407 + 1 THETA=3.141592654
  408 + RETURN
  409 + 2 SQ=SQRT(SQ)
  410 + PHI=ATAN2(Y,X)
  411 + THETA=ATAN2(SQ,Z)
  412 + IF (PHI.LT.0.) PHI=PHI+6.28318531
  413 + RETURN
  414 + 3 SQ=R*SIN(THETA)
  415 + X=SQ*COS(PHI)
  416 + Y=SQ*SIN(PHI)
  417 + Z=R*COS(THETA)
  418 + RETURN
  419 + END
  420 +C
  421 +C===========================================================================
  422 +c
  423 + SUBROUTINE BSPCAR_08 (THETA,PHI,BR,BTHETA,BPHI,BX,BY,BZ)
  424 +C
  425 +C CALCULATES CARTESIAN FIELD COMPONENTS FROM LOCAL SPHERICAL ONES
  426 +C
  427 +C-----INPUT: THETA,PHI - SPHERICAL ANGLES OF THE POINT IN RADIANS
  428 +C BR,BTHETA,BPHI - LOCAL SPHERICAL COMPONENTS OF THE FIELD
  429 +C-----OUTPUT: BX,BY,BZ - CARTESIAN COMPONENTS OF THE FIELD
  430 +C
  431 +C LAST MOFIFICATION: APRIL 1, 2003 (ONLY SOME NOTATION CHANGES)
  432 +C
  433 +C WRITTEN BY: N. A. TSYGANENKO
  434 +C
  435 + S=SIN(THETA)
  436 + C=COS(THETA)
  437 + SF=SIN(PHI)
  438 + CF=COS(PHI)
  439 + BE=BR*S+BTHETA*C
  440 + BX=BE*CF-BPHI*SF
  441 + BY=BE*SF+BPHI*CF
  442 + BZ=BR*C-BTHETA*S
  443 + RETURN
  444 + END
  445 +c
  446 +C==============================================================================
  447 +C
  448 + SUBROUTINE BCARSP_08 (X,Y,Z,BX,BY,BZ,BR,BTHETA,BPHI)
  449 +C
  450 +CALCULATES LOCAL SPHERICAL FIELD COMPONENTS FROM THOSE IN CARTESIAN SYSTEM
  451 +C
  452 +C-----INPUT: X,Y,Z - CARTESIAN COMPONENTS OF THE POSITION VECTOR
  453 +C BX,BY,BZ - CARTESIAN COMPONENTS OF THE FIELD VECTOR
  454 +C-----OUTPUT: BR,BTHETA,BPHI - LOCAL SPHERICAL COMPONENTS OF THE FIELD VECTOR
  455 +C
  456 +C NOTE: AT THE POLES (THETA=0 OR THETA=PI) WE ASSUME PHI=0,
  457 +C AND HENCE BTHETA=BX, BPHI=BY
  458 +C
  459 +C WRITTEN AND ADDED TO THIS PACKAGE: APRIL 1, 2003,
  460 +C AUTHOR: N. A. TSYGANENKO
  461 +C
  462 + RHO2=X**2+Y**2
  463 + R=SQRT(RHO2+Z**2)
  464 + RHO=SQRT(RHO2)
  465 +
  466 + IF (RHO.NE.0.) THEN
  467 + CPHI=X/RHO
  468 + SPHI=Y/RHO
  469 + ELSE
  470 + CPHI=1.
  471 + SPHI=0.
  472 + ENDIF
  473 +
  474 + CT=Z/R
  475 + ST=RHO/R
  476 +
  477 + BR=(X*BX+Y*BY+Z*BZ)/R
  478 + BTHETA=(BX*CPHI+BY*SPHI)*CT-BZ*ST
  479 + BPHI=BY*CPHI-BX*SPHI
  480 +
  481 + RETURN
  482 + END
  483 +C
  484 +c=====================================================================================
  485 +C
  486 + SUBROUTINE RECALC_08 (IYEAR,IDAY,IHOUR,MIN,ISEC,VGSEX,VGSEY,VGSEZ)
  487 +C
  488 +C 1. PREPARES ELEMENTS OF ROTATION MATRICES FOR TRANSFORMATIONS OF VECTORS BETWEEN
  489 +C SEVERAL COORDINATE SYSTEMS, MOST FREQUENTLY USED IN SPACE PHYSICS.
  490 +C
  491 +C 2. PREPARES COEFFICIENTS USED IN THE CALCULATION OF THE MAIN GEOMAGNETIC FIELD
  492 +C (IGRF MODEL)
  493 +C
  494 +C THIS SUBROUTINE SHOULD BE INVOKED BEFORE USING THE FOLLOWING SUBROUTINES:
  495 +C IGRF_GEO_08, IGRF_GSW_08, DIP_08, GEOMAG_08, GEOGSW_08, MAGSW_08, SMGSW_08, GSWGSE_08,
  496 +c GEIGEO_08, TRACE_08, STEP_08, RHAND_08.
  497 +C
  498 +C THERE IS NO NEED TO REPEATEDLY INVOKE RECALC_08, IF MULTIPLE CALCULATIONS ARE MADE
  499 +C FOR THE SAME DATE/TIME AND SOLAR WIND FLOW DIRECTION.
  500 +C
  501 +C-----INPUT PARAMETERS:
  502 +C
  503 +C IYEAR - YEAR NUMBER (FOUR DIGITS)
  504 +C IDAY - DAY OF YEAR (DAY 1 = JAN 1)
  505 +C IHOUR - HOUR OF DAY (00 TO 23)
  506 +C MIN - MINUTE OF HOUR (00 TO 59)
  507 +C ISEC - SECONDS OF MINUTE (00 TO 59)
  508 +C VGSEX,VGSEY,VGSEZ - GSE (GEOCENTRIC SOLAR-ECLIPTIC) COMPONENTS OF THE OBSERVED
  509 +C SOLAR WIND FLOW VELOCITY (IN KM/S)
  510 +C
  511 +C IMPORTANT: IF ONLY QUESTIONABLE INFORMATION (OR NO INFORMATION AT ALL) IS AVAILABLE
  512 +C ON THE SOLAR WIND SPEED, OR, IF THE STANDARD GSM AND/OR SM COORDINATES ARE
  513 +C INTENDED TO BE USED, THEN SET VGSEX=-400.0 AND VGSEY=VGSEZ=0. IN THIS CASE,
  514 +C THE GSW COORDINATE SYSTEM BECOMES IDENTICAL TO THE STANDARD GSM.
  515 +C
  516 +C IF ONLY SCALAR SPEED V OF THE SOLAR WIND IS KNOWN, THEN SETTING
  517 +C VGSEX=-V, VGSEY=29.78, VGSEZ=0.0 WILL TAKE INTO ACCOUNT THE ~4 degs
  518 +C ABERRATION OF THE MAGNETOSPHERE DUE TO EARTH'S ORBITAL MOTION
  519 +C
  520 +C IF ALL THREE GSE COMPONENTS OF THE SOLAR WIND VELOCITY ARE AVAILABLE,
  521 +C PLEASE NOTE THAT IN SOME SOLAR WIND DATABASES THE ABERRATION EFFECT
  522 +C HAS ALREADY BEEN TAKEN INTO ACCOUNT BY SUBTRACTING 29.78 KM/S FROM VYGSE;
  523 +C IN THAT CASE, THE UNABERRATED (OBSERVED) VYGSE VALUES SHOULD BE RESTORED
  524 +C BY ADDING BACK THE 29.78 KM/S CORRECTION. WHETHER OR NOT TO DO THAT, MUST
  525 +C BE EITHER VERIFIED WITH THE DATA ORIGINATOR OR DETERMINED BY AVERAGING
  526 +C VGSEY OVER A SUFFICIENTLY LONG TIME INTERVAL.
  527 +C
  528 +C-----OUTPUT PARAMETERS: NONE (ALL OUTPUT QUANTITIES ARE PLACED
  529 +C INTO THE COMMON BLOCKS /GEOPACK1/ AND /GEOPACK2/)
  530 +C
  531 +C OTHER SUBROUTINES CALLED BY THIS ONE: SUN_08
  532 +C
  533 +C AUTHOR: N.A. TSYGANENKO
  534 +C DATE: DEC.1, 1991
  535 +C
  536 +C REVISION OF JANUARY 30, 2015:
  537 +C
  538 +C The table of IGRF coefficients was extended to include those for the epoch 2015 (igrf-12)
  539 +c (for details, see http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html)
  540 +C
  541 +C REVISION OF NOVEMBER 30, 2010:
  542 +C
  543 +C The table of IGRF coefficients was extended to include those for the epoch 2010
  544 +c (for details, see http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html)
  545 +C
  546 +C REVISION OF NOVEMBER 15, 2007: ADDED THE POSSIBILITY TO TAKE INTO ACCOUNT THE OBSERVED
  547 +C DEFLECTION OF THE SOLAR WIND FLOW FROM STRICTLY RADIAL DIRECTION. TO THAT END, THREE
  548 +C GSE COMPONENTS OF THE SOLAR WIND VELOCITY WERE ADDED TO THE INPUT PARAMETERS.
  549 +C
  550 +c CORRECTION OF MAY 9, 2006: INTERPOLATION OF THE COEFFICIENTS (BETWEEN
  551 +C LABELS 50 AND 105) IS NOW MADE THROUGH THE LAST ELEMENT OF THE ARRAYS
  552 +C G(105) AND H(105) (PREVIOUSLY MADE ONLY THROUGH N=66, WHICH IN SOME
  553 +C CASES CAUSED RUNTIME ERRORS)
  554 +c
  555 +C REVISION OF MAY 3, 2005:
  556 +C The table of IGRF coefficients was extended to include those for the epoch 2005
  557 +c the maximal order of spherical harmonics was also increased up to 13
  558 +c (for details, see http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html)
  559 +c
  560 +C REVISION OF APRIL 3, 2003:
  561 +c The code now includes preparation of the model coefficients for the subroutines
  562 +c IGRF_08 and GEOMAG_08. This eliminates the need for the SAVE statements, used
  563 +c in the old versions, making the codes easier and more compiler-independent.
  564 +C
  565 + SAVE ISW
  566 +C
  567 + COMMON /GEOPACK1/ ST0,CT0,SL0,CL0,CTCL,STCL,CTSL,STSL,SFI,CFI,
  568 + * SPS,CPS,DS3,CGST,SGST,PSI,A11,A21,A31,A12,A22,A32,A13,A23,A33,
  569 + * E11,E21,E31,E12,E22,E32,E13,E23,E33
  570 +C
  571 +C THE COMMON BLOCK /GEOPACK1/ CONTAINS ELEMENTS OF THE ROTATION MATRICES AND OTHER
  572 +C PARAMETERS RELATED TO THE COORDINATE TRANSFORMATIONS PERFORMED BY THIS PACKAGE
  573 +C
  574 + COMMON /GEOPACK2/ G(105),H(105),REC(105)
  575 +C
  576 +C THE COMMON BLOCK /GEOPACK2/ CONTAINS COEFFICIENTS OF THE IGRF FIELD MODEL, CALCULATED
  577 +C FOR A GIVEN YEAR AND DAY FROM THEIR STANDARD EPOCH VALUES. THE ARRAY REC CONTAINS
  578 +C COEFFICIENTS USED IN THE RECURSION RELATIONS FOR LEGENDRE ASSOCIATE POLYNOMIALS.
  579 +C
  580 + DIMENSION G65(105),H65(105),G70(105),H70(105),G75(105),H75(105),
  581 + + G80(105),H80(105),G85(105),H85(105),G90(105),H90(105),G95(105),
  582 + + H95(105),G00(105),H00(105),G05(105),H05(105),G10(105),H10(105),
  583 + + G15(105),H15(105),DG15(45),DH15(45)
  584 +C
  585 + DATA ISW /0/
  586 +c
  587 + DATA G65/0.,-30334.,-2119.,-1662.,2997.,1594.,1297.,-2038.,1292.,
  588 + *856.,957.,804.,479.,-390.,252.,-219.,358.,254.,-31.,-157.,-62.,
  589 + *45.,61.,8.,-228.,4.,1.,-111.,75.,-57.,4.,13.,-26.,-6.,13.,1.,13.,
  590 + *5.,-4.,-14.,0.,8.,-1.,11.,4.,8.,10.,2.,-13.,10.,-1.,-1.,5.,1.,-2.,
  591 + *-2.,-3.,2.,-5.,-2.,4.,4.,0.,2.,2.,0.,39*0./
  592 + DATA H65/0.,0.,5776.,0.,-2016.,114.,0.,-404.,240.,-165.,0.,148.,
  593 + *-269.,13.,-269.,0.,19.,128.,-126.,-97.,81.,0.,-11.,100.,68.,-32.,
  594 + *-8.,-7.,0.,-61.,-27.,-2.,6.,26.,-23.,-12.,0.,7.,-12.,9.,-16.,4.,
  595 + *24.,-3.,-17.,0.,-22.,15.,7.,-4.,-5.,10.,10.,-4.,1.,0.,2.,1.,2.,
  596 + *6.,-4.,0.,-2.,3.,0.,-6.,39*0./
  597 +c
  598 + DATA G70/0.,-30220.,-2068.,-1781.,3000.,1611.,1287.,-2091.,1278.,
  599 + *838.,952.,800.,461.,-395.,234.,-216.,359.,262.,-42.,-160.,-56.,
  600 + *43.,64.,15.,-212.,2.,3.,-112.,72.,-57.,1.,14.,-22.,-2.,13.,-2.,
  601 + *14.,6.,-2.,-13.,-3.,5.,0.,11.,3.,8.,10.,2.,-12.,10.,-1.,0.,3.,
  602 + *1.,-1.,-3.,-3.,2.,-5.,-1.,6.,4.,1.,0.,3.,-1.,39*0./
  603 + DATA H70/0.,0.,5737.,0.,-2047.,25.,0.,-366.,251.,-196.,0.,167.,
  604 + *-266.,26.,-279.,0.,26.,139.,-139.,-91.,83.,0.,-12.,100.,72.,-37.,
  605 + *-6.,1.,0.,-70.,-27.,-4.,8.,23.,-23.,-11.,0.,7.,-15.,6.,-17.,6.,
  606 + *21.,-6.,-16.,0.,-21.,16.,6.,-4.,-5.,10.,11.,-2.,1.,0.,1.,1.,3.,
  607 + *4.,-4.,0.,-1.,3.,1.,-4.,39*0./
  608 +c
  609 + DATA G75/0.,-30100.,-2013.,-1902.,3010.,1632.,1276.,-2144.,1260.,
  610 + *830.,946.,791.,438.,-405.,216.,-218.,356.,264.,-59.,-159.,-49.,
  611 + *45.,66.,28.,-198.,1.,6.,-111.,71.,-56.,1.,16.,-14.,0.,12.,-5.,
  612 + *14.,6.,-1.,-12.,-8.,4.,0.,10.,1.,7.,10.,2.,-12.,10.,-1.,-1.,4.,
  613 + *1.,-2.,-3.,-3.,2.,-5.,-2.,5.,4.,1.,0.,3.,-1.,39*0./
  614 + DATA H75/0.,0.,5675.,0.,-2067.,-68.,0.,-333.,262.,-223.,0.,191.,
  615 + *-265.,39.,-288.,0.,31.,148.,-152.,-83.,88.,0.,-13.,99.,75.,-41.,
  616 + *-4.,11.,0.,-77.,-26.,-5.,10.,22.,-23.,-12.,0.,6.,-16.,4.,-19.,6.,
  617 + *18.,-10.,-17.,0.,-21.,16.,7.,-4.,-5.,10.,11.,-3.,1.,0.,1.,1.,3.,
  618 + *4.,-4.,-1.,-1.,3.,1.,-5.,39*0./
  619 +c
  620 + DATA G80/0.,-29992.,-1956.,-1997.,3027.,1663.,1281.,-2180.,1251.,
  621 + *833.,938.,782.,398.,-419.,199.,-218.,357.,261.,-74.,-162.,-48.,
  622 + *48.,66.,42.,-192.,4.,14.,-108.,72.,-59.,2.,21.,-12.,1.,11.,-2.,
  623 + *18.,6.,0.,-11.,-7.,4.,3.,6.,-1.,5.,10.,1.,-12.,9.,-3.,-1.,7.,2.,
  624 + *-5.,-4.,-4.,2.,-5.,-2.,5.,3.,1.,2.,3.,0.,39*0./
  625 + DATA H80/0.,0.,5604.,0.,-2129.,-200.,0.,-336.,271.,-252.,0.,212.,
  626 + *-257.,53.,-297.,0.,46.,150.,-151.,-78.,92.,0.,-15.,93.,71.,-43.,
  627 + *-2.,17.,0.,-82.,-27.,-5.,16.,18.,-23.,-10.,0.,7.,-18.,4.,-22.,9.,
  628 + *16.,-13.,-15.,0.,-21.,16.,9.,-5.,-6.,9.,10.,-6.,2.,0.,1.,0.,3.,
  629 + *6.,-4.,0.,-1.,4.,0.,-6.,39*0./
  630 +c
  631 + DATA G85/0.,-29873.,-1905.,-2072.,3044.,1687.,1296.,-2208.,1247.,
  632 + *829.,936.,780.,361.,-424.,170.,-214.,355.,253.,-93.,-164.,-46.,
  633 + *53.,65.,51.,-185.,4.,16.,-102.,74.,-62.,3.,24.,-6.,4.,10.,0.,21.,
  634 + *6.,0.,-11.,-9.,4.,4.,4.,-4.,5.,10.,1.,-12.,9.,-3.,-1.,7.,1.,-5.,
  635 + *-4.,-4.,3.,-5.,-2.,5.,3.,1.,2.,3.,0.,39*0./
  636 + DATA H85/0.,0.,5500.,0.,-2197.,-306.,0.,-310.,284.,-297.,0.,232.,
  637 + *-249.,69.,-297.,0.,47.,150.,-154.,-75.,95.,0.,-16.,88.,69.,-48.,
  638 + *-1.,21.,0.,-83.,-27.,-2.,20.,17.,-23.,-7.,0.,8.,-19.,5.,-23.,11.,
  639 + *14.,-15.,-11.,0.,-21.,15.,9.,-6.,-6.,9.,9.,-7.,2.,0.,1.,0.,3.,
  640 + *6.,-4.,0.,-1.,4.,0.,-6.,39*0./
  641 +c
  642 + DATA G90/0., -29775., -1848., -2131., 3059., 1686., 1314.,
  643 + * -2239., 1248., 802., 939., 780., 325., -423.,
  644 + * 141., -214., 353., 245., -109., -165., -36.,
  645 + * 61., 65., 59., -178., 3., 18., -96.,
  646 + * 77., -64., 2., 26., -1., 5., 9.,
  647 + * 0., 23., 5., -1., -10., -12., 3.,
  648 + * 4., 2., -6., 4., 9., 1., -12.,
  649 + * 9., -4., -2., 7., 1., -6., -3.,
  650 + * -4., 2., -5., -2., 4., 3., 1.,
  651 + * 3., 3., 0., 39*0./
  652 +C
  653 + DATA H90/0., 0., 5406., 0., -2279., -373., 0.,
  654 + * -284., 293., -352., 0., 247., -240., 84.,
  655 + * -299., 0., 46., 154., -153., -69., 97.,
  656 + * 0., -16., 82., 69., -52., 1., 24.,
  657 + * 0., -80., -26., 0., 21., 17., -23.,
  658 + * -4., 0., 10., -19., 6., -22., 12.,
  659 + * 12., -16., -10., 0., -20., 15., 11.,
  660 + * -7., -7., 9., 8., -7., 2., 0.,
  661 + * 2., 1., 3., 6., -4., 0., -2.,
  662 + * 3., -1., -6., 39*0./
  663 +C
  664 + DATA G95/0., -29692., -1784., -2200., 3070., 1681., 1335.,
  665 + * -2267., 1249., 759., 940., 780., 290., -418.,
  666 + * 122., -214., 352., 235., -118., -166., -17.,
  667 + * 68., 67., 68., -170., -1., 19., -93.,
  668 + * 77., -72., 1., 28., 5., 4., 8.,
  669 + * -2., 25., 6., -6., -9., -14., 9.,
  670 + * 6., -5., -7., 4., 9., 3., -10.,
  671 + * 8., -8., -1., 10., -2., -8., -3.,
  672 + * -6., 2., -4., -1., 4., 2., 2.,
  673 + * 5., 1., 0., 39*0./
  674 +C
  675 + DATA H95/0., 0., 5306., 0., -2366., -413., 0.,
  676 + * -262., 302., -427., 0., 262., -236., 97.,
  677 + * -306., 0., 46., 165., -143., -55., 107.,
  678 + * 0., -17., 72., 67., -58., 1., 36.,
  679 + * 0., -69., -25., 4., 24., 17., -24.,
  680 + * -6., 0., 11., -21., 8., -23., 15.,
  681 + * 11., -16., -4., 0., -20., 15., 12.,
  682 + * -6., -8., 8., 5., -8., 3., 0.,
  683 + * 1., 0., 4., 5., -5., -1., -2.,
  684 + * 1., -2., -7., 39*0./
  685 +C
  686 + DATA G00/0.,-29619.4, -1728.2, -2267.7, 3068.4, 1670.9, 1339.6,
  687 + * -2288., 1252.1, 714.5, 932.3, 786.8, 250., -403.,
  688 + * 111.3, -218.8, 351.4, 222.3, -130.4, -168.6, -12.9,
  689 + * 72.3, 68.2, 74.2, -160.9, -5.9, 16.9, -90.4,
  690 + * 79.0, -74.0, 0., 33.3, 9.1, 6.9, 7.3,
  691 + * -1.2, 24.4, 6.6, -9.2, -7.9, -16.6, 9.1,
  692 + * 7.0, -7.9, -7., 5., 9.4, 3., - 8.4,
  693 + * 6.3, -8.9, -1.5, 9.3, -4.3, -8.2, -2.6,
  694 + * -6., 1.7, -3.1, -0.5, 3.7, 1., 2.,
  695 + * 4.2, 0.3, -1.1, 2.7, -1.7, -1.9, 1.5,
  696 + * -0.1, 0.1, -0.7, 0.7, 1.7, 0.1, 1.2,
  697 + * 4.0, -2.2, -0.3, 0.2, 0.9, -0.2, 0.9,
  698 + * -0.5, 0.3, -0.3, -0.4, -0.1, -0.2, -0.4,
  699 + * -0.2, -0.9, 0.3, 0.1, -0.4, 1.3, -0.4,
  700 + * 0.7, -0.4, 0.3, -0.1, 0.4, 0., 0.1/
  701 +C
  702 + DATA H00/0., 0., 5186.1, 0., -2481.6, -458.0, 0.,
  703 + * -227.6, 293.4, -491.1, 0., 272.6, -231.9, 119.8,
  704 + * -303.8, 0., 43.8, 171.9, -133.1, -39.3, 106.3,
  705 + * 0., -17.4, 63.7, 65.1, -61.2, 0.7, 43.8,
  706 + * 0., -64.6, -24.2, 6.2, 24., 14.8, -25.4,
  707 + * -5.8, 0.0, 11.9, -21.5, 8.5, -21.5, 15.5,
  708 + * 8.9, -14.9, -2.1, 0.0, -19.7, 13.4, 12.5,
  709 + * -6.2, -8.4, 8.4, 3.8, -8.2, 4.8, 0.0,
  710 + * 1.7, 0.0, 4.0, 4.9, -5.9, -1.2, -2.9,
  711 + * 0.2, -2.2, -7.4, 0.0, 0.1, 1.3, -0.9,
  712 + * -2.6, 0.9, -0.7, -2.8, -0.9, -1.2, -1.9,
  713 + * -0.9, 0.0, -0.4, 0.3, 2.5, -2.6, 0.7,
  714 + * 0.3, 0.0, 0.0, 0.3, -0.9, -0.4, 0.8,
  715 + * 0.0, -0.9, 0.2, 1.8, -0.4, -1.0, -0.1,
  716 + * 0.7, 0.3, 0.6, 0.3, -0.2, -0.5, -0.9/
  717 +C
  718 + DATA G05/0.,-29554.6, -1669.0, -2337.2, 3047.7, 1657.8, 1336.3,
  719 + * -2305.8, 1246.4, 672.5, 920.6, 798.0, 210.7, -379.9,
  720 + * 100.0, -227.0, 354.4, 208.9, -136.5, -168.1, -13.6,
  721 + * 73.6, 69.6, 76.7, -151.3, -14.6, 14.6, -86.4,
  722 + * 79.9, -74.5, -1.7, 38.7, 12.3, 9.4, 5.4,
  723 + * 1.9, 24.8, 7.6, -11.7, -6.9, -18.1, 10.2,
  724 + * 9.4, -11.3, -4.9, 5.6, 9.8, 3.6, -6.9,
  725 + * 5.0, -10.8, -1.3, 8.8, -6.7, -9.2, -2.2,
  726 + * -6.1, 1.4, -2.4, -0.2, 3.1, 0.3, 2.1,
  727 + * 3.8, -0.2, -2.1, 2.9, -1.6, -1.9, 1.4,
  728 + * -0.3, 0.3, -0.8, 0.5, 1.8, 0.2, 1.0,
  729 + * 4.0, -2.2, -0.3, 0.2, 0.9, -0.4, 1.0,
  730 + * -0.3, 0.5, -0.4, -0.4, 0.1, -0.5, -0.1,
  731 + * -0.2, -0.9, 0.3, 0.3, -0.4, 1.2, -0.4,
  732 + * 0.8, -0.3, 0.4, -0.1, 0.4, -0.1, -0.2/
  733 +C
  734 + DATA H05/0., 0.0, 5078.0, 0.0, -2594.5, -515.4, 0.0,
  735 + * -198.9, 269.7, -524.7, 0.0, 282.1, -225.2, 145.2,
  736 + * -305.4, 0.0, 42.7, 180.3, -123.5, -19.6, 103.9,
  737 + * 0.0, -20.3, 54.8, 63.6, -63.5, 0.2, 50.9,
  738 + * 0.0, -61.1, -22.6, 6.8, 25.4, 10.9, -26.3,
  739 + * -4.6, 0.0, 11.2, -20.9, 9.8, -19.7, 16.2,
  740 + * 7.6, -12.8, -0.1, 0.0, -20.1, 12.7, 12.7,
  741 + * -6.7, -8.2, 8.1, 2.9, -7.7, 6.0, 0.0,
  742 + * 2.2, 0.1, 4.5, 4.8, -6.7, -1.0, -3.5,
  743 + * -0.9, -2.3, -7.9, 0.0, 0.3, 1.4, -0.8,
  744 + * -2.3, 0.9, -0.6, -2.7, -1.1, -1.6, -1.9,
  745 + * -1.4, 0.0, -0.6, 0.2, 2.4, -2.6, 0.6,
  746 + * 0.4, 0.0, 0.0, 0.3, -0.9, -0.3, 0.9,
  747 + * 0.0, -0.8, 0.3, 1.7, -0.5, -1.1, 0.0,
  748 + * 0.6, 0.2, 0.5, 0.4, -0.2, -0.6, -0.9/
  749 +C
  750 + DATA G10/0.00,-29496.57,-1586.42,-2396.06,3026.34,1668.17,1339.85,
  751 + * -2326.54, 1232.10, 633.73, 912.66, 808.97, 166.58,-356.83,
  752 + * 89.40, -230.87, 357.29, 200.26,-141.05,-163.17, -8.03,
  753 + * 72.78, 68.69, 75.92, -141.40, -22.83, 13.10, -78.09,
  754 + * 80.44, -75.00, -4.55, 45.24, 14.00, 10.46, 1.64,
  755 + * 4.92, 24.41, 8.21, -14.50, -5.59, -19.34, 11.61,
  756 + * 10.85, -14.05, -3.54, 5.50, 9.45, 3.45, -5.27,
  757 + * 3.13, -12.38, -0.76, 8.43, -8.42, -10.08, -1.94,
  758 + * -6.24, 0.89, -1.07, -0.16, 2.45, -0.33, 2.13,
  759 + * 3.09, -1.03, -2.80, 3.05, -1.48, -2.03, 1.65,
  760 + * -0.51, 0.54, -0.79, 0.37, 1.79, 0.12, 0.75,
  761 + * 3.75, -2.12, -0.21, 0.30, 1.04, -0.63, 0.95,
  762 + * -0.11, 0.52, -0.39, -0.37, 0.21, -0.77, 0.04,
  763 + * -0.09, -0.89, 0.31, 0.42, -0.45, 1.08, -0.31,
  764 + * 0.78, -0.18, 0.38, 0.02, 0.42, -0.26, -0.26/
  765 +C
  766 + DATA H10/0.00, 0.00, 4944.26, 0.00,-2708.54, -575.73, 0.00,
  767 + * -160.40,251.75, -537.03, 0.00, 286.48, -211.03, 164.46,
  768 + * -309.72, 0.00, 44.58, 189.01, -118.06, -0.01, 101.04,
  769 + * 0.00,-20.90, 44.18, 61.54, -66.26, 3.02, 55.40,
  770 + * 0.00,-57.80, -21.20, 6.54, 24.96, 7.03, -27.61,
  771 + * -3.28, 0.00, 10.84, -20.03, 11.83, -17.41, 16.71,
  772 + * 6.96,-10.74, 1.64, 0.00, -20.54, 11.51, 12.75,
  773 + * -7.14, -7.42, 7.97, 2.14, -6.08, 7.01, 0.00,
  774 + * 2.73, -0.10, 4.71, 4.44, -7.22, -0.96, -3.95,
  775 + * -1.99, -1.97, -8.31, 0.00, 0.13, 1.67, -0.66,
  776 + * -1.76, 0.85, -0.39, -2.51, -1.27, -2.11, -1.94,
  777 + * -1.86, 0.00, -0.87, 0.27, 2.13, -2.49, 0.49,
  778 + * 0.59, 0.00, 0.13, 0.27, -0.86, -0.23, 0.87,
  779 + * 0.00, -0.87, 0.30, 1.66, -0.59, -1.14, -0.07,
  780 + * 0.54, 0.10, 0.49, 0.44, -0.25, -0.53, -0.79/
  781 +C
  782 + DATA G15/0.,-29442.0, -1501.0, -2445.1, 3012.9, 1676.7, 1350.7,
  783 + * -2352.3, 1225.6, 582.0, 907.6, 813.7, 120.4, -334.9,
  784 + * 70.4, -232.6, 360.1, 192.4, -140.9, -157.5, 4.1,
  785 + * 70.0, 67.7, 72.7, -129.9, -28.9, 13.2, -70.9,
  786 + * 81.6, -76.1, -6.8, 51.8, 15.0, 9.4, -2.8,
  787 + * 6.8, 24.2, 8.8, -16.9, -3.2, -20.6, 13.4,
  788 + * 11.7, -15.9, -2.0, 5.4, 8.8, 3.1, -3.3,
  789 + * 0.7, -13.3, -0.1, 8.7, -9.1, -10.5, -1.9,
  790 + * -6.3, 0.1, 0.5, -0.5, 1.8, -0.7, 2.1,
  791 + * 2.4, -1.8, -3.6, 3.1, -1.5, -2.3, 2.0,
  792 + * -0.8, 0.6, -0.7, 0.2, 1.7, -0.2, 0.4,
  793 + * 3.5, -1.9, -0.2, 0.4, 1.2, -0.8, 0.9,
  794 + * 0.1, 0.5, -0.3, -0.4, 0.2, -0.9, 0.0,
  795 + * 0.0, -0.9, 0.4, 0.5, -0.5, 1.0, -0.2,
  796 + * 0.8, -0.1, 0.3, 0.1, 0.5, -0.4, -0.3/
  797 +c
  798 + DATA H15/0.0, 0.0, 4797.1, 0.0, -2845.6, -641.9, 0.0,
  799 + * -115.3, 244.9, -538.4, 0.0, 283.3, -188.7, 180.9,
  800 + * -329.5, 0.0, 47.3, 197.0, -119.3, 16.0, 100.2,
  801 + * 0.0, -20.8, 33.2, 58.9, -66.7, 7.3, 62.6,
  802 + * 0.0, -54.1, -19.5, 5.7, 24.4, 3.4, -27.4,
  803 + * -2.2, 0.0, 10.1, -18.3, 13.3, -14.6, 16.2,
  804 + * 5.7, -9.1, 2.1, 0.0, -21.6, 10.8, 11.8,
  805 + * -6.8, -6.9, 7.8, 1.0, -4.0, 8.4, 0.0,
  806 + * 3.2, -0.4, 4.6, 4.4, -7.9, -0.6, -4.2,
  807 + * -2.8, -1.2, -8.7, 0.0, -0.1, 2.0, -0.7,
  808 + * -1.1, 0.8, -0.2, -2.2, -1.4, -2.5, -2.0,
  809 + * -2.4, 0.0, -1.1, 0.4, 1.9, -2.2, 0.3,
  810 + * 0.7, -0.1, 0.3, 0.2, -0.9, -0.1, 0.7,
  811 + * 0.0, -0.9, 0.4, 1.6, -0.5, -1.2, -0.1,
  812 + * 0.4, -0.1, 0.4, 0.5, -0.3, -0.4, -0.8/
  813 +c
  814 + DATA DG15/0.0, 10.3, 18.1, -8.7, -3.3, 2.1, 3.4,
  815 + * -5.5, -0.7, -10.1, -0.7, 0.2, -9.1, 4.1,
  816 + * -4.3, -0.2, 0.5, -1.3, -0.1, 1.4, 3.9,
  817 + * -0.3, -0.1, -0.7, 2.1, -1.2, 0.3, 1.6,
  818 + * 0.3, -0.2, -0.5, 1.3, 0.1, -0.6, -0.8,
  819 + * 0.2, 0.2, 0.0, -0.6, 0.5, -0.2, 0.4,
  820 + * 0.1, -0.4, 0.3/
  821 +c
  822 + DATA DH15/0.0, 0.0, -26.6, 0.0, -27.4, -14.1, 0.0,
  823 + * 8.2, -0.4, 1.8, 0.0, -1.3, 5.3, 2.9,
  824 + * -5.2, 0.0, 0.6, 1.7, -1.2, 3.4, 0.0,
  825 + * 0.0, 0.0, -2.1, -0.7, 0.2, 0.9, 1.0,
  826 + * 0.0, 0.8, 0.4, -0.2, -0.3, -0.6, 0.1,
  827 + * -0.2, 0.0, -0.3, 0.3, 0.1, 0.5, -0.2,
  828 + * -0.3, 0.3, 0.0/
  829 +C
  830 +c IF (VGSEY.EQ.0..AND.VGSEZ.EQ.0..AND.ISW.NE.1) THEN
  831 +c PRINT *, ''
  832 +c PRINT *,
  833 +c *' RECALC_08: RADIAL SOLAR WIND --> GSW SYSTEM IDENTICAL HERE'
  834 +c PRINT *,
  835 +c *' TO STANDARD GSM (I.E., XGSW AXIS COINCIDES WITH EARTH-SUN LINE)'
  836 +c PRINT *, ''
  837 +c ISW=1
  838 +c ENDIF
  839 +
  840 +c IF ((VGSEY.NE.0..OR.VGSEZ.NE.0.).AND.ISW.NE.2) THEN ! CORRECTED NOV.27, 2009
  841 +c PRINT *, ''
  842 +c PRINT *,
  843 +c *' WARNING: NON-RADIAL SOLAR WIND FLOW SPECIFIED IN RECALC_08;'
  844 +c PRINT *,
  845 +c *' HENCE XGSW AXIS IS ASSUMED ORIENTED ANTIPARALLEL TO V_SW VECTOR'
  846 +c PRINT *, ''
  847 +c ISW=2
  848 +c ENDIF
  849 +C
  850 + IY=IYEAR
  851 +C
  852 +C WE ARE RESTRICTED BY THE INTERVAL 1965-2020, FOR WHICH EITHER THE IGRF/DGRF COEFFICIENTS OR SECULAR VELOCITIES
  853 +c ARE KNOWN; IF IYEAR IS OUTSIDE THIS INTERVAL, THEN THE SUBROUTINE USES THE
  854 +C NEAREST LIMITING VALUE AND PRINTS A WARNING:
  855 +C
  856 + IF(IY.LT.1965) THEN
  857 + IY=1965
  858 + WRITE (*,10) IYEAR,IY
  859 + ENDIF
  860 +
  861 + IF(IY.GT.2020) THEN
  862 + IY=2020
  863 + WRITE (*,10) IYEAR,IY
  864 + ENDIF
  865 +C
  866 +C CALCULATE THE ARRAY REC, CONTAINING COEFFICIENTS FOR THE RECURSION RELATIONS,
  867 +C USED IN THE IGRF SUBROUTINE FOR CALCULATING THE ASSOCIATE LEGENDRE POLYNOMIALS
  868 +C AND THEIR DERIVATIVES:
  869 +c
  870 + DO 20 N=1,14
  871 + N2=2*N-1
  872 + N2=N2*(N2-2)
  873 + DO 20 M=1,N
  874 + MN=N*(N-1)/2+M
  875 +20 REC(MN)=FLOAT((N-M)*(N+M-2))/FLOAT(N2)
  876 +C
  877 + IF (IY.LT.1970) GOTO 50 !INTERPOLATE BETWEEN 1965 - 1970
  878 + IF (IY.LT.1975) GOTO 60 !INTERPOLATE BETWEEN 1970 - 1975
  879 + IF (IY.LT.1980) GOTO 70 !INTERPOLATE BETWEEN 1975 - 1980
  880 + IF (IY.LT.1985) GOTO 80 !INTERPOLATE BETWEEN 1980 - 1985
  881 + IF (IY.LT.1990) GOTO 90 !INTERPOLATE BETWEEN 1985 - 1990
  882 + IF (IY.LT.1995) GOTO 100 !INTERPOLATE BETWEEN 1990 - 1995
  883 + IF (IY.LT.2000) GOTO 110 !INTERPOLATE BETWEEN 1995 - 2000
  884 + IF (IY.LT.2005) GOTO 120 !INTERPOLATE BETWEEN 2000 - 2005
  885 + IF (IY.LT.2010) GOTO 130 !INTERPOLATE BETWEEN 2005 - 2010
  886 + IF (IY.LT.2015) GOTO 140 !INTERPOLATE BETWEEN 2010 - 2015
  887 +C
  888 +C EXTRAPOLATE BEYOND 2015:
  889 +C
  890 + DT=FLOAT(IY)+FLOAT(IDAY-1)/365.25-2015.
  891 + DO 40 N=1,105
  892 + G(N)=G15(N)
  893 + H(N)=H15(N)
  894 + IF (N.GT.45) GOTO 40
  895 + G(N)=G(N)+DG15(N)*DT
  896 + H(N)=H(N)+DH15(N)*DT
  897 +40 CONTINUE
  898 + GOTO 300
  899 +C
  900 +C INTERPOLATE BETWEEEN 1965 - 1970:
  901 +C
  902 +50 F2=(FLOAT(IY)+FLOAT(IDAY-1)/365.25-1965)/5.
  903 + F1=1.-F2
  904 + DO 55 N=1,105
  905 + G(N)=G65(N)*F1+G70(N)*F2
  906 +55 H(N)=H65(N)*F1+H70(N)*F2
  907 + GOTO 300
  908 +C
  909 +C INTERPOLATE BETWEEN 1970 - 1975:
  910 +C
  911 +60 F2=(FLOAT(IY)+FLOAT(IDAY-1)/365.25-1970)/5.
  912 + F1=1.-F2
  913 + DO 65 N=1,105
  914 + G(N)=G70(N)*F1+G75(N)*F2
  915 +65 H(N)=H70(N)*F1+H75(N)*F2
  916 + GOTO 300
  917 +C
  918 +C INTERPOLATE BETWEEN 1975 - 1980:
  919 +C
  920 +70 F2=(FLOAT(IY)+FLOAT(IDAY-1)/365.25-1975)/5.
  921 + F1=1.-F2
  922 + DO 75 N=1,105
  923 + G(N)=G75(N)*F1+G80(N)*F2
  924 +75 H(N)=H75(N)*F1+H80(N)*F2
  925 + GOTO 300
  926 +C
  927 +C INTERPOLATE BETWEEN 1980 - 1985:
  928 +C
  929 +80 F2=(FLOAT(IY)+FLOAT(IDAY-1)/365.25-1980)/5.
  930 + F1=1.-F2
  931 + DO 85 N=1,105
  932 + G(N)=G80(N)*F1+G85(N)*F2
  933 +85 H(N)=H80(N)*F1+H85(N)*F2
  934 + GOTO 300
  935 +C
  936 +C INTERPOLATE BETWEEN 1985 - 1990:
  937 +C
  938 +90 F2=(FLOAT(IY)+FLOAT(IDAY-1)/365.25-1985)/5.
  939 + F1=1.-F2
  940 + DO 95 N=1,105
  941 + G(N)=G85(N)*F1+G90(N)*F2
  942 +95 H(N)=H85(N)*F1+H90(N)*F2
  943 + GOTO 300
  944 +C
  945 +C INTERPOLATE BETWEEN 1990 - 1995:
  946 +C
  947 +100 F2=(FLOAT(IY)+FLOAT(IDAY-1)/365.25-1990)/5.
  948 + F1=1.-F2
  949 + DO 105 N=1,105
  950 + G(N)=G90(N)*F1+G95(N)*F2
  951 +105 H(N)=H90(N)*F1+H95(N)*F2
  952 + GOTO 300
  953 +C
  954 +C INTERPOLATE BETWEEN 1995 - 2000:
  955 +C
  956 +110 F2=(FLOAT(IY)+FLOAT(IDAY-1)/365.25-1995)/5.
  957 + F1=1.-F2
  958 + DO 115 N=1,105 ! THE 2000 COEFFICIENTS (G00) GO THROUGH THE ORDER 13, NOT 10
  959 + G(N)=G95(N)*F1+G00(N)*F2
  960 +115 H(N)=H95(N)*F1+H00(N)*F2
  961 + GOTO 300
  962 +C
  963 +C INTERPOLATE BETWEEN 2000 - 2005:
  964 +C
  965 +120 F2=(FLOAT(IY)+FLOAT(IDAY-1)/365.25-2000)/5.
  966 + F1=1.-F2
  967 + DO 125 N=1,105
  968 + G(N)=G00(N)*F1+G05(N)*F2
  969 +125 H(N)=H00(N)*F1+H05(N)*F2
  970 + GOTO 300
  971 +C
  972 +C INTERPOLATE BETWEEN 2005 - 2010:
  973 +C
  974 +130 F2=(FLOAT(IY)+FLOAT(IDAY-1)/365.25-2005)/5.
  975 + F1=1.-F2
  976 + DO 135 N=1,105
  977 + G(N)=G05(N)*F1+G10(N)*F2
  978 +135 H(N)=H05(N)*F1+H10(N)*F2
  979 + GOTO 300
  980 +C
  981 +C INTERPOLATE BETWEEN 2010 - 2015:
  982 +C
  983 +140 F2=(FLOAT(IY)+FLOAT(IDAY-1)/365.25-2010)/5.
  984 + F1=1.-F2
  985 + DO 145 N=1,105
  986 + G(N)=G10(N)*F1+G15(N)*F2
  987 +145 H(N)=H10(N)*F1+H15(N)*F2
  988 + GOTO 300
  989 +C
  990 +C COEFFICIENTS FOR A GIVEN YEAR HAVE BEEN CALCULATED; NOW MULTIPLY
  991 +C THEM BY SCHMIDT NORMALIZATION FACTORS:
  992 +C
  993 +300 S=1.
  994 + DO 150 N=2,14
  995 + MN=N*(N-1)/2+1
  996 + S=S*FLOAT(2*N-3)/FLOAT(N-1)
  997 + G(MN)=G(MN)*S
  998 + H(MN)=H(MN)*S
  999 + P=S
  1000 + DO 150 M=2,N
  1001 + AA=1.
  1002 + IF (M.EQ.2) AA=2.
  1003 + P=P*SQRT(AA*FLOAT(N-M+1)/FLOAT(N+M-2))
  1004 + MNN=MN+M-1
  1005 + G(MNN)=G(MNN)*P
  1006 +150 H(MNN)=H(MNN)*P
  1007 +
  1008 + G_10=-G(2)
  1009 + G_11= G(3)
  1010 + H_11= H(3)
  1011 +C
  1012 +C NOW CALCULATE GEO COMPONENTS OF THE UNIT VECTOR EzMAG, PARALLEL TO GEODIPOLE AXIS:
  1013 +C SIN(TETA0)*COS(LAMBDA0), SIN(TETA0)*SIN(LAMBDA0), AND COS(TETA0)
  1014 +C ST0 * CL0 ST0 * SL0 CT0
  1015 +C
  1016 + SQ=G_11**2+H_11**2
  1017 + SQQ=SQRT(SQ)
  1018 + SQR=SQRT(G_10**2+SQ)
  1019 + SL0=-H_11/SQQ
  1020 + CL0=-G_11/SQQ
  1021 + ST0=SQQ/SQR
  1022 + CT0=G_10/SQR
  1023 + STCL=ST0*CL0
  1024 + STSL=ST0*SL0
  1025 + CTSL=CT0*SL0
  1026 + CTCL=CT0*CL0
  1027 +C
  1028 +C NOW CALCULATE GEI COMPONENTS (S1,S2,S3) OF THE UNIT VECTOR S = EX_GSE
  1029 +C POINTING FROM THE EARTH'S CENTER TO SUN
  1030 +C
  1031 + CALL SUN_08 (IY,IDAY,IHOUR,MIN,ISEC,GST,SLONG,SRASN,SDEC)
  1032 +C
  1033 + S1=COS(SRASN)*COS(SDEC)
  1034 + S2=SIN(SRASN)*COS(SDEC)
  1035 + S3=SIN(SDEC)
  1036 +C
  1037 +C NOW CALCULATE GEI COMPONENTS (DZ1,DZ2,DZ3) OF THE UNIT VECTOR EZGSE
  1038 +C POINTING NORTHWARD AND ORTHOGONAL TO THE ECLIPTIC PLANE, AS
  1039 +C (0,-SIN(OBLIQ),COS(OBLIQ)). FOR THE EPOCH 1978, OBLIQ = 23.44214 DEGS.
  1040 +C HERE WE USE A MORE ACCURATE TIME-DEPENDENT VALUE, DETERMINED AS:
  1041 +C
  1042 + DJ=FLOAT(365*(IY-1900)+(IY-1901)/4 +IDAY)
  1043 + * -0.5+FLOAT(IHOUR*3600+MIN*60+ISEC)/86400.
  1044 + T=DJ/36525.
  1045 + OBLIQ=(23.45229-0.0130125*T)/57.2957795
  1046 + DZ1=0.
  1047 + DZ2=-SIN(OBLIQ)
  1048 + DZ3=COS(OBLIQ)
  1049 +C
  1050 +C NOW WE OBTAIN GEI COMPONENTS OF THE UNIT VECTOR EYGSE=(DY1,DY2,DY3),
  1051 +C COMPLETING THE RIGHT-HANDED SYSTEM. THEY CAN BE FOUND FROM THE VECTOR
  1052 +C PRODUCT EZGSE x EXGSE = (DZ1,DZ2,DZ3) x (S1,S2,S3):
  1053 +C
  1054 + DY1=DZ2*S3-DZ3*S2
  1055 + DY2=DZ3*S1-DZ1*S3
  1056 + DY3=DZ1*S2-DZ2*S1
  1057 +C
  1058 +C NOW LET'S CALCULATE GEI COMPONENTS OF THE UNIT VECTOR X = EXGSW, DIRECTED ANTIPARALLEL
  1059 +C TO THE OBSERVED SOLAR WIND FLOW. FIRST, CALCULATE ITS COMPONENTS IN GSE:
  1060 +C
  1061 + V=SQRT(VGSEX**2+VGSEY**2+VGSEZ**2)
  1062 + DX1=-VGSEX/V
  1063 + DX2=-VGSEY/V
  1064 + DX3=-VGSEZ/V
  1065 +C
  1066 +C THEN IN GEI:
  1067 +C
  1068 + X1=DX1*S1+DX2*DY1+DX3*DZ1
  1069 + X2=DX1*S2+DX2*DY2+DX3*DZ2
  1070 + X3=DX1*S3+DX2*DY3+DX3*DZ3
  1071 +C
  1072 +C NOW CALCULATE GEI COMPONENTS (DIP1,DIP2,DIP3) OF THE UNIT VECTOR DIP = EZ_SM = EZ_MAG,
  1073 +C ALIGNED WITH THE GEODIPOLE AND POINTING NORTHWARD FROM ECLIPTIC PLANE:
  1074 +C
  1075 + CGST=COS(GST)
  1076 + SGST=SIN(GST)
  1077 +C
  1078 + DIP1=STCL*CGST-STSL*SGST
  1079 + DIP2=STCL*SGST+STSL*CGST
  1080 + DIP3=CT0
  1081 +C
  1082 +C THIS ALLOWS US TO CALCULATE GEI COMPONENTS OF THE UNIT VECTOR Y = EYGSW
  1083 +C BY TAKING THE VECTOR PRODUCT DIP x X AND NORMALIZING IT TO UNIT LENGTH:
  1084 +C
  1085 + Y1=DIP2*X3-DIP3*X2
  1086 + Y2=DIP3*X1-DIP1*X3
  1087 + Y3=DIP1*X2-DIP2*X1
  1088 + Y=SQRT(Y1*Y1+Y2*Y2+Y3*Y3)
  1089 + Y1=Y1/Y
  1090 + Y2=Y2/Y
  1091 + Y3=Y3/Y
  1092 +C
  1093 +C AND GEI COMPONENTS OF THE UNIT VECTOR Z = EZGSW = EXGSW x EYGSW = X x Y:
  1094 +C
  1095 + Z1=X2*Y3-X3*Y2
  1096 + Z2=X3*Y1-X1*Y3
  1097 + Z3=X1*Y2-X2*Y1
  1098 +C
  1099 +C ELEMENTS OF THE MATRIX GSE TO GSW ARE THE SCALAR PRODUCTS:
  1100 +C
  1101 +C E11=(EXGSE,EXGSW) E12=(EXGSE,EYGSW) E13=(EXGSE,EZGSW)
  1102 +C E21=(EYGSE,EXGSW) E22=(EYGSE,EYGSW) E23=(EYGSE,EZGSW)
  1103 +C E31=(EZGSE,EXGSW) E32=(EZGSE,EYGSW) E33=(EZGSE,EZGSW)
  1104 +C
  1105 + E11= S1*X1 +S2*X2 +S3*X3
  1106 + E12= S1*Y1 +S2*Y2 +S3*Y3
  1107 + E13= S1*Z1 +S2*Z2 +S3*Z3
  1108 + E21=DY1*X1+DY2*X2+DY3*X3
  1109 + E22=DY1*Y1+DY2*Y2+DY3*Y3
  1110 + E23=DY1*Z1+DY2*Z2+DY3*Z3
  1111 + E31=DZ1*X1+DZ2*X2+DZ3*X3
  1112 + E32=DZ1*Y1+DZ2*Y2+DZ3*Y3
  1113 + E33=DZ1*Z1+DZ2*Z2+DZ3*Z3
  1114 +C
  1115 +C GEODIPOLE TILT ANGLE IN THE GSW SYSTEM: PSI=ARCSIN(DIP,EXGSW)
  1116 +C
  1117 + SPS=DIP1*X1+DIP2*X2+DIP3*X3
  1118 + CPS=SQRT(1.-SPS**2)
  1119 + PSI=ASIN(SPS)
  1120 +C
  1121 +C ELEMENTS OF THE MATRIX GEO TO GSW ARE THE SCALAR PRODUCTS:
  1122 +C
  1123 +C A11=(EXGEO,EXGSW), A12=(EYGEO,EXGSW), A13=(EZGEO,EXGSW),
  1124 +C A21=(EXGEO,EYGSW), A22=(EYGEO,EYGSW), A23=(EZGEO,EYGSW),
  1125 +C A31=(EXGEO,EZGSW), A32=(EYGEO,EZGSW), A33=(EZGEO,EZGSW),
  1126 +C
  1127 +C ALL THE UNIT VECTORS IN BRACKETS ARE ALREADY DEFINED IN GEI:
  1128 +C
  1129 +C EXGEO=(CGST,SGST,0), EYGEO=(-SGST,CGST,0), EZGEO=(0,0,1)
  1130 +C EXGSW=(X1,X2,X3), EYGSW=(Y1,Y2,Y3), EZGSW=(Z1,Z2,Z3)
  1131 +C AND THEREFORE:
  1132 +C
  1133 + A11=X1*CGST+X2*SGST
  1134 + A12=-X1*SGST+X2*CGST
  1135 + A13=X3
  1136 + A21=Y1*CGST+Y2*SGST
  1137 + A22=-Y1*SGST+Y2*CGST
  1138 + A23=Y3
  1139 + A31=Z1*CGST+Z2*SGST
  1140 + A32=-Z1*SGST+Z2*CGST
  1141 + A33=Z3
  1142 +C
  1143 +C NOW CALCULATE ELEMENTS OF THE MATRIX MAG TO SM (ONE ROTATION ABOUT THE GEODIPOLE AXIS);
  1144 +C THEY ARE FOUND AS THE SCALAR PRODUCTS: CFI=GM22=(EYSM,EYMAG)=(EYGSW,EYMAG),
  1145 +C SFI=GM23=(EYSM,EXMAG)=(EYGSW,EXMAG),
  1146 +C DERIVED AS FOLLOWS:
  1147 +C
  1148 +C IN GEO, THE VECTORS EXMAG AND EYMAG HAVE THE COMPONENTS (CT0*CL0,CT0*SL0,-ST0)
  1149 +C AND (-SL0,CL0,0), RESPECTIVELY. HENCE, IN GEI THEIR COMPONENTS ARE:
  1150 +C EXMAG: CT0*CL0*COS(GST)-CT0*SL0*SIN(GST)
  1151 +C CT0*CL0*SIN(GST)+CT0*SL0*COS(GST)
  1152 +C -ST0
  1153 +C EYMAG: -SL0*COS(GST)-CL0*SIN(GST)
  1154 +C -SL0*SIN(GST)+CL0*COS(GST)
  1155 +C 0
  1156 +C NOW, NOTE THAT GEI COMPONENTS OF EYSM=EYGSW WERE FOUND ABOVE AS Y1, Y2, AND Y3,
  1157 +C AND WE ONLY HAVE TO COMBINE THESE QUANTITIES INTO SCALAR PRODUCTS:
  1158 +C
  1159 + EXMAGX=CT0*(CL0*CGST-SL0*SGST)
  1160 + EXMAGY=CT0*(CL0*SGST+SL0*CGST)
  1161 + EXMAGZ=-ST0
  1162 + EYMAGX=-(SL0*CGST+CL0*SGST)
  1163 + EYMAGY=-(SL0*SGST-CL0*CGST)
  1164 + CFI=Y1*EYMAGX+Y2*EYMAGY
  1165 + SFI=Y1*EXMAGX+Y2*EXMAGY+Y3*EXMAGZ
  1166 +C
  1167 + 10 FORMAT(//1X,
  1168 + *'**** RECALC_08 WARNS: YEAR IS OUT OF INTERVAL 1965-2020: IYEAR=',
  1169 + *I4,/,6X,'CALCULATIONS WILL BE DONE FOR IYEAR=',I4,/)
  1170 + RETURN
  1171 + END
  1172 +c
  1173 +c==================================================================================
  1174 +
  1175 + SUBROUTINE GSWGSE_08 (XGSW,YGSW,ZGSW,XGSE,YGSE,ZGSE,J)
  1176 +C
  1177 +C THIS SUBROUTINE TRANSFORMS COMPONENTS OF ANY VECTOR BETWEEN THE STANDARD GSE
  1178 +C COORDINATE SYSTEM AND THE GEOCENTRIC SOLAR-WIND (GSW, aka GSWM), DEFINED AS FOLLOWS
  1179 +C (HONES ET AL., PLANET.SPACE SCI., V.34, P.889, 1986; TSYGANENKO ET AL., JGRA,
  1180 +C V.103(A4), P.6827, 1998):
  1181 +C
  1182 +C IN THE GSW SYSTEM, X AXIS IS ANTIPARALLEL TO THE OBSERVED DIRECTION OF THE SOLAR WIND FLOW.
  1183 +C TWO OTHER AXES, Y AND Z, ARE DEFINED IN THE SAME WAY AS FOR THE STANDARD GSM, THAT IS,
  1184 +C Z AXIS ORTHOGONAL TO X AXIS, POINTS NORTHWARD, AND LIES IN THE PLANE DEFINED BY THE X-
  1185 +C AND GEODIPOLE AXIS. THE Y AXIS COMPLETES THE RIGHT-HANDED SYSTEM.
  1186 +C
  1187 +C THE GSW SYSTEM BECOMES IDENTICAL TO THE STANDARD GSM IN THE CASE OF
  1188 +C A STRICTLY RADIAL SOLAR WIND FLOW.
  1189 +C
  1190 +C AUTHOR: N. A. TSYGANENKO
  1191 +C ADDED TO 2008 VERSION OF GEOPACK: JAN 27, 2008.
  1192 +C
  1193 +C J>0 J<0
  1194 +C-----INPUT: J,XGSW,YGSW,ZGSW J,XGSE,YGSE,ZGSE
  1195 +C-----OUTPUT: XGSE,YGSE,ZGSE XGSW,YGSW,ZGSW
  1196 +C
  1197 +C IMPORTANT THINGS TO REMEMBER:
  1198 +C
  1199 +C (1) BEFORE CALLING GSWGSE_08, BE SURE TO INVOKE SUBROUTINE RECALC_08, IN ORDER
  1200 +C TO DEFINE ALL NECESSARY ELEMENTS OF TRANSFORMATION MATRICES
  1201 +C
  1202 +C (2) IN THE ABSENCE OF INFORMATION ON THE SOLAR WIND DIRECTION, E.G., WITH ONLY SCALAR
  1203 +C SPEED V KNOWN, THIS SUBROUTINE CAN BE USED TO CONVERT VECTORS TO ABERRATED
  1204 +C COORDINATE SYSTEM, TAKING INTO ACCOUNT EARTH'S ORBITAL SPEED OF 29 KM/S.
  1205 +C TO DO THAT, SPECIFY THE LAST 3 PARAMETERS IN RECALC_08 AS FOLLOWS:
  1206 +C VGSEX=-V, VGSEY=29.0, VGSEZ=0.0.
  1207 +C
  1208 +C IT SHOULD ALSO BE KEPT IN MIND THAT IN SOME SOLAR WIND DATABASES THE ABERRATION
  1209 +C EFFECT HAS ALREADY BEEN TAKEN INTO ACCOUNT BY SUBTRACTING 29 KM/S FROM VYGSE;
  1210 +C IN THAT CASE, THE ORIGINAL VYGSE VALUES SHOULD BE RESTORED BY ADDING BACK THE
  1211 +C 29 KM/S CORRECTION. WHETHER OR NOT TO DO THAT, MUST BE VERIFIED WITH THE DATA
  1212 +C ORIGINATOR (OR CAN BE DETERMINED BY CALCULATING THE AVERAGE VGSEY OVER
  1213 +C A SUFFICIENTLY LONG TIME INTERVAL)
  1214 +C
  1215 +C (3) IF NO INFORMATION IS AVAILABLE ON THE SOLAR WIND SPEED, THEN SET VGSEX=-400.0
  1216 +c AND VGSEY=VGSEZ=0. IN THAT CASE, THE GSW COORDINATE SYSTEM BECOMES
  1217 +c IDENTICAL TO THE STANDARD ONE.
  1218 +C
  1219 + COMMON /GEOPACK1/ AAA(25),E11,E21,E31,E12,E22,E32,E13,E23,E33
  1220 +C
  1221 +C DIRECT TRANSFORMATION:
  1222 +C
  1223 + IF (J.GT.0) THEN
  1224 + XGSE=XGSW*E11+YGSW*E12+ZGSW*E13
  1225 + YGSE=XGSW*E21+YGSW*E22+ZGSW*E23
  1226 + ZGSE=XGSW*E31+YGSW*E32+ZGSW*E33
  1227 + ENDIF
  1228 +C
  1229 +C INVERSE TRANSFORMATION: CARRIED OUT USING THE TRANSPOSED MATRIX:
  1230 +C
  1231 + IF (J.LT.0) THEN
  1232 + XGSW=XGSE*E11+YGSE*E21+ZGSE*E31
  1233 + YGSW=XGSE*E12+YGSE*E22+ZGSE*E32
  1234 + ZGSW=XGSE*E13+YGSE*E23+ZGSE*E33
  1235 + ENDIF
  1236 +C
  1237 + RETURN
  1238 + END
  1239 +C
  1240 +C========================================================================================
  1241 +C
  1242 + SUBROUTINE GEOMAG_08 (XGEO,YGEO,ZGEO,XMAG,YMAG,ZMAG,J)
  1243 +C
  1244 +C CONVERTS GEOGRAPHIC (GEO) TO DIPOLE (MAG) COORDINATES OR VICE VERSA.
  1245 +C
  1246 +C J>0 J<0
  1247 +C-----INPUT: J,XGEO,YGEO,ZGEO J,XMAG,YMAG,ZMAG
  1248 +C-----OUTPUT: XMAG,YMAG,ZMAG XGEO,YGEO,ZGEO
  1249 +C
  1250 +C ATTENTION: SUBROUTINE RECALC_08 MUST BE INVOKED BEFORE GEOMAG_08 IN TWO CASES:
  1251 +C /A/ BEFORE THE FIRST TRANSFORMATION OF COORDINATES
  1252 +C /B/ IF THE VALUES OF IYEAR AND/OR IDAY HAVE BEEN CHANGED
  1253 +C
  1254 +C NO INFORMATION IS REQUIRED HERE ON THE SOLAR WIND VELOCITY, SO ONE
  1255 +C CAN SET VGSEX=-400.0, VGSEY=0.0, VGSEZ=0.0 IN RECALC_08.
  1256 +C
  1257 +C LAST MOFIFICATION: JAN 28, 2008
  1258 +C
  1259 +C AUTHOR: N. A. TSYGANENKO
  1260 +C
  1261 + COMMON /GEOPACK1/ ST0,CT0,SL0,CL0,CTCL,STCL,CTSL,STSL,AB(26)
  1262 +
  1263 + IF(J.GT.0) THEN
  1264 + XMAG=XGEO*CTCL+YGEO*CTSL-ZGEO*ST0
  1265 + YMAG=YGEO*CL0-XGEO*SL0
  1266 + ZMAG=XGEO*STCL+YGEO*STSL+ZGEO*CT0
  1267 + ELSE
  1268 + XGEO=XMAG*CTCL-YMAG*SL0+ZMAG*STCL
  1269 + YGEO=XMAG*CTSL+YMAG*CL0+ZMAG*STSL
  1270 + ZGEO=ZMAG*CT0-XMAG*ST0
  1271 + ENDIF
  1272 +
  1273 + RETURN
  1274 + END
  1275 +c
  1276 +c=========================================================================================
  1277 +c
  1278 + SUBROUTINE GEIGEO_08 (XGEI,YGEI,ZGEI,XGEO,YGEO,ZGEO,J)
  1279 +C
  1280 +C CONVERTS EQUATORIAL INERTIAL (GEI) TO GEOGRAPHICAL (GEO) COORDS
  1281 +C OR VICE VERSA.
  1282 +C J>0 J<0
  1283 +C----INPUT: J,XGEI,YGEI,ZGEI J,XGEO,YGEO,ZGEO
  1284 +C----OUTPUT: XGEO,YGEO,ZGEO XGEI,YGEI,ZGEI
  1285 +C
  1286 +C ATTENTION: SUBROUTINE RECALC_08 MUST BE INVOKED BEFORE GEIGEO_08 IN TWO CASES:
  1287 +C /A/ BEFORE THE FIRST TRANSFORMATION OF COORDINATES
  1288 +C /B/ IF THE CURRENT VALUES OF IYEAR,IDAY,IHOUR,MIN,ISEC HAVE BEEN CHANGED
  1289 +C
  1290 +C NO INFORMATION IS REQUIRED HERE ON THE SOLAR WIND VELOCITY, SO ONE
  1291 +C CAN SET VGSEX=-400.0, VGSEY=0.0, VGSEZ=0.0 IN RECALC_08.
  1292 +C
  1293 +C LAST MODIFICATION: JAN 28, 2008
  1294 +
  1295 +C AUTHOR: N. A. TSYGANENKO
  1296 +C
  1297 + COMMON /GEOPACK1/ A(13),CGST,SGST,B(19)
  1298 +C
  1299 + IF(J.GT.0) THEN
  1300 + XGEO=XGEI*CGST+YGEI*SGST
  1301 + YGEO=YGEI*CGST-XGEI*SGST
  1302 + ZGEO=ZGEI
  1303 + ELSE
  1304 + XGEI=XGEO*CGST-YGEO*SGST
  1305 + YGEI=YGEO*CGST+XGEO*SGST
  1306 + ZGEI=ZGEO
  1307 + ENDIF
  1308 +
  1309 + RETURN
  1310 + END
  1311 +C
  1312 +C=======================================================================================
  1313 +C
  1314 + SUBROUTINE MAGSM_08 (XMAG,YMAG,ZMAG,XSM,YSM,ZSM,J)
  1315 +C
  1316 +C CONVERTS DIPOLE (MAG) TO SOLAR MAGNETIC (SM) COORDINATES OR VICE VERSA
  1317 +C
  1318 +C J>0 J<0
  1319 +C----INPUT: J,XMAG,YMAG,ZMAG J,XSM,YSM,ZSM
  1320 +C----OUTPUT: XSM,YSM,ZSM XMAG,YMAG,ZMAG
  1321 +C
  1322 +C ATTENTION: SUBROUTINE RECALC_08 MUST BE INVOKED BEFORE MAGSM_08 IN THREE CASES:
  1323 +C /A/ BEFORE THE FIRST TRANSFORMATION OF COORDINATES, OR
  1324 +C /B/ IF THE VALUES OF IYEAR,IDAY,IHOUR,MIN,ISEC HAVE CHANGED, AND/OR
  1325 +C /C/ IF THE VALUES OF COMPONENTS OF THE SOLAR WIND FLOW VELOCITY HAVE CHANGED
  1326 +C
  1327 +C IMPORTANT NOTE:
  1328 +C
  1329 +C A NON-STANDARD DEFINITION IS IMPLIED HERE FOR THE SOLAR MAGNETIC COORDINATE
  1330 +C SYSTEM: IT IS ASSUMED THAT THE XSM AXIS LIES IN THE PLANE DEFINED BY THE
  1331 +C GEODIPOLE AXIS AND THE OBSERVED VECTOR OF THE SOLAR WIND FLOW (RATHER THAN
  1332 +C THE EARTH-SUN LINE). IN ORDER TO CONVERT MAG COORDINATES TO AND FROM THE
  1333 +C STANDARD SM COORDINATES, INVOKE RECALC_08 WITH VGSEX=-400.0, VGSEY=0.0, VGSEZ=0.0
  1334 +C
  1335 +C LAST MODIFICATION: FEB 07, 2008
  1336 +C
  1337 +C AUTHOR: N. A. TSYGANENKO
  1338 +C
  1339 + COMMON /GEOPACK1/ A(8),SFI,CFI,B(24)
  1340 +C
  1341 + IF (J.GT.0) THEN
  1342 + XSM=XMAG*CFI-YMAG*SFI
  1343 + YSM=XMAG*SFI+YMAG*CFI
  1344 + ZSM=ZMAG
  1345 + ELSE
  1346 + XMAG=XSM*CFI+YSM*SFI
  1347 + YMAG=YSM*CFI-XSM*SFI
  1348 + ZMAG=ZSM
  1349 + ENDIF
  1350 +
  1351 + RETURN
  1352 + END
  1353 +C
  1354 +C=====================================================================================
  1355 +C
  1356 + SUBROUTINE SMGSW_08 (XSM,YSM,ZSM,XGSW,YGSW,ZGSW,J)
  1357 +C
  1358 +C CONVERTS SOLAR MAGNETIC (SM) TO GEOCENTRIC SOLAR-WIND (GSW) COORDINATES OR VICE VERSA.
  1359 +C
  1360 +C J>0 J<0
  1361 +C-----INPUT: J,XSM,YSM,ZSM J,XGSW,YGSW,ZGSW
  1362 +C----OUTPUT: XGSW,YGSW,ZGSW XSM,YSM,ZSM
  1363 +C
  1364 +C ATTENTION: SUBROUTINE RECALC_08 MUST BE INVOKED BEFORE SMGSW_08 IN THREE CASES:
  1365 +C /A/ BEFORE THE FIRST TRANSFORMATION OF COORDINATES
  1366 +C /B/ IF THE VALUES OF IYEAR,IDAY,IHOUR,MIN,ISEC HAVE BEEN CHANGED
  1367 +C /C/ IF THE VALUES OF COMPONENTS OF THE SOLAR WIND FLOW VELOCITY HAVE CHANGED
  1368 +C
  1369 +C IMPORTANT NOTE:
  1370 +C
  1371 +C A NON-STANDARD DEFINITION IS IMPLIED HERE FOR THE SOLAR MAGNETIC (SM) COORDINATE
  1372 +C SYSTEM: IT IS ASSUMED THAT THE XSM AXIS LIES IN THE PLANE DEFINED BY THE
  1373 +C GEODIPOLE AXIS AND THE OBSERVED VECTOR OF THE SOLAR WIND FLOW (RATHER THAN
  1374 +C THE EARTH-SUN LINE). IN ORDER TO CONVERT MAG COORDINATES TO AND FROM THE
  1375 +C STANDARD SM COORDINATES, INVOKE RECALC_08 WITH VGSEX=-400.0, VGSEY=0.0, VGSEZ=0.0
  1376 +C
  1377 +C LAST MODIFICATION: FEB 07, 2008
  1378 +C
  1379 +C AUTHOR: N. A. TSYGANENKO
  1380 +C
  1381 + COMMON /GEOPACK1/ A(10),SPS,CPS,B(22)
  1382 +
  1383 + IF (J.GT.0) THEN
  1384 + XGSW=XSM*CPS+ZSM*SPS
  1385 + YGSW=YSM
  1386 + ZGSW=ZSM*CPS-XSM*SPS
  1387 + ELSE
  1388 + XSM=XGSW*CPS-ZGSW*SPS
  1389 + YSM=YGSW
  1390 + ZSM=XGSW*SPS+ZGSW*CPS
  1391 + ENDIF
  1392 +
  1393 + RETURN
  1394 + END
  1395 +C
  1396 +C==========================================================================================
  1397 +C
  1398 + SUBROUTINE GEOGSW_08 (XGEO,YGEO,ZGEO,XGSW,YGSW,ZGSW,J)
  1399 +C
  1400 +C CONVERTS GEOGRAPHIC (GEO) TO GEOCENTRIC SOLAR-WIND (GSW) COORDINATES OR VICE VERSA.
  1401 +C
  1402 +C J>0 J<0
  1403 +C----- INPUT: J,XGEO,YGEO,ZGEO J,XGSW,YGSW,ZGSW
  1404 +C---- OUTPUT: XGSW,YGSW,ZGSW XGEO,YGEO,ZGEO
  1405 +C
  1406 +C ATTENTION: SUBROUTINE RECALC_08 MUST BE INVOKED BEFORE GEOGSW_08 IN THREE CASES:
  1407 +C /A/ BEFORE THE FIRST TRANSFORMATION OF COORDINATES, OR
  1408 +C /B/ IF THE VALUES OF IYEAR,IDAY,IHOUR,MIN,ISEC HAVE CHANGED, AND/OR
  1409 +C /C/ IF THE VALUES OF COMPONENTS OF THE SOLAR WIND FLOW VELOCITY HAVE CHANGED
  1410 +C
  1411 +C NOTE: THIS SUBROUTINE CONVERTS GEO VECTORS TO AND FROM THE SOLAR-WIND GSW COORDINATE
  1412 +C SYSTEM, TAKING INTO ACCOUNT POSSIBLE DEFLECTIONS OF THE SOLAR WIND DIRECTION FROM
  1413 +C STRICTLY RADIAL. BEFORE CONVERTING TO/FROM STANDARD GSM COORDINATES, INVOKE RECALC_08
  1414 +C WITH VGSEX=-400.0 and VGSEY=0.0, VGSEZ=0.0
  1415 +C
  1416 +C LAST MODIFICATION: FEB 07, 2008
  1417 +C
  1418 +C AUTHOR: N. A. TSYGANENKO
  1419 +C
  1420 + COMMON /GEOPACK1/ AA(16),A11,A21,A31,A12,A22,A32,A13,A23,A33,B(9)
  1421 +C
  1422 + IF (J.GT.0) THEN
  1423 + XGSW=A11*XGEO+A12*YGEO+A13*ZGEO
  1424 + YGSW=A21*XGEO+A22*YGEO+A23*ZGEO
  1425 + ZGSW=A31*XGEO+A32*YGEO+A33*ZGEO
  1426 + ELSE
  1427 + XGEO=A11*XGSW+A21*YGSW+A31*ZGSW
  1428 + YGEO=A12*XGSW+A22*YGSW+A32*ZGSW
  1429 + ZGEO=A13*XGSW+A23*YGSW+A33*ZGSW
  1430 + ENDIF
  1431 +
  1432 + RETURN
  1433 + END
  1434 +C
  1435 +C=====================================================================================
  1436 +C
  1437 + SUBROUTINE GEODGEO_08 (H,XMU,R,THETA,J)
  1438 +C
  1439 +C THIS SUBROUTINE (1) CONVERTS VERTICAL LOCAL HEIGHT (ALTITUDE) H AND GEODETIC
  1440 +C LATITUDE XMU INTO GEOCENTRIC COORDINATES R AND THETA (GEOCENTRIC RADIAL
  1441 +C DISTANCE AND COLATITUDE, RESPECTIVELY; ALSO KNOWN AS ECEF COORDINATES),
  1442 +C AS WELL AS (2) PERFORMS THE INVERSE TRANSFORMATION FROM {R,THETA} TO {H,XMU}.
  1443 +C
  1444 +C THE SUBROUTINE USES WORLD GEODETIC SYSTEM WGS84 PARAMETERS FOR THE EARTH'S
  1445 +C ELLIPSOID. THE ANGULAR QUANTITIES (GEO COLATITUDE THETA AND GEODETIC LATITUDE
  1446 +C XMU) ARE IN RADIANS, AND THE DISTANCES (GEOCENTRIC RADIUS R AND ALTITUDE H
  1447 +C ABOVE THE EARTH'S ELLIPSOID) ARE IN KILOMETERS.
  1448 +C
  1449 +C IF J>0, THE TRANSFORMATION IS MADE FROM GEODETIC TO GEOCENTRIC COORDINATES
  1450 +C USING SIMPLE DIRECT EQUATIONS.
  1451 +C IF J<0, THE INVERSE TRANSFORMATION FROM GEOCENTRIC TO GEODETIC COORDINATES
  1452 +C IS MADE BY MEANS OF A FAST ITERATIVE ALGORITHM.
  1453 +C
  1454 +c-------------------------------------------------------------------------------
  1455 +C J>0 | J<0
  1456 +c-------------------------------------------|-----------------------------------
  1457 +C--INPUT: J H XMU | J R THETA
  1458 +c flag altitude (km) geodetic | flag geocentric spherical
  1459 +c latitude | distance (km) colatitude
  1460 +c (radians) | (radians)
  1461 +c-------------------------------------------|-----------------------------------
  1462 +c |
  1463 +C----OUTPUT: R THETA | H XMU
  1464 +C geocentric spherical | altitude (km) geodetic
  1465 +C distance (km) colatitude | latitude
  1466 +C (radians) | (radians)
  1467 +C-------------------------------------------------------------------------------
  1468 +C
  1469 +C AUTHOR: N. A. TSYGANENKO
  1470 +c DATE: DEC 5, 2007
  1471 +C
  1472 + DATA R_EQ, BETA /6378.137, 6.73949674228E-3/
  1473 +c
  1474 +c R_EQ is the semi-major axis of the Earth's ellipsoid, and BETA is its
  1475 +c second eccentricity squared
  1476 +c
  1477 + DATA TOL /1.E-6/
  1478 +c
  1479 +c Direct transformation (GEOD=>GEO):
  1480 +c
  1481 + IF (J.GT.0) THEN
  1482 + COSXMU=COS(XMU)
  1483 + SINXMU=SIN(XMU)
  1484 + DEN=SQRT(COSXMU**2+(SINXMU/(1.0+BETA))**2)
  1485 + COSLAM=COSXMU/DEN
  1486 + SINLAM=SINXMU/(DEN*(1.0+BETA))
  1487 + RS=R_EQ/SQRT(1.0+BETA*SINLAM**2)
  1488 + X=RS*COSLAM+H*COSXMU
  1489 + Z=RS*SINLAM+H*SINXMU
  1490 + R=SQRT(X**2+Z**2)
  1491 + THETA=ACOS(Z/R)
  1492 + ENDIF
  1493 +
  1494 +c
  1495 +c Inverse transformation (GEO=>GEOD):
  1496 +c
  1497 + IF (J.LT.0) THEN
  1498 + N=0
  1499 + PHI=1.570796327-THETA
  1500 + PHI1=PHI
  1501 + 1 SP=SIN(PHI1)
  1502 + ARG=SP*(1.0D0+BETA)/SQRT(1.0+BETA*(2.0+BETA)*SP**2)
  1503 + XMUS=ASIN(ARG)
  1504 + RS=R_EQ/SQRT(1.0+BETA*SIN(PHI1)**2)
  1505 + COSFIMS=COS(PHI1-XMUS)
  1506 + H=SQRT((RS*COSFIMS)**2+R**2-RS**2)-RS*COSFIMS
  1507 + Z=RS*SIN(PHI1)+H*SIN(XMUS)
  1508 + X=RS*COS(PHI1)+H*COS(XMUS)
  1509 + RR=SQRT(X**2+Z**2)
  1510 + DPHI=ASIN(Z/RR)-PHI
  1511 + PHI1=PHI1-DPHI
  1512 + N=N+1
  1513 + IF (ABS(DPHI).GT.TOL.AND.N.LT.100) GOTO 1
  1514 + XMU=XMUS
  1515 + ENDIF
  1516 +
  1517 + RETURN
  1518 + END
  1519 +C
  1520 +C=====================================================================================
  1521 +C
  1522 + SUBROUTINE RHAND_08 (X,Y,Z,R1,R2,R3,IOPT,PARMOD,EXNAME,INNAME)
  1523 +C
  1524 +C CALCULATES THE COMPONENTS OF THE RIGHT HAND SIDE VECTOR IN THE GEOMAGNETIC FIELD
  1525 +C LINE EQUATION (a subsidiary subroutine for the subroutine STEP_08)
  1526 +C
  1527 +C LAST MODIFICATION: FEB 07, 2008
  1528 +C
  1529 +C AUTHOR: N. A. TSYGANENKO
  1530 +C
  1531 + DIMENSION PARMOD(10)
  1532 +C
  1533 +C EXNAME AND INNAME ARE NAMES OF SUBROUTINES FOR THE EXTERNAL AND INTERNAL
  1534 +C PARTS OF THE TOTAL FIELD, E.G., T96_01 AND IGRF_GSW_08
  1535 +C
  1536 + COMMON /GEOPACK1/ A(12),DS3,BB(2),PSI,CC(18)
  1537 +
  1538 + CALL EXNAME (IOPT,PARMOD,PSI,X,Y,Z,BXGSW,BYGSW,BZGSW)
  1539 + CALL INNAME (X,Y,Z,HXGSW,HYGSW,HZGSW)
  1540 +
  1541 + BX=BXGSW+HXGSW
  1542 + BY=BYGSW+HYGSW
  1543 + BZ=BZGSW+HZGSW
  1544 + B=DS3/SQRT(BX**2+BY**2+BZ**2)
  1545 + R1=BX*B
  1546 + R2=BY*B
  1547 + R3=BZ*B
  1548 + RETURN
  1549 + END
  1550 +C
  1551 +C===================================================================================
  1552 +C
  1553 + SUBROUTINE STEP_08(X,Y,Z,DS,DSMAX,ERRIN,IOPT,PARMOD,EXNAME,INNAME)
  1554 +C
  1555 +C RE-CALCULATES THE INPUT VALUES {X,Y,Z} (IN GSW COORDINATES) FOR ANY POINT ON A FIELD LINE,
  1556 +C BY MAKING A STEP ALONG THAT LINE USING RUNGE-KUTTA-MERSON ALGORITHM (G.N. Lance, Numerical
  1557 +C methods for high-speed computers, Iliffe & Sons, London 1960.)
  1558 +C DS IS A PRESCRIBED VALUE OF THE CURRENT STEP SIZE, DSMAX IS ITS UPPER LIMIT.
  1559 +C ERRIN IS A PERMISSIBLE ERROR (ITS OPTIMAL VALUE SPECIFIED IN THE S/R TRACE_08)
  1560 +C IF THE ACTUAL ERROR (ERRCUR) AT THE CURRENT STEP IS LARGER THAN ERRIN, THE STEP IS REJECTED,
  1561 +C AND THE CALCULATION IS REPEATED ANEW WITH HALVED STEPSIZE DS.
  1562 +C IF ERRCUR IS SMALLER THAN ERRIN, THE STEP IS ACCEPTED, AND THE CURRENT VALUE OF DS IS RETAINED
  1563 +C FOR THE NEXT STEP.
  1564 +C IF ERRCUR IS SMALLER THAN 0.04*ERRIN, THE STEP IS ACCEPTED, AND THE VALUE OF DS FOR THE NEXT STEP
  1565 +C IS INCREASED BY THE FACTOR 1.5, BUT NOT LARGER THAN DSMAX.
  1566 +C IOPT IS A FLAG, RESERVED FOR SPECIFYNG A VERSION OF THE EXTERNAL FIELD MODEL EXNAME.
  1567 +C ARRAY PARMOD(10) CONTAINS INPUT PARAMETERS FOR THE MODEL EXNAME.
  1568 +C EXNAME IS THE NAME OF THE SUBROUTINE FOR THE EXTERNAL FIELD MODEL.
  1569 +C INNAME IS THE NAME OF THE SUBROUTINE FOR THE INTERNAL FIELD MODEL (EITHER DIP_08 OR IGRF_GSW_08)
  1570 +C
  1571 +C ALL THE ABOVE PARAMETERS ARE INPUT ONES; OUTPUT IS THE RECALCULATED VALUES OF X,Y,Z
  1572 +C
  1573 +C LAST MODIFICATION: APR 21, 2008 (SEE ERRATA AS OF THIS DATE)
  1574 +C
  1575 +C AUTHOR: N. A. TSYGANENKO
  1576 +C
  1577 + DIMENSION PARMOD(10)
  1578 + COMMON /GEOPACK1/ A(12),DS3,B(21)
  1579 + EXTERNAL EXNAME,INNAME
  1580 +
  1581 + 1 DS3=-DS/3.
  1582 + CALL RHAND_08 (X,Y,Z,R11,R12,R13,IOPT,PARMOD,EXNAME,INNAME)
  1583 + CALL RHAND_08 (X+R11,Y+R12,Z+R13,R21,R22,R23,IOPT,PARMOD,EXNAME,
  1584 + * INNAME)
  1585 + CALL RHAND_08 (X+.5*(R11+R21),Y+.5*(R12+R22),Z+.5*
  1586 + *(R13+R23),R31,R32,R33,IOPT,PARMOD,EXNAME,INNAME)
  1587 + CALL RHAND_08 (X+.375*(R11+3.*R31),Y+.375*(R12+3.*R32
  1588 + *),Z+.375*(R13+3.*R33),R41,R42,R43,IOPT,PARMOD,EXNAME,INNAME)
  1589 + CALL RHAND_08 (X+1.5*(R11-3.*R31+4.*R41),Y+1.5*(R12-
  1590 + *3.*R32+4.*R42),Z+1.5*(R13-3.*R33+4.*R43),
  1591 + *R51,R52,R53,IOPT,PARMOD,EXNAME,INNAME)
  1592 + ERRCUR=ABS(R11-4.5*R31+4.*R41-.5*R51)+ABS(R12-4.5*R32+4.*R42-.5*
  1593 + *R52)+ABS(R13-4.5*R33+4.*R43-.5*R53)
  1594 +C
  1595 +C READY FOR MAKING THE STEP, BUT CHECK THE ACCURACY; IF INSUFFICIENT,
  1596 +C REPEAT THE STEP WITH HALVED STEPSIZE:
  1597 +C
  1598 + IF (ERRCUR.GT.ERRIN) THEN
  1599 + DS=DS*.5
  1600 + GOTO 1
  1601 + ENDIF
  1602 +C
  1603 +C ACCURACY IS ACCEPTABLE, BUT CHECK IF THE STEPSIZE IS NOT TOO LARGE;
  1604 +C OTHERWISE REPEAT THE STEP WITH DS=DSMAX
  1605 +C
  1606 + IF (ABS(DS).GT.DSMAX) THEN
  1607 + DS=SIGN(DSMAX,DS)
  1608 + GOTO 1
  1609 + ENDIF
  1610 +C
  1611 +C MAKING THE STEP:
  1612 +C
  1613 + 2 X=X+.5*(R11+4.*R41+R51)
  1614 + Y=Y+.5*(R12+4.*R42+R52)
  1615 + Z=Z+.5*(R13+4.*R43+R53)
  1616 +C
  1617 +C IF THE ACTUAL ERROR IS TOO SMALL (LESS THAN 4% OF ERRIN) AND DS SMALLER
  1618 +C THAN DSMAX/1.5, THEN WE INCREASE THE STEPSIZE FOR THE NEXT STEP BY 50%
  1619 +C
  1620 + IF(ERRCUR.LT.ERRIN*.04.AND.DS.LT.DSMAX/1.5) DS=DS*1.5
  1621 + RETURN
  1622 + END
  1623 +C
  1624 +C==============================================================================
  1625 +C
  1626 + SUBROUTINE TRACE_08 (XI,YI,ZI,DIR,DSMAX,ERR,RLIM,R0,IOPT,PARMOD,
  1627 + * EXNAME,INNAME,XF,YF,ZF,XX,YY,ZZ,L,LMAX)
  1628 +C
  1629 +C TRACES A FIELD LINE FROM AN ARBITRARY POINT OF SPACE TO THE EARTH'S
  1630 +C SURFACE OR TO A MODEL LIMITING BOUNDARY.
  1631 +C
  1632 +C THIS SUBROUTINE ALLOWS TWO OPTIONS:
  1633 +C
  1634 +C (1) IF INNAME=IGRF_GSW_08, THEN THE IGRF MODEL WILL BE USED FOR CALCULATING
  1635 +C CONTRIBUTION FROM EARTH'S INTERNAL SOURCES. IN THIS CASE, SUBROUTINE
  1636 +C RECALC_08 MUST BE CALLED BEFORE USING TRACE_08, WITH PROPERLY SPECIFIED DATE,
  1637 +C UNIVERSAL TIME, AND SOLAR WIND VELOCITY COMPONENTS, TO CALCULATE IN ADVANCE
  1638 +C ALL QUANTITIES NEEDED FOR THE MAIN FIELD MODEL AND FOR TRANSFORMATIONS
  1639 +C BETWEEN INVOLVED COORDINATE SYSTEMS.
  1640 +C
  1641 +C (2) IF INNAME=DIP_08, THEN A PURE DIPOLE FIELD WILL BE USED INSTEAD OF THE IGRF MODEL.
  1642 +C IN THIS CASE, THE SUBROUTINE RECALC_08 MUST ALSO BE CALLED BEFORE TRACE_08.
  1643 +C HERE ONE CAN CHOOSE EITHER TO
  1644 +C (a) CALCULATE DIPOLE TILT ANGLE BASED ON DATE, TIME, AND SOLAR WIND DIRECTION,
  1645 +C OR (b) EXPLICITLY SPECIFY THAT ANGLE, WITHOUT ANY REFERENCE TO DATE/UT/SOLAR WIND.
  1646 +C IN THE LAST CASE (b), THE SINE (SPS) AND COSINE (CPS) OF THE DIPOLE TILT
  1647 +C ANGLE MUST BE SPECIFIED IN ADVANCE (BUT AFTER HAVING CALLED RECALC_08) AND FORWARDED
  1648 +C IN THE COMMON BLOCK /GEOPACK1/ (IN ITS 11th AND 12th ELEMENTS, RESPECTIVELY).
  1649 +C IN THIS CASE THE ROLE OF THE SUBROUTINE RECALC_08 IS REDUCED TO ONLY CALCULATING
  1650 +C THE COMPONENTS OF THE EARTH'S DIPOLE MOMENT.
  1651 +C
  1652 +C------------- INPUT PARAMETERS:
  1653 +C
  1654 +C XI,YI,ZI - GSW COORDS OF THE FIELD LINE STARTING POINT (IN EARTH RADII, 1 RE = 6371.2 km),
  1655 +C
  1656 +C DIR - SIGN OF THE TRACING DIRECTION: IF DIR=1.0 THEN THE TRACING IS MADE ANTIPARALLEL
  1657 +C TO THE TOTAL FIELD VECTOR (E.G., FROM NORTHERN TO SOUTHERN CONJUGATE POINT);
  1658 +C IF DIR=-1.0 THEN THE TRACING PROCEEDS IN THE OPPOSITE DIRECTION, THAT IS, PARALLEL TO
  1659 +C THE TOTAL FIELD VECTOR.
  1660 +C
  1661 +C DSMAX - UPPER LIMIT ON THE STEPSIZE (SETS A DESIRED MAXIMAL SPACING BETWEEN
  1662 +C THE FIELD LINE POINTS)
  1663 +C
  1664 +C ERR - PERMISSIBLE STEP ERROR. A REASONABLE ESTIMATE PROVIDING A SUFFICIENT ACCURACY FOR MOST
  1665 +C APPLICATIONS IS ERR=0.0001. SMALLER/LARGER VALUES WILL RESULT IN LARGER/SMALLER NUMBER
  1666 +C OF STEPS AND, HENCE, OF OUTPUT FIELD LINE POINTS. NOTE THAT USING MUCH SMALLER VALUES
  1667 +C OF ERR MAY REQUIRE USING A DOUBLE PRECISION VERSION OF THE ENTIRE PACKAGE.
  1668 +C
  1669 +C R0 - RADIUS OF A SPHERE (IN RE), DEFINING THE INNER BOUNDARY OF THE TRACING REGION
  1670 +C (USUALLY, EARTH'S SURFACE OR THE IONOSPHERE, WHERE R0~1.0)
  1671 +C IF THE FIELD LINE REACHES THAT SPHERE FROM OUTSIDE, ITS INBOUND TRACING IS
  1672 +C TERMINATED AND THE CROSSING POINT COORDINATES XF,YF,ZF ARE CALCULATED.
  1673 +C
  1674 +C RLIM - RADIUS OF A SPHERE (IN RE), DEFINING THE OUTER BOUNDARY OF THE TRACING REGION;
  1675 +C IF THE FIELD LINE REACHES THAT BOUNDARY FROM INSIDE, ITS OUTBOUND TRACING IS
  1676 +C TERMINATED AND THE CROSSING POINT COORDINATES XF,YF,ZF ARE CALCULATED.
  1677 +C
  1678 +C IOPT - A MODEL INDEX; CAN BE USED FOR SPECIFYING A VERSION OF THE EXTERNAL FIELD
  1679 +C MODEL (E.G., A NUMBER OF THE KP-INDEX INTERVAL). ALTERNATIVELY, ONE CAN USE THE ARRAY
  1680 +C PARMOD FOR THAT PURPOSE (SEE BELOW); IN THAT CASE IOPT IS JUST A DUMMY PARAMETER.
  1681 +C
  1682 +C PARMOD - A 10-ELEMENT ARRAY CONTAINING INPUT PARAMETERS NEEDED FOR A UNIQUE
  1683 +C SPECIFICATION OF THE EXTERNAL FIELD MODEL. THE CONCRETE MEANING OF THE COMPONENTS
  1684 +C OF PARMOD DEPENDS ON A SPECIFIC VERSION OF THAT MODEL.
  1685 +C
  1686 +C EXNAME - NAME OF A SUBROUTINE PROVIDING COMPONENTS OF THE EXTERNAL MAGNETIC FIELD
  1687 +C (E.G., T89, OR T96_01, ETC.).
  1688 +C INNAME - NAME OF A SUBROUTINE PROVIDING COMPONENTS OF THE INTERNAL MAGNETIC FIELD
  1689 +C (EITHER DIP_08 OR IGRF_GSW_08).
  1690 +C
  1691 +C LMAX - MAXIMAL LENGTH OF THE ARRAYS XX,YY,ZZ, IN WHICH COORDINATES OF THE FIELD
  1692 +C LINE POINTS ARE STORED. LMAX SHOULD BE SET EQUAL TO THE ACTUAL LENGTH OF
  1693 +C THE ARRAYS, DEFINED IN THE MAIN PROGRAM AS ACTUAL ARGUMENTS OF THIS SUBROUTINE.
  1694 +C
  1695 +C-------------- OUTPUT PARAMETERS:
  1696 +C
  1697 +C XF,YF,ZF - GSW COORDINATES OF THE ENDPOINT OF THE TRACED FIELD LINE.
  1698 +C XX,YY,ZZ - ARRAYS OF LENGTH LMAX, CONTAINING COORDINATES OF THE FIELD LINE POINTS.
  1699 +C L - ACTUAL NUMBER OF FIELD LINE POINTS, GENERATED BY THIS SUBROUTINE.
  1700 +C
  1701 +C ----------------------------------------------------------
  1702 +C
  1703 +C LAST MODIFICATION: JAN 30, 2008.
  1704 +C
  1705 +C AUTHOR: N. A. TSYGANENKO
  1706 +C
  1707 + DIMENSION XX(LMAX),YY(LMAX),ZZ(LMAX), PARMOD(10)
  1708 + COMMON /GEOPACK1/ AA(12),DD,BB(21)
  1709 + EXTERNAL EXNAME,INNAME
  1710 +C
  1711 + L=0
  1712 + NREV=0
  1713 + DD=DIR
  1714 +C
  1715 +C INITIALIZE THE STEP SIZE AND STARTING PONT:
  1716 +C
  1717 + DS=0.5*DIR
  1718 + X=XI
  1719 + Y=YI
  1720 + Z=ZI
  1721 +c
  1722 +c here we call RHAND_08 just to find out the sign of the radial component of the field
  1723 +c vector, and to determine the initial direction of the tracing (i.e., either away
  1724 +c or towards Earth):
  1725 +c
  1726 + CALL RHAND_08 (X,Y,Z,R1,R2,R3,IOPT,PARMOD,EXNAME,INNAME)
  1727 + AD=0.01
  1728 + IF (X*R1+Y*R2+Z*R3.LT.0.) AD=-0.01
  1729 +C
  1730 +c |AD|=0.01 and its sign follows the rule:
  1731 +c (1) if DIR=1 (tracing antiparallel to B vector) then the sign of AD is the same as of Br
  1732 +c (2) if DIR=-1 (tracing parallel to B vector) then the sign of AD is opposite to that of Br
  1733 +c AD is defined in order to initialize the value of RR (radial distance at previous step):
  1734 +
  1735 + RR=SQRT(X**2+Y**2+Z**2)+AD
  1736 +c
  1737 + 1 L=L+1
  1738 + IF(L.GT.LMAX) GOTO 7
  1739 + XX(L)=X
  1740 + YY(L)=Y
  1741 + ZZ(L)=Z
  1742 + RYZ=Y**2+Z**2
  1743 + R2=X**2+RYZ
  1744 + R=SQRT(R2)
  1745 +C
  1746 +c check if the line hit the outer tracing boundary; if yes, then terminate
  1747 +c the tracing (label 8). The outer boundary is assumed reached, when the line
  1748 +c crosses any of the 3 surfaces: (1) a sphere R=RLIM, (2) a cylinder of radius 40Re,
  1749 +c coaxial with the XGSW axis, (3) the plane X=20Re:
  1750 +
  1751 + IF (R.GT.RLIM.OR.RYZ.GT.1600..OR.X.GT.20.) GOTO 8
  1752 +c
  1753 +c check whether or not the inner tracing boundary was crossed from outside,
  1754 +c if yes, then calculate the footpoint position by interpolation (go to label 6):
  1755 +c
  1756 + IF (R.LT.R0.AND.RR.GT.R) GOTO 6
  1757 +
  1758 +c check if we are moving outward, or R is still larger than 3Re; if yes, proceed further:
  1759 +c
  1760 + IF (R.GE.RR.OR.R.GE.3.) GOTO 4
  1761 +c
  1762 +c now we entered inside the sphere R=3: to avoid too large steps (and hence
  1763 +c inaccurate interpolated position of the footpoint), enforce the progressively
  1764 +c smaller stepsize values as we approach the inner boundary R=R0:
  1765 +c
  1766 + FC=0.2
  1767 + IF(R-R0.LT.0.05) FC=0.05
  1768 + AL=FC*(R-R0+0.2)
  1769 + DS=DIR*AL
  1770 +c
  1771 + 4 XR=X
  1772 + YR=Y
  1773 + ZR=Z
  1774 +c
  1775 + DRP=R-RR
  1776 + RR=R
  1777 +c
  1778 + CALL STEP_08 (X,Y,Z,DS,DSMAX,ERR,IOPT,PARMOD,EXNAME,INNAME)
  1779 +c
  1780 +C check the total number NREV of changes in the tracing radial direction; (NREV.GT.2) means
  1781 +c that the line started making multiple loops, in which case we stop the process:
  1782 +C
  1783 + R=SQRT(X**2+Y**2+Z**2)
  1784 + DR=R-RR
  1785 + IF (DRP*DR.LT.0.) NREV=NREV+1
  1786 + IF (NREV.GT.4) GOTO 8
  1787 +C
  1788 + GOTO 1
  1789 +c
  1790 +c find the footpoint position by interpolating between the current and previous
  1791 +c field line points:
  1792 +c
  1793 + 6 R1=(R0-R)/(RR-R)
  1794 + X=X-(X-XR)*R1
  1795 + Y=Y-(Y-YR)*R1
  1796 + Z=Z-(Z-ZR)*R1
  1797 + GOTO 8
  1798 + 7 WRITE (*,10)
  1799 + L=LMAX
  1800 + 8 XF=X
  1801 + YF=Y
  1802 + ZF=Z
  1803 +C
  1804 +C replace the coordinates of the last (L-th) point in the XX,YY,ZZ arrays
  1805 +C so that they correspond to the estimated footpoint position {XF,YF,ZF},
  1806 +c satisfying: sqrt(XF**2+YF**2+ZF**2}=R0
  1807 +C
  1808 + XX(L)=XF
  1809 + YY(L)=YF
  1810 + ZZ(L)=ZF
  1811 +C
  1812 + RETURN
  1813 + 10 FORMAT(//,1X,'**** COMPUTATIONS IN THE SUBROUTINE TRACE_08 ARE',
  1814 + *' TERMINATED: THE NUMBER OF POINTS EXCEEDED LMAX ****'//)
  1815 + END
  1816 +c
  1817 +C====================================================================================
  1818 +C
  1819 + SUBROUTINE SHUETAL_MGNP_08(XN_PD,VEL,BZIMF,XGSW,YGSW,ZGSW,
  1820 + * XMGNP,YMGNP,ZMGNP,DIST,ID)
  1821 +C
  1822 +C FOR ANY POINT OF SPACE WITH COORDINATES (XGSW,YGSW,ZGSW) AND SPECIFIED CONDITIONS
  1823 +C IN THE INCOMING SOLAR WIND, THIS SUBROUTINE:
  1824 +C
  1825 +C (1) DETERMINES IF THE POINT (XGSW,YGSW,ZGSW) LIES INSIDE OR OUTSIDE THE
  1826 +C MODEL MAGNETOPAUSE OF SHUE ET AL. (JGR-A, V.103, P. 17691, 1998).
  1827 +C
  1828 +C (2) CALCULATES THE GSW POSITION OF A POINT {XMGNP,YMGNP,ZMGNP}, LYING AT THE MODEL
  1829 +C MAGNETOPAUSE AND ASYMPTOTICALLY TENDING TO THE NEAREST BOUNDARY POINT WITH
  1830 +C RESPECT TO THE OBSERVATION POINT {XGSW,YGSW,ZGSW}, AS IT APPROACHES THE MAGNETO-
  1831 +C PAUSE.
  1832 +C
  1833 +C INPUT: XN_PD - EITHER SOLAR WIND PROTON NUMBER DENSITY (PER C.C.) (IF VEL>0)
  1834 +C OR THE SOLAR WIND RAM PRESSURE IN NANOPASCALS (IF VEL<0)
  1835 +C BZIMF - IMF BZ IN NANOTESLAS
  1836 +C
  1837 +C VEL - EITHER SOLAR WIND VELOCITY (KM/SEC)
  1838 +C OR ANY NEGATIVE NUMBER, WHICH INDICATES THAT XN_PD STANDS
  1839 +C FOR THE SOLAR WIND PRESSURE, RATHER THAN FOR THE DENSITY
  1840 +C
  1841 +C XGSW,YGSW,ZGSW - GSW POSITION OF THE OBSERVATION POINT IN EARTH RADII
  1842 +C
  1843 +C OUTPUT: XMGNP,YMGNP,ZMGNP - GSW POSITION OF THE BOUNDARY POINT
  1844 +C DIST - DISTANCE (IN RE) BETWEEN THE OBSERVATION POINT (XGSW,YGSW,ZGSW)
  1845 +C AND THE MODEL NAGNETOPAUSE
  1846 +C ID - POSITION FLAG: ID=+1 (-1) MEANS THAT THE OBSERVATION POINT
  1847 +C LIES INSIDE (OUTSIDE) OF THE MODEL MAGNETOPAUSE, RESPECTIVELY.
  1848 +C
  1849 +C OTHER SUBROUTINES USED: T96_MGNP_08
  1850 +C
  1851 +c AUTHOR: N.A. TSYGANENKO,
  1852 +C DATE: APRIL 4, 2003.
  1853 +C
  1854 + IF (VEL.LT.0.) THEN
  1855 + P=XN_PD
  1856 + ELSE
  1857 + P=1.94E-6*XN_PD*VEL**2 ! P IS THE SOLAR WIND DYNAMIC PRESSURE (IN nPa)
  1858 + ENDIF
  1859 +
  1860 +c
  1861 +c DEFINE THE ANGLE PHI, MEASURED DUSKWARD FROM THE NOON-MIDNIGHT MERIDIAN PLANE;
  1862 +C IF THE OBSERVATION POINT LIES ON THE X AXIS, THE ANGLE PHI CANNOT BE UNIQUELY
  1863 +C DEFINED, AND WE SET IT AT ZERO:
  1864 +c
  1865 + IF (YGSW.NE.0..OR.ZGSW.NE.0.) THEN
  1866 + PHI=ATAN2(YGSW,ZGSW)
  1867 + ELSE
  1868 + PHI=0.
  1869 + ENDIF
  1870 +C
  1871 +C FIRST, FIND OUT IF THE OBSERVATION POINT LIES INSIDE THE SHUE ET AL BDRY
  1872 +C AND SET THE VALUE OF THE ID FLAG:
  1873 +C
  1874 + ID=-1
  1875 + R0=(10.22+1.29*TANH(0.184*(BZIMF+8.14)))*P**(-.15151515)
  1876 + ALPHA=(0.58-0.007*BZIMF)*(1.+0.024*ALOG(P))
  1877 + R=SQRT(XGSW**2+YGSW**2+ZGSW**2)
  1878 + RM=R0*(2./(1.+XGSW/R))**ALPHA
  1879 + IF (R.LE.RM) ID=+1
  1880 +C
  1881 +C NOW, FIND THE CORRESPONDING T96 MAGNETOPAUSE POSITION, TO BE USED AS
  1882 +C A STARTING APPROXIMATION IN THE SEARCH OF A CORRESPONDING SHUE ET AL.
  1883 +C BOUNDARY POINT:
  1884 +C
  1885 + CALL T96_MGNP_08(P,-1.,XGSW,YGSW,ZGSW,XMT96,YMT96,ZMT96,DIST,ID96)
  1886 +C
  1887 + RHO2=YMT96**2+ZMT96**2
  1888 + R=SQRT(RHO2+XMT96**2)
  1889 + ST=SQRT(RHO2)/R
  1890 + CT=XMT96/R
  1891 +C
  1892 +C NOW, USE NEWTON'S ITERATIVE METHOD TO FIND THE NEAREST POINT AT THE
  1893 +C SHUE ET AL.'S BOUNDARY:
  1894 +C
  1895 + NIT=0
  1896 +
  1897 + 1 T=ATAN2(ST,CT)
  1898 + RM=R0*(2./(1.+CT))**ALPHA
  1899 +
  1900 + F=R-RM
  1901 + GRADF_R=1.
  1902 + GRADF_T=-ALPHA/R*RM*ST/(1.+CT)
  1903 + GRADF=SQRT(GRADF_R**2+GRADF_T**2)
  1904 +
  1905 + DR=-F/GRADF**2
  1906 + DT= DR/R*GRADF_T
  1907 +
  1908 + R=R+DR
  1909 + T=T+DT
  1910 + ST=SIN(T)
  1911 + CT=COS(T)
  1912 +
  1913 + DS=SQRT(DR**2+(R*DT)**2)
  1914 +
  1915 + NIT=NIT+1
  1916 +
  1917 + IF (NIT.GT.1000) THEN
  1918 + PRINT *,
  1919 + *' BOUNDARY POINT COULD NOT BE FOUND; ITERATIONS DO NOT CONVERGE'
  1920 + ENDIF
  1921 +
  1922 + IF (DS.GT.1.E-4) GOTO 1
  1923 +
  1924 + XMGNP=R*COS(T)
  1925 + RHO= R*SIN(T)
  1926 +
  1927 + YMGNP=RHO*SIN(PHI)
  1928 + ZMGNP=RHO*COS(PHI)
  1929 +
  1930 + DIST=SQRT((XGSW-XMGNP)**2+(YGSW-YMGNP)**2+(ZGSW-ZMGNP)**2)
  1931 +
  1932 + RETURN
  1933 + END
  1934 +C
  1935 +C=======================================================================================
  1936 +C
  1937 + SUBROUTINE T96_MGNP_08(XN_PD,VEL,XGSW,YGSW,ZGSW,XMGNP,YMGNP,ZMGNP,
  1938 + * DIST,ID)
  1939 +C
  1940 +C FOR ANY POINT OF SPACE WITH GIVEN COORDINATES (XGSW,YGSW,ZGSW), THIS SUBROUTINE DEFINES
  1941 +C THE POSITION OF A POINT (XMGNP,YMGNP,ZMGNP) AT THE T96 MODEL MAGNETOPAUSE WITH THE
  1942 +C SAME VALUE OF THE ELLIPSOIDAL TAU-COORDINATE, AND THE DISTANCE BETWEEN THEM. THIS IS
  1943 +C NOT THE SHORTEST DISTANCE D_MIN TO THE BOUNDARY, BUT DIST ASYMPTOTICALLY TENDS TO D_MIN,
  1944 +C AS THE OBSERVATION POINT GETS CLOSER TO THE MAGNETOPAUSE.
  1945 +C
  1946 +C INPUT: XN_PD - EITHER SOLAR WIND PROTON NUMBER DENSITY (PER C.C.) (IF VEL>0)
  1947 +C OR THE SOLAR WIND RAM PRESSURE IN NANOPASCALS (IF VEL<0)
  1948 +C VEL - EITHER SOLAR WIND VELOCITY (KM/SEC)
  1949 +C OR ANY NEGATIVE NUMBER, WHICH INDICATES THAT XN_PD STANDS
  1950 +C FOR THE SOLAR WIND PRESSURE, RATHER THAN FOR THE DENSITY
  1951 +C
  1952 +C XGSW,YGSW,ZGSW - COORDINATES OF THE OBSERVATION POINT IN EARTH RADII
  1953 +C
  1954 +C OUTPUT: XMGNP,YMGNP,ZMGNP - GSW POSITION OF THE BOUNDARY POINT, HAVING THE SAME
  1955 +C VALUE OF TAU-COORDINATE AS THE OBSERVATION POINT (XGSW,YGSW,ZGSW)
  1956 +C DIST - THE DISTANCE BETWEEN THE TWO POINTS, IN RE,
  1957 +C ID - POSITION FLAG; ID=+1 (-1) MEANS THAT THE POINT (XGSW,YGSW,ZGSW)
  1958 +C LIES INSIDE (OUTSIDE) THE MODEL MAGNETOPAUSE, RESPECTIVELY.
  1959 +C
  1960 +C THE PRESSURE-DEPENDENT MAGNETOPAUSE IS THAT USED IN THE T96_01 MODEL
  1961 +C (TSYGANENKO, JGR, V.100, P.5599, 1995; ESA SP-389, P.181, OCT. 1996)
  1962 +C
  1963 +c AUTHOR: N.A. TSYGANENKO
  1964 +C DATE: AUG.1, 1995, REVISED APRIL 3, 2003.
  1965 +C
  1966 +C
  1967 +C DEFINE SOLAR WIND DYNAMIC PRESSURE (NANOPASCALS, ASSUMING 4% OF ALPHA-PARTICLES),
  1968 +C IF NOT EXPLICITLY SPECIFIED IN THE INPUT:
  1969 +
  1970 + IF (VEL.LT.0.) THEN
  1971 + PD=XN_PD
  1972 + ELSE
  1973 + PD=1.94E-6*XN_PD*VEL**2
  1974 +C
  1975 + ENDIF
  1976 +C
  1977 +C RATIO OF PD TO THE AVERAGE PRESSURE, ASSUMED EQUAL TO 2 nPa:
  1978 +
  1979 + RAT=PD/2.0
  1980 + RAT16=RAT**0.14
  1981 +
  1982 +C (THE POWER INDEX 0.14 IN THE SCALING FACTOR IS THE BEST-FIT VALUE OBTAINED FROM DATA
  1983 +C AND USED IN THE T96_01 VERSION)
  1984 +C
  1985 +C VALUES OF THE MAGNETOPAUSE PARAMETERS FOR PD = 2 nPa:
  1986 +C
  1987 + A0=70.
  1988 + S00=1.08
  1989 + X00=5.48
  1990 +C
  1991 +C VALUES OF THE MAGNETOPAUSE PARAMETERS, SCALED BY THE ACTUAL PRESSURE:
  1992 +C
  1993 + A=A0/RAT16
  1994 + S0=S00
  1995 + X0=X00/RAT16
  1996 + XM=X0-A
  1997 +C
  1998 +C (XM IS THE X-COORDINATE OF THE "SEAM" BETWEEN THE ELLIPSOID AND THE CYLINDER)
  1999 +C
  2000 +C (FOR DETAILS OF THE ELLIPSOIDAL COORDINATES, SEE THE PAPER:
  2001 +C N.A.TSYGANENKO, SOLUTION OF CHAPMAN-FERRARO PROBLEM FOR AN
  2002 +C ELLIPSOIDAL MAGNETOPAUSE, PLANET.SPACE SCI., V.37, P.1037, 1989).
  2003 +C
  2004 + IF (YGSW.NE.0..OR.ZGSW.NE.0.) THEN
  2005 + PHI=ATAN2(YGSW,ZGSW)
  2006 + ELSE
  2007 + PHI=0.
  2008 + ENDIF
  2009 +C
  2010 + RHO=SQRT(YGSW**2+ZGSW**2)
  2011 +C
  2012 + IF (XGSW.LT.XM) THEN
  2013 + XMGNP=XGSW
  2014 + RHOMGNP=A*SQRT(S0**2-1)
  2015 + YMGNP=RHOMGNP*SIN(PHI)
  2016 + ZMGNP=RHOMGNP*COS(PHI)
  2017 + DIST=SQRT((XGSW-XMGNP)**2+(YGSW-YMGNP)**2+(ZGSW-ZMGNP)**2)
  2018 + IF (RHOMGNP.GT.RHO) ID=+1
  2019 + IF (RHOMGNP.LE.RHO) ID=-1
  2020 + RETURN
  2021 + ENDIF
  2022 +C
  2023 + XKSI=(XGSW-X0)/A+1.
  2024 + XDZT=RHO/A
  2025 + SQ1=SQRT((1.+XKSI)**2+XDZT**2)
  2026 + SQ2=SQRT((1.-XKSI)**2+XDZT**2)
  2027 + SIGMA=0.5*(SQ1+SQ2)
  2028 + TAU=0.5*(SQ1-SQ2)
  2029 +C
  2030 +C NOW CALCULATE (X,Y,Z) FOR THE CLOSEST POINT AT THE MAGNETOPAUSE
  2031 +C
  2032 + XMGNP=X0-A*(1.-S0*TAU)
  2033 + ARG=(S0**2-1.)*(1.-TAU**2)
  2034 + IF (ARG.LT.0.) ARG=0.
  2035 + RHOMGNP=A*SQRT(ARG)
  2036 + YMGNP=RHOMGNP*SIN(PHI)
  2037 + ZMGNP=RHOMGNP*COS(PHI)
  2038 +C
  2039 +C NOW CALCULATE THE DISTANCE BETWEEN THE POINTS {XGSW,YGSW,ZGSW} AND {XMGNP,YMGNP,ZMGNP}:
  2040 +C (IN GENERAL, THIS IS NOT THE SHORTEST DISTANCE D_MIN, BUT DIST ASYMPTOTICALLY TENDS
  2041 +C TO D_MIN, AS WE ARE GETTING CLOSER TO THE MAGNETOPAUSE):
  2042 +C
  2043 + DIST=SQRT((XGSW-XMGNP)**2+(YGSW-YMGNP)**2+(ZGSW-ZMGNP)**2)
  2044 +C
  2045 + IF (SIGMA.GT.S0) ID=-1 ! ID=-1 MEANS THAT THE POINT LIES OUTSIDE
  2046 + IF (SIGMA.LE.S0) ID=+1 ! ID=+1 MEANS THAT THE POINT LIES INSIDE
  2047 +C THE MAGNETOSPHERE
  2048 + RETURN
  2049 + END
  2050 +C
  2051 +C===================================================================================
  2052 +C
  2053 +c
geopack/src/Modif.for 0 → 100644
  Diff comments   View file @55c4dfc
@@ -0,0 +1,10 @@ @@ -0,0 +1,10 @@
  1 + SUBROUTINE TRACE_08_MODIFIED (XI,YI,ZI,DIR,DSMAX,ERR,RLIM,R0,IOPT,
  2 + * PARMOD,XF,YF,ZF,XX,YY,ZZ,L,LMAX)
  3 +
  4 + DIMENSION XX(LMAX),YY(LMAX),ZZ(LMAX),PARMOD(10)
  5 + EXTERNAL T96_01,IGRF_GSW_08
  6 +
  7 + CALL TRACE_08 (XI,YI,ZI,DIR,DSMAX,ERR,RLIM,R0,IOPT,PARMOD,
  8 + * T96_01,IGRF_GSW_08,XF,YF,ZF,XX,YY,ZZ,L,LMAX)
  9 +
  10 + END
0 \ No newline at end of file 11 \ No newline at end of file
geopack/src/T96.for 0 → 100644
  Diff comments   View file @55c4dfc
@@ -0,0 +1,2571 @@ @@ -0,0 +1,2571 @@
  1 +c
  2 +C----------------------------------------------------------------------
  3 +c
  4 + SUBROUTINE T96_01 (IOPT,PARMOD,PS,X,Y,Z,BX,BY,BZ)
  5 +C
  6 +c RELEASE DATE OF THIS VERSION: JUNE 22, 1996.
  7 +C LAST UPDATE: MAY 01, 2006: IN THE S/R DIPOLE, SPS AND CPS WERE ADDED IN THE SAVE STATEMENT
  8 +
  9 +C----------------------------------------------------------------------
  10 +C
  11 +C WITH TWO CORRECTIONS, SUGGESTED BY T.SOTIRELIS' COMMENTS (APR.7, 1997)
  12 +C
  13 +C (1) A "STRAY " CLOSING PARENTHESIS WAS REMOVED IN THE S/R R2_BIRK
  14 +C (2) A 0/0 PROBLEM ON THE Z-AXIS WAS SIDESTEPPED (LINES 44-46 OF THE
  15 +c DOUBLE PRECISION FUNCTION XKSI)
  16 +c--------------------------------------------------------------------
  17 +C DATA-BASED MODEL CALIBRATED BY (1) SOLAR WIND PRESSURE PDYN (NANOPASCALS),
  18 +C (2) DST (NANOTESLA), (3) BYIMF, AND (4) BZIMF (NANOTESLA).
  19 +c THESE INPUT PARAMETERS SHOULD BE PLACED IN THE FIRST 4 ELEMENTS
  20 +c OF THE ARRAY PARMOD(10).
  21 +C
  22 +C THE REST OF THE INPUT VARIABLES ARE: THE GEODIPOLE TILT ANGLE PS (RADIANS),
  23 +C AND X,Y,Z - GSM POSITION (RE)
  24 +C
  25 +c IOPT IS JUST A DUMMY INPUT PARAMETER, NECESSARY TO MAKE THIS SUBROUTINE
  26 +C COMPATIBLE WITH THE NEW RELEASE (APRIL 1996) OF THE TRACING SOFTWARE
  27 +C PACKAGE (GEOPACK). IOPT VALUE DOES NOT AFFECT THE OUTPUT FIELD.
  28 +c
  29 +C
  30 +c OUTPUT: GSM COMPONENTS OF THE EXTERNAL MAGNETIC FIELD (BX,BY,BZ, nanotesla)
  31 +C COMPUTED AS A SUM OF CONTRIBUTIONS FROM PRINCIPAL FIELD SOURCES
  32 +C
  33 +c (C) Copr. 1995, 1996, Nikolai A. Tsyganenko, Hughes STX, Code 695, NASA GSFC
  34 +c Greenbelt, MD 20771, USA
  35 +c
  36 +C REFERENCES:
  37 +C
  38 +C (1) N.A. TSYGANENKO AND D.P. STERN, A NEW-GENERATION GLOBAL
  39 +C MAGNETOSPHERE FIELD MODEL , BASED ON SPACECRAFT MAGNETOMETER DATA,
  40 +C ISTP NEWSLETTER, V.6, NO.1, P.21, FEB.1996.
  41 +C
  42 +c (2) N.A.TSYGANENKO, MODELING THE EARTH'S MAGNETOSPHERIC
  43 +C MAGNETIC FIELD CONFINED WITHIN A REALISTIC MAGNETOPAUSE,
  44 +C J.GEOPHYS.RES., V.100, P. 5599, 1995.
  45 +C
  46 +C (3) N.A. TSYGANENKO AND M.PEREDO, ANALYTICAL MODELS OF THE
  47 +C MAGNETIC FIELD OF DISK-SHAPED CURRENT SHEETS, J.GEOPHYS.RES.,
  48 +C V.99, P. 199, 1994.
  49 +C
  50 +c----------------------------------------------------------------------
  51 +
  52 + IMPLICIT REAL*8 (A-H,O-Z)
  53 + REAL PDYN,DST,BYIMF,BZIMF,PS,X,Y,Z,BX,BY,BZ,QX,QY,QZ,PARMOD(10),
  54 + * A(9)
  55 +c
  56 + DATA PDYN0,EPS10 /2.,3630.7/
  57 +C
  58 + DATA A/1.162,22.344,18.50,2.602,6.903,5.287,0.5790,0.4462,0.7850/
  59 +C
  60 + DATA AM0,S0,X00,DSIG/70.,1.08,5.48,0.005/
  61 + DATA DELIMFX,DELIMFY /20.,10./
  62 +C
  63 + PDYN=PARMOD(1)
  64 + DST=PARMOD(2)
  65 + BYIMF=PARMOD(3)
  66 + BZIMF=PARMOD(4)
  67 +C
  68 + SPS=SIN(PS)
  69 + PPS=PS
  70 +C
  71 + DEPR=0.8*DST-13.*SQRT(PDYN) ! DEPR is an estimate of total near-Earth
  72 +c depression, based on DST and Pdyn
  73 +c (usually, DEPR &lt; 0 )
  74 +C
  75 +C CALCULATE THE IMF-RELATED QUANTITIES:
  76 +C
  77 + Bt=SQRT(BYIMF**2+BZIMF**2)
  78 +
  79 + IF (BYIMF.EQ.0..AND.BZIMF.EQ.0.) THEN
  80 + THETA=0.
  81 + GOTO 1
  82 + ENDIF
  83 +C
  84 + THETA=ATAN2(BYIMF,BZIMF)
  85 + IF (THETA.LE.0.D0) THETA=THETA+6.2831853
  86 + 1 CT=COS(THETA)
  87 + ST=SIN(THETA)
  88 + EPS=718.5*SQRT(Pdyn)*Bt*SIN(THETA/2.)
  89 +C
  90 + FACTEPS=EPS/EPS10-1.
  91 + FACTPD=SQRT(PDYN/PDYN0)-1.
  92 +C
  93 + RCAMPL=-A(1)*DEPR ! RCAMPL is the amplitude of the ring current
  94 +c (positive and equal to abs.value of RC depression at origin)
  95 +C
  96 + TAMPL2=A(2)+A(3)*FACTPD+A(4)*FACTEPS
  97 + TAMPL3=A(5)+A(6)*FACTPD
  98 + B1AMPL=A(7)+A(8)*FACTEPS
  99 + B2AMPL=20.*B1AMPL ! IT IS EQUIVALENT TO ASSUMING THAT THE TOTAL CURRENT
  100 +C IN THE REGION 2 SYSTEM IS 40% OF THAT IN REGION 1
  101 + RECONN=A(9)
  102 +C
  103 + XAPPA=(PDYN/PDYN0)**0.14
  104 + XAPPA3=XAPPA**3
  105 + YS=Y*CT-Z*ST
  106 + ZS=Z*CT+Y*ST
  107 +C
  108 + FACTIMF=EXP(X/DELIMFX-(YS/DELIMFY)**2)
  109 +C
  110 +C CALCULATE THE "IMF" COMPONENTS OUTSIDE THE LAYER (HENCE BEGIN WITH "O")
  111 +C
  112 + OIMFX=0.
  113 + OIMFY=RECONN*BYIMF*FACTIMF
  114 + OIMFZ=RECONN*BZIMF*FACTIMF
  115 +C
  116 + RIMFAMPL=RECONN*Bt
  117 +C
  118 + PPS=PS
  119 + XX=X*XAPPA
  120 + YY=Y*XAPPA
  121 + ZZ=Z*XAPPA
  122 +C
  123 +C SCALE AND CALCULATE THE MAGNETOPAUSE PARAMETERS FOR THE INTERPOLATION ACROSS
  124 +C THE BOUNDARY LAYER (THE COORDINATES XX,YY,ZZ ARE ALREADY SCALED)
  125 +C
  126 + X0=X00/XAPPA
  127 + AM=AM0/XAPPA
  128 + RHO2=Y**2+Z**2
  129 + ASQ=AM**2
  130 + XMXM=AM+X-X0
  131 + IF (XMXM.LT.0.) XMXM=0. ! THE BOUNDARY IS A CYLINDER TAILWARD OF X=X0-AM
  132 + AXX0=XMXM**2
  133 + ARO=ASQ+RHO2
  134 + SIGMA=SQRT((ARO+AXX0+SQRT((ARO+AXX0)**2-4.*ASQ*AXX0))/(2.*ASQ))
  135 +C
  136 +C NOW, THERE ARE THREE POSSIBLE CASES:
  137 +C (1) INSIDE THE MAGNETOSPHERE
  138 +C (2) IN THE BOUNDARY LAYER
  139 +C (3) OUTSIDE THE MAGNETOSPHERE AND B.LAYER
  140 +C FIRST OF ALL, CONSIDER THE CASES (1) AND (2):
  141 +C
  142 + IF (SIGMA.LT.S0+DSIG) THEN ! CALCULATE THE T95_06 FIELD (WITH THE
  143 +C POTENTIAL "PENETRATED" INTERCONNECTION FIELD):
  144 +
  145 + CALL DIPSHLD(PPS,XX,YY,ZZ,CFX,CFY,CFZ)
  146 + CALL TAILRC96(SPS,XX,YY,ZZ,BXRC,BYRC,BZRC,BXT2,BYT2,BZT2,
  147 + * BXT3,BYT3,BZT3)
  148 + CALL BIRK1TOT_02(PPS,XX,YY,ZZ,R1X,R1Y,R1Z)
  149 + CALL BIRK2TOT_02(PPS,XX,YY,ZZ,R2X,R2Y,R2Z)
  150 + CALL INTERCON(XX,YS*XAPPA,ZS*XAPPA,RIMFX,RIMFYS,RIMFZS)
  151 + RIMFY=RIMFYS*CT+RIMFZS*ST
  152 + RIMFZ=RIMFZS*CT-RIMFYS*ST
  153 +C
  154 + FX=CFX*XAPPA3+RCAMPL*BXRC +TAMPL2*BXT2+TAMPL3*BXT3
  155 + * +B1AMPL*R1X +B2AMPL*R2X +RIMFAMPL*RIMFX
  156 + FY=CFY*XAPPA3+RCAMPL*BYRC +TAMPL2*BYT2+TAMPL3*BYT3
  157 + * +B1AMPL*R1Y +B2AMPL*R2Y +RIMFAMPL*RIMFY
  158 + FZ=CFZ*XAPPA3+RCAMPL*BZRC +TAMPL2*BZT2+TAMPL3*BZT3
  159 + * +B1AMPL*R1Z +B2AMPL*R2Z +RIMFAMPL*RIMFZ
  160 +C
  161 +C NOW, LET US CHECK WHETHER WE HAVE THE CASE (1). IF YES - WE ARE DONE:
  162 +C
  163 + IF (SIGMA.LT.S0-DSIG) THEN
  164 + BX=FX
  165 + BY=FY
  166 + BZ=FZ
  167 + ELSE ! THIS IS THE MOST COMPLEX CASE: WE ARE INSIDE
  168 +C THE INTERPOLATION REGION
  169 + FINT=0.5*(1.-(SIGMA-S0)/DSIG)
  170 + FEXT=0.5*(1.+(SIGMA-S0)/DSIG)
  171 +C
  172 + CALL DIPOLE(PS,X,Y,Z,QX,QY,QZ)
  173 + BX=(FX+QX)*FINT+OIMFX*FEXT -QX
  174 + BY=(FY+QY)*FINT+OIMFY*FEXT -QY
  175 + BZ=(FZ+QZ)*FINT+OIMFZ*FEXT -QZ
  176 +c
  177 + ENDIF ! THE CASES (1) AND (2) ARE EXHAUSTED; THE ONLY REMAINING
  178 +C POSSIBILITY IS NOW THE CASE (3):
  179 + ELSE
  180 + CALL DIPOLE(PS,X,Y,Z,QX,QY,QZ)
  181 + BX=OIMFX-QX
  182 + BY=OIMFY-QY
  183 + BZ=OIMFZ-QZ
  184 + ENDIF
  185 +C
  186 + RETURN
  187 + END
  188 +C=====================================================================
  189 +
  190 + SUBROUTINE DIPSHLD(PS,X,Y,Z,BX,BY,BZ)
  191 +C
  192 +C CALCULATES GSM COMPONENTS OF THE EXTERNAL MAGNETIC FIELD DUE TO
  193 +C SHIELDING OF THE EARTH'S DIPOLE ONLY
  194 +C
  195 + IMPLICIT REAL*8 (A-H,O-Z)
  196 + DIMENSION A1(12),A2(12)
  197 + DATA A1 /.24777,-27.003,-.46815,7.0637,-1.5918,-.90317E-01,57.522,
  198 + * 13.757,2.0100,10.458,4.5798,2.1695/
  199 + DATA A2/-.65385,-18.061,-.40457,-5.0995,1.2846,.78231E-01,39.592,
  200 + * 13.291,1.9970,10.062,4.5140,2.1558/
  201 +C
  202 + CPS=DCOS(PS)
  203 + SPS=DSIN(PS)
  204 + CALL CYLHARM(A1,X,Y,Z,HX,HY,HZ)
  205 + CALL CYLHAR1(A2,X,Y,Z,FX,FY,FZ)
  206 +C
  207 + BX=HX*CPS+FX*SPS
  208 + BY=HY*CPS+FY*SPS
  209 + BZ=HZ*CPS+FZ*SPS
  210 + RETURN
  211 + END
  212 +C
  213 +C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  214 +C
  215 +C THIS CODE YIELDS THE SHIELDING FIELD FOR THE PERPENDICULAR DIPOLE
  216 +C
  217 + SUBROUTINE CYLHARM( A, X,Y,Z, BX,BY,BZ)
  218 +C
  219 +C
  220 +C *** N.A. Tsyganenko *** Sept. 14-18, 1993; revised March 16, 1994 ***
  221 +C
  222 +C An approximation for the Chapman-Ferraro field by a sum of 6 cylin-
  223 +c drical harmonics (see pp. 97-113 in the brown GSFC notebook #1)
  224 +c
  225 +C Description of parameters:
  226 +C
  227 +C A - input vector containing model parameters;
  228 +C X,Y,Z - input GSM coordinates
  229 +C BX,BY,BZ - output GSM components of the shielding field
  230 +C - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  231 +C The 6 linear parameters A(1)-A(6) are amplitudes of the cylindrical harmonic
  232 +c terms.
  233 +c The 6 nonlinear parameters A(7)-A(12) are the corresponding scale lengths
  234 +C for each term (see GSFC brown notebook).
  235 +c
  236 +C - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  237 +C
  238 + IMPLICIT REAL * 8 (A - H, O - Z)
  239 +C
  240 + DIMENSION A(12)
  241 +C
  242 + RHO=DSQRT(Y**2+Z**2)
  243 + IF (RHO.LT.1.D-8) THEN
  244 + SINFI=1.D0
  245 + COSFI=0.D0
  246 + RHO=1.D-8
  247 + GOTO 1
  248 + ENDIF
  249 +C
  250 + SINFI=Z/RHO
  251 + COSFI=Y/RHO
  252 + 1 SINFI2=SINFI**2
  253 + SI2CO2=SINFI2-COSFI**2
  254 +C
  255 + BX=0.D0
  256 + BY=0.D0
  257 + BZ=0.D0
  258 +C
  259 + DO 11 I=1,3
  260 + DZETA=RHO/A(I+6)
  261 + XJ0=BES(DZETA,0)
  262 + XJ1=BES(DZETA,1)
  263 + XEXP=DEXP(X/A(I+6))
  264 + BX=BX-A(I)*XJ1*XEXP*SINFI
  265 + BY=BY+A(I)*(2.D0*XJ1/DZETA-XJ0)*XEXP*SINFI*COSFI
  266 + BZ=BZ+A(I)*(XJ1/DZETA*SI2CO2-XJ0*SINFI2)*XEXP
  267 + 11 CONTINUE
  268 +c
  269 + DO 12 I=4,6
  270 + DZETA=RHO/A(I+6)
  271 + XKSI=X/A(I+6)
  272 + XJ0=BES(DZETA,0)
  273 + XJ1=BES(DZETA,1)
  274 + XEXP=DEXP(XKSI)
  275 + BRHO=(XKSI*XJ0-(DZETA**2+XKSI-1.D0)*XJ1/DZETA)*XEXP*SINFI
  276 + BPHI=(XJ0+XJ1/DZETA*(XKSI-1.D0))*XEXP*COSFI
  277 + BX=BX+A(I)*(DZETA*XJ0+XKSI*XJ1)*XEXP*SINFI
  278 + BY=BY+A(I)*(BRHO*COSFI-BPHI*SINFI)
  279 + BZ=BZ+A(I)*(BRHO*SINFI+BPHI*COSFI)
  280 + 12 CONTINUE
  281 +C
  282 +c
  283 + RETURN
  284 + END
  285 +C
  286 +c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  287 +C
  288 +C THIS CODE YIELDS THE SHIELDING FIELD FOR THE PARALLEL DIPOLE
  289 +C
  290 + SUBROUTINE CYLHAR1(A, X,Y,Z, BX,BY,BZ)
  291 +C
  292 +C
  293 +C *** N.A. Tsyganenko *** Sept. 14-18, 1993; revised March 16, 1994 ***
  294 +C
  295 +C An approximation of the Chapman-Ferraro field by a sum of 6 cylin-
  296 +c drical harmonics (see pages 97-113 in the brown GSFC notebook #1)
  297 +c
  298 +C Description of parameters:
  299 +C
  300 +C A - input vector containing model parameters;
  301 +C X,Y,Z - input GSM coordinates,
  302 +C BX,BY,BZ - output GSM components of the shielding field
  303 +C
  304 +C - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  305 +C
  306 +C The 6 linear parameters A(1)-A(6) are amplitudes of the cylindrical
  307 +c harmonic terms.
  308 +c The 6 nonlinear parameters A(7)-A(12) are the corresponding scale
  309 +c lengths for each term (see GSFC brown notebook).
  310 +c
  311 +C - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  312 +C
  313 + IMPLICIT REAL * 8 (A - H, O - Z)
  314 +C
  315 + DIMENSION A(12)
  316 +C
  317 + RHO=DSQRT(Y**2+Z**2)
  318 + IF (RHO.LT.1.D-10) THEN
  319 + SINFI=1.D0
  320 + COSFI=0.D0
  321 + GOTO 1
  322 + ENDIF
  323 +C
  324 + SINFI=Z/RHO
  325 + COSFI=Y/RHO
  326 +C
  327 + 1 BX=0.D0
  328 + BY=0.D0
  329 + BZ=0.D0
  330 +C
  331 + DO 11 I=1,3
  332 + DZETA=RHO/A(I+6)
  333 + XKSI=X/A(I+6)
  334 + XJ0=BES(DZETA,0)
  335 + XJ1=BES(DZETA,1)
  336 + XEXP=DEXP(XKSI)
  337 + BRHO=XJ1*XEXP
  338 + BX=BX-A(I)*XJ0*XEXP
  339 + BY=BY+A(I)*BRHO*COSFI
  340 + BZ=BZ+A(I)*BRHO*SINFI
  341 + 11 CONTINUE
  342 +c
  343 + DO 12 I=4,6
  344 + DZETA=RHO/A(I+6)
  345 + XKSI=X/A(I+6)
  346 + XJ0=BES(DZETA,0)
  347 + XJ1=BES(DZETA,1)
  348 + XEXP=DEXP(XKSI)
  349 + BRHO=(DZETA*XJ0+XKSI*XJ1)*XEXP
  350 + BX=BX+A(I)*(DZETA*XJ1-XJ0*(XKSI+1.D0))*XEXP
  351 + BY=BY+A(I)*BRHO*COSFI
  352 + BZ=BZ+A(I)*BRHO*SINFI
  353 + 12 CONTINUE
  354 +C
  355 + RETURN
  356 + END
  357 +
  358 +c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  359 +C
  360 + DOUBLE PRECISION FUNCTION BES(X,K)
  361 + IMPLICIT REAL*8 (A-H,O-Z)
  362 +C
  363 + IF (K.EQ.0) THEN
  364 + BES=BES0(X)
  365 + RETURN
  366 + ENDIF
  367 +C
  368 + IF (K.EQ.1) THEN
  369 + BES=BES1(X)
  370 + RETURN
  371 + ENDIF
  372 +C
  373 + IF (X.EQ.0.D0) THEN
  374 + BES=0.D0
  375 + RETURN
  376 + ENDIF
  377 +C
  378 + G=2.D0/X
  379 + IF (X.LE.DFLOAT(K)) GOTO 10
  380 +C
  381 + N=1
  382 + XJN=BES1(X)
  383 + XJNM1=BES0(X)
  384 +C
  385 + 1 XJNP1=G*N*XJN-XJNM1
  386 + N=N+1
  387 + IF (N.LT.K) GOTO 2
  388 + BES=XJNP1
  389 + RETURN
  390 +C
  391 + 2 XJNM1=XJN
  392 + XJN=XJNP1
  393 + GOTO 1
  394 +C
  395 + 10 N=24
  396 + XJN=1.D0
  397 + XJNP1=0.D0
  398 + SUM=0.D0
  399 +C
  400 + 3 IF (MOD(N,2).EQ.0) SUM=SUM+XJN
  401 + XJNM1=G*N*XJN-XJNP1
  402 + N=N-1
  403 +C
  404 + XJNP1=XJN
  405 + XJN=XJNM1
  406 + IF (N.EQ.K) BES=XJN
  407 +C
  408 + IF (DABS(XJN).GT.1.D5) THEN
  409 + XJNP1=XJNP1*1.D-5
  410 + XJN=XJN*1.D-5
  411 + SUM=SUM*1.D-5
  412 + IF (N.LE.K) BES=BES*1.D-5
  413 + ENDIF
  414 +C
  415 + IF (N.EQ.0) GOTO 4
  416 + GOTO 3
  417 +C
  418 + 4 SUM=XJN+2.D0*SUM
  419 + BES=BES/SUM
  420 + RETURN
  421 + END
  422 +c-------------------------------------------------------------------
  423 +c
  424 + DOUBLE PRECISION FUNCTION BES0(X)
  425 +C
  426 + IMPLICIT REAL*8 (A-H,O-Z)
  427 +C
  428 + IF (DABS(X).LT.3.D0) THEN
  429 + X32=(X/3.D0)**2
  430 + BES0=1.D0-X32*(2.2499997D0-X32*(1.2656208D0-X32*
  431 + * (0.3163866D0-X32*(0.0444479D0-X32*(0.0039444D0
  432 + * -X32*0.00021D0)))))
  433 + ELSE
  434 + XD3=3.D0/X
  435 + F0=0.79788456D0-XD3*(0.00000077D0+XD3*(0.00552740D0+XD3*
  436 + * (0.00009512D0-XD3*(0.00137237D0-XD3*(0.00072805D0
  437 + * -XD3*0.00014476D0)))))
  438 + T0=X-0.78539816D0-XD3*(0.04166397D0+XD3*(0.00003954D0-XD3*
  439 + * (0.00262573D0-XD3*(0.00054125D0+XD3*(0.00029333D0
  440 + * -XD3*0.00013558D0)))))
  441 + BES0=F0/DSQRT(X)*DCOS(T0)
  442 + ENDIF
  443 + RETURN
  444 + END
  445 +c
  446 +c--------------------------------------------------------------------------
  447 +c
  448 + DOUBLE PRECISION FUNCTION BES1(X)
  449 +C
  450 + IMPLICIT REAL*8 (A-H,O-Z)
  451 +C
  452 + IF (DABS(X).LT.3.D0) THEN
  453 + X32=(X/3.D0)**2
  454 + BES1XM1=0.5D0-X32*(0.56249985D0-X32*(0.21093573D0-X32*
  455 + * (0.03954289D0-X32*(0.00443319D0-X32*(0.00031761D0
  456 + * -X32*0.00001109D0)))))
  457 + BES1=BES1XM1*X
  458 + ELSE
  459 + XD3=3.D0/X
  460 + F1=0.79788456D0+XD3*(0.00000156D0+XD3*(0.01659667D0+XD3*
  461 + * (0.00017105D0-XD3*(0.00249511D0-XD3*(0.00113653D0
  462 + * -XD3*0.00020033D0)))))
  463 + T1=X-2.35619449D0+XD3*(0.12499612D0+XD3*(0.0000565D0-XD3*
  464 + * (0.00637879D0-XD3*(0.00074348D0+XD3*(0.00079824D0
  465 + * -XD3*0.00029166D0)))))
  466 + BES1=F1/DSQRT(X)*DCOS(T1)
  467 + ENDIF
  468 + RETURN
  469 + END
  470 +C------------------------------------------------------------
  471 +C
  472 + SUBROUTINE INTERCON(X,Y,Z,BX,BY,BZ)
  473 +C
  474 +C Calculates the potential interconnection field inside the magnetosphere,
  475 +c corresponding to DELTA_X = 20Re and DELTA_Y = 10Re (NB#3, p.90, 6/6/1996).
  476 +C The position (X,Y,Z) and field components BX,BY,BZ are given in the rotated
  477 +c coordinate system, in which the Z-axis is always directed along the BzIMF
  478 +c (i.e. rotated by the IMF clock angle Theta)
  479 +C It is also assumed that the IMF Bt=1, so that the components should be
  480 +c (i) multiplied by the actual Bt, and
  481 +c (ii) transformed to standard GSM coords by rotating back around X axis
  482 +c by the angle -Theta.
  483 +c
  484 +C Description of parameters:
  485 +C
  486 +C X,Y,Z - GSM POSITION
  487 +C BX,BY,BZ - INTERCONNECTION FIELD COMPONENTS INSIDE THE MAGNETOSPHERE
  488 +C OF A STANDARD SIZE (TO TAKE INTO ACCOUNT EFFECTS OF PRESSURE CHANGES,
  489 +C APPLY THE SCALING TRANSFORMATION)
  490 +C
  491 +C - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  492 +C
  493 +C The 9 linear parameters are amplitudes of the "cartesian" harmonics
  494 +c The 6 nonlinear parameters are the scales Pi and Ri entering
  495 +c the arguments of exponents, sines, and cosines in the 9 "Cartesian"
  496 +c harmonics (3+3)
  497 +C - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  498 +C
  499 + IMPLICIT REAL * 8 (A - H, O - Z)
  500 +C
  501 + DIMENSION A(15),RP(3),RR(3),P(3),R(3)
  502 + SAVE
  503 +C
  504 + DATA A/-8.411078731,5932254.951,-9073284.93,-11.68794634,
  505 + * 6027598.824,-9218378.368,-6.508798398,-11824.42793,18015.66212,
  506 + * 7.99754043,13.9669886,90.24475036,16.75728834,1015.645781,
  507 + * 1553.493216/
  508 +C
  509 + DATA M/0/
  510 +C
  511 + IF (M.NE.0) GOTO 111
  512 + M=1
  513 +C
  514 + P(1)=A(10)
  515 + P(2)=A(11)
  516 + P(3)=A(12)
  517 + R(1)=A(13)
  518 + R(2)=A(14)
  519 + R(3)=A(15)
  520 +C
  521 +C
  522 + DO 11 I=1,3
  523 + RP(I)=1.D0/P(I)
  524 + 11 RR(I)=1.D0/R(I)
  525 +C
  526 + 111 CONTINUE
  527 +C
  528 + L=0
  529 +C
  530 + BX=0.
  531 + BY=0.
  532 + BZ=0.
  533 +C
  534 +c "PERPENDICULAR" KIND OF SYMMETRY ONLY
  535 +C
  536 + DO 2 I=1,3
  537 + CYPI=DCOS(Y*RP(I))
  538 + SYPI=DSIN(Y*RP(I))
  539 +C
  540 + DO 2 K=1,3
  541 + SZRK=DSIN(Z*RR(K))
  542 + CZRK=DCOS(Z*RR(K))
  543 + SQPR=DSQRT(RP(I)**2+RR(K)**2)
  544 + EPR=DEXP(X*SQPR)
  545 +C
  546 + HX=-SQPR*EPR*CYPI*SZRK
  547 + HY=RP(I)*EPR*SYPI*SZRK
  548 + HZ=-RR(K)*EPR*CYPI*CZRK
  549 + L=L+1
  550 +c
  551 + BX=BX+A(L)*HX
  552 + BY=BY+A(L)*HY
  553 + BZ=BZ+A(L)*HZ
  554 + 2 CONTINUE
  555 +C
  556 + RETURN
  557 + END
  558 +
  559 +C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  560 + SUBROUTINE TAILRC96(SPS,X,Y,Z,BXRC,BYRC,BZRC,BXT2,BYT2,BZT2,
  561 + * BXT3,BYT3,BZT3)
  562 +c
  563 +c COMPUTES THE COMPONENTS OF THE FIELD OF THE MODEL RING CURRENT AND THREE
  564 +c TAIL MODES WITH UNIT AMPLITUDES
  565 +C (FOR THE RING CURRENT, IT MEANS THE DISTURBANCE OF Bz=-1nT AT ORIGIN,
  566 +C AND FOR THE TAIL MODES IT MEANS MAXIMAL BX JUST ABOVE THE SHEET EQUAL 1 nT.
  567 +C
  568 + IMPLICIT REAL*8 (A-H,O-Z)
  569 + DIMENSION ARC(48),ATAIL2(48),ATAIL3(48)
  570 + COMMON /WARP/ CPSS,SPSS,DPSRR,RPS,WARP,D,XS,ZS,DXSX,DXSY,DXSZ,
  571 + * DZSX,DZSY,DZSZ,DZETAS,DDZETADX,DDZETADY,DDZETADZ,ZSWW
  572 +C
  573 + DATA ARC/-3.087699646,3.516259114,18.81380577,-13.95772338,
  574 + * -5.497076303,0.1712890838,2.392629189,-2.728020808,-14.79349936,
  575 + * 11.08738083,4.388174084,0.2492163197E-01,0.7030375685,
  576 + *-.7966023165,-3.835041334,2.642228681,-0.2405352424,-0.7297705678,
  577 + * -0.3680255045,0.1333685557,2.795140897,-1.078379954,0.8014028630,
  578 + * 0.1245825565,0.6149982835,-0.2207267314,-4.424578723,1.730471572,
  579 + * -1.716313926,-0.2306302941,-0.2450342688,0.8617173961E-01,
  580 + * 1.54697858,-0.6569391113,-0.6537525353,0.2079417515,12.75434981,
  581 + * 11.37659788,636.4346279,1.752483754,3.604231143,12.83078674,
  582 + * 7.412066636,9.434625736,676.7557193,1.701162737,3.580307144,
  583 + * 14.64298662/
  584 +C
  585 + DATA ATAIL2/.8747515218,-.9116821411,2.209365387,-2.159059518,
  586 + * -7.059828867,5.924671028,-1.916935691,1.996707344,-3.877101873,
  587 + * 3.947666061,11.38715899,-8.343210833,1.194109867,-1.244316975,
  588 + * 3.73895491,-4.406522465,-20.66884863,3.020952989,.2189908481,
  589 + * -.09942543549,-.927225562,.1555224669,.6994137909,-.08111721003,
  590 + * -.7565493881,.4686588792,4.266058082,-.3717470262,-3.920787807,
  591 + * .02298569870,.7039506341,-.5498352719,-6.675140817,.8279283559,
  592 + * -2.234773608,-1.622656137,5.187666221,6.802472048,39.13543412,
  593 + * 2.784722096,6.979576616,25.71716760,4.495005873,8.068408272,
  594 + * 93.47887103,4.158030104,9.313492566,57.18240483/
  595 +C
  596 + DATA ATAIL3/-19091.95061,-3011.613928,20582.16203,4242.918430,
  597 + * -2377.091102,-1504.820043,19884.04650,2725.150544,-21389.04845,
  598 + * -3990.475093,2401.610097,1548.171792,-946.5493963,490.1528941,
  599 + * 986.9156625,-489.3265930,-67.99278499,8.711175710,-45.15734260,
  600 + * -10.76106500,210.7927312,11.41764141,-178.0262808,.7558830028,
  601 + * 339.3806753,9.904695974,69.50583193,-118.0271581,22.85935896,
  602 + * 45.91014857,-425.6607164,15.47250738,118.2988915,65.58594397,
  603 + * -201.4478068,-14.57062940,19.69877970,20.30095680,86.45407420,
  604 + * 22.50403727,23.41617329,48.48140573,24.61031329,123.5395974,
  605 + * 223.5367692,39.50824342,65.83385762,266.2948657/
  606 +C
  607 + DATA RH,DR,G,D0,DELTADY/9.,4.,10.,2.,10./
  608 +C
  609 +C TO ECONOMIZE THE CODE, WE FIRST CALCULATE COMMON VARIABLES, WHICH ARE
  610 +C THE SAME FOR ALL MODES, AND PUT THEM IN THE COMMON-BLOCK /WARP/
  611 +C
  612 + DR2=DR*DR
  613 + C11=DSQRT((1.D0+RH)**2+DR2)
  614 + C12=DSQRT((1.D0-RH)**2+DR2)
  615 + C1=C11-C12
  616 + SPSC1=SPS/C1
  617 + RPS=0.5*(C11+C12)*SPS ! THIS IS THE SHIFT OF OF THE SHEET WITH RESPECT
  618 +C TO GSM EQ.PLANE FOR THE 3RD (ASYMPTOTIC) TAIL MODE
  619 +C
  620 + R=DSQRT(X*X+Y*Y+Z*Z)
  621 + SQ1=DSQRT((R+RH)**2+DR2)
  622 + SQ2=DSQRT((R-RH)**2+DR2)
  623 + C=SQ1-SQ2
  624 + CS=(R+RH)/SQ1-(R-RH)/SQ2
  625 + SPSS=SPSC1/R*C
  626 + CPSS=DSQRT(1.D0-SPSS**2)
  627 + DPSRR=SPS/(R*R)*(CS*R-C)/DSQRT((R*C1)**2-(C*SPS)**2)
  628 +C
  629 + WFAC=Y/(Y**4+1.D4) ! WARPING
  630 + W=WFAC*Y**3
  631 + WS=4.D4*Y*WFAC**2
  632 + WARP=G*SPS*W
  633 + XS=X*CPSS-Z*SPSS
  634 + ZSWW=Z*CPSS+X*SPSS ! "WW" MEANS "WITHOUT Y-Z WARPING" (IN X-Z ONLY)
  635 + ZS=ZSWW +WARP
  636 +
  637 + DXSX=CPSS-X*ZSWW*DPSRR
  638 + DXSY=-Y*ZSWW*DPSRR
  639 + DXSZ=-SPSS-Z*ZSWW*DPSRR
  640 + DZSX=SPSS+X*XS*DPSRR
  641 + DZSY=XS*Y*DPSRR +G*SPS*WS ! THE LAST TERM IS FOR THE Y-Z WARP
  642 + DZSZ=CPSS+XS*Z*DPSRR ! (TAIL MODES ONLY)
  643 +
  644 + D=D0+DELTADY*(Y/20.D0)**2 ! SHEET HALF-THICKNESS FOR THE TAIL MODES
  645 + DDDY=DELTADY*Y*0.005D0 ! (THICKENS TO FLANKS, BUT NO VARIATION
  646 +C ALONG X, IN CONTRAST TO RING CURRENT)
  647 +C
  648 + DZETAS=DSQRT(ZS**2+D**2) ! THIS IS THE SAME SIMPLE WAY TO SPREAD
  649 +C OUT THE SHEET, AS THAT USED IN T89
  650 + DDZETADX=ZS*DZSX/DZETAS
  651 + DDZETADY=(ZS*DZSY+D*DDDY)/DZETAS
  652 + DDZETADZ=ZS*DZSZ/DZETAS
  653 +C
  654 + CALL SHLCAR3X3(ARC,X,Y,Z,SPS,WX,WY,WZ)
  655 + CALL RINGCURR96(X,Y,Z,HX,HY,HZ)
  656 + BXRC=WX+HX
  657 + BYRC=WY+HY
  658 + BZRC=WZ+HZ
  659 +C
  660 + CALL SHLCAR3X3(ATAIL2,X,Y,Z,SPS,WX,WY,WZ)
  661 + CALL TAILDISK(X,Y,Z,HX,HY,HZ)
  662 + BXT2=WX+HX
  663 + BYT2=WY+HY
  664 + BZT2=WZ+HZ
  665 +C
  666 + CALL SHLCAR3X3(ATAIL3,X,Y,Z,SPS,WX,WY,WZ)
  667 + CALL TAIL87(X,Z,HX,HZ)
  668 + BXT3=WX+HX
  669 + BYT3=WY
  670 + BZT3=WZ+HZ
  671 +C
  672 + RETURN
  673 + END
  674 +C
  675 +c********************************************************************
  676 +C
  677 + SUBROUTINE RINGCURR96(X,Y,Z,BX,BY,BZ)
  678 +c
  679 +c THIS SUBROUTINE COMPUTES THE COMPONENTS OF THE RING CURRENT FIELD,
  680 +C SIMILAR TO THAT DESCRIBED BY TSYGANENKO AND PEREDO (1994). THE
  681 +C DIFFERENCE IS THAT NOW WE USE SPACEWARPING, AS DESCRIBED IN THE
  682 +C PAPER ON MODELING BIRKELAND CURRENTS (TSYGANENKO AND STERN, 1996),
  683 +C INSTEAD OF SHEARING IT IN THE SPIRIT OF THE T89 TAIL MODEL.
  684 +C
  685 +C IN ADDITION, INSTEAD OF 7 TERMS FOR THE RING CURRENT MODEL, WE USE
  686 +C NOW ONLY 2 TERMS; THIS SIMPLIFICATION ALSO GIVES RISE TO AN
  687 +C EASTWARD RING CURRENT LOCATED EARTHWARD FROM THE MAIN ONE,
  688 +C IN LINE WITH WHAT IS ACTUALLY OBSERVED
  689 +C
  690 +C FOR DETAILS, SEE NB #3, PAGES 70-73
  691 +C
  692 + IMPLICIT REAL*8 (A-H,O-Z)
  693 + DIMENSION F(2),BETA(2)
  694 + COMMON /WARP/ CPSS,SPSS,DPSRR, XNEXT(3),XS,ZSWARPED,DXSX,DXSY,
  695 + * DXSZ,DZSX,DZSYWARPED,DZSZ,OTHER(4),ZS ! ZS HERE IS WITHOUT Y-Z WARP
  696 +C
  697 +
  698 + DATA D0,DELTADX,XD,XLDX /2.,0.,0.,4./ ! ACHTUNG !! THE RC IS NOW
  699 +C COMPLETELY SYMMETRIC (DELTADX=0)
  700 +
  701 +C
  702 + DATA F,BETA /569.895366D0,-1603.386993D0,2.722188D0,3.766875D0/
  703 +C
  704 +C THE ORIGINAL VALUES OF F(I) WERE MULTIPLIED BY BETA(I) (TO REDUCE THE
  705 +C NUMBER OF MULTIPLICATIONS BELOW) AND BY THE FACTOR -0.43, NORMALIZING
  706 +C THE DISTURBANCE AT ORIGIN TO B=-1nT
  707 +C
  708 + DZSY=XS*Y*DPSRR ! NO WARPING IN THE Y-Z PLANE (ALONG X ONLY), AND
  709 +C THIS IS WHY WE DO NOT USE DZSY FROM THE COMMON-BLOCK
  710 + XXD=X-XD
  711 + FDX=0.5D0*(1.D0+XXD/DSQRT(XXD**2+XLDX**2))
  712 + DDDX=DELTADX*0.5D0*XLDX**2/DSQRT(XXD**2+XLDX**2)**3
  713 + D=D0+DELTADX*FDX
  714 +
  715 + DZETAS=DSQRT(ZS**2+D**2) ! THIS IS THE SAME SIMPLE WAY TO SPREAD
  716 +C OUT THE SHEET, AS THAT USED IN T89
  717 + RHOS=DSQRT(XS**2+Y**2)
  718 + DDZETADX=(ZS*DZSX+D*DDDX)/DZETAS
  719 + DDZETADY=ZS*DZSY/DZETAS
  720 + DDZETADZ=ZS*DZSZ/DZETAS
  721 + IF (RHOS.LT.1.D-5) THEN
  722 + DRHOSDX=0.D0
  723 + DRHOSDY=DSIGN(1.D0,Y)
  724 + DRHOSDZ=0.D0
  725 + ELSE
  726 + DRHOSDX=XS*DXSX/RHOS
  727 + DRHOSDY=(XS*DXSY+Y)/RHOS
  728 + DRHOSDZ=XS*DXSZ/RHOS
  729 + ENDIF
  730 +C
  731 + BX=0.D0
  732 + BY=0.D0
  733 + BZ=0.D0
  734 +C
  735 + DO 1 I=1,2
  736 +C
  737 + BI=BETA(I)
  738 +C
  739 + S1=DSQRT((DZETAS+BI)**2+(RHOS+BI)**2)
  740 + S2=DSQRT((DZETAS+BI)**2+(RHOS-BI)**2)
  741 + DS1DDZ=(DZETAS+BI)/S1
  742 + DS2DDZ=(DZETAS+BI)/S2
  743 + DS1DRHOS=(RHOS+BI)/S1
  744 + DS2DRHOS=(RHOS-BI)/S2
  745 +C
  746 + DS1DX=DS1DDZ*DDZETADX+DS1DRHOS*DRHOSDX
  747 + DS1DY=DS1DDZ*DDZETADY+DS1DRHOS*DRHOSDY
  748 + DS1DZ=DS1DDZ*DDZETADZ+DS1DRHOS*DRHOSDZ
  749 +C
  750 + DS2DX=DS2DDZ*DDZETADX+DS2DRHOS*DRHOSDX
  751 + DS2DY=DS2DDZ*DDZETADY+DS2DRHOS*DRHOSDY
  752 + DS2DZ=DS2DDZ*DDZETADZ+DS2DRHOS*DRHOSDZ
  753 +C
  754 + S1TS2=S1*S2
  755 + S1PS2=S1+S2
  756 + S1PS2SQ=S1PS2**2
  757 + FAC1=DSQRT(S1PS2SQ-(2.D0*BI)**2)
  758 + AS=FAC1/(S1TS2*S1PS2SQ)
  759 + TERM1=1.D0/(S1TS2*S1PS2*FAC1)
  760 + FAC2=AS/S1PS2SQ
  761 + DASDS1=TERM1-FAC2/S1*(S2*S2+S1*(3.D0*S1+4.D0*S2))
  762 + DASDS2=TERM1-FAC2/S2*(S1*S1+S2*(3.D0*S2+4.D0*S1))
  763 +C
  764 + DASDX=DASDS1*DS1DX+DASDS2*DS2DX
  765 + DASDY=DASDS1*DS1DY+DASDS2*DS2DY
  766 + DASDZ=DASDS1*DS1DZ+DASDS2*DS2DZ
  767 +C
  768 + BX=BX+F(I)*((2.D0*AS+Y*DASDY)*SPSS-XS*DASDZ
  769 + * +AS*DPSRR*(Y**2*CPSS+Z*ZS))
  770 + BY=BY-F(I)*Y*(AS*DPSRR*XS+DASDZ*CPSS+DASDX*SPSS)
  771 + 1 BZ=BZ+F(I)*((2.D0*AS+Y*DASDY)*CPSS+XS*DASDX
  772 + * -AS*DPSRR*(X*ZS+Y**2*SPSS))
  773 +C
  774 + RETURN
  775 + END
  776 +C
  777 +C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  778 +C
  779 + SUBROUTINE TAILDISK(X,Y,Z,BX,BY,BZ)
  780 +C
  781 +c
  782 +c THIS SUBROUTINE COMPUTES THE COMPONENTS OF THE TAIL CURRENT FIELD,
  783 +C SIMILAR TO THAT DESCRIBED BY TSYGANENKO AND PEREDO (1994). THE
  784 +C DIFFERENCE IS THAT NOW WE USE SPACEWARPING, AS DESCRIBED IN OUR
  785 +C PAPER ON MODELING BIRKELAND CURRENTS (TSYGANENKO AND STERN, 1996)
  786 +C INSTEAD OF SHEARING IT IN THE SPIRIT OF T89 TAIL MODEL.
  787 +C
  788 +C IN ADDITION, INSTEAD OF 8 TERMS FOR THE TAIL CURRENT MODEL, WE USE
  789 +C NOW ONLY 4 TERMS
  790 +C
  791 +C FOR DETAILS, SEE NB #3, PAGES 74-
  792 +C
  793 + IMPLICIT REAL*8 (A-H,O-Z)
  794 + DIMENSION F(4),BETA(4)
  795 + COMMON /WARP/ CPSS,SPSS,DPSRR,XNEXT(3),XS,ZS,DXSX,DXSY,DXSZ,
  796 + * OTHER(3),DZETAS,DDZETADX,DDZETADY,DDZETADZ,ZSWW
  797 +C
  798 + DATA XSHIFT /4.5/
  799 +C
  800 + DATA F,BETA
  801 + * / -745796.7338D0,1176470.141D0,-444610.529D0,-57508.01028D0,
  802 + * 7.9250000D0,8.0850000D0,8.4712500D0,27.89500D0/
  803 +c
  804 +c here original F(I) are multiplied by BETA(I), to economize
  805 +c calculations
  806 +C
  807 + RHOS=DSQRT((XS-XSHIFT)**2+Y**2)
  808 + IF (RHOS.LT.1.D-5) THEN
  809 + DRHOSDX=0.D0
  810 + DRHOSDY=DSIGN(1.D0,Y)
  811 + DRHOSDZ=0.D0
  812 + ELSE
  813 + DRHOSDX=(XS-XSHIFT)*DXSX/RHOS
  814 + DRHOSDY=((XS-XSHIFT)*DXSY+Y)/RHOS
  815 + DRHOSDZ=(XS-XSHIFT)*DXSZ/RHOS
  816 + ENDIF
  817 +C
  818 + BX=0.D0
  819 + BY=0.D0
  820 + BZ=0.D0
  821 +C
  822 + DO 1 I=1,4
  823 +C
  824 + BI=BETA(I)
  825 +C
  826 + S1=DSQRT((DZETAS+BI)**2+(RHOS+BI)**2)
  827 + S2=DSQRT((DZETAS+BI)**2+(RHOS-BI)**2)
  828 + DS1DDZ=(DZETAS+BI)/S1
  829 + DS2DDZ=(DZETAS+BI)/S2
  830 + DS1DRHOS=(RHOS+BI)/S1
  831 + DS2DRHOS=(RHOS-BI)/S2
  832 +C
  833 + DS1DX=DS1DDZ*DDZETADX+DS1DRHOS*DRHOSDX
  834 + DS1DY=DS1DDZ*DDZETADY+DS1DRHOS*DRHOSDY
  835 + DS1DZ=DS1DDZ*DDZETADZ+DS1DRHOS*DRHOSDZ
  836 +C
  837 + DS2DX=DS2DDZ*DDZETADX+DS2DRHOS*DRHOSDX
  838 + DS2DY=DS2DDZ*DDZETADY+DS2DRHOS*DRHOSDY
  839 + DS2DZ=DS2DDZ*DDZETADZ+DS2DRHOS*DRHOSDZ
  840 +C
  841 + S1TS2=S1*S2
  842 + S1PS2=S1+S2
  843 + S1PS2SQ=S1PS2**2
  844 + FAC1=DSQRT(S1PS2SQ-(2.D0*BI)**2)
  845 + AS=FAC1/(S1TS2*S1PS2SQ)
  846 + TERM1=1.D0/(S1TS2*S1PS2*FAC1)
  847 + FAC2=AS/S1PS2SQ
  848 + DASDS1=TERM1-FAC2/S1*(S2*S2+S1*(3.D0*S1+4.D0*S2))
  849 + DASDS2=TERM1-FAC2/S2*(S1*S1+S2*(3.D0*S2+4.D0*S1))
  850 +C
  851 + DASDX=DASDS1*DS1DX+DASDS2*DS2DX
  852 + DASDY=DASDS1*DS1DY+DASDS2*DS2DY
  853 + DASDZ=DASDS1*DS1DZ+DASDS2*DS2DZ
  854 +C
  855 + BX=BX+F(I)*((2.D0*AS+Y*DASDY)*SPSS-(XS-XSHIFT)*DASDZ
  856 + * +AS*DPSRR*(Y**2*CPSS+Z*ZSWW))
  857 +C
  858 + BY=BY-F(I)*Y*(AS*DPSRR*XS+DASDZ*CPSS+DASDX*SPSS)
  859 + 1 BZ=BZ+F(I)*((2.D0*AS+Y*DASDY)*CPSS+(XS-XSHIFT)*DASDX
  860 + * -AS*DPSRR*(X*ZSWW+Y**2*SPSS))
  861 +
  862 + RETURN
  863 + END
  864 +
  865 +C-------------------------------------------------------------------------
  866 +C
  867 + SUBROUTINE TAIL87(X,Z,BX,BZ)
  868 +
  869 + IMPLICIT REAL*8 (A-H,O-Z)
  870 +
  871 + COMMON /WARP/ FIRST(3), RPS,WARP,D, OTHER(13)
  872 +C
  873 +C 'LONG' VERSION OF THE 1987 TAIL MAGNETIC FIELD MODEL
  874 +C (N.A.TSYGANENKO, PLANET. SPACE SCI., V.35, P.1347, 1987)
  875 +C
  876 +C D IS THE Y-DEPENDENT SHEET HALF-THICKNESS (INCREASING TOWARDS FLANKS)
  877 +C RPS IS THE TILT-DEPENDENT SHIFT OF THE SHEET IN THE Z-DIRECTION,
  878 +C CORRESPONDING TO THE ASYMPTOTIC HINGING DISTANCE, DEFINED IN THE
  879 +C MAIN SUBROUTINE (TAILRC96) FROM THE PARAMETERS RH AND DR OF THE
  880 +C T96-TYPE MODULE, AND
  881 +C WARP IS THE BENDING OF THE SHEET FLANKS IN THE Z-DIRECTION, DIRECTED
  882 +C OPPOSITE TO RPS, AND INCREASING WITH DIPOLE TILT AND |Y|
  883 +C
  884 +
  885 + DATA DD/3./
  886 +C
  887 + DATA HPI,RT,XN,X1,X2,B0,B1,B2,XN21,XNR,ADLN
  888 + * /1.5707963,40.,-10.,
  889 + * -1.261,-0.663,0.391734,5.89715,24.6833,76.37,-0.1071,0.13238005/
  890 +C !!! THESE ARE NEW VALUES OF X1, X2, B0, B1, B2,
  891 +C CORRESPONDING TO TSCALE=1, INSTEAD OF TSCALE=0.6
  892 +C
  893 +C THE ABOVE QUANTITIES WERE DEFINED AS FOLLOWS:------------------------
  894 +C HPI=PI/2
  895 +C RT=40. ! Z-POSITION OF UPPER AND LOWER ADDITIONAL SHEETS
  896 +C XN=-10. ! INNER EDGE POSITION
  897 +C
  898 +C TSCALE=1 ! SCALING FACTOR, DEFINING THE RATE OF INCREASE OF THE
  899 +C CURRENT DENSITY TAILWARDS
  900 +C
  901 +c ATTENTION ! NOW I HAVE CHANGED TSCALE TO: TSCALE=1.0, INSTEAD OF 0.6
  902 +c OF THE PREVIOUS VERSION
  903 +c
  904 +C B0=0.391734
  905 +C B1=5.89715 *TSCALE
  906 +C B2=24.6833 *TSCALE**2
  907 +C
  908 +C HERE ORIGINAL VALUES OF THE MODE AMPLITUDES (P.77, NB#3) WERE NORMALIZED
  909 +C SO THAT ASYMPTOTIC BX=1 AT X=-200RE
  910 +C
  911 +C X1=(4.589 -5.85) *TSCALE -(TSCALE-1.)*XN ! NONLINEAR PARAMETERS OF THE
  912 +C CURRENT FUNCTION
  913 +C X2=(5.187 -5.85) *TSCALE -(TSCALE-1.)*XN
  914 +c
  915 +c
  916 +C XN21=(XN-X1)**2
  917 +C XNR=1./(XN-X2)
  918 +C ADLN=-DLOG(XNR**2*XN21)
  919 +C
  920 +C---------------------------------------------------------------
  921 +C
  922 + ZS=Z -RPS +WARP
  923 + ZP=Z-RT
  924 + ZM=Z+RT
  925 +C
  926 + XNX=XN-X
  927 + XNX2=XNX**2
  928 + XC1=X-X1
  929 + XC2=X-X2
  930 + XC22=XC2**2
  931 + XR2=XC2*XNR
  932 + XC12=XC1**2
  933 + D2=DD**2 ! SQUARE OF THE TOTAL HALFTHICKNESS (DD=3Re for this mode)
  934 + B20=ZS**2+D2
  935 + B2P=ZP**2+D2
  936 + B2M=ZM**2+D2
  937 + B=DSQRT(B20)
  938 + BP=DSQRT(B2P)
  939 + BM=DSQRT(B2M)
  940 + XA1=XC12+B20
  941 + XAP1=XC12+B2P
  942 + XAM1=XC12+B2M
  943 + XA2=1./(XC22+B20)
  944 + XAP2=1./(XC22+B2P)
  945 + XAM2=1./(XC22+B2M)
  946 + XNA=XNX2+B20
  947 + XNAP=XNX2+B2P
  948 + XNAM=XNX2+B2M
  949 + F=B20-XC22
  950 + FP=B2P-XC22
  951 + FM=B2M-XC22
  952 + XLN1=DLOG(XN21/XNA)
  953 + XLNP1=DLOG(XN21/XNAP)
  954 + XLNM1=DLOG(XN21/XNAM)
  955 + XLN2=XLN1+ADLN
  956 + XLNP2=XLNP1+ADLN
  957 + XLNM2=XLNM1+ADLN
  958 + ALN=0.25*(XLNP1+XLNM1-2.*XLN1)
  959 + S0=(DATAN(XNX/B)+HPI)/B
  960 + S0P=(DATAN(XNX/BP)+HPI)/BP
  961 + S0M=(DATAN(XNX/BM)+HPI)/BM
  962 + S1=(XLN1*.5+XC1*S0)/XA1
  963 + S1P=(XLNP1*.5+XC1*S0P)/XAP1
  964 + S1M=(XLNM1*.5+XC1*S0M)/XAM1
  965 + S2=(XC2*XA2*XLN2-XNR-F*XA2*S0)*XA2
  966 + S2P=(XC2*XAP2*XLNP2-XNR-FP*XAP2*S0P)*XAP2
  967 + S2M=(XC2*XAM2*XLNM2-XNR-FM*XAM2*S0M)*XAM2
  968 + G1=(B20*S0-0.5*XC1*XLN1)/XA1
  969 + G1P=(B2P*S0P-0.5*XC1*XLNP1)/XAP1
  970 + G1M=(B2M*S0M-0.5*XC1*XLNM1)/XAM1
  971 + G2=((0.5*F*XLN2+2.*S0*B20*XC2)*XA2+XR2)*XA2
  972 + G2P=((0.5*FP*XLNP2+2.*S0P*B2P*XC2)*XAP2+XR2)*XAP2
  973 + G2M=((0.5*FM*XLNM2+2.*S0M*B2M*XC2)*XAM2+XR2)*XAM2
  974 + BX=B0*(ZS*S0-0.5*(ZP*S0P+ZM*S0M))
  975 + * +B1*(ZS*S1-0.5*(ZP*S1P+ZM*S1M))+B2*(ZS*S2-0.5*(ZP*S2P+ZM*S2M))
  976 + BZ=B0*ALN+B1*(G1-0.5*(G1P+G1M))+B2*(G2-0.5*(G2P+G2M))
  977 +C
  978 +C CALCULATION OF THE MAGNETOTAIL CURRENT CONTRIBUTION IS FINISHED
  979 +C
  980 + RETURN
  981 + END
  982 +
  983 +C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  984 +C
  985 +C THIS CODE RETURNS THE SHIELDING FIELD REPRESENTED BY 2x3x3=18 "CARTESIAN"
  986 +C HARMONICS
  987 +C
  988 + SUBROUTINE SHLCAR3X3(A,X,Y,Z,SPS,HX,HY,HZ)
  989 +C
  990 +C - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  991 +C The 36 coefficients enter in pairs in the amplitudes of the "cartesian"
  992 +c harmonics (A(1)-A(36).
  993 +c The 12 nonlinear parameters (A(37)-A(48) are the scales Pi,Ri,Qi,and Si
  994 +C entering the arguments of exponents, sines, and cosines in each of the
  995 +C 18 "Cartesian" harmonics
  996 +C - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  997 +C
  998 + IMPLICIT REAL * 8 (A - H, O - Z)
  999 +C
  1000 + DIMENSION A(48)
  1001 +C
  1002 + CPS=DSQRT(1.D0-SPS**2)
  1003 + S3PS=4.D0*CPS**2-1.D0 ! THIS IS SIN(3*PS)/SIN(PS)
  1004 +C
  1005 + HX=0.D0
  1006 + HY=0.D0
  1007 + HZ=0.D0
  1008 + L=0
  1009 +C
  1010 + DO 1 M=1,2 ! M=1 IS FOR THE 1ST SUM ("PERP." SYMMETRY)
  1011 +C AND M=2 IS FOR THE SECOND SUM ("PARALL." SYMMETRY)
  1012 + DO 2 I=1,3
  1013 + P=A(36+I)
  1014 + Q=A(42+I)
  1015 + CYPI=DCOS(Y/P)
  1016 + CYQI=DCOS(Y/Q)
  1017 + SYPI=DSIN(Y/P)
  1018 + SYQI=DSIN(Y/Q)
  1019 +C
  1020 + DO 3 K=1,3
  1021 + R=A(39+K)
  1022 + S=A(45+K)
  1023 + SZRK=DSIN(Z/R)
  1024 + CZSK=DCOS(Z/S)
  1025 + CZRK=DCOS(Z/R)
  1026 + SZSK=DSIN(Z/S)
  1027 + SQPR=DSQRT(1.D0/P**2+1.D0/R**2)
  1028 + SQQS=DSQRT(1.D0/Q**2+1.D0/S**2)
  1029 + EPR=DEXP(X*SQPR)
  1030 + EQS=DEXP(X*SQQS)
  1031 +C
  1032 + DO 4 N=1,2 ! N=1 IS FOR THE FIRST PART OF EACH COEFFICIENT
  1033 +C AND N=2 IS FOR THE SECOND ONE
  1034 +C
  1035 + L=L+1
  1036 + IF (M.EQ.1) THEN
  1037 + IF (N.EQ.1) THEN
  1038 + DX=-SQPR*EPR*CYPI*SZRK
  1039 + DY=EPR/P*SYPI*SZRK
  1040 + DZ=-EPR/R*CYPI*CZRK
  1041 + HX=HX+A(L)*DX
  1042 + HY=HY+A(L)*DY
  1043 + HZ=HZ+A(L)*DZ
  1044 + ELSE
  1045 + DX=DX*CPS
  1046 + DY=DY*CPS
  1047 + DZ=DZ*CPS
  1048 + HX=HX+A(L)*DX
  1049 + HY=HY+A(L)*DY
  1050 + HZ=HZ+A(L)*DZ
  1051 + ENDIF
  1052 + ELSE
  1053 + IF (N.EQ.1) THEN
  1054 + DX=-SPS*SQQS*EQS*CYQI*CZSK
  1055 + DY=SPS*EQS/Q*SYQI*CZSK
  1056 + DZ=SPS*EQS/S*CYQI*SZSK
  1057 + HX=HX+A(L)*DX
  1058 + HY=HY+A(L)*DY
  1059 + HZ=HZ+A(L)*DZ
  1060 + ELSE
  1061 + DX=DX*S3PS
  1062 + DY=DY*S3PS
  1063 + DZ=DZ*S3PS
  1064 + HX=HX+A(L)*DX
  1065 + HY=HY+A(L)*DY
  1066 + HZ=HZ+A(L)*DZ
  1067 + ENDIF
  1068 + ENDIF
  1069 +c
  1070 + 4 CONTINUE
  1071 + 3 CONTINUE
  1072 + 2 CONTINUE
  1073 + 1 CONTINUE
  1074 +C
  1075 + RETURN
  1076 + END
  1077 +
  1078 +C
  1079 +C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
  1080 +C
  1081 + SUBROUTINE BIRK1TOT_02(PS,X,Y,Z,BX,BY,BZ)
  1082 +C
  1083 +C THIS IS THE SECOND VERSION OF THE ANALYTICAL MODEL OF THE REGION 1 FIELD
  1084 +C BASED ON A SEPARATE REPRESENTATION OF THE POTENTIAL FIELD IN THE INNER AND
  1085 +C OUTER SPACE, MAPPED BY MEANS OF A SPHERO-DIPOLAR COORDINATE SYSTEM (NB #3,
  1086 +C P.91). THE DIFFERENCE FROM THE FIRST ONE IS THAT INSTEAD OF OCTAGONAL
  1087 +C CURRENT LOOPS, CIRCULAR ONES ARE USED IN THIS VERSION FOR APPROXIMATING THE
  1088 +C FIELD IN THE OUTER REGION, WHICH IS FASTER.
  1089 +C
  1090 + IMPLICIT REAL*8 (A-H,O-Z)
  1091 +C
  1092 + DIMENSION D1(3,26),D2(3,79),XI(4),C1(26),C2(79)
  1093 +
  1094 + COMMON /COORD11/ XX1(12),YY1(12)
  1095 + COMMON /RHDR/ RH,DR
  1096 + COMMON /LOOPDIP1/ TILT,XCENTRE(2),RADIUS(2), DIPX,DIPY
  1097 +C
  1098 + COMMON /COORD21/ XX2(14),YY2(14),ZZ2(14)
  1099 + COMMON /DX1/ DX,SCALEIN,SCALEOUT
  1100 +C
  1101 + DATA C1/-0.911582E-03,-0.376654E-02,-0.727423E-02,-0.270084E-02,
  1102 + * -0.123899E-02,-0.154387E-02,-0.340040E-02,-0.191858E-01,
  1103 + * -0.518979E-01,0.635061E-01,0.440680,-0.396570,0.561238E-02,
  1104 + * 0.160938E-02,-0.451229E-02,-0.251810E-02,-0.151599E-02,
  1105 + * -0.133665E-02,-0.962089E-03,-0.272085E-01,-0.524319E-01,
  1106 + * 0.717024E-01,0.523439,-0.405015,-89.5587,23.2806/
  1107 +
  1108 +C
  1109 + DATA C2/6.04133,.305415,.606066E-02,.128379E-03,-.179406E-04,
  1110 + * 1.41714,-27.2586,-4.28833,-1.30675,35.5607,8.95792,.961617E-03,
  1111 + * -.801477E-03,-.782795E-03,-1.65242,-16.5242,-5.33798,.424878E-03,
  1112 + * .331787E-03,-.704305E-03,.844342E-03,.953682E-04,.886271E-03,
  1113 + * 25.1120,20.9299,5.14569,-44.1670,-51.0672,-1.87725,20.2998,
  1114 + * 48.7505,-2.97415,3.35184,-54.2921,-.838712,-10.5123,70.7594,
  1115 + * -4.94104,.106166E-03,.465791E-03,-.193719E-03,10.8439,-29.7968,
  1116 + * 8.08068,.463507E-03,-.224475E-04,.177035E-03,-.317581E-03,
  1117 + * -.264487E-03,.102075E-03,7.71390,10.1915,-4.99797,-23.1114,
  1118 + *-29.2043,12.2928,10.9542,33.6671,-9.3851,.174615E-03,-.789777E-06,
  1119 + * .686047E-03,.460104E-04,-.345216E-02,.221871E-02,.110078E-01,
  1120 + * -.661373E-02,.249201E-02,.343978E-01,-.193145E-05,.493963E-05,
  1121 + * -.535748E-04,.191833E-04,-.100496E-03,-.210103E-03,-.232195E-02,
  1122 + * .315335E-02,-.134320E-01,-.263222E-01/
  1123 +c
  1124 + DATA TILT,XCENTRE,RADIUS,DIPX,DIPY /1.00891,2.28397,-5.60831,
  1125 + * 1.86106,7.83281,1.12541,0.945719/
  1126 +
  1127 + DATA DX,SCALEIN,SCALEOUT /-0.16D0,0.08D0,0.4D0/
  1128 + DATA XX1/-11.D0,2*-7.D0,2*-3.D0,3*1.D0,2*5.D0,2*9.D0/
  1129 + DATA YY1/2.D0,0.D0,4.D0,2.D0,6.D0,0.D0,4.D0,8.D0,2.D0,6.D0,0.D0,
  1130 + * 4.D0/
  1131 + DATA XX2/-10.D0,-7.D0,2*-4.D0,0.D0,2*4.D0,7.D0,10.D0,5*0.D0/
  1132 + DATA YY2/3.D0,6.D0,3.D0,9.D0,6.D0,3.D0,9.D0,6.D0,3.D0,5*0.D0/
  1133 + DATA ZZ2/2*20.D0,4.D0,20.D0,2*4.D0,3*20.D0,2.D0,3.D0,4.5D0,
  1134 + * 7.D0,10.D0/
  1135 +C
  1136 + DATA RH,DR /9.D0,4.D0/ ! RH IS THE "HINGING DISTANCE" AND DR IS THE
  1137 +C TRANSITION SCALE LENGTH, DEFINING THE
  1138 +C CURVATURE OF THE WARPING (SEE P.89, NB #2)
  1139 +C
  1140 + DATA XLTDAY,XLTNGHT /78.D0,70.D0/ ! THESE ARE LATITUDES OF THE R-1 OVAL
  1141 +C AT NOON AND AT MIDNIGHT
  1142 + DATA DTET0 /0.034906/ ! THIS IS THE LATITUDINAL HALF-THICKNESS OF THE
  1143 +C R-1 OVAL (THE INTERPOLATION REGION BETWEEN
  1144 +C THE HIGH-LAT. AND THE PLASMA SHEET)
  1145 +C
  1146 + TNOONN=(90.D0-XLTDAY)*0.01745329D0
  1147 + TNOONS=3.141592654D0-TNOONN ! HERE WE ASSUME THAT THE POSITIONS OF
  1148 +C THE NORTHERN AND SOUTHERN R-1 OVALS
  1149 +C ARE SYMMETRIC IN THE SM-COORDINATES
  1150 + DTETDN=(XLTDAY-XLTNGHT)*0.01745329D0
  1151 + DR2=DR**2
  1152 +C
  1153 + SPS=DSIN(PS)
  1154 + R2=X**2+Y**2+Z**2
  1155 + R=DSQRT(R2)
  1156 + R3=R*R2
  1157 +C
  1158 + RMRH=R-RH
  1159 + RPRH=R+RH
  1160 + SQM=DSQRT(RMRH**2+DR2)
  1161 + SQP=DSQRT(RPRH**2+DR2)
  1162 + C=SQP-SQM
  1163 + Q=DSQRT((RH+1.D0)**2+DR2)-DSQRT((RH-1.D0)**2+DR2)
  1164 + SPSAS=SPS/R*C/Q
  1165 + CPSAS=DSQRT(1.D0-SPSAS**2)
  1166 + XAS = X*CPSAS-Z*SPSAS
  1167 + ZAS = X*SPSAS+Z*CPSAS
  1168 + IF (XAS.NE.0.D0.OR.Y.NE.0.D0) THEN
  1169 + PAS = DATAN2(Y,XAS)
  1170 + ELSE
  1171 + PAS=0.D0
  1172 + ENDIF
  1173 +C
  1174 + TAS=DATAN2(DSQRT(XAS**2+Y**2),ZAS)
  1175 + STAS=DSIN(TAS)
  1176 + F=STAS/(STAS**6*(1.D0-R3)+R3)**0.1666666667D0
  1177 +C
  1178 + TET0=DASIN(F)
  1179 + IF (TAS.GT.1.5707963D0) TET0=3.141592654D0-TET0
  1180 + DTET=DTETDN*DSIN(PAS*0.5D0)**2
  1181 + TETR1N=TNOONN+DTET
  1182 + TETR1S=TNOONS-DTET
  1183 +C
  1184 +C NOW LET'S DEFINE WHICH OF THE FOUR REGIONS (HIGH-LAT., NORTHERN PSBL,
  1185 +C PLASMA SHEET, SOUTHERN PSBL) DOES THE POINT (X,Y,Z) BELONG TO:
  1186 +C
  1187 + IF (TET0.LT.TETR1N-DTET0.OR.TET0.GT.TETR1S+DTET0) LOC=1 ! HIGH-LAT.
  1188 + IF (TET0.GT.TETR1N+DTET0.AND.TET0.LT.TETR1S-DTET0) LOC=2 ! PL.SHEET
  1189 + IF (TET0.GE.TETR1N-DTET0.AND.TET0.LE.TETR1N+DTET0) LOC=3 ! NORTH PSBL
  1190 + IF (TET0.GE.TETR1S-DTET0.AND.TET0.LE.TETR1S+DTET0) LOC=4 ! SOUTH PSBL
  1191 +C
  1192 + IF (LOC.EQ.1) THEN ! IN THE HIGH-LAT. REGION USE THE SUBROUTINE DIPOCT
  1193 +C
  1194 +C print *, ' LOC=1 (HIGH-LAT)' ! (test printout; disabled now)
  1195 + XI(1)=X
  1196 + XI(2)=Y
  1197 + XI(3)=Z
  1198 + XI(4)=PS
  1199 + CALL DIPLOOP1(XI,D1)
  1200 + BX=0.D0
  1201 + BY=0.D0
  1202 + BZ=0.D0
  1203 + DO 1 I=1,26
  1204 + BX=BX+C1(I)*D1(1,I)
  1205 + BY=BY+C1(I)*D1(2,I)
  1206 + 1 BZ=BZ+C1(I)*D1(3,I)
  1207 + ENDIF ! END OF THE CASE 1
  1208 +C
  1209 + IF (LOC.EQ.2) THEN
  1210 +C print *, ' LOC=2 (PLASMA SHEET)' ! (test printout; disabled now)
  1211 +C
  1212 + XI(1)=X
  1213 + XI(2)=Y
  1214 + XI(3)=Z
  1215 + XI(4)=PS
  1216 + CALL CONDIP1(XI,D2)
  1217 + BX=0.D0
  1218 + BY=0.D0
  1219 + BZ=0.D0
  1220 + DO 2 I=1,79
  1221 + BX=BX+C2(I)*D2(1,I)
  1222 + BY=BY+C2(I)*D2(2,I)
  1223 + 2 BZ=BZ+C2(I)*D2(3,I)
  1224 + ENDIF ! END OF THE CASE 2
  1225 +C
  1226 + IF (LOC.EQ.3) THEN
  1227 +C print *, ' LOC=3 (north PSBL)' ! (test printout; disabled now)
  1228 +C
  1229 + T01=TETR1N-DTET0
  1230 + T02=TETR1N+DTET0
  1231 + SQR=DSQRT(R)
  1232 + ST01AS=SQR/(R3+1.D0/DSIN(T01)**6-1.D0)**0.1666666667
  1233 + ST02AS=SQR/(R3+1.D0/DSIN(T02)**6-1.D0)**0.1666666667
  1234 + CT01AS=DSQRT(1.D0-ST01AS**2)
  1235 + CT02AS=DSQRT(1.D0-ST02AS**2)
  1236 + XAS1=R*ST01AS*DCOS(PAS)
  1237 + Y1= R*ST01AS*DSIN(PAS)
  1238 + ZAS1=R*CT01AS
  1239 + X1=XAS1*CPSAS+ZAS1*SPSAS
  1240 + Z1=-XAS1*SPSAS+ZAS1*CPSAS ! X1,Y1,Z1 ARE COORDS OF THE NORTHERN
  1241 +c BOUNDARY POINT
  1242 + XI(1)=X1
  1243 + XI(2)=Y1
  1244 + XI(3)=Z1
  1245 + XI(4)=PS
  1246 + CALL DIPLOOP1(XI,D1)
  1247 + BX1=0.D0
  1248 + BY1=0.D0
  1249 + BZ1=0.D0
  1250 + DO 11 I=1,26
  1251 + BX1=BX1+C1(I)*D1(1,I) ! BX1,BY1,BZ1 ARE FIELD COMPONENTS
  1252 + BY1=BY1+C1(I)*D1(2,I) ! IN THE NORTHERN BOUNDARY POINT
  1253 + 11 BZ1=BZ1+C1(I)*D1(3,I) !
  1254 +C
  1255 + XAS2=R*ST02AS*DCOS(PAS)
  1256 + Y2= R*ST02AS*DSIN(PAS)
  1257 + ZAS2=R*CT02AS
  1258 + X2=XAS2*CPSAS+ZAS2*SPSAS
  1259 + Z2=-XAS2*SPSAS+ZAS2*CPSAS ! X2,Y2,Z2 ARE COORDS OF THE SOUTHERN
  1260 +C BOUNDARY POINT
  1261 + XI(1)=X2
  1262 + XI(2)=Y2
  1263 + XI(3)=Z2
  1264 + XI(4)=PS
  1265 + CALL CONDIP1(XI,D2)
  1266 + BX2=0.D0
  1267 + BY2=0.D0
  1268 + BZ2=0.D0
  1269 + DO 12 I=1,79
  1270 + BX2=BX2+C2(I)*D2(1,I)! BX2,BY2,BZ2 ARE FIELD COMPONENTS
  1271 + BY2=BY2+C2(I)*D2(2,I) ! IN THE SOUTHERN BOUNDARY POINT
  1272 + 12 BZ2=BZ2+C2(I)*D2(3,I)
  1273 +C
  1274 +C NOW INTERPOLATE:
  1275 +C
  1276 + SS=DSQRT((X2-X1)**2+(Y2-Y1)**2+(Z2-Z1)**2)
  1277 + DS=DSQRT((X-X1)**2+(Y-Y1)**2+(Z-Z1)**2)
  1278 + FRAC=DS/SS
  1279 + BX=BX1*(1.D0-FRAC)+BX2*FRAC
  1280 + BY=BY1*(1.D0-FRAC)+BY2*FRAC
  1281 + BZ=BZ1*(1.D0-FRAC)+BZ2*FRAC
  1282 +C
  1283 + ENDIF ! END OF THE CASE 3
  1284 +C
  1285 + IF (LOC.EQ.4) THEN
  1286 +C print *, ' LOC=4 (south PSBL)' ! (test printout; disabled now)
  1287 +C
  1288 + T01=TETR1S-DTET0
  1289 + T02=TETR1S+DTET0
  1290 + SQR=DSQRT(R)
  1291 + ST01AS=SQR/(R3+1.D0/DSIN(T01)**6-1.D0)**0.1666666667
  1292 + ST02AS=SQR/(R3+1.D0/DSIN(T02)**6-1.D0)**0.1666666667
  1293 + CT01AS=-DSQRT(1.D0-ST01AS**2)
  1294 + CT02AS=-DSQRT(1.D0-ST02AS**2)
  1295 + XAS1=R*ST01AS*DCOS(PAS)
  1296 + Y1= R*ST01AS*DSIN(PAS)
  1297 + ZAS1=R*CT01AS
  1298 + X1=XAS1*CPSAS+ZAS1*SPSAS
  1299 + Z1=-XAS1*SPSAS+ZAS1*CPSAS ! X1,Y1,Z1 ARE COORDS OF THE NORTHERN
  1300 +C BOUNDARY POINT
  1301 + XI(1)=X1
  1302 + XI(2)=Y1
  1303 + XI(3)=Z1
  1304 + XI(4)=PS
  1305 + CALL CONDIP1(XI,D2)
  1306 + BX1=0.D0
  1307 + BY1=0.D0
  1308 + BZ1=0.D0
  1309 + DO 21 I=1,79
  1310 + BX1=BX1+C2(I)*D2(1,I) ! BX1,BY1,BZ1 ARE FIELD COMPONENTS
  1311 + BY1=BY1+C2(I)*D2(2,I) ! IN THE NORTHERN BOUNDARY POINT
  1312 + 21 BZ1=BZ1+C2(I)*D2(3,I) !
  1313 +C
  1314 + XAS2=R*ST02AS*DCOS(PAS)
  1315 + Y2= R*ST02AS*DSIN(PAS)
  1316 + ZAS2=R*CT02AS
  1317 + X2=XAS2*CPSAS+ZAS2*SPSAS
  1318 + Z2=-XAS2*SPSAS+ZAS2*CPSAS ! X2,Y2,Z2 ARE COORDS OF THE SOUTHERN
  1319 +C BOUNDARY POINT
  1320 + XI(1)=X2
  1321 + XI(2)=Y2
  1322 + XI(3)=Z2
  1323 + XI(4)=PS
  1324 + CALL DIPLOOP1(XI,D1)
  1325 + BX2=0.D0
  1326 + BY2=0.D0
  1327 + BZ2=0.D0
  1328 + DO 22 I=1,26
  1329 + BX2=BX2+C1(I)*D1(1,I) ! BX2,BY2,BZ2 ARE FIELD COMPONENTS
  1330 + BY2=BY2+C1(I)*D1(2,I) ! IN THE SOUTHERN BOUNDARY POINT
  1331 + 22 BZ2=BZ2+C1(I)*D1(3,I)
  1332 +C
  1333 +C NOW INTERPOLATE:
  1334 +C
  1335 + SS=DSQRT((X2-X1)**2+(Y2-Y1)**2+(Z2-Z1)**2)
  1336 + DS=DSQRT((X-X1)**2+(Y-Y1)**2+(Z-Z1)**2)
  1337 + FRAC=DS/SS
  1338 + BX=BX1*(1.D0-FRAC)+BX2*FRAC
  1339 + BY=BY1*(1.D0-FRAC)+BY2*FRAC
  1340 + BZ=BZ1*(1.D0-FRAC)+BZ2*FRAC
  1341 +C
  1342 + ENDIF ! END OF THE CASE 4
  1343 +C
  1344 +C NOW, LET US ADD THE SHIELDING FIELD
  1345 +C
  1346 + CALL BIRK1SHLD(PS,X,Y,Z,BSX,BSY,BSZ)
  1347 + BX=BX+BSX
  1348 + BY=BY+BSY
  1349 + BZ=BZ+BSZ
  1350 + RETURN
  1351 + END
  1352 +C
  1353 +C------------------------------------------------------------------------------
  1354 +C
  1355 +C
  1356 + SUBROUTINE DIPLOOP1(XI,D)
  1357 +C
  1358 +C
  1359 +C Calculates dependent model variables and their deriva-
  1360 +C tives for given independent variables and model parame-
  1361 +C ters. Specifies model functions with free parameters which
  1362 +C must be determined by means of least squares fits (RMS
  1363 +C minimization procedure).
  1364 +C
  1365 +C Description of parameters:
  1366 +C
  1367 +C XI - input vector containing independent variables;
  1368 +C D - output double precision vector containing
  1369 +C calculated values for derivatives of dependent
  1370 +C variables with respect to LINEAR model parameters;
  1371 +C
  1372 +C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  1373 +C
  1374 +c The 26 coefficients are moments (Z- and X-components) of 12 dipoles placed
  1375 +C inside the R1-shell, PLUS amplitudes of two octagonal double loops.
  1376 +C The dipoles with nonzero Yi appear in pairs with equal moments.
  1377 +c (see the notebook #2, pp.102-103, for details)
  1378 +C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  1379 +c
  1380 + IMPLICIT REAL * 8 (A - H, O - Z)
  1381 +C
  1382 + COMMON /COORD11/ XX(12),YY(12)
  1383 + COMMON /LOOPDIP1/ TILT,XCENTRE(2),RADIUS(2), DIPX,DIPY
  1384 + COMMON /RHDR/RH,DR
  1385 + DIMENSION XI(4),D(3,26)
  1386 +C
  1387 + X = XI(1)
  1388 + Y = XI(2)
  1389 + Z = XI(3)
  1390 + PS= XI(4)
  1391 + SPS=DSIN(PS)
  1392 +C
  1393 + DO 1 I=1,12
  1394 + R2=(XX(I)*DIPX)**2+(YY(I)*DIPY)**2
  1395 + R=DSQRT(R2)
  1396 + RMRH=R-RH
  1397 + RPRH=R+RH
  1398 + DR2=DR**2
  1399 + SQM=DSQRT(RMRH**2+DR2)
  1400 + SQP=DSQRT(RPRH**2+DR2)
  1401 + C=SQP-SQM
  1402 + Q=DSQRT((RH+1.D0)**2+DR2)-DSQRT((RH-1.D0)**2+DR2)
  1403 + SPSAS=SPS/R*C/Q
  1404 + CPSAS=DSQRT(1.D0-SPSAS**2)
  1405 + XD= (XX(I)*DIPX)*CPSAS
  1406 + YD= (YY(I)*DIPY)
  1407 + ZD=-(XX(I)*DIPX)*SPSAS
  1408 + CALL DIPXYZ(X-XD,Y-YD,Z-ZD,BX1X,BY1X,BZ1X,BX1Y,BY1Y,BZ1Y,
  1409 + * BX1Z,BY1Z,BZ1Z)
  1410 + IF (DABS(YD).GT.1.D-10) THEN
  1411 + CALL DIPXYZ(X-XD,Y+YD,Z-ZD,BX2X,BY2X,BZ2X,BX2Y,BY2Y,BZ2Y,
  1412 + * BX2Z,BY2Z,BZ2Z)
  1413 + ELSE
  1414 + BX2X=0.D0
  1415 + BY2X=0.D0
  1416 + BZ2X=0.D0
  1417 +C
  1418 + BX2Z=0.D0
  1419 + BY2Z=0.D0
  1420 + BZ2Z=0.D0
  1421 + ENDIF
  1422 +C
  1423 + D(1,I)=BX1Z+BX2Z
  1424 + D(2,I)=BY1Z+BY2Z
  1425 + D(3,I)=BZ1Z+BZ2Z
  1426 + D(1,I+12)=(BX1X+BX2X)*SPS
  1427 + D(2,I+12)=(BY1X+BY2X)*SPS
  1428 + D(3,I+12)=(BZ1X+BZ2X)*SPS
  1429 + 1 CONTINUE
  1430 +c
  1431 + R2=(XCENTRE(1)+RADIUS(1))**2
  1432 + R=DSQRT(R2)
  1433 + RMRH=R-RH
  1434 + RPRH=R+RH
  1435 + DR2=DR**2
  1436 + SQM=DSQRT(RMRH**2+DR2)
  1437 + SQP=DSQRT(RPRH**2+DR2)
  1438 + C=SQP-SQM
  1439 + Q=DSQRT((RH+1.D0)**2+DR2)-DSQRT((RH-1.D0)**2+DR2)
  1440 + SPSAS=SPS/R*C/Q
  1441 + CPSAS=DSQRT(1.D0-SPSAS**2)
  1442 + XOCT1= X*CPSAS-Z*SPSAS
  1443 + YOCT1= Y
  1444 + ZOCT1= X*SPSAS+Z*CPSAS
  1445 +C
  1446 + CALL CROSSLP(XOCT1,YOCT1,ZOCT1,BXOCT1,BYOCT1,BZOCT1,XCENTRE(1),
  1447 + * RADIUS(1),TILT)
  1448 + D(1,25)=BXOCT1*CPSAS+BZOCT1*SPSAS
  1449 + D(2,25)=BYOCT1
  1450 + D(3,25)=-BXOCT1*SPSAS+BZOCT1*CPSAS
  1451 +C
  1452 + R2=(RADIUS(2)-XCENTRE(2))**2
  1453 + R=DSQRT(R2)
  1454 + RMRH=R-RH
  1455 + RPRH=R+RH
  1456 + DR2=DR**2
  1457 + SQM=DSQRT(RMRH**2+DR2)
  1458 + SQP=DSQRT(RPRH**2+DR2)
  1459 + C=SQP-SQM
  1460 + Q=DSQRT((RH+1.D0)**2+DR2)-DSQRT((RH-1.D0)**2+DR2)
  1461 + SPSAS=SPS/R*C/Q
  1462 + CPSAS=DSQRT(1.D0-SPSAS**2)
  1463 + XOCT2= X*CPSAS-Z*SPSAS -XCENTRE(2)
  1464 + YOCT2= Y
  1465 + ZOCT2= X*SPSAS+Z*CPSAS
  1466 + CALL CIRCLE(XOCT2,YOCT2,ZOCT2,RADIUS(2),BX,BY,BZ)
  1467 + D(1,26) = BX*CPSAS+BZ*SPSAS
  1468 + D(2,26) = BY
  1469 + D(3,26) = -BX*SPSAS+BZ*CPSAS
  1470 +C
  1471 + RETURN
  1472 + END
  1473 +c-------------------------------------------------------------------------
  1474 +C
  1475 + SUBROUTINE CIRCLE(X,Y,Z,RL,BX,BY,BZ)
  1476 +C
  1477 +C RETURNS COMPONENTS OF THE FIELD FROM A CIRCULAR CURRENT LOOP OF RADIUS RL
  1478 +C USES THE SECOND (MORE ACCURATE) APPROXIMATION GIVEN IN ABRAMOWITZ AND STEGUN
  1479 +
  1480 + IMPLICIT REAL*8 (A-H,O-Z)
  1481 + REAL*8 K
  1482 + DATA PI/3.141592654D0/
  1483 +C
  1484 + RHO2=X*X+Y*Y
  1485 + RHO=DSQRT(RHO2)
  1486 + R22=Z*Z+(RHO+RL)**2
  1487 + R2=DSQRT(R22)
  1488 + R12=R22-4.D0*RHO*RL
  1489 + R32=0.5D0*(R12+R22)
  1490 + XK2=1.D0-R12/R22
  1491 + XK2S=1.D0-XK2
  1492 + DL=DLOG(1.D0/XK2S)
  1493 + K=1.38629436112d0+XK2S*(0.09666344259D0+XK2S*(0.03590092383+
  1494 + * XK2S*(0.03742563713+XK2S*0.01451196212))) +DL*
  1495 + * (0.5D0+XK2S*(0.12498593597D0+XK2S*(0.06880248576D0+
  1496 + * XK2S*(0.03328355346D0+XK2S*0.00441787012D0))))
  1497 + E=1.D0+XK2S*(0.44325141463D0+XK2S*(0.0626060122D0+XK2S*
  1498 + * (0.04757383546D0+XK2S*0.01736506451D0))) +DL*
  1499 + * XK2S*(0.2499836831D0+XK2S*(0.09200180037D0+XK2S*
  1500 + * (0.04069697526D0+XK2S*0.00526449639D0)))
  1501 +
  1502 + IF (RHO.GT.1.D-6) THEN
  1503 + BRHO=Z/(RHO2*R2)*(R32/R12*E-K) ! THIS IS NOT EXACTLY THE B-RHO COM-
  1504 + ELSE ! PONENT - NOTE THE ADDITIONAL
  1505 + BRHO=PI*RL/R2*(RL-RHO)/R12*Z/(R32-RHO2) ! DIVISION BY RHO
  1506 + ENDIF
  1507 +
  1508 + BX=BRHO*X
  1509 + BY=BRHO*Y
  1510 + BZ=(K-E*(R32-2.D0*RL*RL)/R12)/R2
  1511 + RETURN
  1512 + END
  1513 +C-------------------------------------------------------------
  1514 +C
  1515 + SUBROUTINE CROSSLP(X,Y,Z,BX,BY,BZ,XC,RL,AL)
  1516 +C
  1517 +c RETURNS FIELD COMPONENTS OF A PAIR OF LOOPS WITH A COMMON CENTER AND
  1518 +C DIAMETER, COINCIDING WITH THE X AXIS. THE LOOPS ARE INCLINED TO THE
  1519 +C EQUATORIAL PLANE BY THE ANGLE AL (RADIANS) AND SHIFTED IN THE POSITIVE
  1520 +C X-DIRECTION BY THE DISTANCE XC.
  1521 +c
  1522 + IMPLICIT REAL*8 (A-H,O-Z)
  1523 +C
  1524 + CAL=DCOS(AL)
  1525 + SAL=DSIN(AL)
  1526 +C
  1527 + Y1=Y*CAL-Z*SAL
  1528 + Z1=Y*SAL+Z*CAL
  1529 + Y2=Y*CAL+Z*SAL
  1530 + Z2=-Y*SAL+Z*CAL
  1531 + CALL CIRCLE(X-XC,Y1,Z1,RL,BX1,BY1,BZ1)
  1532 + CALL CIRCLE(X-XC,Y2,Z2,RL,BX2,BY2,BZ2)
  1533 + BX=BX1+BX2
  1534 + BY= (BY1+BY2)*CAL+(BZ1-BZ2)*SAL
  1535 + BZ=-(BY1-BY2)*SAL+(BZ1+BZ2)*CAL
  1536 +C
  1537 + RETURN
  1538 + END
  1539 +
  1540 +C*******************************************************************
  1541 +
  1542 + SUBROUTINE DIPXYZ(X,Y,Z,BXX,BYX,BZX,BXY,BYY,BZY,BXZ,BYZ,BZZ)
  1543 +C
  1544 +C RETURNS THE FIELD COMPONENTS PRODUCED BY THREE DIPOLES, EACH
  1545 +C HAVING M=Me AND ORIENTED PARALLEL TO X,Y, and Z AXIS, RESP.
  1546 +C
  1547 + IMPLICIT REAL*8 (A-H,O-Z)
  1548 +C
  1549 + X2=X**2
  1550 + Y2=Y**2
  1551 + Z2=Z**2
  1552 + R2=X2+Y2+Z2
  1553 +
  1554 + XMR5=30574.D0/(R2*R2*DSQRT(R2))
  1555 + XMR53=3.D0*XMR5
  1556 + BXX=XMR5*(3.D0*X2-R2)
  1557 + BYX=XMR53*X*Y
  1558 + BZX=XMR53*X*Z
  1559 +C
  1560 + BXY=BYX
  1561 + BYY=XMR5*(3.D0*Y2-R2)
  1562 + BZY=XMR53*Y*Z
  1563 +C
  1564 + BXZ=BZX
  1565 + BYZ=BZY
  1566 + BZZ=XMR5*(3.D0*Z2-R2)
  1567 +C
  1568 + RETURN
  1569 + END
  1570 +C
  1571 +C------------------------------------------------------------------------------
  1572 + SUBROUTINE CONDIP1(XI,D)
  1573 +C
  1574 +C Calculates dependent model variables and their derivatives for given
  1575 +C independent variables and model parameters. Specifies model functions with
  1576 +C free parameters which must be determined by means of least squares fits
  1577 +C (RMS minimization procedure).
  1578 +C
  1579 +C Description of parameters:
  1580 +C
  1581 +C XI - input vector containing independent variables;
  1582 +C D - output double precision vector containing
  1583 +C calculated values for derivatives of dependent
  1584 +C variables with respect to LINEAR model parameters;
  1585 +C
  1586 +C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  1587 +C
  1588 +c The 79 coefficients are (1) 5 amplitudes of the conical harmonics, plus
  1589 +c (2) (9x3+5x2)x2=74 components of the dipole moments
  1590 +c (see the notebook #2, pp.113-..., for details)
  1591 +C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  1592 +c
  1593 + IMPLICIT REAL * 8 (A - H, O - Z)
  1594 +C
  1595 + COMMON /DX1/ DX,SCALEIN,SCALEOUT
  1596 + COMMON /COORD21/ XX(14),YY(14),ZZ(14)
  1597 +c
  1598 + DIMENSION XI(4),D(3,79),CF(5),SF(5)
  1599 +C
  1600 + X = XI(1)
  1601 + Y = XI(2)
  1602 + Z = XI(3)
  1603 + PS= XI(4)
  1604 + SPS=DSIN(PS)
  1605 + CPS=DCOS(PS)
  1606 +C
  1607 + XSM=X*CPS-Z*SPS - DX
  1608 + ZSM=Z*CPS+X*SPS
  1609 + RO2=XSM**2+Y**2
  1610 + RO=SQRT(RO2)
  1611 +C
  1612 + CF(1)=XSM/RO
  1613 + SF(1)=Y/RO
  1614 +C
  1615 + CF(2)=CF(1)**2-SF(1)**2
  1616 + SF(2)=2.*SF(1)*CF(1)
  1617 + CF(3)=CF(2)*CF(1)-SF(2)*SF(1)
  1618 + SF(3)=SF(2)*CF(1)+CF(2)*SF(1)
  1619 + CF(4)=CF(3)*CF(1)-SF(3)*SF(1)
  1620 + SF(4)=SF(3)*CF(1)+CF(3)*SF(1)
  1621 + CF(5)=CF(4)*CF(1)-SF(4)*SF(1)
  1622 + SF(5)=SF(4)*CF(1)+CF(4)*SF(1)
  1623 +C
  1624 + R2=RO2+ZSM**2
  1625 + R=DSQRT(R2)
  1626 + C=ZSM/R
  1627 + S=RO/R
  1628 + CH=DSQRT(0.5D0*(1.D0+C))
  1629 + SH=DSQRT(0.5D0*(1.D0-C))
  1630 + TNH=SH/CH
  1631 + CNH=1.D0/TNH
  1632 +C
  1633 + DO 1 M=1,5
  1634 + BT=M*CF(M)/(R*S)*(TNH**M+CNH**M)
  1635 + BF=-0.5D0*M*SF(M)/R*(TNH**(M-1)/CH**2-CNH**(M-1)/SH**2)
  1636 + BXSM=BT*C*CF(1)-BF*SF(1)
  1637 + BY=BT*C*SF(1)+BF*CF(1)
  1638 + BZSM=-BT*S
  1639 +C
  1640 + D(1,M)=BXSM*CPS+BZSM*SPS
  1641 + D(2,M)=BY
  1642 + 1 D(3,M)=-BXSM*SPS+BZSM*CPS
  1643 +C
  1644 + XSM = X*CPS-Z*SPS
  1645 + ZSM = Z*CPS+X*SPS
  1646 +C
  1647 + DO 2 I=1,9
  1648 +C
  1649 + IF (I.EQ.3.OR.I.EQ.5.OR.I.EQ.6) THEN
  1650 + XD = XX(I)*SCALEIN
  1651 + YD = YY(I)*SCALEIN
  1652 + ELSE
  1653 + XD = XX(I)*SCALEOUT
  1654 + YD = YY(I)*SCALEOUT
  1655 + ENDIF
  1656 +C
  1657 + ZD = ZZ(I)
  1658 +C
  1659 + CALL DIPXYZ(XSM-XD,Y-YD,ZSM-ZD,BX1X,BY1X,BZ1X,BX1Y,BY1Y,BZ1Y,
  1660 + * BX1Z,BY1Z,BZ1Z)
  1661 + CALL DIPXYZ(XSM-XD,Y+YD,ZSM-ZD,BX2X,BY2X,BZ2X,BX2Y,BY2Y,BZ2Y,
  1662 + * BX2Z,BY2Z,BZ2Z)
  1663 + CALL DIPXYZ(XSM-XD,Y-YD,ZSM+ZD,BX3X,BY3X,BZ3X,BX3Y,BY3Y,BZ3Y,
  1664 + * BX3Z,BY3Z,BZ3Z)
  1665 + CALL DIPXYZ(XSM-XD,Y+YD,ZSM+ZD,BX4X,BY4X,BZ4X,BX4Y,BY4Y,BZ4Y,
  1666 + * BX4Z,BY4Z,BZ4Z)
  1667 +C
  1668 + IX=I*3+3
  1669 + IY=IX+1
  1670 + IZ=IY+1
  1671 +C
  1672 + D(1,IX)=(BX1X+BX2X-BX3X-BX4X)*CPS+(BZ1X+BZ2X-BZ3X-BZ4X)*SPS
  1673 + D(2,IX)= BY1X+BY2X-BY3X-BY4X
  1674 + D(3,IX)=(BZ1X+BZ2X-BZ3X-BZ4X)*CPS-(BX1X+BX2X-BX3X-BX4X)*SPS
  1675 +C
  1676 + D(1,IY)=(BX1Y-BX2Y-BX3Y+BX4Y)*CPS+(BZ1Y-BZ2Y-BZ3Y+BZ4Y)*SPS
  1677 + D(2,IY)= BY1Y-BY2Y-BY3Y+BY4Y
  1678 + D(3,IY)=(BZ1Y-BZ2Y-BZ3Y+BZ4Y)*CPS-(BX1Y-BX2Y-BX3Y+BX4Y)*SPS
  1679 +C
  1680 + D(1,IZ)=(BX1Z+BX2Z+BX3Z+BX4Z)*CPS+(BZ1Z+BZ2Z+BZ3Z+BZ4Z)*SPS
  1681 + D(2,IZ)= BY1Z+BY2Z+BY3Z+BY4Z
  1682 + D(3,IZ)=(BZ1Z+BZ2Z+BZ3Z+BZ4Z)*CPS-(BX1Z+BX2Z+BX3Z+BX4Z)*SPS
  1683 +C
  1684 + IX=IX+27
  1685 + IY=IY+27
  1686 + IZ=IZ+27
  1687 +C
  1688 + D(1,IX)=SPS*((BX1X+BX2X+BX3X+BX4X)*CPS+(BZ1X+BZ2X+BZ3X+BZ4X)*SPS)
  1689 + D(2,IX)=SPS*(BY1X+BY2X+BY3X+BY4X)
  1690 + D(3,IX)=SPS*((BZ1X+BZ2X+BZ3X+BZ4X)*CPS-(BX1X+BX2X+BX3X+BX4X)*SPS)
  1691 +C
  1692 + D(1,IY)=SPS*((BX1Y-BX2Y+BX3Y-BX4Y)*CPS+(BZ1Y-BZ2Y+BZ3Y-BZ4Y)*SPS)
  1693 + D(2,IY)=SPS*(BY1Y-BY2Y+BY3Y-BY4Y)
  1694 + D(3,IY)=SPS*((BZ1Y-BZ2Y+BZ3Y-BZ4Y)*CPS-(BX1Y-BX2Y+BX3Y-BX4Y)*SPS)
  1695 +C
  1696 + D(1,IZ)=SPS*((BX1Z+BX2Z-BX3Z-BX4Z)*CPS+(BZ1Z+BZ2Z-BZ3Z-BZ4Z)*SPS)
  1697 + D(2,IZ)=SPS*(BY1Z+BY2Z-BY3Z-BY4Z)
  1698 + D(3,IZ)=SPS*((BZ1Z+BZ2Z-BZ3Z-BZ4Z)*CPS-(BX1Z+BX2Z-BX3Z-BX4Z)*SPS)
  1699 + 2 CONTINUE
  1700 +C
  1701 + DO 3 I=1,5
  1702 + ZD=ZZ(I+9)
  1703 + CALL DIPXYZ(XSM,Y,ZSM-ZD,BX1X,BY1X,BZ1X,BX1Y,BY1Y,BZ1Y,BX1Z,BY1Z,
  1704 + * BZ1Z)
  1705 + CALL DIPXYZ(XSM,Y,ZSM+ZD,BX2X,BY2X,BZ2X,BX2Y,BY2Y,BZ2Y,BX2Z,BY2Z,
  1706 + * BZ2Z)
  1707 + IX=58+I*2
  1708 + IZ=IX+1
  1709 + D(1,IX)=(BX1X-BX2X)*CPS+(BZ1X-BZ2X)*SPS
  1710 + D(2,IX)=BY1X-BY2X
  1711 + D(3,IX)=(BZ1X-BZ2X)*CPS-(BX1X-BX2X)*SPS
  1712 +C
  1713 + D(1,IZ)=(BX1Z+BX2Z)*CPS+(BZ1Z+BZ2Z)*SPS
  1714 + D(2,IZ)=BY1Z+BY2Z
  1715 + D(3,IZ)=(BZ1Z+BZ2Z)*CPS-(BX1Z+BX2Z)*SPS
  1716 +C
  1717 + IX=IX+10
  1718 + IZ=IZ+10
  1719 + D(1,IX)=SPS*((BX1X+BX2X)*CPS+(BZ1X+BZ2X)*SPS)
  1720 + D(2,IX)=SPS*(BY1X+BY2X)
  1721 + D(3,IX)=SPS*((BZ1X+BZ2X)*CPS-(BX1X+BX2X)*SPS)
  1722 +C
  1723 + D(1,IZ)=SPS*((BX1Z-BX2Z)*CPS+(BZ1Z-BZ2Z)*SPS)
  1724 + D(2,IZ)=SPS*(BY1Z-BY2Z)
  1725 + 3 D(3,IZ)=SPS*((BZ1Z-BZ2Z)*CPS-(BX1Z-BX2Z)*SPS)
  1726 +C
  1727 + RETURN
  1728 + END
  1729 +C
  1730 +C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  1731 +C
  1732 + SUBROUTINE BIRK1SHLD(PS,X,Y,Z,BX,BY,BZ)
  1733 +C
  1734 +C - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  1735 +C
  1736 +C The 64 linear parameters are amplitudes of the "box" harmonics.
  1737 +c The 16 nonlinear parameters are the scales Pi, and Qk entering the arguments
  1738 +C of sines/cosines and exponents in each of 32 cartesian harmonics
  1739 +c N.A. Tsyganenko, Spring 1994, adjusted for the Birkeland field Aug.22, 1995
  1740 +c Revised June 12, 1996.
  1741 +C - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  1742 +C
  1743 + IMPLICIT REAL * 8 (A - H, O - Z)
  1744 +C
  1745 + DIMENSION A(80)
  1746 + DIMENSION P1(4),R1(4),Q1(4),S1(4),RP(4),RR(4),RQ(4),RS(4)
  1747 +C
  1748 + EQUIVALENCE (P1(1),A(65)),(R1(1),A(69)),(Q1(1),A(73)),
  1749 + * (S1(1),A(77))
  1750 +C
  1751 + DATA A/1.174198045,-1.463820502,4.840161537,-3.674506864,
  1752 + * 82.18368896,-94.94071588,-4122.331796,4670.278676,-21.54975037,
  1753 + * 26.72661293,-72.81365728,44.09887902,40.08073706,-51.23563510,
  1754 + * 1955.348537,-1940.971550,794.0496433,-982.2441344,1889.837171,
  1755 + * -558.9779727,-1260.543238,1260.063802,-293.5942373,344.7250789,
  1756 + * -773.7002492,957.0094135,-1824.143669,520.7994379,1192.484774,
  1757 + * -1192.184565,89.15537624,-98.52042999,-0.8168777675E-01,
  1758 + * 0.4255969908E-01,0.3155237661,-0.3841755213,2.494553332,
  1759 + * -0.6571440817E-01,-2.765661310,0.4331001908,0.1099181537,
  1760 + * -0.6154126980E-01,-0.3258649260,0.6698439193,-5.542735524,
  1761 + * 0.1604203535,5.854456934,-0.8323632049,3.732608869,-3.130002153,
  1762 + * 107.0972607,-32.28483411,-115.2389298,54.45064360,-0.5826853320,
  1763 + * -3.582482231,-4.046544561,3.311978102,-104.0839563,30.26401293,
  1764 + * 97.29109008,-50.62370872,-296.3734955,127.7872523,5.303648988,
  1765 + * 10.40368955,69.65230348,466.5099509,1.645049286,3.825838190,
  1766 + * 11.66675599,558.9781177,1.826531343,2.066018073,25.40971369,
  1767 + * 990.2795225,2.319489258,4.555148484,9.691185703,591.8280358/
  1768 +C
  1769 + BX=0.D0
  1770 + BY=0.D0
  1771 + BZ=0.D0
  1772 + CPS=DCOS(PS)
  1773 + SPS=DSIN(PS)
  1774 + S3PS=4.D0*CPS**2-1.D0
  1775 +C
  1776 + DO 11 I=1,4
  1777 + RP(I)=1.D0/P1(I)
  1778 + RR(I)=1.D0/R1(I)
  1779 + RQ(I)=1.D0/Q1(I)
  1780 + 11 RS(I)=1.D0/S1(I)
  1781 +C
  1782 + L=0
  1783 +C
  1784 + DO 1 M=1,2 ! M=1 IS FOR THE 1ST SUM ("PERP." SYMMETRY)
  1785 +C AND M=2 IS FOR THE SECOND SUM ("PARALL." SYMMETRY)
  1786 + DO 2 I=1,4
  1787 + CYPI=DCOS(Y*RP(I))
  1788 + CYQI=DCOS(Y*RQ(I))
  1789 + SYPI=DSIN(Y*RP(I))
  1790 + SYQI=DSIN(Y*RQ(I))
  1791 +C
  1792 + DO 3 K=1,4
  1793 + SZRK=DSIN(Z*RR(K))
  1794 + CZSK=DCOS(Z*RS(K))
  1795 + CZRK=DCOS(Z*RR(K))
  1796 + SZSK=DSIN(Z*RS(K))
  1797 + SQPR=DSQRT(RP(I)**2+RR(K)**2)
  1798 + SQQS=DSQRT(RQ(I)**2+RS(K)**2)
  1799 + EPR=DEXP(X*SQPR)
  1800 + EQS=DEXP(X*SQQS)
  1801 +C
  1802 + DO 4 N=1,2 ! N=1 IS FOR THE FIRST PART OF EACH COEFFICIENT
  1803 +C AND N=2 IS FOR THE SECOND ONE
  1804 + IF (M.EQ.1) THEN
  1805 + IF (N.EQ.1) THEN
  1806 + HX=-SQPR*EPR*CYPI*SZRK
  1807 + HY=RP(I)*EPR*SYPI*SZRK
  1808 + HZ=-RR(K)*EPR*CYPI*CZRK
  1809 + ELSE
  1810 + HX=HX*CPS
  1811 + HY=HY*CPS
  1812 + HZ=HZ*CPS
  1813 + ENDIF
  1814 + ELSE
  1815 + IF (N.EQ.1) THEN
  1816 + HX=-SPS*SQQS*EQS*CYQI*CZSK
  1817 + HY=SPS*RQ(I)*EQS*SYQI*CZSK
  1818 + HZ=SPS*RS(K)*EQS*CYQI*SZSK
  1819 + ELSE
  1820 + HX=HX*S3PS
  1821 + HY=HY*S3PS
  1822 + HZ=HZ*S3PS
  1823 + ENDIF
  1824 + ENDIF
  1825 + L=L+1
  1826 +c
  1827 + BX=BX+A(L)*HX
  1828 + BY=BY+A(L)*HY
  1829 + 4 BZ=BZ+A(L)*HZ
  1830 + 3 CONTINUE
  1831 + 2 CONTINUE
  1832 + 1 CONTINUE
  1833 +C
  1834 + RETURN
  1835 + END
  1836 +C
  1837 +C##########################################################################
  1838 +C
  1839 + SUBROUTINE BIRK2TOT_02(PS,X,Y,Z,BX,BY,BZ)
  1840 +C
  1841 + IMPLICIT REAL*8 (A-H,O-Z)
  1842 +C
  1843 + CALL BIRK2SHL(X,Y,Z,PS,WX,WY,WZ)
  1844 + CALL R2_BIRK(X,Y,Z,PS,HX,HY,HZ)
  1845 + BX=WX+HX
  1846 + BY=WY+HY
  1847 + BZ=WZ+HZ
  1848 + RETURN
  1849 + END
  1850 +C
  1851 +C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  1852 +C
  1853 +C THIS CODE IS FOR THE FIELD FROM 2x2x2=8 "CARTESIAN" HARMONICS
  1854 +C
  1855 + SUBROUTINE BIRK2SHL(X,Y,Z,PS,HX,HY,HZ)
  1856 +C
  1857 +C - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  1858 +C The model parameters are provided to this module via common-block /A/.
  1859 +C The 16 linear parameters enter in pairs in the amplitudes of the
  1860 +c "cartesian" harmonics.
  1861 +c The 8 nonlinear parameters are the scales Pi,Ri,Qi,and Si entering the
  1862 +c arguments of exponents, sines, and cosines in each of the 8 "Cartesian"
  1863 +c harmonics
  1864 +C - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  1865 +C
  1866 + IMPLICIT REAL * 8 (A - H, O - Z)
  1867 +C
  1868 + DIMENSION P(2),R(2),Q(2),S(2)
  1869 + DIMENSION A(24)
  1870 +C
  1871 + EQUIVALENCE(P(1),A(17)),(R(1),A(19)),(Q(1),A(21)),(S(1),A(23))
  1872 + DATA A/-111.6371348,124.5402702,110.3735178,-122.0095905,
  1873 + * 111.9448247,-129.1957743,-110.7586562,126.5649012,-0.7865034384,
  1874 + * -0.2483462721,0.8026023894,0.2531397188,10.72890902,0.8483902118,
  1875 + * -10.96884315,-0.8583297219,13.85650567,14.90554500,10.21914434,
  1876 + * 10.09021632,6.340382460,14.40432686,12.71023437,12.83966657/
  1877 +C
  1878 + CPS=DCOS(PS)
  1879 + SPS=DSIN(PS)
  1880 + S3PS=4.D0*CPS**2-1.D0 ! THIS IS SIN(3*PS)/SIN(PS)
  1881 +C
  1882 + HX=0.D0
  1883 + HY=0.D0
  1884 + HZ=0.D0
  1885 + L=0
  1886 +C
  1887 + DO 1 M=1,2 ! M=1 IS FOR THE 1ST SUM ("PERP." SYMMETRY)
  1888 +C AND M=2 IS FOR THE SECOND SUM ("PARALL." SYMMETRY)
  1889 + DO 2 I=1,2
  1890 + CYPI=DCOS(Y/P(I))
  1891 + CYQI=DCOS(Y/Q(I))
  1892 + SYPI=DSIN(Y/P(I))
  1893 + SYQI=DSIN(Y/Q(I))
  1894 +C
  1895 + DO 3 K=1,2
  1896 + SZRK=DSIN(Z/R(K))
  1897 + CZSK=DCOS(Z/S(K))
  1898 + CZRK=DCOS(Z/R(K))
  1899 + SZSK=DSIN(Z/S(K))
  1900 + SQPR=DSQRT(1.D0/P(I)**2+1.D0/R(K)**2)
  1901 + SQQS=DSQRT(1.D0/Q(I)**2+1.D0/S(K)**2)
  1902 + EPR=DEXP(X*SQPR)
  1903 + EQS=DEXP(X*SQQS)
  1904 +C
  1905 + DO 4 N=1,2 ! N=1 IS FOR THE FIRST PART OF EACH COEFFICIENT
  1906 +C AND N=2 IS FOR THE SECOND ONE
  1907 +C
  1908 + L=L+1
  1909 + IF (M.EQ.1) THEN
  1910 + IF (N.EQ.1) THEN
  1911 + DX=-SQPR*EPR*CYPI*SZRK
  1912 + DY=EPR/P(I)*SYPI*SZRK
  1913 + DZ=-EPR/R(K)*CYPI*CZRK
  1914 + HX=HX+A(L)*DX
  1915 + HY=HY+A(L)*DY
  1916 + HZ=HZ+A(L)*DZ
  1917 + ELSE
  1918 + DX=DX*CPS
  1919 + DY=DY*CPS
  1920 + DZ=DZ*CPS
  1921 + HX=HX+A(L)*DX
  1922 + HY=HY+A(L)*DY
  1923 + HZ=HZ+A(L)*DZ
  1924 + ENDIF
  1925 + ELSE
  1926 + IF (N.EQ.1) THEN
  1927 + DX=-SPS*SQQS*EQS*CYQI*CZSK
  1928 + DY=SPS*EQS/Q(I)*SYQI*CZSK
  1929 + DZ=SPS*EQS/S(K)*CYQI*SZSK
  1930 + HX=HX+A(L)*DX
  1931 + HY=HY+A(L)*DY
  1932 + HZ=HZ+A(L)*DZ
  1933 + ELSE
  1934 + DX=DX*S3PS
  1935 + DY=DY*S3PS
  1936 + DZ=DZ*S3PS
  1937 + HX=HX+A(L)*DX
  1938 + HY=HY+A(L)*DY
  1939 + HZ=HZ+A(L)*DZ
  1940 + ENDIF
  1941 + ENDIF
  1942 +c
  1943 + 4 CONTINUE
  1944 + 3 CONTINUE
  1945 + 2 CONTINUE
  1946 + 1 CONTINUE
  1947 +C
  1948 + RETURN
  1949 + END
  1950 +
  1951 +c********************************************************************
  1952 +C
  1953 + SUBROUTINE R2_BIRK(X,Y,Z,PS,BX,BY,BZ)
  1954 +C
  1955 +C RETURNS THE MODEL FIELD FOR THE REGION 2 BIRKELAND CURRENT/PARTIAL RC
  1956 +C (WITHOUT SHIELDING FIELD)
  1957 +C
  1958 + IMPLICIT REAL*8 (A-H,O-Z)
  1959 + SAVE PSI,CPS,SPS
  1960 + DATA DELARG/0.030D0/,DELARG1/0.015D0/,PSI/10.D0/
  1961 +C
  1962 + IF (DABS(PSI-PS).GT.1.D-10) THEN
  1963 + PSI=PS
  1964 + CPS=DCOS(PS)
  1965 + SPS=DSIN(PS)
  1966 + ENDIF
  1967 +C
  1968 + XSM=X*CPS-Z*SPS
  1969 + ZSM=Z*CPS+X*SPS
  1970 +C
  1971 + XKS=XKSI(XSM,Y,ZSM)
  1972 + IF (XKS.LT.-(DELARG+DELARG1)) THEN
  1973 + CALL R2OUTER(XSM,Y,ZSM,BXSM,BY,BZSM)
  1974 + BXSM=-BXSM*0.02 ! ALL COMPONENTS ARE MULTIPLIED BY THE
  1975 + BY=-BY*0.02 ! FACTOR -0.02, IN ORDER TO NORMALIZE THE
  1976 + BZSM=-BZSM*0.02 ! FIELD (SO THAT Bz=-1 nT at X=-5.3 RE, Y=Z=0)
  1977 + ENDIF
  1978 + IF (XKS.GE.-(DELARG+DELARG1).AND.XKS.LT.-DELARG+DELARG1) THEN
  1979 + CALL R2OUTER(XSM,Y,ZSM,BXSM1,BY1,BZSM1)
  1980 + CALL R2SHEET(XSM,Y,ZSM,BXSM2,BY2,BZSM2)
  1981 + F2=-0.02*TKSI(XKS,-DELARG,DELARG1)
  1982 + F1=-0.02-F2
  1983 + BXSM=BXSM1*F1+BXSM2*F2
  1984 + BY=BY1*F1+BY2*F2
  1985 + BZSM=BZSM1*F1+BZSM2*F2
  1986 + ENDIF
  1987 +
  1988 + IF (XKS.GE.-DELARG+DELARG1.AND.XKS.LT.DELARG-DELARG1) THEN
  1989 + CALL R2SHEET(XSM,Y,ZSM,BXSM,BY,BZSM)
  1990 + BXSM=-BXSM*0.02
  1991 + BY=-BY*0.02
  1992 + BZSM=-BZSM*0.02
  1993 + ENDIF
  1994 + IF (XKS.GE.DELARG-DELARG1.AND.XKS.LT.DELARG+DELARG1) THEN
  1995 + CALL R2INNER(XSM,Y,ZSM,BXSM1,BY1,BZSM1)
  1996 + CALL R2SHEET(XSM,Y,ZSM,BXSM2,BY2,BZSM2)
  1997 + F1=-0.02*TKSI(XKS,DELARG,DELARG1)
  1998 + F2=-0.02-F1
  1999 + BXSM=BXSM1*F1+BXSM2*F2
  2000 + BY=BY1*F1+BY2*F2
  2001 + BZSM=BZSM1*F1+BZSM2*F2
  2002 + ENDIF
  2003 + IF (XKS.GE.DELARG+DELARG1) THEN
  2004 + CALL R2INNER(XSM,Y,ZSM,BXSM,BY,BZSM)
  2005 + BXSM=-BXSM*0.02
  2006 + BY=-BY*0.02
  2007 + BZSM=-BZSM*0.02
  2008 + ENDIF
  2009 +C
  2010 + BX=BXSM*CPS+BZSM*SPS
  2011 + BZ=BZSM*CPS-BXSM*SPS
  2012 +C
  2013 + RETURN
  2014 + END
  2015 +C
  2016 +C****************************************************************
  2017 +
  2018 +c
  2019 + SUBROUTINE R2INNER (X,Y,Z,BX,BY,BZ)
  2020 +C
  2021 +C
  2022 + IMPLICIT REAL*8 (A-H,O-Z)
  2023 + DIMENSION CBX(5),CBY(5),CBZ(5)
  2024 +C
  2025 + DATA PL1,PL2,PL3,PL4,PL5,PL6,PL7,PL8/154.185,-2.12446,.601735E-01,
  2026 + * -.153954E-02,.355077E-04,29.9996,262.886,99.9132/
  2027 + DATA PN1,PN2,PN3,PN4,PN5,PN6,PN7,PN8/-8.1902,6.5239,5.504,7.7815,
  2028 + * .8573,3.0986,.0774,-.038/
  2029 +C
  2030 + CALL BCONIC(X,Y,Z,CBX,CBY,CBZ,5)
  2031 +C
  2032 +C NOW INTRODUCE ONE 4-LOOP SYSTEM:
  2033 +C
  2034 + CALL LOOPS4(X,Y,Z,DBX8,DBY8,DBZ8,PN1,PN2,PN3,PN4,PN5,PN6)
  2035 +C
  2036 + CALL DIPDISTR(X-PN7,Y,Z,DBX6,DBY6,DBZ6,0)
  2037 + CALL DIPDISTR(X-PN8,Y,Z,DBX7,DBY7,DBZ7,1)
  2038 +
  2039 +C NOW COMPUTE THE FIELD COMPONENTS:
  2040 +
  2041 + BX=PL1*CBX(1)+PL2*CBX(2)+PL3*CBX(3)+PL4*CBX(4)+PL5*CBX(5)
  2042 + * +PL6*DBX6+PL7*DBX7+PL8*DBX8
  2043 + BY=PL1*CBY(1)+PL2*CBY(2)+PL3*CBY(3)+PL4*CBY(4)+PL5*CBY(5)
  2044 + * +PL6*DBY6+PL7*DBY7+PL8*DBY8
  2045 + BZ=PL1*CBZ(1)+PL2*CBZ(2)+PL3*CBZ(3)+PL4*CBZ(4)+PL5*CBZ(5)
  2046 + * +PL6*DBZ6+PL7*DBZ7+PL8*DBZ8
  2047 +C
  2048 + RETURN
  2049 + END
  2050 +C-----------------------------------------------------------------------
  2051 +
  2052 + SUBROUTINE BCONIC(X,Y,Z,CBX,CBY,CBZ,NMAX)
  2053 +C
  2054 +c "CONICAL" HARMONICS
  2055 +c
  2056 + IMPLICIT REAL*8 (A-H,O-Z)
  2057 +C
  2058 + DIMENSION CBX(NMAX),CBY(NMAX),CBZ(NMAX)
  2059 +
  2060 + RO2=X**2+Y**2
  2061 + RO=SQRT(RO2)
  2062 +C
  2063 + CF=X/RO
  2064 + SF=Y/RO
  2065 + CFM1=1.D0
  2066 + SFM1=0.D0
  2067 +C
  2068 + R2=RO2+Z**2
  2069 + R=DSQRT(R2)
  2070 + C=Z/R
  2071 + S=RO/R
  2072 + CH=DSQRT(0.5D0*(1.D0+C))
  2073 + SH=DSQRT(0.5D0*(1.D0-C))
  2074 + TNHM1=1.D0
  2075 + CNHM1=1.D0
  2076 + TNH=SH/CH
  2077 + CNH=1.D0/TNH
  2078 +C
  2079 + DO 1 M=1,NMAX
  2080 + CFM=CFM1*CF-SFM1*SF
  2081 + SFM=CFM1*SF+SFM1*CF
  2082 + CFM1=CFM
  2083 + SFM1=SFM
  2084 + TNHM=TNHM1*TNH
  2085 + CNHM=CNHM1*CNH
  2086 + BT=M*CFM/(R*S)*(TNHM+CNHM)
  2087 + BF=-0.5D0*M*SFM/R*(TNHM1/CH**2-CNHM1/SH**2)
  2088 + TNHM1=TNHM
  2089 + CNHM1=CNHM
  2090 + CBX(M)=BT*C*CF-BF*SF
  2091 + CBY(M)=BT*C*SF+BF*CF
  2092 + 1 CBZ(M)=-BT*S
  2093 +C
  2094 + RETURN
  2095 + END
  2096 +
  2097 +C-------------------------------------------------------------------
  2098 +C
  2099 + SUBROUTINE DIPDISTR(X,Y,Z,BX,BY,BZ,MODE)
  2100 +C
  2101 +C RETURNS FIELD COMPONENTS FROM A LINEAR DISTRIBUTION OF DIPOLAR SOURCES
  2102 +C ON THE Z-AXIS. THE PARAMETER MODE DEFINES HOW THE DIPOLE STRENGTH
  2103 +C VARIES ALONG THE Z-AXIS: MODE=0 IS FOR A STEP-FUNCTION (Mx=const &gt; 0
  2104 +c FOR Z &gt; 0, AND Mx=-const &lt; 0 FOR Z &lt; 0)
  2105 +C WHILE MODE=1 IS FOR A LINEAR VARIATION OF THE DIPOLE MOMENT DENSITY
  2106 +C SEE NB#3, PAGE 53 FOR DETAILS.
  2107 +C
  2108 +C
  2109 +C INPUT: X,Y,Z OF A POINT OF SPACE, AND MODE
  2110 +C
  2111 + IMPLICIT REAL*8 (A-H,O-Z)
  2112 + X2=X*X
  2113 + RHO2=X2+Y*Y
  2114 + R2=RHO2+Z*Z
  2115 + R3=R2*DSQRT(R2)
  2116 +
  2117 + IF (MODE.EQ.0) THEN
  2118 + BX=Z/RHO2**2*(R2*(Y*Y-X2)-RHO2*X2)/R3
  2119 + BY=-X*Y*Z/RHO2**2*(2.D0*R2+RHO2)/R3
  2120 + BZ=X/R3
  2121 + ELSE
  2122 + BX=Z/RHO2**2*(Y*Y-X2)
  2123 + BY=-2.D0*X*Y*Z/RHO2**2
  2124 + BZ=X/RHO2
  2125 + ENDIF
  2126 + RETURN
  2127 + END
  2128 +
  2129 +C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2130 +C
  2131 + SUBROUTINE R2OUTER (X,Y,Z,BX,BY,BZ)
  2132 +C
  2133 + IMPLICIT REAL*8 (A-H,O-Z)
  2134 +C
  2135 + DATA PL1,PL2,PL3,PL4,PL5/-34.105,-2.00019,628.639,73.4847,12.5162/
  2136 + DATA PN1,PN2,PN3,PN4,PN5,PN6,PN7,PN8,PN9,PN10,PN11,PN12,PN13,PN14,
  2137 + * PN15,PN16,PN17 /.55,.694,.0031,1.55,2.8,.1375,-.7,.2,.9625,
  2138 + * -2.994,2.925,-1.775,4.3,-.275,2.7,.4312,1.55/
  2139 +c
  2140 +C THREE PAIRS OF CROSSED LOOPS:
  2141 +C
  2142 + CALL CROSSLP(X,Y,Z,DBX1,DBY1,DBZ1,PN1,PN2,PN3)
  2143 + CALL CROSSLP(X,Y,Z,DBX2,DBY2,DBZ2,PN4,PN5,PN6)
  2144 + CALL CROSSLP(X,Y,Z,DBX3,DBY3,DBZ3,PN7,PN8,PN9)
  2145 +C
  2146 +C NOW AN EQUATORIAL LOOP ON THE NIGHTSIDE
  2147 +C
  2148 + CALL CIRCLE(X-PN10,Y,Z,PN11,DBX4,DBY4,DBZ4)
  2149 +c
  2150 +c NOW A 4-LOOP SYSTEM ON THE NIGHTSIDE
  2151 +c
  2152 +
  2153 + CALL LOOPS4(X,Y,Z,DBX5,DBY5,DBZ5,PN12,PN13,PN14,PN15,PN16,PN17)
  2154 +
  2155 +C---------------------------------------------------------------------
  2156 +
  2157 +C NOW COMPUTE THE FIELD COMPONENTS:
  2158 +
  2159 + BX=PL1*DBX1+PL2*DBX2+PL3*DBX3+PL4*DBX4+PL5*DBX5
  2160 + BY=PL1*DBY1+PL2*DBY2+PL3*DBY3+PL4*DBY4+PL5*DBY5
  2161 + BZ=PL1*DBZ1+PL2*DBZ2+PL3*DBZ3+PL4*DBZ4+PL5*DBZ5
  2162 +
  2163 + RETURN
  2164 + END
  2165 +C
  2166 +C--------------------------------------------------------------------
  2167 +C
  2168 + SUBROUTINE LOOPS4(X,Y,Z,BX,BY,BZ,XC,YC,ZC,R,THETA,PHI)
  2169 +C
  2170 +C RETURNS FIELD COMPONENTS FROM A SYSTEM OF 4 CURRENT LOOPS, POSITIONED
  2171 +C SYMMETRICALLY WITH RESPECT TO NOON-MIDNIGHT MERIDIAN AND EQUATORIAL
  2172 +C PLANES.
  2173 +C INPUT: X,Y,Z OF A POINT OF SPACE
  2174 +C XC,YC,ZC (YC &gt; 0 AND ZC &gt; 0) - POSITION OF THE CENTER OF THE
  2175 +C 1ST-QUADRANT LOOP
  2176 +C R - LOOP RADIUS (THE SAME FOR ALL FOUR)
  2177 +C THETA, PHI - SPECIFY THE ORIENTATION OF THE NORMAL OF THE 1ST LOOP
  2178 +c -----------------------------------------------------------
  2179 +
  2180 + IMPLICIT REAL*8 (A-H,O-Z)
  2181 +C
  2182 + CT=DCOS(THETA)
  2183 + ST=DSIN(THETA)
  2184 + CP=DCOS(PHI)
  2185 + SP=DSIN(PHI)
  2186 +C------------------------------------1ST QUADRANT:
  2187 + XS=(X-XC)*CP+(Y-YC)*SP
  2188 + YSS=(Y-YC)*CP-(X-XC)*SP
  2189 + ZS=Z-ZC
  2190 + XSS=XS*CT-ZS*ST
  2191 + ZSS=ZS*CT+XS*ST
  2192 +
  2193 + CALL CIRCLE(XSS,YSS,ZSS,R,BXSS,BYS,BZSS)
  2194 + BXS=BXSS*CT+BZSS*ST
  2195 + BZ1=BZSS*CT-BXSS*ST
  2196 + BX1=BXS*CP-BYS*SP
  2197 + BY1=BXS*SP+BYS*CP
  2198 +C-------------------------------------2nd QUADRANT:
  2199 + XS=(X-XC)*CP-(Y+YC)*SP
  2200 + YSS=(Y+YC)*CP+(X-XC)*SP
  2201 + ZS=Z-ZC
  2202 + XSS=XS*CT-ZS*ST
  2203 + ZSS=ZS*CT+XS*ST
  2204 +
  2205 + CALL CIRCLE(XSS,YSS,ZSS,R,BXSS,BYS,BZSS)
  2206 + BXS=BXSS*CT+BZSS*ST
  2207 + BZ2=BZSS*CT-BXSS*ST
  2208 + BX2=BXS*CP+BYS*SP
  2209 + BY2=-BXS*SP+BYS*CP
  2210 +C-------------------------------------3RD QUADRANT:
  2211 + XS=-(X-XC)*CP+(Y+YC)*SP
  2212 + YSS=-(Y+YC)*CP-(X-XC)*SP
  2213 + ZS=Z+ZC
  2214 + XSS=XS*CT-ZS*ST
  2215 + ZSS=ZS*CT+XS*ST
  2216 +
  2217 + CALL CIRCLE(XSS,YSS,ZSS,R,BXSS,BYS,BZSS)
  2218 + BXS=BXSS*CT+BZSS*ST
  2219 + BZ3=BZSS*CT-BXSS*ST
  2220 + BX3=-BXS*CP-BYS*SP
  2221 + BY3=BXS*SP-BYS*CP
  2222 +C-------------------------------------4TH QUADRANT:
  2223 + XS=-(X-XC)*CP-(Y-YC)*SP
  2224 + YSS=-(Y-YC)*CP+(X-XC)*SP
  2225 + ZS=Z+ZC
  2226 + XSS=XS*CT-ZS*ST
  2227 + ZSS=ZS*CT+XS*ST
  2228 +
  2229 + CALL CIRCLE(XSS,YSS,ZSS,R,BXSS,BYS,BZSS)
  2230 + BXS=BXSS*CT+BZSS*ST
  2231 + BZ4=BZSS*CT-BXSS*ST
  2232 + BX4=-BXS*CP+BYS*SP
  2233 + BY4=-BXS*SP-BYS*CP
  2234 +
  2235 + BX=BX1+BX2+BX3+BX4
  2236 + BY=BY1+BY2+BY3+BY4
  2237 + BZ=BZ1+BZ2+BZ3+BZ4
  2238 +
  2239 + RETURN
  2240 + END
  2241 +C
  2242 +C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
  2243 +C
  2244 + SUBROUTINE R2SHEET(X,Y,Z,BX,BY,BZ)
  2245 +C
  2246 + IMPLICIT REAL*8 (A-H,O-Z)
  2247 +C
  2248 + DATA PNONX1,PNONX2,PNONX3,PNONX4,PNONX5,PNONX6,PNONX7,PNONX8,
  2249 + * PNONY1,PNONY2,PNONY3,PNONY4,PNONY5,PNONY6,PNONY7,PNONY8,
  2250 + * PNONZ1,PNONZ2,PNONZ3,PNONZ4,PNONZ5,PNONZ6,PNONZ7,PNONZ8
  2251 + */-19.0969D0,-9.28828D0,-0.129687D0,5.58594D0,22.5055D0,
  2252 + * 0.483750D-01,0.396953D-01,0.579023D-01,-13.6750D0,-6.70625D0,
  2253 + * 2.31875D0,11.4062D0,20.4562D0,0.478750D-01,0.363750D-01,
  2254 + * 0.567500D-01,-16.7125D0,-16.4625D0,-0.1625D0,5.1D0,23.7125D0,
  2255 + * 0.355625D-01,0.318750D-01,0.538750D-01/
  2256 +C
  2257 +C
  2258 + DATA A1,A2,A3,A4,A5,A6,A7,A8,A9,A10,A11,A12,A13,A14,A15,A16,A17,
  2259 + * A18,A19,A20,A21,A22,A23,A24,A25,A26,A27,A28,A29,A30,A31,A32,A33,
  2260 + * A34,A35,A36,A37,A38,A39,A40,A41,A42,A43,A44,A45,A46,A47,A48,A49,
  2261 + * A50,A51,A52,A53,A54,A55,A56,A57,A58,A59,A60,A61,A62,A63,A64,A65,
  2262 + * A66,A67,A68,A69,A70,A71,A72,A73,A74,A75,A76,A77,A78,A79,A80
  2263 + * /8.07190D0,-7.39582D0,-7.62341D0,0.684671D0,-13.5672D0,11.6681D0,
  2264 + * 13.1154,-0.890217D0,7.78726D0,-5.38346D0,-8.08738D0,0.609385D0,
  2265 + * -2.70410D0, 3.53741D0,3.15549D0,-1.11069D0,-8.47555D0,0.278122D0,
  2266 + * 2.73514D0,4.55625D0,13.1134D0,1.15848D0,-3.52648D0,-8.24698D0,
  2267 + * -6.85710D0,-2.81369D0, 2.03795D0, 4.64383D0,2.49309D0,-1.22041D0,
  2268 + * -1.67432D0,-0.422526D0,-5.39796D0,7.10326D0,5.53730D0,-13.1918D0,
  2269 + * 4.67853D0,-7.60329D0,-2.53066D0, 7.76338D0, 5.60165D0,5.34816D0,
  2270 + * -4.56441D0,7.05976D0,-2.62723D0,-0.529078D0,1.42019D0,-2.93919D0,
  2271 + * 55.6338D0,-1.55181D0,39.8311D0,-80.6561D0,-46.9655D0,32.8925D0,
  2272 + * -6.32296D0,19.7841D0,124.731D0,10.4347D0,-30.7581D0,102.680D0,
  2273 + * -47.4037D0,-3.31278D0,9.37141D0,-50.0268D0,-533.319D0,110.426D0,
  2274 + * 1000.20D0,-1051.40D0, 1619.48D0,589.855D0,-1462.73D0,1087.10D0,
  2275 + * -1994.73D0,-1654.12D0,1263.33D0,-260.210D0,1424.84D0,1255.71D0,
  2276 + * -956.733D0, 219.946D0/
  2277 +C
  2278 +C
  2279 + DATA B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11,B12,B13,B14,B15,B16,B17,
  2280 + * B18,B19,B20,B21,B22,B23,B24,B25,B26,B27,B28,B29,B30,B31,B32,B33,
  2281 + * B34,B35,B36,B37,B38,B39,B40,B41,B42,B43,B44,B45,B46,B47,B48,B49,
  2282 + * B50,B51,B52,B53,B54,B55,B56,B57,B58,B59,B60,B61,B62,B63,B64,B65,
  2283 + * B66,B67,B68,B69,B70,B71,B72,B73,B74,B75,B76,B77,B78,B79,B80
  2284 + */-9.08427D0,10.6777D0,10.3288D0,-0.969987D0,6.45257D0,-8.42508D0,
  2285 + * -7.97464D0,1.41996D0,-1.92490D0,3.93575D0,2.83283D0,-1.48621D0,
  2286 + *0.244033D0,-0.757941D0,-0.386557D0,0.344566D0,9.56674D0,-2.5365D0,
  2287 + * -3.32916D0,-5.86712D0,-6.19625D0,1.83879D0,2.52772D0,4.34417D0,
  2288 + * 1.87268D0,-2.13213D0,-1.69134D0,-.176379D0,-.261359D0,.566419D0,
  2289 + * 0.3138D0,-0.134699D0,-3.83086D0,-8.4154D0,4.77005D0,-9.31479D0,
  2290 + * 37.5715D0,19.3992D0,-17.9582D0,36.4604D0,-14.9993D0,-3.1442D0,
  2291 + * 6.17409D0,-15.5519D0,2.28621D0,-0.891549D-2,-.462912D0,2.47314D0,
  2292 + * 41.7555D0,208.614D0,-45.7861D0,-77.8687D0,239.357D0,-67.9226D0,
  2293 + * 66.8743D0,238.534D0,-112.136D0,16.2069D0,-40.4706D0,-134.328D0,
  2294 + * 21.56D0,-0.201725D0,2.21D0,32.5855D0,-108.217D0,-1005.98D0,
  2295 + * 585.753D0,323.668D0,-817.056D0,235.750D0,-560.965D0,-576.892D0,
  2296 + * 684.193D0,85.0275D0,168.394D0,477.776D0,-289.253D0,-123.216D0,
  2297 + * 75.6501D0,-178.605D0/
  2298 +C
  2299 + DATA C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12,C13,C14,C15,C16,C17,
  2300 + * C18,C19,C20,C21,C22,C23,C24,C25,C26,C27,C28,C29,C30,C31,C32,C33,
  2301 + * C34,C35,C36,C37,C38,C39,C40,C41,C42,C43,C44,C45,C46,C47,C48,C49,
  2302 + * C50,C51,C52,C53,C54,C55,C56,C57,C58,C59,C60,C61,C62,C63,C64,C65,
  2303 + * C66,C67,C68,C69,C70,C71,C72,C73,C74,C75,C76,C77,C78,C79,C80
  2304 + * / 1167.61D0,-917.782D0,-1253.2D0,-274.128D0,-1538.75D0,1257.62D0,
  2305 + * 1745.07D0,113.479D0,393.326D0,-426.858D0,-641.1D0,190.833D0,
  2306 + * -29.9435D0,-1.04881D0,117.125D0,-25.7663D0,-1168.16D0,910.247D0,
  2307 + * 1239.31D0,289.515D0,1540.56D0,-1248.29D0,-1727.61D0,-131.785D0,
  2308 + * -394.577D0,426.163D0,637.422D0,-187.965D0,30.0348D0,0.221898D0,
  2309 + * -116.68D0,26.0291D0,12.6804D0,4.84091D0,1.18166D0,-2.75946D0,
  2310 + * -17.9822D0,-6.80357D0,-1.47134D0,3.02266D0,4.79648D0,0.665255D0,
  2311 + * -0.256229D0,-0.857282D-1,-0.588997D0,0.634812D-1,0.164303D0,
  2312 + * -0.15285D0,22.2524D0,-22.4376D0,-3.85595D0,6.07625D0,-105.959D0,
  2313 + * -41.6698D0,0.378615D0,1.55958D0,44.3981D0,18.8521D0,3.19466D0,
  2314 + * 5.89142D0,-8.63227D0,-2.36418D0,-1.027D0,-2.31515D0,1035.38D0,
  2315 + * 2040.66D0,-131.881D0,-744.533D0,-3274.93D0,-4845.61D0,482.438D0,
  2316 + * 1567.43D0,1354.02D0,2040.47D0,-151.653D0,-845.012D0,-111.723D0,
  2317 + * -265.343D0,-26.1171D0,216.632D0/
  2318 +C
  2319 +c------------------------------------------------------------------
  2320 +C
  2321 + XKS=XKSI(X,Y,Z) ! variation across the current sheet
  2322 + T1X=XKS/DSQRT(XKS**2+PNONX6**2)
  2323 + T2X=PNONX7**3/DSQRT(XKS**2+PNONX7**2)**3
  2324 + T3X=XKS/DSQRT(XKS**2+PNONX8**2)**5 *3.493856D0*PNONX8**4
  2325 +C
  2326 + T1Y=XKS/DSQRT(XKS**2+PNONY6**2)
  2327 + T2Y=PNONY7**3/DSQRT(XKS**2+PNONY7**2)**3
  2328 + T3Y=XKS/DSQRT(XKS**2+PNONY8**2)**5 *3.493856D0*PNONY8**4
  2329 +C
  2330 + T1Z=XKS/DSQRT(XKS**2+PNONZ6**2)
  2331 + T2Z=PNONZ7**3/DSQRT(XKS**2+PNONZ7**2)**3
  2332 + T3Z=XKS/DSQRT(XKS**2+PNONZ8**2)**5 *3.493856D0*PNONZ8**4
  2333 +C
  2334 + RHO2=X*X+Y*Y
  2335 + R=DSQRT(RHO2+Z*Z)
  2336 + RHO=DSQRT(RHO2)
  2337 +C
  2338 + C1P=X/RHO
  2339 + S1P=Y/RHO
  2340 + S2P=2.D0*S1P*C1P
  2341 + C2P=C1P*C1P-S1P*S1P
  2342 + S3P=S2P*C1P+C2P*S1P
  2343 + C3P=C2P*C1P-S2P*S1P
  2344 + S4P=S3P*C1P+C3P*S1P
  2345 + CT=Z/R
  2346 + ST=RHO/R
  2347 +C
  2348 + S1=FEXP(CT,PNONX1)
  2349 + S2=FEXP(CT,PNONX2)
  2350 + S3=FEXP(CT,PNONX3)
  2351 + S4=FEXP(CT,PNONX4)
  2352 + S5=FEXP(CT,PNONX5)
  2353 +C
  2354 +C NOW COMPUTE THE GSM FIELD COMPONENTS:
  2355 +C
  2356 +C
  2357 + BX=S1*((A1+A2*T1X+A3*T2X+A4*T3X)
  2358 + * +C1P*(A5+A6*T1X+A7*T2X+A8*T3X)
  2359 + * +C2P*(A9+A10*T1X+A11*T2X+A12*T3X)
  2360 + * +C3P*(A13+A14*T1X+A15*T2X+A16*T3X))
  2361 + * +S2*((A17+A18*T1X+A19*T2X+A20*T3X)
  2362 + * +C1P*(A21+A22*T1X+A23*T2X+A24*T3X)
  2363 + * +C2P*(A25+A26*T1X+A27*T2X+A28*T3X)
  2364 + * +C3P*(A29+A30*T1X+A31*T2X+A32*T3X))
  2365 + * +S3*((A33+A34*T1X+A35*T2X+A36*T3X)
  2366 + * +C1P*(A37+A38*T1X+A39*T2X+A40*T3X)
  2367 + * +C2P*(A41+A42*T1X+A43*T2X+A44*T3X)
  2368 + * +C3P*(A45+A46*T1X+A47*T2X+A48*T3X))
  2369 + * +S4*((A49+A50*T1X+A51*T2X+A52*T3X)
  2370 + * +C1P*(A53+A54*T1X+A55*T2X+A56*T3X)
  2371 + * +C2P*(A57+A58*T1X+A59*T2X+A60*T3X)
  2372 + * +C3P*(A61+A62*T1X+A63*T2X+A64*T3X))
  2373 + * +S5*((A65+A66*T1X+A67*T2X+A68*T3X)
  2374 + * +C1P*(A69+A70*T1X+A71*T2X+A72*T3X)
  2375 + * +C2P*(A73+A74*T1X+A75*T2X+A76*T3X)
  2376 + * +C3P*(A77+A78*T1X+A79*T2X+A80*T3X))
  2377 +C
  2378 +C
  2379 + S1=FEXP(CT,PNONY1)
  2380 + S2=FEXP(CT,PNONY2)
  2381 + S3=FEXP(CT,PNONY3)
  2382 + S4=FEXP(CT,PNONY4)
  2383 + S5=FEXP(CT,PNONY5)
  2384 +C
  2385 + BY=S1*(S1P*(B1+B2*T1Y+B3*T2Y+B4*T3Y)
  2386 + * +S2P*(B5+B6*T1Y+B7*T2Y+B8*T3Y)
  2387 + * +S3P*(B9+B10*T1Y+B11*T2Y+B12*T3Y)
  2388 + * +S4P*(B13+B14*T1Y+B15*T2Y+B16*T3Y))
  2389 + * +S2*(S1P*(B17+B18*T1Y+B19*T2Y+B20*T3Y)
  2390 + * +S2P*(B21+B22*T1Y+B23*T2Y+B24*T3Y)
  2391 + * +S3P*(B25+B26*T1Y+B27*T2Y+B28*T3Y)
  2392 + * +S4P*(B29+B30*T1Y+B31*T2Y+B32*T3Y))
  2393 + * +S3*(S1P*(B33+B34*T1Y+B35*T2Y+B36*T3Y)
  2394 + * +S2P*(B37+B38*T1Y+B39*T2Y+B40*T3Y)
  2395 + * +S3P*(B41+B42*T1Y+B43*T2Y+B44*T3Y)
  2396 + * +S4P*(B45+B46*T1Y+B47*T2Y+B48*T3Y))
  2397 + * +S4*(S1P*(B49+B50*T1Y+B51*T2Y+B52*T3Y)
  2398 + * +S2P*(B53+B54*T1Y+B55*T2Y+B56*T3Y)
  2399 + * +S3P*(B57+B58*T1Y+B59*T2Y+B60*T3Y)
  2400 + * +S4P*(B61+B62*T1Y+B63*T2Y+B64*T3Y))
  2401 + * +S5*(S1P*(B65+B66*T1Y+B67*T2Y+B68*T3Y)
  2402 + * +S2P*(B69+B70*T1Y+B71*T2Y+B72*T3Y)
  2403 + * +S3P*(B73+B74*T1Y+B75*T2Y+B76*T3Y)
  2404 + * +S4P*(B77+B78*T1Y+B79*T2Y+B80*T3Y))
  2405 +C
  2406 + S1=FEXP1(CT,PNONZ1)
  2407 + S2=FEXP1(CT,PNONZ2)
  2408 + S3=FEXP1(CT,PNONZ3)
  2409 + S4=FEXP1(CT,PNONZ4)
  2410 + S5=FEXP1(CT,PNONZ5)
  2411 +C
  2412 + BZ=S1*((C1+C2*T1Z+C3*T2Z+C4*T3Z)
  2413 + * +C1P*(C5+C6*T1Z+C7*T2Z+C8*T3Z)
  2414 + * +C2P*(C9+C10*T1Z+C11*T2Z+C12*T3Z)
  2415 + * +C3P*(C13+C14*T1Z+C15*T2Z+C16*T3Z))
  2416 + * +S2*((C17+C18*T1Z+C19*T2Z+C20*T3Z)
  2417 + * +C1P*(C21+C22*T1Z+C23*T2Z+C24*T3Z)
  2418 + * +C2P*(C25+C26*T1Z+C27*T2Z+C28*T3Z)
  2419 + * +C3P*(C29+C30*T1Z+C31*T2Z+C32*T3Z))
  2420 + * +S3*((C33+C34*T1Z+C35*T2Z+C36*T3Z)
  2421 + * +C1P*(C37+C38*T1Z+C39*T2Z+C40*T3Z)
  2422 + * +C2P*(C41+C42*T1Z+C43*T2Z+C44*T3Z)
  2423 + * +C3P*(C45+C46*T1Z+C47*T2Z+C48*T3Z))
  2424 + * +S4*((C49+C50*T1Z+C51*T2Z+C52*T3Z)
  2425 + * +C1P*(C53+C54*T1Z+C55*T2Z+C56*T3Z)
  2426 + * +C2P*(C57+C58*T1Z+C59*T2Z+C60*T3Z)
  2427 + * +C3P*(C61+C62*T1Z+C63*T2Z+C64*T3Z))
  2428 + * +S5*((C65+C66*T1Z+C67*T2Z+C68*T3Z)
  2429 + * +C1P*(C69+C70*T1Z+C71*T2Z+C72*T3Z)
  2430 + * +C2P*(C73+C74*T1Z+C75*T2Z+C76*T3Z)
  2431 + * +C3P*(C77+C78*T1Z+C79*T2Z+C80*T3Z))
  2432 +C
  2433 + RETURN
  2434 + END
  2435 +C
  2436 +C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  2437 +
  2438 + DOUBLE PRECISION FUNCTION XKSI(X,Y,Z)
  2439 +C
  2440 + IMPLICIT REAL*8 (A-H,O-Z)
  2441 +C
  2442 +C A11 - C72, R0, and DR below ARE STRETCH PARAMETERS (P.26-27, NB# 3),
  2443 +C
  2444 + DATA A11A12,A21A22,A41A42,A51A52,A61A62,B11B12,B21B22,C61C62,
  2445 + * C71C72,R0,DR /0.305662,-0.383593,0.2677733,-0.097656,-0.636034,
  2446 + * -0.359862,0.424706,-0.126366,0.292578,1.21563,7.50937/
  2447 +
  2448 + DATA TNOON,DTETA/0.3665191,0.09599309/ ! Correspond to noon and midnight
  2449 +C latitudes 69 and 63.5 degs, resp.
  2450 + DR2=DR*DR
  2451 +C
  2452 + X2=X*X
  2453 + Y2=Y*Y
  2454 + Z2=Z*Z
  2455 + XY=X*Y
  2456 + XYZ=XY*Z
  2457 + R2=X2+Y2+Z2
  2458 + R=DSQRT(R2)
  2459 + R3=R2*R
  2460 + R4=R2*R2
  2461 + XR=X/R
  2462 + YR=Y/R
  2463 + ZR=Z/R
  2464 +C
  2465 + IF (R.LT.R0) THEN
  2466 + PR=0.D0
  2467 + ELSE
  2468 + PR=DSQRT((R-R0)**2+DR2)-DR
  2469 + ENDIF
  2470 +C
  2471 + F=X+PR*(A11A12+A21A22*XR+A41A42*XR*XR+A51A52*YR*YR+
  2472 + * A61A62*ZR*ZR)
  2473 + G=Y+PR*(B11B12*YR+B21B22*XR*YR)
  2474 + H=Z+PR*(C61C62*ZR+C71C72*XR*ZR)
  2475 + G2=G*G
  2476 +C
  2477 + FGH=F**2+G2+H**2
  2478 + FGH32=DSQRT(FGH)**3
  2479 + FCHSG2=F**2+G2
  2480 +
  2481 + IF (FCHSG2.LT.1.D-5) THEN
  2482 + XKSI=-1.D0 ! THIS IS JUST FOR ELIMINATING PROBLEMS
  2483 + RETURN ! ON THE Z-AXIS
  2484 + ENDIF
  2485 +
  2486 + SQFCHSG2=DSQRT(FCHSG2)
  2487 + ALPHA=FCHSG2/FGH32
  2488 + THETA=TNOON+0.5D0*DTETA*(1.D0-F/SQFCHSG2)
  2489 + PHI=DSIN(THETA)**2
  2490 +C
  2491 + XKSI=ALPHA-PHI
  2492 +C
  2493 + RETURN
  2494 + END
  2495 +C
  2496 +C--------------------------------------------------------------------
  2497 +C
  2498 + FUNCTION FEXP(S,A)
  2499 + IMPLICIT REAL*8 (A-H,O-Z)
  2500 + DATA E/2.718281828459D0/
  2501 + IF (A.LT.0.D0) FEXP=DSQRT(-2.D0*A*E)*S*DEXP(A*S*S)
  2502 + IF (A.GE.0.D0) FEXP=S*DEXP(A*(S*S-1.D0))
  2503 + RETURN
  2504 + END
  2505 +C
  2506 +C-----------------------------------------------------------------------
  2507 + FUNCTION FEXP1(S,A)
  2508 + IMPLICIT REAL*8 (A-H,O-Z)
  2509 + IF (A.LE.0.D0) FEXP1=DEXP(A*S*S)
  2510 + IF (A.GT.0.D0) FEXP1=DEXP(A*(S*S-1.D0))
  2511 + RETURN
  2512 + END
  2513 +C
  2514 +C************************************************************************
  2515 +C
  2516 + DOUBLE PRECISION FUNCTION TKSI(XKSI,XKS0,DXKSI)
  2517 + IMPLICIT REAL*8 (A-H,O-Z)
  2518 + SAVE M,TDZ3
  2519 + DATA M/0/
  2520 +C
  2521 + IF (M.EQ.0) THEN
  2522 + TDZ3=2.*DXKSI**3
  2523 + M=1
  2524 + ENDIF
  2525 +C
  2526 + IF (XKSI-XKS0.LT.-DXKSI) TKSII=0.
  2527 + IF (XKSI-XKS0.GE.DXKSI) TKSII=1.
  2528 +C
  2529 + IF (XKSI.GE.XKS0-DXKSI.AND.XKSI.LT.XKS0) THEN
  2530 + BR3=(XKSI-XKS0+DXKSI)**3
  2531 + TKSII=1.5*BR3/(TDZ3+BR3)
  2532 + ENDIF
  2533 +C
  2534 + IF (XKSI.GE.XKS0.AND.XKSI.LT.XKS0+DXKSI) THEN
  2535 + BR3=(XKSI-XKS0-DXKSI)**3
  2536 + TKSII=1.+1.5*BR3/(TDZ3-BR3)
  2537 + ENDIF
  2538 + TKSI=TKSII
  2539 + END
  2540 +C
  2541 +C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  2542 +c
  2543 + SUBROUTINE DIPOLE(PS,X,Y,Z,BX,BY,BZ)
  2544 +C
  2545 +C CALCULATES GSM COMPONENTS OF GEODIPOLE FIELD WITH THE DIPOLE MOMENT
  2546 +C CORRESPONDING TO THE EPOCH OF 1980.
  2547 +C------------INPUT PARAMETERS:
  2548 +C PS - GEODIPOLE TILT ANGLE IN RADIANS, X,Y,Z - GSM COORDINATES IN RE
  2549 +C------------OUTPUT PARAMETERS:
  2550 +C BX,BY,BZ - FIELD COMPONENTS IN GSM SYSTEM, IN NANOTESLA.
  2551 +C
  2552 +C
  2553 +C WRITEN BY: N. A. TSYGANENKO
  2554 +
  2555 + DATA M,PSI/0,5./
  2556 + SAVE M,PSI,SPS,CPS
  2557 + IF(M.EQ.1.AND.ABS(PS-PSI).LT.1.E-5) GOTO 1
  2558 + SPS=SIN(PS)
  2559 + CPS=COS(PS)
  2560 + PSI=PS
  2561 + M=1
  2562 + 1 P=X**2
  2563 + U=Z**2
  2564 + V=3.*Z*X
  2565 + T=Y**2
  2566 + Q=30574./SQRT(P+T+U)**5
  2567 + BX=Q*((T+U-2.*P)*SPS-V*CPS)
  2568 + BY=-3.*Y*Q*(X*SPS+Z*CPS)
  2569 + BZ=Q*((P+T-2.*U)*CPS-V*SPS)
  2570 + RETURN
  2571 + END
0 \ No newline at end of file 2572 \ No newline at end of file
geopack/src/include/geopack.hh 0 → 100644
  Diff comments   View file @55c4dfc
@@ -0,0 +1,29 @@ @@ -0,0 +1,29 @@
  1 +#ifndef GEOPACK_HH
  2 +#define GEOPACK_HH
  3 +
  4 +#ifdef __cplusplus
  5 +extern "C" {
  6 +#endif
  7 +
  8 + void recalc_08_(int* iyr, int *iday, int* ihour, int* min, int* isec,
  9 + float* V_GSE_X_i, float* V_GSE_Y_i, float* V_GSE_Z_i);
  10 +
  11 + void gswgse_08_(float* sat_pos_GSW_X_i, float* sat_pos_GSW_Y_i, float* sat_pos_GSW_Z_i,
  12 + float* sat_pos_GSE_X_i, float* sat_pos_GSE_Y_i, float* sat_pos_GSE_Z_i, int* transform_flag);
  13 +
  14 + void geogsw_08_(float* pos_GEO_X_i, float* pos_GEO_Y_i, float* pos_GEO_Z_i,
  15 + float* pos_GSW_X_i, float* pos_GSW_Y_i, float* pos_GSW_Z_i, int* transform_flag);
  16 +
  17 + void t96_mgnp_08_(float* XN_PD, float* VEL, float* sat_pos_GSW_X_i, float* sat_pos_GSW_Y_i, float* sat_pos_GSW_Z_i,
  18 + float* x_mgnp, float* y_mgnp, float* z_mgnp, float* DIST_MGNP, int* ID_MGNP);
  19 +
  20 + void trace_08_modified_(float* sat_pos_GSW_X_i, float* sat_pos_GSW_Y_i, float* sat_pos_GSW_Z_i, float* DIR,
  21 + float* DSMAX, float* ERR, float* RLIM, float* R0, int* IOPT, float PARMOD[],
  22 + float* FP_GSW_X_i, float* FP_GSW_Y_i, float* FP_GSW_Z_i,
  23 + float XX[], float YY[], float ZZ[], int* L, int* LMAX);
  24 +
  25 +#ifdef __cplusplus
  26 +}
  27 +#endif
  28 +
  29 +#endif
0 \ No newline at end of file 30 \ No newline at end of file