WSJT-X/lib/sun.f90

89 lines
2.8 KiB
Fortran
Raw Normal View History

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 !UT!in hours
real RA,Dec !RA and Dec of sun
! NB: Double caps here are single caps in the writeup.
! 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
! 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
!Ecliptic coords of sun (geocentric)
real xs,ys
!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/
! 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
! 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)
! Ecliptic coordinates of sun (rectangular):
xs = r * cos(lonsun/rad)
ys = r * sin(lonsun/rad)
! Equatorial coordinates of sun (rectangular):
xe = xs
ye = ys * cos(ecl/rad)
ze = ys * sin(ecl/rad)
! 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 subroutine sun