!*************************************************************************************************** ! Copyright 2004 N. Tsyganenko ! ! This file is part of IRBEM-LIB. ! ! IRBEM-LIB is free software: you can redistribute it and/or modify ! it under the terms of the GNU Lesser General Public License as published by ! the Free Software Foundation, either version 3 of the License, or ! (at your option) any later version. ! ! IRBEM-LIB is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU Lesser General Public License for more details. ! ! You should have received a copy of the GNU Lesser General Public License ! along with IRBEM-LIB. If not, see . ! c==================================================================================== c c SUBROUTINE T04_s (PARMOD,X,Y,Z,BX,BY,BZ) c c ASSEMBLED: MARCH 25, 2004; UPDATED: AUGUST 2 & 31, DECEMBER 27, 2004. c A BUG ELIMINATED MARCH 14, 2005 (might cause compilation problems with c some Fortran compilers) C Coefficients in AP_04 DATA, FULL_RC DATA C_SY and T04_S DATA A updated C May 23rd 2012. See http://geo.phys.spbu.ru/~tsyganenko/modeling.html C for details of the bug (corrected by N. Tsyganenko on June 24th 2006). 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_04 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 REAL*8 PARMOD(10),PS,X,Y,Z,BX,BY,BZ,tilt 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 COMMON /dip_ang/tilt 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/ DATA IOPGEN,IOPTT,IOPB,IOPR/0,0,0,0/ C ps=tilt*4.D0*ATAN(1.D0)/180.d0 c PDYN=PARMOD(1) DST_AST=PARMOD(2)*0.8D0-13.D0*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_04 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_tsyg04/ RH0 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.5D0 ! TAIL HINGING DISTANCE c G= 35.0D0 ! TAIL WARPING PARAMETER 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_04 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) IF (DD.GT.1.D-6) GOTO 1 C END OF ITERATIVE SEARCH RHO2=Y**2+ZSS**2 ASQ=AM**2 XMXM=AM+XSS-X0 IF (XMXM.LT.0.D0) XMXM=0.D0 ! 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.D0*ASQ*AXX0))/ *(2.D0*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_04(XX,YY,ZZ,PS,CFX,CFY,CFZ) ! DIPOLE_04 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.D0 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.7D0 CALL DEFORMED_04 (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_04 (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_04(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.5D0*(1.D0-(SIGMA-S0)/DSIG) FEXT=0.5D0*(1.D0+(SIGMA-S0)/DSIG) C CALL DIPOLE_04 (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_04 (PS,X,Y,Z,QX,QY,QZ) BX=OIMFX-QX BY=OIMFY-QY BZ=OIMFZ-QZ ENDIF C END c C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C SUBROUTINE SHLCAR3X3_04(X,Y,Z,PS,BX,BY,BZ) C C THIS S/R RETURNS THE SHIELDING FIELD FOR THE EARTH'S DIPOLE_04, 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.2327248D0,895.8011176D0,817.6208321D0,-845.5880889D0, *-83.73539535D0,86.58542841D0,336.8781402D0,-329.3619944D0, *-311.2947120D0, *308.6011161D0,31.94469304D0,-31.30824526D0,125.8739681D0, *-372.3384278D0, *-235.4720434D0,286.7594095D0,21.86305585D0,-27.42344605D0, *-150.4874688D0, *2.669338538D0,1.395023949D0,-.5540427503D0,-56.85224007D0, *3.681827033D0, *-43.48705106D0,5.103131905D0,1.073551279D0,-.6673083508D0, *12.21404266D0, *4.177465543D0,5.799964188D0,-.3977802319D0,-1.044652977D0, *.5703560010D0, *3.536082962D0,-3.222069852D0,9.620648151D0,6.082014949D0, *27.75216226D0, *12.44199571D0,5.122226936D0,6.982039615D0,20.12149582D0, *6.150973118D0, *4.663639687D0,15.73319647D0,2.303504968D0,5.840511214D0, *.8385953499D-01, *.3477844929D0/ 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_04 (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_04 TILT: C WARPING IN Y-Z (DONE BY THE S/R WARPED_04) AND BENDING IN X-Z (DONE BY THIS SUBROUTINE) C IMPLICIT REAL*8 (A-H,O-Z) COMMON /RH0_tsyg04/ 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_04(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_04 (IOPT,PS,X,Y,Z,BX1,BY1,BZ1,BX2,BY2,BZ2) C C CALCULATES GSM COMPONENTS OF THE WARPED_04 FIELD FOR TWO TAIL UNIT MODES. C THE WARPING DEFORMATION IS IMPOSED ON THE UNWARPED_04 FIELD, COMPUTED C BY THE S/R "UNWARPED_04". 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_04 (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_04 (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.45869857D0,57.35899080D0,317.5501869D0, *-2.626756717D0, *-93.38053698D0,-199.6467926D0,-858.8129729D0,34.09192395D0, *845.4214929D0, *-29.07463068D0,47.10678547D0,-128.9797943D0,-781.7512093D0, *6.165038619D0, *167.8905046D0,492.0680410D0,1654.724031D0,-46.77337920D0, *-1635.922669D0, *40.86186772D0,-.1349775602D0,-.9661991179D-01,-.1662302354D0, *.002810467517D0,.2487355077D0,.1025565237D0,-14.41750229D0, *-.8185333989D0, *11.07693629D0,.7569503173D0,-9.655264745D0,112.2446542D0, *777.5948964D0, *-5.745008536D0,-83.03921993D0,-490.2278695D0,-1155.004209D0, *39.08023320D0, *1172.780574D0,-39.44349797D0,-14.07211198D0,-40.41201127D0, *-313.2277343D0, *2.203920979D0,8.232835341D0,197.7065115D0,391.2733948D0, *-18.57424451D0, *-437.2779053D0,23.04976898D0,11.75673963D0,13.60497313D0, *4.691927060D0, *18.20923547D0,27.59044809D0,6.677425469D0,1.398283308D0, *2.839005878D0, *31.24817706D0,24.53577264D0/ DATA A2/-287187.1962D0,4970.499233D0,410490.1952D0, *-1347.839052D0, *-386370.3240D0,3317.983750D0,-143462.3895D0,5706.513767D0, *171176.2904D0, *250.8882750D0,-506570.8891D0,5733.592632D0,397975.5842D0, *9771.762168D0, *-941834.2436D0,7990.975260D0,54313.10318D0,447.5388060D0, *528046.3449D0, *12751.04453D0,-21920.98301D0,-21.05075617D0,31971.07875D0, *3012.641612D0, *-301822.9103D0,-3601.107387D0,1797.577552D0,-6.315855803D0, *142578.8406D0, *13161.93640D0,804184.8410D0,-14168.99698D0,-851926.6360D0, *-1890.885671D0, *972475.6869D0,-8571.862853D0,26432.49197D0,-2554.752298D0, *-482308.3431D0, *-4391.473324D0,105155.9160D0,-1134.622050D0,-74353.53091D0, *-5382.670711D0, *695055.0788D0,-916.3365144D0,-12111.06667D0,67.20923358D0, *-367200.9285D0, *-21414.14421D0,14.75567902D0,20.75638190D0,59.78601609D0, *16.86431444D0, *32.58482365D0,23.69472951D0,17.24977936D0,13.64902647D0, *68.40989058D0, *11.67828167D0/ 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_04(D0SC1,DELTADX1,DELTADY,XSC1,YSC1, &ZSC1,FX1,FY1,FZ1) CALL SHLCAR5X5_04(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_04(D0SC2,DELTADX2,DELTADY,XSC2,YSC2, &ZSC2,FX2,FY2,FZ2) CALL SHLCAR5X5_04(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_04(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 1 DBZ=DBZ+F(I)*(2.D0*AS+X*DASDX+Y*DASDY) 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_04(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 2 I=1,5 RP=1.D0/A(50+I) CYPI=DCOS(Y*RP) SYPI=DSIN(Y*RP) C DO 2 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 2 CONTINUE HX=DHX HY=DHY HZ=DHZ C RETURN END c c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C SUBROUTINE BIRK_TOT_04 (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.84663D0,-15541.95244D0,-23210.09824D0, *-32625.03856D0, *-109894.4551D0,-71415.32808D0,58168.94612D0,55564.87578D0, *-22890.60626D0, *-6056.763968D0,5091.368100D0,239.7001538D0,-13899.49253D0, *4648.016991D0, *6971.310672D0,9699.351891D0,32633.34599D0,21028.48811D0, *-17395.96190D0, *-16461.11037D0,7447.621471D0,2528.844345D0,-1934.094784D0, *-588.3108359D0, *-32588.88216D0,10894.11453D0,16238.25044D0,22925.60557D0, *77251.11274D0, *50375.97787D0,-40763.78048D0,-39088.60660D0,15546.53559D0, *3559.617561D0, *-3187.730438D0,309.1487975D0,88.22153914D0,-243.0721938D0, *-63.63543051D0, *191.1109142D0,69.94451996D0,-187.9539415D0,-49.89923833D0, *104.0902848D0, *-120.2459738D0,253.5572433D0,89.25456949D0,-205.6516252D0, *-44.93654156D0, *124.7026309D0,32.53005523D0,-98.85321751D0,-36.51904756D0, *98.88241690D0, *24.88493459D0,-55.04058524D0,61.14493565D0,-128.4224895D0, *-45.35023460D0, *105.0548704D0,-43.66748755D0,119.3284161D0,31.38442798D0, *-92.87946767D0, *-33.52716686D0,89.98992001D0,25.87341323D0,-48.86305045D0, *59.69362881D0, *-126.5353789D0,-44.39474251D0,101.5196856D0,59.41537992D0, *41.18892281D0, *80.86101200D0,3.066809418D0,7.893523804D0,30.56212082D0, *10.36861082D0, *8.222335945D0,19.97575641D0,2.050148531D0,4.992657093D0, *2.300564232D0, *.2256245602D0,-.05841594319D0/ DATA SH12/210260.4816D0,-1443587.401D0,-1468919.281D0, *281939.2993D0, *-1131124.839D0,729331.7943D0,2573541.307D0,304616.7457D0, *468887.5847D0, *181554.7517D0,-1300722.650D0,-257012.8601D0,645888.8041D0, *-2048126.412D0, *-2529093.041D0,571093.7972D0,-2115508.353D0,1122035.951D0, *4489168.802D0, *75234.22743D0,823905.6909D0,147926.6121D0,-2276322.876D0, *-155528.5992D0, *-858076.2979D0,3474422.388D0,3986279.931D0,-834613.9747D0, *3250625.781D0, *-1818680.377D0,-7040468.986D0,-414359.6073D0,-1295117.666D0, *-346320.6487D0, *3565527.409D0,430091.9496D0,-.1565573462D0,7.377619826D0, *.4115646037D0, *-6.146078880D0,3.808028815D0,-.5232034932D0,1.454841807D0, *-12.32274869D0, *-4.466974237D0,-2.941184626D0,-.6172620658D0,12.64613490D0, *1.494922012D0, *-21.35489898D0,-1.652256960D0,16.81799898D0,-1.404079922D0, *-24.09369677D0, *-10.99900839D0,45.94237820D0,2.248579894D0,31.91234041D0, *7.575026816D0, *-45.80833339D0,-1.507664976D0,14.60016998D0,1.348516288D0, *-11.05980247D0, *-5.402866968D0,31.69094514D0,12.28261196D0,-37.55354174D0, *4.155626879D0, *-33.70159657D0,-8.437907434D0,36.22672602D0,145.0262164D0, *70.73187036D0, *85.51110098D0,21.47490989D0,24.34554406D0,31.34405345D0, *4.655207476D0, *5.747889264D0,7.802304187D0,1.844169801D0,4.867254550D0, *2.941393119D0, *.1379899178D0,.06607020029D0/ DATA SH21/162294.6224D0,503885.1125D0,-27057.67122D0, *-531450.1339D0, *84747.05678D0,-237142.1712D0,84133.61490D0,259530.0402D0, *69196.05160D0, *-189093.5264D0,-19278.55134D0,195724.5034D0,-263082.6367D0, *-818899.6923D0, *43061.10073D0,863506.6932D0,-139707.9428D0,389984.8850D0, *-135167.5555D0, *-426286.9206D0,-109504.0387D0,295258.3531D0,30415.07087D0, *-305502.9405D0, *100785.3400D0,315010.9567D0,-15999.50673D0,-332052.2548D0, *54964.34639D0, *-152808.3750D0,51024.67566D0,166720.0603D0,40389.67945D0, *-106257.7272D0, *-11126.14442D0,109876.2047D0,2.978695024D0,558.6019011D0, *2.685592939D0, *-338.0004730D0,-81.99724090D0,-444.1102659D0,89.44617716D0, *212.0849592D0, *-32.58562625D0,-982.7336105D0,-35.10860935D0,567.8931751D0, *-1.917212423D0, *-260.2023543D0,-1.023821735D0,157.5533477D0,23.00200055D0, *232.0603673D0, *-36.79100036D0,-111.9110936D0,18.05429984D0,447.0481000D0, *15.10187415D0, *-258.7297813D0,-1.032340149D0,-298.6402478D0,-1.676201415D0, *180.5856487D0, *64.52313024D0,209.0160857D0,-53.85574010D0,-98.52164290D0, *14.35891214D0, *536.7666279D0,20.09318806D0,-309.7349530D0,58.54144539D0, *67.45226850D0, *97.92374406D0,4.752449760D0,10.46824379D0,32.91856110D0, *12.05124381D0, *9.962933904D0,15.91258637D0,1.804233877D0,6.578149088D0, *2.515223491D0, *.1930034238D0,-.02261109942D0/ DATA SH22/-131287.8986D0,-631927.6885D0,-318797.4173D0, *616785.8782D0, *-50027.36189D0,863099.9833D0,47680.20240D0,-1053367.944D0, *-501120.3811D0, *-174400.9476D0,222328.6873D0,333551.7374D0,-389338.7841D0, *-1995527.467D0, *-982971.3024D0,1960434.268D0,297239.7137D0,2676525.168D0, *-147113.4775D0, *-3358059.979D0,-2106979.191D0,-462827.1322D0,1017607.960D0, *1039018.475D0, *520266.9296D0,2627427.473D0,1301981.763D0,-2577171.706D0, *-238071.9956D0, *-3539781.111D0,94628.16420D0,4411304.724D0,2598205.733D0, *637504.9351D0, *-1234794.298D0,-1372562.403D0,-2.646186796D0,-31.10055575D0, *2.295799273D0, *19.20203279D0,30.01931202D0,-302.1028550D0,-14.78310655D0, *162.1561899D0, *.4943938056D0,176.8089129D0,-.2444921680D0,-100.6148929D0, *9.172262228D0, *137.4303440D0,-8.451613443D0,-84.20684224D0,-167.3354083D0, *1321.830393D0, *76.89928813D0,-705.7586223D0,18.28186732D0,-770.1665162D0, *-9.084224422D0, *436.3368157D0,-6.374255638D0,-107.2730177D0,6.080451222D0, *65.53843753D0, *143.2872994D0,-1028.009017D0,-64.22739330D0,547.8536586D0, *-20.58928632D0, *597.3893669D0,10.17964133D0,-337.7800252D0,159.3532209D0, *76.34445954D0, *84.74398828D0,12.76722651D0,27.63870691D0,32.69873634D0, *5.145153451D0, *6.310949163D0,6.996159733D0,1.971629939D0,4.436299219D0, *2.904964304D0, *.1486276863D0,.06859991529D0/ XKAPPA=XKAPPA1 ! FORWARDED IN BIRK_1N2_04 X_SC=XKAPPA1-1.1D0 ! FORWARDED IN BIRK_SHL_04 IF (IOPB.EQ.0.OR.IOPB.EQ.1) THEN CALL BIRK_1N2_04 (1,1,PS,X,Y,Z,FX11,FY11,FZ11) ! REGION 1, MODE 1 CALL BIRK_SHL_04 (SH11,PS,X_SC,X,Y,Z,HX11,HY11,HZ11) BX11=FX11+HX11 BY11=FY11+HY11 BZ11=FZ11+HZ11 CALL BIRK_1N2_04 (1,2,PS,X,Y,Z,FX12,FY12,FZ12) ! REGION 1, MODE 2 CALL BIRK_SHL_04 (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_04 X_SC=XKAPPA2-1.0D0 ! FORWARDED IN BIRK_SHL_04 IF (IOPB.EQ.0.OR.IOPB.EQ.2) THEN CALL BIRK_1N2_04 (2,1,PS,X,Y,Z,FX21,FY21,FZ21) ! REGION 2, MODE 1 CALL BIRK_SHL_04 (SH21,PS,X_SC,X,Y,Z,HX21,HY21,HZ21) BX21=FX21+HX21 BY21=FY21+HY21 BZ21=FZ21+HZ21 CALL BIRK_1N2_04 (2,2,PS,X,Y,Z,FX22,FY22,FZ22) ! REGION 2, MODE 2 CALL BIRK_SHL_04 (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_04 (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/.1618068350D0,-.1797957553D0,2.999642482D0, *-.9322708978D0, *-.6811059760D0,.2099057262D0,-8.358815746D0,-14.86033550D0, *.3838362986D0, *-16.30945494D0,4.537022847D0,2.685836007D0,27.97833029D0, *6.330871059D0, *1.876532361D0,18.95619213D0,.9651528100D0,.4217195118D0, *-.08957770020D0, *-1.823555887D0,.7457045438D0,-.5785916524D0,-1.010200918D0, *.01112389357D0, *.09572927448D0,-.3599292276D0,8.713700514D0,.9763932955D0, *3.834602998D0, *2.492118385D0,.7113544659D0/ DATA A12/.7058026940D0,-.2845938535D0,5.715471266D0, *-2.472820880D0, *-.7738802408D0,.3478293930D0,-11.37653694D0,-38.64768867D0, *.6932927651D0, *-212.4017288D0,4.944204937D0,3.071270411D0,33.05882281D0, *7.387533799D0, *2.366769108D0,79.22572682D0,.6154290178D0,.5592050551D0, *-.1796585105D0, *-1.654932210D0,.7309108776D0,-.4926292779D0,-1.130266095D0, *-.009613974555D0, *.1484586169D0,-.2215347198D0,7.883592948D0,.02768251655D0, *2.950280953D0, *1.212634762D0,.5567714182D0/ DATA A21/.1278764024D0,-.2320034273D0,1.805623266D0, *-32.37241440D0, *-.9931490648D0,.3175085630D0,-2.492465814D0,-16.21600096D0, *.2695393416D0, *-6.752691265D0,3.971794901D0,14.54477563D0,41.10158386D0, *7.912889730D0, *1.258297372D0,9.583547721D0,1.014141963D0,.5104134759D0, *-.1790430468D0, *-1.756358428D0,.7561986717D0,-.6775248254D0,-.04014016420D0, *.01446794851D0, *.1200521731D0,-.2203584559D0,4.508963850D0,.8221623576D0, *1.779933730D0, *1.102649543D0,.8867880020D0/ DATA A22/.4036015198D0,-.3302974212D0,2.827730930D0, *-45.44405830D0, *-1.611103927D0,.4927112073D0,-.003258457559D0,-49.59014949D0, *.3796217108D0, *-233.7884098D0,4.312666980D0,18.05051709D0,28.95320323D0, *11.09948019D0, *.7471649558D0,67.10246193D0,.5667096597D0,.6468519751D0, *-.1560665317D0, *-1.460805289D0,.7719653528D0,-.6658988668D0,.2515179349D-05, *.02426021891D0,.1195003324D0,-.2625739255D0,4.377172556D0, *.2421190547D0, *2.503482679D0,1.071587299D0,.7247997430D0/ B=0.5D0 RHO_0=7.0D0 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_04 (A11,XS,Ysc,ZS,BXS,BYAS,BZS) IF (MODE.EQ.2) CALL TWOCONES_04 (A12,XS,Ysc,ZS,BXS,BYAS,BZS) ELSE IF (MODE.EQ.1) CALL TWOCONES_04 (A21,XS,Ysc,ZS,BXS,BYAS,BZS) IF (MODE.EQ.2) CALL TWOCONES_04 (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_04 (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_04 (A,X,Y,Z,BXN,BYN,BZN) CALL ONE_CONE_04 (A,X,-Y,-Z,BXS,BYS,BZS) BX=BXN-BXS BY=BYN+BYS BZ=BZN+BZS RETURN END c C------------------------------------------------------------------------- C SUBROUTINE ONE_CONE_04(A,X,Y,Z,BX,BY,BZ) c c RETURNS FIELD COMPONENTS FOR A DEFORMED_04 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_04(A,R,THETA) THETAS=THETA_S_04(A,R,THETA) PHIS=PHI C C CALCULATE FIELD COMPONENTS AT THE NEW POSITION (ASTERISKED): C CALL FIALCOS_04 (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_04(A,R+DR,THETA)-R_S_04(A,R-DR,THETA))/(2.D0*DR) DRSDT=(R_S_04(A,R,THETA+DT)-R_S_04(A,R,THETA-DT))/(2.D0*DT) DTSDR=(THETA_S_04(A,R+DR,THETA)-THETA_S_04(A,R-DR,THETA))/ &(2.D0*DR) DTSDT=(THETA_S_04(A,R,THETA+DT)-THETA_S_04(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_04(A,R,THETA) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(31) C R_S_04=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_04(A,R,THETA) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(31) c THETA_S_04=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_04(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 2 BPN(M)=-DTT*SSIN(M)/R BTHETA=BTN(N) *800.D0 BPHI =BPN(N) *800.D0 RETURN END C C------------------------------------------------------------------------- C C SUBROUTINE BIRK_SHL_04 (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) 5 GZ=GZ+HZR*A(L) 4 CONTINUE 3 CONTINUE 2 CONTINUE 1 CONTINUE BX=GX BY=GY BZ=GZ RETURN END C C************************************************************************************ C SUBROUTINE FULL_RC_04 (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, * 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.58481D0,-63965.62048D0,66267.93413D0, *135049.7504D0, *-36.56316878D0,124.6614669D0,56.75637955D0,-87.56841077D0, *5848.631425D0, *4981.097722D0,-6233.712207D0,-10986.40188D0,68716.52057D0, *65682.69473D0, *-69673.32198D0,-138829.3568D0,43.45817708D0,-117.9565488D0, *-62.14836263D0, *79.83651604D0,-6211.451069D0,-5151.633113D0,6544.481271D0, *11353.03491D0, *23.72352603D0,-256.4846331D0,25.77629189D0,145.2377187D0, *-4.472639098D0, *-3.554312754D0,2.936973114D0,2.682302576D0,2.728979958D0, *26.43396781D0, *-9.312348296D0,-29.65427726D0,-247.5855336D0,-206.9111326D0, *74.25277664D0, *106.4069993D0,15.45391072D0,16.35943569D0,-5.965177750D0, *-6.079451700D0, *115.6748385D0,-35.27377307D0,-32.28763497D0,-32.53122151D0, *93.74409310D0, *84.25677504D0,-29.23010465D0,-43.79485175D0,-6.434679514D0, *-6.620247951D0, *2.443524317D0,2.266538956D0,-43.82903825D0,6.904117876D0, *12.24289401D0, *17.62014361D0,152.3078796D0,124.5505289D0,-44.58690290D0, *-63.02382410D0, *-8.999368955D0,-9.693774119D0,3.510930306D0,3.770949738D0, *-77.96705716D0, *22.07730961D0,20.46491655D0,18.67728847D0,9.451290614D0, *9.313661792D0, *644.7620970D0,418.2515954D0,7.183754387D0,35.62128817D0, *19.43180682D0, *39.57218411D0,15.69384715D0,7.123215241D0,2.300635346D0, *21.90881131D0, *-.01775839370D0,.3996346710D0/ CALL SRC_PRC_04 (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_04 (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_04 (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_04 (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_04(XTS,YTS,ZTS,BXS,BYS,BZS) IF (IOPR.EQ.0.OR.IOPR.EQ.2) * CALL PRC_SYMM_04(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_04(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_04 (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_04(R,DS,DC)/DS DARDR=(RP*AP_04(RP,DS,DC)-RM*AP_04(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_04(R,SINTP,COSTP)-SINTM*AP_04(R,SINTM,COSTM)) * /(R*SINT)*DRD BT=(RM*AP_04(RM,SINT,COST)-RP*AP_04(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_04(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=.99994999875D0 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_04) 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.03590092383D0+ * XK2S*(0.03742563713D0+XK2S*0.01451196212D0))) +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.03590092383D0+ * XK2S*(0.03742563713D0+XK2S*0.01451196212D0))) +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_04=A1*APHI1+A2*APHI2 IF (PROX) AP_04=AP_04*SINT/SINT1 ! LINEAR INTERPOLATION, IF TOO CLOSE TO THE Z-AXIS C RETURN END c c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ C SUBROUTINE PRC_SYMM_04 (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_04(R,DS,DC)/DS DARDR=(RP*APPRC_04(RP,DS,DC)-RM*APPRC_04(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_04(R,SINTP,COSTP)-SINTM* * APPRC_04(R,SINTM,COSTM)) * /(R*SINT)*DRD BT=(RM*APPRC_04(RM,SINT,COST)-RP*APPRC_04(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_04(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.11202281D0,12.58246758D0,6.560486035D0,1.930711037D0, *3.827208119D0, *.7789990504D0,.3058309043D0,.1817139853D0,.1257532909D0, *3.422509402D0, *.04742939676D0,-4.800458958D0,-.02845643596D0,.2188114228D0, *2.545944574D0, *.00813272793D0,.35868244D0,103.1601001D0,-.00764731187D0, *.1046487459D0, *2.958863546D0,.01172314188D0,.4382872938D0,.01134908150D0, *14.51339943D0, *.2647095287D0,.07091230197D0,.01512963586D0,6.861329631D0, *.1677400816D0, *.04433648846D0,.05553741389D0,.7665599464D0,.7277854652D0/ 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=.99994999875D0 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_04) GAMMA_S=GAMMA*(1.D0+Q0+Q1*(ALPHA-ALPHA4)*DEXP2 ! GAMMA -> GAMMA_ (DEFORMED_04) * +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.03590092383D0+ * XK2S*(0.03742563713D0+XK2S*0.01451196212D0))) +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.03590092383D0+ * XK2S*(0.03742563713D0+XK2S*0.01451196212D0))) +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_04=A1*APHI1+A2*APHI2 IF (PROX) APPRC_04=APPRC_04*SINT/SINT1 ! LINEAR INTERPOLATION, IF TOO CLOSE TO THE Z-AXIS C RETURN END C C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ C C SUBROUTINE PRC_QUAD_04 (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_04(R,SINT,COST) BT=BT_PRC_Q_04(R,SINT,COST) DBRR=(BR_PRC_Q_04(RP,SINT,COST)-BR_PRC_Q_04(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_04(R,SINTP,COSTP)- & BT_PRC_Q_04(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_04(R,ST,CT) BT=BT_PRC_Q_04(R,ST,CT) DBRR=(BR_PRC_Q_04(RP,ST,CT)-BR_PRC_Q_04(RM,ST,CT))/DD DBTT=(BT_PRC_Q_04(R,SINTP,COSTP)- & BT_PRC_Q_04(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_04 (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.2666329D0, ! SO THAT THEY CORRESPOND TO P_0=1 nPa, *32.24527521D0,-6.062894078D0,7.515660734D0,233.7341288D0, *-227.1195714D0, ! RATHER THAN THE ORIGINAL VALUE OF 10 nPa *8.483233889D0,16.80642754D0,-24.63534184D0,9.067120578D0, *-1.052686913D0, ! ASSUMED IN THE BIOT-SAVART INTEGRAL *-12.08384538D0,18.61969572D0,-12.71686069D0,47017.35679D0, *-50646.71204D0, *7746.058231D0,1.531069371D0,2.318824273D0,.1417519429D0, *.6388013110D-02, *5.303934488D0,4.213397467D0,.7955534018D0,.1401142771D0, *.2306094179D-01, *3.462235072D0,2.568743010D0,3.477425908D0,1.922155110D0, *.1485233485D0, *.2319676273D-01,7.830223587D0,8.492933868D0,.1295221828D0, *.01753008801D0, *.01125504083D0,.1811846095D0,.04841237481D0,.01981805097D0, *6.557801891D0, *6.348576071D0,5.744436687D0,.2265212965D0,.1301957209D0, *.5654023158D0/ SINT2=SINT**2 COST2=COST**2 SC=SINT*COST ALPHA=SINT2/R GAMMA=COST/R**2 CALL FFS_04(ALPHA,AL1,DAL1,F,FA,FS) D1=SC*F**XK1/((R/B1)**BE1+1.D0) D2=D1*COST2 CALL FFS_04(ALPHA,AL2,DAL2,F,FA,FS) D3=SC*FS**XK2/((R/B2)**BE2+1.D0) D4=D3*COST2 CALL FFS_04(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_04(ALPHA,AL6,DAL6,F,FA,FS) D18=SC*FS/(1.D0+((R-1.2D0)/DRM)**2) BR_PRC_Q_04=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_04 (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.74640393D0,-7.516393516D0, ! SO THAT THEY CORRESPOND TO P_0=1 nPa, *-5.476233865D0,3.212704645D0,-59.10926169D0,46.62198189D0, *-.01644280062D0, ! RATHER THAN THE ORIGINAL VALUE OF 10 nPa *.1234229112D0,-.08579198697D0,.01321366966D0,.8970494003D0, *9.136186247D0, ! ASSUMED IN THE BIOT-SAVART INTEGRAL *-38.19301215D0,21.73775846D0,-410.0783424D0,-69.90832690D0, *-848.8543440D0, *1.243288286D0,.2071721360D0,.05030555417D0,7.471332374D0, *3.180533613D0, *1.376743507D0,.1568504222D0,.02092910682D0,1.985148197D0, *.3157139940D0, *1.056309517D0,.1701395257D0,.1019870070D0,6.293740981D0, *5.671824276D0, *.1280772299D0,.02189060799D0,.01040696080D0,.1648265607D0, *.04701592613D0, *.01526400086D0,12.88384229D0,3.361775101D0,23.44173897D0/ SINT2=SINT**2 COST2=COST**2 SC=SINT*COST ALPHA=SINT2/R GAMMA=COST/R**2 CALL FFS_04(ALPHA,AL1,DAL1,F,FA,FS) D1=F**XK1/((R/B1)**BE1+1.D0) D2=D1*COST2 CALL FFS_04(ALPHA,AL2,DAL2,F,FA,FS) D3=FA**XK2/R**BE2 D4=D3*COST2 CALL FFS_04(ALPHA,AL3,DAL3,F,FA,FS) D5=FS**XK3*ALPHA**XK4/((R/B3)**BE3+1.D0) D6=D5*COST2 CALL FFS_04(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_04=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_04(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_04 (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) 5 GZ=GZ+HZR*A(L) 4 CONTINUE 3 CONTINUE 2 CONTINUE 1 CONTINUE BX=GX BY=GY BZ=GZ RETURN END C c=========================================================================== c SUBROUTINE DIPOLE_04 (PS,X,Y,Z,BX,BY,BZ) C C A DOUBLE PRECISION ROUTINE C C CALCULATES GSM COMPONENTS OF A GEODIPOLE_04 FIELD WITH THE DIPOLE_04 MOMENT C CORRESPONDING TO THE EPOCH OF 2000. C C----INPUT PARAMETERS: C PS - GEODIPOLE_04 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(((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((