##// END OF EJS Templates
DATABASES » madrigal
RSS
Fixing plots
Fixing plots
jespinoza - - Load All Authors

File last commit:

r0:b84e1135c2c4
r7:4e0b343b0c61
Show More
Geopack_2008.f
2286 lines | 83.7 KiB | text/x-fortran | FortranFixedLexer
/ source / madf / geolib / Geopack_2008.f
History | Source | Raw |Copy content |Copy permalink
jespinoza
Initial
 r0 c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
C Downloaded from ftp://hanna.ccmc.gsfc.nasa.gov/pub/modelweb/magnetospheric/tsyganenko
C by B. Rideout on Mar. 29, 2013.
C The following changes were made from imported version:
C
C 1. ran script modifyFortran.py to remove
C to avoid conflicts with geo-cgm
C 2. remake_file script (nag_apt + tcl) run
C 3. All print statements commented out
C 4. Some _08 methods renamed to drop _08 due to compiling issues with double underscores
C 5. added common block br_err to signal errors back to trace subroutine
c
c ##########################################################################
c # #
c # GEOPACK-2008 #
c # (MAIN SET OF FORTRAN CODES) #
c # (IGRF-11 coefficients added on Nov 30, 2010) #
c ##########################################################################
C
c
cThis collection of subroutines is a result of several upgrades of the original package
c written by N. A. Tsyganenko in 1978-1979.
c
c PREFATORY NOTE TO THE VERSION OF FEBRUARY 8, 2008:
c
cTo avoid inappropriate use of obsolete subroutines from earlier versions, a suffix 08 was
c added to the name of each subroutine in this release.
c
cA possibility has been added in this version to calculate vector components in the
c"Geocentric Solar Wind" (GSW) coordinate system, which, to our knowledge, was first
cintroduced by Hones et al., Planet. Space Sci., v.34, p.889, 1986 (aka GSWM, see Appendix,
cTsyganenko et al., JGRA, v.103(A4), p.6827, 1998). The GSW system is analogous to the
cstandard GSM, except that its X-axis is antiparallel to the currently observed solar wind
cflow vector, rather than aligned with the Earth-Sun line. The orientation of axes in the
cGSW system can be uniquely defined by specifying three components (VGSEX,VGSEY,VGSEZ) of
cthe solar wind velocity, and in the case of a strictly radial anti-sunward flow (VGSEY=
cVGSEZ=0) the GSW system becomes identical to the standard GSM, which fact was used here
cto minimize the number of subroutines in the package. To that end, instead of the special
ccase of the GSM coordinates, this version uses a more general GSW system, and three more
cinput parameters are added in the subroutine TS_RECALC_08, the observed values (VGSEX,VGSEY,
cVGSEZ) of the solar wind velocity. Invoking TS_RECALC_08 with VGSEY=VGSEZ=0 restores the
cstandard (sunward) orientation of the X axis, which allows one to easily convert vectors
cbetween GSW and GSM, as well as to/from other standard and commonly used systems. For more
c details, see the documentation file GEOPACK-2008.DOC.
c
cAnother modification allows users to have more control over the procedure of field line
cmapping using the subroutine TRACE. To that end, three new input parameters were added
cin that subroutine, allowing one to set (i) an upper limit, DSMAX, on the automatically
cadjusted step size, (ii) a permissible step error, ERR, and (iii) maximal length, LMAX,
cof arrays where field line point coordinates are stored. Minor changes were also made in
cthe tracing subroutine, to make it more compact and easier for understanding, and to
cprevent the algorithm from making uncontrollable large number of multiple loops in some
c cases with plasmoid-like field structures.
c
COne more subroutine, named GEODGEO_08, was added to the package, allowing one to convert
cgeodetic coordinates of a point in space (altitude above the Earth's WGS84 ellipsoid and
cgeodetic latitude) to geocentric radial distance and colatitude, and vice versa.
c
CFor a complete list of modifications made earlier in previous versions, see the
c documentation file GEOPACK-2008.DOC.
c
c----------------------------------------------------------------------------------
c
SUBROUTINE IGRF_GSW_08(XGSW,YGSW,ZGSW,HXGSW,HYGSW,HZGSW)
c
CCALCULATES COMPONENTS OF THE MAIN (INTERNAL) TS_GEOMAGNETIC FIELD IN THE GEOCENTRIC SOLAR-WIND
C(GSW) COORDINATE SYSTEM, USING IAGA INTERNATIONAL TS_GEOMAGNETIC REFERENCE MODEL COEFFICIENTS
C(e.g., http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html, revised 22 March, 2005)
c
CTHE GSW SYSTEM IS ESSENTIALLY SIMILAR TO THE STANDARD GSM (THE TWO SYSTEMS BECOME IDENTICAL
CTO EACH OTHER IN THE CASE OF STRICTLY ANTI-TS_SUNWARD SOLAR WIND FLOW). FOR A DETAILED
C DEFINITION, SEE INTRODUCTORY COMMENTS FOR THE SUBROUTINE GSWGSE .
C
CBEFORE THE FIRST CALL OF THIS SUBROUTINE, OR, IF THE DATE/TIME (IYEAR,IDAY,IHOUR,MIN,ISEC),
COR THE SOLAR WIND VELOCITY COMPONENTS (VGSEX,VGSEY,VGSEZ) HAVE CHANGED, THE MODEL COEFFICIENTS
cAND GEO-GSW ROTATION MATRIX ELEMENTS SHOULD BE UPDATED BY CALLING THE SUBROUTINE TS_RECALC_08.
C
C -----INPUT PARAMETERS:
C
CXGSW,YGSW,ZGSW - CARTESIAN GEOCENTRIC SOLAR-WIND COORDINATES (IN UNITS RE=6371.2 KM)
C
C -----OUTPUT PARAMETERS:
C
CHXGSW,HYGSW,HZGSW - CARTESIAN GEOCENTRIC SOLAR-WIND COMPONENTS OF THE MAIN TS_GEOMAGNETIC
C FIELD IN NANOTESLA
C
C LAST MODIFICATION: FEB 07, 2008.
C THIS VERSION OF THE CODE ACCEPTS DATES FROM 1965 THROUGH 2015.
c
C AUTHOR: N. A. TSYGANENKO
C
C
C .. Scalar Arguments ..
DOUBLE PRECISION HXGSW,HYGSW,HZGSW,XGSW,YGSW,ZGSW
DOUBLE PRECISION XGEO,YGEO,ZGEO,R,THETA,PHI,COLAT,ELONG
DOUBLE PRECISION BN,BE,BZ,B
C .. common block added by B. Rideout to get date from RECALC_08
COMMON /BR_FYEAR/FYEAR
DOUBLE PRECISION FYEAR
C ..
C Rewritten by Bill Rideout. Now calls call igrf13.f methods
C
C convert to position to X,Y,Z
CALL GEOGSW(XGEO,YGEO,ZGEO,XGSW,YGSW,ZGSW,-1)
C convert to spherical position
CALL TS_SPHCAR_08(R,THETA,PHI,XGEO,YGEO,ZGEO,-1)
C convert to R, COLAT, ELONG
R = R * 6378.1
COLAT = THETA * 57.295779
ELONG = PHI * 57.295779
C call igrf
CALL igrf13syn(0,FYEAR,2,R,COLAT,ELONG,BN,BE,BZ,B)
BZ = -1.0 * BZ
C convert to X,Y,Z coordinates
CALL TS_BSPCAR_08(THETA,PHI,BZ,-1.0*BN,BE,XGEO,YGEO,ZGEO)
CALL GEOGSW(XGEO,YGEO,ZGEO,HXGSW,HYGSW,HZGSW,1)
RETURN
END
C
c==========================================================================================
c
SUBROUTINE IGRF_GEO_08(R,THETA,PHI,BR,BTHETA,BPHI)
c
CCALCULATES COMPONENTS OF THE MAIN (INTERNAL) TS_GEOMAGNETIC FIELD IN THE SPHERICAL GEOGRAPHIC
C(GEOCENTRIC) COORDINATE SYSTEM, USING IAGA INTERNATIONAL TS_GEOMAGNETIC REFERENCE MODEL
CCOEFFICIENTS (e.g., http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html, revised: 22 March, 2005)
C
CBEFORE THE FIRST CALL OF THIS SUBROUTINE, OR IF THE DATE (IYEAR AND IDAY) WAS CHANGED,
CTHE MODEL COEFFICIENTS SHOULD BE UPDATED BY CALLING THE SUBROUTINE TS_RECALC_08
C
C -----INPUT PARAMETERS:
C
C R, THETA, PHI - SPHERICAL GEOGRAPHIC (GEOCENTRIC) COORDINATES:
CRADIAL DISTANCE R IN UNITS RE=6371.2 KM, COLATITUDE THETA AND LONGITUDE PHI IN RADIANS
C
C -----OUTPUT PARAMETERS:
C
CBR, BTHETA, BPHI - SPHERICAL COMPONENTS OF THE MAIN TS_GEOMAGNETIC FIELD IN NANOTESLA
C (POSITIVE BR OUTWARD, BTHETA SOUTHWARD, BPHI EASTWARD)
C
C LAST MODIFICATION: MAY 4, 2005.
C THIS VERSION OF THE CODE ACCEPTS DATES FROM 1965 THROUGH 2015.
c
C AUTHOR: N. A. TSYGANENKO
C
C
C .. Scalar Arguments ..
DOUBLE PRECISION BPHI,BR,BTHETA,PHI,R,THETA
C ..
C .. Local Scalars ..
DOUBLE PRECISION AN,BBF,BBR,BBT,BI,C,CF,D,D2,DP,E,HH,P,P2,PM,PP,Q,
* QQ,S,SF,W,X,XK,Y,Z
INTEGER IRP3,K,M,MM,MN,N,NM
C ..
C .. Local Arrays ..
DOUBLE PRECISION A(14),B(14)
C ..
C .. Intrinsic Functions ..
INTRINSIC COS,SIN
C ..
C .. Common blocks ..
COMMON /GEOPACK2/G,H,REC
DOUBLE PRECISION G(105),H(105),REC(105)
C ..
C = COS(THETA)
S = SIN(THETA)
CF = COS(PHI)
SF = SIN(PHI)
C
PP = 1.D0/R
P = PP
C
CIN THIS NEW VERSION, THE OPTIMAL VALUE OF THE PARAMETER NM (MAXIMAL ORDER OF THE SPHERICAL
CHARMONIC EXPANSION) IS NOT USER-PRESCRIBED, BUT CALCULATED INSIDE THE SUBROUTINE, BASED
C ON THE VALUE OF THE RADIAL DISTANCE R:
C
IRP3 = R + 2
NM = 3 + 30/IRP3
IF (NM.GT.13) NM = 13
K = NM + 1
DO 10 N = 1,K
P = P*PP
A(N) = P
B(N) = P*N
10 CONTINUE
P = 1.D0
D = 0.D0
BBR = 0.D0
BBT = 0.D0
BBF = 0.D0
DO 60 M = 1,K
IF (M.EQ.1) GO TO 20
MM = M - 1
W = X
X = W*CF + Y*SF
Y = Y*CF - W*SF
GO TO 30
20 X = 0.D0
Y = 1.D0
30 Q = P
Z = D
BI = 0.D0
P2 = 0.D0
D2 = 0.D0
DO 50 N = M,K
AN = A(N)
MN = N*(N-1)/2 + M
E = G(MN)
HH = H(MN)
W = E*Y + HH*X
BBR = BBR + B(N)*W*Q
BBT = BBT - AN*W*Z
IF (M.EQ.1) GO TO 40
QQ = Q
IF (S.LT.1.D-5) QQ = Z
BI = BI + AN*(E*X-HH*Y)*QQ
40 XK = REC(MN)
DP = C*Z - S*Q - XK*D2
PM = C*Q - XK*P2
D2 = Z
P2 = Q
Z = DP
Q = PM
50 CONTINUE
D = S*D + C*P
P = S*P
IF (M.EQ.1) GO TO 60
BI = BI*MM
BBF = BBF + BI
60 CONTINUE
C
BR = BBR
BTHETA = BBT
IF (S.LT.1.D-5) GO TO 70
BPHI = BBF/S
RETURN
70 IF (C.LT.0.D0) BBF = -BBF
BPHI = BBF
RETURN
END
C
c==========================================================================================
c
SUBROUTINE DIP_08(XGSW,YGSW,ZGSW,BXGSW,BYGSW,BZGSW)
C
CCALCULATES GSW (GEOCENTRIC SOLAR-WIND) COMPONENTS OF GEODIPOLE FIELD WITH THE DIPOLE MOMENT
CCORRESPONDING TO THE EPOCH, SPECIFIED BY CALLING SUBROUTINE TS_RECALC_08 (SHOULD BE
CINVOKED BEFORE THE FIRST USE OF THIS ONE, OR IF THE DATE/TIME, AND/OR THE OBSERVED
C SOLAR WIND DIRECTION, HAVE CHANGED.
C
CTHE GSW COORDINATE SYSTEM IS ESSENTIALLY SIMILAR TO THE STANDARD GSM (THE TWO SYSTEMS BECOME
CIDENTICAL TO EACH OTHER IN THE CASE OF STRICTLY RADIAL ANTI-TS_SUNWARD SOLAR WIND FLOW). ITS
CDETAILED DEFINITION IS GIVEN IN INTRODUCTORY COMMENTS FOR THE SUBROUTINE GSWGSE .
C--INPUT PARAMETERS: XGSW,YGSW,ZGSW - GSW COORDINATES IN RE (1 RE = 6371.2 km)
C
C--OUTPUT PARAMETERS: BXGSW,BYGSW,BZGSW - FIELD COMPONENTS IN GSW SYSTEM, IN NANOTESLA.
C
C LAST MODIFICATION: JAN 28, 2008.
C
C AUTHOR: N. A. TSYGANENKO
C
C
C .. Scalar Arguments ..
DOUBLE PRECISION BXGSW,BYGSW,BZGSW,XGSW,YGSW,ZGSW
C ..
C .. Local Scalars ..
DOUBLE PRECISION DIPMOM,P,Q,T,U,V
C ..
C .. Intrinsic Functions ..
INTRINSIC SQRT
C ..
C .. Common blocks ..
COMMON /GEOPACK1/AA,SPS,CPS,BB
COMMON /GEOPACK2/G,H,REC
DOUBLE PRECISION CPS,SPS
DOUBLE PRECISION AA(10),BB(22),G(105),H(105),REC(105)
C ..
DIPMOM = SQRT(G(2)**2+G(3)**2+H(3)**2)
C
P = XGSW**2
U = ZGSW**2
V = 3.D0*ZGSW*XGSW
T = YGSW**2
Q = DIPMOM/SQRT(P+T+U)**5
BXGSW = Q*((T+U-2.D0*P)*SPS-V*CPS)
BYGSW = -3.D0*YGSW*Q*(XGSW*SPS+ZGSW*CPS)
BZGSW = Q*((P+T-2.D0*U)*CPS-V*SPS)
RETURN
END
C *******************************************************************
c
SUBROUTINE TS_SUN_08(IYEAR,IDAY,IHOUR,MIN,ISEC,GST,SLONG,SRASN,
* SDEC)
C *PT*WARNING* Already double-precision
C
C CALCULATES FOUR QUANTITIES NECESSARY FOR COORDINATE TRANSFORMATIONS
CWHICH DEPEND ON TS_SUN POSITION (AND, HENCE, ON UNIVERSAL TIME AND SEASON)
C
C ------- INPUT PARAMETERS:
CIYR,IDAY,IHOUR,MIN,ISEC - YEAR, DAY, AND UNIVERSAL TIME IN HOURS, MINUTES,
C AND SECONDS (IDAY=1 CORRESPONDS TO JANUARY 1).
C
C ------- OUTPUT PARAMETERS:
C GST - GREENWICH MEAN SIDEREAL TIME, SLONG - LONGITUDE ALONG ECLIPTIC
C SRASN - RIGHT ASCENSION, SDEC - DECLINATION OF THE TS_SUN (RADIANS)
C ORIGINAL VERSION OF THIS SUBROUTINE HAS BEEN COMPILED FROM:
C RUSSELL, C.T., COSMIC ELECTRODYNAMICS, 1971, V.2, PP.184-196.
C
C LAST MODIFICATION: MARCH 31, 2003 (ONLY SOME NOTATION CHANGES)
C
C ORIGINAL VERSION WRITTEN BY: Gilbert D. Mead
C
C .. Scalar Arguments ..
DOUBLE PRECISION GST,SDEC,SLONG,SRASN
INTEGER IDAY,IHOUR,ISEC,IYEAR,MIN
C ..
C .. Local Scalars ..
DOUBLE PRECISION COSD,DJ,FDAY,G,OBLIQ,RAD,SC,SIND,SLP,SOB,T,VL
C ..
C .. External Functions ..
C DOUBLE PRECISION DFLOAT
C EXTERNAL DFLOAT
C ..
C .. Intrinsic Functions ..
INTRINSIC ATAN,ATAN2,COS,DMOD,SIN,SQRT
C ..
C .. Data statements ..
DATA RAD/57.295779513D0/
C ..
C
IF (IYEAR.LT.1901 .OR. IYEAR.GT.2099) RETURN
C *PT*WARNING* Constant already double-precision
FDAY = DBLE(IHOUR*3600+MIN*60+ISEC)/86400.D0
C *PT*WARNING* Constant already double-precision
DJ = 365*(IYEAR-1900) + (IYEAR-1901)/4 + IDAY - 0.5D0 + FDAY
T = DJ/36525.D0
C *PT*WARNING* Constant already double-precision
VL = DMOD(279.696678D0+0.9856473354D0*DJ,360.D0)
C *PT*WARNING* Constant already double-precision
GST = DMOD(279.690983D0+.9856473354D0*DJ+360.D0*FDAY+180.D0,
* 360.D0)/RAD
C *PT*WARNING* Constant already double-precision
G = DMOD(358.475845D0+0.985600267D0*DJ,360.D0)/RAD
SLONG = (VL+(1.91946D0-0.004789D0*T)*SIN(G)+
* 0.020094D0*SIN(2.D0*G))/RAD
IF (SLONG.GT.6.2831853D0) SLONG = SLONG - 6.2831853D0
IF (SLONG.LT.0.D0) SLONG = SLONG + 6.2831853D0
OBLIQ = (23.45229D0-0.0130125D0*T)/RAD
SOB = SIN(OBLIQ)
SLP = SLONG - 9.924D-5
C
C THE LAST CONSTANT IS A CORRECTION FOR THE ANGULAR ABERRATION DUE TO
C EARTH'S ORBITAL MOTION
C
SIND = SOB*SIN(SLP)
COSD = SQRT(1.D0-SIND**2)
SC = SIND/COSD
SDEC = ATAN(SC)
SRASN = 3.141592654D0 - ATAN2(COS(OBLIQ)/SOB*SC,-COS(SLP)/COSD)
RETURN
END
C
C================================================================================
c
SUBROUTINE TS_SPHCAR_08(R,THETA,PHI,X,Y,Z,J)
C
C CONVERTS SPHERICAL COORDS INTO CARTESIAN ONES AND VICE VERSA
C (THETA AND PHI IN RADIANS).
C
C J>0 J<0
C -----INPUT: J,R,THETA,PHI J,X,Y,Z
C ----OUTPUT: X,Y,Z R,THETA,PHI
C
C NOTE: AT THE POLES (X=0 AND Y=0) WE ASSUME PHI=0 WHEN CONVERTING
C FROM CARTESIAN TO SPHERICAL COORDS (I.E., FOR J<0)
C
C LAST MOFIFICATION: APRIL 1, 2003 (ONLY SOME NOTATION CHANGES AND MORE
C COMMENTS ADDED)
C
C AUTHOR: N. A. TSYGANENKO
C
C .. Scalar Arguments ..
DOUBLE PRECISION PHI,R,THETA,X,Y,Z
INTEGER J
C ..
C .. Local Scalars ..
DOUBLE PRECISION SQ
C ..
C .. Intrinsic Functions ..
INTRINSIC ATAN2,COS,SIN,SQRT
C ..
IF (J.GT.0) GO TO 30
SQ = X**2 + Y**2
R = SQRT(SQ+Z**2)
IF (SQ.NE.0.D0) GO TO 20
PHI = 0.D0
IF (Z.LT.0.D0) GO TO 10
THETA = 0.D0
RETURN
10 THETA = 3.141592654D0
RETURN
20 SQ = SQRT(SQ)
PHI = ATAN2(Y,X)
THETA = ATAN2(SQ,Z)
IF (PHI.LT.0.D0) PHI = PHI + 6.28318531D0
RETURN
30 SQ = R*SIN(THETA)
X = SQ*COS(PHI)
Y = SQ*SIN(PHI)
Z = R*COS(THETA)
RETURN
END
C
C===========================================================================
c
SUBROUTINE TS_BSPCAR_08(THETA,PHI,BR,BTHETA,BPHI,BX,BY,BZ)
C
C CALCULATES CARTESIAN FIELD COMPONENTS FROM LOCAL SPHERICAL ONES
C
C -----INPUT: THETA,PHI - SPHERICAL ANGLES OF THE POINT IN RADIANS
C BR,BTHETA,BPHI - LOCAL SPHERICAL COMPONENTS OF THE FIELD
C -----OUTPUT: BX,BY,BZ - CARTESIAN COMPONENTS OF THE FIELD
C
C LAST MOFIFICATION: APRIL 1, 2003 (ONLY SOME NOTATION CHANGES)
C
C WRITTEN BY: N. A. TSYGANENKO
C
C .. Scalar Arguments ..
DOUBLE PRECISION BPHI,BR,BTHETA,BX,BY,BZ,PHI,THETA
C ..
C .. Local Scalars ..
DOUBLE PRECISION BE,C,CF,S,SF
C ..
C .. Intrinsic Functions ..
INTRINSIC COS,SIN
C ..
S = SIN(THETA)
C = COS(THETA)
SF = SIN(PHI)
CF = COS(PHI)
BE = BR*S + BTHETA*C
BX = BE*CF - BPHI*SF
BY = BE*SF + BPHI*CF
BZ = BR*C - BTHETA*S
RETURN
END
c
C==============================================================================
C
SUBROUTINE BCARSP_08(X,Y,Z,BX,BY,BZ,BR,BTHETA,BPHI)
C
CALCULATES LOCAL SPHERICAL FIELD COMPONENTS FROM THOSE IN CARTESIAN SYSTEM
C
C -----INPUT: X,Y,Z - CARTESIAN COMPONENTS OF THE POSITION VECTOR
C BX,BY,BZ - CARTESIAN COMPONENTS OF THE FIELD VECTOR
C-----OUTPUT: BR,BTHETA,BPHI - LOCAL SPHERICAL COMPONENTS OF THE FIELD VECTOR
C
C NOTE: AT THE POLES (THETA=0 OR THETA=PI) WE ASSUME PHI=0,
C AND HENCE BTHETA=BX, BPHI=BY
C
C WRITTEN AND ADDED TO THIS PACKAGE: APRIL 1, 2003,
C AUTHOR: N. A. TSYGANENKO
C
C .. Scalar Arguments ..
DOUBLE PRECISION BPHI,BR,BTHETA,BX,BY,BZ,X,Y,Z
C ..
C .. Local Scalars ..
DOUBLE PRECISION CPHI,CT,R,RHO,RHO2,SPHI,ST
C ..
C .. Intrinsic Functions ..
INTRINSIC SQRT
C ..
RHO2 = X**2 + Y**2
R = SQRT(RHO2+Z**2)
RHO = SQRT(RHO2)
IF (RHO.NE.0.D0) THEN
CPHI = X/RHO
SPHI = Y/RHO
ELSE
CPHI = 1.D0
SPHI = 0.D0
END IF
CT = Z/R
ST = RHO/R
BR = (X*BX+Y*BY+Z*BZ)/R
BTHETA = (BX*CPHI+BY*SPHI)*CT - BZ*ST
BPHI = BY*CPHI - BX*SPHI
RETURN
END
C
c=====================================================================================
C
SUBROUTINE TS_RECALC_08(IYEAR,IDAY,IHOUR,MIN,ISEC,VGSEX,VGSEY,
* VGSEZ)
C
C1. PREPARES ELEMENTS OF ROTATION MATRICES FOR TRANSFORMATIONS OF VECTORS BETWEEN
C SEVERAL COORDINATE SYSTEMS, MOST FREQUENTLY USED IN SPACE PHYSICS.
C
C2. PREPARES COEFFICIENTS USED IN THE CALCULATION OF THE MAIN TS_GEOMAGNETIC FIELD
C (IGRF MODEL)
C
CTHIS SUBROUTINE SHOULD BE INVOKED BEFORE USING THE FOLLOWING SUBROUTINES:
CIGRF_GEO_08, IGRF_GSW_08, DIP_08, TS_GEOMAG_08, GEOGSW, MAGSW_08, SMGSW_08, GSWGSE,
c GEIGEO_08, TRACE, STEP_08, RHAND_08.
C
CTHERE IS NO NEED TO REPEATEDLY INVOKE TS_RECALC_08, IF MULTIPLE CALCULATIONS ARE MADE
C FOR THE SAME DATE/TIME AND SOLAR WIND FLOW DIRECTION.
C
C -----INPUT PARAMETERS:
C
C IYEAR - YEAR NUMBER (FOUR DIGITS)
C IDAY - DAY OF YEAR (DAY 1 = JAN 1)
C IHOUR - HOUR OF DAY (00 TO 23)
C MIN - MINUTE OF HOUR (00 TO 59)
C ISEC - SECONDS OF MINUTE (00 TO 59)
CVGSEX,VGSEY,VGSEZ - GSE (GEOCENTRIC SOLAR-ECLIPTIC) COMPONENTS OF THE OBSERVED
C SOLAR WIND FLOW VELOCITY (IN KM/S)
C
CIMPORTANT: IF ONLY QUESTIONABLE INFORMATION (OR NO INFORMATION AT ALL) IS AVAILABLE
C ON THE SOLAR WIND SPEED, OR, IF THE STANDARD GSM AND/OR SM COORDINATES ARE
C INTENDED TO BE USED, THEN SET VGSEX=-400.0 AND VGSEY=VGSEZ=0. IN THIS CASE,
C THE GSW COORDINATE SYSTEM BECOMES IDENTICAL TO THE STANDARD GSM.
C
C IF ONLY SCALAR SPEED V OF THE SOLAR WIND IS KNOWN, THEN SETTING
C VGSEX=-V, VGSEY=29.78, VGSEZ=0.0 WILL TAKE INTO ACCOUNT THE ~4 degs
C ABERRATION OF THE MAGNETOSPHERE DUE TO EARTH'S ORBITAL MOTION
C
C IF ALL THREE GSE COMPONENTS OF THE SOLAR WIND VELOCITY ARE AVAILABLE,
C PLEASE NOTE THAT IN SOME SOLAR WIND DATABASES THE ABERRATION EFFECT
C HAS ALREADY BEEN TAKEN INTO ACCOUNT BY SUBTRACTING 29.78 KM/S FROM VYGSE;
C IN THAT CASE, THE UNABERRATED (OBSERVED) VYGSE VALUES SHOULD BE RESTORED
C BY ADDING BACK THE 29.78 KM/S CORRECTION. WHETHER OR NOT TO DO THAT, MUST
C BE EITHER VERIFIED WITH THE DATA ORIGINATOR OR DETERMINED BY AVERAGING
C VGSEY OVER A SUFFICIENTLY LONG TIME INTERVAL.
C
C -----OUTPUT PARAMETERS: NONE (ALL OUTPUT QUANTITIES ARE PLACED
C INTO THE COMMON BLOCKS /GEOPACK1/ AND /GEOPACK2/)
C
C OTHER SUBROUTINES CALLED BY THIS ONE: TS_SUN_08
C
C AUTHOR: N.A. TSYGANENKO
C DATE: DEC.1, 1991
C
C REVISION OF NOVEMBER 30, 2010:
C
CThe table of IGRF coefficients was extended to include those for the epoch 2010
c (for details, see http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html)
C
CREVISION OF NOVEMBER 15, 2007: ADDED THE POSSIBILITY TO TAKE INTO ACCOUNT THE OBSERVED
CDEFLECTION OF THE SOLAR WIND FLOW FROM STRICTLY RADIAL DIRECTION. TO THAT END, THREE
CGSE COMPONENTS OF THE SOLAR WIND VELOCITY WERE ADDED TO THE INPUT PARAMETERS.
C
c CORRECTION OF MAY 9, 2006: INTERPOLATION OF THE COEFFICIENTS (BETWEEN
C LABELS 50 AND 105) IS NOW MADE THROUGH THE LAST ELEMENT OF THE ARRAYS
C G(105) AND H(105) (PREVIOUSLY MADE ONLY THROUGH N=66, WHICH IN SOME
C CASES CAUSED RUNTIME ERRORS)
c
C REVISION OF MAY 3, 2005:
CThe table of IGRF coefficients was extended to include those for the epoch 2005
c the maximal order of spherical harmonics was also increased up to 13
c (for details, see http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html)
c
C REVISION OF APRIL 3, 2003:
cThe code now includes preparation of the model coefficients for the subroutines
cIGRF_08 and TS_GEOMAG_08. This eliminates the need for the SAVE statements, used
cin the old versions, making the codes easier and more compiler-independent.
C
C
C
CTHE COMMON BLOCK /GEOPACK1/ CONTAINS ELEMENTS OF THE ROTATION MATRICES AND OTHER
CPARAMETERS RELATED TO THE COORDINATE TRANSFORMATIONS PERFORMED BY THIS PACKAGE
C
C
CTHE COMMON BLOCK /GEOPACK2/ CONTAINS COEFFICIENTS OF THE IGRF FIELD MODEL, CALCULATED
CFOR A GIVEN YEAR AND DAY FROM THEIR STANDARD EPOCH VALUES. THE ARRAY REC CONTAINS
CCOEFFICIENTS USED IN THE RECURSION RELATIONS FOR LEGENDRE ASSOCIATE POLYNOMIALS.
C
C .. Scalar Arguments ..
DOUBLE PRECISION VGSEX,VGSEY,VGSEZ
INTEGER IDAY,IHOUR,ISEC,IYEAR,MIN
C ..
C .. Local Scalars ..
DOUBLE PRECISION AA,DIP1,DIP2,DIP3,DJ,DT,DX1,DX2,DX3,DY1,DY2,DY3,
* DZ1,DZ2,DZ3,EXMAGX,EXMAGY,EXMAGZ,EYMAGX,EYMAGY,
* F1,F2,GST,G_10,G_11,H_11,OBLIQ,P,S,S1,S2,S3,SDEC,
* SLONG,SQ,SQQ,SQR,SRASN,T,V,X1,X2,X3,Y,Y1,Y2,Y3,
* Z1,Z2,Z3
INTEGER ISW,IY,M,MN,MNN,N,N2
C ..
C .. Local Arrays ..
DOUBLE PRECISION DG10(45),DH10(45),G00(105),G05(105),G10(105),
* G65(105),G70(105),G75(105),G80(105),G85(105),
* G90(105),G95(105),H00(105),H05(105),H10(105),
* H65(105),H70(105),H75(105),H80(105),H85(105),
* H90(105),H95(105)
C ..
C .. External Subroutines ..
EXTERNAL TS_SUN_08
C ..
C .. Intrinsic Functions ..
INTRINSIC ASIN,COS,DBLE,SIN,SQRT
C ..
C .. Common blocks ..
COMMON /GEOPACK1/ST0,CT0,SL0,CL0,CTCL,STCL,CTSL,STSL,SFI,CFI,SPS,
* CPS,DS3,CGST,SGST,PSI,A11,A21,A31,A12,A22,A32,A13,A23,A33,
* E11,E21,E31,E12,E22,E32,E13,E23,E33
COMMON /GEOPACK2/G,H,REC
DOUBLE PRECISION A11,A12,A13,A21,A22,A23,A31,A32,A33,CFI,CGST,CL0,
* CPS,CT0,CTCL,CTSL,DS3,E11,E12,E13,E21,E22,E23,
* E31,E32,E33,PSI,SFI,SGST,SL0,SPS,ST0,STCL,STSL
DOUBLE PRECISION G(105),H(105),REC(105)
C added by B. Rideout to allow time to be passed to my igrf call
COMMON /BR_FYEAR/FYEAR
DOUBLE PRECISION FYEAR
C ..
C .. Save statement ..
SAVE ISW
C ..
C .. Data statements ..
C
c
c
c
c
c
c
C
C
C
C
C
C
C
C
C
C
C
DATA ISW/0/
DATA G65/0.D0,-30334.D0,-2119.D0,-1662.D0,2997.D0,1594.D0,1297.D0,
* -2038.D0,1292.D0,856.D0,957.D0,804.D0,479.D0,-390.D0,252.D0,
* -219.D0,358.D0,254.D0,-31.D0,-157.D0,-62.D0,45.D0,61.D0,8.D0,
* -228.D0,4.D0,1.D0,-111.D0,75.D0,-57.D0,4.D0,13.D0,-26.D0,
* -6.D0,13.D0,1.D0,13.D0,5.D0,-4.D0,-14.D0,0.D0,8.D0,-1.D0,
* 11.D0,4.D0,8.D0,10.D0,2.D0,-13.D0,10.D0,-1.D0,-1.D0,5.D0,
* 1.D0,-2.D0,-2.D0,-3.D0,2.D0,-5.D0,-2.D0,4.D0,4.D0,0.D0,2.D0,
* 2.D0,0.D0,39*0.D0/
DATA H65/0.D0,0.D0,5776.D0,0.D0,-2016.D0,114.D0,0.D0,-404.D0,
* 240.D0,-165.D0,0.D0,148.D0,-269.D0,13.D0,-269.D0,0.D0,19.D0,
* 128.D0,-126.D0,-97.D0,81.D0,0.D0,-11.D0,100.D0,68.D0,-32.D0,
* -8.D0,-7.D0,0.D0,-61.D0,-27.D0,-2.D0,6.D0,26.D0,-23.D0,
* -12.D0,0.D0,7.D0,-12.D0,9.D0,-16.D0,4.D0,24.D0,-3.D0,-17.D0,
* 0.D0,-22.D0,15.D0,7.D0,-4.D0,-5.D0,10.D0,10.D0,-4.D0,1.D0,
* 0.D0,2.D0,1.D0,2.D0,6.D0,-4.D0,0.D0,-2.D0,3.D0,0.D0,-6.D0,
* 39*0.D0/
DATA G70/0.D0,-30220.D0,-2068.D0,-1781.D0,3000.D0,1611.D0,1287.D0,
* -2091.D0,1278.D0,838.D0,952.D0,800.D0,461.D0,-395.D0,234.D0,
* -216.D0,359.D0,262.D0,-42.D0,-160.D0,-56.D0,43.D0,64.D0,
* 15.D0,-212.D0,2.D0,3.D0,-112.D0,72.D0,-57.D0,1.D0,14.D0,
* -22.D0,-2.D0,13.D0,-2.D0,14.D0,6.D0,-2.D0,-13.D0,-3.D0,5.D0,
* 0.D0,11.D0,3.D0,8.D0,10.D0,2.D0,-12.D0,10.D0,-1.D0,0.D0,3.D0,
* 1.D0,-1.D0,-3.D0,-3.D0,2.D0,-5.D0,-1.D0,6.D0,4.D0,1.D0,0.D0,
* 3.D0,-1.D0,39*0.D0/
DATA H70/0.D0,0.D0,5737.D0,0.D0,-2047.D0,25.D0,0.D0,-366.D0,
* 251.D0,-196.D0,0.D0,167.D0,-266.D0,26.D0,-279.D0,0.D0,26.D0,
* 139.D0,-139.D0,-91.D0,83.D0,0.D0,-12.D0,100.D0,72.D0,-37.D0,
* -6.D0,1.D0,0.D0,-70.D0,-27.D0,-4.D0,8.D0,23.D0,-23.D0,-11.D0,
* 0.D0,7.D0,-15.D0,6.D0,-17.D0,6.D0,21.D0,-6.D0,-16.D0,0.D0,
* -21.D0,16.D0,6.D0,-4.D0,-5.D0,10.D0,11.D0,-2.D0,1.D0,0.D0,
* 1.D0,1.D0,3.D0,4.D0,-4.D0,0.D0,-1.D0,3.D0,1.D0,-4.D0,39*0.D0/
DATA G75/0.D0,-30100.D0,-2013.D0,-1902.D0,3010.D0,1632.D0,1276.D0,
* -2144.D0,1260.D0,830.D0,946.D0,791.D0,438.D0,-405.D0,216.D0,
* -218.D0,356.D0,264.D0,-59.D0,-159.D0,-49.D0,45.D0,66.D0,
* 28.D0,-198.D0,1.D0,6.D0,-111.D0,71.D0,-56.D0,1.D0,16.D0,
* -14.D0,0.D0,12.D0,-5.D0,14.D0,6.D0,-1.D0,-12.D0,-8.D0,4.D0,
* 0.D0,10.D0,1.D0,7.D0,10.D0,2.D0,-12.D0,10.D0,-1.D0,-1.D0,
* 4.D0,1.D0,-2.D0,-3.D0,-3.D0,2.D0,-5.D0,-2.D0,5.D0,4.D0,1.D0,
* 0.D0,3.D0,-1.D0,39*0.D0/
DATA H75/0.D0,0.D0,5675.D0,0.D0,-2067.D0,-68.D0,0.D0,-333.D0,
* 262.D0,-223.D0,0.D0,191.D0,-265.D0,39.D0,-288.D0,0.D0,31.D0,
* 148.D0,-152.D0,-83.D0,88.D0,0.D0,-13.D0,99.D0,75.D0,-41.D0,
* -4.D0,11.D0,0.D0,-77.D0,-26.D0,-5.D0,10.D0,22.D0,-23.D0,
* -12.D0,0.D0,6.D0,-16.D0,4.D0,-19.D0,6.D0,18.D0,-10.D0,-17.D0,
* 0.D0,-21.D0,16.D0,7.D0,-4.D0,-5.D0,10.D0,11.D0,-3.D0,1.D0,
* 0.D0,1.D0,1.D0,3.D0,4.D0,-4.D0,-1.D0,-1.D0,3.D0,1.D0,-5.D0,
* 39*0.D0/
DATA G80/0.D0,-29992.D0,-1956.D0,-1997.D0,3027.D0,1663.D0,1281.D0,
* -2180.D0,1251.D0,833.D0,938.D0,782.D0,398.D0,-419.D0,199.D0,
* -218.D0,357.D0,261.D0,-74.D0,-162.D0,-48.D0,48.D0,66.D0,
* 42.D0,-192.D0,4.D0,14.D0,-108.D0,72.D0,-59.D0,2.D0,21.D0,
* -12.D0,1.D0,11.D0,-2.D0,18.D0,6.D0,0.D0,-11.D0,-7.D0,4.D0,
* 3.D0,6.D0,-1.D0,5.D0,10.D0,1.D0,-12.D0,9.D0,-3.D0,-1.D0,7.D0,
* 2.D0,-5.D0,-4.D0,-4.D0,2.D0,-5.D0,-2.D0,5.D0,3.D0,1.D0,2.D0,
* 3.D0,0.D0,39*0.D0/
DATA H80/0.D0,0.D0,5604.D0,0.D0,-2129.D0,-200.D0,0.D0,-336.D0,
* 271.D0,-252.D0,0.D0,212.D0,-257.D0,53.D0,-297.D0,0.D0,46.D0,
* 150.D0,-151.D0,-78.D0,92.D0,0.D0,-15.D0,93.D0,71.D0,-43.D0,
* -2.D0,17.D0,0.D0,-82.D0,-27.D0,-5.D0,16.D0,18.D0,-23.D0,
* -10.D0,0.D0,7.D0,-18.D0,4.D0,-22.D0,9.D0,16.D0,-13.D0,-15.D0,
* 0.D0,-21.D0,16.D0,9.D0,-5.D0,-6.D0,9.D0,10.D0,-6.D0,2.D0,
* 0.D0,1.D0,0.D0,3.D0,6.D0,-4.D0,0.D0,-1.D0,4.D0,0.D0,-6.D0,
* 39*0.D0/
DATA G85/0.D0,-29873.D0,-1905.D0,-2072.D0,3044.D0,1687.D0,1296.D0,
* -2208.D0,1247.D0,829.D0,936.D0,780.D0,361.D0,-424.D0,170.D0,
* -214.D0,355.D0,253.D0,-93.D0,-164.D0,-46.D0,53.D0,65.D0,
* 51.D0,-185.D0,4.D0,16.D0,-102.D0,74.D0,-62.D0,3.D0,24.D0,
* -6.D0,4.D0,10.D0,0.D0,21.D0,6.D0,0.D0,-11.D0,-9.D0,4.D0,4.D0,
* 4.D0,-4.D0,5.D0,10.D0,1.D0,-12.D0,9.D0,-3.D0,-1.D0,7.D0,1.D0,
* -5.D0,-4.D0,-4.D0,3.D0,-5.D0,-2.D0,5.D0,3.D0,1.D0,2.D0,3.D0,
* 0.D0,39*0.D0/
DATA H85/0.D0,0.D0,5500.D0,0.D0,-2197.D0,-306.D0,0.D0,-310.D0,
* 284.D0,-297.D0,0.D0,232.D0,-249.D0,69.D0,-297.D0,0.D0,47.D0,
* 150.D0,-154.D0,-75.D0,95.D0,0.D0,-16.D0,88.D0,69.D0,-48.D0,
* -1.D0,21.D0,0.D0,-83.D0,-27.D0,-2.D0,20.D0,17.D0,-23.D0,
* -7.D0,0.D0,8.D0,-19.D0,5.D0,-23.D0,11.D0,14.D0,-15.D0,-11.D0,
* 0.D0,-21.D0,15.D0,9.D0,-6.D0,-6.D0,9.D0,9.D0,-7.D0,2.D0,0.D0,
* 1.D0,0.D0,3.D0,6.D0,-4.D0,0.D0,-1.D0,4.D0,0.D0,-6.D0,39*0.D0/
DATA G90/0.D0,-29775.D0,-1848.D0,-2131.D0,3059.D0,1686.D0,1314.D0,
* -2239.D0,1248.D0,802.D0,939.D0,780.D0,325.D0,-423.D0,141.D0,
* -214.D0,353.D0,245.D0,-109.D0,-165.D0,-36.D0,61.D0,65.D0,
* 59.D0,-178.D0,3.D0,18.D0,-96.D0,77.D0,-64.D0,2.D0,26.D0,
* -1.D0,5.D0,9.D0,0.D0,23.D0,5.D0,-1.D0,-10.D0,-12.D0,3.D0,
* 4.D0,2.D0,-6.D0,4.D0,9.D0,1.D0,-12.D0,9.D0,-4.D0,-2.D0,7.D0,
* 1.D0,-6.D0,-3.D0,-4.D0,2.D0,-5.D0,-2.D0,4.D0,3.D0,1.D0,3.D0,
* 3.D0,0.D0,39*0.D0/
DATA H90/0.D0,0.D0,5406.D0,0.D0,-2279.D0,-373.D0,0.D0,-284.D0,
* 293.D0,-352.D0,0.D0,247.D0,-240.D0,84.D0,-299.D0,0.D0,46.D0,
* 154.D0,-153.D0,-69.D0,97.D0,0.D0,-16.D0,82.D0,69.D0,-52.D0,
* 1.D0,24.D0,0.D0,-80.D0,-26.D0,0.D0,21.D0,17.D0,-23.D0,-4.D0,
* 0.D0,10.D0,-19.D0,6.D0,-22.D0,12.D0,12.D0,-16.D0,-10.D0,0.D0,
* -20.D0,15.D0,11.D0,-7.D0,-7.D0,9.D0,8.D0,-7.D0,2.D0,0.D0,
* 2.D0,1.D0,3.D0,6.D0,-4.D0,0.D0,-2.D0,3.D0,-1.D0,-6.D0,
* 39*0.D0/
DATA G95/0.D0,-29692.D0,-1784.D0,-2200.D0,3070.D0,1681.D0,1335.D0,
* -2267.D0,1249.D0,759.D0,940.D0,780.D0,290.D0,-418.D0,122.D0,
* -214.D0,352.D0,235.D0,-118.D0,-166.D0,-17.D0,68.D0,67.D0,
* 68.D0,-170.D0,-1.D0,19.D0,-93.D0,77.D0,-72.D0,1.D0,28.D0,
* 5.D0,4.D0,8.D0,-2.D0,25.D0,6.D0,-6.D0,-9.D0,-14.D0,9.D0,6.D0,
* -5.D0,-7.D0,4.D0,9.D0,3.D0,-10.D0,8.D0,-8.D0,-1.D0,10.D0,
* -2.D0,-8.D0,-3.D0,-6.D0,2.D0,-4.D0,-1.D0,4.D0,2.D0,2.D0,5.D0,
* 1.D0,0.D0,39*0.D0/
DATA H95/0.D0,0.D0,5306.D0,0.D0,-2366.D0,-413.D0,0.D0,-262.D0,
* 302.D0,-427.D0,0.D0,262.D0,-236.D0,97.D0,-306.D0,0.D0,46.D0,
* 165.D0,-143.D0,-55.D0,107.D0,0.D0,-17.D0,72.D0,67.D0,-58.D0,
* 1.D0,36.D0,0.D0,-69.D0,-25.D0,4.D0,24.D0,17.D0,-24.D0,-6.D0,
* 0.D0,11.D0,-21.D0,8.D0,-23.D0,15.D0,11.D0,-16.D0,-4.D0,0.D0,
* -20.D0,15.D0,12.D0,-6.D0,-8.D0,8.D0,5.D0,-8.D0,3.D0,0.D0,
* 1.D0,0.D0,4.D0,5.D0,-5.D0,-1.D0,-2.D0,1.D0,-2.D0,-7.D0,
* 39*0.D0/
DATA G00/0.D0,-29619.4D0,-1728.2D0,-2267.7D0,3068.4D0,1670.9D0,
* 1339.6D0,-2288.D0,1252.1D0,714.5D0,932.3D0,786.8D0,250.D0,
* -403.D0,111.3D0,-218.8D0,351.4D0,222.3D0,-130.4D0,-168.6D0,
* -12.9D0,72.3D0,68.2D0,74.2D0,-160.9D0,-5.9D0,16.9D0,-90.4D0,
* 79.0D0,-74.0D0,0.D0,33.3D0,9.1D0,6.9D0,7.3D0,-1.2D0,24.4D0,
* 6.6D0,-9.2D0,-7.9D0,-16.6D0,9.1D0,7.0D0,-7.9D0,-7.D0,5.D0,
* 9.4D0,3.D0,-8.4D0,6.3D0,-8.9D0,-1.5D0,9.3D0,-4.3D0,-8.2D0,
* -2.6D0,-6.D0,1.7D0,-3.1D0,-0.5D0,3.7D0,1.D0,2.D0,4.2D0,0.3D0,
* -1.1D0,2.7D0,-1.7D0,-1.9D0,1.5D0,-0.1D0,0.1D0,-0.7D0,0.7D0,
* 1.7D0,0.1D0,1.2D0,4.0D0,-2.2D0,-0.3D0,0.2D0,0.9D0,-0.2D0,
* 0.9D0,-0.5D0,0.3D0,-0.3D0,-0.4D0,-0.1D0,-0.2D0,-0.4D0,-0.2D0,
* -0.9D0,0.3D0,0.1D0,-0.4D0,1.3D0,-0.4D0,0.7D0,-0.4D0,0.3D0,
* -0.1D0,0.4D0,0.D0,0.1D0/
DATA H00/0.D0,0.D0,5186.1D0,0.D0,-2481.6D0,-458.0D0,0.D0,-227.6D0,
* 293.4D0,-491.1D0,0.D0,272.6D0,-231.9D0,119.8D0,-303.8D0,0.D0,
* 43.8D0,171.9D0,-133.1D0,-39.3D0,106.3D0,0.D0,-17.4D0,63.7D0,
* 65.1D0,-61.2D0,0.7D0,43.8D0,0.D0,-64.6D0,-24.2D0,6.2D0,24.D0,
* 14.8D0,-25.4D0,-5.8D0,0.0D0,11.9D0,-21.5D0,8.5D0,-21.5D0,
* 15.5D0,8.9D0,-14.9D0,-2.1D0,0.0D0,-19.7D0,13.4D0,12.5D0,
* -6.2D0,-8.4D0,8.4D0,3.8D0,-8.2D0,4.8D0,0.0D0,1.7D0,0.0D0,
* 4.0D0,4.9D0,-5.9D0,-1.2D0,-2.9D0,0.2D0,-2.2D0,-7.4D0,0.0D0,
* 0.1D0,1.3D0,-0.9D0,-2.6D0,0.9D0,-0.7D0,-2.8D0,-0.9D0,-1.2D0,
* -1.9D0,-0.9D0,0.0D0,-0.4D0,0.3D0,2.5D0,-2.6D0,0.7D0,0.3D0,
* 0.0D0,0.0D0,0.3D0,-0.9D0,-0.4D0,0.8D0,0.0D0,-0.9D0,0.2D0,
* 1.8D0,-0.4D0,-1.0D0,-0.1D0,0.7D0,0.3D0,0.6D0,0.3D0,-0.2D0,
* -0.5D0,-0.9D0/
DATA G05/0.D0,-29554.6D0,-1669.0D0,-2337.2D0,3047.7D0,1657.8D0,
* 1336.3D0,-2305.8D0,1246.4D0,672.5D0,920.6D0,798.0D0,210.7D0,
* -379.9D0,100.0D0,-227.0D0,354.4D0,208.9D0,-136.5D0,-168.1D0,
* -13.6D0,73.6D0,69.6D0,76.7D0,-151.3D0,-14.6D0,14.6D0,-86.4D0,
* 79.9D0,-74.5D0,-1.7D0,38.7D0,12.3D0,9.4D0,5.4D0,1.9D0,24.8D0,
* 7.6D0,-11.7D0,-6.9D0,-18.1D0,10.2D0,9.4D0,-11.3D0,-4.9D0,
* 5.6D0,9.8D0,3.6D0,-6.9D0,5.0D0,-10.8D0,-1.3D0,8.8D0,-6.7D0,
* -9.2D0,-2.2D0,-6.1D0,1.4D0,-2.4D0,-0.2D0,3.1D0,0.3D0,2.1D0,
* 3.8D0,-0.2D0,-2.1D0,2.9D0,-1.6D0,-1.9D0,1.4D0,-0.3D0,0.3D0,
* -0.8D0,0.5D0,1.8D0,0.2D0,1.0D0,4.0D0,-2.2D0,-0.3D0,0.2D0,
* 0.9D0,-0.4D0,1.0D0,-0.3D0,0.5D0,-0.4D0,-0.4D0,0.1D0,-0.5D0,
* -0.1D0,-0.2D0,-0.9D0,0.3D0,0.3D0,-0.4D0,1.2D0,-0.4D0,0.8D0,
* -0.3D0,0.4D0,-0.1D0,0.4D0,-0.1D0,-0.2D0/
DATA H05/0.D0,0.0D0,5078.0D0,0.0D0,-2594.5D0,-515.4D0,0.0D0,
* -198.9D0,269.7D0,-524.7D0,0.0D0,282.1D0,-225.2D0,145.2D0,
* -305.4D0,0.0D0,42.7D0,180.3D0,-123.5D0,-19.6D0,103.9D0,0.0D0,
* -20.3D0,54.8D0,63.6D0,-63.5D0,0.2D0,50.9D0,0.0D0,-61.1D0,
* -22.6D0,6.8D0,25.4D0,10.9D0,-26.3D0,-4.6D0,0.0D0,11.2D0,
* -20.9D0,9.8D0,-19.7D0,16.2D0,7.6D0,-12.8D0,-0.1D0,0.0D0,
* -20.1D0,12.7D0,12.7D0,-6.7D0,-8.2D0,8.1D0,2.9D0,-7.7D0,6.0D0,
* 0.0D0,2.2D0,0.1D0,4.5D0,4.8D0,-6.7D0,-1.0D0,-3.5D0,-0.9D0,
* -2.3D0,-7.9D0,0.0D0,0.3D0,1.4D0,-0.8D0,-2.3D0,0.9D0,-0.6D0,
* -2.7D0,-1.1D0,-1.6D0,-1.9D0,-1.4D0,0.0D0,-0.6D0,0.2D0,2.4D0,
* -2.6D0,0.6D0,0.4D0,0.0D0,0.0D0,0.3D0,-0.9D0,-0.3D0,0.9D0,
* 0.0D0,-0.8D0,0.3D0,1.7D0,-0.5D0,-1.1D0,0.0D0,0.6D0,0.2D0,
* 0.5D0,0.4D0,-0.2D0,-0.6D0,-0.9D0/
DATA G10/0.D0,-29496.5D0,-1585.9D0,-2396.6D0,3026.0D0,1668.6D0,
* 1339.7D0,-2326.3D0,1231.7D0,634.2D0,912.6D0,809.0D0,166.6D0,
* -357.1D0,89.7D0,-231.1D0,357.2D0,200.3D0,-141.2D0,-163.1D0,
* -7.7D0,72.8D0,68.6D0,76.0D0,-141.4D0,-22.9D0,13.1D0,-77.9D0,
* 80.4D0,-75.0D0,-4.7D0,45.3D0,14.0D0,10.4D0,1.6D0,4.9D0,
* 24.3D0,8.2D0,-14.5D0,-5.7D0,-19.3D0,11.6D0,10.9D0,-14.1D0,
* -3.7D0,5.4D0,9.4D0,3.4D0,-5.3D0,3.1D0,-12.4D0,-0.8D0,8.4D0,
* -8.4D0,-10.1D0,-2.0D0,-6.3D0,0.9D0,-1.1D0,-0.2D0,2.5D0,
* -0.3D0,2.2D0,3.1D0,-1.0D0,-2.8D0,3.0D0,-1.5D0,-2.1D0,1.6D0,
* -0.5D0,0.5D0,-0.8D0,0.4D0,1.8D0,0.2D0,0.8D0,3.8D0,-2.1D0,
* -0.2D0,0.3D0,1.0D0,-0.7D0,0.9D0,-0.1D0,0.5D0,-0.4D0,-0.4D0,
* 0.2D0,-0.8D0,0.0D0,-0.2D0,-0.9D0,0.3D0,0.4D0,-0.4D0,1.1D0,
* -0.3D0,0.8D0,-0.2D0,0.4D0,0.0D0,0.4D0,-0.3D0,-0.3D0/
DATA H10/0.0D0,0.0D0,4945.1D0,0.0D0,-2707.7D0,-575.4D0,0.0D0,
* -160.5D0,251.7D0,-536.8D0,0.0D0,286.4D0,-211.2D0,164.4D0,
* -309.2D0,0.0D0,44.7D0,188.9D0,-118.1D0,0.1D0,100.9D0,0.0D0,
* -20.8D0,44.2D0,61.5D0,-66.3D0,3.1D0,54.9D0,0.0D0,-57.8D0,
* -21.2D0,6.6D0,24.9D0,7.0D0,-27.7D0,-3.4D0,0.0D0,10.9D0,
* -20.0D0,11.9D0,-17.4D0,16.7D0,7.1D0,-10.8D0,1.7D0,0.0D0,
* -20.5D0,11.6D0,12.8D0,-7.2D0,-7.4D0,8.0D0,2.2D0,-6.1D0,7.0D0,
* 0.0D0,2.8D0,-0.1D0,4.7D0,4.4D0,-7.2D0,-1.0D0,-4.0D0,-2.0D0,
* -2.0D0,-8.3D0,0.0D0,0.1D0,1.7D0,-0.6D0,-1.8D0,0.9D0,-0.4D0,
* -2.5D0,-1.3D0,-2.1D0,-1.9D0,-1.8D0,0.0D0,-0.8D0,0.3D0,2.2D0,
* -2.5D0,0.5D0,0.6D0,0.0D0,0.1D0,0.3D0,-0.9D0,-0.2D0,0.8D0,
* 0.0D0,-0.8D0,0.3D0,1.7D0,-0.6D0,-1.2D0,-0.1D0,0.5D0,0.1D0,
* 0.5D0,0.4D0,-0.2D0,-0.5D0,-0.8D0/
DATA DG10/0.0D0,11.4D0,16.7D0,-11.3D0,-3.9D0,2.7D0,1.3D0,-3.9D0,
* -2.9D0,-8.1D0,-1.4D0,2.0D0,-8.9D0,4.4D0,-2.3D0,-0.5D0,0.5D0,
* -1.5D0,-0.7D0,1.3D0,1.4D0,-0.3D0,-0.3D0,-0.3D0,1.9D0,-1.6D0,
* -0.2D0,1.8D0,0.2D0,-0.1D0,-0.6D0,1.4D0,0.3D0,0.1D0,-0.8D0,
* 0.4D0,-0.1D0,0.1D0,-0.5D0,0.3D0,-0.3D0,0.3D0,0.2D0,-0.5D0,
* 0.2D0/
DATA DH10/0.0D0,0.0D0,-28.8D0,0.0D0,-23.0D0,-12.9D0,0.0D0,8.6D0,
* -2.9D0,-2.1D0,0.0D0,0.4D0,3.2D0,3.6D0,-0.8D0,0.0D0,0.5D0,
* 1.5D0,0.9D0,3.7D0,-0.6D0,0.0D0,-0.1D0,-2.1D0,-0.4D0,-0.5D0,
* 0.8D0,0.5D0,0.0D0,0.6D0,0.3D0,-0.2D0,-0.1D0,-0.8D0,-0.3D0,
* 0.2D0,0.0D0,0.0D0,0.2D0,0.5D0,0.4D0,0.1D0,-0.1D0,0.4D0,0.4D0/
C ..
C
IF (VGSEY.EQ.0.D0 .AND. VGSEZ.EQ.0.D0 .AND. ISW.NE.1) THEN
C PRINT *, ''
C PRINT *,
C *' TS_RECALC_08: RADIAL SOLAR WIND --> GSW SYSTEM IDENTICAL HERE'
C PRINT *,
C *' TO STANDARD GSM (I.E., XGSW AXIS COINCIDES WITH EARTH-TS_SUN LINE)'
C PRINT *, ''
ISW = 1
END IF
IF ((VGSEY.NE.0.D0.OR.VGSEZ.NE.0.D0) .AND. ISW.NE.2) THEN
C PRINT *, ''
C PRINT *,
C *' WARNING: NON-RADIAL SOLAR WIND FLOW SPECIFIED IN TS_RECALC_08;'
C PRINT *,
C *' HENCE XGSW AXIS IS ASSUMED ORIENTED ANTIPARALLEL TO V_SW VECTOR'
C PRINT *, ''
ISW = 2
END IF
C
IY = IYEAR
FYEAR = DBLE(IYEAR) + (IDAY-1)/366.0
C
CWE ARE RESTRICTED BY THE INTERVAL 1965-2010, FOR WHICH THE IGRF COEFFICIENTS
cARE KNOWN; IF IYEAR IS OUTSIDE THIS INTERVAL, THEN THE SUBROUTINE USES THE
C NEAREST LIMITING VALUE AND PRINTS A WARNING:
C
IF (IY.LT.1965) THEN
IY = 1965
C WRITE (*,FMT=9000) IYEAR,IY
END IF
IF (IY.GT.2015) THEN
IY = 2015
C WRITE (*,FMT=9000) IYEAR,IY
END IF
C
CCALCULATE THE ARRAY REC, CONTAINING COEFFICIENTS FOR THE RECURSION RELATIONS,
CUSED IN THE IGRF SUBROUTINE FOR CALCULATING THE ASSOCIATE LEGENDRE POLYNOMIALS
C AND THEIR DERIVATIVES:
c
DO 20 N = 1,14
N2 = 2*N - 1
N2 = N2*(N2-2)
DO 10 M = 1,N
MN = N*(N-1)/2 + M
REC(MN) = DBLE((N-M)*(N+M-2))/DBLE(N2)
10 CONTINUE
20 CONTINUE
C
IF (IY.LT.1970) GO TO 40
IF (IY.LT.1975) GO TO 60
IF (IY.LT.1980) GO TO 80
IF (IY.LT.1985) GO TO 100
IF (IY.LT.1990) GO TO 120
IF (IY.LT.1995) GO TO 140
IF (IY.LT.2000) GO TO 160
IF (IY.LT.2005) GO TO 180
IF (IY.LT.2010) GO TO 200
C
C EXTRAPOLATE BEYOND 2010:
C
DT = DBLE(IY) + DBLE(IDAY-1)/365.25D0 - 2010.D0
DO 30 N = 1,105
G(N) = G10(N)
H(N) = H10(N)
IF (N.GT.45) GO TO 30
G(N) = G(N) + DG10(N)*DT
H(N) = H(N) + DH10(N)*DT
30 CONTINUE
GO TO 220
C
C INTERPOLATE BETWEEEN 1965 - 1970:
C
40 F2 = (DBLE(IY)+DBLE(IDAY-1)/365.25D0-1965)/5.D0
F1 = 1.D0 - F2
DO 50 N = 1,105
G(N) = G65(N)*F1 + G70(N)*F2
H(N) = H65(N)*F1 + H70(N)*F2
50 CONTINUE
GO TO 220
C
C INTERPOLATE BETWEEN 1970 - 1975:
C
60 F2 = (DBLE(IY)+DBLE(IDAY-1)/365.25D0-1970)/5.D0
F1 = 1.D0 - F2
DO 70 N = 1,105
G(N) = G70(N)*F1 + G75(N)*F2
H(N) = H70(N)*F1 + H75(N)*F2
70 CONTINUE
GO TO 220
C
C INTERPOLATE BETWEEN 1975 - 1980:
C
80 F2 = (DBLE(IY)+DBLE(IDAY-1)/365.25D0-1975)/5.D0
F1 = 1.D0 - F2
DO 90 N = 1,105
G(N) = G75(N)*F1 + G80(N)*F2
H(N) = H75(N)*F1 + H80(N)*F2
90 CONTINUE
GO TO 220
C
C INTERPOLATE BETWEEN 1980 - 1985:
C
100 F2 = (DBLE(IY)+DBLE(IDAY-1)/365.25D0-1980)/5.D0
F1 = 1.D0 - F2
DO 110 N = 1,105
G(N) = G80(N)*F1 + G85(N)*F2
H(N) = H80(N)*F1 + H85(N)*F2
110 CONTINUE
GO TO 220
C
C INTERPOLATE BETWEEN 1985 - 1990:
C
120 F2 = (DBLE(IY)+DBLE(IDAY-1)/365.25D0-1985)/5.D0
F1 = 1.D0 - F2
DO 130 N = 1,105
G(N) = G85(N)*F1 + G90(N)*F2
H(N) = H85(N)*F1 + H90(N)*F2
130 CONTINUE
GO TO 220
C
C INTERPOLATE BETWEEN 1990 - 1995:
C
140 F2 = (DBLE(IY)+DBLE(IDAY-1)/365.25D0-1990)/5.D0
F1 = 1.D0 - F2
DO 150 N = 1,105
G(N) = G90(N)*F1 + G95(N)*F2
H(N) = H90(N)*F1 + H95(N)*F2
150 CONTINUE
GO TO 220
C
C INTERPOLATE BETWEEN 1995 - 2000:
C
160 F2 = (DBLE(IY)+DBLE(IDAY-1)/365.25D0-1995)/5.D0
F1 = 1.D0 - F2
DO 170 N = 1,105
G(N) = G95(N)*F1 + G00(N)*F2
H(N) = H95(N)*F1 + H00(N)*F2
170 CONTINUE
GO TO 220
C
C INTERPOLATE BETWEEN 2000 - 2005:
C
180 F2 = (DBLE(IY)+DBLE(IDAY-1)/365.25D0-2000)/5.D0
F1 = 1.D0 - F2
DO 190 N = 1,105
G(N) = G00(N)*F1 + G05(N)*F2
H(N) = H00(N)*F1 + H05(N)*F2
190 CONTINUE
GO TO 220
C
C INTERPOLATE BETWEEN 2005 - 2010:
C
200 F2 = (DBLE(IY)+DBLE(IDAY-1)/365.25D0-2005)/5.D0
F1 = 1.D0 - F2
DO 210 N = 1,105
G(N) = G05(N)*F1 + G10(N)*F2
H(N) = H05(N)*F1 + H10(N)*F2
210 CONTINUE
GO TO 220
C
C COEFFICIENTS FOR A GIVEN YEAR HAVE BEEN CALCULATED; NOW MULTIPLY
C THEM BY SCHMIDT NORMALIZATION FACTORS:
C
220 S = 1.D0
DO 240 N = 2,14
MN = N*(N-1)/2 + 1
S = S*DBLE(2*N-3)/DBLE(N-1)
G(MN) = G(MN)*S
H(MN) = H(MN)*S
P = S
DO 230 M = 2,N
AA = 1.D0
IF (M.EQ.2) AA = 2.D0
P = P*SQRT(AA*DBLE(N-M+1)/DBLE(N+M-2))
MNN = MN + M - 1
G(MNN) = G(MNN)*P
H(MNN) = H(MNN)*P
230 CONTINUE
240 CONTINUE
G_10 = -G(2)
G_11 = G(3)
H_11 = H(3)
C
CNOW CALCULATE GEO COMPONENTS OF THE UNIT VECTOR EzMAG, PARALLEL TO GEODIPOLE AXIS:
C SIN(TETA0)*COS(LAMBDA0), SIN(TETA0)*SIN(LAMBDA0), AND COS(TETA0)
C ST0 * CL0 ST0 * SL0 CT0
C
SQ = G_11**2 + H_11**2
SQQ = SQRT(SQ)
SQR = SQRT(G_10**2+SQ)
SL0 = -H_11/SQQ
CL0 = -G_11/SQQ
ST0 = SQQ/SQR
CT0 = G_10/SQR
STCL = ST0*CL0
STSL = ST0*SL0
CTSL = CT0*SL0
CTCL = CT0*CL0
C
C NOW CALCULATE GEI COMPONENTS (S1,S2,S3) OF THE UNIT VECTOR S = EX_GSE
C POINTING FROM THE EARTH'S CENTER TO TS_SUN
C
CALL TS_SUN_08(IY,IDAY,IHOUR,MIN,ISEC,GST,SLONG,SRASN,SDEC)
C
S1 = COS(SRASN)*COS(SDEC)
S2 = SIN(SRASN)*COS(SDEC)
S3 = SIN(SDEC)
C
C NOW CALCULATE GEI COMPONENTS (DZ1,DZ2,DZ3) OF THE UNIT VECTOR EZGSE
C POINTING NORTHWARD AND ORTHOGONAL TO THE ECLIPTIC PLANE, AS
C (0,-SIN(OBLIQ),COS(OBLIQ)). FOR THE EPOCH 1978, OBLIQ = 23.44214 DEGS.
C HERE WE USE A MORE ACCURATE TIME-DEPENDENT VALUE, DETERMINED AS:
C
DJ = DBLE(365*(IY-1900)+(IY-1901)/4+IDAY) - 0.5D0 +
* DBLE(IHOUR*3600+MIN*60+ISEC)/86400.D0
T = DJ/36525.D0
OBLIQ = (23.45229D0-0.0130125D0*T)/57.2957795D0
DZ1 = 0.D0
DZ2 = -SIN(OBLIQ)
DZ3 = COS(OBLIQ)
C
C NOW WE OBTAIN GEI COMPONENTS OF THE UNIT VECTOR EYGSE=(DY1,DY2,DY3),
C COMPLETING THE RIGHT-HANDED SYSTEM. THEY CAN BE FOUND FROM THE VECTOR
C PRODUCT EZGSE x EXGSE = (DZ1,DZ2,DZ3) x (S1,S2,S3):
C
DY1 = DZ2*S3 - DZ3*S2
DY2 = DZ3*S1 - DZ1*S3
DY3 = DZ1*S2 - DZ2*S1
C
CNOW LET'S CALCULATE GEI COMPONENTS OF THE UNIT VECTOR X = EXGSW, DIRECTED ANTIPARALLEL
CTO THE OBSERVED SOLAR WIND FLOW. FIRST, CALCULATE ITS COMPONENTS IN GSE:
C
V = SQRT(VGSEX**2+VGSEY**2+VGSEZ**2)
DX1 = -VGSEX/V
DX2 = -VGSEY/V
DX3 = -VGSEZ/V
C
C THEN IN GEI:
C
X1 = DX1*S1 + DX2*DY1 + DX3*DZ1
X2 = DX1*S2 + DX2*DY2 + DX3*DZ2
X3 = DX1*S3 + DX2*DY3 + DX3*DZ3
C
CNOW CALCULATE GEI COMPONENTS (DIP1,DIP2,DIP3) OF THE UNIT VECTOR DIP = EZ_SM = EZ_MAG,
C ALIGNED WITH THE GEODIPOLE AND POINTING NORTHWARD FROM ECLIPTIC PLANE:
C
CGST = COS(GST)
SGST = SIN(GST)
C
DIP1 = STCL*CGST - STSL*SGST
DIP2 = STCL*SGST + STSL*CGST
DIP3 = CT0
C
CTHIS ALLOWS US TO CALCULATE GEI COMPONENTS OF THE UNIT VECTOR Y = EYGSW
CBY TAKING THE VECTOR PRODUCT DIP x X AND NORMALIZING IT TO UNIT LENGTH:
C
Y1 = DIP2*X3 - DIP3*X2
Y2 = DIP3*X1 - DIP1*X3
Y3 = DIP1*X2 - DIP2*X1
Y = SQRT(Y1*Y1+Y2*Y2+Y3*Y3)
Y1 = Y1/Y
Y2 = Y2/Y
Y3 = Y3/Y
C
CAND GEI COMPONENTS OF THE UNIT VECTOR Z = EZGSW = EXGSW x EYGSW = X x Y:
C
Z1 = X2*Y3 - X3*Y2
Z2 = X3*Y1 - X1*Y3
Z3 = X1*Y2 - X2*Y1
C
C ELEMENTS OF THE MATRIX GSE TO GSW ARE THE SCALAR PRODUCTS:
C
C E11=(EXGSE,EXGSW) E12=(EXGSE,EYGSW) E13=(EXGSE,EZGSW)
C E21=(EYGSE,EXGSW) E22=(EYGSE,EYGSW) E23=(EYGSE,EZGSW)
C E31=(EZGSE,EXGSW) E32=(EZGSE,EYGSW) E33=(EZGSE,EZGSW)
C
E11 = S1*X1 + S2*X2 + S3*X3
E12 = S1*Y1 + S2*Y2 + S3*Y3
E13 = S1*Z1 + S2*Z2 + S3*Z3
E21 = DY1*X1 + DY2*X2 + DY3*X3
E22 = DY1*Y1 + DY2*Y2 + DY3*Y3
E23 = DY1*Z1 + DY2*Z2 + DY3*Z3
E31 = DZ1*X1 + DZ2*X2 + DZ3*X3
E32 = DZ1*Y1 + DZ2*Y2 + DZ3*Y3
E33 = DZ1*Z1 + DZ2*Z2 + DZ3*Z3
C
C GEODIPOLE TILT ANGLE IN THE GSW SYSTEM: PSI=ARCSIN(DIP,EXGSW)
C
SPS = DIP1*X1 + DIP2*X2 + DIP3*X3
CPS = SQRT(1.D0-SPS**2)
PSI = ASIN(SPS)
C
C ELEMENTS OF THE MATRIX GEO TO GSW ARE THE SCALAR PRODUCTS:
C
C A11=(EXGEO,EXGSW), A12=(EYGEO,EXGSW), A13=(EZGEO,EXGSW),
C A21=(EXGEO,EYGSW), A22=(EYGEO,EYGSW), A23=(EZGEO,EYGSW),
C A31=(EXGEO,EZGSW), A32=(EYGEO,EZGSW), A33=(EZGEO,EZGSW),
C
C ALL THE UNIT VECTORS IN BRACKETS ARE ALREADY DEFINED IN GEI:
C
C EXGEO=(CGST,SGST,0), EYGEO=(-SGST,CGST,0), EZGEO=(0,0,1)
C EXGSW=(X1,X2,X3), EYGSW=(Y1,Y2,Y3), EZGSW=(Z1,Z2,Z3)
C AND THEREFORE:
C
A11 = X1*CGST + X2*SGST
A12 = -X1*SGST + X2*CGST
A13 = X3
A21 = Y1*CGST + Y2*SGST
A22 = -Y1*SGST + Y2*CGST
A23 = Y3
A31 = Z1*CGST + Z2*SGST
A32 = -Z1*SGST + Z2*CGST
A33 = Z3
C
CNOW CALCULATE ELEMENTS OF THE MATRIX MAG TO SM (ONE ROTATION ABOUT THE GEODIPOLE AXIS);
CTHEY ARE FOUND AS THE SCALAR PRODUCTS: CFI=GM22=(EYSM,EYMAG)=(EYGSW,EYMAG),
C SFI=GM23=(EYSM,EXMAG)=(EYGSW,EXMAG),
C DERIVED AS FOLLOWS:
C
CIN GEO, THE VECTORS EXMAG AND EYMAG HAVE THE COMPONENTS (CT0*CL0,CT0*SL0,-ST0)
C AND (-SL0,CL0,0), RESPECTIVELY. HENCE, IN GEI THEIR COMPONENTS ARE:
C EXMAG: CT0*CL0*COS(GST)-CT0*SL0*SIN(GST)
C CT0*CL0*SIN(GST)+CT0*SL0*COS(GST)
C -ST0
C EYMAG: -SL0*COS(GST)-CL0*SIN(GST)
C -SL0*SIN(GST)+CL0*COS(GST)
C 0
CNOW, NOTE THAT GEI COMPONENTS OF EYSM=EYGSW WERE FOUND ABOVE AS Y1, Y2, AND Y3,
C AND WE ONLY HAVE TO COMBINE THESE QUANTITIES INTO SCALAR PRODUCTS:
C
EXMAGX = CT0*(CL0*CGST-SL0*SGST)
EXMAGY = CT0*(CL0*SGST+SL0*CGST)
EXMAGZ = -ST0
EYMAGX = -(SL0*CGST+CL0*SGST)
EYMAGY = -(SL0*SGST-CL0*CGST)
CFI = Y1*EYMAGX + Y2*EYMAGY
SFI = Y1*EXMAGX + Y2*EXMAGY + Y3*EXMAGZ
C
9000 FORMAT (/,/,1X,
* '****TS_RECALC WARNS: YEAR IS OUT OF INTERVAL 1965-2015: IYEAR='
* ,I4,/,6X,'CALCULATIONS WILL BE DONE FOR IYEAR=',I4,/)
RETURN
END
c
c==================================================================================
SUBROUTINE GSWGSE(XGSW,YGSW,ZGSW,XGSE,YGSE,ZGSE,J)
C
CTHIS SUBROUTINE TRANSFORMS COMPONENTS OF ANY VECTOR BETWEEN THE STANDARD GSE
CCOORDINATE SYSTEM AND THE GEOCENTRIC SOLAR-WIND (GSW, aka GSWM), DEFINED AS FOLLOWS
C(HONES ET AL., PLANET.SPACE SCI., V.34, P.889, 1986; TSYGANENKO ET AL., JGRA,
C V.103(A4), P.6827, 1998):
C
CIN THE GSW SYSTEM, X AXIS IS ANTIPARALLEL TO THE OBSERVED DIRECTION OF THE SOLAR WIND FLOW.
CTWO OTHER AXES, Y AND Z, ARE DEFINED IN THE SAME WAY AS FOR THE STANDARD GSM, THAT IS,
CZ AXIS ORTHOGONAL TO X AXIS, POINTS NORTHWARD, AND LIES IN THE PLANE DEFINED BY THE X-
C AND GEODIPOLE AXIS. THE Y AXIS COMPLETES THE RIGHT-HANDED SYSTEM.
C
C THE GSW SYSTEM BECOMES IDENTICAL TO THE STANDARD GSM IN THE CASE OF
C A STRICTLY RADIAL SOLAR WIND FLOW.
C
C AUTHOR: N. A. TSYGANENKO
C ADDED TO 2008 VERSION OF GEOPACK: JAN 27, 2008.
C
C J>0 J<0
C -----INPUT: J,XGSW,YGSW,ZGSW J,XGSE,YGSE,ZGSE
C -----OUTPUT: XGSE,YGSE,ZGSE XGSW,YGSW,ZGSW
C
C IMPORTANT THINGS TO REMEMBER:
C
C(1) BEFORE CALLING GSWGSE, BE SURE TO INVOKE SUBROUTINE TS_RECALC_08, IN ORDER
C TO DEFINE ALL NECESSARY ELEMENTS OF TRANSFORMATION MATRICES
C
C(2) IN THE ABSENCE OF INFORMATION ON THE SOLAR WIND DIRECTION, E.G., WITH ONLY SCALAR
C SPEED V KNOWN, THIS SUBROUTINE CAN BE USED TO CONVERT VECTORS TO ABERRATED
C COORDINATE SYSTEM, TAKING INTO ACCOUNT EARTH'S ORBITAL SPEED OF 29 KM/S.
C TO DO THAT, SPECIFY THE LAST 3 PARAMETERS IN TS_RECALC_08 AS FOLLOWS:
C VGSEX=-V, VGSEY=29.0, VGSEZ=0.0.
C
C IT SHOULD ALSO BE KEPT IN MIND THAT IN SOME SOLAR WIND DATABASES THE ABERRATION
C EFFECT HAS ALREADY BEEN TAKEN INTO ACCOUNT BY SUBTRACTING 29 KM/S FROM VYGSE;
C IN THAT CASE, THE ORIGINAL VYGSE VALUES SHOULD BE RESTORED BY ADDING BACK THE
C 29 KM/S CORRECTION. WHETHER OR NOT TO DO THAT, MUST BE VERIFIED WITH THE DATA
C ORIGINATOR (OR CAN BE DETERMINED BY CALCULATING THE AVERAGE VGSEY OVER
C A SUFFICIENTLY LONG TIME INTERVAL)
C
C(3) IF NO INFORMATION IS AVAILABLE ON THE SOLAR WIND SPEED, THEN SET VGSEX=-400.0
c AND VGSEY=VGSEZ=0. IN THAT CASE, THE GSW COORDINATE SYSTEM BECOMES
c IDENTICAL TO THE STANDARD ONE.
C
C
C DIRECT TRANSFORMATION:
C
C .. Scalar Arguments ..
DOUBLE PRECISION XGSE,XGSW,YGSE,YGSW,ZGSE,ZGSW
INTEGER J
C ..
C .. Common blocks ..
COMMON /GEOPACK1/AAA,E11,E21,E31,E12,E22,E32,E13,E23,E33
DOUBLE PRECISION E11,E12,E13,E21,E22,E23,E31,E32,E33
DOUBLE PRECISION AAA(25)
C ..
IF (J.GT.0) THEN
XGSE = XGSW*E11 + YGSW*E12 + ZGSW*E13
YGSE = XGSW*E21 + YGSW*E22 + ZGSW*E23
ZGSE = XGSW*E31 + YGSW*E32 + ZGSW*E33
END IF
C
C INVERSE TRANSFORMATION: CARRIED OUT USING THE TRANSPOSED MATRIX:
C
IF (J.LT.0) THEN
XGSW = XGSE*E11 + YGSE*E21 + ZGSE*E31
YGSW = XGSE*E12 + YGSE*E22 + ZGSE*E32
ZGSW = XGSE*E13 + YGSE*E23 + ZGSE*E33
END IF
C
RETURN
END
C
C========================================================================================
C
SUBROUTINE TS_GEOMAG_08(XGEO,YGEO,ZGEO,XMAG,YMAG,ZMAG,J)
C
C CONVERTS GEOGRAPHIC (GEO) TO DIPOLE (MAG) COORDINATES OR VICE VERSA.
C
C J>0 J<0
C -----INPUT: J,XGEO,YGEO,ZGEO J,XMAG,YMAG,ZMAG
C -----OUTPUT: XMAG,YMAG,ZMAG XGEO,YGEO,ZGEO
C
CATTENTION: SUBROUTINE TS_RECALC_08 MUST BE INVOKED BEFORE TS_GEOMAG_08 IN TWO CASES:
C /A/ BEFORE THE FIRST TRANSFORMATION OF COORDINATES
C /B/ IF THE VALUES OF IYEAR AND/OR IDAY HAVE BEEN CHANGED
C
C NO INFORMATION IS REQUIRED HERE ON THE SOLAR WIND VELOCITY, SO ONE
C CAN SET VGSEX=-400.0, VGSEY=0.0, VGSEZ=0.0 IN TS_RECALC_08.
C
C LAST MOFIFICATION: JAN 28, 2008
C
C AUTHOR: N. A. TSYGANENKO
C
C .. Scalar Arguments ..
DOUBLE PRECISION XGEO,XMAG,YGEO,YMAG,ZGEO,ZMAG
INTEGER J
C ..
C .. Common blocks ..
COMMON /GEOPACK1/ST0,CT0,SL0,CL0,CTCL,STCL,CTSL,STSL,AB
DOUBLE PRECISION CL0,CT0,CTCL,CTSL,SL0,ST0,STCL,STSL
DOUBLE PRECISION AB(26)
C ..
IF (J.GT.0) THEN
XMAG = XGEO*CTCL + YGEO*CTSL - ZGEO*ST0
YMAG = YGEO*CL0 - XGEO*SL0
ZMAG = XGEO*STCL + YGEO*STSL + ZGEO*CT0
ELSE
XGEO = XMAG*CTCL - YMAG*SL0 + ZMAG*STCL
YGEO = XMAG*CTSL + YMAG*CL0 + ZMAG*STSL
ZGEO = ZMAG*CT0 - XMAG*ST0
END IF
RETURN
END
c
c=========================================================================================
c
SUBROUTINE GEIGEO_08(XGEI,YGEI,ZGEI,XGEO,YGEO,ZGEO,J)
C
C CONVERTS EQUATORIAL INERTIAL (GEI) TO GEOGRAPHICAL (GEO) COORDS
C OR VICE VERSA.
C J>0 J<0
C ----INPUT: J,XGEI,YGEI,ZGEI J,XGEO,YGEO,ZGEO
C ----OUTPUT: XGEO,YGEO,ZGEO XGEI,YGEI,ZGEI
C
CATTENTION: SUBROUTINE TS_RECALC_08 MUST BE INVOKED BEFORE GEIGEO_08 IN TWO CASES:
C /A/ BEFORE THE FIRST TRANSFORMATION OF COORDINATES
C/B/ IF THE CURRENT VALUES OF IYEAR,IDAY,IHOUR,MIN,ISEC HAVE BEEN CHANGED
C
C NO INFORMATION IS REQUIRED HERE ON THE SOLAR WIND VELOCITY, SO ONE
C CAN SET VGSEX=-400.0, VGSEY=0.0, VGSEZ=0.0 IN TS_RECALC_08.
C
C LAST MODIFICATION: JAN 28, 2008
C AUTHOR: N. A. TSYGANENKO
C
C
C .. Scalar Arguments ..
DOUBLE PRECISION XGEI,XGEO,YGEI,YGEO,ZGEI,ZGEO
INTEGER J
C ..
C .. Common blocks ..
COMMON /GEOPACK1/A,CGST,SGST,B
DOUBLE PRECISION CGST,SGST
DOUBLE PRECISION A(13),B(19)
C ..
IF (J.GT.0) THEN
XGEO = XGEI*CGST + YGEI*SGST
YGEO = YGEI*CGST - XGEI*SGST
ZGEO = ZGEI
ELSE
XGEI = XGEO*CGST - YGEO*SGST
YGEI = YGEO*CGST + XGEO*SGST
ZGEI = ZGEO
END IF
RETURN
END
C
C=======================================================================================
C
SUBROUTINE TS_MAGSM_08(XMAG,YMAG,ZMAG,XSM,YSM,ZSM,J)
C
C CONVERTS DIPOLE (MAG) TO SOLAR MAGNETIC (SM) COORDINATES OR VICE VERSA
C
C J>0 J<0
C ----INPUT: J,XMAG,YMAG,ZMAG J,XSM,YSM,ZSM
C ----OUTPUT: XSM,YSM,ZSM XMAG,YMAG,ZMAG
C
CATTENTION: SUBROUTINE TS_RECALC_08 MUST BE INVOKED BEFORE TS_MAGSM_08 IN THREE CASES:
C /A/ BEFORE THE FIRST TRANSFORMATION OF COORDINATES, OR
C /B/ IF THE VALUES OF IYEAR,IDAY,IHOUR,MIN,ISEC HAVE CHANGED, AND/OR
C/C/ IF THE VALUES OF COMPONENTS OF THE SOLAR WIND FLOW VELOCITY HAVE CHANGED
C
C IMPORTANT NOTE:
C
C A NON-STANDARD DEFINITION IS IMPLIED HERE FOR THE SOLAR MAGNETIC COORDINATE
C SYSTEM: IT IS ASSUMED THAT THE XSM AXIS LIES IN THE PLANE DEFINED BY THE
C GEODIPOLE AXIS AND THE OBSERVED VECTOR OF THE SOLAR WIND FLOW (RATHER THAN
C THE EARTH-TS_SUN LINE). IN ORDER TO CONVERT MAG COORDINATES TO AND FROM THE
C STANDARD SM COORDINATES, INVOKE TS_RECALC_08 WITH VGSEX=-400.0, VGSEY=0.0, VGSEZ=0.0
C
C LAST MODIFICATION: FEB 07, 2008
C
C AUTHOR: N. A. TSYGANENKO
C
C
C .. Scalar Arguments ..
DOUBLE PRECISION XMAG,XSM,YMAG,YSM,ZMAG,ZSM
INTEGER J
C ..
C .. Common blocks ..
COMMON /GEOPACK1/A,SFI,CFI,B
DOUBLE PRECISION CFI,SFI
DOUBLE PRECISION A(8),B(24)
C ..
IF (J.GT.0) THEN
XSM = XMAG*CFI - YMAG*SFI
YSM = XMAG*SFI + YMAG*CFI
ZSM = ZMAG
ELSE
XMAG = XSM*CFI + YSM*SFI
YMAG = YSM*CFI - XSM*SFI
ZMAG = ZSM
END IF
RETURN
END
C
C=====================================================================================
C
SUBROUTINE SMGSW_08(XSM,YSM,ZSM,XGSW,YGSW,ZGSW,J)
C
CCONVERTS SOLAR MAGNETIC (SM) TO GEOCENTRIC SOLAR-WIND (GSW) COORDINATES OR VICE VERSA.
C
C J>0 J<0
C -----INPUT: J,XSM,YSM,ZSM J,XGSW,YGSW,ZGSW
C ----OUTPUT: XGSW,YGSW,ZGSW XSM,YSM,ZSM
C
CATTENTION: SUBROUTINE TS_RECALC_08 MUST BE INVOKED BEFORE SMGSW_08 IN THREE CASES:
C /A/ BEFORE THE FIRST TRANSFORMATION OF COORDINATES
C /B/ IF THE VALUES OF IYEAR,IDAY,IHOUR,MIN,ISEC HAVE BEEN CHANGED
C/C/ IF THE VALUES OF COMPONENTS OF THE SOLAR WIND FLOW VELOCITY HAVE CHANGED
C
C IMPORTANT NOTE:
C
C A NON-STANDARD DEFINITION IS IMPLIED HERE FOR THE SOLAR MAGNETIC (SM) COORDINATE
C SYSTEM: IT IS ASSUMED THAT THE XSM AXIS LIES IN THE PLANE DEFINED BY THE
C GEODIPOLE AXIS AND THE OBSERVED VECTOR OF THE SOLAR WIND FLOW (RATHER THAN
C THE EARTH-TS_SUN LINE). IN ORDER TO CONVERT MAG COORDINATES TO AND FROM THE
C STANDARD SM COORDINATES, INVOKE TS_RECALC_08 WITH VGSEX=-400.0, VGSEY=0.0, VGSEZ=0.0
C
C LAST MODIFICATION: FEB 07, 2008
C
C AUTHOR: N. A. TSYGANENKO
C
C .. Scalar Arguments ..
DOUBLE PRECISION XGSW,XSM,YGSW,YSM,ZGSW,ZSM
INTEGER J
C ..
C .. Common blocks ..
COMMON /GEOPACK1/A,SPS,CPS,B
DOUBLE PRECISION CPS,SPS
DOUBLE PRECISION A(10),B(22)
C ..
IF (J.GT.0) THEN
XGSW = XSM*CPS + ZSM*SPS
YGSW = YSM
ZGSW = ZSM*CPS - XSM*SPS
ELSE
XSM = XGSW*CPS - ZGSW*SPS
YSM = YGSW
ZSM = XGSW*SPS + ZGSW*CPS
END IF
RETURN
END
C
C==========================================================================================
C
SUBROUTINE GEOGSW(XGEO,YGEO,ZGEO,XGSW,YGSW,ZGSW,J)
C
CCONVERTS GEOGRAPHIC (GEO) TO GEOCENTRIC SOLAR-WIND (GSW) COORDINATES OR VICE VERSA.
C
C J>0 J<0
C ----- INPUT: J,XGEO,YGEO,ZGEO J,XGSW,YGSW,ZGSW
C ---- OUTPUT: XGSW,YGSW,ZGSW XGEO,YGEO,ZGEO
C
CATTENTION: SUBROUTINE TS_RECALC_08 MUST BE INVOKED BEFORE GEOGSW IN THREE CASES:
C /A/ BEFORE THE FIRST TRANSFORMATION OF COORDINATES, OR
C /B/ IF THE VALUES OF IYEAR,IDAY,IHOUR,MIN,ISEC HAVE CHANGED, AND/OR
C/C/ IF THE VALUES OF COMPONENTS OF THE SOLAR WIND FLOW VELOCITY HAVE CHANGED
C
CNOTE: THIS SUBROUTINE CONVERTS GEO VECTORS TO AND FROM THE SOLAR-WIND GSW COORDINATE
C SYSTEM, TAKING INTO ACCOUNT POSSIBLE DEFLECTIONS OF THE SOLAR WIND DIRECTION FROM
C STRICTLY RADIAL. BEFORE CONVERTING TO/FROM STANDARD GSM COORDINATES, INVOKE TS_RECALC_08
C WITH VGSEX=-400.0 and VGSEY=0.0, VGSEZ=0.0
C
C LAST MODIFICATION: FEB 07, 2008
C
C AUTHOR: N. A. TSYGANENKO
C
C
C .. Scalar Arguments ..
DOUBLE PRECISION XGEO,XGSW,YGEO,YGSW,ZGEO,ZGSW
INTEGER J
C ..
C .. Common blocks ..
COMMON /GEOPACK1/AA,A11,A21,A31,A12,A22,A32,A13,A23,A33,B
DOUBLE PRECISION A11,A12,A13,A21,A22,A23,A31,A32,A33
DOUBLE PRECISION AA(16),B(9)
C ..
IF (J.GT.0) THEN
XGSW = A11*XGEO + A12*YGEO + A13*ZGEO
YGSW = A21*XGEO + A22*YGEO + A23*ZGEO
ZGSW = A31*XGEO + A32*YGEO + A33*ZGEO
ELSE
XGEO = A11*XGSW + A21*YGSW + A31*ZGSW
YGEO = A12*XGSW + A22*YGSW + A32*ZGSW
ZGEO = A13*XGSW + A23*YGSW + A33*ZGSW
END IF
RETURN
END
C
C=====================================================================================
C
SUBROUTINE GEODGEO_08(H,XMU,R,THETA,J)
C .. Scalar Arguments ..
DOUBLE PRECISION H,R,THETA,XMU
INTEGER J
C ..
C .. Local Scalars ..
DOUBLE PRECISION ARG,BETA,COSFIMS,COSLAM,COSXMU,DEN,DPHI,PHI,PHI1,
* RR,RS,R_EQ,SINLAM,SINXMU,SP,TOL,X,XMUS,Z
INTEGER N
C ..
C .. Intrinsic Functions ..
INTRINSIC ABS,ACOS,ASIN,COS,SIN,SQRT
C ..
C .. Data statements ..
C
CTHIS SUBROUTINE (1) CONVERTS VERTICAL LOCAL HEIGHT (ALTITUDE) H AND GEODETIC
CLATITUDE XMU INTO GEOCENTRIC COORDINATES R AND THETA (GEOCENTRIC RADIAL
CDISTANCE AND COLATITUDE, RESPECTIVELY; ALSO KNOWN AS ECEF COORDINATES),
CAS WELL AS (2) PERFORMS THE INVERSE TRANSFORMATION FROM {R,THETA} TO {H,XMU}.
C
CTHE SUBROUTINE USES WORLD GEODETIC SYSTEM WGS84 PARAMETERS FOR THE EARTH'S
CELLIPSOID. THE ANGULAR QUANTITIES (GEO COLATITUDE THETA AND GEODETIC LATITUDE
CXMU) ARE IN RADIANS, AND THE DISTANCES (GEOCENTRIC RADIUS R AND ALTITUDE H
C ABOVE THE EARTH'S ELLIPSOID) ARE IN KILOMETERS.
C
CIF J>0, THE TRANSFORMATION IS MADE FROM GEODETIC TO GEOCENTRIC COORDINATES
C USING SIMPLE DIRECT EQUATIONS.
CIF J<0, THE INVERSE TRANSFORMATION FROM GEOCENTRIC TO GEODETIC COORDINATES
C IS MADE BY MEANS OF A FAST ITERATIVE ALGORITHM.
C
c-------------------------------------------------------------------------------
C J>0 | J<0
c-------------------------------------------|-----------------------------------
C--INPUT: J H XMU | J R THETA
c flag altitude (km) geodetic | flag geocentric spherical
c latitude | distance (km) colatitude
c (radians) | (radians)
c-------------------------------------------|-----------------------------------
c |
C----OUTPUT: R THETA | H XMU
C geocentric spherical | altitude (km) geodetic
C distance (km) colatitude | latitude
C (radians) | (radians)
C-------------------------------------------------------------------------------
C
C AUTHOR: N. A. TSYGANENKO
c DATE: DEC 5, 2007
C
c
c R_EQ is the semi-major axis of the Earth's ellipsoid, and BETA is its
c second eccentricity squared
c
DATA R_EQ,BETA/6378.137D0,6.73949674228D-3/
DATA TOL/1.D-6/
C ..
c
c Direct transformation (GEOD=>GEO):
c
IF (J.GT.0) THEN
COSXMU = COS(XMU)
SINXMU = SIN(XMU)
DEN = SQRT(COSXMU**2+(SINXMU/(1.0D0+BETA))**2)
COSLAM = COSXMU/DEN
SINLAM = SINXMU/(DEN*(1.0D0+BETA))
RS = R_EQ/SQRT(1.0D0+BETA*SINLAM**2)
X = RS*COSLAM + H*COSXMU
Z = RS*SINLAM + H*SINXMU
R = SQRT(X**2+Z**2)
THETA = ACOS(Z/R)
END IF
c
c Inverse transformation (GEO=>GEOD):
c
IF (J.LT.0) THEN
N = 0
PHI = 1.570796327D0 - THETA
PHI1 = PHI
10 SP = SIN(PHI1)
C *PT*WARNING* Constant already double-precision
ARG = SP*(1.0D0+BETA)/SQRT(1.0D0+BETA*(2.0D0+BETA)*SP**2)
XMUS = ASIN(ARG)
RS = R_EQ/SQRT(1.0D0+BETA*SIN(PHI1)**2)
COSFIMS = COS(PHI1-XMUS)
H = SQRT((RS*COSFIMS)**2+R**2-RS**2) - RS*COSFIMS
Z = RS*SIN(PHI1) + H*SIN(XMUS)
X = RS*COS(PHI1) + H*COS(XMUS)
RR = SQRT(X**2+Z**2)
DPHI = ASIN(Z/RR) - PHI
PHI1 = PHI1 - DPHI
N = N + 1
IF (ABS(DPHI).GT.TOL .AND. N.LT.100) GO TO 10
XMU = XMUS
END IF
RETURN
END
C
C=====================================================================================
C
SUBROUTINE RHAND_08(X,Y,Z,R1,R2,R3,IOPT,PARMOD,EXNAME,INNAME)
C
CCALCULATES THE COMPONENTS OF THE RIGHT HAND SIDE VECTOR IN THE TS_GEOMAGNETIC FIELD
C LINE EQUATION (a subsidiary subroutine for the subroutine STEP_08)
C
C LAST MODIFICATION: FEB 07, 2008
C
C AUTHOR: N. A. TSYGANENKO
C
C
CEXNAME AND INNAME ARE NAMES OF SUBROUTINES FOR THE EXTERNAL AND INTERNAL
C PARTS OF THE TOTAL FIELD, E.G., T96_01 AND IGRF_GSW_08
C
C common block added by B. Rideout to signal error
COMMON /BR_ERR/IERR
INTEGER IERR
C
C .. Scalar Arguments ..
DOUBLE PRECISION R1,R2,R3,X,Y,Z
INTEGER IOPT
C ..
C .. Array Arguments ..
DOUBLE PRECISION PARMOD(10)
C ..
C .. Subroutine Arguments ..
EXTERNAL EXNAME,INNAME
C ..
C .. Local Scalars ..
DOUBLE PRECISION B,BX,BXGSW,BY,BYGSW,BZ,BZGSW,HXGSW,HYGSW,HZGSW
C ..
C .. Intrinsic Functions ..
INTRINSIC SQRT
C ..
C .. Common blocks ..
COMMON /GEOPACK1/A,DS3,BB,PSI,CC
DOUBLE PRECISION DS3,PSI
DOUBLE PRECISION A(12),BB(2),CC(18)
C ..
CALL EXNAME(IOPT,PARMOD,PSI,X,Y,Z,BXGSW,BYGSW,BZGSW)
C error check added by brideout
IF (IERR .EQ. 1) THEN
RETURN
END IF
CALL INNAME(X,Y,Z,HXGSW,HYGSW,HZGSW)
BX = BXGSW + HXGSW
BY = BYGSW + HYGSW
BZ = BZGSW + HZGSW
B = DS3/SQRT(BX**2+BY**2+BZ**2)
R1 = BX*B
R2 = BY*B
R3 = BZ*B
RETURN
END
C
C===================================================================================
C
SUBROUTINE STEP_08(X,Y,Z,DS,DSMAX,ERRIN,IOPT,PARMOD,EXNAME,INNAME)
C
CRE-CALCULATES THE INPUT VALUES {X,Y,Z} (IN GSW COORDINATES) FOR ANY POINT ON A FIELD LINE,
CBY MAKING A STEP ALONG THAT LINE USING RUNGE-KUTTA-MERSON ALGORITHM (G.N. Lance, Numerical
C methods for high-speed computers, Iliffe & Sons, London 1960.)
CDS IS A PRESCRIBED VALUE OF THE CURRENT STEP SIZE, DSMAX IS ITS UPPER LIMIT.
CERRIN IS A PERMISSIBLE ERROR (ITS OPTIMAL VALUE SPECIFIED IN THE S/R TRACE)
CIF THE ACTUAL ERROR (ERRCUR) AT THE CURRENT STEP IS LARGER THAN ERRIN, THE STEP IS REJECTED,
C AND THE CALCULATION IS REPEATED ANEW WITH HALVED STEPSIZE DS.
CIF ERRCUR IS SMALLER THAN ERRIN, THE STEP IS ACCEPTED, AND THE CURRENT VALUE OF DS IS RETAINE
C D
C FOR THE NEXT STEP.
CIF ERRCUR IS SMALLER THAN 0.04*ERRIN, THE STEP IS ACCEPTED, AND THE VALUE OF DS FOR THE NEXT
C STEP
C IS INCREASED BY THE FACTOR 1.5, BUT NOT LARGER THAN DSMAX.
CIOPT IS A FLAG, RESERVED FOR SPECIFYNG A VERSION OF THE EXTERNAL FIELD MODEL EXNAME.
C ARRAY PARMOD(10) CONTAINS INPUT PARAMETERS FOR THE MODEL EXNAME.
C EXNAME IS THE NAME OF THE SUBROUTINE FOR THE EXTERNAL FIELD MODEL.
CINNAME IS THE NAME OF THE SUBROUTINE FOR THE INTERNAL FIELD MODEL (EITHER DIP_08 OR IGRF_GSW_08
C )
C
CALL THE ABOVE PARAMETERS ARE INPUT ONES; OUTPUT IS THE TS_RECALCULATED VALUES OF X,Y,Z
C
C LAST MODIFICATION: APR 21, 2008 (SEE ERRATA AS OF THIS DATE)
C
C common block added by B. Rideout to signal error
COMMON /BR_ERR/IERR
INTEGER IERR
C
C .. Scalar Arguments ..
DOUBLE PRECISION DS,DSMAX,ERRIN,X,Y,Z
INTEGER IOPT
C ..
C .. Array Arguments ..
DOUBLE PRECISION PARMOD(10)
C ..
C .. Subroutine Arguments ..
EXTERNAL EXNAME,INNAME
C ..
C .. Local Scalars ..
DOUBLE PRECISION ERRCUR,R11,R12,R13,R21,R22,R23,R31,R32,R33,R41,
* R42,R43,R51,R52,R53
C ..
C .. External Subroutines ..
EXTERNAL RHAND_08
C ..
C .. Intrinsic Functions ..
INTRINSIC ABS,SIGN
C ..
C .. Common blocks ..
COMMON /GEOPACK1/A,DS3,B
DOUBLE PRECISION DS3
DOUBLE PRECISION A(12),B(21)
C ..
10 DS3 = -DS/3.D0
CALL RHAND_08(X,Y,Z,R11,R12,R13,IOPT,PARMOD,EXNAME,INNAME)
C error check added by brideout
IF (IERR .EQ. 1) THEN
RETURN
END IF
CALL RHAND_08(X+R11,Y+R12,Z+R13,R21,R22,R23,IOPT,PARMOD,EXNAME,
* INNAME)
C error check added by brideout
IF (IERR .EQ. 1) THEN
RETURN
END IF
CALL RHAND_08(X+.5D0*(R11+R21),Y+.5D0*(R12+R22),Z+.5D0*(R13+R23),
* R31,R32,R33,IOPT,PARMOD,EXNAME,INNAME)
C error check added by brideout
IF (IERR .EQ. 1) THEN
RETURN
END IF
CALL RHAND_08(X+.375D0*(R11+3.D0*R31),Y+.375D0*(R12+3.D0*R32),
* Z+.375D0*(R13+3.D0*R33),R41,R42,R43,IOPT,PARMOD,
* EXNAME,INNAME)
C error check added by brideout
IF (IERR .EQ. 1) THEN
RETURN
END IF
CALL RHAND_08(X+1.5D0*(R11-3.D0*R31+4.D0*R41),
* Y+1.5D0*(R12-3.D0*R32+4.D0*R42),
* Z+1.5D0*(R13-3.D0*R33+4.D0*R43),R51,R52,R53,IOPT,
* PARMOD,EXNAME,INNAME)
C error check added by brideout
IF (IERR .EQ. 1) THEN
RETURN
END IF
ERRCUR = ABS(R11-4.5D0*R31+4.D0*R41-.5D0*R51) +
* ABS(R12-4.5D0*R32+4.D0*R42-.5D0*R52) +
* ABS(R13-4.5D0*R33+4.D0*R43-.5D0*R53)
C AUTHOR: N. A. TSYGANENKO
C
C
IF (ERRCUR.GT.ERRIN) THEN
DS = DS*.5D0
GO TO 10
END IF
C READY FOR MAKING THE STEP, BUT CHECK THE ACCURACY; IF INSUFFICIENT,
C REPEAT THE STEP WITH HALVED STEPSIZE:
C
C
IF (ABS(DS).GT.DSMAX) THEN
DS = SIGN(DSMAX,DS)
GO TO 10
END IF
C ACCURACY IS ACCEPTABLE, BUT CHECK IF THE STEPSIZE IS NOT TOO LARGE;
C OTHERWISE REPEAT THE STEP WITH DS=DSMAX
C
20 X = X + .5D0*(R11+4.D0*R41+R51)
Y = Y + .5D0*(R12+4.D0*R42+R52)
Z = Z + .5D0*(R13+4.D0*R43+R53)
C
C MAKING THE STEP:
C
C
IF (ERRCUR.LT.ERRIN*.04D0 .AND. DS.LT.DSMAX/1.5D0) DS = DS*1.5D0
RETURN
END
CIF THE ACTUAL ERROR IS TOO SMALL (LESS THAN 4% OF ERRIN) AND DS SMALLER
C THAN DSMAX/1.5, THEN WE INCREASE THE STEPSIZE FOR THE NEXT STEP BY 50%
C
SUBROUTINE TRACE(XI,YI,ZI,DIR,DSMAX,ERR,RLIM,R0,IOPT,PARMOD,
* EXNAME,INNAME,XF,YF,ZF,XX,YY,ZZ,L,LMAX)
C
C==============================================================================
C
C
C TRACES A FIELD LINE FROM AN ARBITRARY POINT OF SPACE TO THE EARTH'S
C SURFACE OR TO A MODEL LIMITING BOUNDARY.
C
C THIS SUBROUTINE ALLOWS TWO OPTIONS:
C
C(1) IF INNAME=IGRF_GSW_08, THEN THE IGRF MODEL WILL BE USED FOR CALCULATING
C CONTRIBUTION FROM EARTH'S INTERNAL SOURCES. IN THIS CASE, SUBROUTINE
C TS_RECALC_08 MUST BE CALLED BEFORE USING TRACE, WITH PROPERLY SPECIFIED DATE,
C UNIVERSAL TIME, AND SOLAR WIND VELOCITY COMPONENTS, TO CALCULATE IN ADVANCE
C ALL QUANTITIES NEEDED FOR THE MAIN FIELD MODEL AND FOR TRANSFORMATIONS
C BETWEEN INVOLVED COORDINATE SYSTEMS.
C
C(2) IF INNAME=DIP_08, THEN A PURE DIPOLE FIELD WILL BE USED INSTEAD OF THE IGRF MODEL.
C IN THIS CASE, THE SUBROUTINE TS_RECALC_08 MUST ALSO BE CALLED BEFORE TRACE.
C HERE ONE CAN CHOOSE EITHER TO
C (a) CALCULATE DIPOLE TILT ANGLE BASED ON DATE, TIME, AND SOLAR WIND DIRECTION,
COR (b) EXPLICITLY SPECIFY THAT ANGLE, WITHOUT ANY REFERENCE TO DATE/UT/SOLAR WIND.
C IN THE LAST CASE (b), THE SINE (SPS) AND COSINE (CPS) OF THE DIPOLE TILT
C ANGLE MUST BE SPECIFIED IN ADVANCE (BUT AFTER HAVING CALLED TS_RECALC_08) AND FORWARDED
C IN THE COMMON BLOCK /GEOPACK1/ (IN ITS 11th AND 12th ELEMENTS, RESPECTIVELY).
C IN THIS CASE THE ROLE OF THE SUBROUTINE TS_RECALC_08 IS REDUCED TO ONLY CALCULATING
C THE COMPONENTS OF THE EARTH'S DIPOLE MOMENT.
C
C ------------- INPUT PARAMETERS:
C
CXI,YI,ZI - GSW COORDS OF THE FIELD LINE STARTING POINT (IN EARTH RADII, 1 RE = 6371.2 km),
C
CDIR - SIGN OF THE TRACING DIRECTION: IF DIR=1.0 THEN THE TRACING IS MADE ANTIPARALLEL
CTO THE TOTAL FIELD VECTOR (E.G., FROM NORTHERN TO SOUTHERN CONJUGATE POINT);
CIF DIR=-1.0 THEN THE TRACING PROCEEDS IN THE OPPOSITE DIRECTION, THAT IS, PARALLEL TO
C THE TOTAL FIELD VECTOR.
C
CDSMAX - UPPER LIMIT ON THE STEPSIZE (SETS A DESIRED MAXIMAL SPACING BETWEEN
C THE FIELD LINE POINTS)
C
CERR - PERMISSIBLE STEP ERROR. A REASONABLE ESTIMATE PROVIDING A SUFFICIENT ACCURACY FOR MOST
C APPLICATIONS IS ERR=0.0001. SMALLER/LARGER VALUES WILL RESULT IN LARGER/SMALLER NUMBER
C OF STEPS AND, HENCE, OF OUTPUT FIELD LINE POINTS. NOTE THAT USING MUCH SMALLER VALUES
C OF ERR MAY REQUIRE USING A DOUBLE PRECISION VERSION OF THE ENTIRE PACKAGE.
C
CR0 - RADIUS OF A SPHERE (IN RE), DEFINING THE INNER BOUNDARY OF THE TRACING REGION
C (USUALLY, EARTH'S SURFACE OR THE IONOSPHERE, WHERE R0~1.0)
C IF THE FIELD LINE REACHES THAT SPHERE FROM OUTSIDE, ITS INBOUND TRACING IS
C TERMINATED AND THE CROSSING POINT COORDINATES XF,YF,ZF ARE CALCULATED.
C
CRLIM - RADIUS OF A SPHERE (IN RE), DEFINING THE OUTER BOUNDARY OF THE TRACING REGION;
C IF THE FIELD LINE REACHES THAT BOUNDARY FROM INSIDE, ITS OUTBOUND TRACING IS
C TERMINATED AND THE CROSSING POINT COORDINATES XF,YF,ZF ARE CALCULATED.
C
CIOPT - A MODEL INDEX; CAN BE USED FOR SPECIFYING A VERSION OF THE EXTERNAL FIELD
C MODEL (E.G., A NUMBER OF THE KP-INDEX INTERVAL). ALTERNATIVELY, ONE CAN USE THE ARRAY
C PARMOD FOR THAT PURPOSE (SEE BELOW); IN THAT CASE IOPT IS JUST A DUMMY PARAMETER.
C
CPARMOD - A 10-ELEMENT ARRAY CONTAINING INPUT PARAMETERS NEEDED FOR A UNIQUE
C SPECIFICATION OF THE EXTERNAL FIELD MODEL. THE CONCRETE MEANING OF THE COMPONENTS
C OF PARMOD DEPENDS ON A SPECIFIC VERSION OF THAT MODEL.
C
CEXNAME - NAME OF A SUBROUTINE PROVIDING COMPONENTS OF THE EXTERNAL MAGNETIC FIELD
C (E.G., T89, OR T96_01, ETC.).
CINNAME - NAME OF A SUBROUTINE PROVIDING COMPONENTS OF THE INTERNAL MAGNETIC FIELD
C (EITHER DIP_08 OR IGRF_GSW_08).
C
CLMAX - MAXIMAL LENGTH OF THE ARRAYS XX,YY,ZZ, IN WHICH COORDINATES OF THE FIELD
C LINE POINTS ARE STORED. LMAX SHOULD BE SET EQUAL TO THE ACTUAL LENGTH OF
C THE ARRAYS, DEFINED IN THE MAIN PROGRAM AS ACTUAL ARGUMENTS OF THIS SUBROUTINE.
C
C -------------- OUTPUT PARAMETERS:
C
C XF,YF,ZF - GSW COORDINATES OF THE ENDPOINT OF THE TRACED FIELD LINE.
CXX,YY,ZZ - ARRAYS OF LENGTH LMAX, CONTAINING COORDINATES OF THE FIELD LINE POINTS.
C L - ACTUAL NUMBER OF FIELD LINE POINTS, GENERATED BY THIS SUBROUTINE.
C
C ----------------------------------------------------------
C
C LAST MODIFICATION: JAN 30, 2008.
C
C common block added by B. Rideout to signal error
COMMON /BR_ERR/IERR
INTEGER IERR
C
C .. Scalar Arguments ..
DOUBLE PRECISION DIR,DSMAX,ERR,R0,RLIM,XF,XI,YF,YI,ZF,ZI
INTEGER IOPT,L,LMAX
C ..
C .. Array Arguments ..
DOUBLE PRECISION PARMOD(10),XX(LMAX),YY(LMAX),ZZ(LMAX)
C ..
C .. Subroutine Arguments ..
EXTERNAL EXNAME,INNAME
C ..
C .. Local Scalars ..
DOUBLE PRECISION AD,AL,DR,DRP,DS,FC,R,R1,R2,R3,RR,RYZ,X,XR,Y,YR,Z,
* ZR
INTEGER NREV
C ..
C .. External Subroutines ..
EXTERNAL RHAND_08,STEP_08
C ..
C .. Intrinsic Functions ..
INTRINSIC SQRT
C ..
C .. Common blocks ..
COMMON /GEOPACK1/AA,DD,BB
DOUBLE PRECISION DD
DOUBLE PRECISION AA(12),BB(21)
C ..
L = 0
NREV = 0
DD = DIR
C AUTHOR: N. A. TSYGANENKO
C
C
DS = 0.5D0*DIR
X = XI
Y = YI
Z = ZI
C initialize error flag
IERR = 0
C
C INITIALIZE THE STEP SIZE AND STARTING PONT:
C
c
chere we call RHAND_08 just to find out the sign of the radial component of the field
CALL RHAND_08(X,Y,Z,R1,R2,R3,IOPT,PARMOD,EXNAME,INNAME)
C error check added by brideout
IF (IERR .EQ. 1) THEN
L = 999
RETURN
END IF
AD = 0.01D0
IF (X*R1+Y*R2+Z*R3.LT.0.D0) AD = -0.01D0
cvector, and to determine the initial direction of the tracing (i.e., either away
c or towards Earth):
c
C
c |AD|=0.01 and its sign follows the rule:
c(1) if DIR=1 (tracing antiparallel to B vector) then the sign of AD is the same as of Br
RR = SQRT(X**2+Y**2+Z**2) + AD
c(2) if DIR=-1 (tracing parallel to B vector) then the sign of AD is opposite to that of Br
10 L = L + 1
IF (L.GT.LMAX) GO TO 40
XX(L) = X
YY(L) = Y
ZZ(L) = Z
RYZ = Y**2 + Z**2
R2 = X**2 + RYZ
R = SQRT(R2)
cAD is defined in order to initialize the value of RR (radial distance at previous step):
c
C
ccheck if the line hit the outer tracing boundary; if yes, then terminate
cthe tracing (label 8). The outer boundary is assumed reached, when the line
IF (R.GT.RLIM .OR. RYZ.GT.1600.D0 .OR. X.GT.20.D0) GO TO 50
ccrosses any of the 3 surfaces: (1) a sphere R=RLIM, (2) a cylinder of radius 40Re,
c coaxial with the XGSW axis, (3) the plane X=20Re:
c
IF (R.LT.R0 .AND. RR.GT.R) GO TO 30
ccheck whether or not the inner tracing boundary was crossed from outside,
cif yes, then calculate the footpoint position by interpolation (go to label 6):
c
IF (R.GE.RR .OR. R.GE.3.D0) GO TO 20
ccheck if we are moving outward, or R is still larger than 3Re; if yes, proceed further:
c
c
cnow we entered inside the sphere R=3: to avoid too large steps (and hence
FC = 0.2D0
IF (R-R0.LT.0.05D0) FC = 0.05D0
AL = FC*(R-R0+0.2D0)
DS = DIR*AL
cinaccurate interpolated position of the footpoint), enforce the progressively
20 XR = X
YR = Y
ZR = Z
c smaller stepsize values as we approach the inner boundary R=R0:
DRP = R - RR
RR = R
c
CALL STEP_08(X,Y,Z,DS,DSMAX,ERR,IOPT,PARMOD,EXNAME,INNAME)
C error check added by brideout
IF (IERR .EQ. 1) THEN
L = 999
RETURN
END IF
c
c
c
c
R = SQRT(X**2+Y**2+Z**2)
DR = R - RR
IF (DRP*DR.LT.0.D0) NREV = NREV + 1
IF (NREV.GT.2) GO TO 50
Ccheck the total number NREV of changes in the tracing radial direction; (NREV.GT.2) means
GO TO 10
cthat the line started making multiple loops, in which case we stop the process:
C
C
c
30 R1 = (R0-R)/(RR-R)
X = X - (X-XR)*R1
Y = Y - (Y-YR)*R1
Z = Z - (Z-ZR)*R1
GO TO 50
C 40 WRITE (*,FMT=9000)
L = LMAX
40 L = LMAX
50 XF = X
YF = Y
ZF = Z
cfind the footpoint position by interpolating between the current and previous
c field line points:
c
C
Creplace the coordinates of the last (L-th) point in the XX,YY,ZZ arrays
XX(L) = XF
YY(L) = YF
ZZ(L) = ZF
Cso that they correspond to the estimated footpoint position {XF,YF,ZF},
RETURN
9000 FORMAT (/,/,1X,'**** COMPUTATIONS IN THE SUBROUTINE TRACE ARE',
* ' TERMINATED: THE NUMBER OF POINTS EXCEEDED LMAX ****',/,/)
END
c satisfying: sqrt(XF**2+YF**2+ZF**2}=R0
C
C
SUBROUTINE SHUETAL_MGNP_08(XN_PD,VEL,BZIMF,XGSW,YGSW,ZGSW,XMGNP,
* YMGNP,ZMGNP,DIST,ID)
c
C====================================================================================
C
C
CFOR ANY POINT OF SPACE WITH COORDINATES (XGSW,YGSW,ZGSW) AND SPECIFIED CONDITIONS
C IN THE INCOMING SOLAR WIND, THIS SUBROUTINE:
C
C(1) DETERMINES IF THE POINT (XGSW,YGSW,ZGSW) LIES INSIDE OR OUTSIDE THE
C MODEL MAGNETOPAUSE OF SHUE ET AL. (JGR-A, V.103, P. 17691, 1998).
C
C(2) CALCULATES THE GSW POSITION OF A POINT {XMGNP,YMGNP,ZMGNP}, LYING AT THE MODEL
C MAGNETOPAUSE AND ASYMPTOTICALLY TENDING TO THE NEAREST BOUNDARY POINT WITH
C RESPECT TO THE OBSERVATION POINT {XGSW,YGSW,ZGSW}, AS IT APPROACHES THE MAGNETO-
C PAUSE.
C
CINPUT: XN_PD - EITHER SOLAR WIND PROTON NUMBER DENSITY (PER C.C.) (IF VEL>0)
C OR THE SOLAR WIND RAM PRESSURE IN NANOPASCALS (IF VEL<0)
C BZIMF - IMF BZ IN NANOTESLAS
C
C VEL - EITHER SOLAR WIND VELOCITY (KM/SEC)
C OR ANY NEGATIVE NUMBER, WHICH INDICATES THAT XN_PD STANDS
C FOR THE SOLAR WIND PRESSURE, RATHER THAN FOR THE DENSITY
C
C XGSW,YGSW,ZGSW - GSW POSITION OF THE OBSERVATION POINT IN EARTH RADII
C
C OUTPUT: XMGNP,YMGNP,ZMGNP - GSW POSITION OF THE BOUNDARY POINT
C DIST - DISTANCE (IN RE) BETWEEN THE OBSERVATION POINT (XGSW,YGSW,ZGSW)
C AND THE MODEL NAGNETOPAUSE
C ID - POSITION FLAG: ID=+1 (-1) MEANS THAT THE OBSERVATION POINT
C LIES INSIDE (OUTSIDE) OF THE MODEL MAGNETOPAUSE, RESPECTIVELY.
C
C OTHER SUBROUTINES USED: T96_MGNP_08
C
C .. Scalar Arguments ..
DOUBLE PRECISION BZIMF,DIST,VEL,XGSW,XMGNP,XN_PD,YGSW,YMGNP,ZGSW,
* ZMGNP
INTEGER ID
C ..
C .. Local Scalars ..
DOUBLE PRECISION ALPHA,CT,DR,DS,DT,F,GRADF,GRADF_R,GRADF_T,P,PHI,
* R,R0,RHO,RHO2,RM,ST,T,XMT96,YMT96,ZMT96
INTEGER ID96,NIT
C ..
C .. External Subroutines ..
EXTERNAL T96_MGNP_08
C ..
C .. Intrinsic Functions ..
INTRINSIC ATAN2,COS,DLOG,SIN,SQRT,TANH
C ..
IF (VEL.LT.0.D0) THEN
P = XN_PD
ELSE
P = 1.94D-6*XN_PD*VEL**2
END IF
c AUTHOR: N.A. TSYGANENKO,
C DATE: APRIL 4, 2003.
C
c
cDEFINE THE ANGLE PHI, MEASURED DUSKWARD FROM THE NOON-MIDNIGHT MERIDIAN PLANE;
IF (YGSW.NE.0.D0 .OR. ZGSW.NE.0.D0) THEN
PHI = ATAN2(YGSW,ZGSW)
ELSE
PHI = 0.D0
END IF
CIF THE OBSERVATION POINT LIES ON THE X AXIS, THE ANGLE PHI CANNOT BE UNIQUELY
C DEFINED, AND WE SET IT AT ZERO:
c
C
ID = -1
R0 = (10.22D0+1.29D0*TANH(0.184D0*(BZIMF+8.14D0)))*
* P**(-.15151515D0)
ALPHA = (0.58D0-0.007D0*BZIMF)*(1.D0+0.024D0*DLOG(P))
R = SQRT(XGSW**2+YGSW**2+ZGSW**2)
RM = R0*(2.D0/(1.D0+XGSW/R))**ALPHA
IF (R.LE.RM) ID = +1
CFIRST, FIND OUT IF THE OBSERVATION POINT LIES INSIDE THE SHUE ET AL BDRY
C AND SET THE VALUE OF THE ID FLAG:
C
C
C NOW, FIND THE CORRESPONDING T96 MAGNETOPAUSE POSITION, TO BE USED AS
CALL T96_MGNP_08(P,-1.D0,XGSW,YGSW,ZGSW,XMT96,YMT96,ZMT96,DIST,
* ID96)
C A STARTING APPROXIMATION IN THE SEARCH OF A CORRESPONDING SHUE ET AL.
RHO2 = YMT96**2 + ZMT96**2
R = SQRT(RHO2+XMT96**2)
ST = SQRT(RHO2)/R
CT = XMT96/R
C BOUNDARY POINT:
C
C
C
NIT = 0
C NOW, USE NEWTON'S ITERATIVE METHOD TO FIND THE NEAREST POINT AT THE
10 T = ATAN2(ST,CT)
RM = R0*(2.D0/(1.D0+CT))**ALPHA
C SHUE ET AL.'S BOUNDARY:
F = R - RM
GRADF_R = 1.D0
GRADF_T = -ALPHA/R*RM*ST/(1.D0+CT)
GRADF = SQRT(GRADF_R**2+GRADF_T**2)
C
DR = -F/GRADF**2
DT = DR/R*GRADF_T
R = R + DR
T = T + DT
ST = SIN(T)
CT = COS(T)
DS = SQRT(DR**2+(R*DT)**2)
NIT = NIT + 1
C IF (NIT.GT.1000) THEN
C PRINT *,
C * ' BOUNDARY POINT COULD NOT BE FOUND; ITERATIONS DO NOT CONVERGE'
C END IF
IF (DS.GT.1.D-4) GO TO 10
XMGNP = R*COS(T)
RHO = R*SIN(T)
YMGNP = RHO*SIN(PHI)
ZMGNP = RHO*COS(PHI)
DIST = SQRT((XGSW-XMGNP)**2+(YGSW-YMGNP)**2+(ZGSW-ZMGNP)**2)
RETURN
END
SUBROUTINE T96_MGNP_08(XN_PD,VEL,XGSW,YGSW,ZGSW,XMGNP,YMGNP,ZMGNP,
* DIST,ID)
C
C=======================================================================================
C
C
CFOR ANY POINT OF SPACE WITH GIVEN COORDINATES (XGSW,YGSW,ZGSW), THIS SUBROUTINE DEFINES
CTHE POSITION OF A POINT (XMGNP,YMGNP,ZMGNP) AT THE T96 MODEL MAGNETOPAUSE WITH THE
CSAME VALUE OF THE ELLIPSOIDAL TAU-COORDINATE, AND THE DISTANCE BETWEEN THEM. THIS IS
CNOT THE SHORTEST DISTANCE D_MIN TO THE BOUNDARY, BUT DIST ASYMPTOTICALLY TENDS TO D_MIN,
C AS THE OBSERVATION POINT GETS CLOSER TO THE MAGNETOPAUSE.
C
CINPUT: XN_PD - EITHER SOLAR WIND PROTON NUMBER DENSITY (PER C.C.) (IF VEL>0)
C OR THE SOLAR WIND RAM PRESSURE IN NANOPASCALS (IF VEL<0)
C VEL - EITHER SOLAR WIND VELOCITY (KM/SEC)
C OR ANY NEGATIVE NUMBER, WHICH INDICATES THAT XN_PD STANDS
C FOR THE SOLAR WIND PRESSURE, RATHER THAN FOR THE DENSITY
C
C XGSW,YGSW,ZGSW - COORDINATES OF THE OBSERVATION POINT IN EARTH RADII
C
COUTPUT: XMGNP,YMGNP,ZMGNP - GSW POSITION OF THE BOUNDARY POINT, HAVING THE SAME
C VALUE OF TAU-COORDINATE AS THE OBSERVATION POINT (XGSW,YGSW,ZGSW)
C DIST - THE DISTANCE BETWEEN THE TWO POINTS, IN RE,
C ID - POSITION FLAG; ID=+1 (-1) MEANS THAT THE POINT (XGSW,YGSW,ZGSW)
C LIES INSIDE (OUTSIDE) THE MODEL MAGNETOPAUSE, RESPECTIVELY.
C
C THE PRESSURE-DEPENDENT MAGNETOPAUSE IS THAT USED IN THE T96_01 MODEL
C (TSYGANENKO, JGR, V.100, P.5599, 1995; ESA SP-389, P.181, OCT. 1996)
C
c AUTHOR: N.A. TSYGANENKO
C DATE: AUG.1, 1995, REVISED APRIL 3, 2003.
C
C
C .. Scalar Arguments ..
DOUBLE PRECISION DIST,VEL,XGSW,XMGNP,XN_PD,YGSW,YMGNP,ZGSW,ZMGNP
INTEGER ID
C ..
C .. Local Scalars ..
DOUBLE PRECISION A,A0,ARG,PD,PHI,RAT,RAT16,RHO,RHOMGNP,S0,S00,
* SIGMA,SQ1,SQ2,TAU,X0,X00,XDZT,XKSI,XM
C ..
C .. Intrinsic Functions ..
INTRINSIC ATAN2,COS,SIN,SQRT
C ..
IF (VEL.LT.0.D0) THEN
PD = XN_PD
ELSE
PD = 1.94D-6*XN_PD*VEL**2
CDEFINE SOLAR WIND DYNAMIC PRESSURE (NANOPASCALS, ASSUMING 4% OF ALPHA-PARTICLES),
END IF
C IF NOT EXPLICITLY SPECIFIED IN THE INPUT:
C
RAT = PD/2.0D0
RAT16 = RAT**0.14D0
C
C RATIO OF PD TO THE AVERAGE PRESSURE, ASSUMED EQUAL TO 2 nPa:
C(THE POWER INDEX 0.14 IN THE SCALING FACTOR IS THE BEST-FIT VALUE OBTAINED FROM DATA
C AND USED IN THE T96_01 VERSION)
A0 = 70.D0
S00 = 1.08D0
X00 = 5.48D0
C
C VALUES OF THE MAGNETOPAUSE PARAMETERS FOR PD = 2 nPa:
C
A = A0/RAT16
S0 = S00
X0 = X00/RAT16
XM = X0 - A
C
C VALUES OF THE MAGNETOPAUSE PARAMETERS, SCALED BY THE ACTUAL PRESSURE:
C
C
C(XM IS THE X-COORDINATE OF THE "SEAM" BETWEEN THE ELLIPSOID AND THE CYLINDER)
C
C (FOR DETAILS OF THE ELLIPSOIDAL COORDINATES, SEE THE PAPER:
IF (YGSW.NE.0.D0 .OR. ZGSW.NE.0.D0) THEN
PHI = ATAN2(YGSW,ZGSW)
ELSE
PHI = 0.D0
END IF
C N.A.TSYGANENKO, SOLUTION OF CHAPMAN-FERRARO PROBLEM FOR AN
RHO = SQRT(YGSW**2+ZGSW**2)
C ELLIPSOIDAL MAGNETOPAUSE, PLANET.SPACE SCI., V.37, P.1037, 1989).
IF (XGSW.LT.XM) THEN
XMGNP = XGSW
RHOMGNP = A*SQRT(S0**2-1)
YMGNP = RHOMGNP*SIN(PHI)
ZMGNP = RHOMGNP*COS(PHI)
DIST = SQRT((XGSW-XMGNP)**2+(YGSW-YMGNP)**2+(ZGSW-ZMGNP)**2)
IF (RHOMGNP.GT.RHO) ID = +1
IF (RHOMGNP.LE.RHO) ID = -1
RETURN
END IF
C
XKSI = (XGSW-X0)/A + 1.D0
XDZT = RHO/A
SQ1 = SQRT((1.D0+XKSI)**2+XDZT**2)
SQ2 = SQRT((1.D0-XKSI)**2+XDZT**2)
SIGMA = 0.5D0*(SQ1+SQ2)
TAU = 0.5D0*(SQ1-SQ2)
C
C
C
XMGNP = X0 - A*(1.D0-S0*TAU)
ARG = (S0**2-1.D0)*(1.D0-TAU**2)
IF (ARG.LT.0.D0) ARG = 0.D0
RHOMGNP = A*SQRT(ARG)
YMGNP = RHOMGNP*SIN(PHI)
ZMGNP = RHOMGNP*COS(PHI)
C
C NOW CALCULATE (X,Y,Z) FOR THE CLOSEST POINT AT THE MAGNETOPAUSE
C
C
CNOW CALCULATE THE DISTANCE BETWEEN THE POINTS {XGSW,YGSW,ZGSW} AND {XMGNP,YMGNP,ZMGNP}:
DIST = SQRT((XGSW-XMGNP)**2+(YGSW-YMGNP)**2+(ZGSW-ZMGNP)**2)
C(IN GENERAL, THIS IS NOT THE SHORTEST DISTANCE D_MIN, BUT DIST ASYMPTOTICALLY TENDS
IF (SIGMA.GT.S0) ID = -1
IF (SIGMA.LE.S0) ID = +1
C TO D_MIN, AS WE ARE GETTING CLOSER TO THE MAGNETOPAUSE):
RETURN
END