diff --git a/lib/astro.f90 b/lib/astro.f90 new file mode 100644 index 000000000..56711f358 --- /dev/null +++ b/lib/astro.f90 @@ -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 diff --git a/lib/astro0.f90 b/lib/astro0.f90 new file mode 100644 index 000000000..fe047d460 --- /dev/null +++ b/lib/astro0.f90 @@ -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 diff --git a/lib/coord.f90 b/lib/coord.f90 new file mode 100644 index 000000000..b60f71e60 --- /dev/null +++ b/lib/coord.f90 @@ -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 diff --git a/lib/geocentric.f90 b/lib/geocentric.f90 new file mode 100644 index 000000000..11ea7cb60 --- /dev/null +++ b/lib/geocentric.f90 @@ -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 + diff --git a/lib/moon2.f90 b/lib/moon2.f90 new file mode 100644 index 000000000..8498267e5 --- /dev/null +++ b/lib/moon2.f90 @@ -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 diff --git a/lib/moondop.f90 b/lib/moondop.f90 new file mode 100644 index 000000000..e2987d7f5 --- /dev/null +++ b/lib/moondop.f90 @@ -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 diff --git a/lib/sun.f90 b/lib/sun.f90 new file mode 100644 index 000000000..db258bb26 --- /dev/null +++ b/lib/sun.f90 @@ -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 diff --git a/lib/tm2.f90 b/lib/tm2.f90 new file mode 100644 index 000000000..9d7777750 --- /dev/null +++ b/lib/tm2.f90 @@ -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 diff --git a/lib/toxyz.f90 b/lib/toxyz.f90 new file mode 100644 index 000000000..646dd7a2a --- /dev/null +++ b/lib/toxyz.f90 @@ -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