C $Id: icss_tilt_angle.f,v 1.5 2003/09/09 16:29:00 ad Exp ad $
C+
C
C UNIT NAME: ICSS_TILT_ANGLE
C
C PURPOSE: CALCULATE THE TILT_ANGLE
C
C UNIT TYPE: FORTRAN SUBROUTINE
C
C INVOCATION METHOD: CALL ICSS_TILT_ANGLE(ORB_POS_TIME,
C                                         TILT_ANGLE,
C					  DIPOLE_VECTOR,
C                                         RETURN_STATUS)
C
C ARGUMENT/QUALIFIER LIST:
C
C NAME                   TYPE      USE      DESCRIPTION 
C ----                   ----      ---      -----------
C ORB_POS_TIME(2)        I*4        I       REQUEST TIME:
C                                           (1) YYYYDOY
C                                           (2) MILLISECONDS OF DAY
C TILT_ANGLE             R*8        O       TILT ANGLE (RAD)
C DIPOLE_VECTOR(3)       R*8        O       GEO DIPOLE VECTOR
C                                           (1) X-COMPONENT
C                                           (2) Y-COMPONENT
C                                           (3) Z-COMPONENT
C RETURN_STATUS		 I*4        O       RETURN_STATUS
C                                           = SS$_NORMAL
C                                           = SS$_SSFAIL
C
C RETURNED_VALUE: NONE
C
C FILE/RECORD REFERENCES: NONE
C
C INCLUDE FILES: 
C 
C NAME	      DESCRIPTION
C ----        -----------
C $SSDEF      VMS DEFINED STATUS MESSAGE VALUES
C
C EXTERNAL VARIABLES: NONE
C
C EXTERNAL REFERENCES: 
C 
C NAME                    TYPE        DESCRIPTION
C ----                    ----        -----------
C ICSS_POS_OF_SUN         ICSS        RETURNS THE POSITION OF THE SUN FOR
C                                     A GIVEN TIME
C IC_GRNWCH_SIDEREAL      ICSS        RETURNS THE GREENWICH SIDEREAL
C                                     TIME IN RADIANS FOR A GIVEN TIME
C IC_GET_COEF             ICSS        RETURNS DIPOLE CONSTANTS
C
C ABNORMAL TERMINATION CONDITIONS, ERROR MESSAGES: NONE
C
C ERROR HANDLING: NONE
C
C ASSUMPTIONS, CONSTRAINTS, CONDITIONS:
C
C   (1) THIS ROUTINE IS CONSTRAINED BY THE YEARS 1965 TO 2020.  THE ROUTINE
C       USES THE YEAR 2020 FOR ALL CALCULATIONS IF THE THE INPUT YEAR IS 
C       GREATER THEN 2020.  LIKEWISE, THE ROUTINE WILL USE THE YEAR 1965
C       IF THE INPUT YEAR IS LESS THEN 1965,
C
C DEVELOPMENT HISTORY:
C
C AUTHOR       CHANGE ID     RELEASE      DATE     DESCRIPTION OF CHANGE
C ------       ---------     -------      ----     ---------------------
C K. HUGHES     CCR_2171      R6.3      01/30/96   ORIGINAL DEVELOPMENT
C B. SAMUELSON  UNIX PORT               09/18/96   Removed VMS specific
C                                                  code (referneces to
C                                                  SS$NORMAL)
C A.Davis 01/05/2011 Changed limit from 2010 to 2015
C A.Davis 01/05/2016 Changed limit from 2015 to 2020(and updated coeffs )
C A.Davis 01/03/2020 Changed limit from 2020 to 2025(and updated coeffs )
C
C NOTES:
C 
C     (1)  THIS ROUTINE WAS BASED ON THE ROUTINE RECALC IN THE TSYGANENKO
C          MODEL WRITTEN BY NIKOLAI TSYGANENKO AND MAURICIO PEREDO.
C
C-
C
C PDL:
C
C  	EXTRACT DOY AND YEAR FOR ORB_POS_TIME
C       IF YEAR IS LESS THEN 1965 CHANGE IT TO 1965
C       IF YEAR IS GREATER THEN 2020 THEN CHANGE IT TO 2020
C       GET CONSTANTS BASED ON YEAR
C       COMPUTE GEOGRAPHIC DIPOLE MOMENT
C       CALL ICSS_POS_OF_SUN TO GET POSITION OF SUN
C       CONVERT THE SUN_VECT TO UNIT VECTOR
C	CALL IC_GRNWCH_SIDEREAL TO GET GREENWICH SIDEREAL TIME
C       CALCULATE TILT ANGLE
C       END
C
C--+
C	OPTIONS/EXTEND_SOURCE
	SUBROUTINE ICSS_TILT_ANGLE(ORB_POS_TIME,
     .                             TILT_ANGLE,
     .                             DIPOLE_VECTOR,
     .                             RETURN_STATUS)
	IMPLICIT NONE
C
C       -- ARGUMENT DEFINITIONS
C 
	INTEGER*4   ORB_POS_TIME(2) 
	REAL*8	    TILT_ANGLE
	REAL*8	    DIPOLE_VECTOR(3)
	INTEGER*4   RETURN_STATUS
C
C 	-- LOCAL DEFINITIONS
C
	INTEGER*4   IYR, IDAY
	INTEGER*4   STATUS
	INTEGER*4   LAST_DOY_LOADED /0/, LAST_YEAR_LOADED /0/
	INTEGER*4   YEAR, MONTH, DAY   !HOLDS YEAR, MONTH, DAY FOR DATE CONV.

	REAL*8      RAD /57.295779513D0/
	REAL*8	    S1, S2, S3
	REAL*8	    F1, F2
	REAL*8	    G10, G11
	REAL*8	    H11
	REAL*8	    DT, SQ, SQQ, SQR
	REAL*8	    SL0, CL0, ST0, CT0
	REAL*8	    STCL, STSL

	REAL*8	    SUN_VECT(3), SUN_VEL(3)
	REAL*8	    RMAG_SUN_VECT
	REAL*8	    GST, CGST, SGST
	REAL*8	    DIP1, DIP2, DIP3
	REAL*8	    SPS

	REAL*8	    G(66), H(66)
C
C       ***  BEGIN EXECUTABLE CODE  ***
C
C
C       -- EXTRACT YEAR AND DAY OF YEAR FROM ORB_POS_TIME
C
	IDAY  = MOD(ORB_POS_TIME(1), 1000)
	IYR = (ORB_POS_TIME(1) - IDAY) / 1000
C
C       -- YEAR MUST BE BETWEEN 1965 AND 2025
C
	IF (IYR .LT. 1965) THEN
	   IYR = 1965
	   PRINT *, ' YEAR IS LESS THAN 1965. USING 1965 FOR CALCS.'
	ENDIF
	IF (IYR .GT. 2025) THEN 
	   IYR = 2025
           PRINT *, ' YEAR IS GREATER THAN 2025. USING 2025 FOR CALCS.'
	ENDIF
C
C	-- CALL IC_GET_COEF TO GET COEFFICIENTS FOR DIPOLE CALCULATION
C
	IF ((LAST_YEAR_LOADED .NE. IYR) .OR. 
     .      (LAST_DOY_LOADED  .NE. IDAY)) THEN
	    LAST_YEAR_LOADED = IYR
	    LAST_DOY_LOADED  = IDAY
	    CALL IC_GET_COEF(IYR, IDAY, G, H, G10, G11, H11, STATUS)
	ENDIF
C
C	-- BEGIN CALCULATING DIPOLE AND TILT ANGLE
C
        SQ=G11**2+H11**2
        SQQ=dSQRT(SQ)
        SQR=dSQRT(G10**2+SQ)
        SL0=-H11/SQQ
        CL0=-G11/SQQ
        ST0=SQQ/SQR
        CT0=G10/SQR
        STCL=ST0*CL0
        STSL=ST0*SL0

	DIPOLE_VECTOR(1) = STCL
	DIPOLE_VECTOR(2) = STSL
	DIPOLE_VECTOR(3) = CT0
C
C       -- CALCULATE THE POSITION OF THE SUN...USING THE ICSS ROUTINES...
C
        CALL ICSS_POS_OF_SUN (ORB_POS_TIME, SUN_VECT, SUN_VEL, STATUS)
	IF (STATUS .EQ. 1) THEN
	     RETURN_STATUS = 0 
	ELSE
	     RETURN_STATUS = 0 
	ENDIF
C
C	-- CONTINUE IF STATUS IS GOOD
C
	IF (RETURN_STATUS .EQ. 0) THEN
C
C          -- CONVERT THE SUN_VECT INTO A UNIT VECTOR
C
           RMAG_SUN_VECT = DSQRT( (SUN_VECT(1)**2)+
     .                            (SUN_VECT(2)**2)+
     .                            (SUN_VECT(3)**2) )
C
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:
C
           S1 = SUN_VECT(1)/RMAG_SUN_VECT
           S2 = SUN_VECT(2)/RMAG_SUN_VECT
           S3 = SUN_VECT(3)/RMAG_SUN_VECT
C
           CALL IC_GRNWCH_SIDEREAL(ORB_POS_TIME, GST)
           CGST=dCOS(GST)
           SGST=dSIN(GST)
C
C          -- DIP(1), DIP(2), AND DIP(3) ARE THE COMPONENTS OF THE UNIT VECTOR 
C	   -- EZSM=EZMAG IN THE SYSTEM GEI:
C
           DIP1 = STCL*CGST-STSL*SGST
           DIP2 = STCL*SGST+STSL*CGST
           DIP3 = CT0
C
C          -- TILT ANGLE: PSI=ARCSIN(DIP,EXGSM)
C
           SPS=DIP1*S1+DIP2*S2+DIP3*S3
           TILT_ANGLE=dASIN(SPS)
	ENDIF

	RETURN
	END
