      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 <GEOPACK.FOR> 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
