2007-01-11 21:25:52 +00:00
|
|
|
subroutine sun(y,m,DD,UT,lon,lat,RA,Dec,LST,Az,El,mjd)
|
|
|
|
|
|
|
|
|
|
implicit none
|
|
|
|
|
|
|
|
|
|
integer y !Year
|
|
|
|
|
integer m !Month
|
|
|
|
|
integer DD !Day
|
|
|
|
|
integer mjd !Modified Julian Date
|
|
|
|
|
real UT !UTC in hours
|
|
|
|
|
real RA,Dec !RA and Dec of sun
|
|
|
|
|
|
|
|
|
|
C NB: Double caps here are single caps in the writeup.
|
|
|
|
|
|
|
|
|
|
C Orbital elements of the Sun (also N=0, i=0, a=1):
|
|
|
|
|
real w !Argument of perihelion
|
|
|
|
|
real e !Eccentricity
|
|
|
|
|
real MM !Mean anomaly
|
|
|
|
|
real Ls !Mean longitude
|
|
|
|
|
|
|
|
|
|
C Other standard variables:
|
|
|
|
|
real v !True anomaly
|
|
|
|
|
real EE !Eccentric anomaly
|
|
|
|
|
real ecl !Obliquity of the ecliptic
|
|
|
|
|
real d !Ephemeris time argument in days
|
|
|
|
|
real r !Distance to sun, AU
|
|
|
|
|
real xv,yv !x and y coords in ecliptic
|
|
|
|
|
real lonsun !Ecliptic long and lat of sun
|
|
|
|
|
real xs,ys !Ecliptic coords of sun (geocentric)
|
|
|
|
|
real xe,ye,ze !Equatorial coords of sun (geocentric)
|
|
|
|
|
real lon,lat
|
|
|
|
|
real GMST0,LST,HA
|
|
|
|
|
real xx,yy,zz
|
|
|
|
|
real xhor,yhor,zhor
|
|
|
|
|
real Az,El
|
|
|
|
|
|
|
|
|
|
real rad
|
|
|
|
|
data rad/57.2957795/
|
|
|
|
|
|
|
|
|
|
C Time in days, with Jan 0, 2000 equal to 0.0:
|
|
|
|
|
d=367*y - 7*(y+(m+9)/12)/4 + 275*m/9 + DD - 730530 + UT/24.0
|
|
|
|
|
mjd=d + 51543
|
|
|
|
|
ecl = 23.4393 - 3.563e-7 * d
|
|
|
|
|
|
|
|
|
|
C Compute updated orbital elements for Sun:
|
|
|
|
|
w = 282.9404 + 4.70935e-5 * d
|
|
|
|
|
e = 0.016709 - 1.151e-9 * d
|
|
|
|
|
MM = mod(356.0470d0 + 0.9856002585d0 * d + 360000.d0,360.d0)
|
|
|
|
|
Ls = mod(w+MM+720.0,360.0)
|
|
|
|
|
|
|
|
|
|
EE = MM + e*rad*sin(MM/rad) * (1.0 + e*cos(M/rad))
|
|
|
|
|
EE = EE - (EE - e*rad*sin(EE/rad)-MM) / (1.0 - e*cos(EE/rad))
|
|
|
|
|
|
|
|
|
|
xv = cos(EE/rad) - e
|
|
|
|
|
yv = sqrt(1.0-e*e) * sin(EE/rad)
|
|
|
|
|
v = rad*atan2(yv,xv)
|
|
|
|
|
r = sqrt(xv*xv + yv*yv)
|
|
|
|
|
lonsun = mod(v + w + 720.0,360.0)
|
|
|
|
|
C Ecliptic coordinates of sun (rectangular):
|
|
|
|
|
xs = r * cos(lonsun/rad)
|
|
|
|
|
ys = r * sin(lonsun/rad)
|
|
|
|
|
|
|
|
|
|
C Equatorial coordinates of sun (rectangular):
|
|
|
|
|
xe = xs
|
|
|
|
|
ye = ys * cos(ecl/rad)
|
|
|
|
|
ze = ys * sin(ecl/rad)
|
|
|
|
|
|
|
|
|
|
C RA and Dec in degrees:
|
|
|
|
|
RA = rad*atan2(ye,xe)
|
|
|
|
|
Dec = rad*atan2(ze,sqrt(xe*xe + ye*ye))
|
|
|
|
|
|
|
|
|
|
GMST0 = (Ls + 180.0)/15.0
|
|
|
|
|
LST = mod(GMST0+UT+lon/15.0+48.0,24.0) !LST in hours
|
|
|
|
|
HA = 15.0*LST - RA !HA in degrees
|
|
|
|
|
xx = cos(HA/rad)*cos(Dec/rad)
|
|
|
|
|
yy = sin(HA/rad)*cos(Dec/rad)
|
|
|
|
|
zz = sin(Dec/rad)
|
|
|
|
|
xhor = xx*sin(lat/rad) - zz*cos(lat/rad)
|
|
|
|
|
yhor = yy
|
|
|
|
|
zhor = xx*cos(lat/rad) + zz*sin(lat/rad)
|
|
|
|
|
Az = mod(rad*atan2(yhor,xhor) + 180.0 + 360.0,360.0)
|
|
|
|
|
El = rad*asin(zhor)
|
|
|
|
|
|
|
|
|
|
return
|
|
|
|
|
end
|
