C $Id: ic_pos_vel_of_sun.F,v 1.6 2007/03/23 00:52:31 asc Exp $
CCC
C
C  IC_POS_VEL_OF_SUN - return the instantaneous position and velocity of the sun
C 
C  PURPOSE:  Calculate the position and velocity vectors 
C            of the sun in GCI coordinates.
C
C  UNIT TYPE:  SUBROUTINE
C
C  INVOCATION METHOD:  CALL IC_POS_OF_SUN (orb_pos_time, 
C                                          sun_pos_vel,
C					   pv_status)
C
C  ARGUMENT LIST:
C
C  NAME	                  TYPE   USE  DESCRIPTION
C  ----                   ----   ---  -----------
C  ORB_POS_TIME(2)        I*4    I    TIME OF ORB. VECTOR, YEAR-DAY-MILLI OF DAY
C  SUN_POS_VEL(6)         R*8    O    POSITION AND VELOCITY OF THE SUN VECTORS
C  PV_STATUS		  I*4	 O    STATUS OF GETTING POS AND VEL OF SUN
C
C  FILE/RECORD REFERENCES:  
C
C  NAME 	  USE	    DESCRIPTION
C  ----	    	  ---       -----------
C  SLP FILE	  I	    SOLAR/LUNAR/PLANETARY EPHEMERIS (SLP) FILE
C      
C  EXTERNAL VARIABLES:
C
C  SOURCE	NAME	       USE	DESCRIPTION
C  ------       ----           ---      -----------
C  ORLSPR       TSEC			Time of day in seconds
C		PPOLY1(3,20)		Position poly. coeff. for fast non-central body
C		VPOLY1(3,13)		Velocity poly. coeff. for fast body
C		PPOLY2(3,13,7)		
C		APOLY(3,3,10)		Not Used, Poly Coeffs. for A matrix
C		CPOLY(3,3,10)		Polynomial coefficients for the C matrix	
C		PDELH(10)		Coef. for nutation of Greenwich Hour Angle
C		IDAY			Day # within SLP of current record pair
C		IDAY1			1st day of the SLP file
C		IYEAR			1st year of the SLP file
C		ISPAN			Span of the SLP file in days
C		NBEPM(9)		Body index #s of all bodies listed in SLP file
C		NDEGRE(4)		degree of polynomial expansion for 
C					each quantity given in the SLP file
C		NCFDAY		
C		ISLP50			coordinate system of SLP file
C		NBSLP			# of bodies in SLP file (usually 9)
C		IHEADR			Indicates if the SLP header has been read
C
C  EXTERNAL REFERENCES:  
C    IC_FETCH_SLP_COEF - RETRIEVE CHEBYSHEV POLYNOMIAL COEFFICIENTS FOR
C                        A TIME INTERVAL WHICH INCLUDES (ORB_POS_TIME)
C    IC_OPEN_TIME_COEFF_FILE - OPEN TIMING COEF FILE
C    IC_GET_TIME_DIFF - GET TIME DIFFERENCE BETWEEN A1 AND UTC TIMING SYSTEMS
C    ERRSET -           TRAP OVERFLOW AND UNDERFLOW ERRORS
C    OPEN -             FORTRAN OPEN
C    CLOSE -            FORTRAN CLOSE
C
C  ABNORMAL TERMINATION CONDITIONS, ERROR MESSAGES:
C    PV_STATUS is returned as one of the following:
C	0 - Successful completion
C	1 - Error getting a logical unit # for the SLP file
C	2 - Error opening the SLP file
C	3 - Error getting a LUN for the Timing Coef. file
C	4 - SLP file does not contain the body requested
C    	Error status returned from IC_OPEN_TIME_COEFF_FILE
C    	Error status returned from IC_GET_TIME_DIFF
C    	Error status returned from IC_FETCH_SLP_COEF
C
C  ASSUMPTIONS, CONSTRAINTS, RESTRICTIONS:  NONE
C
C  DEVELOPMENT HISTORY
C
C  AUTHOR	CHANGE ID	RELEASE	  DATE	    DESCRIPTION OF CHANGE
C  ------	---------	-------   ----	    ---------------------
C  J. LUBELCZYK                 B1R1      11/29/90  INITIAL PDL
C  J. LUBELCZYK                 B1R1      12/11/90  CODING
C  C. RAYMOND (CSC)                       07/25/91  ALTERED THE BASIC
C                                                    NATURE OF THE ALGORITHM
C                                                    USED TO GET THE SUN
C                                                    POSITION.  SEE NOTES.
C  J. LUBELCZYK		                  09/04/91  Corrected algorithm to
C						    meet SSDM standards
C  J. LUBELCZYK				  09/11/91  Modified to return both
C						    position and velocity. Also,
C						    the pos. is no longer a unit
C						    vector.  Add error handling.
C  J. LUBELCZYK ICCR #83, CCR #'S 130, 137 11/91    B3 update
C
C  M. Harris                     B3R2     12/03/91  open slp file as readonly
C  B. SAMUELSON SPOF PORT        NONE     04/11/94  Change: required to port
C     (CSC)     (see notes)                          icss routines to SPOF
C
C  A. Davis                               08/27/07  Force read of A1UTC each new day
C  A. Davis                               10/16/07  Correct A1UTC after 1997Jun30
C
C  NOTES:
c  1)  The current version of this program was written by C. Raymond, making 
C      heavy use of code taken from the EPHGEN program to read the Solar/
c      Lunar/Planetary Ephemeris file (SLP file).  The SLP file contains
c      coefficients of Chebyshev polynomials which give high-precision
c      Sun, Moon, etc. positions.
C
C  2)  The COMMON block /ORSLPR/ is introduced only as a means for holding
C      the large amount of data that each pair of SLP file records contains
C      for a 10-day span of time.
C  
C  3)  This routine now returns both the positon and velocity of the sun.
C
C  4)  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  5)  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  Convert ORB_POS_TIME to year, day, and seconds of day
C
C  CONVERT resulting time to Julian ephemeris time
C
C  IF ( initialization flag set ) 
C
C  THEN
C     
C     Call lib$get_lun to get a logical unit number for the SLP file
C     If error getting a logical unit number then
C        set PV_STATUS to indicate error getting a LUN for the SLP file
C	 abort to 999
C     Endif
C
C     OPEN SLP FILE
C     If error opening the SLP file then
C        set PV_STATUS to indicate error opening the SLP file
C        abort to 999
C     endif
C 
C     Call lib$get_lun to get a logical unit number for the timing coef. file
C     If error getting a logical unit number then
C        set PV_STATUS to indicate error getting a LUN for the timing coef file
C	 abort to 991
C     Endif
C
C     Call IC_OPEN_TIME_COEFF_FILE to open the timing coef file
C     If error returned from IC_OPEN_TIME_COEFF_FILE abort to 991
C
C     Call IC_GET_TIME_DIFF to get the offset between A1 and UTC timing systems
C     If error returned from IC_GET_TIME_DIFF abort to 990
C
C     RESET initialization flag
C
C  ENDIF
C
C  Add A1UTC offset to TIME
C
C  CALL IC_FETCH_SLP_COEF TO FETCH SLP file coefficients
C  If error return from IC_FETCH_SLP_COEF abort to 990
C
C  CALCULATE SUN position
C
C990 CONTINUE
C  IF ERROR HAS OCCURRED THEN
C     CLOSE TIMING COEF. FILE
C  ENDIF
C
C991 CONTINUE
C  IF ERROR HAS OCCURRED THEN
C     CLOSE SLP FILE
C  ENDIF
C
C999 CONTINUE
C
C  RETURN
C
CCCC
	subroutine IC_POS_VEL_OF_SUN (ORB_POS_TIME, 
     -                                SUN_POS_VEL, PV_STATUS)
C
	IMPLICIT none
C
C*  Calling parameters
C
	integer*4   ORB_POS_TIME(2) !ISTP time format (YYYYDDD, Milli of Day)
        real*8      SUN_POS_VEL(6)  !Position and velocity of sun vectors
	integer*4   PV_STATUS	    !Status of getting the pos and vel of sun
C
C*  Other Variables
C
 	integer*4   YEAR        !Year - YYYY
 	integer*4   DAY	        !day of year
 	integer*4   YRDOYOLD /0/  !previous year/day of year
        real*8      EPT(19)	!successive powers of eptime up to 19
	real*8	    TEMP(3,2)   !Temporary working array for position
	real*8	    TEMPV(3,2)  !Temporary working array for velocity
        real*8      TOLD /0.D0/ !Previous time, old
	integer*4   ICB	/1/     !Body index number of the earth
	integer*4   INCB /3/    !Body index number of the sun
	real*8      SECS        !Seconds
	logical	    INIT_READ /.false./ !Initialized read of SLP file
        integer*2   MAXLIM /15/ !Control error trapping
	real*8	    DXJUL	!Julian ephemeris date at beginning of year
	integer*4   I, J, K	!Loop Variables
	integer*4   IPP		!Position in slow body polynomial array
	integer*4   LUSLP /10/	!Logical Unit number for the SLP file
	integer*4   LUCOEF /11/	!Logical Unit number for the Timing Coef. file
	real*8	    A1UTC	!Offset between UTC time and A1 time in seconds
	real*8	    FDAY	!Fraction of a day
	real*8	    DJ		!Julian Day
	real*8	    TIME	!Seconds from start of the SLP file
	real*8	    TIME_DIFF(2)!Time differences from IC_GET_TIME_DIFF
	real*8	    EPTIME	!Time variable to be used for the evaluation
C				!of the chebyshev polynomial
	integer*4   MAX		!maximum degree of all polynomials used
        integer*4   LENGTH      !record size in longwords for time coeff file
	integer*4   ICHK	!Temporary variable
	integer*4   NDEG2	!Polynomial order for a fast body
	integer*4   NDEG4	!Polynomial order for a slow body
	real*8	    ETK1	!Used w/ EPT
	real*8	    ETK2	!used w/ EPT
	integer*4   STATUS	!General status
        integer*4   LEN 

	integer*4   ieeer
        integer*4   ieee_flags
        character*(15)  out
        character*(100) filenm

C
C*  System routines
C
C
	include	    'orslpr.cmn'
C
C*  START EXECUTABLE STATEMENTS
C
C*  Convert the given millisecond of day [orb_pos_time(2)] to second of day.
C*  Convert the packed form into year and day-of-year
C
	SECS = (real(ORB_POS_TIME(2)))/1000.0D0
        YEAR = ORB_POS_TIME(1)/1000
	DAY  = MOD(ORB_POS_TIME(1),1000)
	FDAY = SECS/86400.0D0
	DJ = DXJUL(YEAR) + real(DAY)+FDAY
C
C TAKE CARE OF UNDERFLOW/OVERFLOW CONDITIONS THAT CAN SOMETIMES RESULT.
C
C
C          	CHECK INITIALIZATION 
C
C	print *, 'Initial PV_STATUS is', PV_STATUS

      IF ( .NOT. INIT_READ ) THEN
        LENGTH  = 248
        LEN  = 2264
#ifdef DATDIR
       filenm = DATDIR // '/m2000.dat'
#else       
       filenm = './m2000.dat'
#endif

C       print *, 'm2000  is ', filenm


	
        OPEN (UNIT=LUSLP,FORM='UNFORMATTED',ACCESS='DIRECT',CONVERT="BIG_ENDIAN",
     1        FILE=filenm,
     2        STATUS='OLD', RECL=LEN, IOSTAT = status)
	if (status .ne. 0) then
C	   print *, 'Whoops!', status
	   pv_status = 2 
	   goto 999
	endif
C
	call IC_OPEN_TIME_COEFF_FILE (LUCOEF, LENGTH, PV_STATUS)
C	print *, 'PV_STATUS is', PV_STATUS
	IF (PV_STATUS .NE. 0) GOTO 991
C
	INIT_READ = .TRUE.
        IHEADR = 0
      ENDIF
C
C*   Get difference between A1 and UTC timing systems...account for leap seconds
C
      IF ( YRDOYOLD .NE. ORB_POS_TIME(1) ) THEN
C
C        print *, 'calling IC_GET_TIME_DIFF'
        call IC_GET_TIME_DIFF(LUCOEF, DJ, TIME_DIFF, PV_STATUS)
C	print *, 'IC_GET_TIME_DIFF status is ', PV_STATUS
        IF (PV_STATUS .NE. 0) GOTO 990
C
	A1UTC = TIME_DIFF(1)
	IF ( ORB_POS_TIME(1) > 1997181 ) A1UTC = 31.0
	IF ( ORB_POS_TIME(1) > 1998365 ) A1UTC = 32.0
	IF ( ORB_POS_TIME(1) > 2005365 ) A1UTC = 33.0
	IF ( ORB_POS_TIME(1) > 2008365 ) A1UTC = 34.0
	IF ( ORB_POS_TIME(1) > 2012182 ) A1UTC = 35.0
	IF ( ORB_POS_TIME(1) > 2015181 ) A1UTC = 36.0
	IF ( ORB_POS_TIME(1) > 2016366 ) A1UTC = 37.0
	YRDOYOLD = ORB_POS_TIME(1)
      ENDIF
C
C      print *, 'timeparts:', ORB_POS_TIME(1), A1UTC
      TIME = DJ + (A1UTC/86400.0D0)
C      print *, 'A1UTC:', A1UTC
C                                                                       
C                                                                       
C     CALL SUBROUTINE IC_FETCH_SLP_COEF TO INSURE APPROPRIATE DATA RECORD
C     IS IN COMMON BLOCK ORSLPR ... AND ALSO TO CALCULATE THE CORRECT
C     TIME INTERVAL WITH WHICH TO EVALUATE THE POLYNOMIAL.
C
C      print *, TIME, TOLD, DABS(TIME-TOLD), EPTIME, LUSLP
      IF(DABS(TIME-TOLD).LE.1.D-6) GO TO 501                            
      CALL IC_FETCH_SLP_COEF(TIME,EPTIME,LUSLP,PV_STATUS)
C      print *, 'pos_vel_of_sun: fetch_slp status:', PV_STATUS
      IF (PV_STATUS.NE.0) GOTO 990
C
C         	LOAD TIME ARRAYS FOR COMPUTATION
C
      TOLD=TIME                                                         
      EPT(1)=EPTIME                                                     
      MAX=19                                                            
      DO 500 I=2,MAX                                                    
      EPT(I)=EPTIME*EPT(I-1)                                            
500   CONTINUE                                                          
501   CONTINUE
C      print *, 'EPTIME etc:', EPTIME, EPT(1)
C
C		BEGIN COMPUTATION OF POLYNOMIAL;  CHECK FOR WHETHER SUN
C           IS LISTED IN "FAST" BODY POSITION OR "SLOW" BODY POSITION
C
      ICHK  =  ICB                                                      
      DO 550 I=1,2                                                      
         IF (ICHK.EQ.NBEPM(1)) GO TO 516                              
         IF (ICHK.EQ.NBEPM(2)) GO TO 518                              
         DO 515 J=3,NBSLP
            IF (ICHK.EQ.NBEPM(J)) GO TO 5210
515         CONTINUE
	 PV_STATUS = 4 
         GO TO 990
516      TEMP(1,I)=0.D0                                          
         TEMP(2,I)=0.D0                                          
         TEMP(3,I)=0.D0                                          
	 TEMPV(1,I)=0.D0
	 TEMPV(2,I)=0.D0
	 TEMPV(3,I)=0.D0
         GO TO 550                                                    
C                                                                       
C     CALCULATE A POSITION VECTOR BY EVALUATING THE 'FAST BODY' POSITION
C     POLYNOMIAL
C                                                                       
518      TEMP(1,I)=PPOLY1(1,1)                                   
         TEMP(2,I)=PPOLY1(2,1)                                   
         TEMP(3,I)=PPOLY1(3,1)                                   
         NDEG2=NDEGRE(2) + 1                                     
         DO 519 K=2,NDEG2                                        
            ETK1=EPT(K-1)                                      
            TEMP(1,I)=ETK1*PPOLY1(1,K)+TEMP(1,I)               
            TEMP(2,I)=ETK1*PPOLY1(2,K)+TEMP(2,I)               
            TEMP(3,I)=ETK1*PPOLY1(3,K)+TEMP(3,I)               
519         CONTINUE                                                
         GO TO 550
5210     IPP= J-2                                                
C                                                                       
C     CALCULATE A POSITION VECTOR BY EVALUATING THE 'SLOW BODY' POSITION
C     POLYNOMIAL                                                        
C                                                                       
         TEMP(1,I)=PPOLY2(1,1,IPP)                               
         TEMP(2,I)=PPOLY2(2,1,IPP)                               
         TEMP(3,I)=PPOLY2(3,1,IPP)                               
         TEMPV(1,I)=PPOLY2(1,2,IPP)                               
         TEMPV(2,I)=PPOLY2(2,2,IPP)                               
         TEMPV(3,I)=PPOLY2(3,2,IPP)                               
         NDEG4=NDEGRE(4) + 1                                     
         DO 522 K=2,NDEG4                                        
            ETK1=EPT(K-1)                                      
            TEMP(1,I)=ETK1*PPOLY2(1,K,IPP)+TEMP(1,I)           
            TEMP(2,I)=ETK1*PPOLY2(2,K,IPP)+TEMP(2,I)           
            TEMP(3,I)=ETK1*PPOLY2(3,K,IPP)+TEMP(3,I)           
	    IF (K .EQ. 2) GOTO 522
	    ETK2=EPT(K-2)*DFLOAT(K-1)
	    TEMPV(1,I)=ETK2*PPOLY2(1,K,IPP)+TEMPV(1,I)
	    TEMPV(2,I)=ETK2*PPOLY2(2,K,IPP)+TEMPV(2,I)
	    TEMPV(3,I)=ETK2*PPOLY2(3,K,IPP)+TEMPV(3,I)
522         CONTINUE                                                
550      ICHK=INCB                                               
C
      SUN_POS_VEL(1)=TEMP(1,2)-TEMP(1,1)
      SUN_POS_VEL(2)=TEMP(2,2)-TEMP(2,1)
      SUN_POS_VEL(3)=TEMP(3,2)-TEMP(3,1)
      SUN_POS_VEL(4)=(TEMPV(1,2)-TEMPV(1,1))/86400.0D0
      SUN_POS_VEL(5)=(TEMPV(2,2)-TEMPV(2,1))/86400.0D0
      SUN_POS_VEL(6)=(TEMPV(3,2)-TEMPV(3,1))/86400.0D0
C
990	CONTINUE
	IF (PV_STATUS .NE. 0) THEN
           print *, 'pos_vel_of_sun error: ', PV_STATUS
           CLOSE (LUCOEF, IOSTAT = STATUS)
	ENDIF
C
991	CONTINUE
	IF (PV_STATUS .NE. 0) THEN
C	   print *, 'Whoa!'
           print *, 'pos_vel_of_sun error: ', PV_STATUS
           CLOSE (LUSLP, IOSTAT = STATUS)
	ENDIF
C
999   CONTINUE
C      ieeer = ieee_flags('clear', 'exception', 'all', out)
      RETURN
C
      END




