1subroutine MoonDop(nyear,month,nday,uth4,lon4,lat4,RAMoon4,DecMoon4,   &
2     LST4,HA4,AzMoon4,ElMoon4,vr4,dist4)
3
4  implicit real*8 (a-h,o-z)
5  real*4 uth4                    !UT in hours
6  real*4 lon4                    !West longitude, degrees
7  real*4 lat4                    !Latitude, degrees
8  real*4 RAMoon4                 !Topocentric RA of moon, hours
9  real*4 DecMoon4                !Topocentric Dec of Moon, degrees
10  real*4 LST4                    !Locat sidereal time, hours
11  real*4 HA4                     !Local Hour angle, degrees
12  real*4 AzMoon4                 !Topocentric Azimuth of moon, degrees
13  real*4 ElMoon4                 !Topocentric Elevation of moon, degrees
14  real*4 vr4                     !Radial velocity of moon wrt obs, km/s
15  real*4 dist4                   !Echo time, seconds
16
17  real*8 LST
18  real*8 RME(6)                  !Vector from Earth center to Moon
19  real*8 RAE(6)                  !Vector from Earth center to Obs
20  real*8 RMA(6)                  !Vector from Obs to Moon
21  real*8 rme0(6)
22  logical km
23
24  data rad/57.2957795130823d0/,twopi/6.28310530717959d0/
25
26  km=.true.
27  dlat=lat4/rad
28  dlong1=lon4/rad
29  elev1=200.d0
30  call geocentric(dlat,elev1,dlat1,erad1)
31
32  dt=100.d0                       !For numerical derivative, in seconds
33  UT=uth4
34
35! NB: geodetic latitude used here, but geocentric latitude used when
36! determining Earth-rotation contribution to Doppler.
37
38  call moon2(nyear,month,nDay,UT-dt/3600.d0,dlong1*rad,dlat*rad,     &
39       RA,Dec,topRA,topDec,LST,HA,Az0,El0,dist)
40  call toxyz(RA/rad,Dec/rad,dist,rme0)      !Convert to rectangular coords
41
42  call moon2(nyear,month,nDay,UT,dlong1*rad,dlat*rad,                &
43       RA,Dec,topRA,topDec,LST,HA,Az,El,dist)
44  call toxyz(RA/rad,Dec/rad,dist,rme)       !Convert to rectangular coords
45
46  phi=LST*twopi/24.d0
47  call toxyz(phi,dlat1,erad1,rae)           !Gencentric numbers used here!
48  radps=twopi/(86400.d0/1.002737909d0)
49  rae(4)=-rae(2)*radps                      !Vel of Obs wrt Earth center
50  rae(5)=rae(1)*radps
51  rae(6)=0.d0
52
53  do i=1,3
54     rme(i+3)=(rme(i)-rme0(i))/dt
55     rma(i)=rme(i)-rae(i)
56     rma(i+3)=rme(i+3)-rae(i+3)
57  enddo
58
59  call fromxyz(rma,alpha1,delta1,dtopo0)     !Get topocentric coords
60  vr=dot(rma(4),rma)/dtopo0
61
62  RAMoon4=topRA
63  DecMoon4=topDec
64  LST4=LST
65  HA4=HA
66  AzMoon4=Az
67  ElMoon4=El
68  vr4=vr
69  dist4=dist
70
71  return
72end subroutine MoonDop
73