OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
l_sun.f
Go to the documentation of this file.
1  subroutine l_sun(iyr,iday,sec,sunr,rs)
2 c
3 c computes unit sun vector in geocentric rotating coodinates, using
4 c subprograms to compute inertial sun vector and greenwich hour angle
5 c
6 c
7 c arguments:
8 c
9 c name Type i/o description
10 c --------------------------------------------------------
11 c iyr i*4 i year, four digits(i.e, 1993)
12 c iday i*4 i day of year(1-366)
13 c sec r*8 i seconds of day
14 c sunr(3) r*4 o unit sun vector in geocentric rotating
15 c coordinates
16 c rs r*4 o earth-to-sun distance(au)
17 c
18 c subprograms referenced:
19 c
20 c sun2000 computes inertial sun vector
21 c gha2000 computes greenwich sidereal angle
22 c
23 c coded by: frederick s. patt, gsc, september 29, 1992
24 c
25 c modification history:
26 c
27 c modifified to use new sun and hour angle routines
28 c frederick s. patt, november 3, 1992
29 c
30 c removed internal jd() function, since it is available as an
31 c independent module. b. a. franz, gsc, november 14, 1997.
32 
33  implicit real*8 (a-h,o-z)
34  real*4 sunr(3),su(3),rs
35  common /gconst/pi,radeg,re,rem,f,omf2,omegae
36 
37 c get unit sun vector in geocentric inertial coordinates
38  call sun2000(iyr,iday,sec,su,rs)
39 
40 c get greenwich mean sideral angle
41  day = iday + sec/864.d2
42  call gha2000(iyr,day,gha)
43  ghar = gha/radeg
44 
45 c transform sun vector into geocentric rotating frame
46  sunr(1) = su(1)*cos(ghar) + su(2)*sin(ghar)
47  sunr(2) = su(2)*cos(ghar) - su(1)*sin(ghar)
48  sunr(3) = su(3)
49 
50  return
51  end
52 
53  subroutine sun2000(iyr,iday,sec,sun,rs)
54 c
55 c this subroutine computes the sun vector in geocentric inertial
56 c(equatorial) coodinates. it uses the model referenced in the
57 c astronomical almanac for 1984, section s(supplement) and documented
58 c in "Exact closed-form geolocation algorithm for Earth survey
59 c sensors", by f.s. patt and w.w. gregg, int. journal of remote
60 c sensing, 1993. the accuracy of the sun vector is approximately 0.1
61 c arcminute.
62 c
63 c arguments:
64 c
65 c name Type i/o description
66 c --------------------------------------------------------
67 c iyr i*4 i year, four digits(i.e, 1993)
68 c iday i*4 i day of year(1-366)
69 c sec r*8 i seconds of day
70 c sun(3) r*4 o unit sun vector in geocentric inertial
71 c coordinates of date
72 c rs r*4 o magnitude of the sun vector(au)
73 c
74 c subprograms referenced:
75 c
76 c jd computes julian day from calendar date
77 c ephparms computes mean solar longitude and anomaly and
78 c mean lunar lontitude and ascending node
79 c nutate compute nutation corrections to lontitude and
80 c obliquity
81 c
82 c coded by: frederick s. patt, gsc, november 2, 1992
83 c modified to include earth constants subroutine by w. gregg,
84 c may 11, 1993.
85 
86 
87  implicit real*8 (a-h,o-z)
88  real*4 sun(3),rs
89  common /nutcm/dpsi,eps,nutime
90  common /gconst/pi,radeg,re,rem,f,omf2,omegae
91  data xk/0.0056932/ !Constant of aberration
92  data imon/1/
93 
94 c compute floating point days since jan 1.5, 2000
95 c note that the julian day starts at noon on the specified date
96  t = jd(iyr,imon,iday) - 2451545.0d0 + (sec-43200.d0)/86400.d0
97 
98 c compute solar ephemeris parameters
99  call ephparms(t,xls,gs,xlm,omega)
100 
101 c check if need to compute nutation corrections for this day
102  nt = t
103  if (nt.ne.nutime) then
104  nutime = nt
105  call nutate(t,xls,gs,xlm,omega,dpsi,eps)
106  end if
107 
108 c compute planet mean anomalies
109 c venus mean anomaly
110  g2 = 50.40828d0 + 1.60213022d0*t
111  g2 = dmod(g2,360.d0)
112 
113 c mars mean anomaly
114  g4 = 19.38816d0 + 0.52402078d0*t
115  g4 = dmod(g4,360.d0)
116 
117 c jupiter mean anomaly
118  g5 = 20.35116d0 + 0.08309121d0*t
119  g5 = dmod(g5,360.d0)
120 
121 c compute solar distance(au)
122  rs = 1.00014d0 - 0.01671d0*cos(gs/radeg)
123  * - 0.00014d0*cos(2.0d0*gs/radeg)
124 
125 c compute geometric solar longitude
126  dls = (6893.0d0 - 4.6543463d-4*t)*sin(gs/radeg)
127  * + 72.0d0*sin(2.0d0*gs/radeg)
128  * - 7.0d0*cos((gs - g5)/radeg)
129  * + 6.0d0*sin((xlm - xls)/radeg)
130  * + 5.0d0*sin((4.0d0*gs - 8.0d0*g4 + 3.0d0*g5)/radeg)
131  * - 5.0d0*cos((2.0d0*gs - 2.0d0*g2)/radeg)
132  * - 4.0d0*sin((gs - g2)/radeg)
133  * + 4.0d0*cos((4.0d0*gs - 8.0d0*g4 + 3.0d0*g5)/radeg)
134  * + 3.0d0*sin((2.0d0*gs - 2.0d0*g2)/radeg)
135  * - 3.0d0*sin(g5/radeg)
136  * - 3.0d0*sin((2.0d0*gs - 2.0d0*g5)/radeg) !arcseconds
137 
138  xlsg = xls + dls/3600.d0
139 
140 c compute apparent solar longitude; includes corrections for nutation
141 c in longitude and velocity aberration
142  xlsa = xlsg + dpsi - xk/rs
143 
144 c compute unit sun vector
145  sun(1) = cos(xlsa/radeg)
146  sun(2) = sin(xlsa/radeg)*cos(eps/radeg)
147  sun(3) = sin(xlsa/radeg)*sin(eps/radeg)
148 c print *,' Sunlon = ',xlsg,xlsa,eps
149 
150  return
151  end
@ floating
float * vector(long nl, long nh)
Definition: nrutil.c:15
short int ephemeris(double tjd, body *cel_obj, short int origin, double *pos, double *vel)
float mean(float *xs, int sample_size)
Definition: numerical.c:81
Definition: l_sun.py:1
subroutine earth(pos, vel, widphse1, widphfl1, widphse2,
Definition: earth.f:2
#define real
Definition: DbAlgOcean.cpp:26
#define re
Definition: l1_czcs_hdf.c:701
Definition: jd.py:1
===========================================================================V5.0.48(Terra) 03/20/2015 Changes shown below are differences from MOD_PR02 V5.0.46(Terra)============================================================================Changes noted for V6.1.20(Terra) below were also instituted for this version.============================================================================V6.1.20(Terra) 03/12/2015 Changes shown below are differences from MOD_PR02 V6.1.18(Terra)============================================================================Changes from v6.1.18 which may affect scientific output:A situation can occur in which a scan which contains sector rotated data has a telemetry value indicating the completeness of the sector rotation. This issue is caused by the timing of the instrument command to perform the sector rotation and the recording of the telemetry point that reports the status of sector rotation. In this case a scan is considered valid by L1B and pass through the calibration - reporting extremely high radiances. Operationally the TEB calibration uses a 40 scan average coefficient, so the 20 scans(one mirror side) after the sector rotation are contaminated with anomalously high radiance values. A similar timing issue appeared before the sector rotation was fixed in V6.1.2. Our analysis indicates the ‘SET_FR_ENC_DELTA’ telemetry correlates well with the sector rotation encoder position. The use of this telemetry point to determine scans that are sector rotated should fix the anomaly occured before and after the sector rotation(usually due to the lunar roll maneuver). The fix related to the sector rotation in V6.1.2 is removed in this version.============================================================================V6.1.18(Terra) 10/01/2014 Changes shown below are differences from MOD_PR02 V6.1.16(Terra)============================================================================Added doi attributes to NRT(Near-Real-Time) product.============================================================================V6.1.16(Terra) 01/27/2014 Changes shown below are differences from MOD_PR02 V6.1.14(Terra)============================================================================Migrate to SDP Toolkit 5.2.17============================================================================V6.1.14(Terra) 06/26/2012 Changes shown below are differences from MOD_PR02 V6.1.12(Terra)============================================================================Added the doi metadata to L1B product============================================================================V6.1.12(Terra) 04/25/2011 Changes shown below are differences from MOD_PR02 V6.1.8(Terra)============================================================================1. The algorithm to calculate uncertainties for reflective solar bands(RSB) is updated. The current uncertainty in L1B code includes 9 terms from prelaunch analysis. The new algorithm regroups them with the new added contributions into 5 terms:u1:the common term(AOI and time independent) and
Definition: HISTORY.txt:126
short int nutate(double tjd, short int fn1, double *pos, double *pos2)
short int aberration(double *pos, double *vel, double lighttime, double *pos2)
#define pi
Definition: vincenty.c:23
int sun2000(size_t sdim, int32_t iyr, int32_t idy, double *sec, orb_array *sun)
#define omf2
Definition: l1_czcs_hdf.c:703
#define f
Definition: l1_czcs_hdf.c:702
for(i=0;i< NROOTS;i++) s[i]
Definition: decode_rs.h:85
integer function julian(DAY, MONTH, YEAR)
Definition: julian.f:2
this program makes no use of any feature of the SDP Toolkit that could generate such a then geolocation is calculated at that and then aggregated up to Resolved feature request Bug by adding three new int8 SDSs for each high resolution offsets between the high resolution geolocation and a bi linear interpolation extrapolation of the positions This can be used to reconstruct the high resolution geolocation Resolved Bug by delaying cumulation of gflags until after validation of derived products Resolved Bug by setting Latitude and Longitude to the correct fill resolving to support Near Real Time because they may be unnecessary if use of entrained ephemeris and attitude data is turned on(as it will be in Near-Real-Time processing).
int ephparms(double t, double &xls, double &gs, double &xlm, double &omega)