SUBROUTINE GEODIP(IYR,SLA,SLO,DLA,DLO,J) C CALL GEODIP(IYEAR,LATI,YLONGI,MLAT,MLONG,JMAG) C Calculates dipole geomagnetic coordinates from geocentric coordinates C or vice versa. C J=0 J=1 C INPUT: J,SLA,SLO J,DLA,DLO C OUTPUT: DLA,DLO SLA,SLO C Last revision: November 2005 (Vladimir Papitashvili) C The code is modifed from GEOCOR written by V.Popov and V.Papitashvili C in mid-1980s. COMMON /CONST/UMR,PI C Earth's radius (km) RE = 6371.2 C The radius of the sphere to compute the coordinates (in Re) C RH = (RE + HI)/RE R = 1. if(j.gt.0) goto 1234 COL = (90.- SLA)*UMR RLO = SLO*UMR CALL SPHCAR(R,COL,RLO,X,Y,Z,1) CALL GEOMAG(X,Y,Z,XM,YM,ZM,1,IYR) CALL SPHCAR(RM,TH,PF,XM,YM,ZM,-1) SZM = ZM DLO = PF/UMR DCO = TH/UMR DLA = 90.- DCO RETURN 1234 continue COL = (90.- DLA)*UMR RLO = DLO*UMR CALL SPHCAR(R,COL,RLO,XM,YM,ZM,1) CALL GEOMAG(X,Y,Z,XM,YM,ZM,-1,IYR) CALL SPHCAR(RM,TH,PF,X,Y,Z,-1) SZM = ZM SLO = PF/UMR SCO = TH/UMR SLA = 90.- SCO 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 J>0 J<0 C-----INPUT: J,R,TETA,PHI J,X,Y,Z C----OUTPUT: X,Y,Z R,TETA,PHI C AUTHOR: NIKOLAI A. TSYGANENKO, INSTITUTE OF PHYSICS, ST.-PETERSBURG C STATE UNIVERSITY, STARY PETERGOF 198904, ST.-PETERSBURG, RUSSIA C (now the NASA Goddard Space Fligth Center, Greenbelt, Maryland) C ********************************************************************* IMPLICIT NONE REAL R,TETA,PHI,X,Y,Z,SQ INTEGER J 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 C C SUBROUTINE GEOMAG(XGEO,YGEO,ZGEO,XMAG,YMAG,ZMAG,J,IYR) C ********************************************************************* C CONVERTS GEOCENTRIC (GEO) TO DIPOLE (MAG) COORDINATES OR VICA VERSA. C IYR IS YEAR NUMBER (FOUR DIGITS). 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 AUTHOR: NIKOLAI A. TSYGANENKO, INSTITUTE OF PHYSICS, ST.-PETERSBURG C STATE UNIVERSITY, STARY PETERGOF 198904, ST.-PETERSBURG, RUSSIA C (now the NASA Goddard Space Fligth Center, Greenbelt, Maryland) C ********************************************************************* IMPLICIT NONE REAL XGEO,YGEO,ZGEO,XMAG,YMAG,ZMAG,ST0,CT0,SL0,CL0,CTCL, * STCL,CTSL,STSL,AB(19),BB(8) INTEGER J,IYR,K,IY,II COMMON/C1/ ST0,CT0,SL0,CL0,CTCL,STCL,CTSL,STSL,AB,K,IY,BB 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 RECALC(IYR,IDAY,IHOUR,MIN,ISEC) C ********************************************************************* C If only IYR is given then CALL RECALC(IYR,0,25,0,0) C THIS IS A MODIFIED VERSION OF THE SUBROUTINE RECOMP WRITTEN BY C N. A. TSYGANENKO. SINCE I WANT TO USE IT IN PLACE OF SUBROUTINE C RECALC, I HAVE RENAMED THIS ROUTINE RECALC AND ELIMINATED THE C ORIGINAL RECALC FROM THIS VERSION OF THE PACKAGE. C THIS WAY ALL ORIGINAL CALLS TO RECALC WILL CONTINUE TO WORK WITHOUT C HAVING TO CHANGE THEM TO CALLS TO RECOMP. C C AN ALTERNATIVE VERSION OF THE SUBROUTINE RECALC FROM THE GEOPACK C PACKAGE BASED ON A DIFFERENT APPROACH TO DERIVATION OF ROTATION C MATRIX ELEMENTS C C THIS SUBROUTINE WORKS BY 20% FASTER THAN RECALC AND IS EASIER TO C UNDERSTAND C ##################################################### C # WRITTEN BY N.A. TSYGANENKO ON DECEMBER 1, 1991 # C ##################################################### C Modified by Mauricio Peredo, Hughes STX at NASA/GSFC Code 695, C September 1992 C C Modified to accept years up to year 2000 and updated IGRF coeficients C from 1945 (updated by V. Papitashvili, February 1995) C C Modified to accept years up to 2005 (V. Papitashvili, January 2001) C C Modified to accept years from 1900 through 2010 using the DGRF & C IGRF-10 coeficients (updated by V. Papitashvili, November 2005) C C Modified to accept years up to 2015 (V. Papitashvili, January 2011) C C Modified to accept years up to 2020 (D. Bilitza, October 2015) C C OTHER SUBROUTINES CALLED BY THIS ONE: SUN C C IYR = YEAR NUMBER (FOUR DIGITS) C IDAY = DAY OF YEAR (DAY 1 = JAN 1) C IHOUR = HOUR OF DAY (00 TO 23) C MIN = MINUTE OF HOUR (00 TO 59) C ISEC = SECONDS OF DAY(00 TO 59) C ********************************************************************* IMPLICIT NONE REAL ST0,CT0,SL0,CL0,CTCL,STCL,CTSL,STSL,SFI,CFI,SPS,CPS, 1 SHI,CHI,HI,PSI,XMUT,A11,A21,A31,A12,A22,A32,A13,A23, 2 A33,DS3,F2,F1,G10,G11,H11,DT,SQ,SQQ,SQR,S1,S2, 3 S3,CGST,SGST,DIP1,DIP2,DIP3,Y1,Y2,Y3,Y,Z1,Z2,Z3,DJ, 4 T,OBLIQ,DZ1,DZ2,DZ3,DY1,DY2,DY3,EXMAGX,EXMAGY,EXMAGZ, 5 EYMAGX,EYMAGY,GST,SLONG,SRASN,SDEC,BA(8),DECARG INTEGER IYR,IDAY,IHOUR,MIN,ISEC,K,IY,IDE,IYE,konsol logical mess COMMON/C1/ ST0,CT0,SL0,CL0,CTCL,STCL,CTSL,STSL,SFI,CFI,SPS,CPS, * SHI,CHI,HI,PSI,XMUT,A11,A21,A31,A12,A22,A32,A13,A23,A33,DS3, * K,IY,BA common/iounit/konsol,mess DATA IYE,IDE/2*0/ IF (IYR.EQ.IYE.AND.IDAY.EQ.IDE) GOTO 5 C IYE AND IDE ARE THE CURRENT VALUES OF YEAR AND DAY NUMBER IY=IYR IDE=IDAY IF(IY.LT.1900) IY=1900 c IF(IY.GT.2015) IY=2015 c IF(IY.GT.2020) IY=2020 IF(IY.GT.2025) IY=2025 C IF IYR IS OUTSIDE THE TIME INTERVAL COVERED BY THE CURRENT IGRF, THE C SUBROUTINE GIVES A WARNING (BUT DOES NOT REPEAT IT AT THE NEXT CALLS) IF(IY.NE.IYR.AND.mess) write(konsol,10) IYR,IY IYE=IY C LINEAR INTERPOLATION OF THE GEODIPOLE MOMENT COMPONENTS BETWEEN THE C VALUES FOR THE NEAREST EPOCHS: IF (IY.LT.1905) THEN !1900-1905 F2=(FLOAT(IY)+FLOAT(IDAY)/365.-1900.)/5. F1=1.D0-F2 G10=31543.*F1+31464.*F2 G11=-2298.*F1-2298.*F2 H11= 5922.*F1+5909.*F2 ELSEIF (IY.LT.1910) THEN !1905-1910 F2=(FLOAT(IY)+FLOAT(IDAY)/365.-1905.)/5. F1=1.D0-F2 G10=31464.*F1+31354.*F2 G11=-2298.*F1-2297.*F2 H11= 5909.*F1+5898.*F2 ELSEIF (IY.LT.1915) THEN !1910-1915 F2=(FLOAT(IY)+FLOAT(IDAY)/365.-1910.)/5. F1=1.D0-F2 G10=31354.*F1+31212.*F2 G11=-2297.*F1-2306.*F2 H11= 5898.*F1+5875.*F2 ELSEIF (IY.LT.1920) THEN !1915-1920 F2=(FLOAT(IY)+FLOAT(IDAY)/365.-1915.)/5. F1=1.D0-F2 G10=31212.*F1+31060.*F2 G11=-2306.*F1-2317.*F2 H11= 5875.*F1+5845.*F2 ELSEIF (IY.LT.1925) THEN !1920-1925 F2=(FLOAT(IY)+FLOAT(IDAY)/365.-1920.)/5. F1=1.D0-F2 G10=31060.*F1+30926.*F2 G11=-2317.*F1-2318.*F2 H11= 5845.*F1+5817.*F2 ELSEIF (IY.LT.1930) THEN !1925-1930 F2=(FLOAT(IY)+FLOAT(IDAY)/365.-1925.)/5. F1=1.D0-F2 G10=30926.*F1+30805.*F2 G11=-2318.*F1-2316.*F2 H11= 5817.*F1+5808.*F2 ELSEIF (IY.LT.1935) THEN !1930-1935 F2=(FLOAT(IY)+FLOAT(IDAY)/365.-1930.)/5. F1=1.D0-F2 G10=30805.*F1+30715.*F2 G11=-2316.*F1-2306.*F2 H11= 5808.*F1+5812.*F2 ELSEIF (IY.LT.1940) THEN !1935-1940 F2=(FLOAT(IY)+FLOAT(IDAY)/365.-1935.)/5. F1=1.D0-F2 G10=30715.*F1+30654.*F2 G11=-2306.*F1-2292.*F2 H11= 5812.*F1+5821.*F2 ELSEIF (IY.LT.1945) THEN !1940-1945 F2=(FLOAT(IY)+FLOAT(IDAY)/365.-1940.)/5. F1=1.D0-F2 G10=30654.*F1+30594.*F2 G11=-2292.*F1-2285.*F2 H11= 5821.*F1+5810.*F2 ELSEIF (IY.LT.1950) THEN !1945-1950 F2=(FLOAT(IY)+FLOAT(IDAY)/365.-1945.)/5. F1=1.D0-F2 G10=30594.*F1+30554.*F2 G11=-2285.*F1-2250.*F2 H11= 5810.*F1+5815.*F2 ELSEIF (IY.LT.1955) THEN !1950-1955 F2=(FLOAT(IY)+FLOAT(IDAY)/365.-1950.)/5. F1=1.D0-F2 G10=30554.*F1+30500.*F2 G11=-2250.*F1-2215.*F2 H11= 5815.*F1+5820.*F2 ELSEIF (IY.LT.1960) THEN !1955-1960 F2=(FLOAT(IY)+FLOAT(IDAY)/365.-1955.)/5. F1=1.D0-F2 G10=30500.*F1+30421.*F2 G11=-2215.*F1-2169.*F2 H11= 5820.*F1+5791.*F2 ELSEIF (IY.LT.1965) THEN !1960-1965 F2=(FLOAT(IY)+FLOAT(IDAY)/365.-1960.)/5. F1=1.D0-F2 G10=30421.*F1+30334.*F2 G11=-2169.*F1-2119.*F2 H11= 5791.*F1+5776.*F2 ELSEIF (IY.LT.1970) THEN !1965-1970 F2=(FLOAT(IY)+FLOAT(IDAY)/365.-1965.)/5. F1=1.D0-F2 G10=30334.*F1+30220.*F2 G11=-2119.*F1-2068.*F2 H11= 5776.*F1+5737.*F2 ELSEIF (IY.LT.1975) THEN !1970-1975 F2=(FLOAT(IY)+FLOAT(IDAY)/365.-1970.)/5. F1=1.D0-F2 G10=30220.*F1+30100.*F2 G11=-2068.*F1-2013.*F2 H11= 5737.*F1+5675.*F2 ELSEIF (IY.LT.1980) THEN !1975-1980 F2=(DFLOAT(IY)+DFLOAT(IDAY)/365.-1975.)/5. F1=1.D0-F2 G10=30100.*F1+29992.*F2 G11=-2013.*F1-1956.*F2 H11= 5675.*F1+5604.*F2 ELSEIF (IY.LT.1985) THEN !1980-1985 F2=(FLOAT(IY)+FLOAT(IDAY)/365.-1980.)/5. F1=1.D0-F2 G10=29992.*F1+29873.*F2 G11=-1956.*F1-1905.*F2 H11= 5604.*F1+5500.*F2 ELSEIF (IY.LT.1990) THEN !1985-1990 F2=(FLOAT(IY)+FLOAT(IDAY)/365.-1985.)/5. F1=1.D0-F2 G10=29873.*F1+29775.*F2 G11=-1905.*F1-1848.*F2 H11= 5500.*F1+5406.*F2 ELSEIF (IY.LT.1995) THEN !1990-1995 F2=(FLOAT(IY)+FLOAT(IDAY)/365.-1990.)/5. F1=1.D0-F2 G10=29775.*F1+29692.*F2 G11=-1848.*F1-1784.*F2 H11= 5406.*F1+5306.*F2 ELSEIF (IY.LT.2000) THEN !1995-2000 F2=(FLOAT(IY)+FLOAT(IDAY)/365.-1995.)/5. F1=1.D0-F2 G10=29692.*F1+29619.4*F2 G11=-1784.*F1-1728.2*F2 H11= 5306.*F1+5186.1*F2 ELSEIF (IY.LT.2005) THEN !2000-2005 F2=(FLOAT(IY)+FLOAT(IDAY)/365.-2000.)/5. F1=1.D0-F2 G10=29619.4*F1+29554.63*F2 G11=-1728.2*F1-1669.05*F2 H11= 5186.1*F1+5077.99*F2 ELSEIF (IY.LT.2010) THEN !2005-2010 F2=(FLOAT(IY)+FLOAT(IDAY)/365.-2005.)/5. F1=1.D0-F2 G10=29554.63*F1+29496.57*F2 G11=-1669.05*F1-1586.42*F2 H11= 5077.99*F1+4944.26*F2 ELSEIF (IY.LT.2015) THEN !2010-2015 F2=(FLOAT(IY)+FLOAT(IDAY)/365.-2010.)/5. F1=1.D0-F2 G10=29496.57*F1+29441.46*F2 G11=-1586.42*F1-1501.77*F2 H11= 4944.26*F1+4795.99*F2 ELSEIF (IY.LT.2020) THEN !2015-2020 F2=(FLOAT(IY)+FLOAT(IDAY)/365.-2010.)/5. F1=1.D0-F2 G10=29441.46*F1+29404.8*F2 G11=-1501.77*F1-1450.9*F2 H11= 4795.99*F1+4652.5*F2 ELSE !2015-2020 DT=FLOAT(IY)+FLOAT(IDAY)/365.-2015. G10=29404.8-5.7*DT G11=-1450.9+7.4*DT H11= 4652.5-25.9*DT ENDIF C NOW CALCULATE THE COMPONENTS OF THE UNIT VECTOR EzMAG IN GEO COORD C SYSTEM: C SIN(TETA0)*COS(LAMBDA0), SIN(TETA0)*SIN(LAMBDA0), AND COS(TETA0) C ST0 * CL0 ST0 * SL0 CT0 SQ=G11**2+H11**2 SQQ=SQRT(SQ) SQR=SQRT(G10**2+SQ) SL0=-H11/SQQ CL0=-G11/SQQ ST0=SQQ/SQR CT0=G10/SQR STCL=ST0*CL0 STSL=ST0*SL0 CTSL=CT0*SL0 CTCL=CT0*CL0 C THE CALCULATIONS ARE TERMINATED IF ONLY GEO-MAG TRANSFORMATION C IS TO BE DONE (IHOUR>24 IS THE AGREED CONDITION FOR THIS CASE): 5 IF (IHOUR.GT.24) RETURN CALL SUN(IY,IDAY,IHOUR,MIN,ISEC,GST,SLONG,SRASN,SDEC) C S1,S2, AND S3 ARE THE COMPONENTS OF THE UNIT VECTOR EXGSM=EXGSE C IN THE SYSTEM GEI POINTING FROM THE EARTH'S CENTER TO THE SUN: S1=COS(SRASN)*COS(SDEC) S2=SIN(SRASN)*COS(SDEC) S3=SIN(SDEC) CGST=COS(GST) SGST=SIN(GST) C DIP1, DIP2, AND DIP3 ARE THE COMPONENTS OF THE UNIT VECTOR C EZSM=EZMAG IN THE SYSTEM GEI: DIP1=STCL*CGST-STSL*SGST DIP2=STCL*SGST+STSL*CGST DIP3=CT0 C NOW CALCULATE THE COMPONENTS OF THE UNIT VECTOR EYGSM IN THE SYSTEM C GEI BY TAKING THE VECTOR PRODUCT D x S AND NORMALIZING IT TO UNIT C LENGTH: Y1=DIP2*S3-DIP3*S2 Y2=DIP3*S1-DIP1*S3 Y3=DIP1*S2-DIP2*S1 Y=SQRT(Y1*Y1+Y2*Y2+Y3*Y3) Y1=Y1/Y Y2=Y2/Y Y3=Y3/Y C THEN IN THE GEI SYSTEM THE UNIT VECTOR Z=EZGSM=EXGSM x EYGSM=S x Y C HAS THE COMPONENTS: Z1=S2*Y3-S3*Y2 Z2=S3*Y1-S1*Y3 Z3=S1*Y2-S2*Y1 C THE VECTOR EZGSE (HERE DZ) IN GEI HAS THE COMPONENTS (0,-SIN(DELTA), C COS(DELTA)) = (0.,-0.397823,0.917462); HERE DELTA = 23.44214 DEG FOR C THE EPOCH 1978 (SEE THE BOOK BY GUREVICH OR OTHER ASTRONOMICAL C HANDBOOKS). HERE THE MOST ACCURATE TIME-DEPENDENT FORMULA IS USED: DJ=FLOAT(365*(IY-1900)+(IY-1901)/4 +IDAY)-0.5+FLOAT(ISEC)/86400. T=DJ/36525. OBLIQ=(23.45229-0.0130125*T)/57.2957795 DZ1=0. DZ2=-SIN(OBLIQ) DZ3=COS(OBLIQ) C THEN THE UNIT VECTOR EYGSE IN GEI SYSTEM IS THE VECTOR PRODUCT DZ x S DY1=DZ2*S3-DZ3*S2 DY2=DZ3*S1-DZ1*S3 DY3=DZ1*S2-DZ2*S1 C THE ELEMENTS OF THE MATRIX GSE TO GSM ARE THE SCALAR PRODUCTS: C CHI=EM22=(EYGSM,EYGSE), SHI=EM23=(EYGSM,EZGSE), C EM32=(EZGSM,EYGSE)=-EM23, AND EM33=(EZGSM,EZGSE)=EM22 CHI=Y1*DY1+Y2*DY2+Y3*DY3 SHI=Y1*DZ1+Y2*DZ2+Y3*DZ3 DECARG=SHI IF(ABS(DECARG).GT.1.) DECARG=SIGN(1.,DECARG) HI=ASIN(DECARG) C TILT ANGLE: PSI=ARCSIN(DIP,EXGSM) SPS=DIP1*S1+DIP2*S2+DIP3*S3 CPS=SQRT(1.-SPS**2) DECARG=SPS IF(ABS(DECARG).GT.1.) DECARG=SIGN(1.,DECARG) PSI=ASIN(DECARG) C THE ELEMENTS OF THE MATRIX MAG TO SM ARE THE SCALAR PRODUCTS: C CFI=GM22=(EYSM,EYMAG), SFI=GM23=(EYSM,EXMAG); THEY CAN BE DERIVED C AS FOLLOWS: C IN GEO THE VECTORS EXMAG AND EYMAG HAVE THE COMPONENTS C (CT0*CL0,CT0*SL0,-ST0) AND (-SL0,CL0,0), RESPECTIVELY. HENCE, IN C GEI SYSTEM THE COMPONENTS ARE: C EXMAG: CT0*CL0*COS(GST)-CT0*SL0*SIN(GST) C CT0*CL0*SIN(GST)+CT0*SL0*COS(GST) C -ST0 C EYMAG: -SL0*COS(GST)-CL0*SIN(GST) C -SL0*SIN(GST)+CL0*COS(GST) C 0 C THE COMPONENTS OF EYSM IN GEI WERE FOUND ABOVE AS Y1, Y2, AND Y3; C NOW WE ONLY HAVE TO COMBINE THE QUANTITIES INTO SCALAR PRODUCTS: EXMAGX=CT0*(CL0*CGST-SL0*SGST) EXMAGY=CT0*(CL0*SGST+SL0*CGST) EXMAGZ=-ST0 EYMAGX=-(SL0*CGST+CL0*SGST) EYMAGY=-(SL0*SGST-CL0*CGST) CFI=Y1*EYMAGX+Y2*EYMAGY SFI=Y1*EXMAGX+Y2*EXMAGY+Y3*EXMAGZ XMUT=(ATAN2(SFI,CFI)+3.1415926536)*3.8197186342 C THE ELEMENTS OF THE MATRIX GEO TO GSM ARE THE SCALAR PRODUCTS: C A11=(EXGEO,EXGSM), A12=(EYGEO,EXGSM), A13=(EZGEO,EXGSM), C A21=(EXGEO,EYGSM), A22=(EYGEO,EYGSM), A23=(EZGEO,EYGSM), C A31=(EXGEO,EZGSM), A32=(EYGEO,EZGSM), A33=(EZGEO,EZGSM), C ALL THE UNIT VECTORS IN BRACKETS ARE ALREADY DEFINED IN GEI: C EXGEO=(CGST,SGST,0), EYGEO=(-SGST,CGST,0), EZGEO=(0,0,1) C EXGSM=(S1,S2,S3), EYGSM=(Y1,Y2,Y3), EZGSM=(Z1,Z2,Z3) C AND THEREFORE: A11=S1*CGST+S2*SGST A12=-S1*SGST+S2*CGST A13=S3 A21=Y1*CGST+Y2*SGST A22=-Y1*SGST+Y2*CGST A23=Y3 A31=Z1*CGST+Z2*SGST A32=-Z1*SGST+Z2*CGST A33=Z3 10 FORMAT(/ * ' RECALC: GIVEN YEAR',I5,' IS OUT OF INTERVAL 1900-2010'/, * ' *** CALCULATIONS WILL BE DONE FOR YEAR =',I5,' ***'/) RETURN END C C C ********************************************************************* C SUBROUTINE SUN (IYEAR,IDAY,IHOUR,MIN,ISEC,GST,SLONG,SRASN,SDEC) C----------------------------------------------------------------------------- C CALCULATES FOUR QUANTITIES NECESSARY FOR COORDINATE TRANSFORMATIONS C WHICH DEPEND ON SUN POSITION (AND, HENCE, ON UNIVERSAL TIME AND SEASON) C C------- INPUT PARAMETERS: C IYR,IDAY,IHOUR,MIN,ISEC - YEAR, DAY, AND UNIVERSAL TIME IN HOURS, MINUTES, C AND SECONDS (IDAY=1 CORRESPONDS TO JANUARY 1). C C------- OUTPUT PARAMETERS: C GST - GREENWICH MEAN SIDEREAL TIME, SLONG - LONGITUDE ALONG ECLIPTIC C SRASN - RIGHT ASCENSION, SDEC - DECLINATION OF THE SUN (RADIANS) C ORIGINAL VERSION OF THIS SUBROUTINE HAS BEEN COMPILED FROM: C RUSSELL, C.T., COSMIC ELECTRODYNAMICS, 1971, V.2, PP.184-196. C C LAST MODIFICATION: MARCH 31, 2003 (ONLY SOME NOTATION CHANGES) C C ORIGINAL VERSION WRITTEN BY: Gilbert D. Mead C----------------------------------------------------------------------------- C DOUBLE PRECISION DJ,FDAY DATA RAD/57.295779513/ C IF(IYEAR.LT.1901.OR.IYEAR.GT.2099) RETURN FDAY=DFLOAT(IHOUR*3600+MIN*60+ISEC)/86400.D0 DJ=365*(IYEAR-1900)+(IYEAR-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) SLP=SLONG-9.924E-5 C C THE LAST CONSTANT IS A CORRECTION FOR THE ANGULAR ABERRATION DUE TO C THE ORBITAL MOTION OF THE EARTH C SIND=SOB*SIN(SLP) COSD=SQRT(1.-SIND**2) SC=SIND/COSD SDEC=ATAN(SC) SRASN=3.141592654-ATAN2(COS(OBLIQ)/SOB*SC,-COS(SLP)/COSD) RETURN END