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

EME Doppler calculations now use JPL Planetary Ephemeris.

Browse Source 
Also new, simplified routines for Doppler spread.
Beware! Not yet fully tested ...



git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@5496 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
This commit is contained in:
Joe Taylor 2015-06-02 16:43:40 +00:00
parent adebd2420a
commit e837b9209e
17 changed files with 4495 additions and 378 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
CMakeLists.txt
View File

@ -268,11 +268,8 @@ set (wsjt_FSRCS
lib/chkss2.f90
lib/coord.f90
lib/db.f90
lib/dcoord.f90
lib/decode65a.f90
lib/decode65b.f90
lib/fftw3mod.f90
lib/jt9fano.f90
lib/decode4.f90
lib/decoder.f90
lib/decjt9.f90
@ -280,18 +277,17 @@ set (wsjt_FSRCS
lib/deg2grid.f90
lib/demod64a.f90
lib/determ.f90
lib/dot.f90
lib/downsam9.f90
lib/encode232.f90
lib/encode4.f90
lib/entail.f90
lib/ephem.f90
lib/extract.F90
lib/extract4.f90
lib/geocentric.f90
lib/getmet4.f90
lib/fano232.f90
lib/fchisq.f90
lib/fchisq65.f90
lib/fftw3mod.f90
lib/fil3.f90
lib/fil3c.f90
lib/fil4.f90
@ -311,6 +307,7 @@ set (wsjt_FSRCS
lib/genwspr.f90
lib/geodist.f90
lib/getlags.f90
lib/getmet4.f90
lib/graycode.f90
lib/graycode65.f90
lib/grayline.f90
@ -323,13 +320,15 @@ set (wsjt_FSRCS
lib/interleave63.f90
lib/interleave9.f90
lib/inter_wspr.f90
lib/jplsubs.f
lib/jt4.f90
lib/jt4a.f90
lib/jt65a.f90
lib/jt9fano.f90
lib/libration.f90
lib/lpf1.f90
lib/mixlpf.f90
lib/moon2.f90
lib/moondop.f90
lib/moondopjpl.f90
lib/morse.f90
lib/move.f90
lib/options.f90
@ -343,6 +342,7 @@ set (wsjt_FSRCS
lib/savec2.f90
lib/sec_midn.f90
lib/setup65.f90
lib/slasubs.f
lib/sleep_msec.f90
lib/slope.f90
lib/smo.f90
@ -358,8 +358,6 @@ set (wsjt_FSRCS
lib/sync9.f90
lib/timer.f90
lib/timf2.f90
lib/tm2.f90
lib/toxyz.f90
lib/twkfreq.f90
lib/twkfreq65.f90
lib/wav11.f90
@ -1002,6 +1000,12 @@ install (FILES
#COMPONENT runtime
)
install (FILES
contrib/Ephemeris/JPLEPH
DESTINATION ${WSJT_SHARE_DESTINATION}/${WSJT_DATA_DESTINATION}
#COMPONENT runtime
)
#
# Mac installer files
#

BIN
contrib/Ephemeris/JPLEPH Normal file
View File

Binary file not shown.

6
lib/astro.f90
View File

@ -42,11 +42,7 @@ subroutine astro(nyear,month,nday,uth,freq8,Mygrid, &
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, &
call MoonDopJPL(nyear,month,nday,uth,lon,lat,RAMoon,DecMoon, &
LST,HA,AzMoon,ElMoon,vr,dist)
! Compute spatial polarization offset

29
lib/astro0.f90
View File

@ -1,14 +1,16 @@
subroutine astro0(nyear,month,nday,uth8,freq8,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)
width1,width2,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 sd8,poloffset8,day8,width1,width2,xlst8
real*8 uth8,techo8,freq8
real*8 xl,b
common/librcom/xl(2),b(2)
data uth8z/0.d0/
save
@ -19,18 +21,19 @@ subroutine astro0(nyear,month,nday,uth8,freq8,mygrid,hisgrid, &
day,xlon2,xlat2,xlst,techo)
AzMoonB8=AzMoon
ElMoonB8=ElMoon
xl2=xl(1)
xl2a=xl(2)
b2=b(1)
b2a=b(2)
call astro(nyear,month,nday,uth,freq8,mygrid,1,1, &
AzSun,ElSun,AzMoon,ElMoon,ntsky,doppler00,doppler, &
dbMoon,RAMoon,DecMoon,HA,Dgrd,sd,poloffset,xnr, &
day,xlon1,xlat1,xlst,techo)
xl1=xl(1)
xl1a=xl(2)
b1=b(1)
b1a=b(2)
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=1.d-9*freq8
dldt1=DEGS*(xl1a-xl1)
dbdt1=DEGS*(b1a-b1)
@ -41,14 +44,6 @@ subroutine astro0(nyear,month,nday,uth8,freq8,mygrid,hisgrid, &
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

6
lib/astrosub.f90
View File

@ -6,12 +6,16 @@ subroutine astrosub(nyear,month,nday,uth8,freq8,mygrid,hisgrid, &
implicit real*8 (a-h,o-z)
character*6 mygrid,hisgrid,c1*1
character*6 AzElFileName*(*),jpleph*(*)
character*80 jpleph_file_name
logical*1 bTx
common/jplcom/jpleph_file_name
jpleph_file_name=jpleph
call astro0(nyear,month,nday,uth8,freq8,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)
width1,width2,xlst8,techo8)
imin=60*uth8
isec=3600*uth8

40
lib/dcoord.f90
View File

@ -1,40 +0,0 @@
SUBROUTINE DCOORD(A0,B0,AP,BP,A1,B1,A2,B2)
implicit real*8 (a-h,o-z)
! 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 ecliptic latitude (eb) and longitude (el) to ra, dec:
! call coord(0.d0,0.d0,-pio2,pio2-23.443*pi/180,ra,dec,el,eb)
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.D0-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.D0) TA2O2=(1.D0-CA2)/SA2
IF(CA2.GT.0.D0) TA2O2=SA2/(1.D0+CA2)
A2=2.D0*atan(TA2O2)
IF(A2.LT.0.D0) A2=A2+6.2831853071795864D0
RETURN
END SUBROUTINE DCOORD

11
lib/dot.f90
View File

@ -1,11 +0,0 @@
real*8 function dot(x,y)
real*8 x(3),y(3)
dot=0.d0
do i=1,3
dot=dot+x(i)*y(i)
enddo
return
end function dot

87
lib/ephem.f90 Normal file
View File

@ -0,0 +1,87 @@
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,jplok !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
! inquire(file='JPLEPH',exist=jplok)
! if(jplok) then
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

17
lib/geocentric.f90
View File

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

899
lib/jplsubs.f Normal file
View File

@ -0,0 +1,899 @@
C++++++++++++++++++++++++
C
SUBROUTINE FSIZER1(NRECL,KSIZE,NRFILE,NAMFIL)
C
C++++++++++++++++++++++++
C
C Version 1.0 uses the INQUIRE statement to find out the the record length
C of the direct access file before opening it. This procedure is non-standard,
C but seems to work for VAX machines.
C
C THE SUBROUTINE ALSO SETS THE VALUES OF NRECL, NRFILE, AND NAMFIL.
C *****************************************************************
C *****************************************************************
C
C THE PARAMETERS NAMFIL, NRECL, AND NRFILE ARE TO BE SET BY THE USER
C
C *****************************************************************
C NAMFIL IS THE EXTERNAL NAME OF THE BINARY EPHEMERIS FILE
CHARACTER*80 NAMFIL
c NAMFIL='JPLEPH'
C *****************************************************************
C NRECL=1 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN S.P. WORDS
C NRECL=4 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN BYTES
C (for a VAX, it is probably 1)
C
NRECL=4
C *****************************************************************
C NRFILE IS THE INTERNAL UNIT NUMBER USED FOR THE EPHEMERIS FILE
c NRFILE=12
C *****************************************************************
C FIND THE RECORD SIZE USING THE INQUIRE STATEMENT
c IRECSZ=0
INQUIRE(FILE=NAMFIL,RECL=IRECSZ)
C IF 'INQUIRE' DOES NOT WORK, USUALLY IRECSZ WILL BE LEFT AT 0
IF(IRECSZ .LE. 0) write(*,*)
. ' INQUIRE STATEMENT PROBABLY DID NOT WORK'
KSIZE=IRECSZ/NRECL
RETURN
END
C++++++++++++++++++++++++
C
SUBROUTINE FSIZER2(NRECL,KSIZE,NRFILE,NAMFIL)
C
C++++++++++++++++++++++++
C THIS SUBROUTINE OPENS THE FILE, 'NAMFIL', WITH A PHONY RECORD LENGTH, READS
C THE FIRST RECORD, AND USES THE INFO TO COMPUTE KSIZE, THE NUMBER OF SINGLE
C PRECISION WORDS IN A RECORD.
C
C THE SUBROUTINE ALSO SETS THE VALUES OF NRECL, NRFILE, AND NAMFIL.
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
SAVE
INTEGER OLDMAX
PARAMETER (OLDMAX = 400)
INTEGER NMAX
PARAMETER (NMAX = 1000)
CHARACTER*6 TTL(14,3),CNAM(NMAX)
CHARACTER*80 NAMFIL,jpleph_file_name
DIMENSION SS(3)
INTEGER IPT(3,13)
common/jplcom/jpleph_file_name
C *****************************************************************
C *****************************************************************
C
C THE PARAMETERS NRECL, NRFILE, AND NAMFIL ARE TO BE SET BY THE USER
C
C *****************************************************************
C NRECL=1 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN S.P. WORDS
C NRECL=4 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN BYTES
C (for UNIX, it is probably 4)
C
NRECL=4
C NRFILE IS THE INTERNAL UNIT NUMBER USED FOR THE EPHEMERIS FILE
NRFILE=12
C NAMFIL IS THE EXTERNAL NAME OF THE BINARY EPHEMERIS FILE
! NAMFIL='JPLEPH'
NAMFIL=jpleph_file_name
C *****************************************************************
C *****************************************************************
C ** OPEN THE DIRECT-ACCESS FILE AND GET THE POINTERS IN ORDER TO
C ** DETERMINE THE SIZE OF THE EPHEMERIS RECORD
MRECL=NRECL*1000
OPEN(NRFILE,
* FILE=NAMFIL,
* ACCESS='DIRECT',
* FORM='UNFORMATTED',
* RECL=MRECL,
* STATUS='OLD')
READ(NRFILE,REC=1)TTL,(CNAM(K),K=1,OLDMAX),SS,NCON,AU,EMRAT,
& ((IPT(I,J),I=1,3),J=1,12),NUMDE,(IPT(I,13),I=1,3)
CLOSE(NRFILE)
C FIND THE NUMBER OF EPHEMERIS COEFFICIENTS FROM THE POINTERS
KMX = 0
KHI = 0
DO I = 1,13
IF (IPT(1,I) .GT. KMX) THEN
KMX = IPT(1,I)
KHI = I
ENDIF
ENDDO
ND = 3
IF (KHI .EQ. 12) ND=2
KSIZE = 2*(IPT(1,KHI)+ND*IPT(2,KHI)*IPT(3,KHI)-1)
RETURN
END
C++++++++++++++++++++++++
C
SUBROUTINE FSIZER3(NRECL,KSIZE,NRFILE,NAMFIL)
C
C++++++++++++++++++++++++
C
C THE SUBROUTINE SETS THE VALUES OF NRECL, KSIZE, NRFILE, AND NAMFIL.
SAVE
CHARACTER*80 NAMFIL,jpleph_file_name
common/jplcom/jpleph_file_name
C *****************************************************************
C *****************************************************************
C
C THE PARAMETERS NRECL, NRFILE, AND NAMFIL ARE TO BE SET BY THE USER
C *****************************************************************
C NRECL=1 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN S.P. WORDS
C NRECL=4 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN BYTES
NRECL=4
C *****************************************************************
C NRFILE IS THE INTERNAL UNIT NUMBER USED FOR THE EPHEMERIS FILE (DEFAULT: 12)
NRFILE=12
C *****************************************************************
C NAMFIL IS THE EXTERNAL NAME OF THE BINARY EPHEMERIS FILE
! NAMFIL='JPLEPH'
NAMFIL=jpleph_file_name
C *****************************************************************
C KSIZE must be set by the user according to the ephemeris to be read
C For de200, set KSIZE to 1652
C For de405, set KSIZE to 2036
C For de406, set KSIZE to 1456
C For de414, set KSIZE to 2036
C For de418, set KSIZE to 2036
C For de421, set KSIZE to 2036
C For de422, set KSIZE to 2036
C For de423, set KSIZE to 2036
C For de424, set KSIZE to 2036
C For de430, set KSIZE to 2036
KSIZE = 2036
C *******************************************************************
RETURN
END
C++++++++++++++++++++++++++
C
SUBROUTINE PLEPH ( ET, NTARG, NCENT, RRD )
C
C++++++++++++++++++++++++++
C NOTE : Over the years, different versions of PLEPH have had a fifth argument:
C sometimes, an error return statement number; sometimes, a logical denoting
C whether or not the requested date is covered by the ephemeris. We apologize
C for this inconsistency; in this present version, we use only the four necessary
C arguments and do the testing outside of the subroutine.
C
C THIS SUBROUTINE READS THE JPL PLANETARY EPHEMERIS
C AND GIVES THE POSITION AND VELOCITY OF THE POINT 'NTARG'
C WITH RESPECT TO 'NCENT'.
C
C CALLING SEQUENCE PARAMETERS:
C
C ET = D.P. JULIAN EPHEMERIS DATE AT WHICH INTERPOLATION
C IS WANTED.
C
C ** NOTE THE ENTRY DPLEPH FOR A DOUBLY-DIMENSIONED TIME **
C THE REASON FOR THIS OPTION IS DISCUSSED IN THE
C SUBROUTINE STATE
C
C NTARG = INTEGER NUMBER OF 'TARGET' POINT.
C
C NCENT = INTEGER NUMBER OF CENTER POINT.
C
C THE NUMBERING CONVENTION FOR 'NTARG' AND 'NCENT' IS:
C
C 1 = MERCURY 8 = NEPTUNE
C 2 = VENUS 9 = PLUTO
C 3 = EARTH 10 = MOON
C 4 = MARS 11 = SUN
C 5 = JUPITER 12 = SOLAR-SYSTEM BARYCENTER
C 6 = SATURN 13 = EARTH-MOON BARYCENTER
C 7 = URANUS 14 = NUTATIONS (LONGITUDE AND OBLIQ)
C 15 = LIBRATIONS, IF ON EPH FILE
C
C (IF NUTATIONS ARE WANTED, SET NTARG = 14. FOR LIBRATIONS,
C SET NTARG = 15. SET NCENT=0.)
C
C RRD = OUTPUT 6-WORD D.P. ARRAY CONTAINING POSITION AND VELOCITY
C OF POINT 'NTARG' RELATIVE TO 'NCENT'. THE UNITS ARE AU AND
C AU/DAY. FOR LIBRATIONS THE UNITS ARE RADIANS AND RADIANS
C PER DAY. IN THE CASE OF NUTATIONS THE FIRST FOUR WORDS OF
C RRD WILL BE SET TO NUTATIONS AND RATES, HAVING UNITS OF
C RADIANS AND RADIANS/DAY.
C
C The option is available to have the units in km and km/sec.
C For this, set km=.true. in the STCOMX common block.
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
INTEGER NMAX
PARAMETER (NMAX = 1000)
DIMENSION RRD(6),ET2Z(2),ET2(2),PV(6,13)
DIMENSION PVST(6,11),PNUT(4)
DIMENSION SS(3),CVAL(NMAX),PVSUN(6),ZIPS(2)
DATA ZIPS/2*0.d0/
LOGICAL BSAVE,KM,BARY
LOGICAL FIRST
DATA FIRST/.TRUE./
INTEGER LIST(12),IPT(39),DENUM
COMMON/EPHHDR/CVAL,SS,AU,EMRAT,DENUM,NCON,IPT
COMMON/STCOMX/KM,BARY,PVSUN
common/jplcom/jpleph_file_name
C INITIALIZE ET2 FOR 'STATE' AND SET UP COMPONENT COUNT
C
ET2(1)=ET
ET2(2)=0.D0
GO TO 11
C ENTRY POINT 'DPLEPH' FOR DOUBLY-DIMENSIONED TIME ARGUMENT
C (SEE THE DISCUSSION IN THE SUBROUTINE STATE)
ENTRY DPLEPH(ET2Z,NTARG,NCENT,RRD)
ET2(1)=ET2Z(1)
ET2(2)=ET2Z(2)
11 IF(FIRST) CALL STATE(ZIPS,LIST,PVST,PNUT)
FIRST=.FALSE.
96 IF(NTARG .EQ. NCENT) RETURN
DO I=1,12
LIST(I)=0
ENDDO
C CHECK FOR NUTATION CALL
IF(NTARG.NE.14) GO TO 97
IF(IPT(35).GT.0) THEN
LIST(11)=2
CALL STATE(ET2,LIST,PVST,PNUT)
DO I=1,4
RRD(I)=PNUT(I)
ENDDO
RRD(5) = 0.d0
RRD(6) = 0.d0
RETURN
ELSE
DO I=1,4
RRD(I)=0.d0
ENDDO
WRITE(6,297)
297 FORMAT(' ***** NO NUTATIONS ON THE EPHEMERIS FILE *****')
STOP
ENDIF
C CHECK FOR LIBRATIONS
97 CONTINUE
DO I=1,6
RRD(I)=0.d0
ENDDO
IF(NTARG.NE.15) GO TO 98
IF(IPT(38).GT.0) THEN
LIST(12)=2
CALL STATE(ET2,LIST,PVST,PNUT)
DO I=1,6
RRD(I)=PVST(I,11)
ENDDO
RETURN
ELSE
WRITE(6,298)
298 FORMAT(' ***** NO LIBRATIONS ON THE EPHEMERIS FILE *****')
STOP
ENDIF
C FORCE BARYCENTRIC OUTPUT BY 'STATE'
98 BSAVE=BARY
BARY=.TRUE.
C SET UP PROPER ENTRIES IN 'LIST' ARRAY FOR STATE CALL
DO I=1,2
K=NTARG
IF(I .EQ. 2) K=NCENT
IF(K .LE. 10) LIST(K)=2
IF(K .EQ. 10) LIST(3)=2
IF(K .EQ. 3) LIST(10)=2
IF(K .EQ. 13) LIST(3)=2
ENDDO
C MAKE CALL TO STATE
CALL STATE(ET2,LIST,PVST,PNUT)
DO I=1,10
DO J = 1,6
PV(J,I) = PVST(J,I)
ENDDO
ENDDO
IF(NTARG .EQ. 11 .OR. NCENT .EQ. 11) THEN
DO I=1,6
PV(I,11)=PVSUN(I)
ENDDO
ENDIF
IF(NTARG .EQ. 12 .OR. NCENT .EQ. 12) THEN
DO I=1,6
PV(I,12)=0.D0
ENDDO
ENDIF
IF(NTARG .EQ. 13 .OR. NCENT .EQ. 13) THEN
DO I=1,6
PV(I,13) = PVST(I,3)
ENDDO
ENDIF
IF(NTARG*NCENT .EQ. 30 .AND. NTARG+NCENT .EQ. 13) THEN
DO I=1,6
PV(I,3)=0.D0
ENDDO
GO TO 99
ENDIF
IF(LIST(3) .EQ. 2) THEN
DO I=1,6
PV(I,3)=PVST(I,3)-PVST(I,10)/(1.D0+EMRAT)
ENDDO
ENDIF
IF(LIST(10) .EQ. 2) THEN
DO I=1,6
PV(I,10) = PV(I,3)+PVST(I,10)
ENDDO
ENDIF
99 DO I=1,6
RRD(I)=PV(I,NTARG)-PV(I,NCENT)
ENDDO
BARY=BSAVE
RETURN
END
C+++++++++++++++++++++++++++++++++
C
SUBROUTINE INTERP(BUF,T,NCF,NCM,NA,IFL,PV)
C
C+++++++++++++++++++++++++++++++++
C
C THIS SUBROUTINE DIFFERENTIATES AND INTERPOLATES A
C SET OF CHEBYSHEV COEFFICIENTS TO GIVE POSITION AND VELOCITY
C
C CALLING SEQUENCE PARAMETERS:
C
C INPUT:
C
C BUF 1ST LOCATION OF ARRAY OF D.P. CHEBYSHEV COEFFICIENTS OF POSITION
C
C T T(1) IS DP FRACTIONAL TIME IN INTERVAL COVERED BY
C COEFFICIENTS AT WHICH INTERPOLATION IS WANTED
C (0 .LE. T(1) .LE. 1). T(2) IS DP LENGTH OF WHOLE
C INTERVAL IN INPUT TIME UNITS.
C
C NCF # OF COEFFICIENTS PER COMPONENT
C
C NCM # OF COMPONENTS PER SET OF COEFFICIENTS
C
C NA # OF SETS OF COEFFICIENTS IN FULL ARRAY
C (I.E., # OF SUB-INTERVALS IN FULL INTERVAL)
C
C IFL INTEGER FLAG: =1 FOR POSITIONS ONLY
C =2 FOR POS AND VEL
C
C
C OUTPUT:
C
C PV INTERPOLATED QUANTITIES REQUESTED. DIMENSION
C EXPECTED IS PV(NCM,IFL), DP.
C
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
SAVE
C
DOUBLE PRECISION BUF(NCF,NCM,*),T(2),PV(NCM,*),PC(18),VC(18)
C
DATA NP/2/
DATA NV/3/
DATA TWOT/0.D0/
DATA PC(1),PC(2)/1.D0,0.D0/
DATA VC(2)/1.D0/
C
C ENTRY POINT. GET CORRECT SUB-INTERVAL NUMBER FOR THIS SET
C OF COEFFICIENTS AND THEN GET NORMALIZED CHEBYSHEV TIME
C WITHIN THAT SUBINTERVAL.
C
DNA=DBLE(NA)
DT1=DINT(T(1))
TEMP=DNA*T(1)
L=IDINT(TEMP-DT1)+1
C TC IS THE NORMALIZED CHEBYSHEV TIME (-1 .LE. TC .LE. 1)
TC=2.D0*(DMOD(TEMP,1.D0)+DT1)-1.D0
C CHECK TO SEE WHETHER CHEBYSHEV TIME HAS CHANGED,
C AND COMPUTE NEW POLYNOMIAL VALUES IF IT HAS.
C (THE ELEMENT PC(2) IS THE VALUE OF T1(TC) AND HENCE
C CONTAINS THE VALUE OF TC ON THE PREVIOUS CALL.)
IF(TC.NE.PC(2)) THEN
NP=2
NV=3
PC(2)=TC
TWOT=TC+TC
ENDIF
C
C BE SURE THAT AT LEAST 'NCF' POLYNOMIALS HAVE BEEN EVALUATED
C AND ARE STORED IN THE ARRAY 'PC'.
C
IF(NP.LT.NCF) THEN
DO 1 I=NP+1,NCF
PC(I)=TWOT*PC(I-1)-PC(I-2)
1 CONTINUE
NP=NCF
ENDIF
C
C INTERPOLATE TO GET POSITION FOR EACH COMPONENT
C
DO 2 I=1,NCM
PV(I,1)=0.D0
DO 3 J=NCF,1,-1
PV(I,1)=PV(I,1)+PC(J)*BUF(J,I,L)
3 CONTINUE
2 CONTINUE
IF(IFL.LE.1) RETURN
C
C IF VELOCITY INTERPOLATION IS WANTED, BE SURE ENOUGH
C DERIVATIVE POLYNOMIALS HAVE BEEN GENERATED AND STORED.
C
VFAC=(DNA+DNA)/T(2)
VC(3)=TWOT+TWOT
IF(NV.LT.NCF) THEN
DO 4 I=NV+1,NCF
VC(I)=TWOT*VC(I-1)+PC(I-1)+PC(I-1)-VC(I-2)
4 CONTINUE
NV=NCF
ENDIF
C
C INTERPOLATE TO GET VELOCITY FOR EACH COMPONENT
C
DO 5 I=1,NCM
PV(I,2)=0.D0
DO 6 J=NCF,2,-1
PV(I,2)=PV(I,2)+VC(J)*BUF(J,I,L)
6 CONTINUE
PV(I,2)=PV(I,2)*VFAC
5 CONTINUE
C
RETURN
C
END
C+++++++++++++++++++++++++
C
SUBROUTINE SPLIT(TT,FR)
C
C+++++++++++++++++++++++++
C
C THIS SUBROUTINE BREAKS A D.P. NUMBER INTO A D.P. INTEGER
C AND A D.P. FRACTIONAL PART.
C
C CALLING SEQUENCE PARAMETERS:
C
C TT = D.P. INPUT NUMBER
C
C FR = D.P. 2-WORD OUTPUT ARRAY.
C FR(1) CONTAINS INTEGER PART
C FR(2) CONTAINS FRACTIONAL PART
C
C FOR NEGATIVE INPUT NUMBERS, FR(1) CONTAINS THE NEXT
C MORE NEGATIVE INTEGER; FR(2) CONTAINS A POSITIVE FRACTION.
C
C CALLING SEQUENCE DECLARATIONS
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION FR(2)
C MAIN ENTRY -- GET INTEGER AND FRACTIONAL PARTS
FR(1)=DINT(TT)
FR(2)=TT-FR(1)
IF(TT.GE.0.D0 .OR. FR(2).EQ.0.D0) RETURN
C MAKE ADJUSTMENTS FOR NEGATIVE INPUT NUMBER
FR(1)=FR(1)-1.D0
FR(2)=FR(2)+1.D0
RETURN
END
C++++++++++++++++++++++++++++++++
C
SUBROUTINE STATE(ET2,LIST,PV,PNUT)
C
C++++++++++++++++++++++++++++++++
C
C THIS SUBROUTINE READS AND INTERPOLATES THE JPL PLANETARY EPHEMERIS FILE
C
C CALLING SEQUENCE PARAMETERS:
C
C INPUT:
C
C ET2 DP 2-WORD JULIAN EPHEMERIS EPOCH AT WHICH INTERPOLATION
C IS WANTED. ANY COMBINATION OF ET2(1)+ET2(2) WHICH FALLS
C WITHIN THE TIME SPAN ON THE FILE IS A PERMISSIBLE EPOCH.
C
C A. FOR EASE IN PROGRAMMING, THE USER MAY PUT THE
C ENTIRE EPOCH IN ET2(1) AND SET ET2(2)=0.
C
C B. FOR MAXIMUM INTERPOLATION ACCURACY, SET ET2(1) =
C THE MOST RECENT MIDNIGHT AT OR BEFORE INTERPOLATION
C EPOCH AND SET ET2(2) = FRACTIONAL PART OF A DAY
C ELAPSED BETWEEN ET2(1) AND EPOCH.
C
C C. AS AN ALTERNATIVE, IT MAY PROVE CONVENIENT TO SET
C ET2(1) = SOME FIXED EPOCH, SUCH AS START OF INTEGRATION,
C AND ET2(2) = ELAPSED INTERVAL BETWEEN THEN AND EPOCH.
C
C LIST 12-WORD INTEGER ARRAY SPECIFYING WHAT INTERPOLATION
C IS WANTED FOR EACH OF THE BODIES ON THE FILE.
C
C LIST(I)=0, NO INTERPOLATION FOR BODY I
C =1, POSITION ONLY
C =2, POSITION AND VELOCITY
C
C THE DESIGNATION OF THE ASTRONOMICAL BODIES BY I IS:
C
C I = 1: MERCURY
C = 2: VENUS
C = 3: EARTH-MOON BARYCENTER
C = 4: MARS
C = 5: JUPITER
C = 6: SATURN
C = 7: URANUS
C = 8: NEPTUNE
C = 9: PLUTO
C =10: GEOCENTRIC MOON
C =11: NUTATIONS IN LONGITUDE AND OBLIQUITY
C =12: LUNAR LIBRATIONS (IF ON FILE)
C
C OUTPUT:
C
C PV DP 6 X 11 ARRAY THAT WILL CONTAIN REQUESTED INTERPOLATED
C QUANTITIES (OTHER THAN NUTATION, STOERD IN PNUT).
C THE BODY SPECIFIED BY LIST(I) WILL HAVE ITS
C STATE IN THE ARRAY STARTING AT PV(1,I).
C (ON ANY GIVEN CALL, ONLY THOSE WORDS IN 'PV' WHICH ARE
C AFFECTED BY THE FIRST 10 'LIST' ENTRIES, AND BY LIST(12)
C IF LIBRATIONS ARE ON THE FILE, ARE SET.
C THE REST OF THE 'PV' ARRAYIS UNTOUCHED.)
C THE ORDER OF COMPONENTS STARTING IN PV(1,I) IS: X,Y,Z,DX,DY,DZ.
C
C ALL OUTPUT VECTORS ARE REFERENCED TO THE EARTH MEAN
C EQUATOR AND EQUINOX OF J2000 IF THE DE NUMBER IS 200 OR
C GREATER; OF B1950 IF THE DE NUMBER IS LESS THAN 200.
C
C THE MOON STATE IS ALWAYS GEOCENTRIC; THE OTHER NINE STATES
C ARE EITHER HELIOCENTRIC OR SOLAR-SYSTEM BARYCENTRIC,
C DEPENDING ON THE SETTING OF COMMON FLAGS (SEE BELOW).
C
C LUNAR LIBRATIONS, IF ON FILE, ARE PUT INTO PV(K,11) IF
C LIST(12) IS 1 OR 2.
C
C NUT DP 4-WORD ARRAY THAT WILL CONTAIN NUTATIONS AND RATES,
C DEPENDING ON THE SETTING OF LIST(11). THE ORDER OF
C QUANTITIES IN NUT IS:
C
C D PSI (NUTATION IN LONGITUDE)
C D EPSILON (NUTATION IN OBLIQUITY)
C D PSI DOT
C D EPSILON DOT
C
C * STATEMENT # FOR ERROR RETURN, IN CASE OF EPOCH OUT OF
C RANGE OR I/O ERRORS.
C
C COMMON AREA STCOMX:
C
C KM LOGICAL FLAG DEFINING PHYSICAL UNITS OF THE OUTPUT
C STATES. KM = .TRUE., KM AND KM/SEC
C = .FALSE., AU AND AU/DAY
C DEFAULT VALUE = .FALSE. (KM DETERMINES TIME UNIT
C FOR NUTATIONS AND LIBRATIONS. ANGLE UNIT IS ALWAYS RADIANS.)
C
C BARY LOGICAL FLAG DEFINING OUTPUT CENTER.
C ONLY THE 9 PLANETS ARE AFFECTED.
C BARY = .TRUE. =\ CENTER IS SOLAR-SYSTEM BARYCENTER
C = .FALSE. =\ CENTER IS SUN
C DEFAULT VALUE = .FALSE.
C
C PVSUN DP 6-WORD ARRAY CONTAINING THE BARYCENTRIC POSITION AND
C VELOCITY OF THE SUN.
C
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
SAVE
INTEGER OLDMAX
PARAMETER ( OLDMAX = 400)
INTEGER NMAX
PARAMETER ( NMAX = 1000)
DIMENSION ET2(2),PV(6,11),PNUT(4),T(2),PJD(4),BUF(1500),
. SS(3),CVAL(NMAX),PVSUN(6)
INTEGER LIST(12),IPT(3,13)
LOGICAL FIRST
DATA FIRST/.TRUE./
CHARACTER*6 TTL(14,3),CNAM(NMAX)
CHARACTER*80 NAMFIL
LOGICAL KM,BARY
COMMON/EPHHDR/CVAL,SS,AU,EMRAT,NUMDE,NCON,IPT
COMMON/CHRHDR/CNAM,TTL
COMMON/STCOMX/KM,BARY,PVSUN
C
C ENTRY POINT - 1ST TIME IN, GET POINTER DATA, ETC., FROM EPH FILE
C
IF(FIRST) THEN
FIRST=.FALSE.
C ************************************************************************
C ************************************************************************
C THE USER MUST SELECT ONE OF THE FOLLOWING BY DELETING THE 'C' IN COLUMN 1
C ************************************************************************
C CALL FSIZER1(NRECL,KSIZE,NRFILE,NAMFIL)
C CALL FSIZER2(NRECL,KSIZE,NRFILE,NAMFIL)
CALL FSIZER3(NRECL,KSIZE,NRFILE,NAMFIL)
IF(NRECL .EQ. 0) WRITE(*,*)' ***** FSIZER IS NOT WORKING *****'
C ************************************************************************
C ************************************************************************
IRECSZ=NRECL*KSIZE
NCOEFFS=KSIZE/2
OPEN(NRFILE,
* FILE=NAMFIL,
* ACCESS='DIRECT',
* FORM='UNFORMATTED',
* RECL=IRECSZ,
* STATUS='OLD')
READ(NRFILE,REC=1)TTL,(CNAM(K),K=1,OLDMAX),SS,NCON,AU,EMRAT,
& ((IPT(I,J),I=1,3),J=1,12),NUMDE,(IPT(I,13),I=1,3)
& ,(CNAM(L),L=OLDMAX+1,NCON)
IF(NCON .LE. OLDMAX)THEN
READ(NRFILE,REC=2)(CVAL(I),I=1,OLDMAX)
ELSE
READ(NRFILE,REC=2)(CVAL(I),I=1,NCON)
ENDIF
NRL=0
ENDIF
C ********** MAIN ENTRY POINT **********
IF(ET2(1) .EQ. 0.D0) RETURN
S=ET2(1)-.5D0
CALL SPLIT(S,PJD(1))
CALL SPLIT(ET2(2),PJD(3))
PJD(1)=PJD(1)+PJD(3)+.5D0
PJD(2)=PJD(2)+PJD(4)
CALL SPLIT(PJD(2),PJD(3))
PJD(1)=PJD(1)+PJD(3)
C ERROR RETURN FOR EPOCH OUT OF RANGE
IF(PJD(1)+PJD(4).LT.SS(1) .OR. PJD(1)+PJD(4).GT.SS(2)) GO TO 98
C CALCULATE RECORD # AND RELATIVE TIME IN INTERVAL
NR=IDINT((PJD(1)-SS(1))/SS(3))+3
IF(PJD(1).EQ.SS(2)) NR=NR-1
tmp1 = DBLE(NR-3)*SS(3) + SS(1)
tmp2 = PJD(1) - tmp1
T(1) = (tmp2 + PJD(4))/SS(3)
C READ CORRECT RECORD IF NOT IN CORE
IF(NR.NE.NRL) THEN
NRL=NR
READ(NRFILE,REC=NR,ERR=99)(BUF(K),K=1,NCOEFFS)
ENDIF
IF(KM) THEN
T(2)=SS(3)*86400.D0
AUFAC=1.D0
ELSE
T(2)=SS(3)
AUFAC=1.D0/AU
ENDIF
C INTERPOLATE SSBARY SUN
CALL INTERP(BUF(IPT(1,11)),T,IPT(2,11),3,IPT(3,11),2,PVSUN)
DO I=1,6
PVSUN(I)=PVSUN(I)*AUFAC
ENDDO
C CHECK AND INTERPOLATE WHICHEVER BODIES ARE REQUESTED
DO 4 I=1,10
IF(LIST(I).EQ.0) GO TO 4
CALL INTERP(BUF(IPT(1,I)),T,IPT(2,I),3,IPT(3,I),
& LIST(I),PV(1,I))
DO J=1,6
IF(I.LE.9 .AND. .NOT.BARY) THEN
PV(J,I)=PV(J,I)*AUFAC-PVSUN(J)
ELSE
PV(J,I)=PV(J,I)*AUFAC
ENDIF
ENDDO
4 CONTINUE
C DO NUTATIONS IF REQUESTED (AND IF ON FILE)
IF(LIST(11).GT.0 .AND. IPT(2,12).GT.0)
* CALL INTERP(BUF(IPT(1,12)),T,IPT(2,12),2,IPT(3,12),
* LIST(11),PNUT)
C GET LIBRATIONS IF REQUESTED (AND IF ON FILE)
IF(LIST(12).GT.0 .AND. IPT(2,13).GT.0)
* CALL INTERP(BUF(IPT(1,13)),T,IPT(2,13),3,IPT(3,13),
* LIST(12),PV(1,11))
RETURN
98 WRITE(*,198)ET2(1)+ET2(2),SS(1),SS(2)
198 FORMAT(' *** Requested JED,',f12.2,
* ' not within ephemeris limits,',2f12.2,' ***')
STOP
99 WRITE(*,'(2F12.2,A80)')ET2,'ERROR RETURN IN STATE'
STOP
END
C+++++++++++++++++++++++++++++
C
SUBROUTINE CONST(NAM,VAL,SSS,N)
C
C+++++++++++++++++++++++++++++
C
C THIS ENTRY OBTAINS THE CONSTANTS FROM THE EPHEMERIS FILE
C
C CALLING SEQEUNCE PARAMETERS (ALL OUTPUT):
C
C NAM = CHARACTER*6 ARRAY OF CONSTANT NAMES
C
C VAL = D.P. ARRAY OF VALUES OF CONSTANTS
C
C SSS = D.P. JD START, JD STOP, STEP OF EPHEMERIS
C
C N = INTEGER NUMBER OF ENTRIES IN 'NAM' AND 'VAL' ARRAYS
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
SAVE
INTEGER NMAX
PARAMETER (NMAX = 1000)
CHARACTER*6 NAM(*),TTL(14,3),CNAM(NMAX)
DOUBLE PRECISION VAL(*),SSS(3),SS(3),CVAL(NMAX),ZIPS(2)
DOUBLE PRECISION PVST(6,11),PNUT(4)
DATA ZIPS/2*0.d0/
INTEGER IPT(3,13),DENUM,LIST(12)
logical first
data first/.true./
COMMON/EPHHDR/CVAL,SS,AU,EMRAT,DENUM,NCON,IPT
COMMON/CHRHDR/CNAM,TTL
C CALL STATE TO INITIALIZE THE EPHEMERIS AND READ IN THE CONSTANTS
IF(FIRST) CALL STATE(ZIPS,LIST,PVST,PNUT)
first=.false.
N=NCON
DO I=1,3
SSS(I)=SS(I)
ENDDO
DO I=1,N
NAM(I)=CNAM(I)
VAL(I)=CVAL(I)
ENDDO
RETURN
END

38
lib/libration.f90 Normal file
View File

@ -0,0 +1,38 @@
subroutine libration(jd,RA,Dec,xl,b)
! Compute optical libration of moon at jd: that is, the sub-observer
! point (xl,b) in selenographic coordinates. RA and Dec are
! topocentric values.
implicit real*8 (a-h,o-z)
parameter (RADS=0.0174532925199433d0)
parameter (TWOPI=6.28318530717959d0)
real*8 jd,j2000,mjd,lambda
j2000=2451545.0d0
RA2000=RA
Dec2000=Dec
year=2000.0d0 + (jd-j2000)/365.25d0
mjd=jd-2400000.d0
call sla_PRECES('FK5',year,2000.d0,RA2000,Dec2000)
call sla_EQECL(RA2000,Dec2000,mjd,lambda,beta)
day=jd - j2000
t = day / 36525.d0
xi = 1.54242 * RADS
ft = 93.2720993 + 483202.0175273 * t - .0034029 * t * t
b= ft / 360
a = 360 * (b - floor(b))
if (a.lt.0.d0) a = 360 + a;
f=a/57.2957795131d0
omega=sla_dranrm(2.182439196d0 - t*33.7570446085d0 + t*t*3.6236526d-5)
w = lambda - omega
y = sin(w) * cos(beta) * cos(xi) - sin(beta) * sin(xi)
x = cos(w) * cos(beta)
a = sla_dranrm(atan2(y, x))
xl = a - f
if(xl.lt.-0.25*TWOPI) xl=xl+TWOPI !Fix 'round the back' angles
if(xl.gt.0.25*TWOPI) xl=xl-TWOPI
b = asin(-sin(w) * cos(beta) * sin(xi) - sin(beta) * cos(xi))
return
end subroutine libration

166
lib/moon2.f90
View File

@ -1,166 +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
! 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

72
lib/moondop.f90
View File

@ -1,72 +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 rme0(6)
logical km
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

43
lib/moondopjpl.f90 Normal file
View File

@ -0,0 +1,43 @@
subroutine MoonDopJPL(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 !East 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
twopi=8.d0*atan(1.d0) !Define some constants
rad=360.d0/twopi
clight=2.99792458d5
call sla_CLDJ(nyear,month,nday,djutc,j)
djutc=djutc + uth4/24.d0
dut=-0.460d0
east_long=lon4/rad
geodetic_lat=lat4/rad
height=40.
nspecial=0
call ephem(djutc,dut,east_long,geodetic_lat,height,nspecial, &
RA,Dec,Az,El,techo,dop,fspread_1GHz,vr)
RAMoon4=RA
DecMoon4=Dec
LST4=LST
HA4=HA
AzMoon4=Az*rad
ElMoon4=El*rad
vr4=vr
dist4=techo
return
end subroutine MoonDopJPL

3396
lib/slasubs.f Normal file
View File

File diff suppressed because it is too large Load Diff

14
lib/tm2.f90
View File

@ -1,14 +0,0 @@
subroutine tm2(day,xlat4,xlon4,xl4,b4)
implicit real*8 (a-h,o-z)
parameter (RADS=0.0174532925199433d0)
real*4 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
View File

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