1 |
jahn |
1.1 |
C $Header$ |
2 |
|
|
C $Name$ |
3 |
|
|
|
4 |
|
|
#include "RADTRANS_OPTIONS.h" |
5 |
|
|
|
6 |
|
|
CBOP |
7 |
|
|
C !ROUTINE: RADTRANS_SUN2000 |
8 |
|
|
|
9 |
|
|
C !INTERFACE: ====================================================== |
10 |
|
|
subroutine radtrans_sun2000 (radeg, iyr, imon, iday, sec, |
11 |
|
|
O sunvec, rs) |
12 |
|
|
|
13 |
|
|
C !DESCRIPTION: |
14 |
|
|
c This subroutine computes the Sun vector in geocentric inertial |
15 |
|
|
c (equatorial) coodinates. It uses the model referenced in The |
16 |
|
|
c Astronomical Almanac for 1984, Section S (Supplement) and documented |
17 |
|
|
c in Exact closed-form geolocation algorithm for Earth survey |
18 |
|
|
c sensors, by F.S. Patt and W.W. Gregg, Int. Journal of Remote |
19 |
|
|
c Sensing, 1993. The accuracy of the Sun vector is approximately 0.1 |
20 |
|
|
c arcminute. |
21 |
|
|
c |
22 |
|
|
c Arguments: |
23 |
|
|
c |
24 |
|
|
c Name Type I/O Description |
25 |
|
|
c -------------------------------------------------------- |
26 |
|
|
c IYR I*4 I Year, four digits (i.e, 1993) |
27 |
|
|
c IDAY I*4 I Day of year (1-366) |
28 |
|
|
c SEC R*8 I Seconds of day |
29 |
|
|
c SUN(3) R*8 O Unit Sun vector in geocentric inertial |
30 |
|
|
c coordinates of date |
31 |
|
|
c RS R*8 O Magnitude of the Sun vector (AU) |
32 |
|
|
c |
33 |
|
|
c Subprograms referenced: |
34 |
|
|
c |
35 |
|
|
c JD Computes Julian day from calendar date |
36 |
|
|
c EPHPARMS Computes mean solar longitude and anomaly and |
37 |
|
|
c mean lunar lontitude and ascending node |
38 |
|
|
c NUTATE Compute nutation corrections to lontitude and |
39 |
|
|
c obliquity |
40 |
|
|
c |
41 |
|
|
c Coded by: Frederick S. Patt, GSC, November 2, 1992 |
42 |
|
|
c Modified to include Earth constants subroutine by W. Gregg, |
43 |
|
|
c May 11, 1993. |
44 |
|
|
|
45 |
|
|
C !USES: =========================================================== |
46 |
|
|
IMPLICIT NONE |
47 |
|
|
#include "RADTRANS_VARS.h" |
48 |
|
|
|
49 |
|
|
C !INPUT PARAMETERS: =============================================== |
50 |
|
|
INTEGER iyr, imon, iday |
51 |
|
|
_RL radeg, sec |
52 |
|
|
c INTEGER myThid |
53 |
|
|
|
54 |
|
|
C !OUTPUT PARAMETERS: ============================================== |
55 |
|
|
_RL sunvec(3), rs |
56 |
|
|
|
57 |
|
|
C !FUNCTIONS: ====================================================== |
58 |
|
|
INTEGER radtrans_jd |
59 |
|
|
EXTERNAL radtrans_jd |
60 |
|
|
|
61 |
|
|
C !LOCAL VARIABLES: ================================================ |
62 |
|
|
INTEGER nt |
63 |
|
|
_RL xk,rjd,t,xls,gs,xlm,omega,g2,g4,g5,dls,xlsg,xlsa |
64 |
|
|
parameter (xk=0.0056932) !Constant of aberration |
65 |
|
|
CEOP |
66 |
|
|
|
67 |
|
|
c Compute floating point days since Jan 1.5, 2000 |
68 |
|
|
c Note that the Julian day starts at noon on the specified date |
69 |
|
|
rjd = float(radtrans_jd(iyr,imon,iday)) |
70 |
|
|
t = rjd - 2451545.0D0 + (sec-43200.0D0)/86400.0D0 |
71 |
|
|
|
72 |
|
|
c Compute solar ephemeris parameters |
73 |
|
|
call radtrans_ephparms (t, xls, gs, xlm, omega) |
74 |
|
|
|
75 |
|
|
c Check if need to compute nutation corrections for this day |
76 |
|
|
nt = int(t) |
77 |
|
|
if (nt.ne.nutime) then |
78 |
|
|
nutime = nt |
79 |
|
|
call radtrans_nutate (radeg, t, xls, gs, xlm, omega, dpsi, eps) |
80 |
|
|
end if |
81 |
|
|
|
82 |
|
|
c Compute planet mean anomalies |
83 |
|
|
c Venus Mean Anomaly |
84 |
|
|
g2 = 50.40828D0 + 1.60213022D0*t |
85 |
|
|
g2 = mod(g2,360.0) |
86 |
|
|
|
87 |
|
|
c Mars Mean Anomaly |
88 |
|
|
g4 = 19.38816D0 + 0.52402078D0*t |
89 |
|
|
g4 = mod(g4,360.0) |
90 |
|
|
|
91 |
|
|
c Jupiter Mean Anomaly |
92 |
|
|
g5 = 20.35116D0 + 0.08309121D0*t |
93 |
|
|
g5 = mod(g5,360.0) |
94 |
|
|
|
95 |
|
|
c Compute solar distance (AU) |
96 |
|
|
rs = 1.00014D0 - 0.01671D0*cos(gs/radeg) |
97 |
|
|
* - 0.00014D0*cos(2.0D0*gs/radeg) |
98 |
|
|
|
99 |
|
|
c Compute Geometric Solar Longitude |
100 |
|
|
dls = (6893.0D0 - 4.6543463D-4*t)*sin(gs/radeg) |
101 |
|
|
* + 72.0D0*sin(2.0D0*gs/radeg) |
102 |
|
|
* - 7.0D0*cos((gs - g5)/radeg) |
103 |
|
|
* + 6.0D0*sin((xlm - xls)/radeg) |
104 |
|
|
* + 5.0D0*sin((4.0D0*gs - 8.0D0*g4 + 3.0D0*g5)/radeg) |
105 |
|
|
* - 5.0D0*cos((2.0D0*gs - 2.0D0*g2)/radeg) |
106 |
|
|
* - 4.0D0*sin((gs - g2)/radeg) |
107 |
|
|
* + 4.0D0*cos((4.0D0*gs - 8.0D0*g4 + 3.0D0*g5)/radeg) |
108 |
|
|
* + 3.0D0*sin((2.0D0*gs - 2.0D0*g2)/radeg) |
109 |
|
|
* - 3.0D0*sin(g5/radeg) |
110 |
|
|
* - 3.0D0*sin((2.0D0*gs - 2.0D0*g5)/radeg) !arcseconds |
111 |
|
|
|
112 |
|
|
xlsg = xls + dls/3600.0D0 |
113 |
|
|
|
114 |
|
|
c Compute Apparent Solar Longitude; includes corrections for nutation |
115 |
|
|
c in longitude and velocity aberration |
116 |
|
|
xlsa = xlsg + dpsi - xk/rs |
117 |
|
|
|
118 |
|
|
c Compute unit Sun vector |
119 |
|
|
sunvec(1) = cos(xlsa/radeg) |
120 |
|
|
sunvec(2) = sin(xlsa/radeg)*cos(eps/radeg) |
121 |
|
|
sunvec(3) = sin(xlsa/radeg)*sin(eps/radeg) |
122 |
|
|
c type *,' Sunlon = ',xlsg,xlsa,eps |
123 |
|
|
|
124 |
|
|
return |
125 |
|
|
end |
126 |
|
|
c |