pkg/sun/sun_sun2000.F

27f9df093b Oliv*0001 #include "SUN_OPTIONS.h"
                0002
                0003 CBOP
                0004 C     !ROUTINE: SUN_SUN2000
                0005
                0006 C     !INTERFACE: ======================================================
                0007       SUBROUTINE SUN_SUN2000(
                0008      I                        t, xls, gs, xlm, asc, dpsi, eps,
                0009      O                        sunvec, rs )
                0010
                0011 C     !DESCRIPTION:
                0012 C  This subroutine computes the Sun vector in geocentric inertial
                0013 C  (equatorial) coodinates.  It uses the model referenced in The
                0014 C  Astronomical Almanac for 1984, Section S (Supplement) and documented
                0015 C  in Exact closed-form geolocation algorithm for Earth survey
                0016 C  sensors, by F.S. Patt and W.W. Gregg, Int. Journal of Remote
                0017 C  Sensing, 1993.  The accuracy of the Sun vector is approximately 0.1
                0018 C  arcminute.
                0019 C
                0020 C       Coded by:  Frederick S. Patt, GSC, November 2, 1992
                0021 C       Modified to include Earth constants subroutine by W. Gregg,
                0022 C               May 11, 1993.
                0023
                0024 C     !USES: ===========================================================
                0025       IMPLICIT NONE
                0026 #include "SIZE.h"
                0027 #include "EEPARAMS.h"
                0028 #include "PARAMS.h"
                0029
                0030 C     !INPUT PARAMETERS: ===============================================
                0031 C     t    :: Time in days since January 1, 2000 at 12 hours UT
                0032 C     xls  :: Mean solar longitude (degrees)
                0033 C     gs   :: Mean solar anomaly   (degrees)
                0034 C     xlm  :: Mean lunar longitude (degrees)
                0035 C     asc  :: Ascending node of mean lunar orbit
                0036 C     dpsi :: Nutation in longitude (degrees)
                0037 C     eps  :: Obliquity of the Ecliptic (degrees)
                0038       _RL t, xls, gs, xlm, asc, dpsi, eps
                0039
                0040 C     !OUTPUT PARAMETERS: ==============================================
                0041 C     sunvec :: Unit Sun vector in geocentric inertial coords of date
                0042 C     rs     :: Magnitude of the Sun vector (AU)
                0043       _RL sunvec(3), rs
                0044 CEOP
                0045
                0046 #ifdef ALLOW_SUN
                0047
                0048 C     !LOCAL VARIABLES: ================================================
                0049       INTEGER nt
                0050       _RL xk,g2,g4,g5,dls,xlsg,xlsa
                0051 C  Constant of aberration
                0052       PARAMETER (xk=0.0056932 _d 0)
                0053
                0054 C  Compute planet mean anomalies
                0055 C   Venus Mean Anomaly
                0056       g2 = 50.40828 _d 0 + 1.60213022 _d 0*t
086a45f245 Oliv*0057       g2 = MOD(g2,360 _d 0)
27f9df093b Oliv*0058
                0059 C   Mars Mean Anomaly
                0060       g4 = 19.38816 _d 0 + 0.52402078 _d 0*t
086a45f245 Oliv*0061       g4 = MOD(g4,360 _d 0)
27f9df093b Oliv*0062
                0063 C  Jupiter Mean Anomaly
                0064       g5 = 20.35116 _d 0 + 0.08309121 _d 0*t
086a45f245 Oliv*0065       g5 = MOD(g5,360 _d 0)
27f9df093b Oliv*0066
                0067 C  Compute solar distance (AU)
                0068       rs = 1.00014 _d 0 - 0.01671 _d 0*COS(gs*deg2rad)
                0069      &                  - 0.00014 _d 0*COS(2.0 _d 0*gs*deg2rad)
                0070
                0071 C  Compute Geometric Solar Longitude
                0072       dls = (6893.0 _d 0 - 4.6543463 _d -4*t)*SIN(gs*deg2rad)
                0073      & + 72.0 _d 0*SIN(2.0 _d 0*gs*deg2rad)
                0074      & - 7.0 _d 0*COS((gs - g5)*deg2rad)
                0075      & + 6.0 _d 0*SIN((xlm - xls)*deg2rad)
                0076      & + 5.0 _d 0*SIN((4.0 _d 0*gs - 8.0 _d 0*g4 + 3.0 _d 0*g5)*deg2rad)
                0077      & - 5.0 _d 0*COS((2.0 _d 0*gs - 2.0 _d 0*g2)*deg2rad)
                0078      & - 4.0 _d 0*SIN((gs - g2)*deg2rad)
                0079      & + 4.0 _d 0*COS((4.0 _d 0*gs - 8.0 _d 0*g4 + 3.0 _d 0*g5)*deg2rad)
                0080      & + 3.0 _d 0*SIN((2.0 _d 0*gs - 2.0 _d 0*g2)*deg2rad)
                0081      & - 3.0 _d 0*SIN(g5*deg2rad)
                0082      & - 3.0 _d 0*SIN((2.0 _d 0*gs - 2.0 _d 0*g5)*deg2rad)
                0083
                0084       xlsg = xls + dls/3600.0 _d 0
                0085
                0086 C  Compute Apparent Solar Longitude; includes corrections for nutation
                0087 C  in longitude and velocity aberration
                0088       xlsa = xlsg + dpsi - xk/rs
                0089
                0090 C  Compute unit Sun vector
                0091       sunvec(1) = COS(xlsa*deg2rad)
                0092       sunvec(2) = SIN(xlsa*deg2rad)*COS(eps*deg2rad)
                0093       sunvec(3) = SIN(xlsa*deg2rad)*SIN(eps*deg2rad)
                0094
                0095 #endif
                0096
                0097       RETURN
                0098       END