C $Id: ic_get_nut_angles.f,v 1.2 1998/07/24 21:54:31 asc Exp $
CCCC
C
C  IC_GET_NUT_ANGLES - Returns angles that are necessary to adjust the
C		       Greenwich Hour angle to true of date
C 
C  PURPOSE : THIS SUBROUTINE CALCULATES, THROUGH APPROPRIATE ANALYTIC
C            EXPRESSIONS, VALUES FOR THE NUTATION
C            ANGLES TO ROTATE FROM MEAN OF
C            JULIAN 2000 TO TRUE OF DATE.
C
C  UNIT TYPE:  SUBROUTINE
C
C  INVOCATION METHOD:  CALL IC_GET_NUT_ANGLES (time,
C					       deleps,
C					       delpsi,
C					       eps)
C
C  ARGUMENT LIST:
C
C  NAME	      TYPE    USE	DESCRIPTION
C  ----       ----    ---       -----------
C  TIME       R*8     I         TIME IN JULIAN CENTURIES OF 36525.0
C                               MEAN SOLAR DAYS FROM J2000. (NOTE: THIS
C                               CAN BE POSITIVE OR NEGATIVE.)
C  DELEPS     R*8     O	        DELTA EPSILON, Nutation in obliquity
C  DELPSI     R*8     O         DELTA PSI, Nutation in longitude
C  EPS        R*8     O         EPSILON
C
C  FILE/RECORD REFERENCES:  NONE
C  
C  EXTERNAL VARIABLES:  NONE
C
C  EXTERNAL REFERENCES:
C	MULMAT - routine that multiplies two matrices
C
C  ABNORMAL TERMINATION CONDITIONS, ERROR MESSAGES:  NONE
C
C  ASSUMPTIONS, CONSTRAINTS, RESTRICTIONS:
C
C  DEVELOPMENT HISTORY
C
C  AUTHOR	CHANGE ID	RELEASE	  DATE	  DESCRIPTION OF CHANGE
C  ------	---------	-------   ----	  ---------------------
C  C. Raymond                                     Original Algorithm, Coding
C  J. Lubelczyk			B3    	  10/91	  Prolog, Coding
C  J. LUBELCZYK ICCR #83, CCR #'S 130, 137 11/91  B3 Update
C  C. Raymond   CCR 568         B4/R2.1   02/92   Updated to increase the
C						   computational efficiency
C                                                  of the coordinate system
C                                                  conversion package.	
C  B. SAMUELSON SPOF PORT       NONE     04/11/94  Change: required to port
C     (CSC)     (see notes)                          icss routines to SPOF
C
C  NOTES:
C
C  1. (CCR 568)  The nutation series in this routine are computationally
C     expensive.  Computation times were reduced by taking advantage of the
C     relatively slowly varying nature of the nutation angles by only 
C     recomputing every 1/2 day.  This gives an error that is bounded by
C     0.1 arc-seconds, which is certainly acceptable. 	
C  2. The changes recorded under ID SPOF-PORT are required to make the
C     ICSS coordinate conversion routines, originally developed under
C     VAX-VMS 5.4 run on the UNIX-based workstations of the SPOF.  (Sun
C     SPARCstations and DEC DECstations).  The changes are as follows:
C        a.  Delete references to ICSS_INC
C        b.  Define Message texts or files to correspond to the messages
C            ICSS_SUCCESSFUL, etc. which are embedded in the error handling.
C        c.  Remove references to the NAG routines F01CRF and F01CKF (matrix
C            transposition and matrix multiplication routines).
C  3. In addition, to successfully run the software packages, copies of the
C     Solar/Lunar/Planetary (SLP) file and timing coefficients file (TCC)
C     must be ported onto the SPOF.
C
CCCC
C
C  PDL:
C
C  Calculate the fundamental arguments for use in nutation calculations
C
C  Calculate the mean obliquity of date
C
C  Calculate the nutation in longitude and nutation in obliquity
C
C  Compute DELPSI and DELEPS
C
C  Calculate the true obliquity of date, EPS
C
C  return
C
CCCCC
        SUBROUTINE IC_GET_NUT_ANGLES (TIME,DELEPS,DELPSI,EPS)
C
        IMPLICIT    NONE
C
C*  Calling Parameters
C
	REAL*8	    TIME    !Time in Julian Centuries of 36525.0 mean solar
C			    !days from J2000.
	REAL*8	    DELEPS  !DELTA EPSILON, Nutation in obliquity
	REAL*8	    DELPSI  !DELTA PSI, Nutation in longitude
	REAL*8	    EPS	    !EPSILON
	REAL*8      TOL
	parameter (TOL = 0.5d0/36525.0d0 )
C
C*  Other Variables
C
        REAL*4      ISINCO(106,2)   !Array used in nutation calculations
	REAL*4	    ICOSCO(106,2)   !Array used in nutation calculations
	REAL*8	    FUNARG(106,5)   !Array used to get the fundamental args.
	REAL*8	    FUND(5,1)	    !The fundamental arguments
	REAL*8	    ARG(106,1)	    !Arguments of sine and cosine
	REAL*8	    T(2,1)	    !Time
	INTEGER*4   IFUNAR(106,5)   !Array used in nutation calculations
	REAL*8	    SINCOF(106,2)   !R*8 Array used in nutation calculations
	REAL*8	    COSCOF(106,2)   !R*8 Array used in nutation calculations
	REAL*8	    COFSIN(106,1)   !Coefficients of sine
	REAL*8	    COFCOS(106,1)   !Coefficients of cosine
	REAL*8	    SINP(106,1)	    !Used to calc. PSI
	REAL*8	    COSE(106,1)	    !Used to calc. EPS
	REAL*8	    T2		    !TIME squared
	REAL*8	    T3		    !TIME cubed
	INTEGER*4   E, F	    !Loop Counters
	REAL*8	    PI		    !3.14159...
	REAL*8	    DTR		    !Conversion factor, Degree to Radians
	REAL*8	    STR		    !Conversion factor, Seconds to Radians
	REAL*8	    R		    !
	REAL*8	    EPSO	    !Mean obliquity of date
	REAL*8	    SUMPSI	    !Used to calc PSI
	REAL*8	    SUMEPS	    !Used to calc EPS
	REAL*8	    PRODPS	    !Used to calc PSI
	REAL*8	    PRODEP	    !Used to calc EPS

        include 'nut_ang.cmn'
C
C                                                                      
C     INITIALIZE VALUES OF IFUNAR, ISINCO, AND ICOSCO ARRAYS FOR USE IN
C     NUTATION CALCULATIONS.                                           
C                                                                      
C     THE 1980 IAU THEORY OF NUTATION,CONTAINED IN JPL
C     DE200 PLANETARY EPHEMERIS.
C
      DATA IFUNAR/0,0,2,-2,2,0,2,-2,2,0,2,2,2,0,2,0,0,2,0,0,2,0,2,
     ;0,0,-2,-2,0,0,2,2,0,2,2,0,2,0,0,0,2,2,2,0,2,2,2,
     ;2,0,0,2,0,2,2,2,0,2,0,2,2,0,0,2,0,-2,0,0,2,2,2,
     ;0,4*2,3*0,2,0,0,2,2,0,3*2,4,0,2,2,0,4,3*2,0,-2,2,0,-2,2,0,-2,0,
     ;2,0,1,2,1,0,2,0,1,1,2,0,2,2,1,0,0,0,1,2,1,1,1,1,1,0,0,
     ;1,0,2,1,0,2,0,1,2,0,2,0,1,1,2,1,2,0,2,2,0,1,1,1,1,
     ;0,2,2,2,0,2,1,1,1,1,0,1,0,0,0,0,0,2,2,1,3*2,2*1,
     ;2,0,2,2,0,2,2,0,2,1,2,2,0,1,2,1,2,2,0,3*1,2,2*0,2*1,2*0,2,0,
     ;0,0,0,0,0,-1,-2,0,0,1,1,-1,0,0,0,2,1,2,-1,0,-1,0,1,
     ;0,1,0,1,1,0,1,0,0,0,0,0,0,0,13*0,1,1,-1,0,0,0,0,0,0,
     ;0,-1,0,1,0,0,1,0,-1,-1,2*0,-1,1,10*0,1,3*0,-1,6*0,1,-1,2*0,1,0,
     ;-1,1,3*0,1,0,0,-2,2,-2,1,0,2,0,0,0,0,0,2,0,0,0,0,0,-2,0,2,0,1,
     ;2,0,0,0,-1,0,0,1,0,1,1,-1,0,1,-1,-1,1,0,2,1,2,0,-1,
     ;-1,1,-1,1,0,0,1,1,2,0,0,1,0,1,2,0,1,0,1,1,1,-1,
     ;-2,3,0,1,-1,2,1,3,0,-1,1,-2,-1,2,1,1,-2,-1,1,2*2,1,0,3,1,0,-1,
     ;3*0,1,0,2*1,2,7*0,-1,-2,0,-2,0,-2,-2,-2,-2,-2,0,0,-2,0,2,-2,-2,
     ;-2,-1,-2,2,2,0,1,-2,0,0,0,0,-2,0,2,0,0,2,0,2,0,-2,0,0,
     ;0,2,-2,2,-2,0,0,2,2,-2,2,2,-2,-2,0,0,-2,0,1,0,0,0,2,
     ;2*0,2,0,-2,3*0,1,0,-4,2,4,-4,-2,2,4,0,2*-2,2*2,
     ;3*-2,0,2,0,-1,2,-2,0,-2,2*2,4,1/
C
      DATA ISINCO/
     ;-171996.,2062.,46.,11.,-3.,-3.,-2.,1.,-13187.,1426.,
     ;-517.,217.,129.,48.,-22.,17.,-15.,-16.,-12.,-6.,-5.,
     ;4.,4.,-4.,1.,1.,-1.,1.,1.,-1.,-2274.,712.,-386.,-301.,
     ;-158.,123.,63.,63.,-58.,-59.,-51.,-38.,29.,29.,-31.,26.,
     ;21.,16.,-13.,-10.,-7.,7.,-7.,-8.,6.,6.,-6.,-7.,6.,-5.,
     ;5.,-5.,-4.,4.,-4.,-3.,3.,-3.,-3.,-2.,-3.,-3.,2.,-2.,
     ;2.,-2.,2*2.,1.,-1.,1.,-2.,-1.,1.,2*-1.,3*1.,2*-1.,
     ;2*1.,-1.,2*1.,7*-1.,1.,-1.,1.,
     ;
     ;-174.2,.2,6*0.,-1.6,-3.4,1.2,-.5,.1,0.,0.,-.1,0.,.1,
     ;12*0.,-.2,.1,-.4,4*0.,.1,-.1,67*0./
C
      DATA ICOSCO/
     ;92025.,-895.,-24.,0.,1.,0.,1.,0.,5736.,54.,224.,-95.,
     ;-70.,1.,0.,0.,9.,7.,6.,3.,3.,-2.,-2.,7*0.,977.,-7.,
     ;200.,129.,-1.,-53.,-2.,-33.,32.,26.,27.,16.,-1.,-12.,
     ;13.,-1.,-10.,-8.,7.,5.,0.,-3.,3.,3.,0.,-3.,3.,3.,-3.,
     ;3.,0.,3.,5*0.,5*1.,-1.,1.,-1.,1.,0.,2*-1.,0.,-1.,
     ;1.,0.,-1.,2*1.,2*0.,-1.,17*0.,
     ;
     ;8.9,.5,6*0.,-3.1,-.1,-.6,.3,18*0.,-.5,0.,0.,-.1,72*0./
C
      DATA R/1296000.0D0/
      IF ( DABS ( TIME - OLDTIM ) .le. tol ) THEN
         DELEPS = OLDDEP 
         DELPSI = OLDDPS 
         EPS    = OLDEPS
         RETURN
      ENDIF
      T2 = TIME*TIME
      T3 = TIME*T2
C
C                                                                       
C     CONVERT IFUNAR, ISINCO, AND ICOSCO ARRAYS TO REAL*8 ARRAYS FUNARG,
C     SINCOF, AND COSCOF, RESPECTIVELY.                                 
C                                                                       
      DO 1 E=1,5                                                       
      DO 1 F=1,106                                                     
    1 FUNARG(F,E)=real(IFUNAR(F,E))                               
      DO 2 E=1,2                                                       
      DO 2 F=1,106                                                     
      SINCOF(F,E)=DBLE(ISINCO(F,E))                                 
    2 COSCOF(F,E)=DBLE(ICOSCO(F,E))                                 
C                                                                       
C     CALCULATE CONVERSION FACTORS: DEGREES TO RADANS (DTR), SECONDS TO 
C     RADIANS (STR)                                                     
C                                                                       
      PI=4.0D0*DATAN(1.D0)                                              
      DTR=PI/180.D0                                                     
      STR=DTR/3600.D0                                                   
C
C     BEGIN COMPUTATION OF NUTATION IN OBLIQUITY AND LONGITUDE          
C
C
C     CALCULATE FUNDAMENTAL ARGUMENTS FOR USE IN NUTATION CALCULATIONS
C     TIME IS REFERENCED TO J2000.0.
C     FUND(1,1)= F
C     FUND(2,1)= OMEGA
C     FUND(3,1)= L PRIME
C     FUND(4,1)= L
C     FUND(5,1)= D
C
      FUND(1,1)=STR*(335778.877D0+(1342.0D0*R+295263.137D0)*TIME
     * -13.257D0*T2+1.1D-2*T3)
      FUND(2,1)=STR*(450160.280D0-(5.D0*R+482890.539D0)*TIME+
     * 7.455D0*T2+8.0D-3*T3)
      FUND(3,1)=STR*(1287099.804D0+(99.0D0*R+1292581.224D0)*TIME-
     * 5.77D-1*T2-1.2D-2*T3)
      FUND(4,1)=STR*(485866.733D0+(1325.0D0*R+715922.633D0)*TIME+
     * 31.310D0*T2+6.4D-2*T3)
      FUND(5,1)=STR*(1072261.307D0+(1236.0D0*R+1105601.328D0)*TIME-
     * 6.891D0*T2+1.9D-2*T3)
C
C     CALCULATE MEAN OBLIQUITY OF DATE (EPSO). WHERE TIME IS MEASURED IN
C     JULIAN CENTURIES FROM 2000.0.
C
      EPSO=(1.813D-3*T3-5.9D-4*T2
     *     -4.6815D+1*TIME+8.4381448D+4)*STR
C                                                                       
C     CALCULATE NUTATION IN LONGITUDE (DELPSI) AND NUTATION IN OBLIQUITY
C     (DELEPS).  THIS IS A THREE STEP PROCESS:                          
C     (1) CALCULATE ARGUMENTS OF SINE (FOR DELPSI) AND COSINE (FOR DELEPS)
C         THESE ARE OF THE FORM                                        
C                                                                      
C         ARG = SUMMATION ( A(I) * FUND(I,1) ), I = 1,5                 
C                                                                       
C         WHERE THE A(I)'S ARE ELEMENTS OF FUNARG.                      
C                                                                       
      CALL MULMAT (FUNARG, FUND, 106, 5, 1, ARG)
C                                                                       
C     (2) CALCULATE COEFFICIENTS OF SINE AND COSINE, WHICH ARE THE PRODUCTS
C         OF SINCOF * T AND COSCOF * T.  THESE COEFFICIENTS ARE IN UNITS
C         OF 0.0001 SECONDS OF ARC.                                     
C                                                                       
      T(1,1)=1.D0                                                       
      T(2,1)=TIME                                                       
      CALL MULMAT (COSCOF, T, 106, 2, 1, COFCOS)
      CALL MULMAT (SINCOF, T, 106, 2, 1, COFSIN)
      DO 4 F=1,106                                                     
      COFCOS(F,1)=COFCOS(F,1)*1.D-4                                   
    4 COFSIN(F,1)=COFSIN(F,1)*1.D-4                                   
C                                                                       
C     (3) CALCULATE THE SINES AND COSINES OF THE ARGUMENTS AND MULTIPLY 
C         BY THEIR COEFFICIENTS, THEN ADD.  COMPUTE DELPSI AND DELEPS.  
C                                                                       
      SUMPSI=0.D0                                                       
      SUMEPS=0.D0                                                       
      DO 5 E=1,106                                                     
      SINP(E,1)=DSIN(ARG(E,1))                                        
      COSE(E,1)=DCOS(ARG(E,1))                                        
      PRODPS=COFSIN(E,1)*SINP(E,1)                                    
      PRODEP=COFCOS(E,1)*COSE(E,1)                                    
      SUMPSI=SUMPSI+PRODPS                                              
    5 SUMEPS=SUMEPS+PRODEP                                              
      DELEPS=SUMEPS*STR                                                 
      DELPSI=SUMPSI*STR                  
C                                         
C     CALCULATE TRUE OBLIQUITY OF DATE (EPS).       
C                                                   
      EPS=EPSO+DELEPS
      OLDDEP = DELEPS
      OLDDPS = DELPSI
      OLDEPS = EPS
      OLDTIM = TIME
C
      RETURN                                                     
C
      END                                                        
