gha2000.f

Go to the documentation of this file.

        subroutine gha2000(iyr,day,gha)
 
c  this subroutine computes the greenwich hour angle in degrees for the
c  input time.  it uses the model referenced in the astronomical almanac
c  for 1984, section s(supplement) and documented in "Exact 
c  closed-form geolocation algorithm for Earth survey sensors", by 
c  f.s. patt and w.w. gregg, int. journal of remote sensing, 1993.
c  it includes the correction to mean sideral time for nutation
c  as well as precession.
 
c  calling arguments
 
c  name         Type    i/o     description
c
c  iyr          i*4      i      year(four digits)
c  day          r*8      i      day(time of day as fraction)
c  gha          r*8      o      greenwich hour angle(degrees)
 
 
c       subprograms referenced:
c
c       jd              computes julian day from calendar date
c       ephparms        computes mean solar longitude and anomaly and
c                        mean lunar lontitude and ascending node
c       nutate          compute nutation corrections to lontitude and 
c                        obliquity
c       
c
c       Program written by:     frederick s. patt
c                               general sciences corporation
c                               november 2, 1992
c
c       modification history:
c
        implicit real*8 (a-h,o-z)
        common /nutcm/dpsi,eps,nutime
        common /gconst/pi,radeg,re,rem,f,omf2,omegae
        data imon/1/,nutime/-9999/
 
c  compute days since j2000
        iday = day
        fday = day - iday
        jday = jd(iyr,imon,iday)
        t = jday - 2451545.5d0 + fday
        
c  compute greenwich mean sidereal time(degrees)
        gmst = 100.4606184d0 + 0.9856473663d0*t + 2.908d-13*t*t
 
c  check if need to compute nutation correction for this day
        nt = t
        if (nt.ne.nutime) then
          nutime = nt
          call ephparms(t,xls,gs,xlm,omega)
          call nutate(t,xls,gs,xlm,omega,dpsi,eps)
        end if
 
c  include apparent time correction and time-of-day
        gha = gmst + dpsi*cos(eps/radeg) + fday*360.d0
        gha = dmod(gha,360.d0)
        if (gha.lt.0.d0) gha = gha + 360.d0
 
        return
        end
 
        subroutine ephparms(t,xls,gs,xlm,omega)
 
c  this subroutine computes ephemeris parameters used by other mission
c  operations routines:  the solar mean longitude and mean anomaly, and
c  the lunar mean longitude and mean ascending node.  it uses the model
c  referenced in the astronomical almanac for 1984, section s 
c(supplement) and documented and documented in "Exact closed-form 
c  geolocation algorithm for Earth survey sensors", by f.s. patt and 
c  w.w. gregg, int. journal of remote sensing, 1993.  these parameters 
c  are used to compute the solar longitude and the nutation in 
c  longitude and obliquity.
 
c  calling arguments
 
c  name         Type    i/o     description
c
c  t            r*8      i      time in days since january 1, 2000 at 
c                                12 hours ut
c  xls          r*8      o      mean solar longitude(degrees)
c  gs           r*8      o      mean solar anomaly(degrees)
c  xlm          r*8      o      mean lunar longitude(degrees)
c  omega        r*8      o      ascending node of mean lunar orbit 
c(degrees)
c
c
c       Program written by:     frederick s. patt
c                               general sciences corporation
c                               november 2, 1992
c
c       modification history:
c
        implicit real*8 (a-h,o-z)
 
c  sun mean longitude           
        xls = 280.46592d0 + 0.9856473516d0*t
        xls = dmod(xls,360.d0)
 
c  sun mean anomaly             
        gs = 357.52772d0 + 0.9856002831d0*t 
        gs = dmod(gs,360.d0)
 
c  moon mean longitude          
        xlm = 218.31643d0 + 13.17639648d0*t 
        xlm = dmod(xlm,360.d0)
 
c  ascending node of moon's Mean Orbit  
        omega = 125.04452d0 - 0.0529537648d0*t 
        omega = dmod(omega,360.d0)
        
        return
        end
 
        subroutine nutate(t,xls,gs,xlm,omega,dpsi,eps)
 
c  This subroutine computes the nutation in longitude and the obliquity
c  of the ecliptic corrected for nutation.  It uses the model referenced
c  in The Astronomical Almanac for 1984, Section S (Supplement) and 
c  documented in "Exact closed-form geolocation algorithm for Earth 
c  survey sensors", by F.S. Patt and W.W. Gregg, Int. Journal of 
c  Remote Sensing, 1993.  These parameters are used to compute the 
c  apparent time correction to the Greenwich Hour Angle and for the 
c  calculation of the geocentric Sun vector.  The input ephemeris 
c  parameters are computed using subroutine ephparms.  Terms are 
c  included to 0.1 arcsecond.
 
c  Calling Arguments
 
c  Name         Type    I/O     Description
c
c  t            R*8      I      Time in days since January 1, 2000 at 
c                                12 hours UT
c  xls          R*8      I      Mean solar longitude (degrees)
c  gs           R*8      I      Mean solar anomaly   (degrees)
c  xlm          R*8      I      Mean lunar longitude (degrees)
c  Omega        R*8      I      Ascending node of mean lunar orbit 
c                                (degrees)
c  dPsi         R*8      O      Nutation in longitude (degrees)
c  Eps          R*8      O      Obliquity of the Ecliptic (degrees)
c                                (includes nutation in obliquity)
c
c
c       Program written by:     Frederick S. Patt
c                               General Sciences Corporation
c                               October 21, 1992
c
c       Modification History:
c
        implicit real*8 (a-h,o-z)
        common /gconst/pi,radeg,re,rem,f,omf2,omegae
        
c  Nutation in Longitude
        dpsi = - 17.1996D0*sin(omega/radeg) 
     *          + 0.2062D0*sin(2.0D0*omega/radeg)
     *          - 1.3187D0*sin(2.0D0*xls/radeg) 
     *          + 0.1426D0*sin(gs/radeg) 
     *          - 0.2274D0*sin(2.0D0*xlm/radeg) 
        
c  Mean Obliquity of the Ecliptic       
        epsm = 23.439291d0 - 3.560d-7*t 
 
c  Nutation in Obliquity 
        deps = 9.2025D0*cos(omega/radeg) + 0.5736D0*cos(2.0D0*xls/radeg)
 
c  True Obliquity of the Ecliptic 
        eps = epsm + deps/3600.d0
 
        dpsi = dpsi/3600.d0
 
        return
        end