c C---------------------------------------------------------------------- c SUBROUTINE T96_01 (IOPT,PARMOD,PS,X,Y,Z,BX,BY,BZ) C c RELEASE DATE OF THIS VERSION: JUNE 22, 1996. C LAST UPDATE: MAY 01, 2006: IN THE S/R DIPOLE, SPS AND CPS WERE ADDED IN THE SAVE STATEMENT C---------------------------------------------------------------------- C C WITH TWO CORRECTIONS, SUGGESTED BY T.SOTIRELIS COMMENTS (APR.7, 1997) C C (1) A "STRAY " CLOSING PARENTHESIS WAS REMOVED IN THE S/R R2_BIRK C (2) A 0/0 PROBLEM ON THE Z-AXIS WAS SIDESTEPPED (LINES 44-46 OF THE c DOUBLE PRECISION FUNCTION XKSI) c-------------------------------------------------------------------- C DATA-BASED MODEL CALIBRATED BY (1) SOLAR WIND PRESSURE PDYN (NANOPASCALS), C (2) DST (NANOTESLA), (3) BYIMF, AND (4) BZIMF (NANOTESLA). c THESE INPUT PARAMETERS SHOULD BE PLACED IN THE FIRST 4 ELEMENTS c OF THE ARRAY PARMOD(10). C C THE REST OF THE INPUT VARIABLES ARE: THE GEODIPOLE TILT ANGLE PS (RADIANS), C AND X,Y,Z - GSM POSITION (RE) C c IOPT IS JUST A DUMMY INPUT PARAMETER, NECESSARY TO MAKE THIS SUBROUTINE C COMPATIBLE WITH THE NEW RELEASE (APRIL 1996) OF THE TRACING SOFTWARE C PACKAGE (GEOPACK). IOPT VALUE DOES NOT AFFECT THE OUTPUT FIELD. 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. 1995, 1996, Nikolai A. Tsyganenko, Hughes STX, Code 695, NASA GSFC c Greenbelt, MD 20771, USA c C REFERENCES: C C (1) N.A. TSYGANENKO AND D.P. STERN, A NEW-GENERATION GLOBAL C MAGNETOSPHERE FIELD MODEL , BASED ON SPACECRAFT MAGNETOMETER DATA, C ISTP NEWSLETTER, V.6, NO.1, P.21, FEB.1996. C c (2) N.A.TSYGANENKO, MODELING THE EARTH S MAGNETOSPHERIC C MAGNETIC FIELD CONFINED WITHIN A REALISTIC MAGNETOPAUSE, C J.GEOPHYS.RES., V.100, P. 5599, 1995. C C (3) N.A. TSYGANENKO AND M.PEREDO, ANALYTICAL MODELS OF THE C MAGNETIC FIELD OF DISK-SHAPED CURRENT SHEETS, J.GEOPHYS.RES., C V.99, P. 199, 1994. C c---------------------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) REAL*8 PDYN,DST,BYIMF,BZIMF,PS,X,Y,Z,BX,BY,BZ,QX,QY,QZ,PARMOD(10), * A(9) c DATA PDYN0,EPS10 /2.,3630.7/ C DATA A/1.162,22.344,18.50,2.602,6.903,5.287,0.5790,0.4462,0.7850/ C DATA AM0,S0,X00,DSIG/70.,1.08,5.48,0.005/ DATA DELIMFX,DELIMFY /20.,10./ C PDYN=PARMOD(1) DST=PARMOD(2) BYIMF=PARMOD(3) BZIMF=PARMOD(4) C SPS=SIN(PS) PPS=PS C DEPR=0.8*DST-13.*SQRT(PDYN) ! DEPR is an estimate of total near-Earth c depression, based on DST and Pdyn c (usually, DEPR < 0 ) C C CALCULATE THE IMF-RELATED QUANTITIES: C Bt=SQRT(BYIMF**2+BZIMF**2) IF (BYIMF.EQ.0..AND.BZIMF.EQ.0.) THEN THETA=0. GOTO 1 ENDIF C THETA=ATAN2(BYIMF,BZIMF) IF (THETA.LE.0.D0) THETA=THETA+6.2831853 1 CT=COS(THETA) ST=SIN(THETA) EPS=718.5*SQRT(Pdyn)*Bt*SIN(THETA/2.) C FACTEPS=EPS/EPS10-1. FACTPD=SQRT(PDYN/PDYN0)-1. C RCAMPL=-A(1)*DEPR ! RCAMPL is the amplitude of the ring current c (positive and equal to abs.value of RC depression at origin) C TAMPL2=A(2)+A(3)*FACTPD+A(4)*FACTEPS TAMPL3=A(5)+A(6)*FACTPD B1AMPL=A(7)+A(8)*FACTEPS B2AMPL=20.*B1AMPL ! IT IS EQUIVALENT TO ASSUMING THAT THE TOTAL CURRENT C IN THE REGION 2 SYSTEM IS 40% OF THAT IN REGION 1 RECONN=A(9) C XAPPA=(PDYN/PDYN0)**0.14 XAPPA3=XAPPA**3 YS=Y*CT-Z*ST ZS=Z*CT+Y*ST C FACTIMF=EXP(X/DELIMFX-(YS/DELIMFY)**2) C C CALCULATE THE "IMF" COMPONENTS OUTSIDE THE LAYER (HENCE BEGIN WITH "O") C OIMFX=0. OIMFY=RECONN*BYIMF*FACTIMF OIMFZ=RECONN*BZIMF*FACTIMF C RIMFAMPL=RECONN*Bt C PPS=PS XX=X*XAPPA YY=Y*XAPPA ZZ=Z*XAPPA C C SCALE AND CALCULATE THE MAGNETOPAUSE PARAMETERS FOR THE INTERPOLATION ACROSS C THE BOUNDARY LAYER (THE COORDINATES XX,YY,ZZ ARE ALREADY SCALED) C X0=X00/XAPPA AM=AM0/XAPPA RHO2=Y**2+Z**2 ASQ=AM**2 XMXM=AM+X-X0 IF (XMXM.LT.0.) XMXM=0. ! THE BOUNDARY IS A CYLINDER TAILWARD OF X=X0-AM AXX0=XMXM**2 ARO=ASQ+RHO2 SIGMA=SQRT((ARO+AXX0+SQRT((ARO+AXX0)**2-4.*ASQ*AXX0))/(2.*ASQ)) C C NOW, THERE ARE THREE POSSIBLE CASES: C (1) INSIDE THE MAGNETOSPHERE 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 IF (SIGMA.LT.S0+DSIG) THEN ! CALCULATE THE T95_06 FIELD (WITH THE C POTENTIAL "PENETRATED" INTERCONNECTION FIELD): CALL DIPSHLD(PPS,XX,YY,ZZ,CFX,CFY,CFZ) CALL TAILRC96(SPS,XX,YY,ZZ,BXRC,BYRC,BZRC,BXT2,BYT2,BZT2, * BXT3,BYT3,BZT3) CALL BIRK1TOT_02(PPS,XX,YY,ZZ,R1X,R1Y,R1Z) CALL BIRK2TOT_02(PPS,XX,YY,ZZ,R2X,R2Y,R2Z) CALL INTERCON(XX,YS*XAPPA,ZS*XAPPA,RIMFX,RIMFYS,RIMFZS) RIMFY=RIMFYS*CT+RIMFZS*ST RIMFZ=RIMFZS*CT-RIMFYS*ST C FX=CFX*XAPPA3+RCAMPL*BXRC +TAMPL2*BXT2+TAMPL3*BXT3 * +B1AMPL*R1X +B2AMPL*R2X +RIMFAMPL*RIMFX FY=CFY*XAPPA3+RCAMPL*BYRC +TAMPL2*BYT2+TAMPL3*BYT3 * +B1AMPL*R1Y +B2AMPL*R2Y +RIMFAMPL*RIMFY FZ=CFZ*XAPPA3+RCAMPL*BZRC +TAMPL2*BZT2+TAMPL3*BZT3 * +B1AMPL*R1Z +B2AMPL*R2Z +RIMFAMPL*RIMFZ C C NOW, LET US CHECK WHETHER WE HAVE THE CASE (1). IF YES - WE ARE DONE: C IF (SIGMA.LT.S0-DSIG) THEN BX=FX BY=FY BZ=FZ 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=(FX+QX)*FINT+OIMFX*FEXT -QX BY=(FY+QY)*FINT+OIMFY*FEXT -QY BZ=(FZ+QZ)*FINT+OIMFZ*FEXT -QZ c ENDIF ! THE CASES (1) AND (2) ARE EXHAUSTED; THE ONLY REMAINING C POSSIBILITY IS NOW THE CASE (3): ELSE CALL DIPOLE(PS,X,Y,Z,QX,QY,QZ) BX=OIMFX-QX BY=OIMFY-QY BZ=OIMFZ-QZ ENDIF C RETURN END C===================================================================== SUBROUTINE DIPSHLD(PS,X,Y,Z,BX,BY,BZ) C C CALCULATES GSM COMPONENTS OF THE EXTERNAL MAGNETIC FIELD DUE TO C SHIELDING OF THE EARTH S DIPOLE ONLY C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A1(12),A2(12) DATA A1 /.24777,-27.003,-.46815,7.0637,-1.5918,-.90317E-01,57.522, * 13.757,2.0100,10.458,4.5798,2.1695/ DATA A2/-.65385,-18.061,-.40457,-5.0995,1.2846,.78231E-01,39.592, * 13.291,1.9970,10.062,4.5140,2.1558/ C CPS=DCOS(PS) SPS=DSIN(PS) CALL CYLHARM(A1,X,Y,Z,HX,HY,HZ) CALL CYLHAR1(A2,X,Y,Z,FX,FY,FZ) C BX=HX*CPS+FX*SPS BY=HY*CPS+FY*SPS BZ=HZ*CPS+FZ*SPS RETURN END C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C THIS CODE YIELDS THE SHIELDING FIELD FOR THE PERPENDICULAR DIPOLE C SUBROUTINE CYLHARM( A, X,Y,Z, BX,BY,BZ) C C C *** N.A. Tsyganenko *** Sept. 14-18, 1993; revised March 16, 1994 *** C C An approximation for the Chapman-Ferraro field by a sum of 6 cylin- c drical harmonics (see pp. 97-113 in the brown GSFC notebook #1) c C Description of parameters: C C A - input vector containing model parameters; C X,Y,Z - input GSM coordinates C BX,BY,BZ - output GSM components of the shielding field C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C The 6 linear parameters A(1)-A(6) are amplitudes of the cylindrical harmonic c terms. c The 6 nonlinear parameters A(7)-A(12) are the corresponding scale lengths C for each term (see GSFC brown notebook). c C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C IMPLICIT REAL * 8 (A - H, O - Z) C DIMENSION A(12) C RHO=DSQRT(Y**2+Z**2) IF (RHO.LT.1.D-8) THEN SINFI=1.D0 COSFI=0.D0 RHO=1.D-8 GOTO 1 ENDIF C SINFI=Z/RHO COSFI=Y/RHO 1 SINFI2=SINFI**2 SI2CO2=SINFI2-COSFI**2 C BX=0.D0 BY=0.D0 BZ=0.D0 C DO 11 I=1,3 DZETA=RHO/A(I+6) XJ0=BES(DZETA,0) XJ1=BES(DZETA,1) XEXP=DEXP(X/A(I+6)) BX=BX-A(I)*XJ1*XEXP*SINFI BY=BY+A(I)*(2.D0*XJ1/DZETA-XJ0)*XEXP*SINFI*COSFI BZ=BZ+A(I)*(XJ1/DZETA*SI2CO2-XJ0*SINFI2)*XEXP 11 CONTINUE c DO 12 I=4,6 DZETA=RHO/A(I+6) XKSI=X/A(I+6) XJ0=BES(DZETA,0) XJ1=BES(DZETA,1) XEXP=DEXP(XKSI) BRHO=(XKSI*XJ0-(DZETA**2+XKSI-1.D0)*XJ1/DZETA)*XEXP*SINFI BPHI=(XJ0+XJ1/DZETA*(XKSI-1.D0))*XEXP*COSFI BX=BX+A(I)*(DZETA*XJ0+XKSI*XJ1)*XEXP*SINFI BY=BY+A(I)*(BRHO*COSFI-BPHI*SINFI) BZ=BZ+A(I)*(BRHO*SINFI+BPHI*COSFI) 12 CONTINUE C c RETURN END C c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C C THIS CODE YIELDS THE SHIELDING FIELD FOR THE PARALLEL DIPOLE C SUBROUTINE CYLHAR1(A, X,Y,Z, BX,BY,BZ) C C C *** N.A. Tsyganenko *** Sept. 14-18, 1993; revised March 16, 1994 *** C C An approximation of the Chapman-Ferraro field by a sum of 6 cylin- c drical harmonics (see pages 97-113 in the brown GSFC notebook #1) c C Description of parameters: C C A - input vector containing model parameters; C X,Y,Z - input GSM coordinates, C BX,BY,BZ - output GSM components of the shielding field C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C C The 6 linear parameters A(1)-A(6) are amplitudes of the cylindrical c harmonic terms. c The 6 nonlinear parameters A(7)-A(12) are the corresponding scale c lengths for each term (see GSFC brown notebook). c C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C IMPLICIT REAL * 8 (A - H, O - Z) C DIMENSION A(12) C RHO=DSQRT(Y**2+Z**2) IF (RHO.LT.1.D-10) THEN SINFI=1.D0 COSFI=0.D0 GOTO 1 ENDIF C SINFI=Z/RHO COSFI=Y/RHO C 1 BX=0.D0 BY=0.D0 BZ=0.D0 C DO 11 I=1,3 DZETA=RHO/A(I+6) XKSI=X/A(I+6) XJ0=BES(DZETA,0) XJ1=BES(DZETA,1) XEXP=DEXP(XKSI) BRHO=XJ1*XEXP BX=BX-A(I)*XJ0*XEXP BY=BY+A(I)*BRHO*COSFI BZ=BZ+A(I)*BRHO*SINFI 11 CONTINUE c DO 12 I=4,6 DZETA=RHO/A(I+6) XKSI=X/A(I+6) XJ0=BES(DZETA,0) XJ1=BES(DZETA,1) XEXP=DEXP(XKSI) BRHO=(DZETA*XJ0+XKSI*XJ1)*XEXP BX=BX+A(I)*(DZETA*XJ1-XJ0*(XKSI+1.D0))*XEXP BY=BY+A(I)*BRHO*COSFI BZ=BZ+A(I)*BRHO*SINFI 12 CONTINUE C RETURN END c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C DOUBLE PRECISION FUNCTION BES(X,K) IMPLICIT REAL*8 (A-H,O-Z) C IF (K.EQ.0) THEN BES=BES0(X) RETURN ENDIF C IF (K.EQ.1) THEN BES=BES1(X) RETURN ENDIF C IF (X.EQ.0.D0) THEN BES=0.D0 RETURN ENDIF C G=2.D0/X IF (X.LE.DFLOAT(K)) GOTO 10 C N=1 XJN=BES1(X) XJNM1=BES0(X) C 1 XJNP1=G*N*XJN-XJNM1 N=N+1 IF (N.LT.K) GOTO 2 BES=XJNP1 RETURN C 2 XJNM1=XJN XJN=XJNP1 GOTO 1 C 10 N=24 XJN=1.D0 XJNP1=0.D0 SUM=0.D0 C 3 IF (MOD(N,2).EQ.0) SUM=SUM+XJN XJNM1=G*N*XJN-XJNP1 N=N-1 C XJNP1=XJN XJN=XJNM1 IF (N.EQ.K) BES=XJN C IF (DABS(XJN).GT.1.D5) THEN XJNP1=XJNP1*1.D-5 XJN=XJN*1.D-5 SUM=SUM*1.D-5 IF (N.LE.K) BES=BES*1.D-5 ENDIF C IF (N.EQ.0) GOTO 4 GOTO 3 C 4 SUM=XJN+2.D0*SUM BES=BES/SUM RETURN END c------------------------------------------------------------------- c DOUBLE PRECISION FUNCTION BES0(X) C IMPLICIT REAL*8 (A-H,O-Z) C IF (DABS(X).LT.3.D0) THEN X32=(X/3.D0)**2 BES0=1.D0-X32*(2.2499997D0-X32*(1.2656208D0-X32* * (0.3163866D0-X32*(0.0444479D0-X32*(0.0039444D0 * -X32*0.00021D0))))) ELSE XD3=3.D0/X F0=0.79788456D0-XD3*(0.00000077D0+XD3*(0.00552740D0+XD3* * (0.00009512D0-XD3*(0.00137237D0-XD3*(0.00072805D0 * -XD3*0.00014476D0))))) T0=X-0.78539816D0-XD3*(0.04166397D0+XD3*(0.00003954D0-XD3* * (0.00262573D0-XD3*(0.00054125D0+XD3*(0.00029333D0 * -XD3*0.00013558D0))))) BES0=F0/DSQRT(X)*DCOS(T0) ENDIF RETURN END c c-------------------------------------------------------------------------- c DOUBLE PRECISION FUNCTION BES1(X) C IMPLICIT REAL*8 (A-H,O-Z) C IF (DABS(X).LT.3.D0) THEN X32=(X/3.D0)**2 BES1XM1=0.5D0-X32*(0.56249985D0-X32*(0.21093573D0-X32* * (0.03954289D0-X32*(0.00443319D0-X32*(0.00031761D0 * -X32*0.00001109D0))))) BES1=BES1XM1*X ELSE XD3=3.D0/X F1=0.79788456D0+XD3*(0.00000156D0+XD3*(0.01659667D0+XD3* * (0.00017105D0-XD3*(0.00249511D0-XD3*(0.00113653D0 * -XD3*0.00020033D0))))) T1=X-2.35619449D0+XD3*(0.12499612D0+XD3*(0.0000565D0-XD3* * (0.00637879D0-XD3*(0.00074348D0+XD3*(0.00079824D0 * -XD3*0.00029166D0))))) BES1=F1/DSQRT(X)*DCOS(T1) ENDIF RETURN END C------------------------------------------------------------ C SUBROUTINE INTERCON(X,Y,Z,BX,BY,BZ) C C Calculates the potential interconnection field inside the magnetosphere, c corresponding to DELTA_X = 20Re and DELTA_Y = 10Re (NB#3, p.90, 6/6/1996). C The position (X,Y,Z) and field components BX,BY,BZ are given in the rotated c coordinate system, in which the Z-axis is always directed along the BzIMF c (i.e. rotated by the IMF clock angle Theta) C It is also assumed that the IMF Bt=1, so that the components should be c (i) multiplied by the actual Bt, and c (ii) transformed to standard GSM coords by rotating back around X axis c by the angle -Theta. c C Description of parameters: C C X,Y,Z - GSM POSITION C BX,BY,BZ - INTERCONNECTION FIELD COMPONENTS INSIDE THE MAGNETOSPHERE C OF A STANDARD SIZE (TO TAKE INTO ACCOUNT EFFECTS OF PRESSURE CHANGES, C APPLY THE SCALING TRANSFORMATION) C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C C The 9 linear parameters are amplitudes of the "cartesian" harmonics c The 6 nonlinear parameters are the scales Pi and Ri entering c the arguments of exponents, sines, and cosines in the 9 "Cartesian" c harmonics (3+3) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C IMPLICIT REAL * 8 (A - H, O - Z) C DIMENSION A(15),RP(3),RR(3),P(3),R(3) SAVE C DATA A/-8.411078731,5932254.951,-9073284.93,-11.68794634, * 6027598.824,-9218378.368,-6.508798398,-11824.42793,18015.66212, * 7.99754043,13.9669886,90.24475036,16.75728834,1015.645781, * 1553.493216/ C DATA M/0/ C IF (M.NE.0) GOTO 111 M=1 C P(1)=A(10) P(2)=A(11) P(3)=A(12) R(1)=A(13) R(2)=A(14) R(3)=A(15) C C DO 11 I=1,3 RP(I)=1.D0/P(I) 11 RR(I)=1.D0/R(I) C 111 CONTINUE C L=0 C BX=0. BY=0. BZ=0. C c "PERPENDICULAR" KIND OF SYMMETRY ONLY C DO 2 I=1,3 CYPI=DCOS(Y*RP(I)) SYPI=DSIN(Y*RP(I)) C DO 2 K=1,3 SZRK=DSIN(Z*RR(K)) CZRK=DCOS(Z*RR(K)) SQPR=DSQRT(RP(I)**2+RR(K)**2) EPR=DEXP(X*SQPR) C HX=-SQPR*EPR*CYPI*SZRK HY=RP(I)*EPR*SYPI*SZRK HZ=-RR(K)*EPR*CYPI*CZRK L=L+1 c BX=BX+A(L)*HX BY=BY+A(L)*HY BZ=BZ+A(L)*HZ 2 CONTINUE C RETURN END C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE TAILRC96(SPS,X,Y,Z,BXRC,BYRC,BZRC,BXT2,BYT2,BZT2, * BXT3,BYT3,BZT3) c c COMPUTES THE COMPONENTS OF THE FIELD OF THE MODEL RING CURRENT AND THREE c TAIL MODES WITH UNIT AMPLITUDES C (FOR THE RING CURRENT, IT MEANS THE DISTURBANCE OF Bz=-1nT AT ORIGIN, C AND FOR THE TAIL MODES IT MEANS MAXIMAL BX JUST ABOVE THE SHEET EQUAL 1 nT. C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION ARC(48),ATAIL2(48),ATAIL3(48) COMMON /WARP/ CPSS,SPSS,DPSRR,RPS,WARP,D,XS,ZS,DXSX,DXSY,DXSZ, * DZSX,DZSY,DZSZ,DZETAS,DDZETADX,DDZETADY,DDZETADZ,ZSWW C DATA ARC/-3.087699646,3.516259114,18.81380577,-13.95772338, * -5.497076303,0.1712890838,2.392629189,-2.728020808,-14.79349936, * 11.08738083,4.388174084,0.2492163197E-01,0.7030375685, *-.7966023165,-3.835041334,2.642228681,-0.2405352424,-0.7297705678, * -0.3680255045,0.1333685557,2.795140897,-1.078379954,0.8014028630, * 0.1245825565,0.6149982835,-0.2207267314,-4.424578723,1.730471572, * -1.716313926,-0.2306302941,-0.2450342688,0.8617173961E-01, * 1.54697858,-0.6569391113,-0.6537525353,0.2079417515,12.75434981, * 11.37659788,636.4346279,1.752483754,3.604231143,12.83078674, * 7.412066636,9.434625736,676.7557193,1.701162737,3.580307144, * 14.64298662/ C DATA ATAIL2/.8747515218,-.9116821411,2.209365387,-2.159059518, * -7.059828867,5.924671028,-1.916935691,1.996707344,-3.877101873, * 3.947666061,11.38715899,-8.343210833,1.194109867,-1.244316975, * 3.73895491,-4.406522465,-20.66884863,3.020952989,.2189908481, * -.09942543549,-.927225562,.1555224669,.6994137909,-.08111721003, * -.7565493881,.4686588792,4.266058082,-.3717470262,-3.920787807, * .02298569870,.7039506341,-.5498352719,-6.675140817,.8279283559, * -2.234773608,-1.622656137,5.187666221,6.802472048,39.13543412, * 2.784722096,6.979576616,25.71716760,4.495005873,8.068408272, * 93.47887103,4.158030104,9.313492566,57.18240483/ C DATA ATAIL3/-19091.95061,-3011.613928,20582.16203,4242.918430, * -2377.091102,-1504.820043,19884.04650,2725.150544,-21389.04845, * -3990.475093,2401.610097,1548.171792,-946.5493963,490.1528941, * 986.9156625,-489.3265930,-67.99278499,8.711175710,-45.15734260, * -10.76106500,210.7927312,11.41764141,-178.0262808,.7558830028, * 339.3806753,9.904695974,69.50583193,-118.0271581,22.85935896, * 45.91014857,-425.6607164,15.47250738,118.2988915,65.58594397, * -201.4478068,-14.57062940,19.69877970,20.30095680,86.45407420, * 22.50403727,23.41617329,48.48140573,24.61031329,123.5395974, * 223.5367692,39.50824342,65.83385762,266.2948657/ C DATA RH,DR,G,D0,DELTADY/9.,4.,10.,2.,10./ C C TO ECONOMIZE THE CODE, WE FIRST CALCULATE COMMON VARIABLES, WHICH ARE C THE SAME FOR ALL MODES, AND PUT THEM IN THE COMMON-BLOCK /WARP/ C DR2=DR*DR C11=DSQRT((1.D0+RH)**2+DR2) C12=DSQRT((1.D0-RH)**2+DR2) C1=C11-C12 SPSC1=SPS/C1 RPS=0.5*(C11+C12)*SPS ! THIS IS THE SHIFT OF OF THE SHEET WITH RESPECT C TO GSM EQ.PLANE FOR THE 3RD (ASYMPTOTIC) TAIL MODE C R=DSQRT(X*X+Y*Y+Z*Z) SQ1=DSQRT((R+RH)**2+DR2) SQ2=DSQRT((R-RH)**2+DR2) C=SQ1-SQ2 CS=(R+RH)/SQ1-(R-RH)/SQ2 SPSS=SPSC1/R*C CPSS=DSQRT(1.D0-SPSS**2) DPSRR=SPS/(R*R)*(CS*R-C)/DSQRT((R*C1)**2-(C*SPS)**2) C WFAC=Y/(Y**4+1.D4) ! WARPING W=WFAC*Y**3 WS=4.D4*Y*WFAC**2 WARP=G*SPS*W XS=X*CPSS-Z*SPSS ZSWW=Z*CPSS+X*SPSS ! "WW" MEANS "WITHOUT Y-Z WARPING" (IN X-Z ONLY) ZS=ZSWW +WARP DXSX=CPSS-X*ZSWW*DPSRR DXSY=-Y*ZSWW*DPSRR DXSZ=-SPSS-Z*ZSWW*DPSRR DZSX=SPSS+X*XS*DPSRR DZSY=XS*Y*DPSRR +G*SPS*WS ! THE LAST TERM IS FOR THE Y-Z WARP DZSZ=CPSS+XS*Z*DPSRR ! (TAIL MODES ONLY) D=D0+DELTADY*(Y/20.D0)**2 ! SHEET HALF-THICKNESS FOR THE TAIL MODES DDDY=DELTADY*Y*0.005D0 ! (THICKENS TO FLANKS, BUT NO VARIATION C ALONG X, IN CONTRAST TO RING CURRENT) C DZETAS=DSQRT(ZS**2+D**2) ! THIS IS THE SAME SIMPLE WAY TO SPREAD C OUT THE SHEET, AS THAT USED IN T89 DDZETADX=ZS*DZSX/DZETAS DDZETADY=(ZS*DZSY+D*DDDY)/DZETAS DDZETADZ=ZS*DZSZ/DZETAS C CALL SHLCAR3X3(ARC,X,Y,Z,SPS,WX,WY,WZ) CALL RINGCURR96(X,Y,Z,HX,HY,HZ) BXRC=WX+HX BYRC=WY+HY BZRC=WZ+HZ C CALL SHLCAR3X3(ATAIL2,X,Y,Z,SPS,WX,WY,WZ) CALL TAILDISK(X,Y,Z,HX,HY,HZ) BXT2=WX+HX BYT2=WY+HY BZT2=WZ+HZ C CALL SHLCAR3X3(ATAIL3,X,Y,Z,SPS,WX,WY,WZ) CALL TAIL87(X,Z,HX,HZ) BXT3=WX+HX BYT3=WY BZT3=WZ+HZ C RETURN END C c******************************************************************** C SUBROUTINE RINGCURR96(X,Y,Z,BX,BY,BZ) c c THIS SUBROUTINE COMPUTES THE COMPONENTS OF THE RING CURRENT FIELD, C SIMILAR TO THAT DESCRIBED BY TSYGANENKO AND PEREDO (1994). THE C DIFFERENCE IS THAT NOW WE USE SPACEWARPING, AS DESCRIBED IN THE C PAPER ON MODELING BIRKELAND CURRENTS (TSYGANENKO AND STERN, 1996), C INSTEAD OF SHEARING IT IN THE SPIRIT OF THE T89 TAIL MODEL. C C IN ADDITION, INSTEAD OF 7 TERMS FOR THE RING CURRENT MODEL, WE USE C NOW ONLY 2 TERMS; THIS SIMPLIFICATION ALSO GIVES RISE TO AN C EASTWARD RING CURRENT LOCATED EARTHWARD FROM THE MAIN ONE, C IN LINE WITH WHAT IS ACTUALLY OBSERVED C C FOR DETAILS, SEE NB #3, PAGES 70-73 C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION F(2),BETA(2) COMMON /WARP/ CPSS,SPSS,DPSRR, XNEXT(3),XS,ZSWARPED,DXSX,DXSY, * DXSZ,DZSX,DZSYWARPED,DZSZ,OTHER(4),ZS C ZS HERE IS WITHOUT Y-Z WARP C DATA D0,DELTADX,XD,XLDX /2.,0.,0.,4./ ! ACHTUNG !! THE RC IS NOW C COMPLETELY SYMMETRIC (DELTADX=0) C DATA F,BETA /569.895366D0,-1603.386993D0,2.722188D0,3.766875D0/ C C THE ORIGINAL VALUES OF F(I) WERE MULTIPLIED BY BETA(I) (TO REDUCE THE C NUMBER OF MULTIPLICATIONS BELOW) AND BY THE FACTOR -0.43, NORMALIZING C THE DISTURBANCE AT ORIGIN TO B=-1nT C DZSY=XS*Y*DPSRR ! NO WARPING IN THE Y-Z PLANE (ALONG X ONLY), AND C THIS IS WHY WE DO NOT USE DZSY FROM THE COMMON-BLOCK XXD=X-XD FDX=0.5D0*(1.D0+XXD/DSQRT(XXD**2+XLDX**2)) DDDX=DELTADX*0.5D0*XLDX**2/DSQRT(XXD**2+XLDX**2)**3 D=D0+DELTADX*FDX DZETAS=DSQRT(ZS**2+D**2) ! THIS IS THE SAME SIMPLE WAY TO SPREAD C OUT THE SHEET, AS THAT USED IN T89 RHOS=DSQRT(XS**2+Y**2) DDZETADX=(ZS*DZSX+D*DDDX)/DZETAS DDZETADY=ZS*DZSY/DZETAS DDZETADZ=ZS*DZSZ/DZETAS IF (RHOS.LT.1.D-5) THEN DRHOSDX=0.D0 DRHOSDY=DSIGN(1.D0,Y) DRHOSDZ=0.D0 ELSE DRHOSDX=XS*DXSX/RHOS DRHOSDY=(XS*DXSY+Y)/RHOS DRHOSDZ=XS*DXSZ/RHOS ENDIF C BX=0.D0 BY=0.D0 BZ=0.D0 C DO 1 I=1,2 C BI=BETA(I) C S1=DSQRT((DZETAS+BI)**2+(RHOS+BI)**2) S2=DSQRT((DZETAS+BI)**2+(RHOS-BI)**2) DS1DDZ=(DZETAS+BI)/S1 DS2DDZ=(DZETAS+BI)/S2 DS1DRHOS=(RHOS+BI)/S1 DS2DRHOS=(RHOS-BI)/S2 C DS1DX=DS1DDZ*DDZETADX+DS1DRHOS*DRHOSDX DS1DY=DS1DDZ*DDZETADY+DS1DRHOS*DRHOSDY DS1DZ=DS1DDZ*DDZETADZ+DS1DRHOS*DRHOSDZ C DS2DX=DS2DDZ*DDZETADX+DS2DRHOS*DRHOSDX DS2DY=DS2DDZ*DDZETADY+DS2DRHOS*DRHOSDY DS2DZ=DS2DDZ*DDZETADZ+DS2DRHOS*DRHOSDZ C S1TS2=S1*S2 S1PS2=S1+S2 S1PS2SQ=S1PS2**2 FAC1=DSQRT(S1PS2SQ-(2.D0*BI)**2) AS=FAC1/(S1TS2*S1PS2SQ) TERM1=1.D0/(S1TS2*S1PS2*FAC1) FAC2=AS/S1PS2SQ DASDS1=TERM1-FAC2/S1*(S2*S2+S1*(3.D0*S1+4.D0*S2)) DASDS2=TERM1-FAC2/S2*(S1*S1+S2*(3.D0*S2+4.D0*S1)) C DASDX=DASDS1*DS1DX+DASDS2*DS2DX DASDY=DASDS1*DS1DY+DASDS2*DS2DY DASDZ=DASDS1*DS1DZ+DASDS2*DS2DZ C BX=BX+F(I)*((2.D0*AS+Y*DASDY)*SPSS-XS*DASDZ * +AS*DPSRR*(Y**2*CPSS+Z*ZS)) BY=BY-F(I)*Y*(AS*DPSRR*XS+DASDZ*CPSS+DASDX*SPSS) 1 BZ=BZ+F(I)*((2.D0*AS+Y*DASDY)*CPSS+XS*DASDX * -AS*DPSRR*(X*ZS+Y**2*SPSS)) C RETURN END C C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C SUBROUTINE TAILDISK(X,Y,Z,BX,BY,BZ) C 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 C IN ADDITION, INSTEAD OF 8 TERMS FOR THE TAIL CURRENT MODEL, WE USE C NOW ONLY 4 TERMS C C FOR DETAILS, SEE NB #3, PAGES 74- C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION F(4),BETA(4) COMMON /WARP/ CPSS,SPSS,DPSRR,XNEXT(3),XS,ZS,DXSX,DXSY,DXSZ, * OTHER(3),DZETAS,DDZETADX,DDZETADY,DDZETADZ,ZSWW C DATA XSHIFT /4.5/ C DATA F,BETA * / -745796.7338D0,1176470.141D0,-444610.529D0,-57508.01028D0, * 7.9250000D0,8.0850000D0,8.4712500D0,27.89500D0/ c c here original F(I) are multiplied by BETA(I), to economize c calculations C RHOS=DSQRT((XS-XSHIFT)**2+Y**2) IF (RHOS.LT.1.D-5) THEN DRHOSDX=0.D0 DRHOSDY=DSIGN(1.D0,Y) DRHOSDZ=0.D0 ELSE DRHOSDX=(XS-XSHIFT)*DXSX/RHOS DRHOSDY=((XS-XSHIFT)*DXSY+Y)/RHOS DRHOSDZ=(XS-XSHIFT)*DXSZ/RHOS ENDIF C BX=0.D0 BY=0.D0 BZ=0.D0 C DO 1 I=1,4 C BI=BETA(I) C S1=DSQRT((DZETAS+BI)**2+(RHOS+BI)**2) S2=DSQRT((DZETAS+BI)**2+(RHOS-BI)**2) DS1DDZ=(DZETAS+BI)/S1 DS2DDZ=(DZETAS+BI)/S2 DS1DRHOS=(RHOS+BI)/S1 DS2DRHOS=(RHOS-BI)/S2 C DS1DX=DS1DDZ*DDZETADX+DS1DRHOS*DRHOSDX DS1DY=DS1DDZ*DDZETADY+DS1DRHOS*DRHOSDY DS1DZ=DS1DDZ*DDZETADZ+DS1DRHOS*DRHOSDZ C DS2DX=DS2DDZ*DDZETADX+DS2DRHOS*DRHOSDX DS2DY=DS2DDZ*DDZETADY+DS2DRHOS*DRHOSDY DS2DZ=DS2DDZ*DDZETADZ+DS2DRHOS*DRHOSDZ C S1TS2=S1*S2 S1PS2=S1+S2 S1PS2SQ=S1PS2**2 FAC1=DSQRT(S1PS2SQ-(2.D0*BI)**2) AS=FAC1/(S1TS2*S1PS2SQ) TERM1=1.D0/(S1TS2*S1PS2*FAC1) FAC2=AS/S1PS2SQ DASDS1=TERM1-FAC2/S1*(S2*S2+S1*(3.D0*S1+4.D0*S2)) DASDS2=TERM1-FAC2/S2*(S1*S1+S2*(3.D0*S2+4.D0*S1)) C DASDX=DASDS1*DS1DX+DASDS2*DS2DX DASDY=DASDS1*DS1DY+DASDS2*DS2DY DASDZ=DASDS1*DS1DZ+DASDS2*DS2DZ C BX=BX+F(I)*((2.D0*AS+Y*DASDY)*SPSS-(XS-XSHIFT)*DASDZ * +AS*DPSRR*(Y**2*CPSS+Z*ZSWW)) C BY=BY-F(I)*Y*(AS*DPSRR*XS+DASDZ*CPSS+DASDX*SPSS) 1 BZ=BZ+F(I)*((2.D0*AS+Y*DASDY)*CPSS+(XS-XSHIFT)*DASDX * -AS*DPSRR*(X*ZSWW+Y**2*SPSS)) RETURN END C------------------------------------------------------------------------- C SUBROUTINE TAIL87(X,Z,BX,BZ) IMPLICIT REAL*8 (A-H,O-Z) COMMON /WARP/ FIRST(3), RPS,WARP,D, OTHER(13) C C 'LONG' VERSION OF THE 1987 TAIL MAGNETIC FIELD MODEL C (N.A.TSYGANENKO, PLANET. SPACE SCI., V.35, P.1347, 1987) C C D IS THE Y-DEPENDENT SHEET HALF-THICKNESS (INCREASING TOWARDS FLANKS) C RPS IS THE TILT-DEPENDENT SHIFT OF THE SHEET IN THE Z-DIRECTION, C CORRESPONDING TO THE ASYMPTOTIC HINGING DISTANCE, DEFINED IN THE C MAIN SUBROUTINE (TAILRC96) FROM THE PARAMETERS RH AND DR OF THE C T96-TYPE MODULE, AND C WARP IS THE BENDING OF THE SHEET FLANKS IN THE Z-DIRECTION, DIRECTED C OPPOSITE TO RPS, AND INCREASING WITH DIPOLE TILT AND |Y| C DATA DD/3./ C DATA HPI,RT,XN,X1,X2,B0,B1,B2,XN21,XNR,ADLN * /1.5707963,40.,-10., * -1.261,-0.663,0.391734,5.89715,24.6833,76.37,-0.1071,0.13238005/ C !!! THESE ARE NEW VALUES OF X1, X2, B0, B1, B2, C CORRESPONDING TO TSCALE=1, INSTEAD OF TSCALE=0.6 C C THE ABOVE QUANTITIES WERE DEFINED AS FOLLOWS:------------------------ C HPI=PI/2 C RT=40. ! Z-POSITION OF UPPER AND LOWER ADDITIONAL SHEETS C XN=-10. ! INNER EDGE POSITION C C TSCALE=1 ! SCALING FACTOR, DEFINING THE RATE OF INCREASE OF THE C CURRENT DENSITY TAILWARDS C c ATTENTION ! NOW I HAVE CHANGED TSCALE TO: TSCALE=1.0, INSTEAD OF 0.6 c OF THE PREVIOUS VERSION c C B0=0.391734 C B1=5.89715 *TSCALE C B2=24.6833 *TSCALE**2 C C HERE ORIGINAL VALUES OF THE MODE AMPLITUDES (P.77, NB#3) WERE NORMALIZED C SO THAT ASYMPTOTIC BX=1 AT X=-200RE C C X1=(4.589 -5.85) *TSCALE -(TSCALE-1.)*XN ! NONLINEAR PARAMETERS OF THE C CURRENT FUNCTION C X2=(5.187 -5.85) *TSCALE -(TSCALE-1.)*XN c c C XN21=(XN-X1)**2 C XNR=1./(XN-X2) C ADLN=-DLOG(XNR**2*XN21) C C--------------------------------------------------------------- C ZS=Z -RPS +WARP ZP=Z-RT ZM=Z+RT C XNX=XN-X XNX2=XNX**2 XC1=X-X1 XC2=X-X2 XC22=XC2**2 XR2=XC2*XNR XC12=XC1**2 D2=DD**2 ! SQUARE OF THE TOTAL HALFTHICKNESS (DD=3Re for this mode) B20=ZS**2+D2 B2P=ZP**2+D2 B2M=ZM**2+D2 B=DSQRT(B20) BP=DSQRT(B2P) BM=DSQRT(B2M) XA1=XC12+B20 XAP1=XC12+B2P XAM1=XC12+B2M XA2=1./(XC22+B20) XAP2=1./(XC22+B2P) XAM2=1./(XC22+B2M) XNA=XNX2+B20 XNAP=XNX2+B2P XNAM=XNX2+B2M F=B20-XC22 FP=B2P-XC22 FM=B2M-XC22 XLN1=DLOG(XN21/XNA) XLNP1=DLOG(XN21/XNAP) XLNM1=DLOG(XN21/XNAM) XLN2=XLN1+ADLN XLNP2=XLNP1+ADLN XLNM2=XLNM1+ADLN ALN=0.25*(XLNP1+XLNM1-2.*XLN1) S0=(DATAN(XNX/B)+HPI)/B S0P=(DATAN(XNX/BP)+HPI)/BP S0M=(DATAN(XNX/BM)+HPI)/BM S1=(XLN1*.5+XC1*S0)/XA1 S1P=(XLNP1*.5+XC1*S0P)/XAP1 S1M=(XLNM1*.5+XC1*S0M)/XAM1 S2=(XC2*XA2*XLN2-XNR-F*XA2*S0)*XA2 S2P=(XC2*XAP2*XLNP2-XNR-FP*XAP2*S0P)*XAP2 S2M=(XC2*XAM2*XLNM2-XNR-FM*XAM2*S0M)*XAM2 G1=(B20*S0-0.5*XC1*XLN1)/XA1 G1P=(B2P*S0P-0.5*XC1*XLNP1)/XAP1 G1M=(B2M*S0M-0.5*XC1*XLNM1)/XAM1 G2=((0.5*F*XLN2+2.*S0*B20*XC2)*XA2+XR2)*XA2 G2P=((0.5*FP*XLNP2+2.*S0P*B2P*XC2)*XAP2+XR2)*XAP2 G2M=((0.5*FM*XLNM2+2.*S0M*B2M*XC2)*XAM2+XR2)*XAM2 BX=B0*(ZS*S0-0.5*(ZP*S0P+ZM*S0M)) * +B1*(ZS*S1-0.5*(ZP*S1P+ZM*S1M))+B2*(ZS*S2-0.5*(ZP*S2P+ZM*S2M)) BZ=B0*ALN+B1*(G1-0.5*(G1P+G1M))+B2*(G2-0.5*(G2P+G2M)) C C CALCULATION OF THE MAGNETOTAIL CURRENT CONTRIBUTION IS FINISHED C RETURN END C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C C THIS CODE RETURNS THE SHIELDING FIELD REPRESENTED BY 2x3x3=18 "CARTESIAN" C HARMONICS C SUBROUTINE SHLCAR3X3(A,X,Y,Z,SPS,HX,HY,HZ) C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C The 36 coefficients enter in pairs in the amplitudes of the "cartesian" c harmonics (A(1)-A(36). c The 12 nonlinear parameters (A(37)-A(48) 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 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C IMPLICIT REAL * 8 (A - H, O - Z) C DIMENSION A(48) C CPS=DSQRT(1.D0-SPS**2) S3PS=4.D0*CPS**2-1.D0 ! THIS IS SIN(3*PS)/SIN(PS) C HX=0.D0 HY=0.D0 HZ=0.D0 L=0 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(36+I) Q=A(42+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(39+K) S=A(45+K) SZRK=DSIN(Z/R) CZSK=DCOS(Z/S) CZRK=DCOS(Z/R) SZSK=DSIN(Z/S) SQPR=DSQRT(1.D0/P**2+1.D0/R**2) SQQS=DSQRT(1.D0/Q**2+1.D0/S**2) EPR=DEXP(X*SQPR) EQS=DEXP(X*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 C L=L+1 IF (M.EQ.1) THEN IF (N.EQ.1) THEN DX=-SQPR*EPR*CYPI*SZRK DY=EPR/P*SYPI*SZRK DZ=-EPR/R*CYPI*CZRK HX=HX+A(L)*DX HY=HY+A(L)*DY HZ=HZ+A(L)*DZ ELSE DX=DX*CPS DY=DY*CPS DZ=DZ*CPS HX=HX+A(L)*DX HY=HY+A(L)*DY HZ=HZ+A(L)*DZ ENDIF ELSE IF (N.EQ.1) THEN DX=-SPS*SQQS*EQS*CYQI*CZSK DY=SPS*EQS/Q*SYQI*CZSK DZ=SPS*EQS/S*CYQI*SZSK HX=HX+A(L)*DX HY=HY+A(L)*DY HZ=HZ+A(L)*DZ ELSE DX=DX*S3PS DY=DY*S3PS DZ=DZ*S3PS HX=HX+A(L)*DX HY=HY+A(L)*DY HZ=HZ+A(L)*DZ ENDIF ENDIF c 4 CONTINUE 3 CONTINUE 2 CONTINUE 1 CONTINUE C RETURN END C C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ C SUBROUTINE BIRK1TOT_02(PS,X,Y,Z,BX,BY,BZ) C C THIS IS THE SECOND VERSION OF THE ANALYTICAL MODEL OF THE REGION 1 FIELD C BASED ON A SEPARATE REPRESENTATION OF THE POTENTIAL FIELD IN THE INNER AND C OUTER SPACE, MAPPED BY MEANS OF A SPHERO-DIPOLAR COORDINATE SYSTEM (NB #3, C P.91). THE DIFFERENCE FROM THE FIRST ONE IS THAT INSTEAD OF OCTAGONAL C CURRENT LOOPS, CIRCULAR ONES ARE USED IN THIS VERSION FOR APPROXIMATING THE C FIELD IN THE OUTER REGION, WHICH IS FASTER. C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION D1(3,26),D2(3,79),XI(4),C1(26),C2(79) COMMON /COORD11/ XX1(12),YY1(12) COMMON /RHDR/ RH,DR COMMON /LOOPDIP1/ TILT,XCENTRE(2),RADIUS(2), DIPX,DIPY C COMMON /COORD21/ XX2(14),YY2(14),ZZ2(14) COMMON /DX1/ DX,SCALEIN,SCALEOUT C DATA C1/-0.911582E-03,-0.376654E-02,-0.727423E-02,-0.270084E-02, * -0.123899E-02,-0.154387E-02,-0.340040E-02,-0.191858E-01, * -0.518979E-01,0.635061E-01,0.440680,-0.396570,0.561238E-02, * 0.160938E-02,-0.451229E-02,-0.251810E-02,-0.151599E-02, * -0.133665E-02,-0.962089E-03,-0.272085E-01,-0.524319E-01, * 0.717024E-01,0.523439,-0.405015,-89.5587,23.2806/ C DATA C2/6.04133,.305415,.606066E-02,.128379E-03,-.179406E-04, * 1.41714,-27.2586,-4.28833,-1.30675,35.5607,8.95792,.961617E-03, * -.801477E-03,-.782795E-03,-1.65242,-16.5242,-5.33798,.424878E-03, * .331787E-03,-.704305E-03,.844342E-03,.953682E-04,.886271E-03, * 25.1120,20.9299,5.14569,-44.1670,-51.0672,-1.87725,20.2998, * 48.7505,-2.97415,3.35184,-54.2921,-.838712,-10.5123,70.7594, * -4.94104,.106166E-03,.465791E-03,-.193719E-03,10.8439,-29.7968, * 8.08068,.463507E-03,-.224475E-04,.177035E-03,-.317581E-03, * -.264487E-03,.102075E-03,7.71390,10.1915,-4.99797,-23.1114, *-29.2043,12.2928,10.9542,33.6671,-9.3851,.174615E-03,-.789777E-06, * .686047E-03,.460104E-04,-.345216E-02,.221871E-02,.110078E-01, * -.661373E-02,.249201E-02,.343978E-01,-.193145E-05,.493963E-05, * -.535748E-04,.191833E-04,-.100496E-03,-.210103E-03,-.232195E-02, * .315335E-02,-.134320E-01,-.263222E-01/ c DATA TILT,XCENTRE,RADIUS,DIPX,DIPY /1.00891,2.28397,-5.60831, * 1.86106,7.83281,1.12541,0.945719/ DATA DX,SCALEIN,SCALEOUT /-0.16D0,0.08D0,0.4D0/ DATA XX1/-11.D0,2*-7.D0,2*-3.D0,3*1.D0,2*5.D0,2*9.D0/ DATA YY1/2.D0,0.D0,4.D0,2.D0,6.D0,0.D0,4.D0,8.D0,2.D0,6.D0,0.D0, * 4.D0/ DATA XX2/-10.D0,-7.D0,2*-4.D0,0.D0,2*4.D0,7.D0,10.D0,5*0.D0/ DATA YY2/3.D0,6.D0,3.D0,9.D0,6.D0,3.D0,9.D0,6.D0,3.D0,5*0.D0/ DATA ZZ2/2*20.D0,4.D0,20.D0,2*4.D0,3*20.D0,2.D0,3.D0,4.5D0, * 7.D0,10.D0/ C DATA RH,DR /9.D0,4.D0/ ! RH IS THE "HINGING DISTANCE" AND DR IS THE C TRANSITION SCALE LENGTH, DEFINING THE C CURVATURE OF THE WARPING (SEE P.89, NB #2) C DATA XLTDAY,XLTNGHT /78.D0,70.D0/ ! THESE ARE LATITUDES OF THE R-1 OVAL C AT NOON AND AT MIDNIGHT DATA DTET0 /0.034906/ ! THIS IS THE LATITUDINAL HALF-THICKNESS OF THE C R-1 OVAL (THE INTERPOLATION REGION BETWEEN C THE HIGH-LAT. AND THE PLASMA SHEET) C TNOONN=(90.D0-XLTDAY)*0.01745329D0 TNOONS=3.141592654D0-TNOONN ! HERE WE ASSUME THAT THE POSITIONS OF C THE NORTHERN AND SOUTHERN R-1 OVALS C ARE SYMMETRIC IN THE SM-COORDINATES DTETDN=(XLTDAY-XLTNGHT)*0.01745329D0 DR2=DR**2 C SPS=DSIN(PS) R2=X**2+Y**2+Z**2 R=DSQRT(R2) R3=R*R2 C RMRH=R-RH RPRH=R+RH SQM=DSQRT(RMRH**2+DR2) SQP=DSQRT(RPRH**2+DR2) C=SQP-SQM Q=DSQRT((RH+1.D0)**2+DR2)-DSQRT((RH-1.D0)**2+DR2) SPSAS=SPS/R*C/Q CPSAS=DSQRT(1.D0-SPSAS**2) XAS = X*CPSAS-Z*SPSAS ZAS = X*SPSAS+Z*CPSAS IF (XAS.NE.0.D0.OR.Y.NE.0.D0) THEN PAS = DATAN2(Y,XAS) ELSE PAS=0.D0 ENDIF C TAS=DATAN2(DSQRT(XAS**2+Y**2),ZAS) STAS=DSIN(TAS) F=STAS/(STAS**6*(1.D0-R3)+R3)**0.1666666667D0 C TET0=DASIN(F) IF (TAS.GT.1.5707963D0) TET0=3.141592654D0-TET0 DTET=DTETDN*DSIN(PAS*0.5D0)**2 TETR1N=TNOONN+DTET TETR1S=TNOONS-DTET C C NOW LET S DEFINE WHICH OF THE FOUR REGIONS (HIGH-LAT., NORTHERN PSBL, C PLASMA SHEET, SOUTHERN PSBL) DOES THE POINT (X,Y,Z) BELONG TO: C IF (TET0.LT.TETR1N-DTET0.OR.TET0.GT.TETR1S+DTET0) LOC=1 ! HIGH-LAT. IF (TET0.GT.TETR1N+DTET0.AND.TET0.LT.TETR1S-DTET0) LOC=2 ! PL.SHEET IF (TET0.GE.TETR1N-DTET0.AND.TET0.LE.TETR1N+DTET0) LOC=3 ! NORTH PSBL IF (TET0.GE.TETR1S-DTET0.AND.TET0.LE.TETR1S+DTET0) LOC=4 ! SOUTH PSBL C IF (LOC.EQ.1) THEN ! IN THE HIGH-LAT. REGION USE THE SUBROUTINE DIPOCT C C print *, ' LOC=1 (HIGH-LAT)' ! (test printout; disabled now) XI(1)=X XI(2)=Y XI(3)=Z XI(4)=PS CALL DIPLOOP1(XI,D1) BX=0.D0 BY=0.D0 BZ=0.D0 DO 1 I=1,26 BX=BX+C1(I)*D1(1,I) BY=BY+C1(I)*D1(2,I) 1 BZ=BZ+C1(I)*D1(3,I) ENDIF ! END OF THE CASE 1 C IF (LOC.EQ.2) THEN C print *, ' LOC=2 (PLASMA SHEET)' ! (test printout; disabled now) C XI(1)=X XI(2)=Y XI(3)=Z XI(4)=PS CALL CONDIP1(XI,D2) BX=0.D0 BY=0.D0 BZ=0.D0 DO 2 I=1,79 BX=BX+C2(I)*D2(1,I) BY=BY+C2(I)*D2(2,I) 2 BZ=BZ+C2(I)*D2(3,I) ENDIF ! END OF THE CASE 2 C IF (LOC.EQ.3) THEN C print *, ' LOC=3 (north PSBL)' ! (test printout; disabled now) C T01=TETR1N-DTET0 T02=TETR1N+DTET0 SQR=DSQRT(R) ST01AS=SQR/(R3+1.D0/DSIN(T01)**6-1.D0)**0.1666666667 ST02AS=SQR/(R3+1.D0/DSIN(T02)**6-1.D0)**0.1666666667 CT01AS=DSQRT(1.D0-ST01AS**2) CT02AS=DSQRT(1.D0-ST02AS**2) XAS1=R*ST01AS*DCOS(PAS) Y1= R*ST01AS*DSIN(PAS) ZAS1=R*CT01AS X1=XAS1*CPSAS+ZAS1*SPSAS Z1=-XAS1*SPSAS+ZAS1*CPSAS ! X1,Y1,Z1 ARE COORDS OF THE NORTHERN c BOUNDARY POINT XI(1)=X1 XI(2)=Y1 XI(3)=Z1 XI(4)=PS CALL DIPLOOP1(XI,D1) BX1=0.D0 BY1=0.D0 BZ1=0.D0 DO 11 I=1,26 BX1=BX1+C1(I)*D1(1,I) ! BX1,BY1,BZ1 ARE FIELD COMPONENTS BY1=BY1+C1(I)*D1(2,I) ! IN THE NORTHERN BOUNDARY POINT 11 BZ1=BZ1+C1(I)*D1(3,I) ! C XAS2=R*ST02AS*DCOS(PAS) Y2= R*ST02AS*DSIN(PAS) ZAS2=R*CT02AS X2=XAS2*CPSAS+ZAS2*SPSAS Z2=-XAS2*SPSAS+ZAS2*CPSAS ! X2,Y2,Z2 ARE COORDS OF THE SOUTHERN C BOUNDARY POINT XI(1)=X2 XI(2)=Y2 XI(3)=Z2 XI(4)=PS CALL CONDIP1(XI,D2) BX2=0.D0 BY2=0.D0 BZ2=0.D0 DO 12 I=1,79 BX2=BX2+C2(I)*D2(1,I)! BX2,BY2,BZ2 ARE FIELD COMPONENTS BY2=BY2+C2(I)*D2(2,I) ! IN THE SOUTHERN BOUNDARY POINT 12 BZ2=BZ2+C2(I)*D2(3,I) C C NOW INTERPOLATE: C SS=DSQRT((X2-X1)**2+(Y2-Y1)**2+(Z2-Z1)**2) DS=DSQRT((X-X1)**2+(Y-Y1)**2+(Z-Z1)**2) FRAC=DS/SS BX=BX1*(1.D0-FRAC)+BX2*FRAC BY=BY1*(1.D0-FRAC)+BY2*FRAC BZ=BZ1*(1.D0-FRAC)+BZ2*FRAC C ENDIF ! END OF THE CASE 3 C IF (LOC.EQ.4) THEN C print *, ' LOC=4 (south PSBL)' ! (test printout; disabled now) C T01=TETR1S-DTET0 T02=TETR1S+DTET0 SQR=DSQRT(R) ST01AS=SQR/(R3+1.D0/DSIN(T01)**6-1.D0)**0.1666666667 ST02AS=SQR/(R3+1.D0/DSIN(T02)**6-1.D0)**0.1666666667 CT01AS=-DSQRT(1.D0-ST01AS**2) CT02AS=-DSQRT(1.D0-ST02AS**2) XAS1=R*ST01AS*DCOS(PAS) Y1= R*ST01AS*DSIN(PAS) ZAS1=R*CT01AS X1=XAS1*CPSAS+ZAS1*SPSAS Z1=-XAS1*SPSAS+ZAS1*CPSAS ! X1,Y1,Z1 ARE COORDS OF THE NORTHERN C BOUNDARY POINT XI(1)=X1 XI(2)=Y1 XI(3)=Z1 XI(4)=PS CALL CONDIP1(XI,D2) BX1=0.D0 BY1=0.D0 BZ1=0.D0 DO 21 I=1,79 BX1=BX1+C2(I)*D2(1,I) ! BX1,BY1,BZ1 ARE FIELD COMPONENTS BY1=BY1+C2(I)*D2(2,I) ! IN THE NORTHERN BOUNDARY POINT 21 BZ1=BZ1+C2(I)*D2(3,I) ! C XAS2=R*ST02AS*DCOS(PAS) Y2= R*ST02AS*DSIN(PAS) ZAS2=R*CT02AS X2=XAS2*CPSAS+ZAS2*SPSAS Z2=-XAS2*SPSAS+ZAS2*CPSAS ! X2,Y2,Z2 ARE COORDS OF THE SOUTHERN C BOUNDARY POINT XI(1)=X2 XI(2)=Y2 XI(3)=Z2 XI(4)=PS CALL DIPLOOP1(XI,D1) BX2=0.D0 BY2=0.D0 BZ2=0.D0 DO 22 I=1,26 BX2=BX2+C1(I)*D1(1,I) ! BX2,BY2,BZ2 ARE FIELD COMPONENTS BY2=BY2+C1(I)*D1(2,I) ! IN THE SOUTHERN BOUNDARY POINT 22 BZ2=BZ2+C1(I)*D1(3,I) C C NOW INTERPOLATE: C SS=DSQRT((X2-X1)**2+(Y2-Y1)**2+(Z2-Z1)**2) DS=DSQRT((X-X1)**2+(Y-Y1)**2+(Z-Z1)**2) FRAC=DS/SS BX=BX1*(1.D0-FRAC)+BX2*FRAC BY=BY1*(1.D0-FRAC)+BY2*FRAC BZ=BZ1*(1.D0-FRAC)+BZ2*FRAC C ENDIF ! END OF THE CASE 4 C C NOW, LET US ADD THE SHIELDING FIELD C CALL BIRK1SHLD(PS,X,Y,Z,BSX,BSY,BSZ) BX=BX+BSX BY=BY+BSY BZ=BZ+BSZ RETURN END C C------------------------------------------------------------------------------ C C SUBROUTINE DIPLOOP1(XI,D) C 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 XI - input vector containing independent variables; C D - output double precision vector containing C calculated values for derivatives of dependent C variables with respect to LINEAR model parameters; C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C c The 26 coefficients are moments (Z- and X-components) of 12 dipoles placed C inside the R1-shell, PLUS amplitudes of two octagonal double loops. C The dipoles with nonzero Yi appear in pairs with equal moments. c (see the notebook #2, pp.102-103, for details) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c IMPLICIT REAL * 8 (A - H, O - Z) C COMMON /COORD11/ XX(12),YY(12) COMMON /LOOPDIP1/ TILT,XCENTRE(2),RADIUS(2), DIPX,DIPY COMMON /RHDR/RH,DR DIMENSION XI(4),D(3,26) C X = XI(1) Y = XI(2) Z = XI(3) PS= XI(4) SPS=DSIN(PS) C DO 1 I=1,12 R2=(XX(I)*DIPX)**2+(YY(I)*DIPY)**2 R=DSQRT(R2) RMRH=R-RH RPRH=R+RH DR2=DR**2 SQM=DSQRT(RMRH**2+DR2) SQP=DSQRT(RPRH**2+DR2) C=SQP-SQM Q=DSQRT((RH+1.D0)**2+DR2)-DSQRT((RH-1.D0)**2+DR2) SPSAS=SPS/R*C/Q CPSAS=DSQRT(1.D0-SPSAS**2) XD= (XX(I)*DIPX)*CPSAS YD= (YY(I)*DIPY) ZD=-(XX(I)*DIPX)*SPSAS CALL DIPXYZ(X-XD,Y-YD,Z-ZD,BX1X,BY1X,BZ1X,BX1Y,BY1Y,BZ1Y, * BX1Z,BY1Z,BZ1Z) IF (DABS(YD).GT.1.D-10) THEN CALL DIPXYZ(X-XD,Y+YD,Z-ZD,BX2X,BY2X,BZ2X,BX2Y,BY2Y,BZ2Y, * BX2Z,BY2Z,BZ2Z) ELSE BX2X=0.D0 BY2X=0.D0 BZ2X=0.D0 C BX2Z=0.D0 BY2Z=0.D0 BZ2Z=0.D0 ENDIF C D(1,I)=BX1Z+BX2Z D(2,I)=BY1Z+BY2Z D(3,I)=BZ1Z+BZ2Z D(1,I+12)=(BX1X+BX2X)*SPS D(2,I+12)=(BY1X+BY2X)*SPS D(3,I+12)=(BZ1X+BZ2X)*SPS 1 CONTINUE c R2=(XCENTRE(1)+RADIUS(1))**2 R=DSQRT(R2) RMRH=R-RH RPRH=R+RH DR2=DR**2 SQM=DSQRT(RMRH**2+DR2) SQP=DSQRT(RPRH**2+DR2) C=SQP-SQM Q=DSQRT((RH+1.D0)**2+DR2)-DSQRT((RH-1.D0)**2+DR2) SPSAS=SPS/R*C/Q CPSAS=DSQRT(1.D0-SPSAS**2) XOCT1= X*CPSAS-Z*SPSAS YOCT1= Y ZOCT1= X*SPSAS+Z*CPSAS C CALL CROSSLP(XOCT1,YOCT1,ZOCT1,BXOCT1,BYOCT1,BZOCT1,XCENTRE(1), * RADIUS(1),TILT) D(1,25)=BXOCT1*CPSAS+BZOCT1*SPSAS D(2,25)=BYOCT1 D(3,25)=-BXOCT1*SPSAS+BZOCT1*CPSAS C R2=(RADIUS(2)-XCENTRE(2))**2 R=DSQRT(R2) RMRH=R-RH RPRH=R+RH DR2=DR**2 SQM=DSQRT(RMRH**2+DR2) SQP=DSQRT(RPRH**2+DR2) C=SQP-SQM Q=DSQRT((RH+1.D0)**2+DR2)-DSQRT((RH-1.D0)**2+DR2) SPSAS=SPS/R*C/Q CPSAS=DSQRT(1.D0-SPSAS**2) XOCT2= X*CPSAS-Z*SPSAS -XCENTRE(2) YOCT2= Y ZOCT2= X*SPSAS+Z*CPSAS CALL CIRCLE(XOCT2,YOCT2,ZOCT2,RADIUS(2),BX,BY,BZ) D(1,26) = BX*CPSAS+BZ*SPSAS D(2,26) = BY D(3,26) = -BX*SPSAS+BZ*CPSAS C RETURN END c------------------------------------------------------------------------- C SUBROUTINE CIRCLE(X,Y,Z,RL,BX,BY,BZ) C C RETURNS COMPONENTS OF THE FIELD FROM A CIRCULAR CURRENT LOOP OF RADIUS RL C USES THE SECOND (MORE ACCURATE) APPROXIMATION GIVEN IN ABRAMOWITZ AND STEGUN IMPLICIT REAL*8 (A-H,O-Z) REAL*8 K DATA PI/3.141592654D0/ C RHO2=X*X+Y*Y RHO=DSQRT(RHO2) R22=Z*Z+(RHO+RL)**2 R2=DSQRT(R22) R12=R22-4.D0*RHO*RL R32=0.5D0*(R12+R22) XK2=1.D0-R12/R22 XK2S=1.D0-XK2 DL=DLOG(1.D0/XK2S) K=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)))) E=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))) IF (RHO.GT.1.D-6) THEN BRHO=Z/(RHO2*R2)*(R32/R12*E-K) ! THIS IS NOT EXACTLY THE B-RHO COM- ELSE ! PONENT - NOTE THE ADDITIONAL BRHO=PI*RL/R2*(RL-RHO)/R12*Z/(R32-RHO2) ! DIVISION BY RHO ENDIF BX=BRHO*X BY=BRHO*Y BZ=(K-E*(R32-2.D0*RL*RL)/R12)/R2 RETURN END C------------------------------------------------------------- C SUBROUTINE CROSSLP(X,Y,Z,BX,BY,BZ,XC,RL,AL) C c RETURNS FIELD COMPONENTS OF A PAIR OF LOOPS WITH A COMMON CENTER AND C DIAMETER, COINCIDING WITH THE X AXIS. THE LOOPS ARE INCLINED TO THE C EQUATORIAL PLANE BY THE ANGLE AL (RADIANS) AND SHIFTED IN THE POSITIVE C X-DIRECTION BY THE DISTANCE XC. c IMPLICIT REAL*8 (A-H,O-Z) C CAL=DCOS(AL) SAL=DSIN(AL) C Y1=Y*CAL-Z*SAL Z1=Y*SAL+Z*CAL Y2=Y*CAL+Z*SAL Z2=-Y*SAL+Z*CAL CALL CIRCLE(X-XC,Y1,Z1,RL,BX1,BY1,BZ1) CALL CIRCLE(X-XC,Y2,Z2,RL,BX2,BY2,BZ2) BX=BX1+BX2 BY= (BY1+BY2)*CAL+(BZ1-BZ2)*SAL BZ=-(BY1-BY2)*SAL+(BZ1+BZ2)*CAL C RETURN END C******************************************************************* SUBROUTINE DIPXYZ(X,Y,Z,BXX,BYX,BZX,BXY,BYY,BZY,BXZ,BYZ,BZZ) C C RETURNS THE FIELD COMPONENTS PRODUCED BY THREE DIPOLES, EACH C HAVING M=Me AND ORIENTED PARALLEL TO X,Y, and Z AXIS, RESP. C IMPLICIT REAL*8 (A-H,O-Z) C X2=X**2 Y2=Y**2 Z2=Z**2 R2=X2+Y2+Z2 XMR5=30574.D0/(R2*R2*DSQRT(R2)) XMR53=3.D0*XMR5 BXX=XMR5*(3.D0*X2-R2) BYX=XMR53*X*Y BZX=XMR53*X*Z C BXY=BYX BYY=XMR5*(3.D0*Y2-R2) BZY=XMR53*Y*Z C BXZ=BZX BYZ=BZY BZZ=XMR5*(3.D0*Z2-R2) C RETURN END C C------------------------------------------------------------------------------ SUBROUTINE CONDIP1(XI,D) C C Calculates dependent model variables and their derivatives for given C independent variables and model parameters. Specifies model functions with C free parameters which must be determined by means of least squares fits C (RMS minimization procedure). C C Description of parameters: C C XI - input vector containing independent variables; C D - output double precision vector containing C calculated values for derivatives of dependent C variables with respect to LINEAR model parameters; C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C c The 79 coefficients are (1) 5 amplitudes of the conical harmonics, plus c (2) (9x3+5x2)x2=74 components of the dipole moments c (see the notebook #2, pp.113-..., for details) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c IMPLICIT REAL * 8 (A - H, O - Z) C COMMON /DX1/ DX,SCALEIN,SCALEOUT COMMON /COORD21/ XX(14),YY(14),ZZ(14) c DIMENSION XI(4),D(3,79),CF(5),SF(5) C X = XI(1) Y = XI(2) Z = XI(3) PS= XI(4) SPS=DSIN(PS) CPS=DCOS(PS) C XSM=X*CPS-Z*SPS - DX ZSM=Z*CPS+X*SPS RO2=XSM**2+Y**2 RO=SQRT(RO2) C CF(1)=XSM/RO SF(1)=Y/RO C CF(2)=CF(1)**2-SF(1)**2 SF(2)=2.*SF(1)*CF(1) CF(3)=CF(2)*CF(1)-SF(2)*SF(1) SF(3)=SF(2)*CF(1)+CF(2)*SF(1) CF(4)=CF(3)*CF(1)-SF(3)*SF(1) SF(4)=SF(3)*CF(1)+CF(3)*SF(1) CF(5)=CF(4)*CF(1)-SF(4)*SF(1) SF(5)=SF(4)*CF(1)+CF(4)*SF(1) C R2=RO2+ZSM**2 R=DSQRT(R2) C=ZSM/R S=RO/R CH=DSQRT(0.5D0*(1.D0+C)) SH=DSQRT(0.5D0*(1.D0-C)) TNH=SH/CH CNH=1.D0/TNH C DO 1 M=1,5 BT=M*CF(M)/(R*S)*(TNH**M+CNH**M) BF=-0.5D0*M*SF(M)/R*(TNH**(M-1)/CH**2-CNH**(M-1)/SH**2) BXSM=BT*C*CF(1)-BF*SF(1) BY=BT*C*SF(1)+BF*CF(1) BZSM=-BT*S C D(1,M)=BXSM*CPS+BZSM*SPS D(2,M)=BY 1 D(3,M)=-BXSM*SPS+BZSM*CPS C XSM = X*CPS-Z*SPS ZSM = Z*CPS+X*SPS C DO 2 I=1,9 C IF (I.EQ.3.OR.I.EQ.5.OR.I.EQ.6) THEN XD = XX(I)*SCALEIN YD = YY(I)*SCALEIN ELSE XD = XX(I)*SCALEOUT YD = YY(I)*SCALEOUT ENDIF C ZD = ZZ(I) C CALL DIPXYZ(XSM-XD,Y-YD,ZSM-ZD,BX1X,BY1X,BZ1X,BX1Y,BY1Y,BZ1Y, * BX1Z,BY1Z,BZ1Z) CALL DIPXYZ(XSM-XD,Y+YD,ZSM-ZD,BX2X,BY2X,BZ2X,BX2Y,BY2Y,BZ2Y, * BX2Z,BY2Z,BZ2Z) CALL DIPXYZ(XSM-XD,Y-YD,ZSM+ZD,BX3X,BY3X,BZ3X,BX3Y,BY3Y,BZ3Y, * BX3Z,BY3Z,BZ3Z) CALL DIPXYZ(XSM-XD,Y+YD,ZSM+ZD,BX4X,BY4X,BZ4X,BX4Y,BY4Y,BZ4Y, * BX4Z,BY4Z,BZ4Z) C IX=I*3+3 IY=IX+1 IZ=IY+1 C D(1,IX)=(BX1X+BX2X-BX3X-BX4X)*CPS+(BZ1X+BZ2X-BZ3X-BZ4X)*SPS D(2,IX)= BY1X+BY2X-BY3X-BY4X D(3,IX)=(BZ1X+BZ2X-BZ3X-BZ4X)*CPS-(BX1X+BX2X-BX3X-BX4X)*SPS C D(1,IY)=(BX1Y-BX2Y-BX3Y+BX4Y)*CPS+(BZ1Y-BZ2Y-BZ3Y+BZ4Y)*SPS D(2,IY)= BY1Y-BY2Y-BY3Y+BY4Y D(3,IY)=(BZ1Y-BZ2Y-BZ3Y+BZ4Y)*CPS-(BX1Y-BX2Y-BX3Y+BX4Y)*SPS C D(1,IZ)=(BX1Z+BX2Z+BX3Z+BX4Z)*CPS+(BZ1Z+BZ2Z+BZ3Z+BZ4Z)*SPS D(2,IZ)= BY1Z+BY2Z+BY3Z+BY4Z D(3,IZ)=(BZ1Z+BZ2Z+BZ3Z+BZ4Z)*CPS-(BX1Z+BX2Z+BX3Z+BX4Z)*SPS C IX=IX+27 IY=IY+27 IZ=IZ+27 C D(1,IX)=SPS*((BX1X+BX2X+BX3X+BX4X)*CPS+(BZ1X+BZ2X+BZ3X+BZ4X)*SPS) D(2,IX)=SPS*(BY1X+BY2X+BY3X+BY4X) D(3,IX)=SPS*((BZ1X+BZ2X+BZ3X+BZ4X)*CPS-(BX1X+BX2X+BX3X+BX4X)*SPS) C D(1,IY)=SPS*((BX1Y-BX2Y+BX3Y-BX4Y)*CPS+(BZ1Y-BZ2Y+BZ3Y-BZ4Y)*SPS) D(2,IY)=SPS*(BY1Y-BY2Y+BY3Y-BY4Y) D(3,IY)=SPS*((BZ1Y-BZ2Y+BZ3Y-BZ4Y)*CPS-(BX1Y-BX2Y+BX3Y-BX4Y)*SPS) C D(1,IZ)=SPS*((BX1Z+BX2Z-BX3Z-BX4Z)*CPS+(BZ1Z+BZ2Z-BZ3Z-BZ4Z)*SPS) D(2,IZ)=SPS*(BY1Z+BY2Z-BY3Z-BY4Z) D(3,IZ)=SPS*((BZ1Z+BZ2Z-BZ3Z-BZ4Z)*CPS-(BX1Z+BX2Z-BX3Z-BX4Z)*SPS) 2 CONTINUE C DO 3 I=1,5 ZD=ZZ(I+9) CALL DIPXYZ(XSM,Y,ZSM-ZD,BX1X,BY1X,BZ1X,BX1Y,BY1Y,BZ1Y,BX1Z,BY1Z, * BZ1Z) CALL DIPXYZ(XSM,Y,ZSM+ZD,BX2X,BY2X,BZ2X,BX2Y,BY2Y,BZ2Y,BX2Z,BY2Z, * BZ2Z) IX=58+I*2 IZ=IX+1 D(1,IX)=(BX1X-BX2X)*CPS+(BZ1X-BZ2X)*SPS D(2,IX)=BY1X-BY2X D(3,IX)=(BZ1X-BZ2X)*CPS-(BX1X-BX2X)*SPS C D(1,IZ)=(BX1Z+BX2Z)*CPS+(BZ1Z+BZ2Z)*SPS D(2,IZ)=BY1Z+BY2Z D(3,IZ)=(BZ1Z+BZ2Z)*CPS-(BX1Z+BX2Z)*SPS C IX=IX+10 IZ=IZ+10 D(1,IX)=SPS*((BX1X+BX2X)*CPS+(BZ1X+BZ2X)*SPS) D(2,IX)=SPS*(BY1X+BY2X) D(3,IX)=SPS*((BZ1X+BZ2X)*CPS-(BX1X+BX2X)*SPS) C D(1,IZ)=SPS*((BX1Z-BX2Z)*CPS+(BZ1Z-BZ2Z)*SPS) D(2,IZ)=SPS*(BY1Z-BY2Z) 3 D(3,IZ)=SPS*((BZ1Z-BZ2Z)*CPS-(BX1Z-BX2Z)*SPS) C RETURN END C C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C SUBROUTINE BIRK1SHLD(PS,X,Y,Z,BX,BY,BZ) C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C C The 64 linear parameters are amplitudes of the "box" harmonics. c The 16 nonlinear parameters are the scales Pi, and Qk entering the arguments C of sines/cosines and exponents in each of 32 cartesian harmonics c N.A. Tsyganenko, Spring 1994, adjusted for the Birkeland field Aug.22, 1995 c Revised June 12, 1996. C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C IMPLICIT REAL * 8 (A - H, O - Z) C DIMENSION A(80) DIMENSION P1(4),R1(4),Q1(4),S1(4),RP(4),RR(4),RQ(4),RS(4) C EQUIVALENCE (P1(1),A(65)),(R1(1),A(69)),(Q1(1),A(73)), * (S1(1),A(77)) C DATA A/1.174198045,-1.463820502,4.840161537,-3.674506864, * 82.18368896,-94.94071588,-4122.331796,4670.278676,-21.54975037, * 26.72661293,-72.81365728,44.09887902,40.08073706,-51.23563510, * 1955.348537,-1940.971550,794.0496433,-982.2441344,1889.837171, * -558.9779727,-1260.543238,1260.063802,-293.5942373,344.7250789, * -773.7002492,957.0094135,-1824.143669,520.7994379,1192.484774, * -1192.184565,89.15537624,-98.52042999,-0.8168777675E-01, * 0.4255969908E-01,0.3155237661,-0.3841755213,2.494553332, * -0.6571440817E-01,-2.765661310,0.4331001908,0.1099181537, * -0.6154126980E-01,-0.3258649260,0.6698439193,-5.542735524, * 0.1604203535,5.854456934,-0.8323632049,3.732608869,-3.130002153, * 107.0972607,-32.28483411,-115.2389298,54.45064360,-0.5826853320, * -3.582482231,-4.046544561,3.311978102,-104.0839563,30.26401293, * 97.29109008,-50.62370872,-296.3734955,127.7872523,5.303648988, * 10.40368955,69.65230348,466.5099509,1.645049286,3.825838190, * 11.66675599,558.9781177,1.826531343,2.066018073,25.40971369, * 990.2795225,2.319489258,4.555148484,9.691185703,591.8280358/ C BX=0.D0 BY=0.D0 BZ=0.D0 CPS=DCOS(PS) SPS=DSIN(PS) S3PS=4.D0*CPS**2-1.D0 C DO 11 I=1,4 RP(I)=1.D0/P1(I) RR(I)=1.D0/R1(I) RQ(I)=1.D0/Q1(I) 11 RS(I)=1.D0/S1(I) C L=0 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,4 CYPI=DCOS(Y*RP(I)) CYQI=DCOS(Y*RQ(I)) SYPI=DSIN(Y*RP(I)) SYQI=DSIN(Y*RQ(I)) C DO 3 K=1,4 SZRK=DSIN(Z*RR(K)) CZSK=DCOS(Z*RS(K)) CZRK=DCOS(Z*RR(K)) SZSK=DSIN(Z*RS(K)) SQPR=DSQRT(RP(I)**2+RR(K)**2) SQQS=DSQRT(RQ(I)**2+RS(K)**2) EPR=DEXP(X*SQPR) EQS=DEXP(X*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 IF (M.EQ.1) THEN IF (N.EQ.1) THEN HX=-SQPR*EPR*CYPI*SZRK HY=RP(I)*EPR*SYPI*SZRK HZ=-RR(K)*EPR*CYPI*CZRK ELSE HX=HX*CPS HY=HY*CPS HZ=HZ*CPS ENDIF ELSE IF (N.EQ.1) THEN HX=-SPS*SQQS*EQS*CYQI*CZSK HY=SPS*RQ(I)*EQS*SYQI*CZSK HZ=SPS*RS(K)*EQS*CYQI*SZSK ELSE HX=HX*S3PS HY=HY*S3PS HZ=HZ*S3PS ENDIF ENDIF L=L+1 c BX=BX+A(L)*HX BY=BY+A(L)*HY 4 BZ=BZ+A(L)*HZ 3 CONTINUE 2 CONTINUE 1 CONTINUE C RETURN END C C########################################################################## C SUBROUTINE BIRK2TOT_02(PS,X,Y,Z,BX,BY,BZ) C IMPLICIT REAL*8 (A-H,O-Z) C CALL BIRK2SHL(X,Y,Z,PS,WX,WY,WZ) CALL R2_BIRK(X,Y,Z,PS,HX,HY,HZ) BX=WX+HX BY=WY+HY BZ=WZ+HZ RETURN END C C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C C THIS CODE IS FOR THE FIELD FROM 2x2x2=8 "CARTESIAN" HARMONICS C SUBROUTINE BIRK2SHL(X,Y,Z,PS,HX,HY,HZ) C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C The model parameters are provided to this module via common-block /A/. C The 16 linear parameters enter in pairs in the amplitudes of the c "cartesian" harmonics. c The 8 nonlinear parameters are the scales Pi,Ri,Qi,and Si entering the c arguments of exponents, sines, and cosines in each of the 8 "Cartesian" c harmonics C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C IMPLICIT REAL * 8 (A - H, O - Z) C DIMENSION P(2),R(2),Q(2),S(2) DIMENSION A(24) C EQUIVALENCE(P(1),A(17)),(R(1),A(19)),(Q(1),A(21)),(S(1),A(23)) DATA A/-111.6371348,124.5402702,110.3735178,-122.0095905, * 111.9448247,-129.1957743,-110.7586562,126.5649012,-0.7865034384, * -0.2483462721,0.8026023894,0.2531397188,10.72890902,0.8483902118, * -10.96884315,-0.8583297219,13.85650567,14.90554500,10.21914434, * 10.09021632,6.340382460,14.40432686,12.71023437,12.83966657/ C CPS=DCOS(PS) SPS=DSIN(PS) S3PS=4.D0*CPS**2-1.D0 ! THIS IS SIN(3*PS)/SIN(PS) C HX=0.D0 HY=0.D0 HZ=0.D0 L=0 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,2 CYPI=DCOS(Y/P(I)) CYQI=DCOS(Y/Q(I)) SYPI=DSIN(Y/P(I)) SYQI=DSIN(Y/Q(I)) C DO 3 K=1,2 SZRK=DSIN(Z/R(K)) CZSK=DCOS(Z/S(K)) CZRK=DCOS(Z/R(K)) SZSK=DSIN(Z/S(K)) SQPR=DSQRT(1.D0/P(I)**2+1.D0/R(K)**2) SQQS=DSQRT(1.D0/Q(I)**2+1.D0/S(K)**2) EPR=DEXP(X*SQPR) EQS=DEXP(X*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 C L=L+1 IF (M.EQ.1) THEN IF (N.EQ.1) THEN DX=-SQPR*EPR*CYPI*SZRK DY=EPR/P(I)*SYPI*SZRK DZ=-EPR/R(K)*CYPI*CZRK HX=HX+A(L)*DX HY=HY+A(L)*DY HZ=HZ+A(L)*DZ ELSE DX=DX*CPS DY=DY*CPS DZ=DZ*CPS HX=HX+A(L)*DX HY=HY+A(L)*DY HZ=HZ+A(L)*DZ ENDIF ELSE IF (N.EQ.1) THEN DX=-SPS*SQQS*EQS*CYQI*CZSK DY=SPS*EQS/Q(I)*SYQI*CZSK DZ=SPS*EQS/S(K)*CYQI*SZSK HX=HX+A(L)*DX HY=HY+A(L)*DY HZ=HZ+A(L)*DZ ELSE DX=DX*S3PS DY=DY*S3PS DZ=DZ*S3PS HX=HX+A(L)*DX HY=HY+A(L)*DY HZ=HZ+A(L)*DZ ENDIF ENDIF c 4 CONTINUE 3 CONTINUE 2 CONTINUE 1 CONTINUE C RETURN END c******************************************************************** C SUBROUTINE R2_BIRK(X,Y,Z,PS,BX,BY,BZ) C C RETURNS THE MODEL FIELD FOR THE REGION 2 BIRKELAND CURRENT/PARTIAL RC C (WITHOUT SHIELDING FIELD) C IMPLICIT REAL*8 (A-H,O-Z) SAVE PSI,CPS,SPS DATA DELARG/0.030D0/,DELARG1/0.015D0/,PSI/10.D0/ C IF (DABS(PSI-PS).GT.1.D-10) THEN PSI=PS CPS=DCOS(PS) SPS=DSIN(PS) ENDIF C XSM=X*CPS-Z*SPS ZSM=Z*CPS+X*SPS C XKS=XKSI(XSM,Y,ZSM) IF (XKS.LT.-(DELARG+DELARG1)) THEN CALL R2OUTER(XSM,Y,ZSM,BXSM,BY,BZSM) BXSM=-BXSM*0.02 ! ALL COMPONENTS ARE MULTIPLIED BY THE BY=-BY*0.02 ! FACTOR -0.02, IN ORDER TO NORMALIZE THE BZSM=-BZSM*0.02 ! FIELD (SO THAT Bz=-1 nT at X=-5.3 RE, Y=Z=0) ENDIF IF (XKS.GE.-(DELARG+DELARG1).AND.XKS.LT.-DELARG+DELARG1) THEN CALL R2OUTER(XSM,Y,ZSM,BXSM1,BY1,BZSM1) CALL R2SHEET(XSM,Y,ZSM,BXSM2,BY2,BZSM2) F2=-0.02*TKSI(XKS,-DELARG,DELARG1) F1=-0.02-F2 BXSM=BXSM1*F1+BXSM2*F2 BY=BY1*F1+BY2*F2 BZSM=BZSM1*F1+BZSM2*F2 ENDIF IF (XKS.GE.-DELARG+DELARG1.AND.XKS.LT.DELARG-DELARG1) THEN CALL R2SHEET(XSM,Y,ZSM,BXSM,BY,BZSM) BXSM=-BXSM*0.02 BY=-BY*0.02 BZSM=-BZSM*0.02 ENDIF IF (XKS.GE.DELARG-DELARG1.AND.XKS.LT.DELARG+DELARG1) THEN CALL R2INNER(XSM,Y,ZSM,BXSM1,BY1,BZSM1) CALL R2SHEET(XSM,Y,ZSM,BXSM2,BY2,BZSM2) F1=-0.02*TKSI(XKS,DELARG,DELARG1) F2=-0.02-F1 BXSM=BXSM1*F1+BXSM2*F2 BY=BY1*F1+BY2*F2 BZSM=BZSM1*F1+BZSM2*F2 ENDIF IF (XKS.GE.DELARG+DELARG1) THEN CALL R2INNER(XSM,Y,ZSM,BXSM,BY,BZSM) BXSM=-BXSM*0.02 BY=-BY*0.02 BZSM=-BZSM*0.02 ENDIF C BX=BXSM*CPS+BZSM*SPS BZ=BZSM*CPS-BXSM*SPS C RETURN END C C**************************************************************** c SUBROUTINE R2INNER (X,Y,Z,BX,BY,BZ) C C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION CBX(5),CBY(5),CBZ(5) C DATA PL1,PL2,PL3,PL4,PL5,PL6,PL7,PL8/154.185,-2.12446,.601735E-01, * -.153954E-02,.355077E-04,29.9996,262.886,99.9132/ DATA PN1,PN2,PN3,PN4,PN5,PN6,PN7,PN8/-8.1902,6.5239,5.504,7.7815, * .8573,3.0986,.0774,-.038/ C CALL BCONIC(X,Y,Z,CBX,CBY,CBZ,5) C C NOW INTRODUCE ONE 4-LOOP SYSTEM: C CALL LOOPS4(X,Y,Z,DBX8,DBY8,DBZ8,PN1,PN2,PN3,PN4,PN5,PN6) C CALL DIPDISTR(X-PN7,Y,Z,DBX6,DBY6,DBZ6,0) CALL DIPDISTR(X-PN8,Y,Z,DBX7,DBY7,DBZ7,1) C NOW COMPUTE THE FIELD COMPONENTS: BX=PL1*CBX(1)+PL2*CBX(2)+PL3*CBX(3)+PL4*CBX(4)+PL5*CBX(5) * +PL6*DBX6+PL7*DBX7+PL8*DBX8 BY=PL1*CBY(1)+PL2*CBY(2)+PL3*CBY(3)+PL4*CBY(4)+PL5*CBY(5) * +PL6*DBY6+PL7*DBY7+PL8*DBY8 BZ=PL1*CBZ(1)+PL2*CBZ(2)+PL3*CBZ(3)+PL4*CBZ(4)+PL5*CBZ(5) * +PL6*DBZ6+PL7*DBZ7+PL8*DBZ8 C RETURN END C----------------------------------------------------------------------- SUBROUTINE BCONIC(X,Y,Z,CBX,CBY,CBZ,NMAX) C c "CONICAL" HARMONICS c IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION CBX(NMAX),CBY(NMAX),CBZ(NMAX) RO2=X**2+Y**2 RO=SQRT(RO2) C CF=X/RO SF=Y/RO CFM1=1.D0 SFM1=0.D0 C R2=RO2+Z**2 R=DSQRT(R2) C=Z/R S=RO/R CH=DSQRT(0.5D0*(1.D0+C)) SH=DSQRT(0.5D0*(1.D0-C)) TNHM1=1.D0 CNHM1=1.D0 TNH=SH/CH CNH=1.D0/TNH C DO 1 M=1,NMAX CFM=CFM1*CF-SFM1*SF SFM=CFM1*SF+SFM1*CF CFM1=CFM SFM1=SFM TNHM=TNHM1*TNH CNHM=CNHM1*CNH BT=M*CFM/(R*S)*(TNHM+CNHM) BF=-0.5D0*M*SFM/R*(TNHM1/CH**2-CNHM1/SH**2) TNHM1=TNHM CNHM1=CNHM CBX(M)=BT*C*CF-BF*SF CBY(M)=BT*C*SF+BF*CF 1 CBZ(M)=-BT*S C RETURN END C------------------------------------------------------------------- C SUBROUTINE DIPDISTR(X,Y,Z,BX,BY,BZ,MODE) C C RETURNS FIELD COMPONENTS FROM A LINEAR DISTRIBUTION OF DIPOLAR SOURCES C ON THE Z-AXIS. THE PARAMETER MODE DEFINES HOW THE DIPOLE STRENGTH C VARIES ALONG THE Z-AXIS: MODE=0 IS FOR A STEP-FUNCTION (Mx=const > 0 c FOR Z > 0, AND Mx=-const < 0 FOR Z < 0) C WHILE MODE=1 IS FOR A LINEAR VARIATION OF THE DIPOLE MOMENT DENSITY C SEE NB#3, PAGE 53 FOR DETAILS. C C C INPUT: X,Y,Z OF A POINT OF SPACE, AND MODE C IMPLICIT REAL*8 (A-H,O-Z) X2=X*X RHO2=X2+Y*Y R2=RHO2+Z*Z R3=R2*DSQRT(R2) IF (MODE.EQ.0) THEN BX=Z/RHO2**2*(R2*(Y*Y-X2)-RHO2*X2)/R3 BY=-X*Y*Z/RHO2**2*(2.D0*R2+RHO2)/R3 BZ=X/R3 ELSE BX=Z/RHO2**2*(Y*Y-X2) BY=-2.D0*X*Y*Z/RHO2**2 BZ=X/RHO2 ENDIF RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C SUBROUTINE R2OUTER (X,Y,Z,BX,BY,BZ) C IMPLICIT REAL*8 (A-H,O-Z) C DATA PL1,PL2,PL3,PL4,PL5/-34.105,-2.00019,628.639,73.4847,12.5162/ DATA PN1,PN2,PN3,PN4,PN5,PN6,PN7,PN8,PN9,PN10,PN11,PN12,PN13,PN14, * PN15,PN16,PN17 /.55,.694,.0031,1.55,2.8,.1375,-.7,.2,.9625, * -2.994,2.925,-1.775,4.3,-.275,2.7,.4312,1.55/ c C THREE PAIRS OF CROSSED LOOPS: C CALL CROSSLP(X,Y,Z,DBX1,DBY1,DBZ1,PN1,PN2,PN3) CALL CROSSLP(X,Y,Z,DBX2,DBY2,DBZ2,PN4,PN5,PN6) CALL CROSSLP(X,Y,Z,DBX3,DBY3,DBZ3,PN7,PN8,PN9) C C NOW AN EQUATORIAL LOOP ON THE NIGHTSIDE C CALL CIRCLE(X-PN10,Y,Z,PN11,DBX4,DBY4,DBZ4) c c NOW A 4-LOOP SYSTEM ON THE NIGHTSIDE c CALL LOOPS4(X,Y,Z,DBX5,DBY5,DBZ5,PN12,PN13,PN14,PN15,PN16,PN17) C--------------------------------------------------------------------- C NOW COMPUTE THE FIELD COMPONENTS: BX=PL1*DBX1+PL2*DBX2+PL3*DBX3+PL4*DBX4+PL5*DBX5 BY=PL1*DBY1+PL2*DBY2+PL3*DBY3+PL4*DBY4+PL5*DBY5 BZ=PL1*DBZ1+PL2*DBZ2+PL3*DBZ3+PL4*DBZ4+PL5*DBZ5 RETURN END C C-------------------------------------------------------------------- C SUBROUTINE LOOPS4(X,Y,Z,BX,BY,BZ,XC,YC,ZC,R,THETA,PHI) C C RETURNS FIELD COMPONENTS FROM A SYSTEM OF 4 CURRENT LOOPS, POSITIONED C SYMMETRICALLY WITH RESPECT TO NOON-MIDNIGHT MERIDIAN AND EQUATORIAL C PLANES. C INPUT: X,Y,Z OF A POINT OF SPACE C XC,YC,ZC (YC > 0 AND ZC > 0) - POSITION OF THE CENTER OF THE C 1ST-QUADRANT LOOP C R - LOOP RADIUS (THE SAME FOR ALL FOUR) C THETA, PHI - SPECIFY THE ORIENTATION OF THE NORMAL OF THE 1ST LOOP c ----------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) C CT=DCOS(THETA) ST=DSIN(THETA) CP=DCOS(PHI) SP=DSIN(PHI) C------------------------------------1ST QUADRANT: XS=(X-XC)*CP+(Y-YC)*SP YSS=(Y-YC)*CP-(X-XC)*SP ZS=Z-ZC XSS=XS*CT-ZS*ST ZSS=ZS*CT+XS*ST CALL CIRCLE(XSS,YSS,ZSS,R,BXSS,BYS,BZSS) BXS=BXSS*CT+BZSS*ST BZ1=BZSS*CT-BXSS*ST BX1=BXS*CP-BYS*SP BY1=BXS*SP+BYS*CP C-------------------------------------2nd QUADRANT: XS=(X-XC)*CP-(Y+YC)*SP YSS=(Y+YC)*CP+(X-XC)*SP ZS=Z-ZC XSS=XS*CT-ZS*ST ZSS=ZS*CT+XS*ST CALL CIRCLE(XSS,YSS,ZSS,R,BXSS,BYS,BZSS) BXS=BXSS*CT+BZSS*ST BZ2=BZSS*CT-BXSS*ST BX2=BXS*CP+BYS*SP BY2=-BXS*SP+BYS*CP C-------------------------------------3RD QUADRANT: XS=-(X-XC)*CP+(Y+YC)*SP YSS=-(Y+YC)*CP-(X-XC)*SP ZS=Z+ZC XSS=XS*CT-ZS*ST ZSS=ZS*CT+XS*ST CALL CIRCLE(XSS,YSS,ZSS,R,BXSS,BYS,BZSS) BXS=BXSS*CT+BZSS*ST BZ3=BZSS*CT-BXSS*ST BX3=-BXS*CP-BYS*SP BY3=BXS*SP-BYS*CP C-------------------------------------4TH QUADRANT: XS=-(X-XC)*CP-(Y-YC)*SP YSS=-(Y-YC)*CP+(X-XC)*SP ZS=Z+ZC XSS=XS*CT-ZS*ST ZSS=ZS*CT+XS*ST CALL CIRCLE(XSS,YSS,ZSS,R,BXSS,BYS,BZSS) BXS=BXSS*CT+BZSS*ST BZ4=BZSS*CT-BXSS*ST BX4=-BXS*CP+BYS*SP BY4=-BXS*SP-BYS*CP BX=BX1+BX2+BX3+BX4 BY=BY1+BY2+BY3+BY4 BZ=BZ1+BZ2+BZ3+BZ4 RETURN END C C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ C SUBROUTINE R2SHEET(X,Y,Z,BX,BY,BZ) C IMPLICIT REAL*8 (A-H,O-Z) C DATA PNONX1,PNONX2,PNONX3,PNONX4,PNONX5,PNONX6,PNONX7,PNONX8, * PNONY1,PNONY2,PNONY3,PNONY4,PNONY5,PNONY6,PNONY7,PNONY8, * PNONZ1,PNONZ2,PNONZ3,PNONZ4,PNONZ5,PNONZ6,PNONZ7,PNONZ8 */-19.0969D0,-9.28828D0,-0.129687D0,5.58594D0,22.5055D0, * 0.483750D-01,0.396953D-01,0.579023D-01,-13.6750D0,-6.70625D0, * 2.31875D0,11.4062D0,20.4562D0,0.478750D-01,0.363750D-01, * 0.567500D-01,-16.7125D0,-16.4625D0,-0.1625D0,5.1D0,23.7125D0, * 0.355625D-01,0.318750D-01,0.538750D-01/ C C DATA A1,A2,A3,A4,A5,A6,A7,A8,A9,A10,A11,A12,A13,A14,A15,A16,A17, * A18,A19,A20,A21,A22,A23,A24,A25,A26,A27,A28,A29,A30,A31,A32,A33, * A34,A35,A36,A37,A38,A39,A40,A41,A42,A43,A44,A45,A46,A47,A48,A49, * A50,A51,A52,A53,A54,A55,A56,A57,A58,A59,A60,A61,A62,A63,A64,A65, * A66,A67,A68,A69,A70,A71,A72,A73,A74,A75,A76,A77,A78,A79,A80 * /8.07190D0,-7.39582D0,-7.62341D0,0.684671D0,-13.5672D0,11.6681D0, * 13.1154,-0.890217D0,7.78726D0,-5.38346D0,-8.08738D0,0.609385D0, * -2.70410D0, 3.53741D0,3.15549D0,-1.11069D0,-8.47555D0,0.278122D0, * 2.73514D0,4.55625D0,13.1134D0,1.15848D0,-3.52648D0,-8.24698D0, * -6.85710D0,-2.81369D0, 2.03795D0, 4.64383D0,2.49309D0,-1.22041D0, * -1.67432D0,-0.422526D0,-5.39796D0,7.10326D0,5.53730D0,-13.1918D0, * 4.67853D0,-7.60329D0,-2.53066D0, 7.76338D0, 5.60165D0,5.34816D0, * -4.56441D0,7.05976D0,-2.62723D0,-0.529078D0,1.42019D0,-2.93919D0, * 55.6338D0,-1.55181D0,39.8311D0,-80.6561D0,-46.9655D0,32.8925D0, * -6.32296D0,19.7841D0,124.731D0,10.4347D0,-30.7581D0,102.680D0, * -47.4037D0,-3.31278D0,9.37141D0,-50.0268D0,-533.319D0,110.426D0, * 1000.20D0,-1051.40D0, 1619.48D0,589.855D0,-1462.73D0,1087.10D0, * -1994.73D0,-1654.12D0,1263.33D0,-260.210D0,1424.84D0,1255.71D0, * -956.733D0, 219.946D0/ C C DATA B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11,B12,B13,B14,B15,B16,B17, * B18,B19,B20,B21,B22,B23,B24,B25,B26,B27,B28,B29,B30,B31,B32,B33, * B34,B35,B36,B37,B38,B39,B40,B41,B42,B43,B44,B45,B46,B47,B48,B49, * B50,B51,B52,B53,B54,B55,B56,B57,B58,B59,B60,B61,B62,B63,B64,B65, * B66,B67,B68,B69,B70,B71,B72,B73,B74,B75,B76,B77,B78,B79,B80 */-9.08427D0,10.6777D0,10.3288D0,-0.969987D0,6.45257D0,-8.42508D0, * -7.97464D0,1.41996D0,-1.92490D0,3.93575D0,2.83283D0,-1.48621D0, *0.244033D0,-0.757941D0,-0.386557D0,0.344566D0,9.56674D0,-2.5365D0, * -3.32916D0,-5.86712D0,-6.19625D0,1.83879D0,2.52772D0,4.34417D0, * 1.87268D0,-2.13213D0,-1.69134D0,-.176379D0,-.261359D0,.566419D0, * 0.3138D0,-0.134699D0,-3.83086D0,-8.4154D0,4.77005D0,-9.31479D0, * 37.5715D0,19.3992D0,-17.9582D0,36.4604D0,-14.9993D0,-3.1442D0, * 6.17409D0,-15.5519D0,2.28621D0,-0.891549D-2,-.462912D0,2.47314D0, * 41.7555D0,208.614D0,-45.7861D0,-77.8687D0,239.357D0,-67.9226D0, * 66.8743D0,238.534D0,-112.136D0,16.2069D0,-40.4706D0,-134.328D0, * 21.56D0,-0.201725D0,2.21D0,32.5855D0,-108.217D0,-1005.98D0, * 585.753D0,323.668D0,-817.056D0,235.750D0,-560.965D0,-576.892D0, * 684.193D0,85.0275D0,168.394D0,477.776D0,-289.253D0,-123.216D0, * 75.6501D0,-178.605D0/ C DATA C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12,C13,C14,C15,C16,C17, * C18,C19,C20,C21,C22,C23,C24,C25,C26,C27,C28,C29,C30,C31,C32,C33, * C34,C35,C36,C37,C38,C39,C40,C41,C42,C43,C44,C45,C46,C47,C48,C49, * C50,C51,C52,C53,C54,C55,C56,C57,C58,C59,C60,C61,C62,C63,C64,C65, * C66,C67,C68,C69,C70,C71,C72,C73,C74,C75,C76,C77,C78,C79,C80 * / 1167.61D0,-917.782D0,-1253.2D0,-274.128D0,-1538.75D0,1257.62D0, * 1745.07D0,113.479D0,393.326D0,-426.858D0,-641.1D0,190.833D0, * -29.9435D0,-1.04881D0,117.125D0,-25.7663D0,-1168.16D0,910.247D0, * 1239.31D0,289.515D0,1540.56D0,-1248.29D0,-1727.61D0,-131.785D0, * -394.577D0,426.163D0,637.422D0,-187.965D0,30.0348D0,0.221898D0, * -116.68D0,26.0291D0,12.6804D0,4.84091D0,1.18166D0,-2.75946D0, * -17.9822D0,-6.80357D0,-1.47134D0,3.02266D0,4.79648D0,0.665255D0, * -0.256229D0,-0.857282D-1,-0.588997D0,0.634812D-1,0.164303D0, * -0.15285D0,22.2524D0,-22.4376D0,-3.85595D0,6.07625D0,-105.959D0, * -41.6698D0,0.378615D0,1.55958D0,44.3981D0,18.8521D0,3.19466D0, * 5.89142D0,-8.63227D0,-2.36418D0,-1.027D0,-2.31515D0,1035.38D0, * 2040.66D0,-131.881D0,-744.533D0,-3274.93D0,-4845.61D0,482.438D0, * 1567.43D0,1354.02D0,2040.47D0,-151.653D0,-845.012D0,-111.723D0, * -265.343D0,-26.1171D0,216.632D0/ C c------------------------------------------------------------------ C XKS=XKSI(X,Y,Z) ! variation across the current sheet T1X=XKS/DSQRT(XKS**2+PNONX6**2) T2X=PNONX7**3/DSQRT(XKS**2+PNONX7**2)**3 T3X=XKS/DSQRT(XKS**2+PNONX8**2)**5 *3.493856D0*PNONX8**4 C T1Y=XKS/DSQRT(XKS**2+PNONY6**2) T2Y=PNONY7**3/DSQRT(XKS**2+PNONY7**2)**3 T3Y=XKS/DSQRT(XKS**2+PNONY8**2)**5 *3.493856D0*PNONY8**4 C T1Z=XKS/DSQRT(XKS**2+PNONZ6**2) T2Z=PNONZ7**3/DSQRT(XKS**2+PNONZ7**2)**3 T3Z=XKS/DSQRT(XKS**2+PNONZ8**2)**5 *3.493856D0*PNONZ8**4 C RHO2=X*X+Y*Y R=DSQRT(RHO2+Z*Z) RHO=DSQRT(RHO2) C C1P=X/RHO S1P=Y/RHO S2P=2.D0*S1P*C1P C2P=C1P*C1P-S1P*S1P S3P=S2P*C1P+C2P*S1P C3P=C2P*C1P-S2P*S1P S4P=S3P*C1P+C3P*S1P CT=Z/R ST=RHO/R C S1=FEXP(CT,PNONX1) S2=FEXP(CT,PNONX2) S3=FEXP(CT,PNONX3) S4=FEXP(CT,PNONX4) S5=FEXP(CT,PNONX5) C C NOW COMPUTE THE GSM FIELD COMPONENTS: C C BX=S1*((A1+A2*T1X+A3*T2X+A4*T3X) * +C1P*(A5+A6*T1X+A7*T2X+A8*T3X) * +C2P*(A9+A10*T1X+A11*T2X+A12*T3X) * +C3P*(A13+A14*T1X+A15*T2X+A16*T3X)) * +S2*((A17+A18*T1X+A19*T2X+A20*T3X) * +C1P*(A21+A22*T1X+A23*T2X+A24*T3X) * +C2P*(A25+A26*T1X+A27*T2X+A28*T3X) * +C3P*(A29+A30*T1X+A31*T2X+A32*T3X)) * +S3*((A33+A34*T1X+A35*T2X+A36*T3X) * +C1P*(A37+A38*T1X+A39*T2X+A40*T3X) * +C2P*(A41+A42*T1X+A43*T2X+A44*T3X) * +C3P*(A45+A46*T1X+A47*T2X+A48*T3X)) * +S4*((A49+A50*T1X+A51*T2X+A52*T3X) * +C1P*(A53+A54*T1X+A55*T2X+A56*T3X) * +C2P*(A57+A58*T1X+A59*T2X+A60*T3X) * +C3P*(A61+A62*T1X+A63*T2X+A64*T3X)) * +S5*((A65+A66*T1X+A67*T2X+A68*T3X) * +C1P*(A69+A70*T1X+A71*T2X+A72*T3X) * +C2P*(A73+A74*T1X+A75*T2X+A76*T3X) * +C3P*(A77+A78*T1X+A79*T2X+A80*T3X)) C C S1=FEXP(CT,PNONY1) S2=FEXP(CT,PNONY2) S3=FEXP(CT,PNONY3) S4=FEXP(CT,PNONY4) S5=FEXP(CT,PNONY5) C BY=S1*(S1P*(B1+B2*T1Y+B3*T2Y+B4*T3Y) * +S2P*(B5+B6*T1Y+B7*T2Y+B8*T3Y) * +S3P*(B9+B10*T1Y+B11*T2Y+B12*T3Y) * +S4P*(B13+B14*T1Y+B15*T2Y+B16*T3Y)) * +S2*(S1P*(B17+B18*T1Y+B19*T2Y+B20*T3Y) * +S2P*(B21+B22*T1Y+B23*T2Y+B24*T3Y) * +S3P*(B25+B26*T1Y+B27*T2Y+B28*T3Y) * +S4P*(B29+B30*T1Y+B31*T2Y+B32*T3Y)) * +S3*(S1P*(B33+B34*T1Y+B35*T2Y+B36*T3Y) * +S2P*(B37+B38*T1Y+B39*T2Y+B40*T3Y) * +S3P*(B41+B42*T1Y+B43*T2Y+B44*T3Y) * +S4P*(B45+B46*T1Y+B47*T2Y+B48*T3Y)) * +S4*(S1P*(B49+B50*T1Y+B51*T2Y+B52*T3Y) * +S2P*(B53+B54*T1Y+B55*T2Y+B56*T3Y) * +S3P*(B57+B58*T1Y+B59*T2Y+B60*T3Y) * +S4P*(B61+B62*T1Y+B63*T2Y+B64*T3Y)) * +S5*(S1P*(B65+B66*T1Y+B67*T2Y+B68*T3Y) * +S2P*(B69+B70*T1Y+B71*T2Y+B72*T3Y) * +S3P*(B73+B74*T1Y+B75*T2Y+B76*T3Y) * +S4P*(B77+B78*T1Y+B79*T2Y+B80*T3Y)) C S1=FEXP1(CT,PNONZ1) S2=FEXP1(CT,PNONZ2) S3=FEXP1(CT,PNONZ3) S4=FEXP1(CT,PNONZ4) S5=FEXP1(CT,PNONZ5) C BZ=S1*((C1+C2*T1Z+C3*T2Z+C4*T3Z) * +C1P*(C5+C6*T1Z+C7*T2Z+C8*T3Z) * +C2P*(C9+C10*T1Z+C11*T2Z+C12*T3Z) * +C3P*(C13+C14*T1Z+C15*T2Z+C16*T3Z)) * +S2*((C17+C18*T1Z+C19*T2Z+C20*T3Z) * +C1P*(C21+C22*T1Z+C23*T2Z+C24*T3Z) * +C2P*(C25+C26*T1Z+C27*T2Z+C28*T3Z) * +C3P*(C29+C30*T1Z+C31*T2Z+C32*T3Z)) * +S3*((C33+C34*T1Z+C35*T2Z+C36*T3Z) * +C1P*(C37+C38*T1Z+C39*T2Z+C40*T3Z) * +C2P*(C41+C42*T1Z+C43*T2Z+C44*T3Z) * +C3P*(C45+C46*T1Z+C47*T2Z+C48*T3Z)) * +S4*((C49+C50*T1Z+C51*T2Z+C52*T3Z) * +C1P*(C53+C54*T1Z+C55*T2Z+C56*T3Z) * +C2P*(C57+C58*T1Z+C59*T2Z+C60*T3Z) * +C3P*(C61+C62*T1Z+C63*T2Z+C64*T3Z)) * +S5*((C65+C66*T1Z+C67*T2Z+C68*T3Z) * +C1P*(C69+C70*T1Z+C71*T2Z+C72*T3Z) * +C2P*(C73+C74*T1Z+C75*T2Z+C76*T3Z) * +C3P*(C77+C78*T1Z+C79*T2Z+C80*T3Z)) C RETURN END C C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ DOUBLE PRECISION FUNCTION XKSI(X,Y,Z) C IMPLICIT REAL*8 (A-H,O-Z) C C A11 - C72, R0, and DR below ARE STRETCH PARAMETERS (P.26-27, NB# 3), C DATA A11A12,A21A22,A41A42,A51A52,A61A62,B11B12,B21B22,C61C62, * C71C72,R0,DR /0.305662,-0.383593,0.2677733,-0.097656,-0.636034, * -0.359862,0.424706,-0.126366,0.292578,1.21563,7.50937/ DATA TNOON,DTETA/0.3665191,0.09599309/ ! Correspond to noon and midnight C latitudes 69 and 63.5 degs, resp. DR2=DR*DR C X2=X*X Y2=Y*Y Z2=Z*Z XY=X*Y XYZ=XY*Z R2=X2+Y2+Z2 R=DSQRT(R2) R3=R2*R R4=R2*R2 XR=X/R YR=Y/R ZR=Z/R C IF (R.LT.R0) THEN PR=0.D0 ELSE PR=DSQRT((R-R0)**2+DR2)-DR ENDIF C F=X+PR*(A11A12+A21A22*XR+A41A42*XR*XR+A51A52*YR*YR+ * A61A62*ZR*ZR) G=Y+PR*(B11B12*YR+B21B22*XR*YR) H=Z+PR*(C61C62*ZR+C71C72*XR*ZR) G2=G*G C FGH=F**2+G2+H**2 FGH32=DSQRT(FGH)**3 FCHSG2=F**2+G2 IF (FCHSG2.LT.1.D-5) THEN XKSI=-1.D0 ! THIS IS JUST FOR ELIMINATING PROBLEMS RETURN ! ON THE Z-AXIS ENDIF SQFCHSG2=DSQRT(FCHSG2) ALPHA=FCHSG2/FGH32 THETA=TNOON+0.5D0*DTETA*(1.D0-F/SQFCHSG2) PHI=DSIN(THETA)**2 C XKSI=ALPHA-PHI C RETURN END C C-------------------------------------------------------------------- C FUNCTION FEXP(S,A) IMPLICIT REAL*8 (A-H,O-Z) DATA E/2.718281828459D0/ IF (A.LT.0.D0) FEXP=DSQRT(-2.D0*A*E)*S*DEXP(A*S*S) IF (A.GE.0.D0) FEXP=S*DEXP(A*(S*S-1.D0)) RETURN END C C----------------------------------------------------------------------- FUNCTION FEXP1(S,A) IMPLICIT REAL*8 (A-H,O-Z) IF (A.LE.0.D0) FEXP1=DEXP(A*S*S) IF (A.GT.0.D0) FEXP1=DEXP(A*(S*S-1.D0)) RETURN END C C************************************************************************ C DOUBLE PRECISION FUNCTION TKSI(XKSI,XKS0,DXKSI) IMPLICIT REAL*8 (A-H,O-Z) SAVE M,TDZ3 DATA M/0/ C IF (M.EQ.0) THEN TDZ3=2.*DXKSI**3 M=1 ENDIF C IF (XKSI-XKS0.LT.-DXKSI) TKSII=0. IF (XKSI-XKS0.GE.DXKSI) TKSII=1. C IF (XKSI.GE.XKS0-DXKSI.AND.XKSI.LT.XKS0) THEN BR3=(XKSI-XKS0+DXKSI)**3 TKSII=1.5*BR3/(TDZ3+BR3) ENDIF C IF (XKSI.GE.XKS0.AND.XKSI.LT.XKS0+DXKSI) THEN BR3=(XKSI-XKS0-DXKSI)**3 TKSII=1.+1.5*BR3/(TDZ3-BR3) ENDIF TKSI=TKSII END C C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ c SUBROUTINE DIPOLE(PS,X,Y,Z,BX,BY,BZ) C C CALCULATES GSM COMPONENTS OF GEODIPOLE FIELD WITH THE DIPOLE MOMENT C CORRESPONDING TO THE EPOCH OF 1980. C------------INPUT PARAMETERS: C PS - GEODIPOLE TILT ANGLE IN RADIANS, X,Y,Z - GSM COORDINATES IN RE C------------OUTPUT PARAMETERS: C BX,BY,BZ - FIELD COMPONENTS IN GSM SYSTEM, IN NANOTESLA. C C C WRITEN BY: N. A. TSYGANENKO IMPLICIT REAL*8 (A-H,O-Z) DATA M,PSI/0,5./ SAVE M,PSI,SPS,CPS IF(M.EQ.1.AND.ABS(PS-PSI).LT.1.E-5) GOTO 1 SPS=SIN(PS) CPS=COS(PS) PSI=PS M=1 1 P=X**2 U=Z**2 V=3.*Z*X T=Y**2 Q=30574./SQRT(P+T+U)**5 BX=Q*((T+U-2.*P)*SPS-V*CPS) BY=-3.*Y*Q*(X*SPS+Z*CPS) BZ=Q*((P+T-2.*U)*CPS-V*SPS) RETURN END