diff --git a/libm65/CMakeLists.txt b/libm65/CMakeLists.txt index 12b33e697..c424af52e 100644 --- a/libm65/CMakeLists.txt +++ b/libm65/CMakeLists.txt @@ -88,11 +88,19 @@ set (FSRCS getdphi.f90 getpfx1.f90 getpfx2.f90 + graycode.f90 graycode65.f90 + grid2deg.f90 + grid2k.f90 + indexx.f90 + interleave63.f90 iqcal.f90 iqfix.f90 jt65code.f90 + k2grid.f90 map65a.f90 + moon2.f90 + moondop.f90 noisegen.f90 pfxdump.f90 recvpkt.f90 @@ -107,14 +115,6 @@ set (FSRCS zplot.f90 f77_wisdom.f - graycode.f - grid2deg.f - grid2k.f - indexx.f - interleave63.f - k2grid.f - moon2.f - moondop.f nchar.f packcall.f packdxcc.f diff --git a/libm65/graycode.f b/libm65/graycode.f deleted file mode 100644 index 24c03e943..000000000 --- a/libm65/graycode.f +++ /dev/null @@ -1,10 +0,0 @@ - subroutine graycode(dat,n,idir) - - integer dat(n) - do i=1,n - dat(i)=igray(dat(i),idir) - enddo - - return - end - diff --git a/libm65/graycode.f90 b/libm65/graycode.f90 new file mode 100644 index 000000000..2074241a6 --- /dev/null +++ b/libm65/graycode.f90 @@ -0,0 +1,10 @@ +subroutine graycode(dat,n,idir) + + integer dat(n) + do i=1,n + dat(i)=igray(dat(i),idir) + enddo + + return +end subroutine graycode + diff --git a/libm65/grid2deg.f b/libm65/grid2deg.f deleted file mode 100644 index 7947073d8..000000000 --- a/libm65/grid2deg.f +++ /dev/null @@ -1,38 +0,0 @@ - subroutine grid2deg(grid0,dlong,dlat) - -C Converts Maidenhead grid locator to degrees of West longitude -C and North latitude. - - character*6 grid0,grid - character*1 g1,g2,g3,g4,g5,g6 - - grid=grid0 - i=ichar(grid(5:5)) - if(grid(5:5).eq.' ' .or. i.le.64 .or. i.ge.128) grid(5:6)='mm' - - if(grid(1:1).ge.'a' .and. grid(1:1).le.'z') grid(1:1)= - + char(ichar(grid(1:1))+ichar('A')-ichar('a')) - if(grid(2:2).ge.'a' .and. grid(2:2).le.'z') grid(2:2)= - + char(ichar(grid(2:2))+ichar('A')-ichar('a')) - if(grid(5:5).ge.'A' .and. grid(5:5).le.'Z') grid(5:5)= - + char(ichar(grid(5:5))-ichar('A')+ichar('a')) - if(grid(6:6).ge.'A' .and. grid(6:6).le.'Z') grid(6:6)= - + char(ichar(grid(6:6))-ichar('A')+ichar('a')) - - g1=grid(1:1) - g2=grid(2:2) - g3=grid(3:3) - g4=grid(4:4) - g5=grid(5:5) - g6=grid(6:6) - - nlong = 180 - 20*(ichar(g1)-ichar('A')) - n20d = 2*(ichar(g3)-ichar('0')) - xminlong = 5*(ichar(g5)-ichar('a')+0.5) - dlong = nlong - n20d - xminlong/60.0 - nlat = -90+10*(ichar(g2)-ichar('A')) + ichar(g4)-ichar('0') - xminlat = 2.5*(ichar(g6)-ichar('a')+0.5) - dlat = nlat + xminlat/60.0 - - return - end diff --git a/libm65/grid2deg.f90 b/libm65/grid2deg.f90 new file mode 100644 index 000000000..344351dd7 --- /dev/null +++ b/libm65/grid2deg.f90 @@ -0,0 +1,38 @@ +subroutine grid2deg(grid0,dlong,dlat) + +! Converts Maidenhead grid locator to degrees of West longitude +! and North latitude. + + character*6 grid0,grid + character*1 g1,g2,g3,g4,g5,g6 + + grid=grid0 + i=ichar(grid(5:5)) + if(grid(5:5).eq.' ' .or. i.le.64 .or. i.ge.128) grid(5:6)='mm' + + if(grid(1:1).ge.'a' .and. grid(1:1).le.'z') grid(1:1)= & + char(ichar(grid(1:1))+ichar('A')-ichar('a')) + if(grid(2:2).ge.'a' .and. grid(2:2).le.'z') grid(2:2)= & + char(ichar(grid(2:2))+ichar('A')-ichar('a')) + if(grid(5:5).ge.'A' .and. grid(5:5).le.'Z') grid(5:5)= & + char(ichar(grid(5:5))-ichar('A')+ichar('a')) + if(grid(6:6).ge.'A' .and. grid(6:6).le.'Z') grid(6:6)= & + char(ichar(grid(6:6))-ichar('A')+ichar('a')) + + g1=grid(1:1) + g2=grid(2:2) + g3=grid(3:3) + g4=grid(4:4) + g5=grid(5:5) + g6=grid(6:6) + + nlong = 180 - 20*(ichar(g1)-ichar('A')) + n20d = 2*(ichar(g3)-ichar('0')) + xminlong = 5*(ichar(g5)-ichar('a')+0.5) + dlong = nlong - n20d - xminlong/60.0 + nlat = -90+10*(ichar(g2)-ichar('A')) + ichar(g4)-ichar('0') + xminlat = 2.5*(ichar(g6)-ichar('a')+0.5) + dlat = nlat + xminlat/60.0 + + return +end subroutine grid2deg diff --git a/libm65/grid2k.f b/libm65/grid2k.f deleted file mode 100644 index 1306a95a2..000000000 --- a/libm65/grid2k.f +++ /dev/null @@ -1,12 +0,0 @@ - subroutine grid2k(grid,k) - - character*6 grid - - call grid2deg(grid,xlong,xlat) - nlong=nint(xlong) - nlat=nint(xlat) - k=0 - if(nlat.ge.85) k=5*(nlong+179)/2 + nlat-84 - - return - end diff --git a/libm65/grid2k.f90 b/libm65/grid2k.f90 new file mode 100644 index 000000000..f68b1409e --- /dev/null +++ b/libm65/grid2k.f90 @@ -0,0 +1,12 @@ +subroutine grid2k(grid,k) + + character*6 grid + + call grid2deg(grid,xlong,xlat) + nlong=nint(xlong) + nlat=nint(xlat) + k=0 + if(nlat.ge.85) k=5*(nlong+179)/2 + nlat-84 + + return +end subroutine grid2k diff --git a/libm65/indexx.f b/libm65/indexx.f deleted file mode 100644 index df2a5330e..000000000 --- a/libm65/indexx.f +++ /dev/null @@ -1,19 +0,0 @@ - subroutine indexx(n,arr,indx) - - parameter (NMAX=3000) - integer indx(n) - real arr(n) - real brr(NMAX) - if(n.gt.NMAX) then - print*,'n=',n,' too big in indexx.' - stop - endif - do i=1,n - brr(i)=arr(i) - indx(i)=i - enddo - call ssort(brr,indx,n,2) - - return - end - diff --git a/libm65/indexx.f90 b/libm65/indexx.f90 new file mode 100644 index 000000000..57c1ec075 --- /dev/null +++ b/libm65/indexx.f90 @@ -0,0 +1,19 @@ +subroutine indexx(n,arr,indx) + + parameter (NMAX=3000) + integer indx(n) + real arr(n) + real brr(NMAX) + if(n.gt.NMAX) then + print*,'n=',n,' too big in indexx.' + stop + endif + do i=1,n + brr(i)=arr(i) + indx(i)=i + enddo + call ssort(brr,indx,n,2) + + return +end subroutine indexx + diff --git a/libm65/interleave63.f b/libm65/interleave63.f deleted file mode 100644 index f6793023a..000000000 --- a/libm65/interleave63.f +++ /dev/null @@ -1,25 +0,0 @@ - subroutine interleave63(d1,idir) - -C Interleave (idir=1) or de-interleave (idir=-1) the array d1. - - integer d1(0:6,0:8) - integer d2(0:8,0:6) - - if(idir.ge.0) then - do i=0,6 - do j=0,8 - d2(j,i)=d1(i,j) - enddo - enddo - call move(d2,d1,63) - else - call move(d1,d2,63) - do i=0,6 - do j=0,8 - d1(i,j)=d2(j,i) - enddo - enddo - endif - - return - end diff --git a/libm65/interleave63.f90 b/libm65/interleave63.f90 new file mode 100644 index 000000000..a32ef34cd --- /dev/null +++ b/libm65/interleave63.f90 @@ -0,0 +1,25 @@ +subroutine interleave63(d1,idir) + +! Interleave (idir=1) or de-interleave (idir=-1) the array d1. + + integer d1(0:6,0:8) + integer d2(0:8,0:6) + + if(idir.ge.0) then + do i=0,6 + do j=0,8 + d2(j,i)=d1(i,j) + enddo + enddo + call move(d2,d1,63) + else + call move(d1,d2,63) + do i=0,6 + do j=0,8 + d1(i,j)=d2(j,i) + enddo + enddo + endif + + return +end subroutine interleave63 diff --git a/libm65/k2grid.f b/libm65/k2grid.f deleted file mode 100644 index 6fcd7f3e4..000000000 --- a/libm65/k2grid.f +++ /dev/null @@ -1,12 +0,0 @@ - subroutine k2grid(k,grid) - character grid*6 - - nlong=2*mod((k-1)/5,90)-179 - if(k.gt.450) nlong=nlong+180 - nlat=mod(k-1,5)+ 85 - dlat=nlat - dlong=nlong - call deg2grid(dlong,dlat,grid) - - return - end diff --git a/libm65/k2grid.f90 b/libm65/k2grid.f90 new file mode 100644 index 000000000..aa7631579 --- /dev/null +++ b/libm65/k2grid.f90 @@ -0,0 +1,12 @@ +subroutine k2grid(k,grid) + character grid*6 + + nlong=2*mod((k-1)/5,90)-179 + if(k.gt.450) nlong=nlong+180 + nlat=mod(k-1,5)+ 85 + dlat=nlat + dlong=nlong + call deg2grid(dlong,dlat,grid) + + return +end subroutine k2grid diff --git a/libm65/moon2.f b/libm65/moon2.f deleted file mode 100644 index 8144b675f..000000000 --- a/libm65/moon2.f +++ /dev/null @@ -1,167 +0,0 @@ - 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 - -C 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 - -C 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) - -C 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)) - -C 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) - -C 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 - -C Geocentric coordinates: -C 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) - -C 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) - -C Right Ascension, Declination: - RA = mod(rad*atan2(ye,xe)+360.d0,360.d0) - Dec = rad*atan2(ze,sqrt(xe*xe + ye*ye)) - -C Now convert to topocentric system: - mpar=rad*asin(1.d0/r) -C 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 diff --git a/libm65/moon2.f90 b/libm65/moon2.f90 new file mode 100644 index 000000000..9fa72f1df --- /dev/null +++ b/libm65/moon2.f90 @@ -0,0 +1,163 @@ +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/libm65/moondop.f b/libm65/moondop.f deleted file mode 100644 index ebac631c4..000000000 --- a/libm65/moondop.f +++ /dev/null @@ -1,73 +0,0 @@ - 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 - -C NB: geodetic latitude used here, but geocentric latitude used when -C 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 diff --git a/libm65/moondop.f90 b/libm65/moondop.f90 new file mode 100644 index 000000000..4808949f4 --- /dev/null +++ b/libm65/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