; ; 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