;
; IC_GCI_TRANSF
;
; Routines to transform from GCI coordinates to GEO coordinates.
; The following ISTP ICSS routines are appended in this file.
;
; IC_CONV_MATRIX.FOR
; IC_GCI_TO_GEO.FOR
; IC_GET_NUT_ANGLES.FOR
; IC_GRNWCH_SIDEREAL.FOR
;
; In addition, the NAG routine F01CKF is emulated by a subroutine
; included in this file. Consequently, this code should not be
; linked with the NAG libraries if available. The IC_ routines
; have not been modified in any fashion.
;
; A good reference for the routines below is 'An Explanatory
; Supplement to the Astronomical Almanac,' P. Kenneth Seidelmann, ed.
; University Science Books, 1992. ISBN 0-935702-68-7
;
; G. GERMANY 8/9/95
;
;------------------------------------------------------------------
;
; IC_GCI_TO_GEO - return a transformation matrix
;
; PURPOSE: Calculate the transformation matrix from GCI
; coordinates to GEO coordinates at a given date and time.
;
; UNIT TYPE: SUBROUTINE
;
; INVOCATION METHOD: CALL IC_GCI_TO_GEO (orb_pos_time,
; transform_matrix)
;
; argUMENT LIST:
;
; NAME TYPE USE DESCRIPTION
; ---- ---- --- -----------
; orb_pos_time(2) I*4 I time OF ORB. VECTOR, year-day-MILLI OF day
; transform_matrix(3,3) R*8 O TRANSFORMATION MATRIX
;
; FILE/RECORD REFERENCES: NONE
;
; EXTERNAL VARIABLES: NONE
;
; EXTERNAL REFERENCES:
; IC_GRNWCH_SIDEREAL - Returns the Greenwich sidereal time in radians
; F01CKF - Multiplies two matrices
; IC_CONV_MATRIX - Returns the conversion matrix to rotate
; from mean of date to true of date
;
; ABNORMAL TERMINATION CONDITIONS, ERROR MESSAGES: NONE
;
; ASSUMPTIONS, CONstrAINTS, REstrICTIONS: NONE
;
; DEVELOPMENT HISTORY
;
; AUTHOR CHANGE ID RELEASE DATE DESCRIPTION OF CHANGE
; ------ --------- ------- ---- ---------------------
; J. LUBELCZYK B1R1 11/21/90 INITIAL PDL
; J. Lubelczyk B1R1 12/10/90 CODING
; J. Lubelczyk 09/18/91 Updated to return true
; of date trans matrix
; J. LUBELCZYK ICCR #83, CCR #'S 130, 137 11/91 B3 update
;
; NOTES:
;
PRO ic_gci_to_geo, orb_pos_time, transform_matrix
mean_matrix = DBLARR(3,3)
transform_matrix = DBLARR(3,3) ;transformation matrix
cmatrix = DBLARR(3,3) ;the conversion matrix to rotate from
;mean of date to true of date
ic_grnwch_sidereal, orb_pos_time, grnwch_sidereal_time
;
; calculate the sin and cos of the greenwich mean sidereal time
;
sin_gst = sin(grnwch_sidereal_time)
cos_gst = cos(grnwch_sidereal_time)
;
; Fill the mean of date transformation matrix using the sin and cos
; of the greenwich mean sidereal time
;
mean_matrix(0,0) = cos_gst
mean_matrix(0,1) = -sin_gst
mean_matrix(0,2) = 0
mean_matrix(1,0) = sin_gst
mean_matrix(1,1) = cos_gst
mean_matrix(1,2) = 0
mean_matrix(2,0) = 0
mean_matrix(2,1) = 0
mean_matrix(2,2) = 1
ic_conv_matrix,orb_pos_time, cmatrix
transform_matrix = TRANSPOSE( $
TRANSPOSE(mean_matrix) # TRANSPOSE(cmatrix) )
END
; STATEMENT FUNCTION DEFINITION FOR dxjul -- JULIAN EPHEMERIS
; DATE AT BEGINNING OF year.
;
FUNCTION dxjul,i
RETURN,DOUBLE((-32075+1461*(i+4800-13/12)/4 $
+367*(-1+13/12*12)/12-3 $
*((i+4900-13/12)/100)/4)-0.5 )
END
; IC_CONV_MATRIX - Returns the conversion matrix that is necessary to rotate
; from mean of date to true of date
;
; PURPOSE: THIS SUBROUTINE CALCULATES, THROUGH APPROPRIATE ANALYTIC
; EXPRESSIONS, VALUES FOR THE PRECESSION AND NUTATION
; ANGLES AND THE MATRIX REQUIRED TO ROTATE FROM MEAN OF
; JULIAN 2000 TO TRUE OF DATE.
; NAME TYPE USE DESCRIPTION
; ---- ---- --- -----------
; orb_pos_time(2) I*4 I time OF ORB. VECTOR, year-day-MILLI OF day
; CMATRIX(3,3) R*8 O Matrix to rotate from J2000 to true of date
;
; EXTERNAL REFERENCES:
; F01CKF - NAG routine that multiplies two matrices
; IC_GET_NUT_ANGLES Routine to compute the nutation angles
;
PRO ic_conv_matrix,orb_pos_time, cmatrix
cmatrix = DBLARR(3,3) ;conversion matrix
nutmat = DBLARR(3,3) ;Nutation matrix
premat = DBLARR(3,3) ;precession matrix
Jdj2000 = DOUBLE (2451545.0)
R = DOUBLE (1296000.0)
; Convert the given millisecond of day [orb_pos_time(1)] to second of day.
; Convert the packed form into year and day-of-year
secs = (DOUBLE(orb_pos_time(1)))/DOUBLE(1000.0)
year = orb_pos_time(0)/1000
day = orb_pos_time(0) MOD 1000
;
; Calculate the julian date and the time in Julian centuries from J2000
;
fday = secs/DOUBLE(86400.00)
jul_day = dxjul(year) + DOUBLE(day)+fday
time = (jul_day - Jdj2000)/DOUBLE(36525.0)
T2 = time*time
T3 = time*T2
; CALCULATE CONVERSION FACTORS: DEGREES TO RADANS (dtr), SECONDS TO
; RADIANS (str)
;
PI= DOUBLE(4.0 * ATAN(1.0))
dtr=PI/DOUBLE(180.0)
str=dtr/DOUBLE(3600.0)
; CALCULATE PRECESSION ANGLES
zeta = DOUBLE( $
0.11180860865024398D-01*time $
+ 0.14635555405334670D-05*T2 $
+ 0.87256766326094274D-07*T3 )
theta = DOUBLE( $
0.97171734551696701D-02*time $
- 0.20684575704538352D-05*T2 $
- 0.20281210721855218D-06*T3 )
zee = DOUBLE( $
0.11180860865024398D-01*time $
+ 0.53071584043698687D-05*T2 $
+ 0.88250634372368822D-07*T3 )
sinzet = sin(zeta)
coszet = cos(zeta)
sinzee = sin(zee)
coszee = cos(zee)
sinthe = sin(theta)
costhe = cos(theta)
;
; COMPUTE THE TRANSFORMATION MATRIX BETWEEN MEAN EQUATOR AND
; EQUINOX OF 1950 AND MEAN EQUATOR AND EQUINOX OF DATE. THIS
; MATRIX IS CALLED premat.
;
premat(0,0) = -sinzet*sinzee + coszet*coszee*costhe
premat(0,1) = coszee*sinzet + sinzee*costhe*coszet
premat(0,2) = sinthe*coszet
premat(1,0) = -sinzee*coszet - coszee*costhe*sinzet
premat(1,1) = coszee*coszet - sinzee*costhe*sinzet
premat(1,2) = -sinthe*sinzet
premat(2,0) = -coszee*sinthe
premat(2,1) = -sinzee*sinthe
premat(2,2) = costhe
; CALCULATE MEAN OBLIQUITY OF DATE (epso). WHERE TIME IS MEASURED IN
; JULIAN CENTURIES FROM 2000.0.
epso=DOUBLE( (1.813E-3*T3-5.9E-4*T2 $
-4.6815E+1*time+8.4381448E+4)*str )
; CALL IC_GET_NUT_ANGLES TO COMPUTE NUTATION IN OBLIQUITY AND LONGITUDE
ic_get_nut_angles,time,deleps,delpsi,eps
cosep=cos(eps)
cosepO=cos(epso)
cospsi=cos(delpsi)
sinep=sin(eps)
sinepO=sin(epso)
sinpsi=sin(delpsi)
nutmat(0,0)=cospsi
nutmat(1,0)=-sinpsi*cosepO
nutmat(2,0)=-sinpsi*sinepO
nutmat(0,1)=sinpsi*cosep
nutmat(1,1)=cospsi*cosep*cosepO+sinep*sinepO
nutmat(2,1)=cospsi*cosep*sinepO-sinep*cosepO
nutmat(0,2)=sinpsi*sinep
nutmat(1,2)=cospsi*sinep*cosepO-cosep*sinepO
nutmat(2,2)=cospsi*sinep*sinepO+cosep*cosepO
; CALCULATE ELEMENTS OF nutmat * premat. THIS MATRIX IS THE
; ANALYTICALLY CALCULATED TRANSFORMATION MATRIX, WHICH WILL
; TRANSFORM THE MEAN EARTH EQUATOR AND EQUINOX OF J2000 INTO
; THE TRUE EARTH EQUATOR AND EXQUINOX OF DATE.
; cmatrix = nutmat # premat
cmatrix = TRANSPOSE(TRANSPOSE(premat) # TRANSPOSE(nutmat))
END
; IC_GET_NUT_ANGLES - Returns angles that are necessary to adjust the
; Greenwich Hour angle to true of date
;
; PURPOSE : THIS SUBROUTINE CALCULATES, THROUGH APPROPRIATE ANALYTIC
; EXPRESSIONS, VALUES FOR THE NUTATION
; ANGLES TO ROTATE FROM MEAN OF
; JULIAN 2000 TO TRUE OF DATE.
;
; NAME TYPE USE DESCRIPTION
; ---- ---- --- -----------
; time R*8 I time IN JULIAN CENTURIES OF 36525.0
; MEAN SOLAR dayS FROM J2000. (NOTE: THIS
; CAN BE POSITIVE OR NEGATIVE.)
; deleps R*8 O DELTA epsILON, Nutation in obliquity
; delpsi R*8 O DELTA PSI, Nutation in longitude
; eps R*8 O epsILON
PRO ic_get_nut_angles,time,deleps,delpsi,eps
TOL = DOUBLE(0.5)/DOUBLE(36525.0)
isinco = FLTARR(2,106) ;Array used in nutation calculations
icosco = FLTARR(2,106) ;Array used in nutation calculations
fund = DBLARR(1,5) ;The fundamental arguments
T = DBLARR(1,2) ;Time
ifunar = intarr(5,106) ;array used in nutation calculations
oldtim = DOUBLE(0.0)
;
; INITIALIZE VALUES OF IFUNAR, isinco, AND icosco ARRAYS FOR USE IN
; NUTATION CALCULATIONS.
;
; THE 1980 IAU THEORY OF NUTATION,CONTAINED IN JPL
; DE200 PLANETARY EPHEMERIS.
ifunar(0,0:79) = [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, 2,2,2,2,0, 0,0,2,0,0]
ifunar(0,80:105) = [2,2,0,2,2, 2,4,0,2,2, 0,4,2,2,2, 0,-2,2,0,-2, $
2,0,-2,0,2, 0]
ifunar(1,0:79) = [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, 2,2,2,1,1, 2,0,2,2,0]
ifunar(1,80:105) = [2,2,0,2,1, 2,2,0,1,2, 1,2,2,0,1, 1,1,2,0,0, $
1,1,0,0,2, 0]
ifunar(2,0:79) = [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,0,0,0, $
0,0,0,0,0, 0,0,0,0,0, 1,1,-1,0,0, 0,0,0,0,0, $
-1,0,1,0,0, 1,0,-1,-1,0, 0,-1,1,0,0, 0,0,0,0,0]
ifunar(2,80:105) = [0,0,0,1,0, 0,0,-1,0,0, 0,0,0,0,1, -1,0,0,1,0, $
-1,1,0,0,0, 1]
ifunar(3,0:79) = [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]
ifunar(3,80:105) = [-2,-1,2,1,1, -2,-1,1,2,2, 1,0,3,1,0, -1,0,0,0,1, $
0,1,1,2,0, 0]
ifunar(4,0:79) = [0,0,0,0,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,0, 0,2,0,-2,0, 0,0,1,0,-4]
ifunar(4,80:105) = [2,4,-4,-2,2, 4,0,-2,-2,2, 2,-2,-2,-2,0, 2,0,-1,2,-2, $
0,-2,2,2,4, 1]
isinco(0,0:79) = [-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.]
isinco(0,80:105) = [1.,-2.,-1.,1.,-1., -1.,1.,1.,1.,-1., $
-1.,1.,1.,-1.,1., 1.,-1.,-1.,-1.,-1., $
-1.,-1.,-1.,1.,-1., 1.]
isinco(1,0:38) = [-174.2,.2,0.,0.,0., 0.,0.,0.,-1.6,-3.4, $
1.2,-.5,.1,0.,0., -.1,0.,.1,0.,0., $
0.,0.,0.,0.,0., 0.,0.,0.,0.,0., $
-.2,.1,-.4,0.,0., 0.,0.,.1,-.1]
icosco(0,0:79) = [92025.,-895.,-24.,0.,1., 0.,1.,0.,5736.,54., $
224.,-95.,-70.,1.,0., 0.,9.,7.,6.,3., $
3.,-2.,-2.,0.,0., 0.,0.,0.,0.,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.,0.,0.,0., 0.,0.,1.,1.,1., $
1.,1.,-1.,1.,-1., 1.,0.,-1.,-1.,0.]
icosco(0,80:105) = [-1.,1.,0.,-1.,1., 1.,0.,0.,-1.,0., $
0.,0.,0.,0.,0., 0.,0.,0.,0.,0., $
0.,0.,0.,0.,0., 0.]
icosco(1,0:33) = [8.9,.5,0.,0.,0., 0.,0.,0.,-3.1,-.1, $
-.6,.3,0.,0.,0., 0.,0.,0.,0.,0., $
0.,0.,0.,0.,0., 0.,0.,0.,0.,0., $
-.5,0.,0.,-.1]
R = DOUBLE(1296000.0)
IF ( ABS ( time - oldtim ) LE TOL ) THEN BEGIN
deleps = olddep
delpsi = olddps
eps = oldeps
RETURN
ENDIF
T2 = time*time
T3 = time*T2
; CONVERT IFUNAR, isinco, AND icosco ARRAYS TO REAL*8 ARRAYS FUNarg,
; sincof, AND coscof, RESPECTIVELY.
;
funarg = DOUBLE(ifunar)
sincof = DOUBLE(isinco)
coscof = DOUBLE(icosco)
; CALCULATE CONVERSION FACTORS: DEGREES TO RADANS (dtr), SECONDS TO
; RADIANS (str)
PI= DOUBLE(4.0 * ATAN(1.0))
dtr=PI/DOUBLE(180.0)
str=dtr/DOUBLE(3600.0)
; BEGIN COMPUTATION OF NUTATION IN OBLIQUITY AND LONGITUDE
; CALCULATE FUNDAMENTAL argUMENTS FOR USE IN NUTATION CALCULATIONS
; time IS REFERENCED TO J2000.0.
; fund(1,1)= F
; fund(2,1)= OMEGA
; fund(3,1)= L PRIME
; fund(4,1)= L
; fund(5,1)= D
fund(0,0)=str*(335778.877E0+(1342.0E0*R+295263.137E0)*time $
-13.257E0*T2+1.1E-2*T3)
fund(0,1)=str*(450160.280E0-(5.E0*R+482890.539E0)*time+ $
7.455E0*T2+8.0E-3*T3)
fund(0,2)=str*(1287099.804E0+(99.0E0*R+1292581.224E0)*time- $
5.77E-1*T2-1.2E-2*T3)
fund(0,3)=str*(485866.733E0+(1325.0E0*R+715922.633E0)*time+ $
31.310E0*T2+6.4E-2*T3)
fund(0,4)=str*(1072261.307E0+(1236.0E0*R+1105601.328E0)*time- $
6.891E0*T2+1.9E-2*T3)
; CALCULATE MEAN OBLIQUITY OF DATE (epso). WHERE time IS MEASURED IN
; JULIAN CENTURIES FROM 2000.0.
epso=(1.813E-3*T3-5.9E-4*T2-4.6815E+1*time+8.4381448E+4)*str
;
; CALCULATE NUTATION IN LONGITUDE (delpsi) AND NUTATION IN OBLIQUITY
; (deleps). THIS IS A THREE STEP PROCESS:
; (1) CALCULATE argUMENTS OF sinE (FOR delpsi) AND COsinE (FOR deleps)
; THESE ARE OF THE FORM
;
; arg = SUMMATION ( A(I) * fund(I,1) ), I = 1,5
;
; WHERE THE A(I)'S ARE ELEMENTS OF FUNarg.
;
; arg = funarg # fund
; arg = fund # funarg
arg = TRANSPOSE(TRANSPOSE(funarg) # TRANSPOSE(fund))
;
; (2) CALCULATE COEFFICIENTS OF sinE AND COsinE, WHICH ARE THE PRODUCTS
; OF sincof * T AND coscof * T. THESE COEFFICIENTS ARE IN UNITS
; OF 0.0001 SECONDS OF ARC.
;
T(0,0)=DOUBLE(1.0)
T(0,1)=time
; cofcos = coscof # T
; cofsin = sincof # T
; cofcos = T # coscof
; cofsin = T # sincof
cofcos = TRANSPOSE(TRANSPOSE(coscof) # TRANSPOSE(T))
cofsin = TRANSPOSE(TRANSPOSE(sincof) # TRANSPOSE(T))
cofcos=cofcos*DOUBLE(1.E-4)
cofsin=cofsin*DOUBLE(1.E-4)
;
; (3) CALCULATE THE sinES AND COsinES OF THE argUMENTS AND MULTIPLY
; BY THEIR COEFFICIENTS, THEN ADD. COMPUTE delpsi AND deleps.
;
sumpsi=DOUBLE(0.0)
sumeps=DOUBLE(0.0)
sinp=sin(arg)
cose=cos(arg)
FOR E=0,105 DO BEGIN
prodps=cofsin(0,E)*sinp(0,E)
prodep=cofcos(0,E)*cose(0,E)
sumpsi=sumpsi+prodps
sumeps=sumeps+prodep
ENDFOR
deleps=sumeps*str
delpsi=sumpsi*str
;
; CALCULATE TRUE OBLIQUITY OF DATE (eps).
;
eps=epso+deleps
olddep = deleps
olddps = delpsi
oldeps = eps
oldtim = time
END
; IC_GRNWCH_SIDEREAL - return the greenwich true sidereal time in radians
;
; PURPOSE: Calculate the true of date greenwich sidereal time in radians.
;
; NAME TYPE USE DESCRIPTION
; ---- ---- --- -----------
; orb_pos_time(2) I*4 I time OF ORB. VECTOR, year-day-MILLI OF day
; gst R*8 O GREENWICH MEAN SIDEREAL time
; EXTERNAL REFERENCES:
; IC_GET_NUT_ANGLES - Returns angles necessary to adjust the Greenwich
; hour angle to true of date
; NOTES:
; 1) THE ORIGINAL ALGORITHM USED WAS COPIED FROM A SHORT PROGRAM BY
; G. D. MEAD, INCLUDED IN 'GEOPHYSICAL COORDINATE
; TRANSFORMATIONS' BY CHRISTOPHER T. RUSSELL
; 2) THIS VERSION INCORPORATES SEVERAL CHANGES TO CALCULATE THE GREENWICH
; MEAN SIDEREAL time CORRECTLY ON THE J2000 COORDINATE SYSTEM. THE
; PREVIOUS VERSION WAS ONLY CORRECT IN THE B1950 COORDINATE SYSTEM.
; 3) NOW RETURNS THE TRUE OF DATE GREENWICH SIDEREAL time ON THE J2000 SYS
PRO ic_grnwch_sidereal,orb_pos_time, gst
half = DOUBLE(0.50)
C0 = DOUBLE(1.7533685592332653) ;Polynomial Coef.
C1 = DOUBLE(628.33197068884084) ;Polynomial Coef.
C2 = DOUBLE(0.67707139449033354E-05) ;Polynomial Coef.
C3 = DOUBLE(6.3003880989848915) ;Polynomial Coef.
C4 = DOUBLE(-0.45087672343186841E-09) ;Polynomial Coef.
TWOPI = DOUBLE(6.283185307179586) ;Two PI
year = orb_pos_time(0)/1000
day = orb_pos_time(0) MOD 1000
;
; Convert the given millisecond of day [orb_pos_time(1)] to second of day.
;
secs = (DOUBLE(orb_pos_time(1)))/DOUBLE(1000.0)
;
; Begin calculating the greenwich mean sidereal time **
;
fday = secs/86400.00
dj = DOUBLE(365*(year-1900)+(year-1901)/4+day-half)
;
; THE NEXT STATEMENT CAUSES THE REFERENCE EPOCH TO BE SHIFTED
; TO THE J2000 REFERENCE EPOCH.
;
T = (dj-DOUBLE(36525.0))/DOUBLE(36525.0)
gst = DOUBLE(C0 + T*(C1 + T*(C2 + C4*T)) + C3*fday)
gst = DOUBLE(gst MOD TWOPI)
IF (gst LT DOUBLE(0.0)) THEN gst = gst + TWOPI
;
; Convert gst to true of date by adjusting for nutation
;
ic_get_nut_angles, T, deleps, delpsi, eps
gst = gst + delpsi*cos(eps+deleps)
END