|
|
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
|
|
|
|
