;+NAME/ONE LINE DESCRIPTION OF ROUTINE: ; LONLAT2SOLZ calculates the solar and sensor view angles from geodetic ; longitude and latitude and observation time. ; ; ----------------------------------------------------------------- ; pro lonlat2solz ; ; Calculates solar and sensor view angles from geodetic longitude ; and latitude, and observation time. ; ; Inputs: ; lon - longitude of pixel in degrees ; lat - longitude of pixel in degrees ; year - observation year ; day - observation day of year ; sec - observation seconds of day ; ; Outputs: ; solz - solar zenith angle of pixel in degrees ; sola - solar azimuth angle of pixel in degrees ; ; Notes: ; Inputs can be scalar or vector, but vectors must be of ; equal length. ; ; Written By: ; Bryan Franz, SAIC GSC, NASA/SIMBIOS Project, April 1998. ; (with much help from geonav.pro by Fred Patt) ; ; ----------------------------------------------------------------- pro lonlat2solz,lon,lat,year,day,sec,solz,sola ;lon=[-120,0,120] ;lat=[0,0,0] ;year=[2008,2008,2008] ;day=[180,180,180] ;for msec=0UL,3600000*24,3600000 do begin Re = 6378.137 ; Earth radius in km f = 1/298.257 ; Earth flattening factor ;l_sun,year,day,msec/1000.D0,usun l_sun,year,day,double(sec),usun n = n_elements(lon) solz = fltarr(n) sola = fltarr(n) rmat = fltarr(3,3) ; Rotation matrix for i=0L,n-1 do begin rlon = lon(i)*!pi/180. rlat = lat(i)*!pi/180. ; ; First, we must define the axes (in ECR space) of a ; pixel-local coordinate system, with the z-axis along ; the geodetic pixel normal, x-axis pointing east, and ; y-axis pointing north. ; up = [cos(rlat)*cos(rlon),cos(rlat)*sin(rlon),sin(rlat)] upxy = sqrt(up(0)*up(0)+up(1)*up(1)) ea = [-up(1)/upxy,up(0)/upxy,0.0] no = crossp(up,ea) ; ; Compute geocentric pixel location vector. ; phi = atan(tan(rlat)*(1-f)^2) ; geocentric latitude R = Re*(1-f)/sqrt(1-(2-f)*f*(cos(phi)^2)) ; dist to Earth surface gvec = R*[cos(phi)*cos(rlon),cos(phi)*sin(rlon),sin(phi)] ; ; Now we can transform Sun vectors into the local frame. ; rmat(0,*) = ea rmat(1,*) = no rmat(2,*) = up sunl = rmat # usun(*,i) ; ; Compute the solar zenith and azimuth ; solz(i) = atan(sqrt(sunl(0)*sunl(0)+sunl(1)*sunl(1)),sunl(2)) * !radeg if (solz(i) gt 0.05) then begin sola(i) = atan(sunl(0),sunl(1)) * !radeg endif else begin sola(i) = 0.0 endelse if (sola(i) lt 0.0) then $ sola(i) = sola(i) + 360.0d0 endfor if (n eq 1) then begin solz = solz(0) sola = sola(0) endif return ;print,'solz',msec/3600000,solz ;print,'sola',sola ;endfor end ;+NAME/ONE LINE DESCRIPTION OF ROUTINE: ; L_SUN computes unit sun vector in geocentric rotating coordinates. ; ; Computes unit Sun vector in geocentric rotating coodinates, using ; subprograms to compute inertial Sun vector and Greenwich hour angle pro l_sun,iyr,iday,sec,sunr,rs ; implicit real*8 (a-h,o-z) ; real*4 sunr(3),su(3),rs ; radeg = 180.d0/3.14159265359d0 ; Get unit Sun vector in geocentric inertial coordinates sun2000,iyr,iday,sec,su,rs ; Get Greenwich mean sideral angle day = iday + sec/864.d2 gha2000,iyr,day,gha ghar = gha/!radeg ; Transform Sun vector into geocentric rotating frame n = n_elements(day) sunr = fltarr(3,n) sunr(0,*) = su(0,*)*cos(ghar) + su(1,*)*sin(ghar) sunr(1,*) = su(1,*)*cos(ghar) - su(0,*)*sin(ghar) sunr(2,*) = su(2,*) return end ;+NAME/ONE LINE DESCRIPTION OF ROUTINE: ; SUN2000 computes the sun vector in geocentric inertial (equatorial) ; coordinates. ; ; This subroutine computes the Sun vector in geocentric inertial ; (equatorial) coodinates. It uses the model referenced in The ; Astronomical Almanac for 1984, Section S (Supplement) and documented ; for the SeaWiFS Project in "Constants and Parameters for SeaWiFS ; Mission Operations", in TBD. The accuracy of the Sun vector is ; approximately 0.1 arcminute. ; pro sun2000,iyr,iday,sec,sun,rs ; implicit real*8 (a-h,o-z) ; real*4 sun(3),rs xk=0.0056932 ;Constant of aberration imon=1 ; common nutcm,dpsi,eps,nutime ; radeg = 180.d0/3.14159265359d0 ; Compute floating point days since Jan 1.5, 2000 ; Note that the Julian day starts at noon on the specified date jd,iyr,imon,iday,jday t = jday - 2451545.0d0 + (sec-43200.d0)/86400.d0 n = n_elements(t) sun=fltarr(3,n) ; Compute solar ephemeris parameters ephparms,t,xls,gs,xlm,omega ; Check if need to compute nutation corrections for this day ; nt = t ; if (nt ne nutime) then begin ; nutime = nt nutate,t,xls,gs,xlm,omega,dpsi,eps ; endif ; Compute planet mean anomalies ; Venus Mean Anomaly g2 = 50.40828 + 1.60213022*t g2 = g2 mod 360.d0 ; Mars Mean Anomaly g4 = 19.38816 + 0.52402078*t g4 = g4 mod 360.d0 ; Jupiter Mean Anomaly g5 = 20.35116 + 0.08309121*t g5 = g5 mod 360.d0 ; Compute solar distance (AU) rs = 1.00014-0.01671*cos(gs/!radeg)-0.00014*cos(2.*gs/!radeg) ; Compute Geometric Solar Longitude dls = (6893. - 4.6543463D-4*t)*sin(gs/!radeg) $ + 72.*sin(2.*gs/!radeg) $ - 7.*cos((gs - g5)/!radeg) $ + 6.*sin((xlm - xls)/!radeg) $ + 5.*sin((4.*gs - 8.*g4 + 3.*g5)/!radeg) $ - 5.*cos((2.*gs - 2.*g2)/!radeg) $ - 4.*sin((gs - g2)/!radeg) $ + 4.*cos((4.*gs - 8.*g4 + 3.*g5)/!radeg) $ + 3.*sin((2.*gs - 2.*g2)/!radeg) $ - 3.*sin(g5/!radeg) $ - 3.*sin((2.*gs - 2.*g5)/!radeg) xlsg = xls + dls/3600.d0 ; Compute Apparent Solar Longitude; includes corrections for nutation ; in longitude and velocity aberration xlsa = xlsg + dpsi - xk/rs ; Compute unit Sun vector sun(0,*) = cos(xlsa/!radeg) sun(1,*) = sin(xlsa/!radeg)*cos(eps/!radeg) sun(2,*) = sin(xlsa/!radeg)*sin(eps/!radeg) return end ;+NAME/ONE LINE DESCRIPTION OF ROUTINE: ; JD converts a calendar date to a julian day from noon on the calendar day. ; ; This function converts a calendar date to the corresponding Julian ; day starting at noon on the calendar date. The algorithm used is ; from Van Flandern and Pulkkinen, Ap. J. Supplement Series 41, ; November 1979, p. 400. ; ; ; ; Arguments ; ; Name Type I/O Description ; ---- ---- --- ----------- ; i I*4 I Year - e.g. 1970 ; j I*4 I Month - (1-12) ; k I*4 I Day - (1-31) ; jday I*4 O Julian day ; ; external references ; ------------------- ; none ; ; ; Written by Frederick S. Patt, GSC, November 4, 1992 ; ; pro jd,i,j,k,jday i = long(i) j = long(j) k = long(k) jday = 367*i - 7*(i+(j+9)/12)/4 + 275*j/9 + k + 1721014 ; This additional calculation is needed only for dates outside of the ; period March 1, 1900 to February 28, 2100 ; jday = jday + 15 - 3*((i+(j-9)/7)/100+1)/4 return end ;+NAME/ONE LINE DESCRIPTION OF ROUTINE: ; GHA2000 computes the greenwich hour angle in degrees for the input time. ; ; This subroutine computes the Greenwich hour angle in degrees for the ; input time. It uses the model referenced in The Astronomical Almanac ; for 1984, Section S (Supplement) and documented for the SeaWiFS ; Project in "Constants and Parameters for SeaWiFS Mission Operations", ; in TBD. It includes the correction to mean sideral time for nutation ; as well as precession. ; pro gha2000,iyr,day,gha ; implicit real*8 (a-h,o-z) imon=1 ; common nutcm,dpsi,eps,nutime ; radeg = 180.d0/3.14159265359d0 ; Compute days since J2000 iday = long(day) fday = day - iday jday,iyr,iday,jd t = jd - 2451545.5d0 + fday ; Compute Greenwich Mean Sidereal Time (degrees) gmst = 100.4606184d0 + 0.9856473663d0*t + 2.908d-13*t*t ; Check if need to compute nutation correction for this day ; nt = t ; if (nt ne nutime) then begin ; nutime = nt ephparms,t,xls,gs,xlm,omega nutate,t,xls,gs,xlm,omega,dpsi,eps ; endif ; Include apparent time correction and time-of-day gha = gmst + dpsi*cos(eps/!radeg) + fday*360.d0 gha = gha mod 360.d0 neg = where(gha lt 0.d0) if (neg(0) gt -1) then gha(neg) = gha(neg) + 360.d0 ; return end pro ephparms,t,xls,gs,xlm,omega ; This subroutine computes ephemeris parameters used by other Mission ; Operations routines: the solar mean longitude and mean anomaly, and ; the lunar mean longitude and mean ascending node. It uses the model ; referenced in The Astronomical Almanac for 1984, Section S ; (Supplement) and documented for the SeaWiFS Project in "Constants ; and Parameters for SeaWiFS Mission Operations", in TBD. These ; parameters are used to compute the solar longitude and the nutation ; in longitude and obliquity. ; implicit real*8 (a-h,o-z) ; radeg = 180.d0/3.14159265359d0 ; Sun Mean Longitude xls = 280.46592d0 + 0.9856473516d0*t xls = xls mod 360.d0 ; Sun Mean Anomaly gs = 357.52772d0 + 0.9856002831d0*t gs = gs mod 360.d0 ; Moon Mean Longitude xlm = 218.31643d0 + 13.17639648d0*t xlm = xlm mod 360.d0 ; Ascending Node of Moon's Mean Orbit omega = 125.04452d0 - 0.0529537648d0*t omega = omega mod 360.d0 return end pro nutate,t,xls,gs,xlm,omega,dpsi,eps ; This subroutine computes the nutation in longitude and the obliquity ; of the ecliptic corrected for nutation. It uses the model referenced ; in The Astronomical Almanac for 1984, Section S (Supplement) and ; documented for the SeaWiFS Project in "Constants and Parameters for ; SeaWiFS Mission Operations", in TBD. These parameters are used to ; compute the apparent time correction to the Greenwich Hour Angle and ; for the calculation of the geocentric Sun vector. The input ; ephemeris parameters are computed using subroutine ephparms. Terms ; are included to 0.1 arcsecond. ; ; implicit real*8 (a-h,o-z) ; radeg = 180.d0/3.14159265359d0 ; Nutation in Longitude dpsi = - 17.1996*sin(omega/!radeg) $ + 0.2062*sin(2.*omega/!radeg) $ - 1.3187*sin(2.*xls/!radeg) $ + 0.1426*sin(gs/!radeg) $ - 0.2274*sin(2.*xlm/!radeg) ; Mean Obliquity of the Ecliptic epsm = 23.439291d0 - 3.560d-7*t ; Nutation in Obliquity deps = 9.2025*cos(omega/!radeg) + 0.5736*cos(2.*xls/!radeg) ; True Obliquity of the Ecliptic eps = epsm + deps/3600.d0 dpsi = dpsi/3600.d0 return end ;+NAME/ONE LINE DESCRIPTION OF ROUTINE: ; JDAY returns the julian day. ; pro jday,i,k,jd j = 1 l = long((j-14)/12) jd = k - long(32075) + long(1461)*long(i+4800+l)/4 jd = jd + long(367)*long(j-2-l*12)/12 - 3*((i+4900+l)/100)/4 return end