w4kek/WSJT-X
1
0
Fork 0
mirror of https://github.com/saitohirga/WSJT-X.git synced 2024-11-28 23:28:49 -05:00
Code Issues Packages Projects Releases Wiki Activity

Another bunch needesd for astro display.

Browse Source 
git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@3839 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
This commit is contained in:
Joe Taylor 2014-03-05 20:14:36 +00:00
parent 2e6611dbb5
commit f7d4a71454
9 changed files with 614 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

108
lib/astro.f90 Normal file
View File

@ -0,0 +1,108 @@
subroutine astro(nyear,month,nday,uth,nfreq,Mygrid, &
NStation,MoonDX,AzSun,ElSun,AzMoon0,ElMoon0, &
ntsky,doppler00,doppler,dbMoon,RAMoon,DecMoon,HA,Dgrd,sd, &
poloffset,xnr,day,lon,lat,LST,techo)
! Computes astronomical quantities for display and tracking.
! NB: may want to smooth the Tsky map to 10 degrees or so.
character*6 MyGrid,HisGrid
real LST
real lat,lon
integer*2 nt144(180)
! common/echo/xdop(2),techo,AzMoon,ElMoon,mjd
real xdop(2)
data rad/57.2957795/
data nt144/ &
234, 246, 257, 267, 275, 280, 283, 286, 291, 298, &
305, 313, 322, 331, 341, 351, 361, 369, 376, 381, &
383, 382, 379, 374, 370, 366, 363, 361, 363, 368, &
376, 388, 401, 415, 428, 440, 453, 467, 487, 512, &
544, 579, 607, 618, 609, 588, 563, 539, 512, 482, &
450, 422, 398, 379, 363, 349, 334, 319, 302, 282, &
262, 242, 226, 213, 205, 200, 198, 197, 196, 197, &
200, 202, 204, 205, 204, 203, 202, 201, 203, 206, &
212, 218, 223, 227, 231, 236, 240, 243, 247, 257, &
276, 301, 324, 339, 346, 344, 339, 331, 323, 316, &
312, 310, 312, 317, 327, 341, 358, 375, 392, 407, &
422, 437, 451, 466, 480, 494, 511, 530, 552, 579, &
612, 653, 702, 768, 863,1008,1232,1557,1966,2385, &
2719,2924,3018,3038,2986,2836,2570,2213,1823,1461, &
1163, 939, 783, 677, 602, 543, 494, 452, 419, 392, &
373, 360, 353, 350, 350, 350, 350, 350, 350, 348, &
344, 337, 329, 319, 307, 295, 284, 276, 272, 272, &
273, 274, 274, 271, 266, 260, 252, 245, 238, 231/
save
call grid2deg(MyGrid,elon,lat)
lon=-elon
call sun(nyear,month,nday,uth,lon,lat,RASun,DecSun,LST, &
AzSun,ElSun,mjd,day)
freq=nfreq*1.e6
if(nfreq.eq.2) freq=1.8e6
if(nfreq.eq.4) freq=3.5e6
call MoonDop(nyear,month,nday,uth,lon,lat,RAMoon,DecMoon, &
LST,HA,AzMoon,ElMoon,vr,dist)
! Compute spatial polarization offset
xx=sin(lat/rad)*cos(ElMoon/rad) - cos(lat/rad)* &
cos(AzMoon/rad)*sin(ElMoon/rad)
yy=cos(lat/rad)*sin(AzMoon/rad)
if(NStation.eq.1) poloffset1=rad*atan2(yy,xx)
if(NStation.eq.2) poloffset2=rad*atan2(yy,xx)
techo=2.0 * dist/2.99792458e5 !Echo delay time
doppler=-freq*vr/2.99792458e5 !One-way Doppler
call coord(0.,0.,-1.570796,1.161639,RAMoon/rad,DecMoon/rad,el,eb)
longecl_half=nint(rad*el/2.0)
if(longecl_half.lt.1 .or. longecl_half.gt.180) longecl_half=180
t144=nt144(longecl_half)
tsky=(t144-2.7)*(144.0/nfreq)**2.6 + 2.7 !Tsky for obs freq
xdop(NStation)=doppler
if(NStation.eq.2) then
HisGrid=MyGrid
go to 900
endif
doppler00=2.0*xdop(1)
doppler=xdop(1)+xdop(2)
! if(mode.eq.3) doppler=2.0*xdop(1)
dBMoon=-40.0*log10(dist/356903.)
sd=16.23*370152.0/dist
! if(NStation.eq.1 .and. MoonDX.ne.0 .and.
! + (mode.eq.2 .or. mode.eq.5)) then
if(NStation.eq.1 .and. MoonDX.ne.0) then
poloffset=mod(poloffset2-poloffset1+720.0,180.0)
if(poloffset.gt.90.0) poloffset=poloffset-180.0
x1=abs(cos(2*poloffset/rad))
if(x1.lt.0.056234) x1=0.056234
xnr=-20.0*log10(x1)
if(HisGrid(1:1).lt.'A' .or. HisGrid(1:1).gt.'R') xnr=0
endif
tr=80.0 !Good preamp
tskymin=13.0*(408.0/nfreq)**2.6 !Cold sky temperature
tsysmin=tskymin+tr
tsys=tsky+tr
dgrd=-10.0*log10(tsys/tsysmin) + dbMoon
900 AzMoon0=Azmoon
ElMoon0=Elmoon
ntsky=nint(tsky)
! auxHA = 15.0*(LST-auxra) !HA in degrees
! pi=3.14159265
! pio2=0.5*pi
! call coord(pi,pio2-lat/rad,0.0,lat/rad,auxha*pi/180.0,
! + auxdec/rad,azaux,elaux)
! AzAux=azaux*rad
! ElAux=ElAux*rad
return
end subroutine astro

82
lib/astro0.f90 Normal file
View File

@ -0,0 +1,82 @@
subroutine astro0(nyear,month,nday,uth8,nfreq,mygrid,hisgrid, &
AzSun8,ElSun8,AzMoon8,ElMoon8,AzMoonB8,ElMoonB8,ntsky,ndop,ndop00, &
dbMoon8,RAMoon8,DecMoon8,HA8,Dgrd8,sd8,poloffset8,xnr8,dfdt,dfdt0, &
width1,width2,w501,w502,xlst8,techo8)
parameter (DEGS=57.2957795130823d0)
character*6 mygrid,hisgrid
real*8 AzSun8,ElSun8,AzMoon8,ElMoon8,AzMoonB8,ElMoonB8
real*8 dbMoon8,RAMoon8,DecMoon8,HA8,Dgrd8,xnr8,dfdt,dfdt0,dt
real*8 sd8,poloffset8,day8,width1,width2,w501,w502,xlst8
real*8 uth8,techo8
data uth8z/0.d0/
save
uth=uth8
call astro(nyear,month,nday,uth,nfreq,hisgrid,2,1, &
AzSun,ElSun,AzMoon,ElMoon,ntsky,doppler00,doppler, &
dbMoon,RAMoon,DecMoon,HA,Dgrd,sd,poloffset,xnr, &
day,xlon2,xlat2,xlst,techo)
AzMoonB8=AzMoon
ElMoonB8=ElMoon
call astro(nyear,month,nday,uth,nfreq,mygrid,1,1, &
AzSun,ElSun,AzMoon,ElMoon,ntsky,doppler00,doppler, &
dbMoon,RAMoon,DecMoon,HA,Dgrd,sd,poloffset,xnr, &
day,xlon1,xlat1,xlst,techo)
day8=day
xlst8=xlst
techo8=techo
call tm2(day8,xlat1,xlon1,xl1,b1)
call tm2(day8,xlat2,xlon2,xl2,b2)
call tm2(day8+1.d0/1440.0,xlat1,xlon1,xl1a,b1a)
call tm2(day8+1.d0/1440.0,xlat2,xlon2,xl2a,b2a)
fghz=0.001*nfreq
dldt1=DEGS*(xl1a-xl1)
dbdt1=DEGS*(b1a-b1)
dldt2=DEGS*(xl2a-xl2)
dbdt2=DEGS*(b2a-b2)
rate1=2.0*sqrt(dldt1**2 + dbdt1**2)
width1=0.5*6741*fghz*rate1
rate2=sqrt((dldt1+dldt2)**2 + (dbdt1+dbdt2)**2)
width2=0.5*6741*fghz*rate2
fbend=0.7
a2=0.0045*log(fghz/fbend)/log(1.05)
if(fghz.lt.fbend) a2=0.0
f50=0.19 * (fghz/fbend)**a2
if(f50.gt.1.0) f50=1.0
w501=f50*width1
w502=f50*width2
AzSun8=AzSun
ElSun8=ElSun
AzMoon8=AzMoon
ElMoon8=ElMoon
dbMoon8=dbMoon
RAMoon8=RAMoon/15.0
DecMoon8=DecMoon
HA8=HA
Dgrd8=Dgrd
sd8=sd
poloffset8=poloffset
xnr8=xnr
ndop=nint(doppler)
ndop00=nint(doppler00)
if(uth8z.eq.0.d0) then
uth8z=uth8-1.d0/3600.d0
dopplerz=doppler
doppler00z=doppler00
endif
dt=60.0*(uth8-uth8z)
if(dt.le.0) dt=1.d0/60.d0
dfdt=(doppler-dopplerz)/dt
dfdt0=(doppler00-doppler00z)/dt
uth8z=uth8
dopplerz=doppler
doppler00z=doppler00
return
end subroutine astro0

41
lib/coord.f90 Normal file
View File

@ -0,0 +1,41 @@
SUBROUTINE COORD(A0,B0,AP,BP,A1,B1,A2,B2)
! Examples:
! 1. From ha,dec to az,el:
! call coord(pi,pio2-lat,0.,lat,ha,dec,az,el)
! 2. From az,el to ha,dec:
! call coord(pi,pio2-lat,0.,lat,az,el,ha,dec)
! 3. From ra,dec to l,b
! call coord(4.635594495,-0.504691042,3.355395488,0.478220215,
! ra,dec,l,b)
! 4. From l,b to ra,dec
! call coord(1.705981071d0,-1.050357016d0,2.146800277d0,
! 0.478220215d0,l,b,ra,dec)
! 5. From ra,dec to ecliptic latitude (eb) and longitude (el):
! call coord(0.d0,0.d0,-pio2,pio2-23.443*pi/180,ra,dec,el,eb)
! 6. From ecliptic latitude (eb) and longitude (el) to ra,dec:
! call coord(0.d0,0.d0,-pio2,pio2-23.443*pi/180,el,eb,ra,dec)
SB0=sin(B0)
CB0=cos(B0)
SBP=sin(BP)
CBP=cos(BP)
SB1=sin(B1)
CB1=cos(B1)
SB2=SBP*SB1 + CBP*CB1*cos(AP-A1)
CB2=SQRT(1.e0-SB2**2)
B2=atan(SB2/CB2)
SAA=sin(AP-A1)*CB1/CB2
CAA=(SB1-SB2*SBP)/(CB2*CBP)
CBB=SB0/CBP
SBB=sin(AP-A0)*CB0
SA2=SAA*CBB-CAA*SBB
CA2=CAA*CBB+SAA*SBB
TA2O2=0.0 !Shut up compiler warnings. -db
IF(CA2.LE.0.e0) TA2O2=(1.e0-CA2)/SA2
IF(CA2.GT.0.e0) TA2O2=SA2/(1.e0+CA2)
A2=2.e0*atan(TA2O2)
IF(A2.LT.0.e0) A2=A2+6.2831853
RETURN
END SUBROUTINE COORD

17
lib/geocentric.f90 Normal file
View File

@ -0,0 +1,17 @@
subroutine geocentric(alat,elev,hlt,erad)
implicit real*8 (a-h,o-z)
! IAU 1976 flattening f, equatorial radius a
f = 1.d0/298.257d0
a = 6378140.d0
c = 1.d0/sqrt(1.d0 + (-2.d0 + f)*f*sin(alat)*sin(alat))
arcf = (a*c + elev)*cos(alat)
arsf = (a*(1.d0 - f)*(1.d0 - f)*c + elev)*sin(alat)
hlt = datan2(arsf,arcf)
erad = sqrt(arcf*arcf + arsf*arsf)
erad = 0.001d0*erad
return
end subroutine geocentric

166
lib/moon2.f90 Normal file
View File

@ -0,0 +1,166 @@
subroutine moon2(y,m,Day,UT,lon,lat,RA,Dec,topRA,topDec, &
LST,HA,Az,El,dist)
implicit none
integer y !Year
integer m !Month
integer Day !Day
real*8 UT !UTC in hours
real*8 RA,Dec !RA and Dec of moon
! NB: Double caps are single caps in the writeup.
real*8 NN !Longitude of ascending node
real*8 i !Inclination to the ecliptic
real*8 w !Argument of perigee
real*8 a !Semi-major axis
real*8 e !Eccentricity
real*8 MM !Mean anomaly
real*8 v !True anomaly
real*8 EE !Eccentric anomaly
real*8 ecl !Obliquity of the ecliptic
real*8 d !Ephemeris time argument in days
real*8 r !Distance to sun, AU
real*8 xv,yv !x and y coords in ecliptic
real*8 lonecl,latecl !Ecliptic long and lat of moon
real*8 xg,yg,zg !Ecliptic rectangular coords
real*8 Ms !Mean anomaly of sun
real*8 ws !Argument of perihelion of sun
real*8 Ls !Mean longitude of sun (Ns=0)
real*8 Lm !Mean longitude of moon
real*8 DD !Mean elongation of moon
real*8 FF !Argument of latitude for moon
real*8 xe,ye,ze !Equatorial geocentric coords of moon
real*8 mpar !Parallax of moon (r_E / d)
real*8 lat,lon !Station coordinates on earth
real*8 gclat !Geocentric latitude
real*8 rho !Earth radius factor
real*8 GMST0,LST,HA
real*8 g
real*8 topRA,topDec !Topocentric coordinates of Moon
real*8 Az,El
real*8 dist
real*8 rad,twopi,pi,pio2
data rad/57.2957795131d0/,twopi/6.283185307d0/
d=367*y - 7*(y+(m+9)/12)/4 + 275*m/9 + Day - 730530 + UT/24.d0
ecl = 23.4393d0 - 3.563d-7 * d
! Orbital elements for Moon:
NN = 125.1228d0 - 0.0529538083d0 * d
i = 5.1454d0
w = mod(318.0634d0 + 0.1643573223d0 * d + 360000.d0,360.d0)
a = 60.2666d0
e = 0.054900d0
MM = mod(115.3654d0 + 13.0649929509d0 * d + 360000.d0,360.d0)
EE = MM + e*rad*sin(MM/rad) * (1.d0 + e*cos(MM/rad))
EE = EE - (EE - e*rad*sin(EE/rad)-MM) / (1.d0 - e*cos(EE/rad))
EE = EE - (EE - e*rad*sin(EE/rad)-MM) / (1.d0 - e*cos(EE/rad))
xv = a * (cos(EE/rad) - e)
yv = a * (sqrt(1.d0-e*e) * sin(EE/rad))
v = mod(rad*atan2(yv,xv)+720.d0,360.d0)
r = sqrt(xv*xv + yv*yv)
! Get geocentric position in ecliptic rectangular coordinates:
xg = r * (cos(NN/rad)*cos((v+w)/rad) - &
sin(NN/rad)*sin((v+w)/rad)*cos(i/rad))
yg = r * (sin(NN/rad)*cos((v+w)/rad) + &
cos(NN/rad)*sin((v+w)/rad)*cos(i/rad))
zg = r * (sin((v+w)/rad)*sin(i/rad))
! Ecliptic longitude and latitude of moon:
lonecl = mod(rad*atan2(yg/rad,xg/rad)+720.d0,360.d0)
latecl = rad*atan2(zg/rad,sqrt(xg*xg + yg*yg)/rad)
! Now include orbital perturbations:
Ms = mod(356.0470d0 + 0.9856002585d0 * d + 3600000.d0,360.d0)
ws = 282.9404d0 + 4.70935d-5*d
Ls = mod(Ms + ws + 720.d0,360.d0)
Lm = mod(MM + w + NN+720.d0,360.d0)
DD = mod(Lm - Ls + 360.d0,360.d0)
FF = mod(Lm - NN + 360.d0,360.d0)
lonecl = lonecl &
-1.274d0 * sin((MM-2.d0*DD)/rad) &
+0.658d0 * sin(2.d0*DD/rad) &
-0.186d0 * sin(Ms/rad) &
-0.059d0 * sin((2.d0*MM-2.d0*DD)/rad) &
-0.057d0 * sin((MM-2.d0*DD+Ms)/rad) &
+0.053d0 * sin((MM+2.d0*DD)/rad) &
+0.046d0 * sin((2.d0*DD-Ms)/rad) &
+0.041d0 * sin((MM-Ms)/rad) &
-0.035d0 * sin(DD/rad) &
-0.031d0 * sin((MM+Ms)/rad) &
-0.015d0 * sin((2.d0*FF-2.d0*DD)/rad) &
+0.011d0 * sin((MM-4.d0*DD)/rad)
latecl = latecl &
-0.173d0 * sin((FF-2.d0*DD)/rad) &
-0.055d0 * sin((MM-FF-2.d0*DD)/rad) &
-0.046d0 * sin((MM+FF-2.d0*DD)/rad) &
+0.033d0 * sin((FF+2.d0*DD)/rad) &
+0.017d0 * sin((2.d0*MM+FF)/rad)
r = 60.36298d0 &
- 3.27746d0*cos(MM/rad) &
- 0.57994d0*cos((MM-2.d0*DD)/rad) &
- 0.46357d0*cos(2.d0*DD/rad) &
- 0.08904d0*cos(2.d0*MM/rad) &
+ 0.03865d0*cos((2.d0*MM-2.d0*DD)/rad) &
- 0.03237d0*cos((2.d0*DD-Ms)/rad) &
- 0.02688d0*cos((MM+2.d0*DD)/rad) &
- 0.02358d0*cos((MM-2.d0*DD+Ms)/rad) &
- 0.02030d0*cos((MM-Ms)/rad) &
+ 0.01719d0*cos(DD/rad) &
+ 0.01671d0*cos((MM+Ms)/rad)
dist=r*6378.140d0
! Geocentric coordinates:
! Rectangular ecliptic coordinates of the moon:
xg = r * cos(lonecl/rad)*cos(latecl/rad)
yg = r * sin(lonecl/rad)*cos(latecl/rad)
zg = r * sin(latecl/rad)
! Rectangular equatorial coordinates of the moon:
xe = xg
ye = yg*cos(ecl/rad) - zg*sin(ecl/rad)
ze = yg*sin(ecl/rad) + zg*cos(ecl/rad)
! Right Ascension, Declination:
RA = mod(rad*atan2(ye,xe)+360.d0,360.d0)
Dec = rad*atan2(ze,sqrt(xe*xe + ye*ye))
! Now convert to topocentric system:
mpar=rad*asin(1.d0/r)
! alt_topoc = alt_geoc - mpar*cos(alt_geoc)
gclat = lat - 0.1924d0*sin(2.d0*lat/rad)
rho = 0.99883d0 + 0.00167d0*cos(2.d0*lat/rad)
GMST0 = (Ls + 180.d0)/15.d0
LST = mod(GMST0+UT+lon/15.d0+48.d0,24.d0) !LST in hours
HA = 15.d0*LST - RA !HA in degrees
g = rad*atan(tan(gclat/rad)/cos(HA/rad))
topRA = RA - mpar*rho*cos(gclat/rad)*sin(HA/rad)/cos(Dec/rad)
topDec = Dec - mpar*rho*sin(gclat/rad)*sin((g-Dec)/rad)/sin(g/rad)
HA = 15.d0*LST - topRA !HA in degrees
if(HA.gt.180.d0) HA=HA-360.d0
if(HA.lt.-180.d0) HA=HA+360.d0
pi=0.5d0*twopi
pio2=0.5d0*pi
call dcoord(pi,pio2-lat/rad,0.d0,lat/rad,ha*twopi/360,topDec/rad,az,el)
Az=az*rad
El=El*rad
return
end subroutine moon2

73
lib/moondop.f90 Normal file
View File

@ -0,0 +1,73 @@
subroutine MoonDop(nyear,month,nday,uth4,lon4,lat4,RAMoon4, &
DecMoon4,LST4,HA4,AzMoon4,ElMoon4,vr4,dist4)
implicit real*8 (a-h,o-z)
real*4 uth4 !UT in hours
real*4 lon4 !West longitude, degrees
real*4 lat4 !Latitude, degrees
real*4 RAMoon4 !Topocentric RA of moon, hours
real*4 DecMoon4 !Topocentric Dec of Moon, degrees
real*4 LST4 !Locat sidereal time, hours
real*4 HA4 !Local Hour angle, degrees
real*4 AzMoon4 !Topocentric Azimuth of moon, degrees
real*4 ElMoon4 !Topocentric Elevation of moon, degrees
real*4 vr4 !Radial velocity of moon wrt obs, km/s
real*4 dist4 !Echo time, seconds
real*8 LST
real*8 RME(6) !Vector from Earth center to Moon
real*8 RAE(6) !Vector from Earth center to Obs
real*8 RMA(6) !Vector from Obs to Moon
real*8 pvsun(6)
real*8 rme0(6)
logical km,bary
data rad/57.2957795130823d0/,twopi/6.28310530717959d0/
km=.true.
dlat=lat4/rad
dlong1=lon4/rad
elev1=200.d0
call geocentric(dlat,elev1,dlat1,erad1)
dt=100.d0 !For numerical derivative, in seconds
UT=uth4
! NB: geodetic latitude used here, but geocentric latitude used when
! determining Earth-rotation contribution to Doppler.
call moon2(nyear,month,nDay,UT-dt/3600.d0,dlong1*rad,dlat*rad, &
RA,Dec,topRA,topDec,LST,HA,Az0,El0,dist)
call toxyz(RA/rad,Dec/rad,dist,rme0) !Convert to rectangular coords
call moon2(nyear,month,nDay,UT,dlong1*rad,dlat*rad, &
RA,Dec,topRA,topDec,LST,HA,Az,El,dist)
call toxyz(RA/rad,Dec/rad,dist,rme) !Convert to rectangular coords
phi=LST*twopi/24.d0
call toxyz(phi,dlat1,erad1,rae) !Gencentric numbers used here!
radps=twopi/(86400.d0/1.002737909d0)
rae(4)=-rae(2)*radps !Vel of Obs wrt Earth center
rae(5)=rae(1)*radps
rae(6)=0.d0
do i=1,3
rme(i+3)=(rme(i)-rme0(i))/dt
rma(i)=rme(i)-rae(i)
rma(i+3)=rme(i+3)-rae(i+3)
enddo
call fromxyz(rma,alpha1,delta1,dtopo0) !Get topocentric coords
vr=dot(rma(4),rma)/dtopo0
RAMoon4=topRA
DecMoon4=topDec
LST4=LST
HA4=HA
AzMoon4=Az
ElMoon4=El
vr4=vr
dist4=dist
return
end subroutine MoonDop

88
lib/sun.f90 Normal file
View File

@ -0,0 +1,88 @@
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

14
lib/tm2.f90 Normal file
View File

@ -0,0 +1,14 @@
subroutine tm2(day,xlat4,xlon4,xl4,b4)
implicit real*8 (a-h,o-z)
parameter (RADS=0.0174532925199433d0)
real*4 day4,xlat4,xlon4,xl4,b4
glat=xlat4*RADS
glong=xlon4*RADS
call tmoonsub(day,glat,glong,el,rv,xl,b,pax)
xl4=xl
b4=b
end subroutine tm2

25
lib/toxyz.f90 Normal file
View File

@ -0,0 +1,25 @@
subroutine toxyz(alpha,delta,r,vec)
implicit real*8 (a-h,o-z)
real*8 vec(3)
vec(1)=r*cos(delta)*cos(alpha)
vec(2)=r*cos(delta)*sin(alpha)
vec(3)=r*sin(delta)
return
end subroutine toxyz
subroutine fromxyz(vec,alpha,delta,r)
implicit real*8 (a-h,o-z)
real*8 vec(3)
data twopi/6.283185307d0/
r=sqrt(vec(1)**2 + vec(2)**2 + vec(3)**2)
alpha=atan2(vec(2),vec(1))
if(alpha.lt.0.d0) alpha=alpha+twopi
delta=asin(vec(3)/r)
return
end subroutine fromxyz