c
c The small main program below is an example of how to compute field
c components with T89C.
c See GEOPACK.DOC for an example of field line tracing.
c
dimension parmod(10)
1 print *, ' enter x,y,z,ps,iopt'
read*, x,y,z,ps,iopt
call t89c(iopt,parmod,ps,x,y,z,bx,by,bz)
print *, bx,by,bz
goto 1
end
C
SUBROUTINE T89C(IOPT,PARMOD,PS,X,Y,Z,BX,BY,BZ)
C
C
C COMPUTES GSM COMPONENTS OF THE MAGNETIC FIELD PRODUCED BY EXTRA-
C TERRESTRIAL CURRENT SYSTEMS IN THE GEOMAGNETOSPHERE. THE MODEL IS
C VALID UP TO GEOCENTRIC DISTANCES OF 70 RE AND IS BASED ON THE MER-
C GED IMP-A,C,D,E,F,G,H,I,J (1966-1974), HEOS-1 AND -2 (1969-1974),
C AND ISEE-1 AND -2 SPACECRAFT DATA SET.
C
C THIS IS A MODIFIED VERSION (T89c), WHICH REPLACED THE ORIGINAL ONE
C IN 1992 AND DIFFERS FROM IT IN THE FOLLOWING:
C
C (1) ISEE-1,2 DATA WERE ADDED TO THE ORIGINAL IMP-HEOS DATASET
C (2) TWO TERMS WERE ADDED TO THE ORIGINAL TAIL FIELD MODES, ALLOWING
C A MODULATION OF THE CURRENT BY THE GEODIPOLE TILT ANGLE
C
C
C REFERENCE FOR THE ORIGINAL MODEL: N.A. TSYGANENKO, A MAGNETOSPHERIC MAGNETIC
C FIELD MODEL WITH A WARPED TAIL CURRENT SHEET: PLANET.SPACE SCI., V.37,
C PP.5-20, 1989.
C
C----INPUT PARAMETERS: IOPT - SPECIFIES THE GROUND DISTURBANCE LEVEL:
C
C IOPT= 1 2 3 4 5 6 7
C CORRESPOND TO:
C KP= 0,0+ 1-,1,1+ 2-,2,2+ 3-,3,3+ 4-,4,4+ 5-,5,5+ > =6-
C
C PS - GEODIPOLE TILT ANGLE IN RADIANS
C X, Y, Z - GSM COORDINATES OF THE POINT IN EARTH RADII
C
C----OUTPUT PARAMETERS: BX,BY,BZ - GSM COMPONENTS OF THE MODEL MAGNETIC
C FIELD IN NANOTESLAS
c
c THE PARAMETER PARMOD(10) IS A DUMMY ARRAY. IT IS NOT USED IN THIS
C SUBROUTINE AND IS PROVIDED JUST FOR MAKING IT COMPATIBLE WITH THE
C NEW VERSION (4/16/96) OF THE GEOPACK SOFTWARE.
C
C THIS RELEASE OF T89C IS DATED FEB 12, 1996;
C--------------------------------------------------------------------------
C
C
C AUTHOR: NIKOLAI A. TSYGANENKO
C HSTX CORP./NASA GSFC
C
DIMENSION XI(4),F(3),DER(3,30),PARAM(30,7),A(30),PARMOD(10)
DOUBLE PRECISION F,DER
DATA PARAM/-116.53,-10719.,42.375,59.753,-11363.,1.7844,30.268,
* -0.35372E-01,-0.66832E-01,0.16456E-01,-1.3024,0.16529E-02,
* 0.20293E-02,20.289,-0.25203E-01,224.91,-9234.8,22.788,7.8813,
* 1.8362,-0.27228,8.8184,2.8714,14.468,32.177,0.01,0.0,
* 7.0459,4.0,20.0,-55.553,-13198.,60.647,61.072,-16064.,
* 2.2534,34.407,-0.38887E-01,-0.94571E-01,0.27154E-01,-1.3901,
* 0.13460E-02,0.13238E-02,23.005,-0.30565E-01,55.047,-3875.7,
* 20.178,7.9693,1.4575,0.89471,9.4039,3.5215,14.474,36.555,
* 0.01,0.0,7.0787,4.0,20.0,-101.34,-13480.,111.35,12.386,-24699.,
* 2.6459,38.948,-0.34080E-01,-0.12404,0.29702E-01,-1.4052,
* 0.12103E-02,0.16381E-02,24.49,-0.37705E-01,-298.32,4400.9,18.692,
* 7.9064,1.3047,2.4541,9.7012,7.1624,14.288,33.822,0.01,0.0,6.7442,
* 4.0,20.0,-181.69,-12320.,173.79,-96.664,-39051.,3.2633,44.968,
* -0.46377E-01,-0.16686,0.048298,-1.5473,0.10277E-02,0.31632E-02,
* 27.341,-0.50655E-01,-514.10,12482.,16.257,8.5834,1.0194,3.6148,
* 8.6042,5.5057,13.778,32.373,0.01,0.0,7.3195,4.0,20.0,-436.54,
* -9001.0,323.66,-410.08,-50340.,3.9932,58.524,-0.38519E-01,
* -0.26822,0.74528E-01,-1.4268,-0.10985E-02,0.96613E-02,27.557,
* -0.56522E-01,-867.03,20652.,14.101,8.3501,0.72996,3.8149,9.2908,
* 6.4674,13.729,28.353,0.01,0.0,7.4237,4.0,20.0,-707.77,-4471.9,
* 432.81,-435.51,-60400.,4.6229,68.178,-0.88245E-01,-0.21002,
* 0.11846,-2.6711,0.22305E-02,0.10910E-01,27.547,-0.54080E-01,
* -424.23,1100.2,13.954,7.5337,0.89714,3.7813,8.2945,5.174,14.213,
* 25.237,0.01,0.0,7.0037,4.0,20.0,-1190.4,2749.9,742.56,-1110.3,
* -77193.,7.6727,102.05,-0.96015E-01,-0.74507,0.11214,-1.3614,
* 0.15157E-02,0.22283E-01,23.164,-0.74146E-01,-2219.1,48253.,
* 12.714,7.6777,0.57138,2.9633,9.3909,9.7263,11.123,21.558,0.01,
* 0.0,4.4518,4.0,20.0/
DATA IOP/10/
C
IF (IOP.NE.IOPT) THEN
C
ID=1
IOP=IOPT
DO 1 I=1,30
1 A(I)=PARAM(I,IOPT)
C
ENDIF
C
XI(1)=X
XI(2)=Y
XI(3)=Z
XI(4)=PS
CALL T89(ID,A,XI,F,DER)
IF (ID.EQ.1) ID=2
BX=F(1)
BY=F(2)
BZ=F(3)
RETURN
END
C-------------------------------------------------------------------
C
SUBROUTINE T89 (ID, A, XI, F, DER)
C
C *** N.A. Tsyganenko *** 8-10.12.1991 ***
C
C Calculates dependent model variables and their deriva-
C tives for given independent variables and model parame-
C ters. Specifies model functions with free parameters which
C must be determined by means of least squares fits (RMS
C minimization procedure).
C
C Description of parameters:
C
C ID - number of the data point in a set (initial assignments are performed
c only for ID=1, saving thus CPU time)
C A - input vector containing model parameters;
C XI - input vector containing independent variables;
C F - output double precision vector containing
C calculated values of dependent variables;
C DER - output double precision vector containing
C calculated values for derivatives of dependent
C variables with respect to model parameters;
C
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C T89 represents external magnetospheric magnetic field
C in Cartesian SOLAR MAGNETOSPHERIC coordinates (Tsyganenko N.A.,
C Planet. Space Sci., 1989, v.37, p.5-20; the "T89 model" with the warped
c tail current sheet) + A MODIFICATION ADDED IN APRIL 1992 (SEE BELOW)
C
C Model formulas for the magnetic field components contain in total
c 30 free parameters (17 linear and 13 nonlinear parameters).
C First 2 independent linear parameters A(1)-A(2) correspond to contribu-
c tion from the tail current system, then follow A(3) and A(4) which are the
c amplitudes of symmetric and antisymmetric terms in the contribution from
c the closure currents; A(5) is the ring current amplitude. Then follow the
c coefficients A(6)-A(15) which define Chapman-Ferraro+Birkeland current field.
c The coefficients c16-c19 (see Formula 20 in the original paper),
c due to DivB=0 condition, are expressed through A(6)-A(15) and hence are not
c independent ones.
c A(16) AND A(17) CORRESPOND TO THE TERMS WHICH YIELD THE TILT ANGLE DEPEN-
C DENCE OF THE TAIL CURRENT INTENSITY (ADDED ON APRIL 9, 1992)
C
C Nonlinear parameters:
C
C A(18) : DX - Characteristic scale of the Chapman-Ferraro field along the
c X-axis
C A(19) : ADR (aRC) - Characteristic radius of the ring current
c A(20) : D0 - Basic half-thickness of the tail current sheet
C A(21) : DD (GamRC)- defines rate of thickening of the ring current, as
c we go from night- to dayside
C A(22) : Rc - an analog of "hinging distance" entering formula (11)
C A(23) : G - amplitude of tail current warping in the Y-direction
C A(24) : aT - Characteristic radius of the tail current
c A(25) : Dy - characteristic scale distance in the Y direction entering
c in W(x,y) in (13)
c A(26) : Delta - defines the rate of thickening of the tail current sheet
c in the Y-direction (in T89 it was fixed at 0.01)
c A(27) : Q - this parameter was fixed at 0 in the final version of T89;
c initially it was introduced for making Dy to depend on X
c A(28) : Sx (Xo) - enters in W(x,y) ; see (13)
c A(29) : Gam (GamT) - enters in DT in (13) and defines rate of tail sheet
c thickening on going from night to dayside; in T89 fixed at 4.0
c A(30) : Dyc - the Dy parameter for closure current system; in T89 fixed
c at 20.0
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
IMPLICIT REAL * 8 (A - H, O - Z)
C
REAL A(1), XI(1)
C
DIMENSION F(3), DER(3,30)
C
INTEGER ID, I, L
DATA A02,XLW2,YN,RPI,RT/25.D0,170.D0,30.D0,0.31830989D0,30.D0/
DATA XD,XLD2/0.D0,40.D0/
C
C The last four quantities define variation of tail sheet thickness along X
C
DATA SXC,XLWC2/4.D0,50.D0/
C
C The two quantities belong to the function WC which confines tail closure
c current in X- and Y- direction
C
DATA DXL/20.D0/
C
C
IF (ID.NE.1) GOTO 3
DO 2 I = 1, 30
DO 1 L = 1, 3
1 DER(L,I) = 0.0D0
2 CONTINUE
C
DYC=A(30)
DYC2=DYC**2
DX=A(18)
HA02=0.5D0*A02
RDX2M=-1.D0/DX**2
RDX2=-RDX2M
RDYC2=1.D0/DYC2
HLWC2M=-0.5D0*XLWC2
DRDYC2=-2.D0*RDYC2
DRDYC3=2.D0*RDYC2*DSQRT(RDYC2)
HXLW2M=-0.5D0*XLW2
ADR=A(19)
D0=A(20)
DD=A(21)
RC=A(22)
G=A(23)
AT=A(24)
DT=D0
DEL=A(26)
P=A(25)
Q=A(27)
SX=A(28)
GAM=A(29)
HXLD2M=-0.5D0*XLD2
ADSL=0.D0
XGHS=0.D0
H=0.D0
HS=0.D0
GAMH=0.D0
W1=-0.5D0/DX
DBLDEL=2.D0*DEL
W2=W1*2.D0
W4=-1.D0/3.D0
W3=W4/DX
W5=-0.5D0
W6=-3.D0
AK1=A(1)
AK2=A(2)
AK3=A(3)
AK4=A(4)
AK5=A(5)
AK6=A(6)
AK7=A(7)
AK8=A(8)
AK9=A(9)
AK10=A(10)
AK11=A(11)
AK12=A(12)
AK13=A(13)
AK14=A(14)
AK15=A(15)
AK16=A(16)
AK17=A(17)
SXA=0.D0
SYA=0.D0
SZA=0.D0
AK610=AK6*W1+AK10*W5
AK711=AK7*W2-AK11
AK812=AK8*W2+AK12*W6
AK913=AK9*W3+AK13*W4
RDXL=1.D0/DXL
HRDXL=0.5D0*RDXL
A6H=AK6*0.5D0
A9T=AK9/3.D0
YNP=RPI/YN*0.5D0
YND=2.D0*YN
C
3 CONTINUE
C
X = XI(1)
Y = XI(2)
Z = XI(3)
TILT=XI(4)
TLT2=TILT**2
SPS = DSIN(TILT)
CPS = DSQRT (1.0D0 - SPS ** 2)
C
X2=X*X
Y2=Y*Y
Z2=Z*Z
TPS=SPS/CPS
HTP=TPS*0.5D0
GSP=G*SPS
XSM=X*CPS-Z*SPS
ZSM=X*SPS+Z*CPS
C
C CALCULATE THE FUNCTION ZS DEFINING THE SHAPE OF THE TAIL CURRENT SHEET
C AND ITS SPATIAL DERIVATIVES:
C
XRC=XSM+RC
XRC16=XRC**2+16.D0
SXRC=DSQRT(XRC16)
Y4=Y2*Y2
Y410=Y4+1.D4
SY4=SPS/Y410
GSY4=G*SY4
ZS1=HTP*(XRC-SXRC)
DZSX=-ZS1/SXRC
ZS=ZS1-GSY4*Y4
D2ZSGY=-SY4/Y410*4.D4*Y2*Y
DZSY=G*D2ZSGY
C
C CALCULATE THE COMPONENTS OF THE RING CURRENT CONTRIBUTION:
C
XSM2=XSM**2
DSQT=DSQRT(XSM2+A02)
FA0=0.5D0*(1.D0+XSM/DSQT)
DDR=D0+DD*FA0
DFA0=HA02/DSQT**3
ZR=ZSM-ZS
TR=DSQRT(ZR**2+DDR**2)
RTR=1.D0/TR
RO2=XSM2+Y2
ADRT=ADR+TR
ADRT2=ADRT**2
FK=1.D0/(ADRT2+RO2)
DSFC=DSQRT(FK)
FC=FK**2*DSFC
FACXY=3.0D0*ADRT*FC*RTR
XZR=XSM*ZR
YZR=Y*ZR
DBXDP=FACXY*XZR
DER(2,5)=FACXY*YZR
XZYZ=XSM*DZSX+Y*DZSY
FAQ=ZR*XZYZ-DDR*DD*DFA0*XSM
DBZDP=FC*(2.D0*ADRT2-RO2)+FACXY*FAQ
DER(1,5)=DBXDP*CPS+DBZDP*SPS
DER(3,5)=DBZDP*CPS-DBXDP*SPS
C
C CALCULATE THE TAIL CURRENT SHEET CONTRIBUTION:
C
DELY2=DEL*Y2
D=DT+DELY2
IF (DABS(GAM).LT.1.D-6) GOTO 8
XXD=XSM-XD
RQD=1.D0/(XXD**2+XLD2)
RQDS=DSQRT(RQD)
H=0.5D0*(1.D0+XXD*RQDS)
HS=-HXLD2M*RQD*RQDS
GAMH=GAM*H
D=D+GAMH
XGHS=XSM*GAM*HS
ADSL=-D*XGHS
8 D2=D**2
T=DSQRT(ZR**2+D2)
XSMX=XSM-SX
RDSQ2=1.D0/(XSMX**2+XLW2)
RDSQ=DSQRT(RDSQ2)
V=0.5D0*(1.D0-XSMX*RDSQ)
DVX=HXLW2M*RDSQ*RDSQ2
OM=DSQRT(DSQRT(XSM2+16.D0)-XSM)
OMS=-OM/(OM*OM+XSM)*0.5D0
RDY=1.D0/(P+Q*OM)
OMSV=OMS*V
RDY2=RDY**2
FY=1.D0/(1.D0+Y2*RDY2)
W=V*FY
YFY1=2.D0*FY*Y2*RDY2
FYPR=YFY1*RDY
FYDY=FYPR*FY
DWX=DVX*FY+FYDY*Q*OMSV
YDWY=-V*YFY1*FY
DDY=DBLDEL*Y
ATT=AT+T
S1=DSQRT(ATT**2+RO2)
F5=1.D0/S1
F7=1.D0/(S1+ATT)
F1=F5*F7
F3=F5**3
F9=ATT*F3
FS=ZR*XZYZ-D*Y*DDY+ADSL
XDWX=XSM*DWX+YDWY
RTT=1.D0/T
WT=W*RTT
BRRZ1=WT*F1
BRRZ2=WT*F3
DBXC1=BRRZ1*XZR
DBXC2=BRRZ2*XZR
DER(2,1)=BRRZ1*YZR
DER(2,2)=BRRZ2*YZR
DER(2,16)=DER(2,1)*TLT2
DER(2,17)=DER(2,2)*TLT2
WTFS=WT*FS
DBZC1=W*F5+XDWX*F7+WTFS*F1
DBZC2=W*F9+XDWX*F1+WTFS*F3
DER(1,1)=DBXC1*CPS+DBZC1*SPS
DER(1,2)=DBXC2*CPS+DBZC2*SPS
DER(3,1)=DBZC1*CPS-DBXC1*SPS
DER(3,2)=DBZC2*CPS-DBXC2*SPS
DER(1,16)=DER(1,1)*TLT2
DER(1,17)=DER(1,2)*TLT2
DER(3,16)=DER(3,1)*TLT2
DER(3,17)=DER(3,2)*TLT2
C
C CALCULATE CONTRIBUTION FROM THE CLOSURE CURRENTS
C
ZPL=Z+RT
ZMN=Z-RT
ROGSM2=X2+Y2
SPL=DSQRT(ZPL**2+ROGSM2)
SMN=DSQRT(ZMN**2+ROGSM2)
XSXC=X-SXC
RQC2=1.D0/(XSXC**2+XLWC2)
RQC=DSQRT(RQC2)
FYC=1.D0/(1.D0+Y2*RDYC2)
WC=0.5D0*(1.D0-XSXC*RQC)*FYC
DWCX=HLWC2M*RQC2*RQC*FYC
DWCY=DRDYC2*WC*FYC*Y
SZRP=1.D0/(SPL+ZPL)
SZRM=1.D0/(SMN-ZMN)
XYWC=X*DWCX+Y*DWCY
WCSP=WC/SPL
WCSM=WC/SMN
FXYP=WCSP*SZRP
FXYM=WCSM*SZRM
FXPL=X*FXYP
FXMN=-X*FXYM
FYPL=Y*FXYP
FYMN=-Y*FXYM
FZPL=WCSP+XYWC*SZRP
FZMN=WCSM+XYWC*SZRM
DER(1,3)=FXPL+FXMN
DER(1,4)=(FXPL-FXMN)*SPS
DER(2,3)=FYPL+FYMN
DER(2,4)=(FYPL-FYMN)*SPS
DER(3,3)=FZPL+FZMN
DER(3,4)=(FZPL-FZMN)*SPS
C
C NOW CALCULATE CONTRIBUTION FROM CHAPMAN-FERRARO SOURCES + ALL OTHER
C
EX=DEXP(X/DX)
EC=EX*CPS
ES=EX*SPS
ECZ=EC*Z
ESZ=ES*Z
ESZY2=ESZ*Y2
ESZZ2=ESZ*Z2
ECZ2=ECZ*Z
ESY=ES*Y
C
DER(1,6)=ECZ
DER(1,7)=ES
DER(1,8)=ESY*Y
DER(1,9)=ESZ*Z
DER(2,10)=ECZ*Y
DER(2,11)=ESY
DER(2,12)=ESY*Y2
DER(2,13)=ESY*Z2
DER(3,14)=EC
DER(3,15)=EC*Y2
DER(3,6)=ECZ2*W1
DER(3,10)=ECZ2*W5
DER(3,7)=ESZ*W2
DER(3,11)=-ESZ
DER(3,8)=ESZY2*W2
DER(3,12)=ESZY2*W6
DER(3,9)=ESZZ2*W3
DER(3,13)=ESZZ2*W4
C
C FINALLY, CALCULATE NET EXTERNAL MAGNETIC FIELD COMPONENTS,
C BUT FIRST OF ALL THOSE FOR C.-F. FIELD:
C
SX1=AK6*DER(1,6)+AK7*DER(1,7)+AK8*DER(1,8)+AK9*DER(1,9)
SY1=AK10*DER(2,10)+AK11*DER(2,11)+AK12*DER(2,12)+AK13*DER(2,13)
SZ1=AK14*DER(3,14)+AK15*DER(3,15)+AK610*ECZ2+AK711*ESZ+AK812
* *ESZY2+AK913*ESZZ2
BXCL=AK3*DER(1,3)+AK4*DER(1,4)
BYCL=AK3*DER(2,3)+AK4*DER(2,4)
BZCL=AK3*DER(3,3)+AK4*DER(3,4)
BXT=AK1*DER(1,1)+AK2*DER(1,2)+BXCL +AK16*DER(1,16)+AK17*DER(1,17)
BYT=AK1*DER(2,1)+AK2*DER(2,2)+BYCL +AK16*DER(2,16)+AK17*DER(2,17)
BZT=AK1*DER(3,1)+AK2*DER(3,2)+BZCL +AK16*DER(3,16)+AK17*DER(3,17)
F(1)=BXT+AK5*DER(1,5)+SX1+SXA
F(2)=BYT+AK5*DER(2,5)+SY1+SYA
F(3)=BZT+AK5*DER(3,5)+SZ1+SZA
C
RETURN
END
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%