mirror of
https://github.com/saitohirga/WSJT-X.git
synced 2024-11-22 20:28:42 -05:00
More conversions of .f to .f90.
git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/map65@7474 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
This commit is contained in:
Joe Taylor
parent
9dd0cddffa
commit
e5c1c14543
@ -88,11 +88,19 @@ set (FSRCS
|
|||||||
getdphi.f90
|
getdphi.f90
|
||||||
getpfx1.f90
|
getpfx1.f90
|
||||||
getpfx2.f90
|
getpfx2.f90
|
||||||
|
graycode.f90
|
||||||
graycode65.f90
|
graycode65.f90
|
||||||
|
grid2deg.f90
|
||||||
|
grid2k.f90
|
||||||
|
indexx.f90
|
||||||
|
interleave63.f90
|
||||||
iqcal.f90
|
iqcal.f90
|
||||||
iqfix.f90
|
iqfix.f90
|
||||||
jt65code.f90
|
jt65code.f90
|
||||||
|
k2grid.f90
|
||||||
map65a.f90
|
map65a.f90
|
||||||
|
moon2.f90
|
||||||
|
moondop.f90
|
||||||
noisegen.f90
|
noisegen.f90
|
||||||
pfxdump.f90
|
pfxdump.f90
|
||||||
recvpkt.f90
|
recvpkt.f90
|
||||||
@ -107,14 +115,6 @@ set (FSRCS
|
|||||||
zplot.f90
|
zplot.f90
|
||||||
|
|
||||||
f77_wisdom.f
|
f77_wisdom.f
|
||||||
graycode.f
|
|
||||||
grid2deg.f
|
|
||||||
grid2k.f
|
|
||||||
indexx.f
|
|
||||||
interleave63.f
|
|
||||||
k2grid.f
|
|
||||||
moon2.f
|
|
||||||
moondop.f
|
|
||||||
nchar.f
|
nchar.f
|
||||||
packcall.f
|
packcall.f
|
||||||
packdxcc.f
|
packdxcc.f
|
||||||
|
|||||||
@ -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
|
|
||||||
|
|
||||||
10
libm65/graycode.f90
Normal file
10
libm65/graycode.f90
Normal file
@ -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
|
||||||
|
|
||||||
@ -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
|
|
||||||
38
libm65/grid2deg.f90
Normal file
38
libm65/grid2deg.f90
Normal file
@ -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
|
||||||
@ -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
|
|
||||||
12
libm65/grid2k.f90
Normal file
12
libm65/grid2k.f90
Normal file
@ -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
|
||||||
@ -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
|
|
||||||
|
|
||||||
19
libm65/indexx.f90
Normal file
19
libm65/indexx.f90
Normal file
@ -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
|
||||||
|
|
||||||
@ -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
|
|
||||||
25
libm65/interleave63.f90
Normal file
25
libm65/interleave63.f90
Normal file
@ -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
|
||||||
@ -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
|
|
||||||
12
libm65/k2grid.f90
Normal file
12
libm65/k2grid.f90
Normal file
@ -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
|
||||||
167
libm65/moon2.f
167
libm65/moon2.f
@ -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
|
|
||||||
163
libm65/moon2.f90
Normal file
163
libm65/moon2.f90
Normal file
@ -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
|
||||||
@ -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
|
|
||||||
73
libm65/moondop.f90
Normal file
73
libm65/moondop.f90
Normal file
@ -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
|
||||||
Loading…
Reference in New Issue
Block a user