c-----this program transforms the MF from the GSE r.f.to GSM c-----with a use of the GEOPACK software integer UT_mf(8) real Bgse(3),Bgsm(3) c character*12 inputd,output open(1,file=' ') open(3,file=' ') do while (.not.EOF(1)) c-------read MF in the GSE system------------------ read(1,678) UT_mf,Bgse 678 format(1x,i4.4,5i2.2,i4.4,i3.3,3(1x,f6.1)) call trnsfrm(UT_mf,Bgse,Bgsm) write(3,678) UT_mf,Bgsm END DO stop end c-------SUB--------------------------------- SUBROUTINE TRNSFRM(UT_s,BinE,BoutM) COMMON ST0,CT0,SL0,CL0,CTCL,STCL,CTSL,STSL,SFI,CFI,SPS,CPS, *SHI,CHI,HI,PSI,XMUT,A(3,3),AB,K,IY,BA(8) integer UT_s(8) real BinE(3),BoutM(3) Bxgse=BinE(1) Bygse=BinE(2) Bzgse=BinE(3) call RECALC(UT_s(1),UT_s(8), * UT_s(4),UT_s(5),UT_s(6)) call GSMGSE(Bxgsm,Bygsm,Bzgsm,Bxgse,Bygse,Bzgse,-1) BoutM(1)=Bxgsm BoutM(2)=Bygsm BoutM(3)=Bzgsm return end c++++++++++++++++++++++++++++++++++++++++++++++++ c=======MAGNETIC FIELD MODEL===================== c++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE SUN(IYR,IDAY,IHOUR,MIN,ISEC,GST,SLONG,SRASN,SDEC) C C CALCULATES FOUR QUANTITIES NECESSARY FOR COORDINATE TRANSFORMATIONS C DEPENDENT ON SUN POSITION (AND, HENCE, ON UNIVERSAL TIME AND SEASON) C------- INPUT PARAMETERS: C IYR,IDAY,IHOUR,MIN,ISEC - YEAR, DAY, AND UNIVERSAL TIME C (IDAY=1 CORRESPONDS TO JANUARY 1). C------- OUTPUT PARAMETERS: C GST - GREENWICH MEAN SIDEREAL TIME, SLONG - LONGITUDE ALONG ECLIPTIC C SRASN - RIGHT ASCENSION, SDEC - DECLINATION OF THE SUN (RADIANS) C THIS SUBROUTINE HAS BEEN COMPILED FROM: RUSSELL C.T., COSM.ELECTRO- C DYN., 1971, V.2,PP.184-196. C C C AUTHOR: CHRISTOPHER T.RUSSELL C INSTITUTE OF GEOPHYSICS AND PLANETARY PHYSICS C UNIVERSITY OF CALIFORNIA, LOS ANGELES C LOS ANGELES, CALIFORNIA 90024 C U.S.A. C DOUBLE PRECISION DJ,FDAY DATA RAD/57.29578/ IF(IYR.LT.1901.OR.IYR.GT.2099) RETURN FDAY=DBLE((IHOUR*3600.+MIN*60.+ISEC)/86400.) DJ=365*(IYR-1900)+(IYR-1901)/4+IDAY-0.5D0+FDAY T=DJ/36525. VL=DMOD(279.696678+0.9856473354*DJ,360.D0) GST=DMOD(279.690983+.9856473354*DJ+360.*FDAY+180.,360.D0)/RAD G=DMOD(358.475845+0.985600267*DJ,360.D0)/RAD SLONG=(VL+(1.91946-0.004789*T)*SIN(G)+0.020094*SIN(2.*G))/RAD IF(SLONG.GT.6.2831853) SLONG=SLONG-6.2831853 IF (SLONG.LT.0.) SLONG=SLONG+6.2831853 OBLIQ=(23.45229-0.0130125*T)/RAD SOB=SIN(OBLIQ) SIND=SOB*SIN(SLONG-9.924E-5) COSD=SQRT(1.-SIND**2) SC=SIND/COSD SDEC=ATAN(SC) SRASN=3.141592654-ATAN2(COS(OBLIQ)/SOB*SC,-COS(SLONG)/COSD) RETURN END SUBROUTINE SPHCAR(R,TETA,PHI,X,Y,Z,J) C C CONVERTS SPHERICAL COORDS INTO CARTESIAN ONES AND VICA VERSA C (TETA AND PHI IN RADIANS). C C J>0 J<0 C-----INPUT: J,R,TETA,PHI J,X,Y,Z C----OUTPUT: X,Y,Z R,TETA,PHI C C C AUTHOR: NIKOLAI A. TSYGANENKO C INSTITUTE OF PHYSICS C LENINGRAD STATE UNIVERSITY C STARY PETERGOF 198904 C LENINGRAD C USSR C IF(J.GT.0) GOTO 3 SQ=X**2+Y**2 R=SQRT(SQ+Z**2) IF (SQ.NE.0.) GOTO 2 PHI=0. IF (Z.LT.0.) GOTO 1 TETA=0. RETURN 1 TETA=3.141592654 RETURN 2 SQ=SQRT(SQ) PHI=ATAN2(Y,X) TETA=ATAN2(SQ,Z) IF (PHI.LT.0.) PHI=PHI+6.28318531 RETURN 3 SQ=R*SIN(TETA) X=SQ*COS(PHI) Y=SQ*SIN(PHI) Z=R*COS(TETA) RETURN END SUBROUTINE BSPCAR(TETA,PHI,BR,BTET,BPHI,BX,BY,BZ) C CALCULATES CARTESIAN FIELD COMPONENTS FROM SPHERICAL ONES C-----INPUT: TETA,PHI - SPHERICAL ANGLES OF THE POINT IN RADIANS C BR,BTET,BPHI - SPHERICAL COMPONENTS OF THE FIELD C-----OUTPUT: BX,BY,BZ - CARTESIAN COMPONENTS OF THE FIELD C C C AUTHOR: NIKOLAI A. TSYGANENKO C INSTITUTE OF PHYSICS C LENINGRAD STATE UNIVERSITY C STARY PETERGOF 198904 C LENINGRAD C USSR C S=SIN(TETA) C=COS(TETA) SF=SIN(PHI) CF=COS(PHI) BE=BR*S+BTET*C BX=BE*CF-BPHI*SF BY=BE*SF+BPHI*CF BZ=BR*C-BTET*S RETURN END SUBROUTINE RECALC(IYR,IDAY,IHOUR,MIN,ISEC) C C CALCULATES GEODIPOLE AXIS ORIENTATION FOR A GIVEN YEAR, DAY AND UT; C PREPARES MATRICES FOR THE FIVE TRANSFORMATION SUBROUTINES: geomag,geogsm, c magsm,smgsm,gsmgse C ALL THE PARAMETERS ARE INPUT ONES: C IYR - YEAR (FROM 1965 UP TO 1990), C IDAY - DAY NUMBER (JANUARY 1 IS IDAY=1), C IHOUR,MIN,ISEC - UT. IF IHOUR>24 THEN ONLY GEOMAG MATRIX IS PRERARED. C C C AUTHOR: NIKOLAI A. TSYGANENKO C INSTITUTE OF PHYSICS C LENINGRAD STATE UNIVERSITY C STARY PETERGOF 198904 C LENINGRAD C USSR C COMMON ST0,CT0,SL0,CL0,CTCL,STCL,CTSL,STSL,SFI,CFI,SPS,CPS, *SHI,CHI,HI,PSI,XMUT,A(3,3),AB,K,IY,BA(8) DIMENSION X(3,3) DATA IYE,IPR/0,0/ IF(IYR.EQ.IYE) GOTO 5 IY=IYR IF (IY.LT.1965) IY=1965 C IF (IY.GT.1990) IY=1990 IF (IY.NE.IYR.AND.IPR.EQ.0) write(*,10)IY,IYR IF (IY.NE.IYR) IPR=1 IYE=IY IF (IY.LT.1970) GOTO 1 IF (IY.LT.1975) GOTO 2 IF (IY.LT.1980) GOTO 3 IF (IY.LT.1985) GOTO 30 DT=FLOAT(IY)-1985. G10=29877.-23.2*DT G11=-1903.+10.0*DT H11=5497.-24.5*DT GOTO 4 1 F2=(IY-1965)/5. F1=1.-F2 G10=30334.*F1+30220.*F2 G11=-2119.*F1-2068.*F2 H11=5776.*F1+5737.*F2 GOTO 4 2 F2=(IY-1970)/5. F1=1.-F2 G10=30220.*F1+30100.*F2 G11=-2068.*F1-2013.*F2 H11=5737.*F1+5675.*F2 GOTO 4 3 F2=(IY-1975)/5. F1=1.-F2 G10=30100.*F1+29992.*F2 G11=-2013.*F1-1956.*F2 H11=5675.*F1+5604.*F2 30 F2=(IY-1980)/5. F1=1.-F2 G10=29992.*F1+29877.*F2 G11=-1956.*F1-1903.*F2 H11=5604.*F1+5497.*F2 4 XL0=ATAN(H11/G11) SQ=G11**2+H11**2 SQQ=SQRT(SQ) SQ=SQRT(SQ+G10**2) ST0=SQQ/SQ CT0=G10/SQ SL0=-H11/SQQ CL0=-G11/SQQ CTCL= CT0*CL0 CTSL=CT0*SL0 STCL=ST0*CL0 STSL=ST0*SL0 5 IF(IHOUR.GT.24) GOTO 8 CALL SUN(IYR,IDAY,IHOUR,MIN,ISEC,GST,SLONG,SRASN,SDEC) AL=GST-SRASN SAL=SIN(AL) CAL=COS(AL) DEL=AL+XL0 SDL=SIN(DEL) CDL=COS(DEL) CSL=COS(SLONG) SSD=SIN(SDEC) CSD=COS(SDEC) STCD=ST0*CDL SPS=CT0*SSD+STCD*CSD Y1=ST0*SDL Z1=CT0*CSD-STCD*SSD TB=.4336072*CSL CB=1./SQRT(1.+TB**2) SB=TB*CB THI=(Y1*CB+Z1*SB)/(Y1*SB-Z1*CB) CHI=1./SQRT(1.+THI**2) SHI=CHI*THI CPS=SQRT(1.-SPS**2) HI=ATAN(THI) PSI=ATAN(SPS/CPS) X1=CSD*CAL Y1=-CSD*SAL Z1=(X1*CL0+Y1*SL0)*CT0-SSD*ST0 Y1=Y1*CL0-X1*SL0 XMUT=12.-3.8197*ATAN2(Y1,Z1) FI=.2618*(XMUT-12.) SFI=SIN(FI) CFI=COS(FI) DO 6 I=1,3 DO 6 L=1,3 X(I,L)=0. IF (I.EQ.L) X(I,L)=1. 6 CONTINUE DO 7 I=1,3 Y1=X(I,1)*CTCL+X(I,2)*CTSL-X(I,3)*ST0 Y2=X(I,2)*CL0-X(I,1)*SL0 Y3=X(I,1)*STCL+X(I,2)*STSL+X(I,3)*CT0 CALL MAGSM(Y1,Y2,Y3,Z1,Z2,Z3,1) 7 CALL SMGSM(Z1,Z2,Z3,A(1,I),A(2,I),A(3,I),1) 8 CONTINUE 10 FORMAT(//1X,'**** YEAR IS OUT OF INTERVAL 1965-1990: IY =',I5, *', CALCULATIONS WILL BE DONE FOR IYR =',I5,' ****'//) RETURN END SUBROUTINE GEOMAG(XGEO,YGEO,ZGEO,XMAG,YMAG,ZMAG,J,IYR) C C CONVERTS GEOGRAPHIC (GEO) TO DIPOLE (MAG) COORDINATES OR VICA VERSA. C IYR IS YEAR NUMBER (FOUR DIGITS). C C J>0 J<0 C-----INPUT: J,XGEO,YGEO,ZGEO,IYR J,XMAG,YMAG,ZMAG,IYR C-----OUTPUT: XMAG,YMAG,ZMAG XGEO,YGEO,ZGEO C C C AUTHOR: NIKOLAI A. TSYGANENKO C INSTITUTE OF PHYSICS C LENINGRAD STATE UNIVERSITY C STARY PETERGOF 198904 C LENINGRAD C USSR C COMMON ST0,CT0,SL0,CL0,CTCL,STCL,CTSL,STSL,AB(19),K,IY,BB(8) DATA II/1/ IF(IYR.EQ.II) GOTO 1 II=IYR CALL RECALC(II,0,25,0,0) 1 CONTINUE IF(J.LT.0) GOTO 2 XMAG=XGEO*CTCL+YGEO*CTSL-ZGEO*ST0 YMAG=YGEO*CL0-XGEO*SL0 ZMAG=XGEO*STCL+YGEO*STSL+ZGEO*CT0 RETURN 2 XGEO=XMAG*CTCL-YMAG*SL0+ZMAG*STCL YGEO=XMAG*CTSL+YMAG*CL0+ZMAG*STSL ZGEO=ZMAG*CT0-XMAG*ST0 RETURN END SUBROUTINE MAGSM(XMAG,YMAG,ZMAG,XSM,YSM,ZSM,J) C C CONVERTS DIPOLE (MAG) TO SOLAR MAGNETIC (SM) COORDINATES OR VICA VERSA C C J>0 J<0 C-----INPUT: J,XMAG,YMAG,ZMAG J,XSM,YSM,ZSM C----OUTPUT: XSM,YSM,ZSM XMAG,YMAG,ZMAG C ATTENTION: SUBROUTINE RECALC MUST BE CALLED BEFORE MAGSM IN TWO CASES: C /A/ BEFORE THE FIRST USE OF MAGSM C /B/ IF THE CURRENT VALUES OF IYEAR,IDAY,IHOUR,MIN,ISEC ARE DIFFERENT C FROM THOSE IN THE PRECEDING CALL OF MAGSM C C AUTHOR: NIKOLAI A. TSYGANENKO C INSTITUTE OF PHYSICS C LENINGRAD STATE UNIVERSITY C STARY PETERGOF 198904 C LENINGRAD C USSR C COMMON A(8),SFI,CFI,B(7),AB(10),K,IY,BA(8) IF (J.LT.0) GOTO 1 XSM=XMAG*CFI-YMAG*SFI YSM=XMAG*SFI+YMAG*CFI ZSM=ZMAG RETURN 1 XMAG=XSM*CFI+YSM*SFI YMAG=YSM*CFI-XSM*SFI ZMAG=ZSM RETURN END SUBROUTINE GSMGSE(XGSM,YGSM,ZGSM,XGSE,YGSE,ZGSE,J) C C CONVERTS SOLAR MAGNETOSPHERIC (GSM) TO SOLAR ECLIPTICAL (GSE) COORDS C OR VICA VERSA. C J>0 J<0 C-----INPUT: J,XGSM,YGSM,ZGSM J,XGSE,YGSE,ZGSE C----OUTPUT: XGSE,YGSE,ZGSE XGSM,YGSM,ZGSM C ATTENTION: SUBROUTINE RECALC MUST BE CALLED BEFORE GSMGSE IN TWO CASES: C /A/ BEFORE THE FIRST CALL OF GSMGSE C /B/ IF THE CURRENT VALUES OF IYEAR,IDAY,IHOUR,MIN,ISEC ARE DIFFERENT C FROM THOSE IN THE PRECEDING CALL OF GSMGSE C C C AUTHOR: NIKOLAI A. TSYGANENKO C INSTITUTE OF PHYSICS C LENINGRAD STATE UNIVERSITY C STARY PETERGOF 198904 C LENINGRAD C USSR C COMMON A(12),SHI,CHI,AB(13),K,IY,BA(8) IF(J.LT.0) GOTO 1 XGSE=XGSM YGSE=YGSM*CHI-ZGSM*SHI ZGSE=YGSM*SHI+ZGSM*CHI RETURN 1 XGSM=XGSE YGSM=YGSE*CHI+ZGSE*SHI ZGSM=ZGSE*CHI-YGSE*SHI RETURN END SUBROUTINE SMGSM(XSM,YSM,ZSM,XGSM,YGSM,ZGSM,J) C C CONVERTS SOLAR MAGNETIC (SM) TO SOLAR MAGNETOSPHERIC (GSM) COORDINATES C OR VICA VERSA. C J>0 J<0 C-----INPUT: J,XSM,YSM,ZSM J,XGSM,YGSM,ZGSM C----OUTPUT: XGSM,YGSM,ZGSM XSM,YSM,ZSM C C ATTENTION: SUBROUTINE RECALC MUST BE CALLED BEFORE SMGSM IN TWO CASES: C /A/ BEFORE THE FIRST USE OF SMGSM C /B/ IF THE CURRENT VALUES OF IYEAR,IDAY,IHOUR,MIN,ISEC ARE DIFFERENT C FROM THOSE IN THE PRECEDING CALL OF SMGSM C C C AUTHOR: NIKOLAI A. TSYGANENKO C INSTITUTE OF PHYSICS C LENINGRAD UNIVERSITY C STARY PETERGOF 198904 C LENINGRAD C U.S.S.R. C COMMON A(10),SPS,CPS,B(15),K,IY,AB(8) IF (J.LT.0) GOTO 1 XGSM=XSM*CPS+ZSM*SPS YGSM=YSM ZGSM=ZSM*CPS-XSM*SPS RETURN 1 XSM=XGSM*CPS-ZGSM*SPS YSM=YGSM ZSM=XGSM*SPS+ZGSM*CPS RETURN END SUBROUTINE GEOGSM(XGEO,YGEO,ZGEO,XGSM,YGSM,ZGSM,J) C C CONVERTS GEOGRAPHIC TO SOLAR MAGNETOSPHERIC COORDINATES OR VICA VERSA. C C J>0 J<0 C----- INPUT: J,XGEO,YGEO,ZGEO J,XGSM,YGSM,ZGSM C---- OUTPUT: XGSM,YGSM,ZGSM XGEO,YGEO,ZGEO C ATTENTION: SUBROUTINE RECALC MUST BE CALLED BEFORE GEOGSM IN TWO CASES: C /A/ BEFORE THE FIRST USE OF GEOGSM C /B/ IF THE CURRENT VALUES OF IYEAR,IDAY,IHOUR,MIN,ISEC ARE DIFFERENT C FROM THOSE IN THE PREVIOUS CALL OF THIS SUBROUTINE C C C AUTHOR: NIKOLAI A. TSYGANENKO C INSTITUTE OF PHYSICS C LENINGRAD STATE UNIVERSITY C STARY PETERGOF 198904 C LENINGRAD C USSR C COMMON AA(17),A11,A21,A31,A12,A22,A32,A13,A23,A33,D,K,IY,B(8) IF (J.LT.0) GOTO 1 XGSM=A11*XGEO+A12*YGEO+A13*ZGEO YGSM=A21*XGEO+A22*YGEO+A23*ZGEO ZGSM=A31*XGEO+A32*YGEO+A33*ZGEO RETURN 1 XGEO=A11*XGSM+A21*YGSM+A31*ZGSM YGEO=A12*XGSM+A22*YGSM+A32*ZGSM ZGEO=A13*XGSM+A23*YGSM+A33*ZGSM RETURN END