86 lines
3.0 KiB
Fortran
86 lines
3.0 KiB
Fortran
subroutine ephem(mjd0,dut,east_long,geodetic_lat,height,nspecial, &
|
|
RA,Dec,Az,El,techo,dop,fspread_1GHz,vr)
|
|
|
|
implicit real*8 (a-h,o-z)
|
|
real*8 jd !Time of observationa as a Julian Date
|
|
real*8 mjd,mjd0 !Modified Julian Date
|
|
real*8 prec(3,3) !Precession matrix, J2000 to Date
|
|
real*8 rmatn(3,3) !Nutation matrix
|
|
real*8 rme2000(6) !Vector from Earth center to Moon, JD2000
|
|
real*8 rmeDate(6) !Vector from Earth center to Moon at Date
|
|
real*8 rmeTrue(6) !Include nutation
|
|
real*8 raeTrue(6) !Vector from Earth center to Obs at Date
|
|
real*8 rmaTrue(6) !Vector from Obs to Moon at Date
|
|
logical km,bary !Set km=.true. to get km, km/s from ephemeris
|
|
common/stcomx/km,bary,pvsun(6) !Common used in JPL subroutines
|
|
common/librcom/xl(2),b(2)
|
|
|
|
twopi=8.d0*atan(1.d0) !Define some constants
|
|
rad=360.d0/twopi
|
|
clight=2.99792458d5
|
|
au2km=0.1495978706910000d9
|
|
pi=0.5d0*twopi
|
|
pio2=0.5d0*pi
|
|
km=.true.
|
|
freq=1000.0d6
|
|
|
|
do jj=1,2
|
|
mjd=mjd0
|
|
if(jj.eq.1) mjd=mjd - 1.d0/1440.d0
|
|
djutc=mjd
|
|
jd=2400000.5d0 + mjd
|
|
djtt=mjd + sla_DTT(jd)/86400.d0
|
|
ttjd=jd + sla_DTT(jd)/86400.d0
|
|
|
|
if(nspecial.ne.8) then
|
|
call pleph(ttjd,10,3,rme2000) !RME (J2000) from JPL ephemeris
|
|
|
|
year=2000.d0 + (jd-2451545.d0)/365.25d0
|
|
call sla_PREC (2000.0d0, year, prec) !Get precession matrix
|
|
rmeDate(1:3)=matmul(prec,rme2000(1:3)) !Moon geocentric xyz at Date
|
|
rmeDate(4:6)=matmul(prec,rme2000(4:6)) !Moon geocentric vel at Date
|
|
else
|
|
call sla_DMOON(djtt,rmeDate) !No JPL ephemeris, use DMOON
|
|
rmeDate=rmeDate*au2km
|
|
endif
|
|
|
|
if(nspecial.eq.7) then
|
|
rmeTrue=rmeDate
|
|
else
|
|
!Nutation to true equinox of Date
|
|
call sla_NUT(djtt,rmatn)
|
|
call sla_DMXV(rmatn,rmeDate,rmeTrue)
|
|
call sla_DMXV(rmatn,rmeDate(4),rmeTrue(4))
|
|
endif
|
|
|
|
! Local Apparent Sidereal Time:
|
|
djut1=djutc + dut/86400.d0
|
|
if(nspecial.eq.6) djut1=djutc
|
|
xlast=sla_DRANRM(sla_GMST(djut1) + sla_EQEQX(djtt) + east_long)
|
|
call sla_PVOBS(geodetic_lat,height,xlast,raeTrue)
|
|
rmaTrue=rmeTrue - raeTrue*au2km
|
|
|
|
if(nspecial.ne.2) then
|
|
! Allow for planetary aberration
|
|
tl=499.004782D0*SQRT(rmaTrue(1)**2 + rmaTrue(2)**2 + rmaTrue(3)**2)
|
|
rmaTrue(1:3)=rmaTrue(1:3)-tl*rmaTrue(4:6)/au2km
|
|
endif
|
|
|
|
!Topocentric RA, Dec, dist, velocity
|
|
call sla_DC62S(rmaTrue,RA,Dec,dist,RAdot,DECdot,vr)
|
|
dop=-2.d0 * freq * vr/clight !EME doppler shift
|
|
techo=2.d0*dist/clight !Echo delay time (s)
|
|
call libration(jd,RA,Dec,xl(jj),b(jj))
|
|
enddo
|
|
|
|
fspread_1GHz=0.0d0
|
|
dldt=57.2957795131*(xl(2)-xl(1))
|
|
dbdt=57.2957795131*(b(2)-b(1))
|
|
rate=sqrt((2*dldt)**2 + (2*dbdt)**2)
|
|
fspread_1GHz=0.5*6741*rate
|
|
|
|
call sla_DE2H(xlast-RA,Dec,geodetic_lat,Az,El)
|
|
|
|
return
|
|
end subroutine ephem
|