##// END OF EJS Templates
DATABASES » madrigal
RSS
Review updata app
Review updata app
jespinoza - - Load All Authors

File last commit:

r0:b84e1135c2c4
r8:9c1b6919a8a9
Show More
TS04.f
2562 lines | 83.4 KiB | text/x-fortran | FortranFixedLexer
/ source / madf / geolib / TS04.f
History | Source | Raw |Copy content |Copy permalink
jespinoza
Initial
 r0 c
c====================================================================================
c
c
SUBROUTINE T04_s (IOPT,PARMOD,PS,X,Y,Z,BX,BY,BZ)
c
c ASSEMBLED: MARCH 25, 2004;
C UPDATED: AUGUST 2 & 31, DECEMBER 27, 2004.
c LATEST MODIFICATIONS/BUGS REMOVED:
c
C (1) MARCH 14, 2005: 79 -> 69 (LINE 94; might cause compilation problems with some Fortran compilers)
c
C (2) JUNE 24, 2006: REPLACED COEFFICIENTS IN
c (i) DATA statement in FUNCTION AP,
C (ii) DATA C_SY statement in SUBROUTINE FULL_RC, and
c (iii) DATA A statement in SUBROUTINE T04_s.
C This correction was needed due to a bug found in the symmetric ring current module.
c Its impact can be significant (up to ~20 nT) only in the innermost magnetosphere (R<=2)
C and only for strongly disturbed conditions; otherwise, the change in the model field
c does not exceed a few percent.
c
c--------------------------------------------------------------------
C A DATA-BASED MODEL OF THE EXTERNAL (I.E., WITHOUT EARTH'S CONTRIBUTION) PART OF THE
C MAGNETOSPHERIC MAGNETIC FIELD, CALIBRATED BY
C (1) SOLAR WIND PRESSURE PDYN (NANOPASCALS),
C (2) DST (NANOTESLA),
C (3) BYIMF,
C (4) BZIMF (NANOTESLA)
C (5-10) INDICES W1 - W6, CALCULATED AS TIME INTEGRALS FROM THE BEGINNING OF A STORM
c SEE THE REFERENCE (3) BELOW, FOR A DETAILED DEFINITION OF THOSE VARIABLES
C
c THE ABOVE 10 INPUT PARAMETERS SHOULD BE PLACED IN THE ELEMENTS
c OF THE ARRAY PARMOD(10).
C
C THE REST OF THE INPUT VARIABLES ARE: THE GEODIPOLE TILT ANGLE PS (RADIANS),
C X,Y,Z - GSM POSITION (RE)
C
c IOPT IS A DUMMY INPUT PARAMETER, INCLUDED TO MAKE THIS SUBROUTINE
C COMPATIBLE WITH THE TRACING SOFTWARE PACKAGE (GEOPACK). IN THIS MODEL,
C THE PARAMETER IOPT DOES NOT AFFECT THE OUTPUT FIELD.
c
C*******************************************************************************************
c** ATTENTION: THE MODEL IS BASED ON DATA TAKEN SUNWARD FROM X=-15Re, AND HENCE BECOMES *
C** INVALID AT LARGER TAILWARD DISTANCES !!! *
C*******************************************************************************************
C
c OUTPUT: GSM COMPONENTS OF THE EXTERNAL MAGNETIC FIELD (BX,BY,BZ, nanotesla)
C COMPUTED AS A SUM OF CONTRIBUTIONS FROM PRINCIPAL FIELD SOURCES
C
c (C) Copr. 2004, Nikolai A. Tsyganenko, USRA/Code 612.3, NASA GSFC
c Greenbelt, MD 20771, USA
c
C REFERENCES:
C
C (1) N. A. Tsyganenko, A new data-based model of the near magnetosphere magnetic field:
c 1. Mathematical structure.
c 2. Parameterization and fitting to observations. JGR v. 107(A8), 1176/1179, doi:10.1029/2001JA000219/220, 2002.
c
c (2) N. A. Tsyganenko, H. J. Singer, J. C. Kasper, Storm-time distortion of the
c inner magnetosphere: How severe can it get ? JGR v. 108(A5), 1209, doi:10.1029/2002JA009808, 2003.
c (3) N. A. Tsyganenko and M. I. Sitnov, Modeling the dynamics of the inner magnetosphere during
c strong geomagnetic storms, J. Geophys. Res., v. 110 (A3), A03208, doi: 10.1029/2004JA010798, 2005.
c----------------------------------------------------------------------
c
DOUBLE PRECISION PARMOD(10),PS,X,Y,Z,BX,BY,BZ
REAL*8 A(69),PDYN,DST_AST,BXIMF,BYIMF,BZIMF,W1,W2,W3,W4,W5,W6,
* PSS,XX,YY,ZZ,BXCF,BYCF,BZCF,BXT1,BYT1,BZT1,BXT2,BYT2,BZT2,
* BXSRC,BYSRC,BZSRC,BXPRC,BYPRC,BZPRC, BXR11,BYR11,BZR11,
* BXR12,BYR12,BZR12,BXR21,BYR21,BZR21,BXR22,BYR22,BZR22,HXIMF,
* HYIMF,HZIMF,BBX,BBY,BBZ
C
DATA A/1.00000,5.44118,0.891995,9.09684,0.00000,-7.18972,12.2700,
* -4.89408,0.00000,0.870536,1.36081,0.00000,0.688650,0.602330,
* 0.00000,0.316346,1.22728,-0.363620E-01,-0.405821,0.452536,
* 0.755831,0.215662,0.152759,5.96235,23.2036,11.2994,69.9596,
* 0.989596,-0.132131E-01,0.985681,0.344212E-01,1.02389,0.207867,
* 1.51220,0.682715E-01,1.84714,1.76977,1.37690,0.696350,0.343280,
* 3.28846,111.293,5.82287,4.39664,0.383403,0.648176,0.318752E-01,
* 0.581168,1.15070,0.843004,0.394732,0.846509,0.916555,0.550920,
* 0.180725,0.898772,0.387365,2.26596,1.29123,0.436819,1.28211,
* 1.33199,.405553,1.6229,.699074,1.26131,2.42297,.537116,.619441/
c
DATA IOPGEN,IOPTT,IOPB,IOPR/0,0,0,0/
C
PDYN=PARMOD(1)
DST_AST=PARMOD(2)*0.8-13.*SQRT(PDYN)
BYIMF=PARMOD(3)
BZIMF=PARMOD(4)
C
W1=PARMOD (5)
W2=PARMOD (6)
W3=PARMOD (7)
W4=PARMOD (8)
W5=PARMOD (9)
W6=PARMOD(10)
PSS=PS
XX=X
YY=Y
ZZ=Z
C
CALL EXTERN (IOPGEN,IOPTT,IOPB,IOPR,A,69,PDYN,DST_AST,BXIMF,BYIMF,
+ BZIMF,W1,W2,W3,W4,W5,W6,PSS,XX,YY,ZZ,BXCF,BYCF,BZCF,BXT1,BYT1,
+ BZT1,BXT2,BYT2,BZT2,BXSRC,BYSRC,BZSRC,BXPRC,BYPRC,BZPRC, BXR11,
+ BYR11,BZR11,BXR12,BYR12,BZR12,BXR21,BYR21,BZR21,BXR22,BYR22,
+ BZR22,HXIMF,HYIMF,HZIMF,BBX,BBY,BBZ)
C
BX=BBX
BY=BBY
BZ=BBZ
C
RETURN
END
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
SUBROUTINE EXTERN (IOPGEN,IOPT,IOPB,IOPR,A,NTOT,
* PDYN,DST,BXIMF,BYIMF,BZIMF,W1,W2,W3,W4,W5,W6,PS,X,Y,Z,
* BXCF,BYCF,BZCF,BXT1,BYT1,BZT1,BXT2,BYT2,BZT2,
* BXSRC,BYSRC,BZSRC,BXPRC,BYPRC,BZPRC, BXR11,BYR11,BZR11,
* BXR12,BYR12,BZR12,BXR21,BYR21,BZR21,BXR22,BYR22,BZR22,HXIMF,
* HYIMF,HZIMF,BX,BY,BZ)
C
C IOPGEN - GENERAL OPTION FLAG: IOPGEN=0 - CALCULATE TOTAL FIELD
C IOPGEN=1 - DIPOLE SHIELDING ONLY
C IOPGEN=2 - TAIL FIELD ONLY
C IOPGEN=3 - BIRKELAND FIELD ONLY
C IOPGEN=4 - RING CURRENT FIELD ONLY
C IOPGEN=5 - INTERCONNECTION FIELD ONLY
C
C IOPT - TAIL FIELD FLAG: IOPT=0 - BOTH MODES
C IOPT=1 - MODE 1 ONLY
C IOPT=2 - MODE 2 ONLY
C
C IOPB - BIRKELAND FIELD FLAG: IOPB=0 - ALL 4 TERMS
C IOPB=1 - REGION 1, MODES 1 AND 2
C IOPB=2 - REGION 2, MODES 1 AND 2
C
C IOPR - RING CURRENT FLAG: IOPR=0 - BOTH SRC AND PRC
C IOPR=1 - SRC ONLY
C IOPR=2 - PRC ONLY
C
IMPLICIT REAL * 8 (A - H, O - Z)
C
DIMENSION A(NTOT)
C
COMMON /TAIL/ DXSHIFT1,DXSHIFT2,D,DELTADY ! THE COMMON BLOCKS FORWARD NONLINEAR PARAMETERS
COMMON /BIRKPAR/ XKAPPA1,XKAPPA2
COMMON /RCPAR/ SC_SY,SC_AS,PHI
COMMON /G/ G
COMMON /RH0/ RH0
C common block added by B. Rideout to signal error
COMMON /BR_ERR/IERR
INTEGER IERR
INTEGER LOOPCNT
C
C
DATA A0_A,A0_S0,A0_X0 /34.586D0,1.1960D0,3.4397D0/ ! SHUE ET AL. PARAMETERS
DATA DSIG /0.005D0/, RH0,RH2 /7.5D0,-5.2D0/
c
XAPPA=(PDYN/2.)**A(23) ! OVERALL SCALING PARAMETER
RH0=7.5 ! TAIL HINGING DISTANCE
c
G= 35.0 ! TAIL WARPING PARAMETER
LOOPCNT = 0
XAPPA3=XAPPA**3
XX=X*XAPPA
YY=Y*XAPPA
ZZ=Z*XAPPA
C
SPS=DSIN(PS)
c
X0=A0_X0/XAPPA
AM=A0_A/XAPPA
S0=A0_S0
c
C CALCULATE "IMF" COMPONENTS OUTSIDE THE MAGNETOPAUSE LAYER (HENCE BEGIN WITH "O")
C THEY ARE NEEDED ONLY IF THE POINT (X,Y,Z) IS WITHIN THE TRANSITION MAGNETOPAUSE LAYER
C OR OUTSIDE THE MAGNETOSPHERE:
C
FACTIMF=A(20)
c
OIMFX=0.D0
OIMFY=BYIMF*FACTIMF
OIMFZ=BZIMF*FACTIMF
c
R=DSQRT(X**2+Y**2+Z**2)
XSS=X
ZSS=Z
1 XSOLD=XSS ! BEGIN ITERATIVE SEARCH OF UNWARPED COORDS (TO FIND SIGMA)
ZSOLD=ZSS
RH=RH0+RH2*(ZSS/R)**2
SINPSAS=SPS/(1.D0+(R/RH)**3)**0.33333333D0
COSPSAS=DSQRT(1.D0-SINPSAS**2)
ZSS=X*SINPSAS+Z*COSPSAS
XSS=X*COSPSAS-Z*SINPSAS
DD=DABS(XSS-XSOLD)+DABS(ZSS-ZSOLD)
C error check added by B. Rideout
LOOPCNT = LOOPCNT + 1
IF (LOOPCNT.GT.500) THEN
IERR = 1
C print*,'loopcnt error at ',DD
RETURN
ENDIF
IF (DD.GT.1.D-6) GOTO 1
C END OF ITERATIVE SEARCH
C print*,'loopcnt success at ',DD
RHO2=Y**2+ZSS**2
ASQ=AM**2
XMXM=AM+XSS-X0
IF (XMXM.LT.0.) XMXM=0. ! THE BOUNDARY IS A CYLINDER TAILWARD OF X=X0-AM
AXX0=XMXM**2
ARO=ASQ+RHO2
SIGMA=DSQRT((ARO+AXX0+SQRT((ARO+AXX0)**2-4.*ASQ*AXX0))/(2.*ASQ))
C
C NOW, THERE ARE THREE POSSIBLE CASES:
C (1) INSIDE THE MAGNETOSPHERE (SIGMA
C (2) IN THE BOUNDARY LAYER
C (3) OUTSIDE THE MAGNETOSPHERE AND B.LAYER
C FIRST OF ALL, CONSIDER THE CASES (1) AND (2):
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
IF (SIGMA.LT.S0+DSIG) THEN ! CASES (1) OR (2); CALCULATE THE MODEL FIELD
C (WITH THE POTENTIAL "PENETRATED" INTERCONNECTION FIELD):
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
IF (IOPGEN.LE.1) THEN
CALL SHLCAR3X3(XX,YY,ZZ,PS,CFX,CFY,CFZ) ! DIPOLE SHIELDING FIELD
BXCF=CFX*XAPPA3
BYCF=CFY*XAPPA3
BZCF=CFZ*XAPPA3
ELSE
BXCF=0.D0
BYCF=0.D0
BZCF=0.D0
ENDIF ! DONE
IF (IOPGEN.EQ.0.OR.IOPGEN.EQ.2) THEN
DSTT=-20.
IF (DST.LT.DSTT) DSTT=DST
ZNAM=DABS(DSTT)**0.37
DXSHIFT1=A(24)-A(25)/ZNAM
DXSHIFT2=A(26)-A(27)/ZNAM
D=A(36)*DEXP(-W1/A(37)) +A(69)
DELTADY=4.7
CALL DEFORMED (IOPT,PS,XX,YY,ZZ, ! TAIL FIELD (THREE MODES)
* BXT1,BYT1,BZT1,BXT2,BYT2,BZT2)
ELSE
BXT1=0.D0
BYT1=0.D0
BZT1=0.D0
BXT2=0.D0
BYT2=0.D0
BZT2=0.D0
ENDIF
IF (IOPGEN.EQ.0.OR.IOPGEN.EQ.3) THEN
ZNAM=DABS(DST)
IF (DST.GE.-20.D0) ZNAM=20.D0
XKAPPA1=A(32)*(ZNAM/20.D0)**A(33)
XKAPPA2=A(34)*(ZNAM/20.D0)**A(35)
CALL BIRK_TOT (IOPB,PS,XX,YY,ZZ,BXR11,BYR11,BZR11,BXR12,BYR12,
* BZR12,BXR21,BYR21,BZR21,BXR22,BYR22,BZR22) ! BIRKELAND FIELD (TWO MODES FOR R1 AND TWO MODES FOR R2)
ELSE
BXR11=0.D0
BYR11=0.D0
BZR11=0.D0
BXR21=0.D0
BYR21=0.D0
BZR21=0.D0
ENDIF
IF (IOPGEN.EQ.0.OR.IOPGEN.EQ.4) THEN
PHI=A(38)
ZNAM=DABS(DST)
IF (DST.GE.-20.D0) ZNAM=20.D0
SC_SY=A(28)*(20.D0/ZNAM)**A(29) *XAPPA !
SC_AS=A(30)*(20.D0/ZNAM)**A(31) *XAPPA ! MULTIPLICATION BY XAPPA IS MADE IN ORDER TO MAKE THE SRC AND PRC
! SCALING COMPLETELY INDEPENDENT OF THE GENERAL SCALING DUE TO THE
C MAGNETOPAUSE COMPRESSION/EXPANSION !
C
CALL FULL_RC(IOPR,PS,XX,YY,ZZ,BXSRC,BYSRC,BZSRC,BXPRC,BYPRC,
* BZPRC) ! SHIELDED RING CURRENT (SRC AND PRC)
ELSE
BXSRC=0.D0
BYSRC=0.D0
BZSRC=0.D0
BXPRC=0.D0
BYPRC=0.D0
BZPRC=0.D0
ENDIF
C
IF (IOPGEN.EQ.0.OR.IOPGEN.EQ.5) THEN
HXIMF=0.D0
HYIMF=BYIMF
HZIMF=BZIMF ! THESE ARE COMPONENTS OF THE PENETRATED FIELD PER UNIT OF THE PENETRATION COEFFICIENT.
C IN OTHER WORDS, THESE ARE DERIVATIVES OF THE PENETRATION FIELD COMPONENTS WITH RESPECT
C TO THE PENETRATION COEFFICIENT. WE ASSUME THAT ONLY TRANSVERSE COMPONENT OF THE
C FIELD PENETRATES INSIDE.
ELSE
HXIMF=0.D0
HYIMF=0.D0
HZIMF=0.D0
ENDIF
C
C-----------------------------------------------------------
C
C NOW, ADD UP ALL THE COMPONENTS:
c
DLP1=(PDYN/2.D0)**A(21)
DLP2=(PDYN/2.D0)**A(22)
TAMP1=A(2)+A(3)*DLP1+A(4)*A(39)*W1/DSQRT(W1**2+A(39)**2)+A(5)*DST
TAMP2=A(6)+A(7)*DLP2+A(8)*A(40)*W2/DSQRT(W2**2+A(40)**2)+A(9)*DST
A_SRC=A(10)+A(11)*A(41)*W3/DSQRT(W3**2+A(41)**2)
* +A(12)*DST
A_PRC=A(13)+A(14)*A(42)*W4/DSQRT(W4**2+A(42)**2)
* +A(15)*DST
A_R11=A(16)+A(17)*A(43)*W5/DSQRT(W5**2+A(43)**2)
A_R21=A(18)+A(19)*A(44)*W6/DSQRT(W6**2+A(44)**2)
BBX=A(1)*BXCF+TAMP1*BXT1+TAMP2*BXT2+A_SRC*BXSRC+A_PRC*BXPRC
* +A_R11*BXR11+A_R21*BXR21+A(20)*HXIMF
BBY=A(1)*BYCF+TAMP1*BYT1+TAMP2*BYT2+A_SRC*BYSRC+A_PRC*BYPRC
* +A_R11*BYR11+A_R21*BYR21+A(20)*HYIMF
BBZ=A(1)*BZCF+TAMP1*BZT1+TAMP2*BZT2+A_SRC*BZSRC+A_PRC*BZPRC
* +A_R11*BZR11+A_R21*BZR21+A(20)*HZIMF
C
C AND WE HAVE THE TOTAL EXTERNAL FIELD.
C
C
C NOW, LET US CHECK WHETHER WE HAVE THE CASE (1). IF YES - ALL DONE:
C
IF (SIGMA.LT.S0-DSIG) THEN ! (X,Y,Z) IS INSIDE THE MAGNETOSPHERE
BX=BBX
BY=BBY
BZ=BBZ
ELSE ! THIS IS THE MOST COMPLEX CASE: WE ARE INSIDE
C THE INTERPOLATION REGION
FINT=0.5*(1.-(SIGMA-S0)/DSIG)
FEXT=0.5*(1.+(SIGMA-S0)/DSIG)
C
CALL DIPOLE (PS,X,Y,Z,QX,QY,QZ)
BX=(BBX+QX)*FINT+OIMFX*FEXT -QX
BY=(BBY+QY)*FINT+OIMFY*FEXT -QY
BZ=(BBZ+QZ)*FINT+OIMFZ*FEXT -QZ
c
ENDIF ! THE CASES (1) AND (2) ARE EXHAUSTED; THE ONLY REMAINING
C POSSIBILITY IS NOW THE CASE (3):
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
ELSE
CALL DIPOLE (PS,X,Y,Z,QX,QY,QZ)
BX=OIMFX-QX
BY=OIMFY-QY
BZ=OIMFZ-QZ
ENDIF
C
END
c
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
SUBROUTINE SHLCAR3X3(X,Y,Z,PS,BX,BY,BZ)
C
C THIS S/R RETURNS THE SHIELDING FIELD FOR THE EARTH'S DIPOLE,
C REPRESENTED BY 2x3x3=18 "CARTESIAN" HARMONICS, tilted with respect
C to the z=0 plane
C
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C The 36 coefficients enter in pairs in the amplitudes of the "cartesian"
c harmonics (A(1)-A(36).
c The 14 nonlinear parameters (A(37)-A(50) are the scales Pi,Ri,Qi,and Si
C entering the arguments of exponents, sines, and cosines in each of the
C 18 "Cartesian" harmonics PLUS TWO TILT ANGLES FOR THE CARTESIAN HARMONICS
C (ONE FOR THE PSI=0 MODE AND ANOTHER FOR THE PSI=90 MODE)
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
IMPLICIT REAL * 8 (A - H, O - Z)
C
DIMENSION A(50)
DATA A/-901.2327248,895.8011176,817.6208321,-845.5880889,
*-83.73539535,86.58542841,336.8781402,-329.3619944,-311.2947120,
*308.6011161,31.94469304,-31.30824526,125.8739681,-372.3384278,
*-235.4720434,286.7594095,21.86305585,-27.42344605,-150.4874688,
*2.669338538,1.395023949,-.5540427503,-56.85224007,3.681827033,
*-43.48705106,5.103131905,1.073551279,-.6673083508,12.21404266,
*4.177465543,5.799964188,-.3977802319,-1.044652977,.5703560010,
*3.536082962,-3.222069852,9.620648151,6.082014949,27.75216226,
*12.44199571,5.122226936,6.982039615,20.12149582,6.150973118,
*4.663639687,15.73319647,2.303504968,5.840511214,.8385953499E-01,
*.3477844929/
C
P1=A(37)
P2=A(38)
P3=A(39)
R1=A(40)
R2=A(41)
R3=A(42)
Q1=A(43)
Q2=A(44)
Q3=A(45)
S1=A(46)
S2=A(47)
S3=A(48)
T1 =A(49)
T2 =A(50)
C
CPS=DCOS(PS)
SPS=DSIN(PS)
S2PS=2.D0*CPS
C
ST1=DSIN(PS*T1)
CT1=DCOS(PS*T1)
ST2=DSIN(PS*T2)
CT2=DCOS(PS*T2)
X1=X*CT1-Z*ST1
Z1=X*ST1+Z*CT1
X2=X*CT2-Z*ST2
Z2=X*ST2+Z*CT2
C
C
c MAKE THE TERMS IN THE 1ST SUM ("PERPENDICULAR" SYMMETRY):
C
C I=1
C
SQPR= DSQRT(1.D0/P1**2+1.D0/R1**2)
CYP = DCOS(Y/P1)
SYP = DSIN(Y/P1)
CZR = DCOS(Z1/R1)
SZR = DSIN(Z1/R1)
EXPR= DEXP(SQPR*X1)
FX1 =-SQPR*EXPR*CYP*SZR
HY1 = EXPR/P1*SYP*SZR
FZ1 =-EXPR*CYP/R1*CZR
HX1 = FX1*CT1+FZ1*ST1
HZ1 =-FX1*ST1+FZ1*CT1
SQPR= DSQRT(1.D0/P1**2+1.D0/R2**2)
CYP = DCOS(Y/P1)
SYP = DSIN(Y/P1)
CZR = DCOS(Z1/R2)
SZR = DSIN(Z1/R2)
EXPR= DEXP(SQPR*X1)
FX2 =-SQPR*EXPR*CYP*SZR
HY2 = EXPR/P1*SYP*SZR
FZ2 =-EXPR*CYP/R2*CZR
HX2 = FX2*CT1+FZ2*ST1
HZ2 =-FX2*ST1+FZ2*CT1
SQPR= DSQRT(1.D0/P1**2+1.D0/R3**2)
CYP = DCOS(Y/P1)
SYP = DSIN(Y/P1)
CZR = DCOS(Z1/R3)
SZR = DSIN(Z1/R3)
EXPR= DEXP(SQPR*X1)
FX3 =-EXPR*CYP*(SQPR*Z1*CZR+SZR/R3*(X1+1.D0/SQPR))
HY3 = EXPR/P1*SYP*(Z1*CZR+X1/R3*SZR/SQPR)
FZ3 =-EXPR*CYP*(CZR*(1.D0+X1/R3**2/SQPR)-Z1/R3*SZR)
HX3 = FX3*CT1+FZ3*ST1
HZ3 =-FX3*ST1+FZ3*CT1
C
C I=2:
C
SQPR= DSQRT(1.D0/P2**2+1.D0/R1**2)
CYP = DCOS(Y/P2)
SYP = DSIN(Y/P2)
CZR = DCOS(Z1/R1)
SZR = DSIN(Z1/R1)
EXPR= DEXP(SQPR*X1)
FX4 =-SQPR*EXPR*CYP*SZR
HY4 = EXPR/P2*SYP*SZR
FZ4 =-EXPR*CYP/R1*CZR
HX4 = FX4*CT1+FZ4*ST1
HZ4 =-FX4*ST1+FZ4*CT1
SQPR= DSQRT(1.D0/P2**2+1.D0/R2**2)
CYP = DCOS(Y/P2)
SYP = DSIN(Y/P2)
CZR = DCOS(Z1/R2)
SZR = DSIN(Z1/R2)
EXPR= DEXP(SQPR*X1)
FX5 =-SQPR*EXPR*CYP*SZR
HY5 = EXPR/P2*SYP*SZR
FZ5 =-EXPR*CYP/R2*CZR
HX5 = FX5*CT1+FZ5*ST1
HZ5 =-FX5*ST1+FZ5*CT1
SQPR= DSQRT(1.D0/P2**2+1.D0/R3**2)
CYP = DCOS(Y/P2)
SYP = DSIN(Y/P2)
CZR = DCOS(Z1/R3)
SZR = DSIN(Z1/R3)
EXPR= DEXP(SQPR*X1)
FX6 =-EXPR*CYP*(SQPR*Z1*CZR+SZR/R3*(X1+1.D0/SQPR))
HY6 = EXPR/P2*SYP*(Z1*CZR+X1/R3*SZR/SQPR)
FZ6 =-EXPR*CYP*(CZR*(1.D0+X1/R3**2/SQPR)-Z1/R3*SZR)
HX6 = FX6*CT1+FZ6*ST1
HZ6 =-FX6*ST1+FZ6*CT1
C
C I=3:
C
SQPR= DSQRT(1.D0/P3**2+1.D0/R1**2)
CYP = DCOS(Y/P3)
SYP = DSIN(Y/P3)
CZR = DCOS(Z1/R1)
SZR = DSIN(Z1/R1)
EXPR= DEXP(SQPR*X1)
FX7 =-SQPR*EXPR*CYP*SZR
HY7 = EXPR/P3*SYP*SZR
FZ7 =-EXPR*CYP/R1*CZR
HX7 = FX7*CT1+FZ7*ST1
HZ7 =-FX7*ST1+FZ7*CT1
SQPR= DSQRT(1.D0/P3**2+1.D0/R2**2)
CYP = DCOS(Y/P3)
SYP = DSIN(Y/P3)
CZR = DCOS(Z1/R2)
SZR = DSIN(Z1/R2)
EXPR= DEXP(SQPR*X1)
FX8 =-SQPR*EXPR*CYP*SZR
HY8 = EXPR/P3*SYP*SZR
FZ8 =-EXPR*CYP/R2*CZR
HX8 = FX8*CT1+FZ8*ST1
HZ8 =-FX8*ST1+FZ8*CT1
SQPR= DSQRT(1.D0/P3**2+1.D0/R3**2)
CYP = DCOS(Y/P3)
SYP = DSIN(Y/P3)
CZR = DCOS(Z1/R3)
SZR = DSIN(Z1/R3)
EXPR= DEXP(SQPR*X1)
FX9 =-EXPR*CYP*(SQPR*Z1*CZR+SZR/R3*(X1+1.D0/SQPR))
HY9 = EXPR/P3*SYP*(Z1*CZR+X1/R3*SZR/SQPR)
FZ9 =-EXPR*CYP*(CZR*(1.D0+X1/R3**2/SQPR)-Z1/R3*SZR)
HX9 = FX9*CT1+FZ9*ST1
HZ9 =-FX9*ST1+FZ9*CT1
A1=A(1)+A(2)*CPS
A2=A(3)+A(4)*CPS
A3=A(5)+A(6)*CPS
A4=A(7)+A(8)*CPS
A5=A(9)+A(10)*CPS
A6=A(11)+A(12)*CPS
A7=A(13)+A(14)*CPS
A8=A(15)+A(16)*CPS
A9=A(17)+A(18)*CPS
BX=A1*HX1+A2*HX2+A3*HX3+A4*HX4+A5*HX5+A6*HX6+A7*HX7+A8*HX8+A9*HX9
BY=A1*HY1+A2*HY2+A3*HY3+A4*HY4+A5*HY5+A6*HY6+A7*HY7+A8*HY8+A9*HY9
BZ=A1*HZ1+A2*HZ2+A3*HZ3+A4*HZ4+A5*HZ5+A6*HZ6+A7*HZ7+A8*HZ8+A9*HZ9
c MAKE THE TERMS IN THE 2ND SUM ("PARALLEL" SYMMETRY):
C
C I=1
C
SQQS= DSQRT(1.D0/Q1**2+1.D0/S1**2)
CYQ = DCOS(Y/Q1)
SYQ = DSIN(Y/Q1)
CZS = DCOS(Z2/S1)
SZS = DSIN(Z2/S1)
EXQS= DEXP(SQQS*X2)
FX1 =-SQQS*EXQS*CYQ*CZS *SPS
HY1 = EXQS/Q1*SYQ*CZS *SPS
FZ1 = EXQS*CYQ/S1*SZS *SPS
HX1 = FX1*CT2+FZ1*ST2
HZ1 =-FX1*ST2+FZ1*CT2
SQQS= DSQRT(1.D0/Q1**2+1.D0/S2**2)
CYQ = DCOS(Y/Q1)
SYQ = DSIN(Y/Q1)
CZS = DCOS(Z2/S2)
SZS = DSIN(Z2/S2)
EXQS= DEXP(SQQS*X2)
FX2 =-SQQS*EXQS*CYQ*CZS *SPS
HY2 = EXQS/Q1*SYQ*CZS *SPS
FZ2 = EXQS*CYQ/S2*SZS *SPS
HX2 = FX2*CT2+FZ2*ST2
HZ2 =-FX2*ST2+FZ2*CT2
SQQS= DSQRT(1.D0/Q1**2+1.D0/S3**2)
CYQ = DCOS(Y/Q1)
SYQ = DSIN(Y/Q1)
CZS = DCOS(Z2/S3)
SZS = DSIN(Z2/S3)
EXQS= DEXP(SQQS*X2)
FX3 =-SQQS*EXQS*CYQ*CZS *SPS
HY3 = EXQS/Q1*SYQ*CZS *SPS
FZ3 = EXQS*CYQ/S3*SZS *SPS
HX3 = FX3*CT2+FZ3*ST2
HZ3 =-FX3*ST2+FZ3*CT2
C
C I=2
C
SQQS= DSQRT(1.D0/Q2**2+1.D0/S1**2)
CYQ = DCOS(Y/Q2)
SYQ = DSIN(Y/Q2)
CZS = DCOS(Z2/S1)
SZS = DSIN(Z2/S1)
EXQS= DEXP(SQQS*X2)
FX4 =-SQQS*EXQS*CYQ*CZS *SPS
HY4 = EXQS/Q2*SYQ*CZS *SPS
FZ4 = EXQS*CYQ/S1*SZS *SPS
HX4 = FX4*CT2+FZ4*ST2
HZ4 =-FX4*ST2+FZ4*CT2
SQQS= DSQRT(1.D0/Q2**2+1.D0/S2**2)
CYQ = DCOS(Y/Q2)
SYQ = DSIN(Y/Q2)
CZS = DCOS(Z2/S2)
SZS = DSIN(Z2/S2)
EXQS= DEXP(SQQS*X2)
FX5 =-SQQS*EXQS*CYQ*CZS *SPS
HY5 = EXQS/Q2*SYQ*CZS *SPS
FZ5 = EXQS*CYQ/S2*SZS *SPS
HX5 = FX5*CT2+FZ5*ST2
HZ5 =-FX5*ST2+FZ5*CT2
SQQS= DSQRT(1.D0/Q2**2+1.D0/S3**2)
CYQ = DCOS(Y/Q2)
SYQ = DSIN(Y/Q2)
CZS = DCOS(Z2/S3)
SZS = DSIN(Z2/S3)
EXQS= DEXP(SQQS*X2)
FX6 =-SQQS*EXQS*CYQ*CZS *SPS
HY6 = EXQS/Q2*SYQ*CZS *SPS
FZ6 = EXQS*CYQ/S3*SZS *SPS
HX6 = FX6*CT2+FZ6*ST2
HZ6 =-FX6*ST2+FZ6*CT2
C
C I=3
C
SQQS= DSQRT(1.D0/Q3**2+1.D0/S1**2)
CYQ = DCOS(Y/Q3)
SYQ = DSIN(Y/Q3)
CZS = DCOS(Z2/S1)
SZS = DSIN(Z2/S1)
EXQS= DEXP(SQQS*X2)
FX7 =-SQQS*EXQS*CYQ*CZS *SPS
HY7 = EXQS/Q3*SYQ*CZS *SPS
FZ7 = EXQS*CYQ/S1*SZS *SPS
HX7 = FX7*CT2+FZ7*ST2
HZ7 =-FX7*ST2+FZ7*CT2
SQQS= DSQRT(1.D0/Q3**2+1.D0/S2**2)
CYQ = DCOS(Y/Q3)
SYQ = DSIN(Y/Q3)
CZS = DCOS(Z2/S2)
SZS = DSIN(Z2/S2)
EXQS= DEXP(SQQS*X2)
FX8 =-SQQS*EXQS*CYQ*CZS *SPS
HY8 = EXQS/Q3*SYQ*CZS *SPS
FZ8 = EXQS*CYQ/S2*SZS *SPS
HX8 = FX8*CT2+FZ8*ST2
HZ8 =-FX8*ST2+FZ8*CT2
SQQS= DSQRT(1.D0/Q3**2+1.D0/S3**2)
CYQ = DCOS(Y/Q3)
SYQ = DSIN(Y/Q3)
CZS = DCOS(Z2/S3)
SZS = DSIN(Z2/S3)
EXQS= DEXP(SQQS*X2)
FX9 =-SQQS*EXQS*CYQ*CZS *SPS
HY9 = EXQS/Q3*SYQ*CZS *SPS
FZ9 = EXQS*CYQ/S3*SZS *SPS
HX9 = FX9*CT2+FZ9*ST2
HZ9 =-FX9*ST2+FZ9*CT2
A1=A(19)+A(20)*S2PS
A2=A(21)+A(22)*S2PS
A3=A(23)+A(24)*S2PS
A4=A(25)+A(26)*S2PS
A5=A(27)+A(28)*S2PS
A6=A(29)+A(30)*S2PS
A7=A(31)+A(32)*S2PS
A8=A(33)+A(34)*S2PS
A9=A(35)+A(36)*S2PS
BX=BX+A1*HX1+A2*HX2+A3*HX3+A4*HX4+A5*HX5+A6*HX6+A7*HX7+A8*HX8
* +A9*HX9
BY=BY+A1*HY1+A2*HY2+A3*HY3+A4*HY4+A5*HY5+A6*HY6+A7*HY7+A8*HY8
* +A9*HY9
BZ=BZ+A1*HZ1+A2*HZ2+A3*HZ3+A4*HZ4+A5*HZ5+A6*HZ6+A7*HZ7+A8*HZ8
* +A9*HZ9
C
RETURN
END
c
c############################################################################
c
C
SUBROUTINE DEFORMED (IOPT,PS,X,Y,Z,BX1,BY1,BZ1,BX2,BY2,BZ2)
C
C IOPT - TAIL FIELD MODE FLAG: IOPT=0 - THE TWO TAIL MODES ARE ADDED UP
C IOPT=1 - MODE 1 ONLY
C IOPT=2 - MODE 2 ONLY
C
C CALCULATES GSM COMPONENTS OF TWO UNIT-AMPLITUDE TAIL FIELD MODES,
C TAKING INTO ACCOUNT BOTH EFFECTS OF DIPOLE TILT:
C WARPING IN Y-Z (DONE BY THE S/R WARPED) AND BENDING IN X-Z (DONE BY THIS SUBROUTINE)
C
IMPLICIT REAL*8 (A-H,O-Z)
COMMON /RH0/ RH0
DATA RH2,IEPS /-5.2D0,3/
C
C RH0,RH1,RH2, AND IEPS CONTROL THE TILT-RELATED DEFORMATION OF THE TAIL FIELD
C
SPS=DSIN(PS)
CPS=DSQRT(1.D0-SPS**2)
R2=X**2+Y**2+Z**2
R=SQRT(R2)
ZR=Z/R
RH=RH0+RH2*ZR**2
DRHDR=-ZR/R*2.D0*RH2*ZR
DRHDZ= 2.D0*RH2*ZR/R
C
RRH=R/RH
F=1.D0/(1.D0+RRH**IEPS)**(1.D0/IEPS)
DFDR=-RRH**(IEPS-1)*F**(IEPS+1)/RH
DFDRH=-RRH*DFDR
c
SPSAS=SPS*F
CPSAS=DSQRT(1.D0-SPSAS**2)
C
XAS=X*CPSAS-Z*SPSAS
ZAS=X*SPSAS+Z*CPSAS
C
FACPS=SPS/CPSAS*(DFDR+DFDRH*DRHDR)/R
PSASX=FACPS*X
PSASY=FACPS*Y
PSASZ=FACPS*Z+SPS/CPSAS*DFDRH*DRHDZ
C
DXASDX=CPSAS-ZAS*PSASX
DXASDY=-ZAS*PSASY
DXASDZ=-SPSAS-ZAS*PSASZ
DZASDX=SPSAS+XAS*PSASX
DZASDY=XAS*PSASY
DZASDZ=CPSAS+XAS*PSASZ
FAC1=DXASDZ*DZASDY-DXASDY*DZASDZ
FAC2=DXASDX*DZASDZ-DXASDZ*DZASDX
FAC3=DZASDX*DXASDY-DXASDX*DZASDY
C
C DEFORM:
C
CALL WARPED(IOPT,PS,XAS,Y,ZAS,BXAS1,BYAS1,BZAS1,BXAS2,BYAS2,BZAS2)
C
BX1=BXAS1*DZASDZ-BZAS1*DXASDZ +BYAS1*FAC1
BY1=BYAS1*FAC2
BZ1=BZAS1*DXASDX-BXAS1*DZASDX +BYAS1*FAC3
BX2=BXAS2*DZASDZ-BZAS2*DXASDZ +BYAS2*FAC1
BY2=BYAS2*FAC2
BZ2=BZAS2*DXASDX-BXAS2*DZASDX +BYAS2*FAC3
RETURN
END
C
C------------------------------------------------------------------
c
C
SUBROUTINE WARPED (IOPT,PS,X,Y,Z,BX1,BY1,BZ1,BX2,BY2,BZ2)
C
C CALCULATES GSM COMPONENTS OF THE WARPED FIELD FOR TWO TAIL UNIT MODES.
C THE WARPING DEFORMATION IS IMPOSED ON THE UNWARPED FIELD, COMPUTED
C BY THE S/R "UNWARPED". THE WARPING PARAMETERS WERE TAKEN FROM THE
C RESULTS OF GEOTAIL OBSERVATIONS (TSYGANENKO ET AL. [1998]).
C NB # 6, P.106, OCT 12, 2000.
C
C IOPT - TAIL FIELD MODE FLAG: IOPT=0 - THE TWO TAIL MODES ARE ADDED UP
C IOPT=1 - MODE 1 ONLY
C IOPT=2 - MODE 2 ONLY
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /G/ G
DGDX=0.D0
XL=20.D0
DXLDX=0.D0
SPS=DSIN(PS)
RHO2=Y**2+Z**2
RHO=DSQRT(RHO2)
IF (Y.EQ.0.D0.AND.Z.EQ.0.D0) THEN
PHI=0.D0
CPHI=1.D0
SPHI=0.D0
ELSE
PHI=DATAN2(Z,Y)
CPHI=Y/RHO
SPHI=Z/RHO
ENDIF
RR4L4=RHO/(RHO2**2+XL**4)
F=PHI+G*RHO2*RR4L4*CPHI*SPS
DFDPHI=1.D0-G*RHO2*RR4L4*SPHI*SPS
DFDRHO=G*RR4L4**2*(3.D0*XL**4-RHO2**2)*CPHI*SPS
DFDX=RR4L4*CPHI*SPS*(DGDX*RHO2-G*RHO*RR4L4*4.D0*XL**3*DXLDX)
CF=DCOS(F)
SF=DSIN(F)
YAS=RHO*CF
ZAS=RHO*SF
CALL UNWARPED (IOPT,X,YAS,ZAS,BX_AS1,BY_AS1,BZ_AS1,
* BX_AS2,BY_AS2,BZ_AS2)
BRHO_AS = BY_AS1*CF+BZ_AS1*SF ! DEFORM THE 1ST MODE
BPHI_AS = -BY_AS1*SF+BZ_AS1*CF
BRHO_S = BRHO_AS*DFDPHI
BPHI_S = BPHI_AS-RHO*(BX_AS1*DFDX+BRHO_AS*DFDRHO)
BX1 = BX_AS1*DFDPHI
BY1 = BRHO_S*CPHI-BPHI_S*SPHI
BZ1 = BRHO_S*SPHI+BPHI_S*CPHI ! DONE
BRHO_AS = BY_AS2*CF+BZ_AS2*SF ! DEFORM THE 2ND MODE
BPHI_AS = -BY_AS2*SF+BZ_AS2*CF
BRHO_S = BRHO_AS*DFDPHI
BPHI_S = BPHI_AS-RHO*(BX_AS2*DFDX+BRHO_AS*DFDRHO)
BX2 = BX_AS2*DFDPHI
BY2 = BRHO_S*CPHI-BPHI_S*SPHI
BZ2 = BRHO_S*SPHI+BPHI_S*CPHI ! DONE
RETURN
END
C
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C
SUBROUTINE UNWARPED (IOPT,X,Y,Z,BX1,BY1,BZ1,BX2,BY2,BZ2)
C IOPT - TAIL FIELD MODE FLAG: IOPT=0 - THE TWO TAIL MODES ARE ADDED UP
C IOPT=1 - MODE 1 ONLY
C IOPT=2 - MODE 2 ONLY
C
C CALCULATES GSM COMPONENTS OF THE SHIELDED FIELD OF TWO TAIL MODES WITH UNIT
C AMPLITUDES, WITHOUT ANY WARPING OR BENDING. NONLINEAR PARAMETERS OF THE MODES
C ARE FORWARDED HERE VIA A COMMON BLOCK /TAIL/.
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION A1(60),A2(60) ! TAIL SHIELDING FIELD PARAMETERS FOR THE MODES #1 & #2
COMMON /TAIL/ DXSHIFT1,DXSHIFT2,D0,DELTADY ! ATTENTION: HERE D0 & DELTADY ARE INCLUDED IN /TAIL/
C AND EXCLUDED FROM DATA
DATA DELTADX1,ALPHA1,XSHIFT1
* /1.D0,1.1D0,6.D0/
DATA DELTADX2,ALPHA2,XSHIFT2
* /0.D0,.25D0,4.D0/
DATA A1/-25.45869857,57.35899080,317.5501869,-2.626756717,
*-93.38053698,-199.6467926,-858.8129729,34.09192395,845.4214929,
*-29.07463068,47.10678547,-128.9797943,-781.7512093,6.165038619,
*167.8905046,492.0680410,1654.724031,-46.77337920,-1635.922669,
*40.86186772,-.1349775602,-.9661991179E-01,-.1662302354,
*.002810467517,.2487355077,.1025565237,-14.41750229,-.8185333989,
*11.07693629,.7569503173,-9.655264745,112.2446542,777.5948964,
*-5.745008536,-83.03921993,-490.2278695,-1155.004209,39.08023320,
*1172.780574,-39.44349797,-14.07211198,-40.41201127,-313.2277343,
*2.203920979,8.232835341,197.7065115,391.2733948,-18.57424451,
*-437.2779053,23.04976898,11.75673963,13.60497313,4.691927060,
*18.20923547,27.59044809,6.677425469,1.398283308,2.839005878,
*31.24817706,24.53577264/
DATA A2/-287187.1962,4970.499233,410490.1952,-1347.839052,
*-386370.3240,3317.983750,-143462.3895,5706.513767,171176.2904,
*250.8882750,-506570.8891,5733.592632,397975.5842,9771.762168,
*-941834.2436,7990.975260,54313.10318,447.5388060,528046.3449,
*12751.04453,-21920.98301,-21.05075617,31971.07875,3012.641612,
*-301822.9103,-3601.107387,1797.577552,-6.315855803,142578.8406,
*13161.93640,804184.8410,-14168.99698,-851926.6360,-1890.885671,
*972475.6869,-8571.862853,26432.49197,-2554.752298,-482308.3431,
*-4391.473324,105155.9160,-1134.622050,-74353.53091,-5382.670711,
*695055.0788,-916.3365144,-12111.06667,67.20923358,-367200.9285,
*-21414.14421,14.75567902,20.75638190,59.78601609,16.86431444,
*32.58482365,23.69472951,17.24977936,13.64902647,68.40989058,
*11.67828167/
DATA XM1,XM2/2*-12.D0/
IF (IOPT.EQ.2) GOTO 1
XSC1=(X-XSHIFT1-DXSHIFT1)*ALPHA1-XM1*(ALPHA1-1.D0)
YSC1=Y*ALPHA1
ZSC1=Z*ALPHA1
D0SC1=D0*ALPHA1 ! HERE WE USE A SINGLE VALUE D0 OF THE THICKNESS FOR BOTH MODES
CALL TAILDISK(D0SC1,DELTADX1,DELTADY,XSC1,YSC1,ZSC1,FX1,FY1,FZ1)
CALL SHLCAR5X5(A1,X,Y,Z,DXSHIFT1,HX1,HY1,HZ1)
BX1=FX1+HX1
BY1=FY1+HY1
BZ1=FZ1+HZ1
IF (IOPT.EQ.1) THEN
BX2=0.D0
BY2=0.D0
BZ2=0.D0
RETURN
ENDIF
1 XSC2=(X-XSHIFT2-DXSHIFT2)*ALPHA2-XM2*(ALPHA2-1.D0)
YSC2=Y*ALPHA2
ZSC2=Z*ALPHA2
D0SC2=D0*ALPHA2 ! HERE WE USE A SINGLE VALUE D0 OF THE THICKNESS FOR BOTH MODES
CALL TAILDISK(D0SC2,DELTADX2,DELTADY,XSC2,YSC2,ZSC2,FX2,FY2,FZ2)
CALL SHLCAR5X5(A2,X,Y,Z,DXSHIFT2,HX2,HY2,HZ2)
BX2=FX2+HX2
BY2=FY2+HY2
BZ2=FZ2+HZ2
IF (IOPT.EQ.2) THEN
BX1=0.D0
BY1=0.D0
BZ1=0.D0
RETURN
ENDIF
RETURN
END
C
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
SUBROUTINE TAILDISK(D0,DELTADX,DELTADY,X,Y,Z,BX,BY,BZ)
c
c THIS SUBROUTINE COMPUTES THE COMPONENTS OF THE TAIL CURRENT FIELD,
C SIMILAR TO THAT DESCRIBED BY TSYGANENKO AND PEREDO (1994). THE
C DIFFERENCE IS THAT NOW WE USE SPACEWARPING, AS DESCRIBED IN OUR
C PAPER ON MODELING BIRKELAND CURRENTS (TSYGANENKO AND STERN, 1996)
C INSTEAD OF SHEARING IT IN THE SPIRIT OF T89 TAIL MODEL.
C
IMPLICIT REAL*8 (A-H,O-Z)
c
DIMENSION F(5),B(5),C(5)
C
DATA F /-71.09346626D0,-1014.308601D0,-1272.939359D0,
* -3224.935936D0,-44546.86232D0/
DATA B /10.90101242D0,12.68393898D0,13.51791954D0,14.86775017D0,
* 15.12306404D0/
DATA C /.7954069972D0,.6716601849D0,1.174866319D0,2.565249920D0,
* 10.01986790D0/
C
RHO=DSQRT(X**2+Y**2)
DRHODX=X/RHO
DRHODY=Y/RHO
DEX=DEXP(X/7.D0)
D=D0+DELTADY*(Y/20.D0)**2 +DELTADX*DEX ! THE LAST TERM (INTRODUCED 10/11/2000) MAKES THE SHEET
DDDY=DELTADY*Y*0.005D0 ! THICKEN SUNWARD, TO AVOID PROBLEMS IN THE SUBSOLAR REGION
DDDX=DELTADX/7.D0*DEX
C
DZETA=DSQRT(Z**2+D**2) ! THIS IS THE SAME SIMPLE WAY TO SPREAD
C OUT THE SHEET, AS THAT USED IN T89
DDZETADX=D*DDDX/DZETA
DDZETADY=D*DDDY/DZETA
DDZETADZ=Z/DZETA
C
DBX=0.D0
DBY=0.D0
DBZ=0.D0
C
DO 1 I=1,5
C
BI=B(I)
CI=C(I)
C
S1=DSQRT((RHO+BI)**2+(DZETA+CI)**2)
S2=DSQRT((RHO-BI)**2+(DZETA+CI)**2)
DS1DRHO=(RHO+BI)/S1
DS2DRHO=(RHO-BI)/S2
DS1DDZ=(DZETA+CI)/S1
DS2DDZ=(DZETA+CI)/S2
C
DS1DX=DS1DRHO*DRHODX +DS1DDZ*DDZETADX
DS1DY=DS1DRHO*DRHODY + DS1DDZ*DDZETADY
DS1DZ= DS1DDZ*DDZETADZ
C
DS2DX=DS2DRHO*DRHODX +DS2DDZ*DDZETADX
DS2DY=DS2DRHO*DRHODY + DS2DDZ*DDZETADY
DS2DZ= DS2DDZ*DDZETADZ
C
S1TS2=S1*S2
S1PS2=S1+S2
S1PS2SQ=S1PS2**2
FAC1=DSQRT(S1PS2SQ-(2.D0*BI)**2)
AS=FAC1/(S1TS2*S1PS2SQ)
DASDS1=(1.D0/(FAC1*S2)-AS/S1PS2*(S2*S2+S1*(3.D0*S1+4.D0*S2)))
* /(S1*S1PS2)
DASDS2=(1.D0/(FAC1*S1)-AS/S1PS2*(S1*S1+S2*(3.D0*S2+4.D0*S1)))
* /(S2*S1PS2)
C
DASDX=DASDS1*DS1DX+DASDS2*DS2DX
DASDY=DASDS1*DS1DY+DASDS2*DS2DY
DASDZ=DASDS1*DS1DZ+DASDS2*DS2DZ
C
DBX=DBX-F(I)*X*DASDZ
DBY=DBY-F(I)*Y*DASDZ
DBZ=DBZ+F(I)*(2.D0*AS+X*DASDX+Y*DASDY)
1 CONTINUE
BX=DBX
BY=DBY
BZ=DBZ
RETURN
END
C
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C THIS CODE RETURNS THE SHIELDING FIELD REPRESENTED BY 5x5=25 "CARTESIAN"
C HARMONICS
C
SUBROUTINE SHLCAR5X5(A,X,Y,Z,DSHIFT,HX,HY,HZ)
C
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C The NLIN coefficients are the amplitudes of the "cartesian"
c harmonics (A(1)-A(NLIN).
c The NNP nonlinear parameters (A(NLIN+1)-A(NTOT) are the scales Pi and Ri
C entering the arguments of exponents, sines, and cosines in each of the
C NLIN "Cartesian" harmonics
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
IMPLICIT REAL * 8 (A - H, O - Z)
C
DIMENSION A(60)
C
DHX=0.D0
DHY=0.D0
DHZ=0.D0
L=0
C
DO I=1,5
RP=1.D0/A(50+I)
CYPI=DCOS(Y*RP)
SYPI=DSIN(Y*RP)
C
DO K=1,5
RR=1.D0/A(55+K)
SZRK=DSIN(Z*RR)
CZRK=DCOS(Z*RR)
SQPR=DSQRT(RP**2+RR**2)
EPR=DEXP(X*SQPR)
C
DBX=-SQPR*EPR*CYPI*SZRK
DBY= RP*EPR*SYPI*SZRK
DBZ=-RR*EPR*CYPI*CZRK
L=L+2
COEF=A(L-1)+A(L)*DSHIFT
DHX=DHX+COEF*DBX
DHY=DHY+COEF*DBY
DHZ=DHZ+COEF*DBZ
c
END DO
END DO
HX=DHX
HY=DHY
HZ=DHZ
C
RETURN
END
c
c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C
SUBROUTINE BIRK_TOT (IOPB,PS,X,Y,Z,BX11,BY11,BZ11,BX12,BY12,BZ12,
* BX21,BY21,BZ21,BX22,BY22,BZ22)
C
C IOPB - BIRKELAND FIELD MODE FLAG:
C IOPB=0 - ALL COMPONENTS
C IOPB=1 - REGION 1, MODES 1 & 2
C IOPB=2 - REGION 2, MODES 1 & 2
C
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION SH11(86),SH12(86),SH21(86),SH22(86)
COMMON /BIRKPAR/ XKAPPA1,XKAPPA2 ! INPUT PARAMETERS, SPECIFIED FROM A MAIN PROGRAM
COMMON /DPHI_B_RHO0/ DPHI,B,RHO_0,XKAPPA ! PARAMETERS, CONTROL DAY-NIGHT ASYMMETRY OF F.A.C.
DATA SH11/46488.84663,-15541.95244,-23210.09824,-32625.03856,
*-109894.4551,-71415.32808,58168.94612,55564.87578,-22890.60626,
*-6056.763968,5091.368100,239.7001538,-13899.49253,4648.016991,
*6971.310672,9699.351891,32633.34599,21028.48811,-17395.96190,
*-16461.11037,7447.621471,2528.844345,-1934.094784,-588.3108359,
*-32588.88216,10894.11453,16238.25044,22925.60557,77251.11274,
*50375.97787,-40763.78048,-39088.60660,15546.53559,3559.617561,
*-3187.730438,309.1487975,88.22153914,-243.0721938,-63.63543051,
*191.1109142,69.94451996,-187.9539415,-49.89923833,104.0902848,
*-120.2459738,253.5572433,89.25456949,-205.6516252,-44.93654156,
*124.7026309,32.53005523,-98.85321751,-36.51904756,98.88241690,
*24.88493459,-55.04058524,61.14493565,-128.4224895,-45.35023460,
*105.0548704,-43.66748755,119.3284161,31.38442798,-92.87946767,
*-33.52716686,89.98992001,25.87341323,-48.86305045,59.69362881,
*-126.5353789,-44.39474251,101.5196856,59.41537992,41.18892281,
*80.86101200,3.066809418,7.893523804,30.56212082,10.36861082,
*8.222335945,19.97575641,2.050148531,4.992657093,2.300564232,
*.2256245602,-.05841594319/
DATA SH12/210260.4816,-1443587.401,-1468919.281,281939.2993,
*-1131124.839,729331.7943,2573541.307,304616.7457,468887.5847,
*181554.7517,-1300722.650,-257012.8601,645888.8041,-2048126.412,
*-2529093.041,571093.7972,-2115508.353,1122035.951,4489168.802,
*75234.22743,823905.6909,147926.6121,-2276322.876,-155528.5992,
*-858076.2979,3474422.388,3986279.931,-834613.9747,3250625.781,
*-1818680.377,-7040468.986,-414359.6073,-1295117.666,-346320.6487,
*3565527.409,430091.9496,-.1565573462,7.377619826,.4115646037,
*-6.146078880,3.808028815,-.5232034932,1.454841807,-12.32274869,
*-4.466974237,-2.941184626,-.6172620658,12.64613490,1.494922012,
*-21.35489898,-1.652256960,16.81799898,-1.404079922,-24.09369677,
*-10.99900839,45.94237820,2.248579894,31.91234041,7.575026816,
*-45.80833339,-1.507664976,14.60016998,1.348516288,-11.05980247,
*-5.402866968,31.69094514,12.28261196,-37.55354174,4.155626879,
*-33.70159657,-8.437907434,36.22672602,145.0262164,70.73187036,
*85.51110098,21.47490989,24.34554406,31.34405345,4.655207476,
*5.747889264,7.802304187,1.844169801,4.867254550,2.941393119,
*.1379899178,.06607020029/
DATA SH21/162294.6224,503885.1125,-27057.67122,-531450.1339,
*84747.05678,-237142.1712,84133.61490,259530.0402,69196.05160,
*-189093.5264,-19278.55134,195724.5034,-263082.6367,-818899.6923,
*43061.10073,863506.6932,-139707.9428,389984.8850,-135167.5555,
*-426286.9206,-109504.0387,295258.3531,30415.07087,-305502.9405,
*100785.3400,315010.9567,-15999.50673,-332052.2548,54964.34639,
*-152808.3750,51024.67566,166720.0603,40389.67945,-106257.7272,
*-11126.14442,109876.2047,2.978695024,558.6019011,2.685592939,
*-338.0004730,-81.99724090,-444.1102659,89.44617716,212.0849592,
*-32.58562625,-982.7336105,-35.10860935,567.8931751,-1.917212423,
*-260.2023543,-1.023821735,157.5533477,23.00200055,232.0603673,
*-36.79100036,-111.9110936,18.05429984,447.0481000,15.10187415,
*-258.7297813,-1.032340149,-298.6402478,-1.676201415,180.5856487,
*64.52313024,209.0160857,-53.85574010,-98.52164290,14.35891214,
*536.7666279,20.09318806,-309.7349530,58.54144539,67.45226850,
*97.92374406,4.752449760,10.46824379,32.91856110,12.05124381,
*9.962933904,15.91258637,1.804233877,6.578149088,2.515223491,
*.1930034238,-.02261109942/
DATA SH22/-131287.8986,-631927.6885,-318797.4173,616785.8782,
*-50027.36189,863099.9833,47680.20240,-1053367.944,-501120.3811,
*-174400.9476,222328.6873,333551.7374,-389338.7841,-1995527.467,
*-982971.3024,1960434.268,297239.7137,2676525.168,-147113.4775,
*-3358059.979,-2106979.191,-462827.1322,1017607.960,1039018.475,
*520266.9296,2627427.473,1301981.763,-2577171.706,-238071.9956,
*-3539781.111,94628.16420,4411304.724,2598205.733,637504.9351,
*-1234794.298,-1372562.403,-2.646186796,-31.10055575,2.295799273,
*19.20203279,30.01931202,-302.1028550,-14.78310655,162.1561899,
*.4943938056,176.8089129,-.2444921680,-100.6148929,9.172262228,
*137.4303440,-8.451613443,-84.20684224,-167.3354083,1321.830393,
*76.89928813,-705.7586223,18.28186732,-770.1665162,-9.084224422,
*436.3368157,-6.374255638,-107.2730177,6.080451222,65.53843753,
*143.2872994,-1028.009017,-64.22739330,547.8536586,-20.58928632,
*597.3893669,10.17964133,-337.7800252,159.3532209,76.34445954,
*84.74398828,12.76722651,27.63870691,32.69873634,5.145153451,
*6.310949163,6.996159733,1.971629939,4.436299219,2.904964304,
*.1486276863,.06859991529/
XKAPPA=XKAPPA1 ! FORWARDED IN BIRK_1N2
X_SC=XKAPPA1-1.1D0 ! FORWARDED IN BIRK_SHL
IF (IOPB.EQ.0.OR.IOPB.EQ.1) THEN
CALL BIRK_1N2 (1,1,PS,X,Y,Z,FX11,FY11,FZ11) ! REGION 1, MODE 1
CALL BIRK_SHL (SH11,PS,X_SC,X,Y,Z,HX11,HY11,HZ11)
BX11=FX11+HX11
BY11=FY11+HY11
BZ11=FZ11+HZ11
CALL BIRK_1N2 (1,2,PS,X,Y,Z,FX12,FY12,FZ12) ! REGION 1, MODE 2
CALL BIRK_SHL (SH12,PS,X_SC,X,Y,Z,HX12,HY12,HZ12)
BX12=FX12+HX12
BY12=FY12+HY12
BZ12=FZ12+HZ12
ENDIF
XKAPPA=XKAPPA2 ! FORWARDED IN BIRK_1N2
X_SC=XKAPPA2-1.0D0 ! FORWARDED IN BIRK_SHL
IF (IOPB.EQ.0.OR.IOPB.EQ.2) THEN
CALL BIRK_1N2 (2,1,PS,X,Y,Z,FX21,FY21,FZ21) ! REGION 2, MODE 1
CALL BIRK_SHL (SH21,PS,X_SC,X,Y,Z,HX21,HY21,HZ21)
BX21=FX21+HX21
BY21=FY21+HY21
BZ21=FZ21+HZ21
CALL BIRK_1N2 (2,2,PS,X,Y,Z,FX22,FY22,FZ22) ! REGION 2, MODE 2
CALL BIRK_SHL (SH22,PS,X_SC,X,Y,Z,HX22,HY22,HZ22)
BX22=FX22+HX22
BY22=FY22+HY22
BZ22=FZ22+HZ22
ENDIF
RETURN
END
C
c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c
SUBROUTINE BIRK_1N2 (NUMB,MODE,PS,X,Y,Z,BX,BY,BZ)
C
C CALCULATES COMPONENTS OF REGION 1/2 FIELD IN SPHERICAL COORDS. DERIVED FROM THE S/R DIPDEF2C (WHICH
C DOES THE SAME JOB, BUT INPUT/OUTPUT THERE WAS IN SPHERICAL COORDS, WHILE HERE WE USE CARTESIAN ONES)
C
C INPUT: NUMB=1 (2) FOR REGION 1 (2) CURRENTS
C MODE=1 YIELDS SIMPLE SINUSOIDAL MLT VARIATION, WITH MAXIMUM CURRENT AT DAWN/DUSK MERIDIAN
C WHILE MODE=2 YIELDS THE SECOND HARMONIC.
C
C
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION A11(31),A12(31),A21(31),A22(31)
COMMON /MODENUM/ M
COMMON /DTHETA/ DTHETA
COMMON /DPHI_B_RHO0/ DPHI,B,RHO_0,XKAPPA ! THESE PARAMETERS CONTROL DAY-NIGHT ASYMMETRY OF F.A.C., AS FOLLOWS:
C (1) DPHI: HALF-DIFFERENCE (IN RADIANS) BETWEEN DAY AND NIGHT LATITUDE OF FAC OVAL AT IONOSPHERIC ALTITUDE;
C TYPICAL VALUE: 0.06
C (2) B: AN ASYMMETRY FACTOR AT HIGH-ALTITUDES; FOR B=0, THE ONLY ASYMMETRY IS THAT FROM DPHI
C TYPICAL VALUES: 0.35-0.70
C (3) RHO_0: A FIXED PARAMETER, DEFINING THE DISTANCE RHO, AT WHICH THE LATITUDE SHIFT GRADUALLY SATURATES AND
C STOPS INCREASING
C ITS VALUE WAS ASSUMED FIXED, EQUAL TO 7.0.
C (4) XKAPPA: AN OVERALL SCALING FACTOR, WHICH CAN BE USED FOR CHANGING THE SIZE OF THE F.A.C. OVAL
C
DATA BETA,RH,EPS/0.9D0,10.D0,3.D0/ ! parameters of the tilt-dependent deformation of the untilted F.A.C. field
DATA A11/.1618068350,-.1797957553,2.999642482,-.9322708978,
*-.6811059760,.2099057262,-8.358815746,-14.86033550,.3838362986,
*-16.30945494,4.537022847,2.685836007,27.97833029,6.330871059,
*1.876532361,18.95619213,.9651528100,.4217195118,-.08957770020,
*-1.823555887,.7457045438,-.5785916524,-1.010200918,.01112389357,
*.09572927448,-.3599292276,8.713700514,.9763932955,3.834602998,
*2.492118385,.7113544659/
DATA A12/.7058026940,-.2845938535,5.715471266,-2.472820880,
*-.7738802408,.3478293930,-11.37653694,-38.64768867,.6932927651,
*-212.4017288,4.944204937,3.071270411,33.05882281,7.387533799,
*2.366769108,79.22572682,.6154290178,.5592050551,-.1796585105,
*-1.654932210,.7309108776,-.4926292779,-1.130266095,-.009613974555,
*.1484586169,-.2215347198,7.883592948,.02768251655,2.950280953,
*1.212634762,.5567714182/
DATA A21/.1278764024,-.2320034273,1.805623266,-32.37241440,
*-.9931490648,.3175085630,-2.492465814,-16.21600096,.2695393416,
*-6.752691265,3.971794901,14.54477563,41.10158386,7.912889730,
*1.258297372,9.583547721,1.014141963,.5104134759,-.1790430468,
*-1.756358428,.7561986717,-.6775248254,-.04014016420,.01446794851,
*.1200521731,-.2203584559,4.508963850,.8221623576,1.779933730,
*1.102649543,.8867880020/
DATA A22/.4036015198,-.3302974212,2.827730930,-45.44405830,
*-1.611103927,.4927112073,-.003258457559,-49.59014949,.3796217108,
*-233.7884098,4.312666980,18.05051709,28.95320323,11.09948019,
*.7471649558,67.10246193,.5667096597,.6468519751,-.1560665317,
*-1.460805289,.7719653528,-.6658988668,.2515179349E-05,
*.02426021891,.1195003324,-.2625739255,4.377172556,.2421190547,
*2.503482679,1.071587299,.7247997430/
B=0.5
RHO_0=7.0
M=MODE
IF (NUMB.EQ.1) THEN
DPHI=0.055D0
DTHETA=0.06D0
ENDIF
IF (NUMB.EQ.2) THEN
DPHI=0.030D0
DTHETA=0.09D0
ENDIF
Xsc=X*XKAPPA
Ysc=Y*XKAPPA
Zsc=Z*XKAPPA
RHO=DSQRT(Xsc**2+Zsc**2)
Rsc=DSQRT(Xsc**2+Ysc**2+Zsc**2) ! SCALED
RHO2=RHO_0**2
IF (Xsc.EQ.0.D0.AND.Zsc.EQ.0.D0) THEN
PHI=0.D0
ELSE
PHI=DATAN2(-Zsc,Xsc) ! FROM CARTESIAN TO CYLINDRICAL (RHO,PHI,Y)
ENDIF
SPHIC=DSIN(PHI)
CPHIC=DCOS(PHI) ! "C" means "CYLINDRICAL", TO DISTINGUISH FROM SPHERICAL PHI
BRACK=DPHI+B*RHO2/(RHO2+1.D0)*(RHO**2-1.D0)/(RHO2+RHO**2)
R1RH=(Rsc-1.D0)/RH
PSIAS=BETA*PS/(1.D0+R1RH**EPS)**(1.D0/EPS)
PHIS=PHI-BRACK*DSIN(PHI) -PSIAS
DPHISPHI=1.D0-BRACK*DCOS(PHI)
DPHISRHO=-2.D0*B*RHO2*RHO/(RHO2+RHO**2)**2 *DSIN(PHI)
* +BETA*PS*R1RH**(EPS-1.D0)*RHO/(RH*Rsc*
* (1.D0+R1RH**EPS)**(1.D0/EPS+1.D0))
DPHISDY= BETA*PS*R1RH**(EPS-1.D0)*Ysc/(RH*Rsc*
* (1.D0+R1RH**EPS)**(1.D0/EPS+1.D0))
SPHICS=DSIN(PHIS)
CPHICS=DCOS(PHIS)
XS= RHO*CPHICS
ZS=-RHO*SPHICS
IF (NUMB.EQ.1) THEN
IF (MODE.EQ.1) CALL TWOCONES (A11,XS,Ysc,ZS,BXS,BYAS,BZS)
IF (MODE.EQ.2) CALL TWOCONES (A12,XS,Ysc,ZS,BXS,BYAS,BZS)
ELSE
IF (MODE.EQ.1) CALL TWOCONES (A21,XS,Ysc,ZS,BXS,BYAS,BZS)
IF (MODE.EQ.2) CALL TWOCONES (A22,XS,Ysc,ZS,BXS,BYAS,BZS)
ENDIF
BRHOAS=BXS*CPHICS-BZS*SPHICS
BPHIAS=-BXS*SPHICS-BZS*CPHICS
BRHO_S=BRHOAS*DPHISPHI *XKAPPA ! SCALING
BPHI_S=(BPHIAS-RHO*(BYAS*DPHISDY+BRHOAS*DPHISRHO)) *XKAPPA
BY_S=BYAS*DPHISPHI *XKAPPA
BX=BRHO_S*CPHIC-BPHI_S*SPHIC
BY=BY_S
BZ=-BRHO_S*SPHIC-BPHI_S*CPHIC
RETURN
END
c
C=========================================================================
c
SUBROUTINE TWOCONES (A,X,Y,Z,BX,BY,BZ)
C
C ADDS FIELDS FROM TWO CONES (NORTHERN AND SOUTHERN), WITH A PROPER SYMMETRY OF THE CURRENT AND FIELD,
C CORRESPONDING TO THE REGION 1 BIRKELAND CURRENTS.
C
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION A(31)
CALL ONE_CONE (A,X,Y,Z,BXN,BYN,BZN)
CALL ONE_CONE (A,X,-Y,-Z,BXS,BYS,BZS)
BX=BXN-BXS
BY=BYN+BYS
BZ=BZN+BZS
RETURN
END
c
C-------------------------------------------------------------------------
C
SUBROUTINE ONE_CONE(A,X,Y,Z,BX,BY,BZ)
c
c RETURNS FIELD COMPONENTS FOR A DEFORMED CONICAL CURRENT SYSTEM, FITTED TO A BIOSAVART FIELD
c BY SIM_14.FOR. HERE ONLY THE NORTHERN CONE IS TAKEN INTO ACCOUNT.
c
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION A(31)
COMMON /DTHETA/ DTHETA
COMMON /MODENUM/ M
DATA DR,DT/1.D-6,1.D-6/ ! JUST FOR NUMERICAL DIFFERENTIATION
THETA0=A(31)
RHO2=X**2+Y**2
RHO=DSQRT(RHO2)
R=DSQRT(RHO2+Z**2)
THETA=DATAN2(RHO,Z)
PHI=DATAN2(Y,X)
C
C MAKE THE DEFORMATION OF COORDINATES:
C
RS=R_S(A,R,THETA)
THETAS=THETA_S(A,R,THETA)
PHIS=PHI
C
C CALCULATE FIELD COMPONENTS AT THE NEW POSITION (ASTERISKED):
C
CALL FIALCOS (RS,THETAS,PHIS,BTAST,BFAST,M,THETA0,DTHETA) ! MODE #M
C
C NOW TRANSFORM B{R,T,F}_AST BY THE DEFORMATION TENSOR:
C
C FIRST OF ALL, FIND THE DERIVATIVES:
C
DRSDR=(R_S(A,R+DR,THETA)-R_S(A,R-DR,THETA))/(2.D0*DR)
DRSDT=(R_S(A,R,THETA+DT)-R_S(A,R,THETA-DT))/(2.D0*DT)
DTSDR=(THETA_S(A,R+DR,THETA)-THETA_S(A,R-DR,THETA))/(2.D0*DR)
DTSDT=(THETA_S(A,R,THETA+DT)-THETA_S(A,R,THETA-DT))/(2.D0*DT)
STSST=DSIN(THETAS)/DSIN(THETA)
RSR=RS/R
BR =-RSR/R*STSST*BTAST*DRSDT
BTHETA = RSR*STSST*BTAST*DRSDR
BPHI = RSR*BFAST*(DRSDR*DTSDT-DRSDT*DTSDR)
S=RHO/R
C=Z/R
SF=Y/RHO
CF=X/RHO
BE=BR*S+BTHETA*C
BX=A(1)*(BE*CF-BPHI*SF)
BY=A(1)*(BE*SF+BPHI*CF)
BZ=A(1)*(BR*C-BTHETA*S)
RETURN
END
C
C=====================================================================================
DOUBLE PRECISION FUNCTION R_S(A,R,THETA)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION A(31)
C
R_S=R+A(2)/R+A(3)*R/DSQRT(R**2+A(11)**2)+A(4)*R/(R**2+A(12)**2)
*+(A(5)+A(6)/R+A(7)*R/DSQRT(R**2+A(13)**2)+A(8)*R/(R**2+A(14)**2))*
* DCOS(THETA)
*+(A(9)*R/DSQRT(R**2+A(15)**2)+A(10)*R/(R**2+A(16)**2)**2)
* *DCOS(2.D0*THETA)
C
RETURN
END
C
C-----------------------------------------------------------------------------
C
DOUBLE PRECISION FUNCTION THETA_S(A,R,THETA)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION A(31)
c
THETA_S=THETA+(A(17)+A(18)/R+A(19)/R**2
* +A(20)*R/DSQRT(R**2+A(27)**2))*DSIN(THETA)
* +(A(21)+A(22)*R/DSQRT(R**2+A(28)**2)
* +A(23)*R/(R**2+A(29)**2))*DSIN(2.D0*THETA)
* +(A(24)+A(25)/R+A(26)*R/(R**2+A(30)**2))*DSIN(3.D0*THETA)
C
RETURN
END
C
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
SUBROUTINE FIALCOS(R,THETA,PHI,BTHETA,BPHI,N,THETA0,DT)
C
C CONICAL MODEL OF BIRKELAND CURRENT FIELD; BASED ON THE OLD S/R FIALCO (OF 1990-91)
C BTN, AND BPN ARE THE ARRAYS OF BTHETA AND BPHI (BTN(i), BPN(i) CORRESPOND TO i-th MODE).
C ONLY FIRST N MODE AMPLITUDES ARE COMPUTED (N<=10).
C THETA0 IS THE ANGULAR HALF-WIDTH OF THE CONE, DT IS THE ANGULAR H.-W. OF THE CURRENT LAYER
C NOTE: BR=0 (BECAUSE ONLY RADIAL CURRENTS ARE PRESENT IN THIS MODEL)
C
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION BTN(10),BPN(10),CCOS(10),SSIN(10)
SINTE=DSIN(THETA)
RO=R*SINTE
COSTE=DCOS(THETA)
SINFI=DSIN(PHI)
COSFI=DCOS(PHI)
TG=SINTE/(1.D0+COSTE) ! TAN(THETA/2)
CTG=SINTE/(1.D0-COSTE) ! CTG(THETA/2)
C
C
TETANP=THETA0+DT
TETANM=THETA0-DT
IF(THETA.LT.TETANM) GOTO 1
TGP=DTAN(TETANP*0.5D0)
TGM=DTAN(TETANM*0.5D0)
TGM2=TGM*TGM
TGP2=TGP*TGP
1 CONTINUE
COSM1=1.D0
SINM1=0.D0
TM=1.D0
TGM2M=1.D0
TGP2M=1.D0
DO 2 M=1,N
TM=TM*TG
CCOS(M)=COSM1*COSFI-SINM1*SINFI
SSIN(M)=SINM1*COSFI+COSM1*SINFI
COSM1=CCOS(M)
SINM1=SSIN(M)
IF(THETA.LT.TETANM) THEN
T=TM
DTT=0.5D0*M*TM*(TG+CTG)
DTT0=0.D0
ELSE IF(THETA.LT.TETANP) THEN
TGM2M=TGM2M*TGM2
FC=1.D0/(TGP-TGM)
FC1=1.D0/(2*M+1)
TGM2M1=TGM2M*TGM
TG21=1.D0+TG*TG
T=FC*(TM*(TGP-TG)+FC1*(TM*TG-TGM2M1/TM))
DTT=0.5D0*M*FC*TG21*(TM/TG*(TGP-TG)-FC1*(TM-TGM2M1/(TM*TG)))
DTT0=0.5D0*FC*((TGP+TGM)*(TM*TG-FC1*(TM*TG-TGM2M1/TM))+
* TM*(1.D0-TGP*TGM)-(1.D0+TGM2)*TGM2M/TM)
ELSE
TGP2M=TGP2M*TGP2
TGM2M=TGM2M*TGM2
FC=1.D0/(TGP-TGM)
FC1=1.D0/(2*M+1)
T=FC*FC1*(TGP2M*TGP-TGM2M*TGM)/TM
DTT=-T*M*0.5D0*(TG+CTG)
ENDIF
BTN(M)=M*T*CCOS(M)/RO
BPN(M)=-DTT*SSIN(M)/R
2 CONTINUE
BTHETA=BTN(N) *800.
BPHI =BPN(N) *800.
RETURN
END
C
C-------------------------------------------------------------------------
C
C
SUBROUTINE BIRK_SHL (A,PS,X_SC,X,Y,Z,BX,BY,BZ)
C
IMPLICIT REAL * 8 (A - H, O - Z)
DIMENSION A(86)
C
CPS=DCOS(PS)
SPS=DSIN(PS)
S3PS=2.D0*CPS
C
PST1=PS*A(85)
PST2=PS*A(86)
ST1=DSIN(PST1)
CT1=DCOS(PST1)
ST2=DSIN(PST2)
CT2=DCOS(PST2)
X1=X*CT1-Z*ST1
Z1=X*ST1+Z*CT1
X2=X*CT2-Z*ST2
Z2=X*ST2+Z*CT2
C
L=0
GX=0.D0
GY=0.D0
GZ=0.D0
C
DO 1 M=1,2 ! M=1 IS FOR THE 1ST SUM ("PERP." SYMMETRY)
C AND M=2 IS FOR THE SECOND SUM ("PARALL." SYMMETRY)
DO 2 I=1,3
P=A(72+I)
Q=A(78+I)
CYPI=DCOS(Y/P)
CYQI=DCOS(Y/Q)
SYPI=DSIN(Y/P)
SYQI=DSIN(Y/Q)
C
DO 3 K=1,3
R=A(75+K)
S=A(81+K)
SZRK=DSIN(Z1/R)
CZSK=DCOS(Z2/S)
CZRK=DCOS(Z1/R)
SZSK=DSIN(Z2/S)
SQPR=DSQRT(1.D0/P**2+1.D0/R**2)
SQQS=DSQRT(1.D0/Q**2+1.D0/S**2)
EPR=DEXP(X1*SQPR)
EQS=DEXP(X2*SQQS)
C
DO 4 N=1,2 ! N=1 IS FOR THE FIRST PART OF EACH COEFFICIENT
C AND N=2 IS FOR THE SECOND ONE
DO 5 NN=1,2 ! NN = 1,2 FURTHER SPLITS THE COEFFICIENTS INTO 2 PARTS,
C TO TAKE INTO ACCOUNT THE SCALE FACTOR DEPENDENCE
IF (M.EQ.1) THEN
FX=-SQPR*EPR*CYPI*SZRK
FY=EPR*SYPI*SZRK/P
FZ=-EPR*CYPI*CZRK/R
IF (N.EQ.1) THEN
IF (NN.EQ.1) THEN
HX=FX
HY=FY
HZ=FZ
ELSE
HX=FX*X_SC
HY=FY*X_SC
HZ=FZ*X_SC
ENDIF
ELSE
IF (NN.EQ.1) THEN
HX=FX*CPS
HY=FY*CPS
HZ=FZ*CPS
ELSE
HX=FX*CPS*X_SC
HY=FY*CPS*X_SC
HZ=FZ*CPS*X_SC
ENDIF
ENDIF
ELSE ! M.EQ.2
FX=-SPS*SQQS*EQS*CYQI*CZSK
FY=SPS/Q*EQS*SYQI*CZSK
FZ=SPS/S*EQS*CYQI*SZSK
IF (N.EQ.1) THEN
IF (NN.EQ.1) THEN
HX=FX
HY=FY
HZ=FZ
ELSE
HX=FX*X_SC
HY=FY*X_SC
HZ=FZ*X_SC
ENDIF
ELSE
IF (NN.EQ.1) THEN
HX=FX*S3PS
HY=FY*S3PS
HZ=FZ*S3PS
ELSE
HX=FX*S3PS*X_SC
HY=FY*S3PS*X_SC
HZ=FZ*S3PS*X_SC
ENDIF
ENDIF
ENDIF
L=L+1
IF (M.EQ.1) THEN
HXR=HX*CT1+HZ*ST1
HZR=-HX*ST1+HZ*CT1
ELSE
HXR=HX*CT2+HZ*ST2
HZR=-HX*ST2+HZ*CT2
ENDIF
GX=GX+HXR*A(L)
GY=GY+HY *A(L)
GZ=GZ+HZR*A(L)
5 CONTINUE
4 CONTINUE
3 CONTINUE
2 CONTINUE
1 CONTINUE
BX=GX
BY=GY
BZ=GZ
RETURN
END
C
C************************************************************************************
C
SUBROUTINE FULL_RC (IOPR,PS,X,Y,Z,BXSRC,BYSRC,BZSRC,BXPRC,BYPRC,
* BZPRC)
C
C CALCULATES GSM FIELD COMPONENTS OF THE SYMMETRIC (SRC) AND PARTIAL (PRC) COMPONENTS OF THE RING CURRENT
C SRC PROVIDES A DEPRESSION OF -28 nT AT EARTH
C PRC CORRESPONDS TO THE PRESSURE DIFFERENCE OF 2 nPa BETWEEN MIDNIGHT AND NOON RING CURRENT
C PARTICLE PRESSURE AND YIELDS A DEPRESSION OF -17 nT AT X=-6Re
C
C SC_SY AND SC_PR ARE SCALING FACTORS FOR THE SYMMETRIC AND PARTIAL COMPONENTS:
C VALUES LARGER THAN 1 RESULT IN SPATIALLY LARGER CURRENTS
C
C PHI IS THE ROTATION ANGLE IN RADIANS OF THE PARTIAL RING CURRENT (MEASURED FROM MIDNIGHT TOWARD DUSK)
C
C IOPR - A RING CURRENT CALCULATION FLAG (FOR LEAST-SQUARES FITTING ONLY):
C IOPR=0 - BOTH SRC AND PRC FIELDS ARE CALCULATED
C IOPR=1 - SRC ONLY
C IOPR=2 - PRC ONLY
C
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION C_SY(86),C_PR(86)
COMMON /RCPAR/ SC_SY,SC_PR,PHI
C
DATA C_SY/-957.2534900,-817.5450246,583.2991249,758.8568270, ! CORRECTED VALUES (AS OF MAY 2006)
*13.17029064,68.94173502,-15.29764089,-53.43151590,27.34311724,
*149.5252826,-11.00696044,-179.7031814,953.0914774,817.2340042,
*-581.0791366,-757.5387665,-13.10602697,-68.58155678,15.22447386,
*53.15535633,-27.07982637,-149.1413391,10.91433279,179.3251739,
*-6.028703251,1.303196101,-1.345909343,-1.138296330,-0.06642634348,
*-0.3795246458,.07487833559,.2891156371,-.5506314391,-.4443105812,
*0.2273682152,0.01086886655,-9.130025352,1.118684840,1.110838825,
*.1219761512,-.06263009645,-.1896093743,.03434321042,.01523060688,
*-.4913171541,-.2264814165,-.04791374574,.1981955976,-68.32678140,
*-48.72036263,14.03247808,16.56233733,2.369921099,6.200577111,
*-1.415841250,-0.8184867835,-3.401307527,-8.490692287,3.217860767,
*-9.037752107,66.09298105,48.23198578,-13.67277141,-16.27028909,
*-2.309299411,-6.016572391,1.381468849,0.7935312553,3.436934845,
* 8.260038635,-3.136213782,8.833214943,8.041075485,8.024818618,
* 35.54861873,12.55415215,1.738167799,3.721685353,23.06768025,
* 6.871230562,6.806229878,21.35990364,1.687412298,3.500885177,
* 0.3498952546,0.6595919814/
DATA C_PR/-64820.58481,-63965.62048,66267.93413,135049.7504,
*-36.56316878,124.6614669,56.75637955,-87.56841077,5848.631425,
*4981.097722,-6233.712207,-10986.40188,68716.52057,65682.69473,
*-69673.32198,-138829.3568,43.45817708,-117.9565488,-62.14836263,
*79.83651604,-6211.451069,-5151.633113,6544.481271,11353.03491,
*23.72352603,-256.4846331,25.77629189,145.2377187,-4.472639098,
*-3.554312754,2.936973114,2.682302576,2.728979958,26.43396781,
*-9.312348296,-29.65427726,-247.5855336,-206.9111326,74.25277664,
*106.4069993,15.45391072,16.35943569,-5.965177750,-6.079451700,
*115.6748385,-35.27377307,-32.28763497,-32.53122151,93.74409310,
*84.25677504,-29.23010465,-43.79485175,-6.434679514,-6.620247951,
*2.443524317,2.266538956,-43.82903825,6.904117876,12.24289401,
*17.62014361,152.3078796,124.5505289,-44.58690290,-63.02382410,
*-8.999368955,-9.693774119,3.510930306,3.770949738,-77.96705716,
*22.07730961,20.46491655,18.67728847,9.451290614,9.313661792,
*644.7620970,418.2515954,7.183754387,35.62128817,19.43180682,
*39.57218411,15.69384715,7.123215241,2.300635346,21.90881131,
*-.01775839370,.3996346710/
CALL SRC_PRC (IOPR,SC_SY,SC_PR,PHI,PS,X,Y,Z,HXSRC,HYSRC,HZSRC,
* HXPRC,HYPRC,HZPRC)
X_SC=SC_SY-1.D0
IF (IOPR.EQ.0.OR.IOPR.EQ.1) THEN
CALL RC_SHIELD (C_SY,PS,X_SC,X,Y,Z,FSX,FSY,FSZ)
ELSE
FSX=0.D0
FSY=0.D0
FSZ=0.D0
ENDIF
X_SC=SC_PR-1.D0
IF (IOPR.EQ.0.OR.IOPR.EQ.2) THEN
CALL RC_SHIELD (C_PR,PS,X_SC,X,Y,Z,FPX,FPY,FPZ)
ELSE
FPX=0.D0
FPY=0.D0
FPZ=0.D0
ENDIF
BXSRC=HXSRC+FSX
BYSRC=HYSRC+FSY
BZSRC=HZSRC+FSZ
BXPRC=HXPRC+FPX
BYPRC=HYPRC+FPY
BZPRC=HZPRC+FPZ
RETURN
END
C---------------------------------------------------------------------------------------
C
SUBROUTINE SRC_PRC (IOPR,SC_SY,SC_PR,PHI,PS,X,Y,Z,BXSRC,BYSRC,
* BZSRC,BXPRC,BYPRC,BZPRC)
C
C RETURNS FIELD COMPONENTS FROM A MODEL RING CURRENT, INCLUDING ITS SYMMETRIC PART
C AND A PARTIAL RING CURRENT, CLOSED VIA BIRKELAND CURRENTS. BASED ON RESULTS, DESCRIBED
C IN A PAPER "MODELING THE INNER MAGNETOSPHERE: ASYMMETRIC RING CURRENT AND REGION 2
C BIRKELAND CURRENTS REVISITED" (JGR, DEC.2000).
C
C IOPR - A RING CURRENT CALCULATION FLAG (FOR LEAST-SQUARES FITTING ONLY):
C IOPR=0 - BOTH SRC AND PRC FIELDS ARE CALCULATED
C IOPR=1 - SRC ONLY
C IOPR=2 - PRC ONLY
C
C SC_SY & SC_PR ARE SCALE FACTORS FOR THE ABOVE COMPONENTS; TAKING SC<1 OR SC>1 MAKES THE CURRENTS
C SHRINK OR EXPAND, RESPECTIVELY.
C
C PHI IS THE ROTATION ANGLE (RADIANS) OF THE PARTIAL RING CURRENT (MEASURED FROM MIDNIGHT TOWARD DUSK)
C
IMPLICIT REAL*8 (A-H,O-Z)
c
c 1. TRANSFORM TO TILTED COORDINATES (i.e., SM coordinates):
C
CPS=DCOS(PS)
SPS=DSIN(PS)
XT=X*CPS-Z*SPS
ZT=Z*CPS+X*SPS
C
C 2. SCALE THE COORDINATES FOR THE SYMMETRIC AND PARTIAL RC COMPONENTS:
C
XTS=XT/SC_SY ! SYMMETRIC
YTS=Y /SC_SY
ZTS=ZT/SC_SY
XTA=XT/SC_PR ! PARTIAL
YTA=Y /SC_PR
ZTA=ZT/SC_PR
C
C 3. CALCULATE COMPONENTS OF THE TOTAL FIELD IN THE TILTED (SOLAR-MAGNETIC) COORDINATE SYSTEM:
C
C
C 3a. SYMMETRIC FIELD:
C
IF (IOPR.LE.1) CALL RC_SYMM(XTS,YTS,ZTS,BXS,BYS,BZS)
IF (IOPR.EQ.0.OR.IOPR.EQ.2)
* CALL PRC_SYMM(XTA,YTA,ZTA,BXA_S,BYA_S,BZA_S)
C 3b. ROTATE THE SCALED SM COORDINATES BY PHI AROUND ZSM AXIS AND CALCULATE QUADRUPOLE PRC FIELD
C IN THOSE COORDS:
CP=DCOS(PHI)
SP=DSIN(PHI)
XR=XTA*CP-YTA*SP
YR=XTA*SP+YTA*CP
IF (IOPR.EQ.0.OR.IOPR.EQ.2)
* CALL PRC_QUAD(XR,YR,ZTA,BXA_QR,BYA_QR,BZA_Q)
C 3c. TRANSFORM THE QUADRUPOLE FIELD COMPONENTS BACK TO THE SM COORDS:
C
BXA_Q= BXA_QR*CP+BYA_QR*SP
BYA_Q=-BXA_QR*SP+BYA_QR*CP
C 3d. FIND THE TOTAL FIELD OF PRC (SYMM.+QUADR.) IN THE SM COORDS:
C
BXP=BXA_S+BXA_Q
BYP=BYA_S+BYA_Q
BZP=BZA_S+BZA_Q
C
C 4. TRANSFORM THE FIELDS OF BOTH PARTS OF THE RING CURRENT BACK TO THE GSM SYSTEM:
C
BXSRC=BXS*CPS+BZS*SPS ! SYMMETRIC RC
BYSRC=BYS
BZSRC=BZS*CPS-BXS*SPS
C
BXPRC=BXP*CPS+BZP*SPS ! PARTIAL RC
BYPRC=BYP
BZPRC=BZP*CPS-BXP*SPS
C
RETURN
END
C
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
C
SUBROUTINE RC_SYMM (X,Y,Z,BX,BY,BZ)
IMPLICIT REAL * 8 (A - H, O - Z)
DATA DS,DC/1.D-2,0.99994999875D0/, D/1.D-4/,DRD/5.D3/ ! DS=SIN(THETA) AT THE BOUNDARY OF THE LINEARITY
REGION; DC=SQRT(1-DS**2); DRD=1/(2*D)
RHO2=X**2+Y**2
R2=RHO2+Z**2
R=DSQRT(R2)
RP=R+D
RM=R-D
SINT=DSQRT(RHO2)/R
COST=Z/R
IF (SINT.LT.DS) THEN ! TOO CLOSE TO THE Z-AXIS; USING A LINEAR APPROXIMATION A_PHI~SINT,
C TO AVOID THE SINGULARITY PROBLEM
A=AP(R,DS,DC)/DS
DARDR=(RP*AP(RP,DS,DC)-RM*AP(RM,DS,DC))*DRD
FXY=Z*(2.D0*A-DARDR)/(R*R2)
BX=FXY*X
BY=FXY*Y
BZ=(2.D0*A*COST**2+DARDR*SINT**2)/R
ELSE
THETA=DATAN2(SINT,COST)
TP=THETA+D
TM=THETA-D
SINTP=DSIN(TP)
SINTM=DSIN(TM)
COSTP=DCOS(TP)
COSTM=DCOS(TM)
BR=(SINTP*AP(R,SINTP,COSTP)-SINTM*AP(R,SINTM,COSTM))
* /(R*SINT)*DRD
BT=(RM*AP(RM,SINT,COST)-RP*AP(RP,SINT,COST))/R*DRD
FXY=(BR+BT*COST/SINT)/R
BX=FXY*X
BY=FXY*Y
BZ=BR*COST-BT*SINT
ENDIF
RETURN
END
c
c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
C
DOUBLE PRECISION FUNCTION AP(R,SINT,COST)
C
C Calculates azimuthal component of the vector potential of the symmetric
c part of the model ring current.
C
IMPLICIT REAL * 8 (A - H, O - Z)
LOGICAL PROX ! INDICATES WHETHER WE ARE TOO CLOSE TO THE AXIS OF SYMMETRY, WHERE THE INVERSION
C OF DIPOLAR COORDINATES BECOMES INACCURATE
DATA A1,A2,RRC1,DD1,RRC2,DD2,P1,R1,DR1,DLA1,P2,R2,DR2,DLA2,P3,
*R3,DR3/
*-456.5289941,375.9055332,4.274684950,2.439528329,3.367557287, ! CORRECTED VALUES
*3.146382545,-0.2291904607,3.746064740,1.508802177,0.5873525737, ! (UPDATED 04/20/06 (NB#9, P.37))
*0.1556236119,4.993638842,3.324180497,0.4368407663,0.1855957207,
*2.969226745,2.243367377/
PROX=.FALSE.
SINT1=SINT
COST1=COST
IF (SINT1.LT.1.D-2) THEN ! TOO CLOSE TO Z-AXIS; USE LINEAR INTERPOLATION BETWEEN SINT=0 & SINT=0.01
SINT1=1.D-2
COST1=.99994999875
PROX=.TRUE.
ENDIF
ALPHA=SINT1**2/R ! R,THETA -> ALPHA,GAMMA
GAMMA=COST1/R**2
ARG1=-((R-R1)/DR1)**2-(COST1/DLA1)**2
ARG2=-((R-R2)/DR2)**2-(COST1/DLA2)**2
ARG3=-((R-R3)/DR3)**2
IF (ARG1.LT.-500.D0) THEN ! TO PREVENT "FLOATING UNDERFLOW" CRASHES
DEXP1=0.D0
ELSE
DEXP1=DEXP(ARG1)
ENDIF
IF (ARG2.LT.-500.D0) THEN
DEXP2=0.D0
ELSE
DEXP2=DEXP(ARG2)
ENDIF
IF (ARG3.LT.-500.D0) THEN
DEXP3=0.D0
ELSE
DEXP3=DEXP(ARG3)
ENDIF
ALPHA_S=ALPHA*(1.D0+P1*DEXP1+P2*DEXP2+P3*DEXP3) ! ALPHA -> ALPHA_S (DEFORMED)
GAMMA_S=GAMMA
GAMMAS2=GAMMA_S**2
ALSQH=ALPHA_S**2/2.D0 ! ALPHA_S,GAMMA_S -> RS,SINTS,COSTS
F=64.D0/27.D0*GAMMAS2+ALSQH**2
Q=(DSQRT(F)+ALSQH)**(1.D0/3.D0)
C=Q-4.D0*GAMMAS2**(1.D0/3.D0)/(3.D0*Q)
IF (C.LT.0.D0) C=0.D0
G=DSQRT(C**2+4.D0*GAMMAS2**(1.D0/3.D0))
RS=4.D0/((DSQRT(2.D0*G-C)+DSQRT(C))*(G+C))
COSTS=GAMMA_S*RS**2
SINTS=DSQRT(1.D0-COSTS**2)
RHOS=RS*SINTS
RHOS2=RHOS**2
ZS=RS*COSTS
C
c 1st loop:
P=(RRC1+RHOS)**2+ZS**2+DD1**2
XK2=4.D0*RRC1*RHOS/P
XK=SQRT(XK2)
XKRHO12=XK*SQRT(RHOS)
C
XK2S=1.D0-XK2
DL=DLOG(1.D0/XK2S)
ELK=1.38629436112d0+XK2S*(0.09666344259D0+XK2S*(0.03590092383+
* XK2S*(0.03742563713+XK2S*0.01451196212))) +DL*
* (0.5D0+XK2S*(0.12498593597D0+XK2S*(0.06880248576D0+
* XK2S*(0.03328355346D0+XK2S*0.00441787012D0))))
ELE=1.D0+XK2S*(0.44325141463D0+XK2S*(0.0626060122D0+XK2S*
* (0.04757383546D0+XK2S*0.01736506451D0))) +DL*
* XK2S*(0.2499836831D0+XK2S*(0.09200180037D0+XK2S*
* (0.04069697526D0+XK2S*0.00526449639D0)))
C
APHI1=((1.D0-XK2*0.5D0)*ELK-ELE)/XKRHO12
c
c 2nd loop:
P=(RRC2+RHOS)**2+ZS**2+DD2**2
XK2=4.D0*RRC2*RHOS/P
XK=SQRT(XK2)
XKRHO12=XK*SQRT(RHOS)
C
XK2S=1.D0-XK2
DL=DLOG(1.D0/XK2S)
ELK=1.38629436112d0+XK2S*(0.09666344259D0+XK2S*(0.03590092383+
* XK2S*(0.03742563713+XK2S*0.01451196212))) +DL*
* (0.5D0+XK2S*(0.12498593597D0+XK2S*(0.06880248576D0+
* XK2S*(0.03328355346D0+XK2S*0.00441787012D0))))
ELE=1.D0+XK2S*(0.44325141463D0+XK2S*(0.0626060122D0+XK2S*
* (0.04757383546D0+XK2S*0.01736506451D0))) +DL*
* XK2S*(0.2499836831D0+XK2S*(0.09200180037D0+XK2S*
* (0.04069697526D0+XK2S*0.00526449639D0)))
C
APHI2=((1.D0-XK2*0.5D0)*ELK-ELE)/XKRHO12
AP=A1*APHI1+A2*APHI2
IF (PROX) AP=AP*SINT/SINT1 ! LINEAR INTERPOLATION, IF TOO CLOSE TO THE Z-AXIS
C
RETURN
END
c
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
SUBROUTINE PRC_SYMM (X,Y,Z,BX,BY,BZ)
IMPLICIT REAL * 8 (A - H, O - Z)
DATA DS,DC/1.D-2,0.99994999875D0/, D/1.D-4/,DRD/5.D3/ ! DS=SIN(THETA) AT THE BOUNDARY OF THE LINEARITY
REGION; DC=SQRT(1-DS**2); DRD=1/(2*D)
RHO2=X**2+Y**2
R2=RHO2+Z**2
R=DSQRT(R2)
RP=R+D
RM=R-D
SINT=DSQRT(RHO2)/R
COST=Z/R
IF (SINT.LT.DS) THEN ! TOO CLOSE TO THE Z-AXIS; USING A LINEAR APPROXIMATION A_PHI~SINT,
C TO AVOID THE SINGULARITY PROBLEM
A=APPRC(R,DS,DC)/DS
DARDR=(RP*APPRC(RP,DS,DC)-RM*APPRC(RM,DS,DC))*DRD
FXY=Z*(2.D0*A-DARDR)/(R*R2)
BX=FXY*X
BY=FXY*Y
BZ=(2.D0*A*COST**2+DARDR*SINT**2)/R
ELSE
THETA=DATAN2(SINT,COST)
TP=THETA+D
TM=THETA-D
SINTP=DSIN(TP)
SINTM=DSIN(TM)
COSTP=DCOS(TP)
COSTM=DCOS(TM)
BR=(SINTP*APPRC(R,SINTP,COSTP)-SINTM*APPRC(R,SINTM,COSTM))
* /(R*SINT)*DRD
BT=(RM*APPRC(RM,SINT,COST)-RP*APPRC(RP,SINT,COST))/R*DRD
FXY=(BR+BT*COST/SINT)/R
BX=FXY*X
BY=FXY*Y
BZ=BR*COST-BT*SINT
ENDIF
RETURN
END
c
c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
C
C
DOUBLE PRECISION FUNCTION APPRC(R,SINT,COST)
C
C Calculates azimuthal component of the vector potential of the symmetric
c part of the model PARTIAL ring current.
C
IMPLICIT REAL * 8 (A - H, O - Z)
LOGICAL PROX
DATA A1,A2,RRC1,DD1,RRC2,DD2,P1,ALPHA1,DAL1,BETA1,DG1,P2,ALPHA2,
* DAL2,BETA2,DG2,BETA3,P3,ALPHA3,DAL3,BETA4,DG3,BETA5,Q0,Q1,ALPHA4,
* DAL4,DG4,Q2,ALPHA5,DAL5,DG5,BETA6,BETA7
* /-80.11202281,12.58246758,6.560486035,1.930711037,3.827208119,
*.7789990504,.3058309043,.1817139853,.1257532909,3.422509402,
*.04742939676,-4.800458958,-.02845643596,.2188114228,2.545944574,
*.00813272793,.35868244,103.1601001,-.00764731187,.1046487459,
*2.958863546,.01172314188,.4382872938,.01134908150,14.51339943,
*.2647095287,.07091230197,.01512963586,6.861329631,.1677400816,
*.04433648846,.05553741389,.7665599464,.7277854652/
PROX=.FALSE.
SINT1=SINT
COST1=COST
IF (SINT1.LT.1.D-2) THEN ! TOO CLOSE TO Z-AXIS; USE LINEAR INTERPOLATION BETWEEN SINT=0 & SINT=0.01
SINT1=1.D-2
COST1=.99994999875
PROX=.TRUE.
ENDIF
ALPHA=SINT1**2/R ! R,THETA -> ALPHA,GAMMA
GAMMA=COST1/R**2
ARG1=-(GAMMA/DG1)**2
ARG2=-((ALPHA-ALPHA4)/DAL4)**2-(GAMMA/DG4)**2
IF (ARG1.LT.-500.D0) THEN ! TO PREVENT "FLOATING UNDERFLOW" CRASHES
DEXP1=0.D0
ELSE
DEXP1=DEXP(ARG1)
ENDIF
IF (ARG2.LT.-500.D0) THEN ! TO PREVENT "FLOATING UNDERFLOW" CRASHES
DEXP2=0.D0
ELSE
DEXP2=DEXP(ARG2)
ENDIF
ALPHA_S=ALPHA*(1.D0+P1/(1.D0+((ALPHA-ALPHA1)/DAL1)**2)**BETA1
* *DEXP1+P2*(ALPHA-ALPHA2)/(1.D0+((ALPHA-ALPHA2)/DAL2)**2)**BETA2
*/(1.D0+(GAMMA/DG2)**2)**BETA3
*+P3*(ALPHA-ALPHA3)**2/(1.D0+((ALPHA-ALPHA3)/DAL3)**2)**BETA4
*/(1.D0+(GAMMA/DG3)**2)**BETA5) ! ALPHA -> ALPHA_S (DEFORMED)
GAMMA_S=GAMMA*(1.D0+Q0+Q1*(ALPHA-ALPHA4)*DEXP2 ! GAMMA -> GAMMA_ (DEFORMED)
* +Q2*(ALPHA-ALPHA5)/(1.D0+((ALPHA-ALPHA5)/DAL5)**2)**BETA6
* /(1.D0+(GAMMA/DG5)**2)**BETA7)
GAMMAS2=GAMMA_S**2
ALSQH=ALPHA_S**2/2.D0 ! ALPHA_S,GAMMA_S -> RS,SINTS,COSTS
F=64.D0/27.D0*GAMMAS2+ALSQH**2
Q=(DSQRT(F)+ALSQH)**(1.D0/3.D0)
C=Q-4.D0*GAMMAS2**(1.D0/3.D0)/(3.D0*Q)
IF (C.LT.0.D0) C=0.D0
G=DSQRT(C**2+4.D0*GAMMAS2**(1.D0/3.D0))
RS=4.D0/((DSQRT(2.D0*G-C)+DSQRT(C))*(G+C))
COSTS=GAMMA_S*RS**2
SINTS=DSQRT(1.D0-COSTS**2)
RHOS=RS*SINTS
RHOS2=RHOS**2
ZS=RS*COSTS
C
c 1st loop:
P=(RRC1+RHOS)**2+ZS**2+DD1**2
XK2=4.D0*RRC1*RHOS/P
XK=SQRT(XK2)
XKRHO12=XK*SQRT(RHOS)
C
XK2S=1.D0-XK2
DL=DLOG(1.D0/XK2S)
ELK=1.38629436112d0+XK2S*(0.09666344259D0+XK2S*(0.03590092383+
* XK2S*(0.03742563713+XK2S*0.01451196212))) +DL*
* (0.5D0+XK2S*(0.12498593597D0+XK2S*(0.06880248576D0+
* XK2S*(0.03328355346D0+XK2S*0.00441787012D0))))
ELE=1.D0+XK2S*(0.44325141463D0+XK2S*(0.0626060122D0+XK2S*
* (0.04757383546D0+XK2S*0.01736506451D0))) +DL*
* XK2S*(0.2499836831D0+XK2S*(0.09200180037D0+XK2S*
* (0.04069697526D0+XK2S*0.00526449639D0)))
C
APHI1=((1.D0-XK2*0.5D0)*ELK-ELE)/XKRHO12
c
c 2nd loop:
P=(RRC2+RHOS)**2+ZS**2+DD2**2
XK2=4.D0*RRC2*RHOS/P
XK=SQRT(XK2)
XKRHO12=XK*SQRT(RHOS)
C
XK2S=1.D0-XK2
DL=DLOG(1.D0/XK2S)
ELK=1.38629436112d0+XK2S*(0.09666344259D0+XK2S*(0.03590092383+
* XK2S*(0.03742563713+XK2S*0.01451196212))) +DL*
* (0.5D0+XK2S*(0.12498593597D0+XK2S*(0.06880248576D0+
* XK2S*(0.03328355346D0+XK2S*0.00441787012D0))))
ELE=1.D0+XK2S*(0.44325141463D0+XK2S*(0.0626060122D0+XK2S*
* (0.04757383546D0+XK2S*0.01736506451D0))) +DL*
* XK2S*(0.2499836831D0+XK2S*(0.09200180037D0+XK2S*
* (0.04069697526D0+XK2S*0.00526449639D0)))
C
APHI2=((1.D0-XK2*0.5D0)*ELK-ELE)/XKRHO12
APPRC=A1*APHI1+A2*APHI2
IF (PROX) APPRC=APPRC*SINT/SINT1 ! LINEAR INTERPOLATION, IF TOO CLOSE TO THE Z-AXIS
C
RETURN
END
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C
SUBROUTINE PRC_QUAD (X,Y,Z,BX,BY,BZ)
C
IMPLICIT REAL * 8 (A - H, O - Z)
DATA D,DD/1.D-4,2.D-4/, DS/1.D-2/,DC/0.99994999875D0/
RHO2=X**2+Y**2
R=DSQRT(RHO2+Z**2)
RHO=DSQRT(RHO2)
SINT=RHO/R
COST=Z/R
RP=R+D
RM=R-D
IF (SINT.GT.DS) THEN
CPHI=X/RHO
SPHI=Y/RHO
BR=BR_PRC_Q(R,SINT,COST)
BT=BT_PRC_Q(R,SINT,COST)
DBRR=(BR_PRC_Q(RP,SINT,COST)-BR_PRC_Q(RM,SINT,COST))/DD
THETA=DATAN2(SINT,COST)
TP=THETA+D
TM=THETA-D
SINTP=DSIN(TP)
COSTP=DCOS(TP)
SINTM=DSIN(TM)
COSTM=DCOS(TM)
DBTT=(BT_PRC_Q(R,SINTP,COSTP)-BT_PRC_Q(R,SINTM,COSTM))/DD
BX=SINT*(BR+(BR+R*DBRR+DBTT)*SPHI**2)+COST*BT
BY=-SINT*SPHI*CPHI*(BR+R*DBRR+DBTT)
BZ=(BR*COST-BT*SINT)*CPHI
ELSE
ST=DS
CT=DC
IF (Z.LT.0.D0) CT=-DC
THETA=DATAN2(ST,CT)
TP=THETA+D
TM=THETA-D
SINTP=DSIN(TP)
COSTP=DCOS(TP)
SINTM=DSIN(TM)
COSTM=DCOS(TM)
BR=BR_PRC_Q(R,ST,CT)
BT=BT_PRC_Q(R,ST,CT)
DBRR=(BR_PRC_Q(RP,ST,CT)-BR_PRC_Q(RM,ST,CT))/DD
DBTT=(BT_PRC_Q(R,SINTP,COSTP)-BT_PRC_Q(R,SINTM,COSTM))/DD
FCXY=R*DBRR+DBTT
BX=(BR*(X**2+2.D0*Y**2)+FCXY*Y**2)/(R*ST)**2+BT*COST
BY=-(BR+FCXY)*X*Y/(R*ST)**2
BZ=(BR*COST/ST-BT)*X/R
ENDIF
RETURN
END
c
c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
C
DOUBLE PRECISION FUNCTION BR_PRC_Q (R,SINT,COST)
C
Calculates the radial component of the "quadrupole" part of the model partial ring current.
C
IMPLICIT REAL * 8 (A - H, O - Z)
DATA A1,A2,A3,A4,A5,A6,A7,A8,A9,A10,A11,A12,A13,A14,A15,A16,A17, ! ALL LINEAR PARAMETERS HERE
* A18,XK1,AL1,DAL1,B1,BE1,XK2,AL2,DAL2,B2,BE2,XK3,XK4,AL3,DAL3,B3, ! WERE MULTIPLIED BY 0.1,
* BE3,AL4,DAL4,DG1,AL5,DAL5,DG2,C1,C2,C3,AL6,DAL6,DRM/-21.2666329, ! SO THAT THEY CORRESPOND TO P_0=1 nPa,
*32.24527521,-6.062894078,7.515660734,233.7341288,-227.1195714, ! RATHER THAN THE ORIGINAL VALUE OF 10 nPa
*8.483233889,16.80642754,-24.63534184,9.067120578,-1.052686913, ! ASSUMED IN THE BIOT-SAVART INTEGRAL
*-12.08384538,18.61969572,-12.71686069,47017.35679,-50646.71204,
*7746.058231,1.531069371,2.318824273,.1417519429,.6388013110E-02,
*5.303934488,4.213397467,.7955534018,.1401142771,.2306094179E-01,
*3.462235072,2.568743010,3.477425908,1.922155110,.1485233485,
*.2319676273E-01,7.830223587,8.492933868,.1295221828,.01753008801,
*.01125504083,.1811846095,.04841237481,.01981805097,6.557801891,
*6.348576071,5.744436687,.2265212965,.1301957209,.5654023158/
SINT2=SINT**2
COST2=COST**2
SC=SINT*COST
ALPHA=SINT2/R
GAMMA=COST/R**2
CALL FFS(ALPHA,AL1,DAL1,F,FA,FS)
D1=SC*F**XK1/((R/B1)**BE1+1.D0)
D2=D1*COST2
CALL FFS(ALPHA,AL2,DAL2,F,FA,FS)
D3=SC*FS**XK2/((R/B2)**BE2+1.D0)
D4=D3*COST2
CALL FFS(ALPHA,AL3,DAL3,F,FA,FS)
D5=SC*(ALPHA**XK3)*(FS**XK4)/((R/B3)**BE3+1.D0)
D6=D5*COST2
ARGA=((ALPHA-AL4)/DAL4)**2+1.D0
ARGG=1.D0+(GAMMA/DG1)**2
D7=SC/ARGA/ARGG
D8=D7/ARGA
D9=D8/ARGA
D10=D9/ARGA
ARGA=((ALPHA-AL5)/DAL5)**2+1.D0
ARGG=1.D0+(GAMMA/DG2)**2
D11=SC/ARGA/ARGG
D12=D11/ARGA
D13=D12/ARGA
D14=D13/ARGA
D15=SC/(R**4+C1**4)
D16=SC/(R**4+C2**4)*COST2
D17=SC/(R**4+C3**4)*COST2**2
CALL FFS(ALPHA,AL6,DAL6,F,FA,FS)
D18=SC*FS/(1.D0+((R-1.2D0)/DRM)**2)
BR_PRC_Q=A1*D1+A2*D2+A3*D3+A4*D4+A5*D5+A6*D6+A7*D7+A8*D8+A9*D9+
* A10*D10+A11*D11+A12*D12+A13*D13+A14*D14+A15*D15+A16*D16+A17*D17+
* A18*D18
C
RETURN
END
c
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C
DOUBLE PRECISION FUNCTION BT_PRC_Q (R,SINT,COST)
C
Calculates the Theta component of the "quadrupole" part of the model partial ring current.
C
IMPLICIT REAL * 8 (A - H, O - Z)
DATA A1,A2,A3,A4,A5,A6,A7,A8,A9,A10,A11,A12,A13,A14,A15,A16,A17, ! ALL LINEAR PARAMETERS HERE
*XK1,AL1,DAL1,B1,BE1,XK2,AL2,DAL2,BE2,XK3,XK4,AL3,DAL3,B3,BE3,AL4, ! WERE MULTIPLIED BY 0.1,
*DAL4,DG1,AL5,DAL5,DG2,C1,C2,C3/12.74640393,-7.516393516, ! SO THAT THEY CORRESPOND TO P_0=1 nPa,
*-5.476233865,3.212704645,-59.10926169,46.62198189,-.01644280062, ! RATHER THAN THE ORIGINAL VALUE OF 10 nPa
*.1234229112,-.08579198697,.01321366966,.8970494003,9.136186247, ! ASSUMED IN THE BIOT-SAVART INTEGRAL
*-38.19301215,21.73775846,-410.0783424,-69.90832690,-848.8543440,
*1.243288286,.2071721360,.05030555417,7.471332374,3.180533613,
*1.376743507,.1568504222,.02092910682,1.985148197,.3157139940,
*1.056309517,.1701395257,.1019870070,6.293740981,5.671824276,
*.1280772299,.02189060799,.01040696080,.1648265607,.04701592613,
*.01526400086,12.88384229,3.361775101,23.44173897/
SINT2=SINT**2
COST2=COST**2
SC=SINT*COST
ALPHA=SINT2/R
GAMMA=COST/R**2
CALL FFS(ALPHA,AL1,DAL1,F,FA,FS)
D1=F**XK1/((R/B1)**BE1+1.D0)
D2=D1*COST2
CALL FFS(ALPHA,AL2,DAL2,F,FA,FS)
D3=FA**XK2/R**BE2
D4=D3*COST2
CALL FFS(ALPHA,AL3,DAL3,F,FA,FS)
D5=FS**XK3*ALPHA**XK4/((R/B3)**BE3+1.D0)
D6=D5*COST2
CALL FFS(GAMMA,0.D0,DG1,F,FA,FS)
FCC=(1.D0+((ALPHA-AL4)/DAL4)**2)
D7 =1.D0/FCC*FS
D8 =D7/FCC
D9 =D8/FCC
D10=D9/FCC
ARG=1.D0+((ALPHA-AL5)/DAL5)**2
D11=1.D0/ARG/(1.D0+(GAMMA/DG2)**2)
D12=D11/ARG
D13=D12/ARG
D14=D13/ARG
D15=1.D0/(R**4+C1**2)
D16=COST2/(R**4+C2**2)
D17=COST2**2/(R**4+C3**2)
C
BT_PRC_Q=A1*D1+A2*D2+A3*D3+A4*D4+A5*D5+A6*D6+A7*D7+A8*D8+A9*D9+
* A10*D10+A11*D11+A12*D12+A13*D13+A14*D14+A15*D15+A16*D16+A17*D17
C
RETURN
END
c
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
SUBROUTINE FFS(A,A0,DA,F,FA,FS)
IMPLICIT REAL * 8 (A - H, O - Z)
SQ1=DSQRT((A+A0)**2+DA**2)
SQ2=DSQRT((A-A0)**2+DA**2)
FA=2.D0/(SQ1+SQ2)
F=FA*A
FS=0.5D0*(SQ1+SQ2)/(SQ1*SQ2)*(1.D0-F*F)
RETURN
END
C
C||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
C
C-------------------------------------------------------------------------
C
C
SUBROUTINE RC_SHIELD (A,PS,X_SC,X,Y,Z,BX,BY,BZ)
C
IMPLICIT REAL * 8 (A - H, O - Z)
DIMENSION A(86)
C
FAC_SC=(X_SC+1.D0)**3
C
CPS=DCOS(PS)
SPS=DSIN(PS)
S3PS=2.D0*CPS
C
PST1=PS*A(85)
PST2=PS*A(86)
ST1=DSIN(PST1)
CT1=DCOS(PST1)
ST2=DSIN(PST2)
CT2=DCOS(PST2)
X1=X*CT1-Z*ST1
Z1=X*ST1+Z*CT1
X2=X*CT2-Z*ST2
Z2=X*ST2+Z*CT2
C
L=0
GX=0.D0
GY=0.D0
GZ=0.D0
C
DO 1 M=1,2 ! M=1 IS FOR THE 1ST SUM ("PERP." SYMMETRY)
C AND M=2 IS FOR THE SECOND SUM ("PARALL." SYMMETRY)
DO 2 I=1,3
P=A(72+I)
Q=A(78+I)
CYPI=DCOS(Y/P)
CYQI=DCOS(Y/Q)
SYPI=DSIN(Y/P)
SYQI=DSIN(Y/Q)
C
DO 3 K=1,3
R=A(75+K)
S=A(81+K)
SZRK=DSIN(Z1/R)
CZSK=DCOS(Z2/S)
CZRK=DCOS(Z1/R)
SZSK=DSIN(Z2/S)
SQPR=DSQRT(1.D0/P**2+1.D0/R**2)
SQQS=DSQRT(1.D0/Q**2+1.D0/S**2)
EPR=DEXP(X1*SQPR)
EQS=DEXP(X2*SQQS)
C
DO 4 N=1,2 ! N=1 IS FOR THE FIRST PART OF EACH COEFFICIENT
C AND N=2 IS FOR THE SECOND ONE
DO 5 NN=1,2 ! NN = 1,2 FURTHER SPLITS THE COEFFICIENTS INTO 2 PARTS,
C TO TAKE INTO ACCOUNT THE SCALE FACTOR DEPENDENCE
IF (M.EQ.1) THEN
FX=-SQPR*EPR*CYPI*SZRK *FAC_SC
FY=EPR*SYPI*SZRK/P *FAC_SC
FZ=-EPR*CYPI*CZRK/R *FAC_SC
IF (N.EQ.1) THEN
IF (NN.EQ.1) THEN
HX=FX
HY=FY
HZ=FZ
ELSE
HX=FX*X_SC
HY=FY*X_SC
HZ=FZ*X_SC
ENDIF
ELSE
IF (NN.EQ.1) THEN
HX=FX*CPS
HY=FY*CPS
HZ=FZ*CPS
ELSE
HX=FX*CPS*X_SC
HY=FY*CPS*X_SC
HZ=FZ*CPS*X_SC
ENDIF
ENDIF
ELSE ! M.EQ.2
FX=-SPS*SQQS*EQS*CYQI*CZSK *FAC_SC
FY=SPS/Q*EQS*SYQI*CZSK *FAC_SC
FZ=SPS/S*EQS*CYQI*SZSK *FAC_SC
IF (N.EQ.1) THEN
IF (NN.EQ.1) THEN
HX=FX
HY=FY
HZ=FZ
ELSE
HX=FX*X_SC
HY=FY*X_SC
HZ=FZ*X_SC
ENDIF
ELSE
IF (NN.EQ.1) THEN
HX=FX*S3PS
HY=FY*S3PS
HZ=FZ*S3PS
ELSE
HX=FX*S3PS*X_SC
HY=FY*S3PS*X_SC
HZ=FZ*S3PS*X_SC
ENDIF
ENDIF
ENDIF
L=L+1
IF (M.EQ.1) THEN
HXR=HX*CT1+HZ*ST1
HZR=-HX*ST1+HZ*CT1
ELSE
HXR=HX*CT2+HZ*ST2
HZR=-HX*ST2+HZ*CT2
ENDIF
GX=GX+HXR*A(L)
GY=GY+HY *A(L)
GZ=GZ+HZR*A(L)
5 CONTINUE
4 CONTINUE
3 CONTINUE
2 CONTINUE
1 CONTINUE
BX=GX
BY=GY
BZ=GZ
RETURN
END
C
c===========================================================================
c
SUBROUTINE DIPOLE (PS,X,Y,Z,BX,BY,BZ)
C
C A DOUBLE PRECISION ROUTINE
C
C CALCULATES GSM COMPONENTS OF A GEODIPOLE FIELD WITH THE DIPOLE MOMENT
C CORRESPONDING TO THE EPOCH OF 2000.
C
C----INPUT PARAMETERS:
C PS - GEODIPOLE TILT ANGLE IN RADIANS,
C X,Y,Z - GSM COORDINATES IN RE (1 RE = 6371.2 km)
C
C----OUTPUT PARAMETERS:
C BX,BY,BZ - FIELD COMPONENTS IN GSM SYSTEM, IN NANOTESLA.
C
IMPLICIT REAL*8 (A-H,O-Z)
SPS=DSIN(PS)
CPS=DCOS(PS)
P=X**2
U=Z**2
V=3.D0*Z*X
T=Y**2
Q=30115.D0/DSQRT(P+T+U)**5
BX=Q*((T+U-2.D0*P)*SPS-V*CPS)
BY=-3.D0*Y*Q*(X*SPS+Z*CPS)
BZ=Q*((P+T-2.D0*U)*CPS-V*SPS)
RETURN
END
C
C(((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((