WSJT-X/libm65/sun.f
Joe Taylor 9e73f87218 Major changes to the MAP65 branch. This branch is now a copy of
MAP65 v2.3.0, r631, as checked out from the SVN repository on
pulsar.princeton.edu.  If all goes well with this commit, subsequent
MAP65 development will use the Berlios repository.


git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/map65@2461 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
2012-05-22 14:28:39 +00:00

89 lines
3.0 KiB
Fortran

subroutine sun(y,m,DD,UT,lon,lat,RA,Dec,LST,Az,El,mjd,day)
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
C Ecliptic coords of sun (geocentric)
real xs,ys
C Equatorial coords of sun (geocentric)
real xe,ye,ze
real lon,lat
real GMST0,LST,HA
real xx,yy,zz
real xhor,yhor,zhor
real Az,El
real day
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)
day=d-1.5
return
end