w4kek
/
WSJT-X
mirror of https://github.com/saitohirga/WSJT-X.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Remove stuff not used in WSJT-X. 

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@2634 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
This commit is contained in:
Joe Taylor 2012-10-03 14:17:08 +00:00
parent 7d28a1883b
commit f435b4028b
62 changed files with 2 additions and 4662 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

47
libm65/JT65code.f90
View File

@ -1,47 +0,0 @@
program JT65code
! Provides examples of message packing, bit and symbol ordering,
! Reed Solomon encoding, and other necessary details of the JT65
! protocol.
character*22 msg0,msg,decoded,cok*3
integer dgen(12),sent(63),recd(12),era(51)
nargs=iargc()
if(nargs.ne.1) then
print*,'Usage: JT65code "message"'
go to 999
endif
call getarg(1,msg0) !Get message from command line
msg=msg0
call chkmsg(msg,cok,nspecial,flip) !See if it includes "OOO" report
if(nspecial.gt.0) then !or is a shorthand message
write(*,1010)
1010 format('Shorthand message.')
go to 999
endif
call packmsg(msg,dgen) !Pack message into 72 bits
write(*,1020) msg0
1020 format('Message: ',a22) !Echo input message
if(iand(dgen(10),8).ne.0) write(*,1030) !Is plain text bit set?
1030 format('Plain text.')
write(*,1040) dgen
1040 format('Packed message, 6-bit symbols: ',12i3) !Display packed symbols
call rs_encode(dgen,sent) !RS encode
call interleave63(sent,1) !Interleave channel symbols
call graycode(sent,63,1) !Apply Gray code
write(*,1050) sent
1050 format('Channel symbols, including FEC:'/(i5,20i3))
call graycode(sent,63,-1)
call interleave63(sent,-1)
call rs_decode(sent,era,0,recd,nerr)
call unpackmsg(recd,decoded) !Unpack the user message
write(*,1060) decoded,cok
1060 format('Decoded message: ',a22,2x,a3)
999 end program JT65code

68
libm65/afc65b.f
View File

@ -1,68 +0,0 @@
subroutine afc65b(cx,cy,npts,fsample,nflip,ipol,xpol,a,
+ ccfbest,dtbest)
logical xpol
complex cx(npts)
complex cy(npts)
real a(5),deltaa(5)
a(1)=0.
a(2)=0.
a(3)=0.
a(4)=45.0*(ipol-1.0)
deltaa(1)=2.0
deltaa(2)=2.0
deltaa(3)=2.0
deltaa(4)=22.5
deltaa(5)=0.05
nterms=3
if(xpol) nterms=4
chisqr=0.
C Start the iteration
chisqr0=1.e6
do iter=1,3 !One iteration is enough?
do j=1,nterms
chisq1=fchisq(cx,cy,npts,fsample,nflip,a,ccfmax,dtmax)
fn=0.
delta=deltaa(j)
10 a(j)=a(j)+delta
chisq2=fchisq(cx,cy,npts,fsample,nflip,a,ccfmax,dtmax)
if(chisq2.eq.chisq1) go to 10
if(chisq2.gt.chisq1) then
delta=-delta !Reverse direction
a(j)=a(j)+delta
tmp=chisq1
chisq1=chisq2
chisq2=tmp
endif
20 fn=fn+1.0
a(j)=a(j)+delta
chisq3=fchisq(cx,cy,npts,fsample,nflip,a,ccfmax,dtmax)
if(chisq3.lt.chisq2) then
chisq1=chisq2
chisq2=chisq3
go to 20
endif
C Find minimum of parabola defined by last three points
delta=delta*(1./(1.+(chisq1-chisq2)/(chisq3-chisq2))+0.5)
a(j)=a(j)-delta
deltaa(j)=deltaa(j)*fn/3.
enddo
chisqr=fchisq(cx,cy,npts,fsample,nflip,a,ccfmax,dtmax)
if(chisqr/chisqr0.gt.0.9999) go to 30
chisqr0=chisqr
enddo
30 ccfbest=ccfmax * (1378.125/fsample)**2
dtbest=dtmax
if(a(4).lt.0.0) a(4)=a(4)+180.0
if(a(4).ge.180.0) a(4)=a(4)-180.0
if(nint(a(4)).eq.180) a(4)=0.
ipol=nint(a(4)/45.0) + 1
if(ipol.gt.4) ipol=ipol-4
return
end

38
libm65/alignmsg.f90
View File

@ -1,38 +0,0 @@
subroutine alignmsg(word0,nmin,msg,msglen,idone)
character*(*) word0
character*29 msg,word
word=word0//' '
idone=0
! Test for two (or more) <space> characters
if(word(1:2).eq.' ' .and. len(word).eq.2) then
i2=index(msg,' ')
if((i2.ge.1.and.i2.lt.msglen) .or. &
(msg(1:1).eq.' '.and.msg(msglen:msglen).eq.' ')) then
if(i2.eq.1) msg=msg(i2+2:msglen) !Align on EOM
if(i2.ge.2) msg=msg(i2+2:msglen)//msg(1:i2-1)
idone=1
endif
! Align on single <space> (as last resort)
else if(word(1:1).eq.' ' .and. len(word).eq.1) then
i3=index(msg,' ')
if(i3.ge.1 .and. i3.lt.msglen) msg=msg(i3+1:msglen)//msg(1:i3)
if(i3.eq.msglen) msg=msg(1:msglen)
msg=msg(1:msglen)//msg(1:msglen)
idone=1
! Align on specified word
else
call match(word,msg(1:msglen),nstart,nmatch)
if(nmatch.ge.nmin) then
if(nstart.eq.1) msg=msg(nstart:msglen)
if(nstart.gt.1) msg=msg(nstart:msglen)//msg(1:nstart-1)
idone=1
endif
endif
return
end subroutine alignmsg

109
libm65/astro.f
View File

@ -1,109 +0,0 @@
subroutine astro(nyear,month,nday,uth,nfreq,Mygrid,
+ NStation,MoonDX,AzSun,ElSun,AzMoon0,ElMoon0,
+ ntsky,doppler00,doppler,dbMoon,RAMoon,DecMoon,HA,Dgrd,sd,
+ poloffset,xnr,day,lon,lat,LST)
C Computes astronomical quantities for display and tracking.
C NB: may want to smooth the Tsky map to 10 degrees or so.
character*6 MyGrid,HisGrid
real LST
real lat,lon
integer*2 nt144(180)
! common/echo/xdop(2),techo,AzMoon,ElMoon,mjd
real xdop(2)
data rad/57.2957795/
data nt144/
+ 234, 246, 257, 267, 275, 280, 283, 286, 291, 298,
+ 305, 313, 322, 331, 341, 351, 361, 369, 376, 381,
+ 383, 382, 379, 374, 370, 366, 363, 361, 363, 368,
+ 376, 388, 401, 415, 428, 440, 453, 467, 487, 512,
+ 544, 579, 607, 618, 609, 588, 563, 539, 512, 482,
+ 450, 422, 398, 379, 363, 349, 334, 319, 302, 282,
+ 262, 242, 226, 213, 205, 200, 198, 197, 196, 197,
+ 200, 202, 204, 205, 204, 203, 202, 201, 203, 206,
+ 212, 218, 223, 227, 231, 236, 240, 243, 247, 257,
+ 276, 301, 324, 339, 346, 344, 339, 331, 323, 316,
+ 312, 310, 312, 317, 327, 341, 358, 375, 392, 407,
+ 422, 437, 451, 466, 480, 494, 511, 530, 552, 579,
+ 612, 653, 702, 768, 863,1008,1232,1557,1966,2385,
+ 2719,2924,3018,3038,2986,2836,2570,2213,1823,1461,
+ 1163, 939, 783, 677, 602, 543, 494, 452, 419, 392,
+ 373, 360, 353, 350, 350, 350, 350, 350, 350, 348,
+ 344, 337, 329, 319, 307, 295, 284, 276, 272, 272,
+ 273, 274, 274, 271, 266, 260, 252, 245, 238, 231/
save
call grid2deg(MyGrid,elon,lat)
lon=-elon
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,
+ LST,HA,AzMoon,ElMoon,vr,dist)
C Compute spatial polarization offset
xx=sin(lat/rad)*cos(ElMoon/rad) - cos(lat/rad)*
+ cos(AzMoon/rad)*sin(ElMoon/rad)
yy=cos(lat/rad)*sin(AzMoon/rad)
if(NStation.eq.1) poloffset1=rad*atan2(yy,xx)
if(NStation.eq.2) poloffset2=rad*atan2(yy,xx)
techo=2.0 * dist/2.99792458e5 !Echo delay time
doppler=-freq*vr/2.99792458e5 !One-way Doppler
call coord(0.,0.,-1.570796,1.161639,RAMoon/rad,DecMoon/rad,el,eb)
longecl_half=nint(rad*el/2.0)
if(longecl_half.lt.1 .or. longecl_half.gt.180) longecl_half=180
t144=nt144(longecl_half)
tsky=(t144-2.7)*(144.0/nfreq)**2.6 + 2.7 !Tsky for obs freq
xdop(NStation)=doppler
if(NStation.eq.2) then
HisGrid=MyGrid
go to 900
endif
doppler00=2.0*xdop(1)
doppler=xdop(1)+xdop(2)
! if(mode.eq.3) doppler=2.0*xdop(1)
dBMoon=-40.0*log10(dist/356903.)
sd=16.23*370152.0/dist
! if(NStation.eq.1 .and. MoonDX.ne.0 .and.
! + (mode.eq.2 .or. mode.eq.5)) then
if(NStation.eq.1 .and. MoonDX.ne.0) then
poloffset=mod(poloffset2-poloffset1+720.0,180.0)
if(poloffset.gt.90.0) poloffset=poloffset-180.0
x1=abs(cos(2*poloffset/rad))
if(x1.lt.0.056234) x1=0.056234
xnr=-20.0*log10(x1)
if(HisGrid(1:1).lt.'A' .or. HisGrid(1:1).gt.'R') xnr=0
endif
tr=80.0 !Good preamp
tskymin=13.0*(408.0/nfreq)**2.6 !Cold sky temperature
tsysmin=tskymin+tr
tsys=tsky+tr
dgrd=-10.0*log10(tsys/tsysmin) + dbMoon
900 AzMoon0=Azmoon
ElMoon0=Elmoon
ntsky=nint(tsky)
! auxHA = 15.0*(LST-auxra) !HA in degrees
! pi=3.14159265
! pio2=0.5*pi
! call coord(pi,pio2-lat/rad,0.0,lat/rad,auxha*pi/180.0,
! + auxdec/rad,azaux,elaux)
! AzAux=azaux*rad
! ElAux=ElAux*rad
return
end

81
libm65/astro0.f90
View File

@ -1,81 +0,0 @@
subroutine astro0(nyear,month,nday,uth8,nfreq,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)
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 uth8
data uth8z/0.d0/
save
uth=uth8
call astro(nyear,month,nday,uth,nfreq,hisgrid,2,1, &
AzSun,ElSun,AzMoon,ElMoon,ntsky,doppler00,doppler, &
dbMoon,RAMoon,DecMoon,HA,Dgrd,sd,poloffset,xnr, &
day,xlon2,xlat2,xlst)
AzMoonB8=AzMoon
ElMoonB8=ElMoon
call astro(nyear,month,nday,uth,nfreq,mygrid,1,1, &
AzSun,ElSun,AzMoon,ElMoon,ntsky,doppler00,doppler, &
dbMoon,RAMoon,DecMoon,HA,Dgrd,sd,poloffset,xnr, &
day,xlon1,xlat1,xlst)
day8=day
xlst8=xlst
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=0.001*nfreq
dldt1=DEGS*(xl1a-xl1)
dbdt1=DEGS*(b1a-b1)
dldt2=DEGS*(xl2a-xl2)
dbdt2=DEGS*(b2a-b2)
rate1=2.0*sqrt(dldt1**2 + dbdt1**2)
width1=0.5*6741*fghz*rate1
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
ElMoon8=ElMoon
dbMoon8=dbMoon
RAMoon8=RAMoon/15.0
DecMoon8=DecMoon
HA8=HA
Dgrd8=Dgrd
sd8=sd
poloffset8=poloffset
xnr8=xnr
ndop=nint(doppler)
ndop00=nint(doppler00)
if(uth8z.eq.0.d0) then
uth8z=uth8-1.d0/3600.d0
dopplerz=doppler
doppler00z=doppler00
endif
dt=60.0*(uth8-uth8z)
if(dt.le.0) dt=1.d0/60.d0
dfdt=(doppler-dopplerz)/dt
dfdt0=(doppler00-doppler00z)/dt
uth8z=uth8
dopplerz=doppler
doppler00z=doppler00
return
end subroutine astro0

14
libm65/astrosub.f90
View File

@ -1,14 +0,0 @@
subroutine astrosub(nyear,month,nday,uth8,nfreq,mygrid,hisgrid, &
AzSun8,ElSun8,AzMoon8,ElMoon8,AzMoonB8,ElMoonB8,ntsky,ndop,ndop00, &
RAMoon8,DecMoon8,Dgrd8,poloffset8,xnr8)
implicit real*8 (a-h,o-z)
character*6 mygrid,hisgrid
call astro0(nyear,month,nday,uth8,nfreq,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)
return
end subroutine astrosub

45
libm65/ccf2.f
View File

@ -1,45 +0,0 @@
subroutine ccf2(ss,nz,nflip,ccfbest,lagpk)
parameter (LAGMAX=60)
! parameter (LAGMAX=200)
real ss(nz)
real ccf(-LAGMAX:LAGMAX)
integer npr(126)
C The JT65 pseudo-random sync pattern:
data npr/
+ 1,0,0,1,1,0,0,0,1,1,1,1,1,1,0,1,0,1,0,0,
+ 0,1,0,1,1,0,0,1,0,0,0,1,1,1,0,0,1,1,1,1,
+ 0,1,1,0,1,1,1,1,0,0,0,1,1,0,1,0,1,0,1,1,
+ 0,0,1,1,0,1,0,1,0,1,0,0,1,0,0,0,0,0,0,1,
+ 1,0,0,0,0,0,0,0,1,1,0,1,0,0,1,0,1,1,0,1,
+ 0,1,0,1,0,0,1,1,0,0,1,0,0,1,0,0,0,0,1,1,
+ 1,1,1,1,1,1/
save
ccfbest=0.
lag1=-LAGMAX
lag2=LAGMAX
do lag=lag1,lag2
s0=0.
s1=0.
do i=1,126
j=2*(8*i + 43) + lag
if(j.ge.1 .and. j.le.nz-8) then
x=ss(j)+ss(j+8) !Add two half-symbol contributions
if(npr(i).eq.0) then
s0=s0 + x
else
s1=s1 + x
endif
endif
enddo
ccf(lag)=nflip*(s1-s0)
if(ccf(lag).gt.ccfbest) then
ccfbest=ccf(lag)
lagpk=lag
endif
enddo
return
end

118
libm65/ccf65.f90
View File

@ -1,118 +0,0 @@
subroutine ccf65(ss,nhsym,ssmax,sync1,ipol1,jpz,dt1,flipk,syncshort, &
snr2,ipol2,dt2)
parameter (NFFT=512,NH=NFFT/2)
real ss(4,322) !Input: half-symbol powers, 4 pol'ns
real s(NFFT) !CCF = ss*pr
complex cs(0:NH) !Complex FT of s
real s2(NFFT) !CCF = ss*pr2
complex cs2(0:NH) !Complex FT of s2
real pr(NFFT) !JT65 pseudo-random sync pattern
complex cpr(0:NH) !Complex FT of pr
real pr2(NFFT) !JT65 shorthand pattern
complex cpr2(0:NH) !Complex FT of pr2
real tmp1(322)
real tmp2(322)
real ccf(-27:27,4)
logical first
integer npr(126)
data first/.true./
equivalence (s,cs),(pr,cpr),(s2,cs2),(pr2,cpr2)
save
! The JT65 pseudo-random sync pattern:
data npr/ &
1,0,0,1,1,0,0,0,1,1,1,1,1,1,0,1,0,1,0,0, &
0,1,0,1,1,0,0,1,0,0,0,1,1,1,0,0,1,1,1,1, &
0,1,1,0,1,1,1,1,0,0,0,1,1,0,1,0,1,0,1,1, &
0,0,1,1,0,1,0,1,0,1,0,0,1,0,0,0,0,0,0,1, &
1,0,0,0,0,0,0,0,1,1,0,1,0,0,1,0,1,1,0,1, &
0,1,0,1,0,0,1,1,0,0,1,0,0,1,0,0,0,0,1,1, &
1,1,1,1,1,1/
if(first) then
! Initialize pr, pr2; compute cpr, cpr2.
fac=1.0/NFFT
do i=1,NFFT
pr(i)=0.
pr2(i)=0.
k=2*mod((i-1)/8,2)-1
if(i.le.NH) pr2(i)=fac*k
enddo
do i=1,126
j=2*i
pr(j)=fac*(2*npr(i)-1)
! Not sure why, but it works significantly better without the following line:
! pr(j-1)=pr(j)
enddo
call four2a(pr,NFFT,1,-1,0)
call four2a(pr2,NFFT,1,-1,0)
first=.false.
endif
! Look for JT65 sync pattern and shorthand square-wave pattern.
ccfbest=0.
ccfbest2=0.
ipol1=1
ipol2=1
do ip=1,jpz !Do jpz polarizations
do i=1,nhsym-1
! s(i)=ss(ip,i)+ss(ip,i+1)
s(i)=min(ssmax,ss(ip,i)+ss(ip,i+1))
enddo
s(nhsym:NFFT)=0.
call four2a(s,NFFT,1,-1,0) !Real-to-complex FFT
do i=0,NH
cs2(i)=cs(i)*conjg(cpr2(i)) !Mult by complex FFT of pr2
cs(i)=cs(i)*conjg(cpr(i)) !Mult by complex FFT of pr
enddo
call four2a(cs,NFFT,1,1,-1) !Complex-to-real inv-FFT
call four2a(cs2,NFFT,1,1,-1) !Complex-to-real inv-FFT
do lag=-27,27 !Check for best JT65 sync
ccf(lag,ip)=s(lag+28)
if(abs(ccf(lag,ip)).gt.ccfbest) then
ccfbest=abs(ccf(lag,ip))
lagpk=lag
ipol1=ip
flipk=1.0
if(ccf(lag,ip).lt.0.0) flipk=-1.0
endif
enddo
do lag=-8,7 !Check for best shorthand
ccf2=s2(lag+28)
if(ccf2.gt.ccfbest2) then
ccfbest2=ccf2
lagpk2=lag
ipol2=ip
endif
enddo
enddo
! Find rms level on baseline of "ccfblue", for normalization.
sum=0.
do lag=-26,26
if(abs(lag-lagpk).gt.1) sum=sum + ccf(lag,ipol1)
enddo
base=sum/50.0
sq=0.
do lag=-26,26
if(abs(lag-lagpk).gt.1) sq=sq + (ccf(lag,ipol1)-base)**2
enddo
rms=sqrt(sq/49.0)
sync1=ccfbest/rms - 4.0
dt1=2.5 + lagpk*(2048.0/11025.0)
! Find base level for normalizing snr2.
do i=1,nhsym
tmp1(i)=ss(ipol2,i)
enddo
call pctile(tmp1,tmp2,nhsym,40,base)
snr2=0.398107*ccfbest2/base !### empirical
syncshort=0.5*ccfbest2/rms - 4.0 !### better normalizer than rms?
dt2=2.5 + lagpk2*(2048.0/11025.0)
return
end subroutine ccf65

23
libm65/chkhist.f
View File

@ -1,23 +0,0 @@
subroutine chkhist(mrsym,nmax,ipk)
integer mrsym(63)
integer hist(0:63)
do i=0,63
hist(i)=0
enddo
do j=1,63
i=mrsym(j)
hist(i)=hist(i)+1
enddo
nmax=0
do i=0,63
if(hist(i).gt.nmax) then
nmax=hist(i)
ipk=i+1
endif
enddo
return
end

32
libm65/chkmsg.f
View File

@ -1,32 +0,0 @@
subroutine chkmsg(message,cok,nspecial,flip)
character message*22,cok*3
nspecial=0
flip=1.0
cok=" "
do i=22,1,-1
if(message(i:i).ne.' ') go to 10
enddo
i=22
10 if(i.ge.11) then
if ((message(i-3:i).eq.' OOO') .or.
+ (message(20:22).eq.' OO')) then
cok='OOO'
flip=-1.0
if(message(20:22).eq.' OO') then
message=message(1:19)
else
message=message(1:i-4)
endif
endif
endif
if(message(1:3).eq.'RO ') nspecial=2
if(message(1:4).eq.'RRR ') nspecial=3
if(message(1:3).eq.'73 ') nspecial=4
return
end

41
libm65/coord.f
View File

@ -1,41 +0,0 @@
SUBROUTINE COORD(A0,B0,AP,BP,A1,B1,A2,B2)
C Examples:
C 1. From ha,dec to az,el:
C call coord(pi,pio2-lat,0.,lat,ha,dec,az,el)
C 2. From az,el to ha,dec:
C call coord(pi,pio2-lat,0.,lat,az,el,ha,dec)
C 3. From ra,dec to l,b
C call coord(4.635594495,-0.504691042,3.355395488,0.478220215,
C ra,dec,l,b)
C 4. From l,b to ra,dec
C call coord(1.705981071d0,-1.050357016d0,2.146800277d0,
C 0.478220215d0,l,b,ra,dec)
C 5. From ra,dec to ecliptic latitude (eb) and longitude (el):
C call coord(0.d0,0.d0,-pio2,pio2-23.443*pi/180,ra,dec,el,eb)
C 6. From ecliptic latitude (eb) and longitude (el) to ra,dec:
C call coord(0.d0,0.d0,-pio2,pio2-23.443*pi/180,el,eb,ra,dec)
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.e0-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.e0) TA2O2=(1.e0-CA2)/SA2
IF(CA2.GT.0.e0) TA2O2=SA2/(1.e0+CA2)
A2=2.e0*atan(TA2O2)
IF(A2.LT.0.e0) A2=A2+6.2831853
RETURN
END

40
libm65/dcoord.f
View File

@ -1,40 +0,0 @@
SUBROUTINE DCOORD(A0,B0,AP,BP,A1,B1,A2,B2)
implicit real*8 (a-h,o-z)
C Examples:
C 1. From ha,dec to az,el:
C call coord(pi,pio2-lat,0.,lat,ha,dec,az,el)
C 2. From az,el to ha,dec:
C call coord(pi,pio2-lat,0.,lat,az,el,ha,dec)
C 3. From ra,dec to l,b
C call coord(4.635594495,-0.504691042,3.355395488,0.478220215,
C ra,dec,l,b)
C 4. From l,b to ra,dec
C call coord(1.705981071d0,-1.050357016d0,2.146800277d0,
C 0.478220215d0,l,b,ra,dec)
C 5. From ecliptic latitude (eb) and longitude (el) to ra, dec:
C 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

149
libm65/decode1a.f
View File

@ -1,149 +0,0 @@
subroutine decode1a(dd,newdat,f0,nflip,mode65,nfsample,xpol,
+ mycall,hiscall,hisgrid,neme,ndepth,nqd,dphi,
+ nutc,nkhz,ndf,ipol,sync2,a,dt,pol,nkv,nhist,qual,decoded)
! Apply AFC corrections to a candidate JT65 signal, then decode it.
parameter (NMAX=60*96000) !Samples per 60 s
real*4 dd(4,NMAX) !92 MB: raw data from Linrad timf2
complex cx(NMAX/64), cy(NMAX/64) !Data at 1378.125 samples/s
complex c5x(NMAX/256),c5y(NMAX/256),c5tmp(NMAX/256) !Data at 344.53125 Hz
complex c5a(512)
complex z
real s2(66,126)
real s3(64,63),sy(63)
real a(5)
logical first,xpol
character decoded*22
character mycall*12,hiscall*12,hisgrid*6
data first/.true./,jjjmin/1000/,jjjmax/-1000/
data nutc0/-999/,nkhz0/-999/
save
! Mix sync tone to baseband, low-pass filter, downsample to 1378.125 Hz
dt00=dt
call timer('filbig ',0)
call filbig(dd,NMAX,f0,newdat,nfsample,xpol,cx,cy,n5)
! NB: cx, cy have sample rate 96000*77125/5376000 = 1378.125 Hz
call timer('filbig ',1)
joff=0
sqa=0.
sqb=0.
do i=1,n5
sqa=sqa + real(cx(i))**2 + aimag(cx(i))**2
sqb=sqb + real(cy(i))**2 + aimag(cy(i))**2
enddo
sqa=sqa/n5
sqb=sqb/n5
! Find best DF, f1, f2, DT, and pol. Start by downsampling to 344.53125 Hz
z=cmplx(cos(dphi),sin(dphi))
cy(:n5)=z*cy(:n5) !Adjust for cable length difference
call timer('fil6521 ',0)
call fil6521(cx,n5,c5x,n6)
call fil6521(cy,n5,c5y,n6)
call timer('fil6521 ',1)
! Add some zeros at start of c5 arrays -- empirical fix for negative DT's
! NB: might be better to add zeros to cx and cy, rather than here.
! Q: is the DT search range big enough?
nadd=200
c5tmp(1:nadd)=0.
c5tmp(1+nadd:n6+nadd)=c5x(1:n6)
c5x(1:n6+nadd)=c5tmp(1:n6+nadd)
c5tmp(1+nadd:n6+nadd)=c5y(1:n6)
c5y(1:n6+nadd)=c5tmp(1:n6+nadd)
n6=n6+nadd
fsample=1378.125/4.
a(5)=dt00
i0=nint((a(5)+0.5)*fsample) - 2 + 200
if(i0.lt.1) then
write(13,*) 'i0 too small in decode1a:',i0,f0
flush(13)
i0=1
endif
nz=n6+1-i0
! We're looking only at sync tone here... so why not downsample by another
! factor of 1/8, say? Should be a significant execution speed-up.
call timer('afc65b ',0)
! Best fit for DF, f1, f2, pol
call afc65b(c5x(i0),c5y(i0),nz,fsample,nflip,ipol,xpol,a,
+ ccfbest,dtbest)
call timer('afc65b ',1)
pol=a(4)/57.2957795
aa=cos(pol)
bb=sin(pol)
sq0=aa*aa*sqa + bb*bb*sqb
sync2=3.7*ccfbest/sq0
! Apply AFC corrections to the time-domain signal
! Now we are back to using the 1378.125 Hz sample rate, enough to
! accommodate the full JT65C bandwidth.
call timer('twkfreq ',0)
call twkfreq(cx,cy,n5,a)
call timer('twkfreq ',1)
! Compute spectrum at best polarization for each half symbol.
! Adding or subtracting a small number (e.g., 5) to j may make it decode.\
! NB: might want to try computing full-symbol spectra (nfft=512, even for
! submodes B and C).
nsym=126
nfft=512/mode65
j=(dt00+dtbest+2.685)*1378.125 + joff
if(j.lt.0) j=0
call timer('sh_ffts ',0)
! Perhaps should try full-symbol-length FFTs even in B, C sub-modes?
! (Tried this, found no significant difference in decodes.)
do k=1,nsym
do n=1,mode65
do i=1,nfft
j=j+1
c5a(i)=aa*cx(j) + bb*cy(j)
enddo
call four2a(c5a,nfft,1,1,1)
if(n.eq.1) then
do i=1,66
s2(i,k)=real(c5a(i))**2 + aimag(c5a(i))**2
enddo
else
do i=1,66
s2(i,k)=s2(i,k) + real(c5a(i))**2 + aimag(c5a(i))**2
enddo
endif
enddo
enddo
call timer('sh_ffts ',1)
flip=nflip
call timer('dec65b ',0)
call decode65b(s2,flip,mycall,hiscall,hisgrid,mode65,neme,ndepth,
+ nqd,nkv,nhist,qual,decoded,s3,sy)
dt=dt00 + dtbest
call timer('dec65b ',1)
if(nqd.eq.1 .and. nkv.eq.0) then
if(nutc.ne.nutc0) syncbest=0.
if(sync2.gt.syncbest) then
if(nutc.eq.nutc0) nsave=nsave-1
if(nkhz.ne.nkhz0) nsave=0
nkhz0=nkhz
nsave=min(32,nsave+1)
npol=nint(57.296*pol)
call s3avg(nsave,mode65,nutc,ndf,dt+0.8,npol,s3,nkv,decoded)
syncbest=sync2
endif
nutc0=nutc
endif
return
end

49
libm65/decode65b.f
View File

@ -1,49 +0,0 @@
subroutine decode65b(s2,flip,mycall,hiscall,hisgrid,mode65,
+ neme,ndepth,nqd,nkv,nhist,qual,decoded,s3,sy)
real s2(66,126)
real s3(64,63),sy(63)
logical first,ltext
character decoded*22,deepmsg*22
character mycall*12,hiscall*12,hisgrid*6
common/prcom/pr(126),mdat(126),mref(126,2),mdat2(126),mref2(126,2)
data first/.true./
save
if(first) call setup65
first=.false.
do j=1,63
k=mdat(j) !Points to data symbol
if(flip.lt.0.0) k=mdat2(j)
do i=1,64
s3(i,j)=s2(i+2,k)
enddo
k=mdat2(j) !Points to data symbol
if(flip.lt.0.0) k=mdat(j)
sy(j)=s2(1,k)
enddo
nadd=mode65
call extract(s3,nadd,ncount,nhist,decoded,ltext) !Extract the message
C Suppress "birdie messages" and other garbage decodes:
if(decoded(1:7).eq.'000AAA ') ncount=-1
if(decoded(1:7).eq.'0L6MWK ') ncount=-1
if(flip.lt.0.0 .and. ltext) ncount=-1
nkv=1
if(ncount.lt.0) then
nkv=0
decoded=' '
endif
qual=0.
if(ndepth.ge.1 .and. (nqd.eq.1 .or. flip.eq.1.0)) then
call deep65(s3,mode65,neme,flip,mycall,hiscall,
+ hisgrid,deepmsg,qual)
if(nqd.ne.1 .and. qual.lt.10.0) qual=0.0
if(ndepth.lt.2 .and. qual.lt.6.0) qual=0.0
endif
if(nkv.eq.0 .and. qual.ge.1.0) decoded=deepmsg
return
end

39
libm65/decodemsk.f90
View File

@ -1,39 +0,0 @@
subroutine decodemsk(cdat,npts,cw,i1,nchar,s2,msg)
! DF snd sync have been established, now decode the message
complex cdat(npts)
complex cw(168,0:63) !Complex waveforms for codewords
real s2(0:63,400)
character msg*400
complex z,zmax
character cc*64
! 1 2 3 4 5 6
! 0123456789012345678901234567890123456789012345678901234567890123
data cc/'0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ./?- _ @'/
msg=' '
do j=1,nchar !Find best match for each character
ia=i1 + (j-1)*168
smax=0.
do k=0,40
kk=k
if(k.eq.40) kk=44
z=0.
do i=1,168
z=z + cdat(ia+i)*conjg(cw(i,kk))
enddo
ss=abs(z)
s2(k,j)=ss
if(ss.gt.smax) then
smax=ss
zmax=z
kpk=kk
endif
enddo
msg(j:j)=cc(kpk+1:kpk+1)
if(kpk.eq.44) msg(j:j)=' '
enddo
return
end subroutine decodemsk

169
libm65/deep65.f90
View File

@ -1,169 +0,0 @@
subroutine deep65(s3,mode65,neme,flip,mycall,hiscall,hisgrid,decoded,qual)
parameter (MAXCALLS=7000,MAXRPT=63)
real s3(64,63)
character callsign*12,grid*4,message*22,hisgrid*6,c*1,ceme*3
character*12 mycall,hiscall
character*22 decoded
character*22 testmsg(2*MAXCALLS + 2 + MAXRPT)
character*15 callgrid(MAXCALLS)
character*180 line
character*4 rpt(MAXRPT)
integer ncode(63,2*MAXCALLS + 2 + MAXRPT)
real pp(2*MAXCALLS + 2 + MAXRPT)
common/mrscom/ mrs(63),mrs2(63)
common/c3com/ mcall3a
data rpt/'-01','-02','-03','-04','-05', &
'-06','-07','-08','-09','-10', &
'-11','-12','-13','-14','-15', &
'-16','-17','-18','-19','-20', &
'-21','-22','-23','-24','-25', &
'-26','-27','-28','-29','-30', &
'R-01','R-02','R-03','R-04','R-05', &
'R-06','R-07','R-08','R-09','R-10', &
'R-11','R-12','R-13','R-14','R-15', &
'R-16','R-17','R-18','R-19','R-20', &
'R-21','R-22','R-23','R-24','R-25', &
'R-26','R-27','R-28','R-29','R-30', &
'RO','RRR','73'/
save
if(mcall3a.eq.0) go to 30
call timer('deep65a ',0)
mcall3a=0
rewind 23
k=0
icall=0
do n=1,MAXCALLS
if(n.eq.1) then
callsign=hiscall
do i=4,12
if(ichar(callsign(i:i)).eq.0) callsign(i:i)=' '
enddo
grid=hisgrid(1:4)
if(ichar(grid(3:3)).eq.0) grid(3:3)=' '
if(ichar(grid(4:4)).eq.0) grid(4:4)=' '
else
read(23,1002,end=20) line
1002 format (A80)
if(line(1:4).eq.'ZZZZ') go to 20
if(line(1:2).eq.'//') go to 10
i1=index(line,',')
if(i1.lt.4) go to 10
i2=index(line(i1+1:),',')
if(i2.lt.5) go to 10
i2=i2+i1
i3=index(line(i2+1:),',')
if(i3.lt.1) i3=index(line(i2+1:),' ')
i3=i2+i3
callsign=line(1:i1-1)
grid=line(i1+1:i2-1)
ceme=line(i2+1:i3-1)
if(neme.eq.1 .and. ceme.ne.'EME') go to 10
endif
icall=icall+1
j1=index(mycall,' ') - 1
if(j1.le.-1) j1=12
if(j1.lt.3) j1=6
j2=index(callsign,' ') - 1
if(j2.le.-1) j2=12
if(j2.lt.3) j2=6
j3=index(mycall,'/') ! j3>0 means compound mycall
j4=index(callsign,'/') ! j4>0 means compound hiscall
callgrid(icall)=callsign(1:j2)
mz=1
! Allow MyCall + HisCall + rpt (?)
if(n.eq.1 .and. j3.lt.1 .and. j4.lt.1 .and. &
flip.gt.0.0 .and. callsign(1:6).ne.' ') mz=MAXRPT+1
do m=1,mz
if(m.gt.1) grid=rpt(m-1)
if(j3.lt.1 .and.j4.lt.1) callgrid(icall)=callsign(1:j2)//' '//grid
message=mycall(1:j1)//' '//callgrid(icall)
k=k+1
testmsg(k)=message
call encode65(message,ncode(1,k))
! Insert CQ message
if(j4.lt.1) callgrid(icall)=callsign(1:j2)//' '//grid
message='CQ '//callgrid(icall)
k=k+1
testmsg(k)=message
call encode65(message,ncode(1,k))
enddo
10 continue
enddo
20 continue
ntot=k
call timer('deep65a ',1)
30 continue
call timer('deep65b ',0)
ref0=0.
do j=1,63
ref0=ref0 + s3(mrs(j),j)
enddo
p1=-1.e30
p2=-1.e30
do k=1,ntot
pp(k)=0.
! Test all messages if flip=+1; skip the CQ messages if flip=-1.
if(flip.gt.0.0 .or. testmsg(k)(1:3).ne.'CQ ') then
sum=0.
ref=ref0
do j=1,63
i=ncode(j,k)+1
sum=sum + s3(i,j)
if(i.eq.mrs(j)) ref=ref - s3(i,j) + s3(mrs2(j),j)
enddo
p=sum/ref
pp(k)=p
if(p.gt.p1) then
p1=p
ip1=k
endif
endif
enddo
do i=1,ntot
if(pp(i).gt.p2 .and. pp(i).ne.p1) p2=pp(i)
enddo
! ### DO NOT REMOVE ###
rewind 77
write(77,*) p1,p2
! ### Works OK without it (in both Windows and Linux) if compiled
! ### without optimization. However, in Windows this is a colossal
! ### pain because of the way F2PY wants to run the compile step.
if(mode65.eq.1) bias=max(1.12*p2,0.335)
if(mode65.eq.2) bias=max(1.08*p2,0.405)
if(mode65.ge.4) bias=max(1.04*p2,0.505)
if(p2.eq.p1 .and. p1.ne.-1.e30) stop 'Error in deep65'
qual=100.0*(p1-bias)
decoded=' '
c=' '
if(qual.gt.1.0) then
if(qual.lt.6.0) c='?'
decoded=testmsg(ip1)
else
qual=0.
endif
decoded(22:22)=c
! Make sure everything is upper case.
do i=1,22
if(decoded(i:i).ge.'a' .and. decoded(i:i).le.'z') &
decoded(i:i)=char(ichar(decoded(i:i))-32)
enddo
call timer('deep65b ',1)
return
end subroutine deep65

73
libm65/demod64a.f
View File

@ -1,73 +0,0 @@
subroutine demod64a(s3,nadd,mrsym,mrprob,
+ mr2sym,mr2prob,ntest,nlow)
C Demodulate the 64-bin spectra for each of 63 symbols in a frame.
C Parameters
C nadd number of spectra already summed
C mrsym most reliable symbol value
C mr2sym second most likely symbol value
C mrprob probability that mrsym was the transmitted value
C mr2prob probability that mr2sym was the transmitted value
implicit real*8 (a-h,o-z)
real*4 s3(64,63)
real*8 fs(64)
integer mrsym(63),mrprob(63),mr2sym(63),mr2prob(63)
common/mrscom/ mrs(63),mrs2(63)
afac=1.1 * float(nadd)**0.64
scale=255.999
C Compute average spectral value
sum=0.
do j=1,63
do i=1,64
sum=sum+s3(i,j)
enddo
enddo
ave=sum/(64.*63.)
i1=1 !Silence warning
i2=1
C Compute probabilities for most reliable symbol values
do j=1,63
s1=-1.e30
fsum=0.
do i=1,64
x=min(afac*s3(i,j)/ave,50.d0)
fs(i)=exp(x)
fsum=fsum+fs(i)
if(s3(i,j).gt.s1) then
s1=s3(i,j)
i1=i !Most reliable
endif
enddo
s2=-1.e30
do i=1,64
if(i.ne.i1 .and. s3(i,j).gt.s2) then
s2=s3(i,j)
i2=i !Second most reliable
endif
enddo
p1=fs(i1)/fsum !Normalized probabilities
p2=fs(i2)/fsum
mrsym(j)=i1-1
mr2sym(j)=i2-1
mrprob(j)=scale*p1
mr2prob(j)=scale*p2
mrs(j)=i1
mrs2(j)=i2
enddo
sum=0.
nlow=0
do j=1,63
sum=sum+mrprob(j)
if(mrprob(j).le.5) nlow=nlow+1
enddo
ntest=sum/63
return
end

169
libm65/display.f90
View File

@ -1,169 +0,0 @@
subroutine display(nkeep,ftol)
parameter (MAXLINES=400,MX=400)
integer indx(MAXLINES),indx2(MX)
character*81 line(MAXLINES),line2(MX),line3(MAXLINES)
character out*50,cfreq0*3,cqlive*52
character*6 callsign,callsign0
character*12 freqcall(100)
real freqkHz(MAXLINES)
integer utc(MAXLINES),utc2(MX),utcz
real*8 f0
rewind 26
do i=1,MAXLINES
read(26,1010,end=10) line(i)
1010 format(a80)
read(line(i),1020) f0,ndf,nh,nm
1020 format(f8.3,i5,25x,i3,i2)
utc(i)=60*nh + nm
freqkHz(i)=1000.d0*(f0-144.d0) + 0.001d0*ndf
enddo
10 nz=i-1
utcz=utc(nz)
nz=nz-1
if(nz.lt.1) go to 999
nquad=max(nkeep/4,3)
do i=1,nz
nage=utcz-utc(i)
if(nage.lt.0) nage=nage+1440
iage=nage/nquad
write(line(i)(80:81),1021) iage
1021 format(i2)
enddo
nage=utcz-utc(1)
if(nage.lt.0) nage=nage+1440
if(nage.gt.nkeep) then
do i=1,nz
nage=utcz-utc(i)
if(nage.lt.0) nage=nage+1440
if(nage.le.nkeep) go to 20
enddo
20 i0=i
nz=nz-i0+1
rewind 26
if(nz.lt.1) go to 999
do i=1,nz
j=i+i0-1
line(i)=line(j)
utc(i)=utc(j)
freqkHz(i)=freqkHz(j)
write(26,1010) line(i)
enddo
endif
call flush(26)
call indexx(nz,freqkHz,indx)
nstart=1
k3=0
k=1
m=indx(1)
if(m.lt.1 .or. m.gt.MAXLINES) then
print*,'Error in display.f90: ',nz,m
m=1
endif
line2(1)=line(m)
utc2(1)=utc(m)
do i=2,nz
j0=indx(i-1)
j=indx(i)
if(freqkHz(j)-freqkHz(j0).gt.2.0*ftol) then
if(nstart.eq.0) then
k=k+1
line2(k)=""
utc2(k)=-1
endif
kz=k
if(nstart.eq.1) then
call indexx(kz,utc2,indx2)
k3=0
do k=1,kz
k3=min(k3+1,400)
line3(k3)=line2(indx2(k))
enddo
nstart=0
else
call indexx(kz,utc2,indx2)
do k=1,kz
k3=min(k3+1,400)
line3(k3)=line2(indx2(k))
enddo
endif
k=0
endif
if(i.eq.nz) then
k=k+1
line2(k)=""
utc2(k)=-1
endif
k=k+1
line2(k)=line(j)
utc2(k)=utc(j)
j0=j
enddo
kz=k
call indexx(kz,utc2,indx2)
do k=1,kz
k3=min(k3+1,400)
line3(k3)=line2(indx2(k))
enddo
rewind 19
rewind 20
cfreq0=' '
nc=0
callsign0=' '
do k=1,k3
out=line3(k)(6:13)//line3(k)(28:31)//line3(k)(39:43)// &
line3(k)(35:38)//line3(k)(44:67)//line3(k)(77:81)
if(out(1:3).ne.' ') then
cfreq0=out(1:3)
if(iw.lt.MAXLINES-1) iw=iw+1
cqlive=line3(k)(6:13)//line3(k)(28:31)//line3(k)(39:43)// &
line3(k)(23:27)//line3(k)(35:38)//line3(k)(44:67)// &
line3(k)(80:81)
if(index(cqlive,' CQ ').gt.0 .or. index(cqlive,' QRZ ').gt.0 .or. &
index(cqlive,' QRT ').gt.0 .or. index(cqlive,' CQV ').gt.0 .or. &
index(cqlive,' CQH ').gt.0) write(19,1029) cqlive
1029 format(a52)
write(*,1030) out
1030 format('@',a50)
i1=index(out(24:),' ')
callsign=out(i1+24:)
i2=index(callsign,' ')
if(i2.gt.1) callsign(i2:)=' '
if(callsign.ne.' ' .and. callsign.ne.callsign0) then
len=i2-1
if(len.lt.0) len=6
if(len.ge.4) then !Omit short "callsigns"
nc=nc+1
freqcall(nc)=cfreq0//' '//callsign//line3(k)(80:81)
callsign0=callsign
endif
endif
if(callsign.ne.' ' .and. callsign.eq.callsign0) then
freqcall(nc)=cfreq0//' '//callsign//line3(k)(80:81)
endif
endif
enddo
flush(19)
nc=nc+1
freqcall(nc)=' '
nc=nc+1
freqcall(nc)=' '
freqcall(nc+1)=' '
freqcall(nc+2)=' '
iz=(nc+2)/3
do i=1,nc
write(*,1042) freqcall(i)
1042 format('&',a12)
enddo
999 continue
return
end subroutine display

11
libm65/dot.f
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

41
libm65/dpol.f90
View File

@ -1,41 +0,0 @@
real function dpol(mygrid,hisgrid)
! Compute spatial polartzation offset in degrees for the present
! time, between two specified grid locators.
character*6 MyGrid,HisGrid
real lat,lon,LST
character cdate*8,ctime2*10,czone*5,fnamedate*6
integer it(8)
data rad/57.2957795/
call date_and_time(cdate,ctime2,czone,it)
nyear=it(1)
month=it(2)
nday=it(3)
nh=it(5)-it(4)/60
nm=it(6)
ns=it(7)
uth=nh + nm/60.0 + ns/3600.0
call grid2deg(MyGrid,lon,lat)
call MoonDop(nyear,month,nday,uth,-lon,lat,RAMoon,DecMoon, &
LST,HA,AzMoon,ElMoon,vr,dist)
xx=sin(lat/rad)*cos(ElMoon/rad) - cos(lat/rad)* &
cos(AzMoon/rad)*sin(ElMoon/rad)
yy=cos(lat/rad)*sin(AzMoon/rad)
poloffset1=rad*atan2(yy,xx)
call grid2deg(hisGrid,lon,lat)
call MoonDop(nyear,month,nday,uth,-lon,lat,RAMoon,DecMoon, &
LST,HA,AzMoon,ElMoon,vr,dist)
xx=sin(lat/rad)*cos(ElMoon/rad) - cos(lat/rad)* &
cos(AzMoon/rad)*sin(ElMoon/rad)
yy=cos(lat/rad)*sin(AzMoon/rad)
poloffset2=rad*atan2(yy,xx)
dpol=mod(poloffset2-poloffset1+720.0,180.0)
if(dpol.gt.90.0) dpol=dpol-180.0
return
end function dpol

13
libm65/encode65.f
View File

@ -1,13 +0,0 @@
subroutine encode65(message,sent)
character message*22
integer dgen(12)
integer sent(63)
call packmsg(message,dgen)
call rs_encode(dgen,sent)
call interleave63(sent,1)
call graycode(sent,63,1)
return
end

109
libm65/extract.F
View File

@ -1,109 +0,0 @@
subroutine extract(s3,nadd,ncount,nhist,decoded,ltext)
real s3(64,63)
real tmp(4032)
character decoded*22
integer era(51),dat4(12),indx(64)
integer mrsym(63),mr2sym(63),mrprob(63),mr2prob(63)
logical first,ltext
data first/.true./,nsec1/0/
save
nfail=0
1 continue
! call timer('demod64a',0)
call demod64a(s3,nadd,mrsym,mrprob,mr2sym,mr2prob,ntest,nlow)
! call timer('demod64a',1)
if(ntest.lt.50 .or. nlow.gt.20) then
ncount=-999 !Flag bad data
go to 900
endif
call chkhist(mrsym,nhist,ipk)
if(nhist.ge.20) then
nfail=nfail+1
call pctile(s3,tmp,4032,50,base) ! ### or, use ave from demod64a
do j=1,63
s3(ipk,j)=base
enddo
if(nfail.gt.30) then
decoded=' '
ncount=-1
go to 900
endif
go to 1
endif
call graycode(mrsym,63,-1)
call interleave63(mrsym,-1)
call interleave63(mrprob,-1)
ndec=1
nemax=30 !Was 200 (30)
maxe=8
xlambda=13.0 !Was 12
if(ndec.eq.1) then
call graycode(mr2sym,63,-1)
call interleave63(mr2sym,-1)
call interleave63(mr2prob,-1)
nsec1=nsec1+1
write(22,rec=1) nsec1,xlambda,maxe,200,
+ mrsym,mrprob,mr2sym,mr2prob
call flush(22)
! call timer('kvasd ',0)
#ifdef UNIX
iret=system('./kvasd -q > dev_null')
#else
iret=system('kvasd -q > dev_null')
#endif
! call timer('kvasd ',1)
if(iret.ne.0) then
if(first) write(*,1000) iret
1000 format('Error in KV decoder, or no KV decoder present.'/
+ 'Return code:',i8,'. Will use BM algorithm.')
ndec=0
first=.false.
go to 20
endif
read(22,rec=2) nsec2,ncount,dat4
j=nsec2 !Silence compiler warning
decoded=' '
ltext=.false.
if(ncount.ge.0) then
call unpackmsg(dat4,decoded) !Unpack the user message
if(iand(dat4(10),8).ne.0) ltext=.true.
do i=2,12
if(dat4(i).ne.dat4(1)) go to 20
enddo
write(13,*) 'Bad decode?',nhist,nfail,ipk,
+ ' ',dat4,decoded
ncount=-1 !Suppress supposedly bogus decodes
decoded=' '
endif
endif
20 if(ndec.eq.0) then
call indexx(63,mrprob,indx)
do i=1,nemax
j=indx(i)
if(mrprob(j).gt.120) then
ne2=i-1
go to 2
endif
era(i)=j-1
enddo
ne2=nemax
2 decoded=' '
do nerase=0,ne2,2
call rs_decode(mrsym,era,nerase,dat4,ncount)
if(ncount.ge.0) then
call unpackmsg(dat4,decoded)
go to 900
endif
enddo
endif
900 return
end

76
libm65/fchisq.f
View File

@ -1,76 +0,0 @@
real function fchisq(cx,cy,npts,fsample,nflip,a,ccfmax,dtmax)
parameter (NMAX=60*96000) !Samples per 60 s
complex cx(npts),cy(npts)
real a(5)
complex w,wstep,za,zb,z
real ss(2600)
complex csx(0:NMAX/64),csy(0:NMAX/64)
data twopi/6.283185307/a1,a2,a3/99.,99.,99./
save
call timer('fchisq ',0)
baud=11025.0/4096.0
if(a(1).ne.a1 .or. a(2).ne.a2 .or. a(3).ne.a3) then
a1=a(1)
a2=a(2)
a3=a(3)
C Mix and integrate the complex X and Y signals
csx(0)=0.
csy(0)=0.
w=1.0
x0=0.5*(npts+1)
s=2.0/npts
do i=1,npts
x=s*(i-x0)
if(mod(i,100).eq.1) then
p2=1.5*x*x - 0.5
! p3=2.5*(x**3) - 1.5*x
! p4=4.375*(x**4) - 3.75*(x**2) + 0.375
dphi=(a(1) + x*a(2) + p2*a(3)) * (twopi/fsample)
wstep=cmplx(cos(dphi),sin(dphi))
endif
w=w*wstep
csx(i)=csx(i-1) + w*cx(i)
csy(i)=csy(i-1) + w*cy(i)
enddo
endif
C Compute 1/2-symbol powers at 1/16-symbol steps.
fac=1.e-4
pol=a(4)/57.2957795
aa=cos(pol)
bb=sin(pol)
nsps=nint(fsample/baud) !Samples per symbol
nsph=nsps/2 !Samples per half-symbol
ndiv=16 !Output ss() steps per symbol
nout=ndiv*npts/nsps
dtstep=1.0/(ndiv*baud) !Time per output step
do i=1,nout
j=i*nsps/ndiv
k=j-nsph
ss(i)=0.
if(k.ge.1) then
za=csx(j)-csx(k)
zb=csy(j)-csy(k)
z=aa*za + bb*zb
ss(i)=fac*(real(z)**2 + aimag(z)**2)
endif
enddo
ccfmax=0.
call timer('ccf2 ',0)
call ccf2(ss,nout,nflip,ccf,lagpk)
call timer('ccf2 ',1)
if(ccf.gt.ccfmax) then
ccfmax=ccf
dtmax=lagpk*dtstep
endif
fchisq=-ccfmax
call timer('fchisq ',1)
return
end

44
libm65/fil6521.f
View File

@ -1,44 +0,0 @@
subroutine fil6521(c1,n1,c2,n2)
C FIR lowpass filter designed using ScopeFIR
C Pass #1 Pass #2
C-----------------------------------------------
C fsample (Hz) 1378.125 Input sample rate
C Ntaps 21 Number of filter taps
C fc (Hz) 40 Cutoff frequency
C fstop (Hz) 172.266 Lower limit of stopband
C Ripple (dB) 0.1 Ripple in passband
C Stop Atten (dB) 38 Stopband attenuation
C fout (Hz) 344.531 Output sample rate
parameter (NTAPS=21)
parameter (NH=NTAPS/2)
parameter (NDOWN=4) !Downsample ratio = 1/4
complex c1(n1)
complex c2(n1/NDOWN)
C Filter coefficients:
real a(-NH:NH)
data a/
+ -0.011958606980,-0.013888627387,-0.015601306443,-0.010602249570,
+ 0.003804023436, 0.028320058273, 0.060903935217, 0.096841904411,
+ 0.129639871228, 0.152644580853, 0.160917511283, 0.152644580853,
+ 0.129639871228, 0.096841904411, 0.060903935217, 0.028320058273,
+ 0.003804023436,-0.010602249570,-0.015601306443,-0.013888627387,
+ -0.011958606980/
n2=(n1-NTAPS+NDOWN)/NDOWN
k0=NH-NDOWN+1
C Loop over all output samples
do i=1,n2
c2(i)=0.
k=k0 + NDOWN*i
do j=-NH,NH
c2(i)=c2(i) + c1(j+k)*a(j)
enddo
enddo
return
end

134
libm65/filbig.f
View File

@ -1,134 +0,0 @@
subroutine filbig(dd,nmax,f0,newdat,nfsample,xpol,c4a,c4b,n4)
C Filter and downsample complex data stored in array dd(4,nmax).
C Output is downsampled from 96000 Hz to 1375.125 Hz.
parameter (MAXFFT1=5376000,MAXFFT2=77175)
real*4 dd(4,nmax) !Input data
complex ca(MAXFFT1),cb(MAXFFT1) !FFTs of input
complex c4a(MAXFFT2),c4b(MAXFFT2) !Output data
real*8 df
real halfpulse(8) !Impulse response of filter (one sided)
complex cfilt(MAXFFT2) !Filter (complex; imag = 0)
real rfilt(MAXFFT2) !Filter (real)
integer*8 plan1,plan2,plan3,plan4,plan5
logical first,xpol
include 'fftw3.f'
equivalence (rfilt,cfilt)
data first/.true./,npatience/1/
data halfpulse/114.97547150,36.57879257,-20.93789101,
+ 5.89886379,1.59355187,-2.49138308,0.60910773,-0.04248129/
save
nfft1=MAXFFT1
nfft2=MAXFFT2
if(nfsample.eq.95238) then
nfft1=5120000
nfft2=74088
endif
if(nmax.lt.0) go to 900
if(first) then
nflags=FFTW_ESTIMATE
if(npatience.eq.1) nflags=FFTW_ESTIMATE_PATIENT
if(npatience.eq.2) nflags=FFTW_MEASURE
if(npatience.eq.3) nflags=FFTW_PATIENT
if(npatience.eq.4) nflags=FFTW_EXHAUSTIVE
C Plan the FFTs just once
call timer('FFTplans ',0)
call sfftw_plan_dft_1d(plan1,nfft1,ca,ca,
+ FFTW_BACKWARD,nflags)
call sfftw_plan_dft_1d(plan2,nfft1,cb,cb,
+ FFTW_BACKWARD,nflags)
call sfftw_plan_dft_1d(plan3,nfft2,c4a,c4a,
+ FFTW_FORWARD,nflags)
call sfftw_plan_dft_1d(plan4,nfft2,c4b,c4b,
+ FFTW_FORWARD,nflags)
call sfftw_plan_dft_1d(plan5,nfft2,cfilt,cfilt,
+ FFTW_BACKWARD,nflags)
call timer('FFTplans ',1)
C Convert impulse response to filter function
do i=1,nfft2
cfilt(i)=0.
enddo
fac=0.00625/nfft1
cfilt(1)=fac*halfpulse(1)
do i=2,8
cfilt(i)=fac*halfpulse(i)
cfilt(nfft2+2-i)=fac*halfpulse(i)
enddo
call timer('FFTfilt ',0)
call sfftw_execute(plan5)
call timer('FFTfilt ',1)
base=cfilt(nfft2/2+1)
do i=1,nfft2
rfilt(i)=real(cfilt(i))-base
enddo
df=96000.d0/nfft1
if(nfsample.eq.95238) df=95238.1d0/nfft1
first=.false.
endif
C When new data comes along, we need to compute a new "big FFT"
C If we just have a new f0, continue with the existing ca and cb.
if(newdat.ne.0) then
nz=min(nmax,nfft1)
do i=1,nz
ca(i)=cmplx(dd(1,i),dd(2,i))
if(xpol) cb(i)=cmplx(dd(3,i),dd(4,i))
enddo
if(nmax.lt.nfft1) then
do i=nmax+1,nfft1
ca(i)=0.
if(xpol) cb(i)=0.
enddo
endif
call timer('FFTbig ',0)
call sfftw_execute(plan1)
if(xpol) call sfftw_execute(plan2)
call timer('FFTbig ',1)
newdat=0
endif
C NB: f0 is the frequency at which we want our filter centered.
C i0 is the bin number in ca and cb closest to f0.
i0=nint(f0/df) + 1
nh=nfft2/2
do i=1,nh !Copy data into c4a and c4b,
j=i0+i-1 !and apply the filter function
if(j.ge.1 .and. j.le.nfft1) then
c4a(i)=rfilt(i)*ca(j)
if(xpol) c4b(i)=rfilt(i)*cb(j)
else
c4a(i)=0.
if(xpol) c4b(i)=0.
endif
enddo
do i=nh+1,nfft2
j=i0+i-1-nfft2
if(j.lt.1) j=j+nfft1 !nfft1 was nfft2
c4a(i)=rfilt(i)*ca(j)
if(xpol) c4b(i)=rfilt(i)*cb(j)
enddo
C Do the short reverse transform, to go back to time domain.
call timer('FFTsmall',0)
call sfftw_execute(plan3)
if(xpol) call sfftw_execute(plan4)
call timer('FFTsmall',1)
n4=min(nmax/64,nfft2)
go to 999
900 call sfftw_destroy_plan(plan1)
call sfftw_destroy_plan(plan2)
call sfftw_destroy_plan(plan3)
call sfftw_destroy_plan(plan4)
call sfftw_destroy_plan(plan5)
999 return
end

50
libm65/foldmsk.f90
View File

@ -1,50 +0,0 @@
subroutine foldmsk(s2,msglen,nchar,mycall,msg,msg29)
! Fold the 2-d "goodness of fit" array s2 modulo message length,
! then decode the folded message.
real s2(0:63,400)
real fs2(0:63,29)
integer nfs2(29)
character mycall*12
character msg*400,msg29*29
character cc*64
! 1 2 3 4 5 6
! 0123456789012345678901234567890123456789012345678901234567890123
data cc/'0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ./?- _ @'/
fs2=0.
nfs2=0
do j=1,nchar !Fold s2 into fs2, modulo msglen
jj=mod(j-1,msglen)+1
nfs2(jj)=nfs2(jj)+1
do i=0,40
fs2(i,jj)=fs2(i,jj) + s2(i,j)
enddo
enddo
msg=' '
do j=1,msglen
smax=0.
do k=0,40
if(fs2(k,j).gt.smax) then
smax=fs2(k,j)
kpk=k
endif
enddo
if(kpk.eq.40) kpk=57
msg(j:j)=cc(kpk+1:kpk+1)
if(kpk.eq.57) msg(j:j)=' '
enddo
msg29=msg(1:msglen)
call alignmsg(' ',2,msg29,msglen,idone)
if(idone.eq.0) call alignmsg('CQ', 3,msg29,msglen,idone)
if(idone.eq.0) call alignmsg('QRZ', 3,msg29,msglen,idone)
if(idone.eq.0) call alignmsg(mycall,4,msg29,msglen,idone)
if(idone.eq.0) call alignmsg(' ', 1,msg29,msglen,idone)
msg29=adjustl(msg29)
return
end subroutine foldmsk

93
libm65/gen65.f90
View File

@ -1,93 +0,0 @@
subroutine gen65(message,mode65,samfac,nsendingsh,msgsent,iwave,nwave)
! Encodes a JT65 message into a wavefile.
! Executes in 17 ms on opti-745.
parameter (NMAX=60*11025) !Max length of wave file
character*22 message !Message to be generated
character*22 msgsent !Message as it will be received
character*3 cok !' ' or 'OOO'
real*8 dt,phi,f,f0,dfgen,dphi,twopi,samfac
integer*2 iwave(NMAX) !Generated wave file
integer dgen(12)
integer sent(63)
logical first
integer nprc(126)
real pr(126)
data nprc/1,0,0,1,1,0,0,0,1,1,1,1,1,1,0,1,0,1,0,0, &
0,1,0,1,1,0,0,1,0,0,0,1,1,1,0,0,1,1,1,1, &
0,1,1,0,1,1,1,1,0,0,0,1,1,0,1,0,1,0,1,1, &
0,0,1,1,0,1,0,1,0,1,0,0,1,0,0,0,0,0,0,1, &
1,0,0,0,0,0,0,0,1,1,0,1,0,0,1,0,1,1,0,1, &
0,1,0,1,0,0,1,1,0,0,1,0,0,1,0,0,0,0,1,1, &
1,1,1,1,1,1/
data twopi/6.283185307179586476d0/,first/.true./
save
if(first) then
do i=1,126
pr(i)=2*nprc(i)-1
enddo
first=.false.
endif
call chkmsg(message,cok,nspecial,flip)
if(nspecial.eq.0) then
call packmsg(message,dgen) !Pack message into 72 bits
nsendingsh=0
if(iand(dgen(10),8).ne.0) nsendingsh=-1 !Plain text flag
call rs_encode(dgen,sent)
call interleave63(sent,1) !Apply interleaving
call graycode(sent,63,1) !Apply Gray code
nsym=126 !Symbols per transmission
nsps=4096
else
nsym=32
nsps=16384
nsendingsh=1 !Flag for shorthand message
endif
if(mode65.eq.0) go to 900
! Set up necessary constants
dt=1.d0/(samfac*11025.d0)
f0=118*11025.d0/1024
dfgen=mode65*11025.d0/4096.d0
phi=0.d0
i=0
k=0
do j=1,nsym
f=f0
if(nspecial.ne.0 .and. mod(j,2).eq.0) f=f0+10*nspecial*dfgen
if(nspecial.eq.0 .and. flip*pr(j).lt.0.0) then
k=k+1
f=f0+(sent(k)+2)*dfgen
endif
dphi=twopi*dt*f
do ii=1,nsps
phi=phi+dphi
if(phi.gt.twopi) phi=phi-twopi
xphi=phi
i=i+1
iwave(i)=32767.0*sin(xphi)
enddo
enddo
iwave(nsym*nsps+1:)=0
nwave=nsym*nsps + 5512
call unpackmsg(dgen,msgsent)
if(flip.lt.0.0) then
do i=22,1,-1
if(msgsent(i:i).ne.' ') goto 10
enddo
10 msgsent=msgsent(1:i)//' OOO'
endif
if(nsendingsh.eq.1) then
if(nspecial.eq.2) msgsent='RO'
if(nspecial.eq.3) msgsent='RRR'
if(nspecial.eq.4) msgsent='73'
endif
900 return
end subroutine gen65

158
libm65/genjtms3.f90
View File

@ -1,158 +0,0 @@
subroutine genjtms3(msg28,iwave,nwave)
!subroutine genjtms3(msg28,iwave,cwave,isrch,nwave)
! Generate a JTMS3 wavefile.
parameter (NMAX=30*48000) !Max length of wave file
integer*2 iwave(NMAX) !Generated wave file
complex cwave(NMAX) !Alternative for searchms
character*28 msg28 !User message
character*29 msg
character cc*64
integer sentsym(203) !Transmitted symbols (0/1)
real sentsam(4872) !Transmitted waveform
real*8 dt,phi,f0,dphi,pi,twopi,samfac
real p(0:420)
real carrier(4872)
real dat(4872),bb(4872),wave(4872)
complex cdat(0:2436)
logical first
integer np(9)
data np/5,7,9,11,13,17,19,23,29/ !Permissible message lengths
! 1 2 3 4 5 6
! 0123456789012345678901234567890123456789012345678901234567890123
data cc/'0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ./?- _ @'/
data samfac/1.d0/,first/.true./
equivalence (dat,cdat)
save
sinc(x)=sin(pi*x)/(pi*x)
if(first) then
pi=4.d0*atan(1.d0)
twopi=2.d0*pi
x=0.
dx=1.0/24.0
width=3.0
do i=1,420 !Generate the BPSK pulse shape
x=x+dx
fac=0.0
if(x/width.lt.0.5*pi) then
fac=(cos(x/width))**2
ipz=i
endif
p(i)=fac*sinc(x)
enddo
p(0)=1.0
f0=10000.d0/7.d0
dt=1.d0/48000.d0
dphi=twopi*f0*dt
phi=0.d0
do i=1,4872 !Generate the carrier
phi=phi+dphi
if(phi.gt.twopi)phi=phi-twopi
xphi=phi
carrier(i)=sin(xphi)
enddo
first=.false.
endif
msg=msg28//' ' !Extend to 29 characters
do i=28,1,-1 !Find user's message length
if(msg(i:i).ne.' ') go to 1
enddo
1 iz=i+1 !Add one for space at EOM
msglen=iz
if(isrch.ne.0) go to 3
do i=1,9
if(np(i).ge.iz) go to 2
enddo
i=8
2 msglen=np(i)
! Convert message to a bit sequence, 7 bits per character (6 + even parity)
3 sentsym=0
k=0
do j=1,msglen
if(msg(j:j).eq.' ') then
i=58
go to 5
else
do i=1,64
if(msg(j:j).eq.cc(i:i)) go to 5
enddo
endif
5 m=0
do n=5,0,-1 !Each character gets 6 bits
k=k+1
sentsym(k)=iand(1,ishft(i-1,-n))
m=m+sentsym(k)
enddo
k=k+1
sentsym(k)=iand(m,1) !Insert bit for even parity
enddo
nsym=7*msglen !# symbols in message
nsam=24*nsym !# samples in message
bb(1:nsam)=0.
do j=1,nsym
fac=1.0
if(sentsym(j).eq.0) fac=-1.0
k0=24*j - 23
do i=0,ipz
k=k0+i
if(k.gt.nsam) k=k-nsam
bb(k)=bb(k) + fac*p(i)
if(i.gt.0) then
k=k0-i
if(k.lt.1) k=k+nsam
bb(k)=bb(k) + fac*p(i)
endif
enddo
enddo
sq=0.
wmax=0.
do i=1,nsam
wave(i)=carrier(i)*bb(i)
sq=sq + wave(i)**2
wmax=max(wmax,abs(wave(i)))
! write(15,3002) i*dt,bb(i),wave(i)
!3002 format(f12.6,2f12.3)
enddo
rms=sqrt(sq/nsam)
! print*,rms,wmax,wmax/rms
fac=32767.0/wmax
iwave(1:nsam)=fac*wave(1:nsam)
! nwave=nsam
nblk=30*48000/nsam
do n=2,nblk
i0=(n-1)*nsam
iwave(i0+1:i0+nsam)=iwave(1:nsam)
enddo
nwave=i0+nsam
! Compute the spectrum
! nfft=nsam
! df=48000.0/nfft
! ib=4000.0/df
! fac=10.0/nfft
! dat(1:nfft)=fac*bb(1:nfft)
! call four2a(dat,nfft,1,-1,0)
! do i=0,ib
! sq=real(cdat(i))**2 + aimag(cdat(i))**2
! write(14,3010) i*df,sq,10.0*log10(sq)
!3010 format(3f12.3)
! enddo
! if(isrch.eq.0) iwave(k+1:)=0
! nwave=k
return
end subroutine genjtms3

17
libm65/geocentric.f
View File

@ -1,17 +0,0 @@
subroutine geocentric(alat,elev,hlt,erad)
implicit real*8 (a-h,o-z)
C 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

18
libm65/getdphi.f90
View File

@ -1,18 +0,0 @@
subroutine getdphi(qphi)
real qphi(12)
s=0.
c=0.
do i=1,12
th=i*30/57.2957795
s=s+qphi(i)*sin(th)
c=c+qphi(i)*cos(th)
enddo
dphi=57.2957795*atan2(s,c)
write(*,1010) nint(dphi)
1010 format('!Best-fit Dphi =',i4,' deg')
return
end

23
libm65/hipass.f90
View File

@ -1,23 +0,0 @@
subroutine hipass(y,npts,nwidth)
! Hipass filter for time-domain data. Removes an RC-type running
! mean (time constant nwidth) from array y(1:npts).
real y(npts)
c1=1.0/nwidth
c2=1.0-c1
s=0.
do i=1,nwidth !Get initial average
s=s+y(i)
enddo
ave=c1*s
do i=1,npts !Do the filtering
y0=y(i)
y(i)=y0-ave !Remove the mean
ave=c1*y0 + c2*ave !Update the mean
enddo
return
end subroutine hipass

25
libm65/interleave63.f
View File

@ -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

30
libm65/iqcal.f90
View File

@ -1,30 +0,0 @@
subroutine iqcal(nn,c,nfft,gain,phase,zsum,ipk,reject)
complex c(0:nfft-1)
complex z,zsum,zave
if(nn.eq.0) then
zsum=0.
endif
nn=nn+1
smax=0.
ipk=1
do i=1,nfft-1 !Find strongest signal
s=real(c(i))**2 + aimag(c(i))**2
if(s.gt.smax) then
smax=s
ipk=i
endif
enddo
pimage=real(c(nfft-ipk))**2 + aimag(c(nfft-ipk))**2
p=smax + pimage
z=c(ipk)*c(nfft-ipk)/p !Synchronous detection of image
zsum=zsum+z
zave=zsum/nn
tmp=sqrt(1.0 - (2.0*real(zave))**2)
phase=asin(2.0*aimag(zave)/tmp) !Estimate phase
gain=tmp/(1.0-2.0*real(zave)) !Estimate gain
reject=10.0*log10(pimage/smax)
return
end subroutine iqcal

29
libm65/iqfix.f90
View File

@ -1,29 +0,0 @@
subroutine iqfix(c,nfft,gain,phase)
complex c(0:nfft-1)
complex z,h,u,v
real*8 sq1,sq2
nh=nfft/2
h=gain*cmplx(cos(phase),sin(phase))
do i=1,nh-1
u=c(i)
v=c(nfft-i)
x=real(u) + real(v) - (aimag(u) + aimag(v))*aimag(h) + &
(real(u) - real(v))*real(h)
y=aimag(u) - aimag(v) + (aimag(u) + aimag(v))*real(h) + &
(real(u) - real(v))*aimag(h)
c(i)=0.5*cmplx(x,y)
z=u
u=v
v=z
x=real(u) + real(v) - (aimag(u) + aimag(v))*aimag(h) + &
(real(u) - real(v))*real(h)
y=aimag(u) - aimag(v) + (aimag(u) + aimag(v))*real(h) + &
(real(u) - real(v))*aimag(h)
c(nfft-i)=0.5*cmplx(x,y)
enddo
return
end subroutine iqfix

72
libm65/jtmsk.f90
View File

@ -1,72 +0,0 @@
subroutine jtmsk(dat,npts,cfile6,tpk,mswidth,ndb,nrpt,Nfreeze, &
ntol,MouseDF,pick,mycall)
! Decode a JTMSK ping
parameter (NZ=2*48000)
real dat(npts) !Raw data
complex cdat(NZ) !Analytic form of signal
character*6 cfile6 !FileID
logical pick
character*12 mycall
real s(NZ) !Power spectrum
real s2(0:63,400)
real r(60000)
complex cw(168,0:63) !Complex waveforms for all codewords
complex cwb(168) !Complex waveform for <space>
logical first
character msg*400,msg29*29
data first/.true./
save first,cw,cwb
save cdat !Fix its address, for four2
if(first) call setupmsk(cw,cwb) !Calculate waveforms for codewords
first=.false.
nsps=24 !Samples per symbol
f0=1000.0 !Nominal frequency for bit=0
n=log(float(npts))/log(2.0) + 1.0
nfft1=2**n !FFT length
call analytic(dat,npts,nfft1,s,cdat) !Convert to analytic signal
call mskdf(cdat,npts,tpk,nfft1,f0,nfreeze,mousedf,ntol, &
dfx,snrsq2) !Get DF
sq2lim=7.0
if(pick) sq2lim=5.0
if(snrsq2.lt.sq2lim) go to 900 !Reject non-JTMS signals
call tweak1(cdat,npts,-dfx,cdat) !Mix to standard frequency
call syncmsk(cdat,npts,cwb,r,i1) !Get character sync
call lenmsk(r,npts,msglen) !Find message length
s2=0.
nchar=(npts-168+1-i1)/168
if(nchar.gt.400) nchar=400
call decodemsk(cdat,npts,cw,i1,nchar,s2,msg) !Decode the message
ia=1
if(nchar.ge.40) ia=min(nchar/3,nchar-28)
ib=min(ia+28,nchar) !Can better limits ia, ib be found?
msg29=adjustl(msg(ia:ib))
msg=adjustl(msg)
ib=min(nchar,45)
ndf=nint(dfx)
nchk=max(20,nint(1.5*msglen))
if(msglen.eq.0 .or. nchar.lt.nchk .or. pick) then
write(*,1110) cfile6,tpk,mswidth,ndb,nrpt,ndf,msg(1:45)
1110 format(a6,f5.1,i5,i3,1x,i2.2,i5,5x,a45)
endif
if(msglen.gt.0 .and. nchar.ge.nchk) then
call foldmsk(s2,msglen,nchar,mycall,msg,msg29) !Decode folded message
write(*,1120) cfile6,tpk,mswidth,ndb,nrpt,ndf,msg29
1120 format(a6,f5.1,i5,i3,1x,i2.2,i5,5x,a29,11x,'*')
endif
900 continue
return
end subroutine jtmsk

56
libm65/lenmsk.f90
View File

@ -1,56 +0,0 @@
subroutine lenmsk(r,npts,msglen)
! Determine length of the user message in a JTMS ping.
real r(60000)
real acf(4872)
integer np(9)
data np/5,7,9,11,13,17,19,23,29/ !Permissible message lengths
save acf !Why necessary? (But don't remove!)
msglen=0 !Use ACF to find msg length
if(npts.ge.8*168) then
r=r-sum(r(1:npts))/npts
acfmax=0.
acf0=dot_product(r(1:npts),r(1:npts))
kz=min(nint(0.75*npts),29*168)
do k=8,kz
fac=float(npts)/(npts-k)
acf(k)=fac*dot_product(r(1:npts),r(1+k:npts+k))/acf0
enddo
call hipass(acf(8),kz-7,50)
do k=8,kz !Find acfmax, kpk
if(acf(k).gt.acfmax) then
acfmax=acf(k)
kpk=k
endif
enddo
sumsq=0.
n=0
do k=8,kz !Find rms, skipping around kpk
if(abs(k-kpk).gt.10) then
sumsq=sumsq+acf(k)**2
n=n+1
endif
enddo
rms=sqrt(sumsq/n)
acf=acf/rms !Normalize the acf
amax2=0.
acflim=3.5
do i=1,9
k=168*np(i) !Check only the permitted lengths
if(k.gt.kz) go to 10
if(acf(k).gt.acflim .and. acf(k).gt.amax2) then
amax2=acf(k) !Save best value >3.5 sigma
msglen=np(i) !Save message length
kpk2=k
endif
enddo
10 continue
endif
return
end subroutine lenmsk

130
libm65/m65.f90
View File

@ -1,130 +0,0 @@
program m65
! Decoder for map65. Can run stand-alone, reading data from *.tf2 files;
! or as the back end of map65, with data placed in a shared memory region.
parameter (NSMAX=60*96000)
parameter (NFFT=32768)
integer*2 i2(4,87)
real*8 hsym
real*4 ssz5a(NFFT)
logical*1 lstrong(0:1023)
common/tracer/limtrace,lu
real*8 fc0,fcenter
character*80 arg,infile
character mycall*12,hiscall*12,mygrid*6,hisgrid*6,datetime*20
common/datcom/dd(4,5760000),ss(4,322,NFFT),savg(4,NFFT),fc0,nutc0,junk(34)
common/npar/fcenter,nutc,idphi,mousedf,mousefqso,nagain, &
ndepth,ndiskdat,neme,newdat,nfa,nfb,nfcal,nfshift, &
mcall3,nkeep,ntol,nxant,nrxlog,nfsample,nxpol,mode65, &
mycall,mygrid,hiscall,hisgrid,datetime
nargs=iargc()
if(nargs.lt.1) then
print*,'Usage: m65 [95238] file1 [file2 ...]'
print*,' Reads data from *.tf2 files.'
print*,''
print*,' m65 -s'
print*,' Gets data from shared memory region.'
go to 999
endif
call getarg(1,arg)
if(arg(1:2).eq.'-s') then
call m65a
call ftnquit
go to 999
endif
nfsample=96000
nxpol=1
mode65=2
ifile1=1
if(arg.eq.'95238') then
nfsample=95238
call getarg(2,arg)
ifile1=2
endif
limtrace=0
lu=12
nfa=100
nfb=162
nfshift=6
ndepth=2
nfcal=344
idphi=-50
ntol=500
nkeep=10
call ftninit('.')
do ifile=ifile1,nargs
call getarg(ifile,infile)
open(10,file=infile,access='stream',status='old',err=998)
i1=index(infile,'.tf2')
read(infile(i1-4:i1-1),*,err=1) nutc0
go to 2
1 nutc0=0
2 hsym=2048.d0*96000.d0/11025.d0 !Samples per half symbol
nhsym0=-999
k=0
fcenter=144.125d0
mousedf=0
mousefqso=125
newdat=1
mycall='K1JT'
if(ifile.eq.ifile1) call timer('m65 ',0)
do irec=1,9999999
call timer('read_tf2',0)
read(10) i2
call timer('read_tf2',1)
call timer('float ',0)
do i=1,87
k=k+1
dd(1,k)=i2(1,i)
dd(2,k)=i2(2,i)
dd(3,k)=i2(3,i)
dd(4,k)=i2(4,i)
enddo
call timer('float ',1)
nhsym=(k-2048)/hsym
if(nhsym.ge.1 .and. nhsym.ne.nhsym0) then
ndiskdat=1
nb=0
! Emit signal readyForFFT
call timer('symspec ',0)
fgreen=-13.0
iqadjust=1
iqapply=1
nbslider=100
gainx=0.9962
gainy=1.0265
phasex=0.01426
phasey=-0.01195
call symspec(k,nxpol,ndiskdat,nb,nbslider,idphi,nfsample,fgreen, &
iqadjust,iqapply,gainx,gainy,phasex,phasey,rejectx,rejecty, &
pxdb,pydb,ssz5a,nkhz,ihsym,nzap,slimit,lstrong)
call timer('symspec ',1)
nhsym0=nhsym
if(ihsym.ge.278) go to 10
endif
enddo
10 continue
if(iqadjust.ne.0) write(*,3002) rejectx,rejecty
3002 format('Image rejection:',2f7.1,' dB')
nutc=nutc0
nstandalone=1
call decode0(dd,ss,savg,nstandalone,nfsample)
enddo
call timer('m65 ',1)
call timer('m65 ',101)
call ftnquit
go to 999
998 print*,'Cannot open file:'
print*,infile
999 end program m65

97
libm65/m65a.F90
View File

@ -1,97 +0,0 @@
subroutine m65a
! NB: this interface block is required by g95, but must be omitted
! for gfortran. (????)
#ifndef UNIX
interface
function address_m65()
end function address_m65
end interface
#endif
integer*1 attach_m65,lock_m65,unlock_m65
integer size_m65
integer*1, pointer :: address_m65,p_m65
character*80 cwd
logical fileExists
common/tracer/limtrace,lu
call getcwd(cwd)
call ftninit(trim(cwd))
limtrace=0
lu=12
i1=attach_m65()
10 inquire(file=trim(cwd)//'/.lock',exist=fileExists)
if(fileExists) then
call sleep_msec(100)
go to 10
endif
inquire(file=trim(cwd)//'/.quit',exist=fileExists)
if(fileExists) then
call ftnquit
i=detach_m65()
go to 999
endif
nbytes=size_m65()
if(nbytes.le.0) then
print*,'m65a: Shared memory mem_m65 does not exist.'
print*,'Program m65a should be started automatically from within map65.'
go to 999
endif
p_m65=>address_m65()
call m65b(p_m65,nbytes)
write(*,1010)
1010 format('<m65aFinished>')
flush(6)
100 inquire(file=trim(cwd)//'/.lock',exist=fileExists)
if(fileExists) go to 10
call sleep_msec(100)
go to 100
999 return
end subroutine m65a
subroutine m65b(m65com,nbytes)
integer*1 m65com(0:nbytes-1)
kss=4*4*60*96000
ksavg=kss+4*4*322*32768
kfcenter=ksavg+4*4*32768
call m65c(m65com(0),m65com(kss),m65com(ksavg),m65com(kfcenter))
return
end subroutine m65b
subroutine m65c(dd,ss,savg,nparams0)
integer*1 detach_m65
real*4 dd(4,5760000),ss(4,322,32768),savg(4,32768)
real*8 fcenter
integer nparams0(37),nparams(37)
character*12 mycall,hiscall
character*6 mygrid,hisgrid
character*20 datetime
common/npar/fcenter,nutc,idphi,mousedf,mousefqso,nagain, &
ndepth,ndiskdat,neme,newdat,nfa,nfb,nfcal,nfshift, &
mcall3,nkeep,ntol,nxant,nrxlog,nfsample,nxpol,mode65, &
mycall,mygrid,hiscall,hisgrid,datetime
equivalence (nparams,fcenter)
nparams=nparams0 !Copy parameters into common/npar/
npatience=1
if(iand(nrxlog,1).ne.0) then
write(21,1000) datetime(:17)
1000 format(/'UTC Date: 'a17/78('-'))
flush(21)
endif
if(iand(nrxlog,2).ne.0) rewind 21
if(iand(nrxlog,4).ne.0) rewind 26
nstandalone=0
if(sum(nparams).ne.0) call decode0(dd,ss,savg,nstandalone)
return
end subroutine m65c

29
libm65/makepings.f90
View File

@ -1,29 +0,0 @@
subroutine makepings(iwave,nwave)
parameter (NFSAMPLE=48000)
integer*2 iwave(nwave)
real*8 t
iping0=-999
dt=1.0/NFSAMPLE
do i=1,nwave
iping=i/(3*NFSAMPLE)
if(iping.ne.iping0) then
ip=mod(iping,3)
w=0.015*(4-ip)
ig=(iping-1)/3
amp=sqrt((3.0-ig)/3.0)
t0=dt*(iping+0.5)*(3*NFSAMPLE)
iping0=iping
endif
t=(i*dt-t0)/w
if(t.lt.0.d0 .and. t.lt.10.0) then
fac=0.
else
fac=2.718*t*dexp(-t)
endif
iwave(i)=nint(fac*amp*iwave(i))
enddo
return
end subroutine makepings

438
libm65/map65a.f90
View File

@ -1,438 +0,0 @@
subroutine map65a(dd,ss,savg,newdat,nutc,fcenter,ntol,idphi,nfa,nfb, &
mousedf,mousefqso,nagain,ndecdone,ndiskdat,nfshift,ndphi, &
nfcal,nkeep,mcall3b,nsave,nxant,rmsdd,mycall,mygrid, &
neme,ndepth,hiscall,hisgrid,nhsym,nfsample,nxpol,mode65)
! Processes timf2 data from Linrad to find and decode JT65 signals.
parameter (MAXMSG=1000) !Size of decoded message list
parameter (NSMAX=60*96000)
parameter (NFFT=32768)
real dd(4,NSMAX)
real*4 ss(4,322,NFFT),savg(4,NFFT)
real tavg(-50:50) !Temp for finding local base level
real base(4) !Local basel level at 4 pol'ns
real tmp (200) !Temp storage for pctile sorting
real sig(MAXMSG,30) !Parameters of detected signals
real a(5)
real*8 fcenter
character*22 msg(MAXMSG)
character*3 shmsg0(4)
character mycall*12,hiscall*12,mygrid*6,hisgrid*6,grid*6,cp*1
integer indx(MAXMSG),nsiz(MAXMSG)
logical done(MAXMSG)
logical xpol
character decoded*22,blank*22
real short(3,NFFT) !SNR dt ipol for potential shorthands
real qphi(12)
common/c3com/ mcall3a
common/testcom/ifreq
data blank/' '/
data shmsg0/'ATT','RO ','RRR','73 '/
data nfile/0/,nutc0/-999/,nid/0/,ip000/1/,ip001/1/,mousefqso0/-999/
save
mcall3a=mcall3b
mousefqso0=mousefqso
xpol=(nxpol.ne.0)
if(.not.xpol) ndphi=0
!### Should use AppDir! ###
! open(23,file='release/CALL3.TXT',status='unknown')
open(23,file='CALL3.TXT',status='unknown')
if(nutc.ne.nutc0) nfile=nfile+1
nutc0=nutc
df=96000.0/NFFT !df = 96000/NFFT = 2.930 Hz
if(nfsample.eq.95238) df=95238.1/NFFT
ftol=0.010 !Frequency tolerance (kHz)
dphi=idphi/57.2957795
foffset=0.001*(1270 + nfcal) !Offset from sync tone, plus CAL
fqso=mousefqso + foffset - 0.5*(nfa+nfb) + nfshift !fqso at baseband (khz)
iloop=0
2 if(ndphi.eq.1) dphi=30*iloop/57.2957795
do nqd=1,0,-1
if(nqd.eq.1) then !Quick decode, at fQSO
fa=1000.0*(fqso+0.001*mousedf) - ntol
fb=1000.0*(fqso+0.001*mousedf) + ntol + 4*53.8330078
else !Wideband decode at all freqs
fa=-1000*0.5*(nfb-nfa) + 1000*nfshift
fb= 1000*0.5*(nfb-nfa) + 1000*nfshift
endif
ia=nint(fa/df) + 16385
ib=nint(fb/df) + 16385
ia=max(51,ia)
ib=min(32768-51,ib)
km=0
nkm=1
nz=n/8
freq0=-999.
sync10=-999.
fshort0=-999.
syncshort0=-999.
ntry=0
short=0. !Zero the whole short array
jpz=1
if(xpol) jpz=4
do i=ia,ib !Search over freq range
freq=0.001*(i-16385)*df
! Find the local base level for each polarization; update every 10 bins.
if(mod(i-ia,10).eq.0) then
do jp=1,jpz
do ii=-50,50
iii=i+ii
if(iii.ge.1 .and. iii.le.32768) then
tavg(ii)=savg(jp,iii)
else
write(13,*) ,'Error in iii:',iii,ia,ib,fa,fb
flush(13)
go to 999
endif
enddo
call pctile(tavg,tmp,101,50,base(jp))
enddo
endif
! Find max signal at this frequency
smax=0.
do jp=1,jpz
if(savg(jp,i)/base(jp).gt.smax) then
smax=savg(jp,i)/base(jp)
jpmax=jp
endif
enddo
if(smax.gt.1.1) then
! Look for JT65 sync patterns and shorthand square-wave patterns.
call timer('ccf65 ',0)
! ssmax=4.0*(rmsdd/22.5)**2
ssmax=savg(jpmax,i)
call ccf65(ss(1,1,i),nhsym,ssmax,sync1,ipol,jpz,dt,flipk, &
syncshort,snr2,ipol2,dt2)
call timer('ccf65 ',1)
! ########################### Search for Shorthand Messages #################
! Is there a shorthand tone above threshold?
thresh0=1.0
! Use lower thresh0 at fQSO
if(nqd.eq.1 .and. ntol.le.100) thresh0=0.
if(syncshort.gt.thresh0) then
! ### Do shorthand AFC here (or maybe after finding a pair?) ###
short(1,i)=syncshort
short(2,i)=dt2
short(3,i)=ipol2
! Check to see if lower tone of shorthand pair was found.
do j=2,4
i0=i-nint(j*mode65*10.0*(11025.0/4096.0)/df)
! Should this be i0 +/- 1, or just i0?
! Should we also insist that difference in DT be either 1.5 or -1.5 s?
if(short(1,i0).gt.thresh0) then
fshort=0.001*(i0-16385)*df
noffset=0
if(nqd.eq.1) noffset=nint(1000.0*(fshort-fqso)-mousedf)
if(abs(noffset).le.ntol) then
! Keep only the best candidate within ftol.
!### NB: sync2 was not defined here!
! sync2=syncshort !### try this ???
if(fshort-fshort0.le.ftol .and. syncshort.gt.syncshort0 &
.and. nkm.eq.2) km=km-1
if(fshort-fshort0.gt.ftol .or. &
syncshort.gt.syncshort0) then
if(km.lt.MAXMSG) km=km+1
sig(km,1)=nfile
sig(km,2)=nutc
sig(km,3)=fshort + 0.5*(nfa+nfb)
sig(km,4)=syncshort
sig(km,5)=dt2
sig(km,6)=45*(ipol2-1)/57.2957795
sig(km,7)=0
sig(km,8)=snr2
sig(km,9)=0
sig(km,10)=0
! sig(km,11)=rms0
sig(km,12)=savg(ipol2,i)
sig(km,13)=0
sig(km,14)=0
sig(km,15)=0
sig(km,16)=0
! sig(km,17)=0
sig(km,18)=0
msg(km)=shmsg0(j)
fshort0=fshort
syncshort0=syncshort
nkm=2
endif
endif
endif
enddo
endif
! ########################### Search for Normal Messages ###########
! Is sync1 above threshold?
thresh1=1.0
! Use lower thresh1 at fQSO
if(nqd.eq.1 .and. ntol.le.100) thresh1=0.
noffset=0
if(nqd.eq.1) noffset=nint(1000.0*(freq-fqso)-mousedf)
if(sync1.gt.thresh1 .and. abs(noffset).le.ntol) then
! Keep only the best candidate within ftol.
! (Am I deleting any good decodes by doing this?)
if(freq-freq0.le.ftol .and. sync1.gt.sync10 .and. &
nkm.eq.1) km=km-1
if(freq-freq0.gt.ftol .or. sync1.gt.sync10) then
nflip=nint(flipk)
f00=(i-1)*df !Freq of detected sync tone (0-96000 Hz)
ntry=ntry+1
if((nqd.eq.1 .and. ntry.ge.40) .or. &
(nqd.eq.0 .and. ntry.ge.400)) then
! Too many calls to decode1a!
write(*,*) '! Signal too strong, or suspect data? Decoding aborted.'
write(13,*) 'Signal too strong, or suspect data? Decoding aborted.'
call flush(13)
go to 999
endif
call timer('decode1a',0)
ifreq=i
ikHz=nint(freq+0.5*(nfa+nfb)-foffset)-nfshift
idf=nint(1000.0*(freq+0.5*(nfa+nfb)-foffset-(ikHz+nfshift)))
call decode1a(dd,newdat,f00,nflip,mode65,nfsample,xpol, &
mycall,hiscall,hisgrid,neme,ndepth,nqd,dphi, &
nutc,ikHz,idf,ipol,sync2,a,dt,pol,nkv,nhist,qual,decoded)
dt=dt+0.8 !### empirical tweak
call timer('decode1a',1)
if(km.lt.MAXMSG) km=km+1
sig(km,1)=nfile
sig(km,2)=nutc
sig(km,3)=freq + 0.5*(nfa+nfb)
sig(km,4)=sync1
sig(km,5)=dt
sig(km,6)=pol
sig(km,7)=flipk
sig(km,8)=sync2
sig(km,9)=nkv
sig(km,10)=qual
! sig(km,11)=idphi
sig(km,12)=savg(ipol,i)
sig(km,13)=a(1)
sig(km,14)=a(2)
sig(km,15)=a(3)
sig(km,16)=a(4)
! sig(km,17)=a(5)
sig(km,18)=nhist
msg(km)=decoded
freq0=freq
sync10=sync1
nkm=1
endif
endif
endif
!70 continue
enddo
if(nqd.eq.1) then
nwrite=0
do k=1,km
decoded=msg(k)
if(decoded.ne.' ') then
nutc=sig(k,2)
freq=sig(k,3)
sync1=sig(k,4)
dt=sig(k,5)
npol=nint(57.2957795*sig(k,6))
flip=sig(k,7)
sync2=sig(k,8)
nkv=sig(k,9)
nqual=sig(k,10)
! idphi=nint(sig(k,11))
if(flip.lt.0.0) then
do i=22,1,-1
if(decoded(i:i).ne.' ') go to 8
enddo
stop 'Error in message format'
8 if(i.le.18) decoded(i+2:i+4)='OOO'
endif
nkHz=nint(freq-foffset)-nfshift
mhz=fcenter ! ... +fadd ???
f0=mhz+0.001*nkHz
ndf=nint(1000.0*(freq-foffset-(nkHz+nfshift)))
nsync1=sync1
nsync2=nint(10.0*log10(sync2)) - 40 !### empirical ###
if(decoded(1:4).eq.'RO ' .or. decoded(1:4).eq.'RRR ' .or. &
decoded(1:4).eq.'73 ') nsync2=nsync2-6
nwrite=nwrite+1
if(nxant.ne.0) then
npol=npol-45
if(npol.lt.0) npol=npol+180
endif
! If Tx station's grid is in decoded message, compute optimum TxPol
i1=index(decoded,' ')
i2=index(decoded(i1+1:),' ') + i1
grid=' '
if(i2.ge.8 .and. i2.le.18) grid=decoded(i2+1:i2+4)//'mm'
ntxpol=0
cp=' '
if(xpol) then
if(grid(1:1).ge.'A' .and. grid(1:1).le.'R' .and. &
grid(2:2).ge.'A' .and. grid(2:2).le.'R' .and. &
grid(3:3).ge.'0' .and. grid(3:3).le.'9' .and. &
grid(4:4).ge.'0' .and. grid(4:4).le.'9') then
ntxpol=mod(npol-nint(2.0*dpol(mygrid,grid))+720,180)
if(nxant.eq.0) then
cp='H'
if(ntxpol.gt.45 .and. ntxpol.le.135) cp='V'
else
cp='/'
if(ntxpol.ge.90 .and. ntxpol.lt.180) cp='\\'
endif
endif
endif
if(ndphi.eq.0) then
write(*,1010) nkHz,ndf,npol,nutc,dt,nsync2, &
decoded,nkv,nqual,ntxpol,cp
1010 format('!',i3,i5,i4,i5.4,f5.1,i4,2x,a22,i4,i5,i5,1x,a1)
else
if(iloop.ge.1) qphi(iloop)=sig(k,10)
write(*,1010) nkHz,ndf,npol,nutc,dt,nsync2, &
decoded,nkv,nqual,30*iloop
write(27,1011) 30*iloop,nkHz,ndf,npol,nutc, &
dt,sync2,nkv,nqual,decoded
1011 format(i3,i4,i5,i4,i5.4,f5.1,f7.1,i3,i5,2x,a22)
endif
endif
enddo
if(nwrite.eq.0) then
write(*,1012) mousefqso,nutc
1012 format('!',i3,9x,i5.4,' ')
endif
endif
if(ndphi.eq.1 .and.iloop.lt.12) then
iloop=iloop+1
go to 2
endif
if(ndphi.eq.1 .and.iloop.eq.12) call getdphi(qphi)
if(nagain.eq.1) go to 999
enddo
! Trim the list and produce a sorted index and sizes of groups.
! (Should trimlist remove all but best SNR for given UTC and message content?)
call trimlist(sig,km,ftol,indx,nsiz,nz)
do i=1,km
done(i)=.false.
enddo
j=0
ilatest=-1
do n=1,nz
ifile0=0
do m=1,nsiz(n)
i=indx(j+m)
ifile=sig(i,1)
if(ifile.gt.ifile0 .and.msg(i).ne.blank) then
ilatest=i
ifile0=ifile
endif
enddo
i=ilatest
if(i.ge.1) then
if(.not.done(i)) then
done(i)=.true.
nutc=sig(i,2)
freq=sig(i,3)
sync1=sig(i,4)
dt=sig(i,5)
npol=nint(57.2957795*sig(i,6))
flip=sig(i,7)
sync2=sig(i,8)
nkv=sig(i,9)
nqual=min(sig(i,10),10.0)
! rms0=sig(i,11)
do k=1,5
a(k)=sig(i,12+k)
enddo
nhist=sig(i,18)
decoded=msg(i)
if(flip.lt.0.0) then
do i=22,1,-1
if(decoded(i:i).ne.' ') go to 10
enddo
stop 'Error in message format'
10 if(i.le.18) decoded(i+2:i+4)='OOO'
endif
mhz=fcenter !... +fadd ???
nkHz=nint(freq-foffset)-nfshift
f0=mhz+0.001*nkHz
ndf=nint(1000.0*(freq-foffset-(nkHz+nfshift)))
ndf0=nint(a(1))
ndf1=nint(a(2))
ndf2=nint(a(3))
nsync1=sync1
nsync2=nint(10.0*log10(sync2)) - 40 !### empirical ###
if(decoded(1:4).eq.'RO ' .or. decoded(1:4).eq.'RRR ' .or. &
decoded(1:4).eq.'73 ') nsync2=nsync2-6
if(nxant.ne.0) then
npol=npol-45
if(npol.lt.0) npol=npol+180
endif
! If Tx station's grid is in decoded message, compute optimum TxPol
i1=index(decoded,' ')
i2=index(decoded(i1+1:),' ') + i1
grid=' '
if(i2.ge.8 .and. i2.le.18) grid=decoded(i2+1:i2+4)//'mm'
ntxpol=0
cp=' '
if(xpol) then
if(grid(1:1).ge.'A' .and. grid(1:1).le.'R' .and. &
grid(2:2).ge.'A' .and. grid(2:2).le.'R' .and. &
grid(3:3).ge.'0' .and. grid(3:3).le.'9' .and. &
grid(4:4).ge.'0' .and. grid(4:4).le.'9') then
ntxpol=mod(npol-nint(2.0*dpol(mygrid,grid))+720,180)
if(nxant.eq.0) then
cp='H'
if(ntxpol.gt.45 .and. ntxpol.le.135) cp='V'
else
cp='/'
if(ntxpol.ge.90 .and. ntxpol.lt.180) cp='\\'
endif
endif
endif
write(26,1014) f0,ndf,ndf0,ndf1,ndf2,dt,npol,nsync1, &
nsync2,nutc,decoded,nkv,nqual,nhist,cp
write(21,1014) f0,ndf,ndf0,ndf1,ndf2,dt,npol,nsync1, &
nsync2,nutc,decoded,nkv,nqual,nhist
1014 format(f8.3,i5,3i3,f5.1,i4,i3,i4,i5.4,2x,a22,3i3,1x,a1)
endif
endif
j=j+nsiz(n)
enddo
write(26,1015) nutc
1015 format(39x,i4.4)
call flush(21)
call flush(26)
call display(nkeep,ftol)
ndecdone=2
999 close(23)
ndphi=0
nagain=0
mcall3b=mcall3a
return
end subroutine map65a

28
libm65/match.f90
View File

@ -1,28 +0,0 @@
subroutine match(s1,s2,nstart,nmatch)
character*(*) s1,s2
character s1a*29
nstart=-1
nmatch=0
n1=len_trim(s1)+1
n2=len(s2)
s1a=s1//' '
if(n2.ge.n1) then
do j=1,n2
n=0
do i=1,n1
k=j+i-1
if(k.gt.n2) k=k-n2
if(s2(k:k).eq.s1a(i:i)) n=n+1
enddo
if(n.gt.nmatch) then
nmatch=n
nstart=j
endif
enddo
endif
return
end subroutine match

167
libm65/moon2.f
View File

@ -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

73
libm65/moondop.f
View File

@ -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

75
libm65/mskdf.f90
View File

@ -1,75 +0,0 @@
subroutine mskdf(cdat,npts,t2,nfft1,f0,nfreeze,mousedf,ntol,dfx,snrsq2)
! Determine DF for a JTMS signal. Also find ferr, a measure of
! frequency differerence between 1st and 2nd harmonic.
! (Should be 0.000)
parameter (NZ=128*1024)
complex cdat(npts)
integer ntol
real sq(NZ)
real ccf(-2600:2600) !Correct limits?
real tmp(NZ)
complex c(NZ)
data nsps/8/
save c
df1=48000.0/nfft1
nh=nfft1/2
fac=1.0/(nfft1**2)
do i=1,npts
c(i)=fac*cdat(i)**2
enddo
c(npts+1:nfft1)=0.
call four2a(c,nfft1,1,-1,1)
! In the "doubled-frequencies" spectrum of squared cdat:
fa=2.0*(f0-400)
fb=2.0*(f0+400)
j0=nint(2.0*f0/df1)
ja=nint(fa/df1)
jb=nint(fb/df1)
jd=nfft1/nsps
do j=1,nh+1
sq(j)=real(c(j))**2 + aimag(c(j))**2
! if(j*df1.lt.6000.0) write(54,3009) j*df1,sq(j),db(sq(j))
!3009 format(3f12.3)
enddo
ccf=0.
do j=ja,jb
ccf(j-j0-1)=sq(j)+sq(j+jd)
enddo
call pctile(ccf(ja-j0-1),tmp,jb-ja+1,50,base)
ccf=ccf/base
if(NFreeze.gt.0) then
fa=2.0*(f0+MouseDF-ntol)
fb=2.0*(f0+MouseDF+ntol)
endif
ja=nint(fa/df1)
jb=nint(fb/df1)
! rewind 51
smax=0.
do j=ja,jb
k=j-j0-1
if(ccf(k).gt.smax) then
smax=ccf(k)
jpk=j
endif
f=0.5*k*df1
! write(51,3002) f,ccf(k)
!3002 format(2f12.3)
enddo
! call flush(51)
fpk=(jpk-1)*df1
dfx=0.5*fpk-f0
snrsq2=smax
return
end subroutine mskdf

42
libm65/ping.f90
View File

@ -1,42 +0,0 @@
subroutine ping(s,nz,dtbuf,slim,wmin,pingdat,nping)
! Detect pings and make note of their start time, duration, and peak strength.
real*4 s(nz)
real*4 pingdat(3,100)
logical inside
nping=0
peak=0.
inside=.false.
snrlim=10.0**(0.1*slim) - 1.0
sdown=10.0*log10(0.25*snrlim+1.0)
i0=0 !Silence compiler warnings.
tstart=0.0
do i=2,nz
if(s(i).ge.slim .and. .not.inside) then
i0=i
tstart=i0*dtbuf
inside=.true.
peak=0.
endif
if(inside .and. s(i).gt.peak) then
peak=s(i)
endif
if(inside .and. (s(i).lt.sdown .or. i.eq.nz)) then
if(i.gt.i0) then
if(dtbuf*(i-i0).ge.wmin) then
if(nping.le.99) nping=nping+1
pingdat(1,nping)=tstart
pingdat(2,nping)=dtbuf*(i-i0)
pingdat(3,nping)=peak
endif
inside=.false.
peak=0.
endif
endif
enddo
return
end subroutine ping

70
libm65/recvpkt.f90
View File

@ -1,70 +0,0 @@
subroutine recvpkt(nsam,nblock2,userx_no,k,buf4,buf8,buf16)
! Reformat timf2 data from Linrad and stuff data into r*4 array dd().
parameter (NSMAX=60*96000) !Total sample intervals per minute
parameter (NFFT=32768)
integer*1 userx_no
real*4 d4,buf4(*) !(348)
real*8 d8,buf8(*) !(174)
complex*16 c16,buf16(*) !(87)
integer*2 jd(4),kd(2),nblock2
real*4 xd(4),yd(2)
real*8 fcenter
common/datcom/dd(4,5760000),ss(4,322,NFFT),savg(4,NFFT),fcenter,nutc,junk(34)
equivalence (kd,d4)
equivalence (jd,d8,yd)
equivalence (xd,c16)
if(nsam.eq.-1) then
! Move data from the UDP packet buffer into array dd().
if(userx_no.eq.-1) then
do i=1,174 !One RF channel, r*4 data
k=k+1
d8=buf8(i)
dd(1,k)=yd(1)
dd(2,k)=yd(2)
enddo
else if(userx_no.eq.1) then
do i=1,348 !One RF channel, i*2 data
k=k+1
d4=buf4(i)
dd(1,k)=kd(1)
dd(2,k)=kd(2)
enddo
else if(userx_no.eq.-2) then
do i=1,87 !Two RF channels, r*4 data
k=k+1
c16=buf16(i)
dd(1,k)=xd(1)
dd(2,k)=xd(2)
dd(3,k)=xd(3)
dd(4,k)=xd(4)
enddo
else if(userx_no.eq.2) then
do i=1,174 !Two RF channels, i*2 data
k=k+1
d8=buf8(i)
dd(1,k)=jd(1)
dd(2,k)=jd(2)
dd(3,k)=jd(3)
dd(4,k)=jd(4)
enddo
endif
else
if(userx_no.eq.1) then
do i=1,nsam !One RF channel, r*4 data
k=k+1
d4=buf4(i)
dd(1,k)=kd(1)
dd(2,k)=kd(2)
k=k+1
dd(1,k)=kd(1)
dd(2,k)=kd(2)
enddo
endif
endif
return
end subroutine recvpkt

14
libm65/rfile3a.f90
View File

@ -1,14 +0,0 @@
subroutine rfile3a(infile,ibuf,n,fcenter,ierr)
character*(*) infile
integer*8 ibuf(n)
real*8 fcenter
open(10,file=infile,access='stream',status='old',err=998)
read(10,end=998) (ibuf(i),i=1,n/8),fcenter
ierr=0
go to 999
998 ierr=1002
999 close(10)
return
end subroutine rfile3a

92
libm65/rtping.f90
View File

@ -1,92 +0,0 @@
subroutine rtping(dat,k,cfile6,MinSigdB,MouseDF,ntol,mycall)
! Called from datasink() every 2048 sample intervals (approx 43 ms).
! Detects pings (signal level above MinSigdB). When a ping ends, its
! MSK signal is synchronized and decoded.
parameter (NSMAX=30*48000)
parameter (NZMAX=703) !703 = NSMAX/2048
real dat(NSMAX) !Raw audio data
character*6 cfile6 !Time hhmmss at start of this Rx interval
character*12 mycall
real sig(NZMAX) !Sq-law detected signal, sampled at 43 ms
real sigdb(NZMAX) !Signal in dB, sampled at 43 ms
real tmp(NZMAX)
logical inside,pingFound
data k0/9999999/
save
if(k.lt.k0) then
inside=.false.
pingFound=.false.
j0=0
t1=0.
width=0.
peak=0.
tpk=0.
dt=1.0/48000.0
kstep=2048
sdt=dt*kstep
wmin=0.043
endif
k0=k
slim=MinSigdB
snrlim=10.0**(0.1*slim) - 1.0
sdown=10.0*log10(0.25*snrlim+1.0)
! Find signal power
j=k/kstep
sig(j)=dot_product(dat(k-kstep+1:k),dat(k-kstep+1:k))/kstep
if(j.lt.20) return
! Determine baseline noise level
if(mod(j,20).eq.0) call pctile (sig,tmp,j,50,base)
sigdb(j)=db(sig(j)/base) ! (S+N)/N in dB
! write(72,3001) j*sdt,base,sig(j),sigdb(j)
!3001 format(f10.3,3f12.6)
if(sigdb(j).ge.slim .and. .not.inside) then
j0=j !Mark the start of a ping
t1=j0*sdt
inside=.true.
peak=0.
endif
if(inside .and. sigdb(j).gt.peak) then
peak=sigdb(j) !Save peak strength
tpk=j*sdt ! and time of peak
endif
if(inside .and. (sigdb(j).lt.sdown .or. j.eq.NZMAX)) then
width=(j-j0)*sdt !Save ping width
if(width.ge.wmin) pingFound=.true.
endif
if(.not.pingFound) return
! A ping was found! Assemble a signal report:
nwidth=0
if(width.ge.0.04) nwidth=1
if(width.ge.0.12) nwidth=2
if(width.gt.1.00) nwidth=3
nstrength=6
if(peak.ge.11.0) nstrength=7
if(peak.ge.17.0) nstrength=8
if(peak.ge.23.0) nstrength=9
nrpt=10*nwidth + nstrength
mswidth=10*nint(100.0*width)
i1=(t1-0.02)/dt
if(i1.lt.1) i1=1
iz=nint((width+0.02)/dt) + 1
iz=min(iz,k+1-i1,2*48000) !Length of ping in samples
call jtmsk(dat(i1),iz,cfile6,tpk,mswidth,int(peak),nrpt, &
nfreeze,DFTolerance,MouseDF,pick,mycall)
pingFound=.false.
inside=.false.
return
end subroutine rtping

42
libm65/s3avg.f90
View File

@ -1,42 +0,0 @@
subroutine s3avg(nsave,mode65,nutc,ndf,xdt,npol,s3,nkv,decoded)
real s3(64,63),s3b(64,63)
real s3a(64,63,32)
integer iutc(32),idf(32),ipol(32)
real dt(32)
character*22 decoded
logical ltext
save
n=nsave
iutc(n)=nutc
idf(n)=ndf
ipol(n)=npol
dt(n)=xdt
s3a(1:64,1:63,n)=s3
s3b=0.
nsum=0
idfdiff=100
dtdiff=0.2
do i=1,n
if(mod(iutc(i),2).ne.mod(nutc,2)) cycle
if(abs(ndf-idf(i)).gt.idfdiff) cycle
if(abs(xdt-dt(i)).gt.dtdiff) cycle
s3b=s3b + s3a(1:64,1:63,i)
nsum=nsum+1
enddo
decoded=' '
if(nsum.ge.2) then
nadd=mode65*nsum
call extract(s3b,nadd,ncount,nhist,decoded,ltext) !Extract the message
nkv=nsum
if(ncount.lt.0) then
nkv=0
decoded=' '
endif
endif
return
end subroutine s3avg

96
libm65/setup65.f
View File

@ -1,96 +0,0 @@
subroutine setup65
C Defines arrays related to the JT65 pseudo-random synchronizing pattern.
C Executed at program start.
integer nprc(126)
common/prcom/pr(126),mdat(126),mref(126,2),mdat2(126),mref2(126,2)
C JT65
data nprc/
+ 1,0,0,1,1,0,0,0,1,1,1,1,1,1,0,1,0,1,0,0,
+ 0,1,0,1,1,0,0,1,0,0,0,1,1,1,0,0,1,1,1,1,
+ 0,1,1,0,1,1,1,1,0,0,0,1,1,0,1,0,1,0,1,1,
+ 0,0,1,1,0,1,0,1,0,1,0,0,1,0,0,0,0,0,0,1,
+ 1,0,0,0,0,0,0,0,1,1,0,1,0,0,1,0,1,1,0,1,
+ 0,1,0,1,0,0,1,1,0,0,1,0,0,1,0,0,0,0,1,1,
+ 1,1,1,1,1,1/
data mr2/0/ !Silence compiler warning
C Put the appropriate pseudo-random sequence into pr
nsym=126
do i=1,nsym
pr(i)=2*nprc(i)-1
enddo
C Determine locations of data and reference symbols
k=0
mr1=0
do i=1,nsym
if(pr(i).lt.0.0) then
k=k+1
mdat(k)=i
else
mr2=i
if(mr1.eq.0) mr1=i
endif
enddo
nsig=k
C Determine the reference symbols for each data symbol.
do k=1,nsig
m=mdat(k)
mref(k,1)=mr1
do n=1,10 !Get ref symbol before data
if((m-n).gt.0) then
if (pr(m-n).gt.0.0) go to 10
endif
enddo
go to 12
10 mref(k,1)=m-n
12 mref(k,2)=mr2
do n=1,10 !Get ref symbol after data
if((m+n).le.nsym) then
if (pr(m+n).gt.0.0) go to 20
endif
enddo
go to 22
20 mref(k,2)=m+n
22 enddo
C Now do it all again, using opposite logic on pr(i)
k=0
mr1=0
do i=1,nsym
if(pr(i).gt.0.0) then
k=k+1
mdat2(k)=i
else
mr2=i
if(mr1.eq.0) mr1=i
endif
enddo
nsig=k
do k=1,nsig
m=mdat2(k)
mref2(k,1)=mr1
do n=1,10
if((m-n).gt.0) then
if (pr(m-n).lt.0.0) go to 110
endif
enddo
go to 112
110 mref2(k,1)=m-n
112 mref2(k,2)=mr2
do n=1,10
if((m+n).le.nsym) then
if (pr(m+n).lt.0.0) go to 120
endif
enddo
go to 122
120 mref2(k,2)=m+n
122 enddo
return
end

51
libm65/setupmsk.f90
View File

@ -1,51 +0,0 @@
subroutine setupmsk(cw,cwb)
! Calculate the JTMS character waveforms.
complex cw(168,0:63)
complex cwb(168)
integer nb(7)
! real*8 twopi,dt,f0,f1
character cc*64
! 1 2 3 4 5 6
! 0123456789012345678901234567890123456789012345678901234567890123
data cc/'0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ./?- _ @'/
nsps=24
twopi=8.d0*atan(1.d0)
dt=1.d0/48000.d0 !Sample interval
f0=1000.d0
f1=2000.d0
dphi0=twopi*dt*f0
dphi1=twopi*dt*f1
do i=0,63
k=0
m=0
do n=5,0,-1 !Each character gets 6+1 bits
k=k+1
nb(k)=iand(1,ishft(i,-n))
m=m+nb(k)
enddo
k=k+1
nb(k) = 1 - iand(m,1) !Insert odd parity bit
phi=0.
j=0
do k=1,7 !Generate the waveform
if(nb(k).eq.0) then
dphi=dphi0
else
dphi=dphi1
endif
do ii=1,nsps
j=j+1
phi=phi+dphi
cw(j,i)=cmplx(cos(phi),sin(phi))
enddo
enddo
enddo
cwb=cw(1:168,44)
return
end subroutine setupmsk

106
libm65/specjtms.f90
View File

@ -1,106 +0,0 @@
subroutine specjtms(k,px,pxsmo,spk0,f0)
! Compute noise level and 2D spectra, for GUI display.
parameter (NSMAX=30*48000)
parameter (MAXFFT=8192)
integer*2 id
real x(MAXFFT)
complex cx(MAXFFT),cx2(MAXFFT)
logical first
common/mscom/id(1440000),s1(215),s2(215),kin,ndiskdat,kline,nutc
common/jt8com/nn(1800*12000),ss(184,32768),savg(32768)
data first/.true./
save
if(first) then
pi=4.0*atan(1.0)
twopi=2.0*pi
kstep=2048
first=.false.
sqsmo=0.
endif
if(k.lt.kstep .or. k.gt.NSMAX) return
t=k/48000.0
nfft=4096
df=48000.0/nfft
nh=nfft/2
ib=k
ia=k-kstep+1
i0=k-nfft+1
sq=0.
do i=ia,ib
d=id(i)
sq=sq + d*d
enddo
sq=sq/2048.0
sqsmo=0.95*sqsmo + 0.05*sq
rms=sqrt(sq)
px=db(sq) - 23.0
pxsmo=db(sqsmo) - 23.0
if(i0.gt.0) then
do i=i0,ia-1
d=id(i)
sq=sq + d*d
enddo
endif
sq0=sq
! write(13,1010) t,rms,sq,px,pxsmo
!1010 format(5f12.3)
if(k.lt.nfft) return
x(1:nfft)=id(i0:ib)
fac=0.002/nfft
cx=fac*x
call four2a(cx,nfft,1,-1,1) !Forward c2c FFT
iz=nint(2500.0/df)
do i=1,iz !Save spectrum for plot
s1(i)=real(cx(i+1))**2 + aimag(cx(i+1))**2
enddo
cx(1)=0.5*cx(1)
cx(nh+2:nfft)=0.
call four2a(cx,nfft,1,1,1) !Inverse c2c FFT
cx2(1:nfft)=cx(1:nfft)*cx(1:nfft)
cx2(nfft+1:)=0.0
nfft=8192
df=48000.0/nfft
call four2a(cx2,nfft,1,-1,1) !Forward c2c FFT of cx2
! j0=nint(2.0*1428.57/df)
j0=nint(2.0*1500.0/df)
ja=j0-107
jb=j0+107
do j=ja,jb
s2(j-ja+1)=1.e-4*(real(cx2(j))**2 + aimag(cx2(j))**2)
enddo
spk0=0.
fac=(5e8/sq0)**2
s2=fac*s2
do j=1,215
if(s2(j).gt.spk0) then
spk0=s2(j)
f=(j+ja-1)*df
f0=0.5*(f-3000.0)
phi0=0.5*atan2(aimag(cx2(j)),real(cx2(j)))
endif
enddo
spk0=0.5*db(spk0) - 7.0
kline=k/2048
! write(14,3001) k/2048,spk0,f0,phi0
!3001 format(i8,3f12.3)
! call flush(14)
return
end subroutine specjtms

88
libm65/sun.f
View File

@ -1,88 +0,0 @@
subroutine sun(y,m,DD,UT,lon,lat,RA,Dec,LST,Az,El,mjd,day)
implicit none
integer y !Year
integer m !Month
integer DD !Day
integer mjd !Modified Julian Date
real UT !UTC in hours
real RA,Dec !RA and Dec of sun
C NB: Double caps here are single caps in the writeup.
C Orbital elements of the Sun (also N=0, i=0, a=1):
real w !Argument of perihelion
real e !Eccentricity
real MM !Mean anomaly
real Ls !Mean longitude
C Other standard variables:
real v !True anomaly
real EE !Eccentric anomaly
real ecl !Obliquity of the ecliptic
real d !Ephemeris time argument in days
real r !Distance to sun, AU
real xv,yv !x and y coords in ecliptic
real lonsun !Ecliptic long and lat of sun
C Ecliptic coords of sun (geocentric)
real xs,ys
C Equatorial coords of sun (geocentric)
real xe,ye,ze
real lon,lat
real GMST0,LST,HA
real xx,yy,zz
real xhor,yhor,zhor
real Az,El
real day
real rad
data rad/57.2957795/
C Time in days, with Jan 0, 2000 equal to 0.0:
d=367*y - 7*(y+(m+9)/12)/4 + 275*m/9 + DD - 730530 + UT/24.0
mjd=d + 51543
ecl = 23.4393 - 3.563e-7 * d
C Compute updated orbital elements for Sun:
w = 282.9404 + 4.70935e-5 * d
e = 0.016709 - 1.151e-9 * d
MM = mod(356.0470d0 + 0.9856002585d0 * d + 360000.d0,360.d0)
Ls = mod(w+MM+720.0,360.0)
EE = MM + e*rad*sin(MM/rad) * (1.0 + e*cos(M/rad))
EE = EE - (EE - e*rad*sin(EE/rad)-MM) / (1.0 - e*cos(EE/rad))
xv = cos(EE/rad) - e
yv = sqrt(1.0-e*e) * sin(EE/rad)
v = rad*atan2(yv,xv)
r = sqrt(xv*xv + yv*yv)
lonsun = mod(v + w + 720.0,360.0)
C Ecliptic coordinates of sun (rectangular):
xs = r * cos(lonsun/rad)
ys = r * sin(lonsun/rad)
C Equatorial coordinates of sun (rectangular):
xe = xs
ye = ys * cos(ecl/rad)
ze = ys * sin(ecl/rad)
C RA and Dec in degrees:
RA = rad*atan2(ye,xe)
Dec = rad*atan2(ze,sqrt(xe*xe + ye*ye))
GMST0 = (Ls + 180.0)/15.0
LST = mod(GMST0+UT+lon/15.0+48.0,24.0) !LST in hours
HA = 15.0*LST - RA !HA in degrees
xx = cos(HA/rad)*cos(Dec/rad)
yy = sin(HA/rad)*cos(Dec/rad)
zz = sin(Dec/rad)
xhor = xx*sin(lat/rad) - zz*cos(lat/rad)
yhor = yy
zhor = xx*cos(lat/rad) + zz*sin(lat/rad)
Az = mod(rad*atan2(yhor,xhor) + 180.0 + 360.0,360.0)
El = rad*asin(zhor)
day=d-1.5
return
end

5
libm65/sync9.f90
View File

@ -34,8 +34,6 @@ subroutine sync9(ss,tstep,f0a,df3,lagpk,fpk)
enddo
fpk=f0a + (npk-1)*df3
write(*,1010) lagpk,npk,fpk
1010 format('lagpk:',i4,' npk:',i6,' fpk:',f8.2)
do lag=-lagmax,lagmax
sum=0.
@ -43,10 +41,7 @@ subroutine sync9(ss,tstep,f0a,df3,lagpk,fpk)
k=ii(i) + lag
if(k.ge.1) sum=sum + ss(k,npk)
enddo
! write(73,3000) lag,sum
!3000 format(i3,f12.3)
enddo
flush(73)
return
end subroutine sync9

38
libm65/syncmsk.f90
View File

@ -1,38 +0,0 @@
subroutine syncmsk(cdat,npts,cwb,r,i1)
! Establish character sync within a JTMS ping.
complex cdat(npts) !Analytic signal
complex cwb(168) !Complex waveform for 'space'
real r(60000)
real tmp(60000)
integer hist(168),hmax(1)
complex z
r=0.
jz=npts-168+1
do j=1,jz
z=0.
ss=0.
do i=1,168
ss=ss + abs(cdat(i+j-1)) !Total power
z=z + cdat(i+j-1)*conjg(cwb(i)) !Signal matching <space>
enddo
r(j)=abs(z)/ss !Goodness-of-fit to <space>
! write(52,3001) j/168.0,r(j),cdat(j)
!3001 format(4f12.3)
enddo
ncut=99.0*float(jz-10)/float(jz)
call pctile(r,tmp,jz,ncut,rlim)
hist=0
do j=1,jz
k=mod(j-1,168)+1
if(r(j).gt.rlim) hist(k)=hist(k)+1
enddo
hmax=maxloc(hist)
i1=hmax(1)
return
end subroutine syncmsk

14
libm65/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 day4,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

514
libm65/tmoonsub.c
View File

@ -1,514 +0,0 @@
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define RADS 0.0174532925199433
#define DEGS 57.2957795130823
#define TPI 6.28318530717959
#define PI 3.1415927
/* ratio of earth radius to astronomical unit */
#define ER_OVER_AU 0.0000426352325194252
/* all prototypes here */
double getcoord(int coord);
void getargs(int argc, char *argv[], int *y, int *m, double *tz, double *glong, double *glat);
double range(double y);
double rangerad(double y);
double days(int y, int m, int dn, double hour);
double days_(int *y, int *m, int *dn, double *hour);
void moonpos(double, double *, double *, double *);
void sunpos(double , double *, double *, double *);
double moontransit(int y, int m, int d, double timezone, double glat, double glong, int *nt);
double atan22(double y, double x);
double epsilon(double d);
void equatorial(double d, double *lon, double *lat, double *r);
void ecliptic(double d, double *lon, double *lat, double *r);
double gst(double d);
void topo(double lst, double glat, double *alp, double *dec, double *r);
double alt(double glat, double ha, double dec);
void libration(double day, double lambda, double beta, double alpha, double *l, double *b, double *p);
void illumination(double day, double lra, double ldec, double dr, double sra, double sdec, double *pabl, double *ill);
int daysinmonth(int y, int m);
int isleap(int y);
void tmoonsub_(double *day, double *glat, double *glong, double *moonalt,
double *mrv, double *l, double *b, double *paxis);
static const char
*usage = " Usage: tmoon date[yyyymm] timz[+/-h.hh] long[+/-dddmm] lat[+/-ddmm]\n"
"example: tmoon 200009 0 -00155 5230\n";
/*
getargs() gets the arguments from the command line, does some basic error
checking, and converts arguments into numerical form. Arguments are passed
back in pointers. Error messages print to stderr so re-direction of output
to file won't leave users blind. Error checking prints list of all errors
in a command line before quitting.
*/
void getargs(int argc, char *argv[], int *y, int *m, double *tz,
double *glong, double *glat) {
int date, latitude, longitude;
int mflag = 0, yflag = 0, longflag = 0, latflag = 0, tzflag = 0;
int longminflag = 0, latminflag = 0, dflag = 0;
/* if not right number of arguments, then print example command line */
if (argc !=5) {
fprintf(stderr, usage);
exit(EXIT_FAILURE);
}
date = atoi(argv[1]);
*y = date / 100;
*m = date - *y * 100;
*tz = (double) atof(argv[2]);
longitude = atoi(argv[3]);
latitude = atoi(argv[4]);
*glong = RADS * getcoord(longitude);
*glat = RADS * getcoord(latitude);
/* set a flag for each error found */
if (*m > 12 || *m < 1) mflag = 1;
if (*y > 2500) yflag = 1;
if (date < 150001) dflag = 1;
if (fabs((float) *glong) > 180 * RADS) longflag = 1;
if (abs(longitude) % 100 > 59) longminflag = 1;
if (fabs((float) *glat) > 90 * RADS) latflag = 1;
if (abs(latitude) % 100 > 59) latminflag = 1;
if (fabs((float) *tz) > 12) tzflag = 1;
/* print all the errors found */
if (dflag == 1) {
fprintf(stderr, "date: dates must be in form yyyymm, gregorian, and later than 1500 AD\n");
}
if (yflag == 1) {
fprintf(stderr, "date: too far in future - accurate from 1500 to 2500\n");
}
if (mflag == 1) {
fprintf(stderr, "date: month must be in range 0 to 12, eg - August 2000 is entered as 200008\n");
}
if (tzflag == 1) {
fprintf(stderr, "timz: must be in range +/- 12 hours, eg -6 for Chicago\n");
}
if (longflag == 1) {
fprintf(stderr, "long: must be in range +/- 180 degrees\n");
}
if (longminflag == 1) {
fprintf(stderr, "long: last two digits are arcmin - max 59\n");
}
if (latflag == 1) {
fprintf(stderr, " lat: must be in range +/- 90 degrees\n");
}
if (latminflag == 1) {
fprintf(stderr, " lat: last two digits are arcmin - max 59\n");
}
/* quits if one or more flags set */
if (dflag + mflag + yflag + longflag + latflag + tzflag + longminflag + latminflag > 0) {
exit(EXIT_FAILURE);
}
}
/*
returns coordinates in decimal degrees given the
coord as a ddmm value stored in an integer.
*/
double getcoord(int coord) {
int west = 1;
double glg, deg;
if (coord < 0) west = -1;
glg = fabs((double) coord/100);
deg = floor(glg);
glg = west* (deg + (glg - deg)*100 / 60);
return(glg);
}
/*
days() takes the year, month, day in the month and decimal hours
in the day and returns the number of days since J2000.0.
Assumes Gregorian calendar.
*/
double days(int y, int m, int d, double h) {
int a, b;
double day;
/*
The lines below work from 1900 march to feb 2100
a = 367 * y - 7 * (y + (m + 9) / 12) / 4 + 275 * m / 9 + d;
day = (double)a - 730531.5 + hour / 24;
*/
/* These lines work for any Gregorian date since 0 AD */
if (m ==1 || m==2) {
m +=12;
y -= 1;
}
a = y / 100;
b = 2 - a + a/4;
day = floor(365.25*(y + 4716)) + floor(30.6001*(m + 1))
+ d + b - 1524.5 - 2451545 + h/24;
return(day);
}
double days_(int *y0, int *m0, int *d0, double *h0)
{
return days(*y0,*m0,*d0,*h0);
}
/*
Returns 1 if y a leap year, and 0 otherwise, according
to the Gregorian calendar
*/
int isleap(int y) {
int a = 0;
if(y % 4 == 0) a = 1;
if(y % 100 == 0) a = 0;
if(y % 400 == 0) a = 1;
return(a);
}
/*
Given the year and the month, function returns the
number of days in the month. Valid for Gregorian
calendar.
*/
int daysinmonth(int y, int m) {
int b = 31;
if(m == 2) {
if(isleap(y) == 1) b= 29;
else b = 28;
}
if(m == 4 || m == 6 || m == 9 || m == 11) b = 30;
return(b);
}
/*
moonpos() takes days from J2000.0 and returns ecliptic coordinates
of moon in the pointers. Note call by reference.
This function is within a couple of arcminutes most of the time,
and is truncated from the Meeus Ch45 series, themselves truncations of
ELP-2000. Returns moon distance in earth radii.
Terms have been written out explicitly rather than using the
table based method as only a small number of terms is
retained.
*/
void moonpos(double d, double *lambda, double *beta, double *rvec) {
double dl, dB, dR, L, D, M, M1, F, e, lm, bm, rm, t;
t = d / 36525;
L = range(218.3164591 + 481267.88134236 * t) * RADS;
D = range(297.8502042 + 445267.1115168 * t) * RADS;
M = range(357.5291092 + 35999.0502909 * t) * RADS;
M1 = range(134.9634114 + 477198.8676313 * t - .008997 * t * t) * RADS;
F = range(93.27209929999999 + 483202.0175273 * t - .0034029*t*t)*RADS;
e = 1 - .002516 * t;
dl = 6288774 * sin(M1);
dl += 1274027 * sin(2 * D - M1);
dl += 658314 * sin(2 * D);
dl += 213618 * sin(2 * M1);
dl -= e * 185116 * sin(M);
dl -= 114332 * sin(2 * F) ;
dl += 58793 * sin(2 * D - 2 * M1);
dl += e * 57066 * sin(2 * D - M - M1) ;
dl += 53322 * sin(2 * D + M1);
dl += e * 45758 * sin(2 * D - M);
dl -= e * 40923 * sin(M - M1);
dl -= 34720 * sin(D) ;
dl -= e * 30383 * sin(M + M1) ;
dl += 15327 * sin(2 * D - 2 * F) ;
dl -= 12528 * sin(M1 + 2 * F);
dl += 10980 * sin(M1 - 2 * F);
lm = rangerad(L + dl / 1000000 * RADS);
dB = 5128122 * sin(F);
dB += 280602 * sin(M1 + F);
dB += 277693 * sin(M1 - F);
dB += 173237 * sin(2 * D - F);
dB += 55413 * sin(2 * D - M1 + F);
dB += 46271 * sin(2 * D - M1 - F);
dB += 32573 * sin(2 * D + F);
dB += 17198 * sin(2 * M1 + F);
dB += 9266 * sin(2 * D + M1 - F);
dB += 8822 * sin(2 * M1 - F);
dB += e * 8216 * sin(2 * D - M - F);
dB += 4324 * sin(2 * D - 2 * M1 - F);
bm = dB / 1000000 * RADS;
dR = -20905355 * cos(M1);
dR -= 3699111 * cos(2 * D - M1);
dR -= 2955968 * cos(2 * D);
dR -= 569925 * cos(2 * M1);
dR += e * 48888 * cos(M);
dR -= 3149 * cos(2 * F);
dR += 246158 * cos(2 * D - 2 * M1);
dR -= e * 152138 * cos(2 * D - M - M1) ;
dR -= 170733 * cos(2 * D + M1);
dR -= e * 204586 * cos(2 * D - M);
dR -= e * 129620 * cos(M - M1);
dR += 108743 * cos(D);
dR += e * 104755 * cos(M + M1);
dR += 79661 * cos(M1 - 2 * F);
rm = 385000.56 + dR / 1000;
*lambda = lm;
*beta = bm;
/* distance to Moon must be in Earth radii */
*rvec = rm / 6378.14;
}
/*
topomoon() takes the local siderial time, the geographical
latitude of the observer, and pointers to the geocentric
equatorial coordinates. The function overwrites the geocentric
coordinates with topocentric coordinates on a simple spherical
earth model (no polar flattening). Expects Moon-Earth distance in
Earth radii. Formulas scavenged from Astronomical Almanac 'low
precision formulae for Moon position' page D46.
*/
void topo(double lst, double glat, double *alp, double *dec, double *r) {
double x, y, z, r1;
x = *r * cos(*dec) * cos(*alp) - cos(glat) * cos(lst);
y = *r * cos(*dec) * sin(*alp) - cos(glat) * sin(lst);
z = *r * sin(*dec) - sin(glat);
r1 = sqrt(x*x + y*y + z*z);
*alp = atan22(y, x);
*dec = asin(z / r1);
*r = r1;
}
/*
moontransit() takes date, the time zone and geographic longitude
of observer and returns the time (decimal hours) of lunar transit
on that day if there is one, and sets the notransit flag if there
isn't. See Explanatory Supplement to Astronomical Almanac
section 9.32 and 9.31 for the method.
*/
double moontransit(int y, int m, int d, double tz, double glat, double glong, int *notransit) {
double hm, ht, ht1, lon, lat, rv, dnew, lst;
int itcount;
ht1 = 180 * RADS;
ht = 0;
itcount = 0;
*notransit = 0;
do {
ht = ht1;
itcount++;
dnew = days(y, m, d, ht * DEGS/15) - tz/24;
lst = gst(dnew) + glong;
/* find the topocentric Moon ra (hence hour angle) and dec */
moonpos(dnew, &lon, &lat, &rv);
equatorial(dnew, &lon, &lat, &rv);
topo(lst, glat, &lon, &lat, &rv);
hm = rangerad(lst - lon);
ht1 = rangerad(ht - hm);
/* if no convergence, then no transit on that day */
if (itcount > 30) {
*notransit = 1;
break;
}
}
while (fabs(ht - ht1) > 0.04 * RADS);
return(ht1);
}
/*
Calculates the selenographic coordinates of either the sub Earth point
(optical libration) or the sub-solar point (selen. coords of centre of
bright hemisphere). Based on Meeus chapter 51 but neglects physical
libration and nutation, with some simplification of the formulas.
*/
void libration(double day, double lambda, double beta, double alpha, double *l, double *b, double *p) {
double i, f, omega, w, y, x, a, t, eps;
t = day / 36525;
i = 1.54242 * RADS;
eps = epsilon(day);
f = range(93.2720993 + 483202.0175273 * t - .0034029 * t * t) * RADS;
omega = range(125.044555 - 1934.1361849 * t + .0020762 * t * t) * RADS;
w = lambda - omega;
y = sin(w) * cos(beta) * cos(i) - sin(beta) * sin(i);
x = cos(w) * cos(beta);
a = atan22(y, x);
*l = a - f;
/* kludge to catch cases of 'round the back' angles */
if (*l < -90 * RADS) *l += TPI;
if (*l > 90 * RADS) *l -= TPI;
*b = asin(-sin(w) * cos(beta) * sin(i) - sin(beta) * cos(i));
/* pa pole axis - not used for Sun stuff */
x = sin(i) * sin(omega);
y = sin(i) * cos(omega) * cos(eps) - cos(i) * sin(eps);
w = atan22(x, y);
*p = rangerad(asin(sqrt(x*x + y*y) * cos(alpha - w) / cos(*b)));
}
/*
Takes: days since J2000.0, eq coords Moon, ratio of moon to sun distance,
eq coords Sun
Returns: position angle of bright limb wrt NCP, percentage illumination
of Sun
*/
void illumination(double day, double lra, double ldec, double dr, double sra, double sdec, double *pabl, double *ill) {
double x, y, phi, i;
y = cos(sdec) * sin(sra - lra);
x = sin(sdec) * cos(ldec) - cos(sdec) * sin(ldec) * cos (sra - lra);
*pabl = atan22(y, x);
phi = acos(sin(sdec) * sin(ldec) + cos(sdec) * cos(ldec) * cos(sra-lra));
i = atan22(sin(phi) , (dr - cos(phi)));
*ill = 0.5*(1 + cos(i));
}
/*
sunpos() takes days from J2000.0 and returns ecliptic longitude
of Sun in the pointers. Latitude is zero at this level of precision,
but pointer left in for consistency in number of arguments.
This function is within 0.01 degree (1 arcmin) almost all the time
for a century either side of J2000.0. This is from the 'low precision
fomulas for the Sun' from C24 of Astronomical Alamanac
*/
void sunpos(double d, double *lambda, double *beta, double *rvec) {
double L, g, ls, bs, rs;
L = range(280.461 + .9856474 * d) * RADS;
g = range(357.528 + .9856003 * d) * RADS;
ls = L + (1.915 * sin(g) + .02 * sin(2 * g)) * RADS;
bs = 0;
rs = 1.00014 - .01671 * cos(g) - .00014 * cos(2 * g);
*lambda = ls;
*beta = bs;
*rvec = rs;
}
/*
this routine returns the altitude given the days since J2000.0
the hour angle and declination of the object and the latitude
of the observer. Used to find the Sun's altitude to put a letter
code on the transit time, and to find the Moon's altitude at
transit just to make sure that the Moon is visible.
*/
double alt(double glat, double ha, double dec) {
return(asin(sin(dec) * sin(glat) + cos(dec) * cos(glat) * cos(ha)));
}
/* returns an angle in degrees in the range 0 to 360 */
double range(double x) {
double a, b;
b = x / 360;
a = 360 * (b - floor(b));
if (a < 0)
a = 360 + a;
return(a);
}
/* returns an angle in rads in the range 0 to two pi */
double rangerad(double x) {
double a, b;
b = x / TPI;
a = TPI * (b - floor(b));
if (a < 0)
a = TPI + a;
return(a);
}
/*
gets the atan2 function returning angles in the right
order and range
*/
double atan22(double y, double x) {
double a;
a = atan2(y, x);
if (a < 0) a += TPI;
return(a);
}
/*
returns mean obliquity of ecliptic in radians given days since
J2000.0.
*/
double epsilon(double d) {
double t = d/ 36525;
return((23.4392911111111 - (t* (46.8150 + 0.00059*t)/3600)) *RADS);
}
/*
replaces ecliptic coordinates with equatorial coordinates
note: call by reference destroys original values
R is unchanged.
*/
void equatorial(double d, double *lon, double *lat, double *r) {
double eps, ceps, seps, l, b;
l = *lon;
b = * lat;
eps = epsilon(d);
ceps = cos(eps);
seps = sin(eps);
*lon = atan22(sin(l)*ceps - tan(b)*seps, cos(l));
*lat = asin(sin(b)*ceps + cos(b)*seps*sin(l));
}
/*
replaces equatorial coordinates with ecliptic ones. Inverse
of above, but used to find topocentric ecliptic coords.
*/
void ecliptic(double d, double *lon, double *lat, double *r) {
double eps, ceps, seps, alp, dec;
alp = *lon;
dec = *lat;
eps = epsilon(d);
ceps = cos(eps);
seps = sin(eps);
*lon = atan22(sin(alp)*ceps + tan(dec)*seps, cos(alp));
*lat = asin(sin(dec)*ceps - cos(dec)*seps*sin(alp));
}
/*
returns the siderial time at greenwich meridian as
an angle in radians given the days since J2000.0
*/
double gst( double d) {
double t = d / 36525;
double theta;
theta = range(280.46061837 + 360.98564736629 * d + 0.000387933 * t * t);
return(theta * RADS);
}
void tmoonsub_(double *day, double *glat, double *glong, double *moonalt,
double *mrv, double *l, double *b, double *paxis)
{
double mlambda, mbeta;
double malpha, mdelta;
double lst, mhr;
double tlambda, tbeta, trv;
lst = gst(*day) + *glong;
/* find Moon topocentric coordinates for libration calculations */
moonpos(*day, &mlambda, &mbeta, mrv);
malpha = mlambda;
mdelta = mbeta;
equatorial(*day, &malpha, &mdelta, mrv);
topo(lst, *glat, &malpha, &mdelta, mrv);
mhr = rangerad(lst - malpha);
*moonalt = alt(*glat, mhr, mdelta);
/* Optical libration and Position angle of the Pole */
tlambda = malpha;
tbeta = mdelta;
trv = *mrv;
ecliptic(*day, &tlambda, &tbeta, &trv);
libration(*day, tlambda, tbeta, malpha, l, b, paxis);
}

25
libm65/toxyz.f
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 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

24
libm65/tweak1.f90
View File

@ -1,24 +0,0 @@
subroutine tweak1(ca,jz,f0,cb)
! Shift frequency of analytic signal ca, with output to cb
complex ca(jz),cb(jz)
real*8 twopi
complex*16 w,wstep
data twopi/0.d0/
save twopi
if(twopi.eq.0.d0) twopi=8.d0*atan(1.d0)
w=1.d0
dphi=twopi*f0/11025.d0
wstep=cmplx(cos(dphi),sin(dphi))
x0=0.5*(jz+1)
s=2.0/jz
do i=1,jz
x=s*(i-x0)
w=w*wstep
cb(i)=w*ca(i)
enddo
return
end subroutine tweak1

29
libm65/twkfreq.f
View File

@ -1,29 +0,0 @@
subroutine twkfreq(c4aa,c4bb,n5,a)
complex c4aa(n5)
complex c4bb(n5)
real a(5)
complex w,wstep
data twopi/6.283185307/
C Apply AFC corrections to the c4aa and c4bb data
w=1.0
wstep=1.0
x0=0.5*(n5+1)
s=2.0/n5
do i=1,n5
x=s*(i-x0)
if(mod(i,1000).eq.1) then
p2=1.5*x*x - 0.5
! p3=2.5*(x**3) - 1.5*x
! p4=4.375*(x**4) - 3.75*(x**2) + 0.375
dphi=(a(1) + x*a(2) + p2*a(3)) * (twopi/1378.125)
wstep=cmplx(cos(dphi),sin(dphi))
endif
w=w*wstep
c4aa(i)=w*c4aa(i)
c4bb(i)=w*c4bb(i)
enddo
return
end

70
libm65/wrapkarn.c
View File

@ -1,70 +0,0 @@
#include <math.h>
#include <stdio.h>
#include <float.h>
#include <limits.h>
#include <stdlib.h>
#include "rs.h"
static void *rs;
static int first=1;
void rs_encode_(int *dgen, int *sent)
// Encode JT65 data dgen[12], producing sent[63].
{
int dat1[12];
int b[51];
int i;
if(first) {
// Initialize the JT65 codec
rs=init_rs_int(6,0x43,3,1,51,0);
first=0;
}
// Reverse data order for the Karn codec.
for(i=0; i<12; i++) {
dat1[i]=dgen[11-i];
}
// Compute the parity symbols
encode_rs_int(rs,dat1,b);
// Move parity symbols and data into sent[] array, in reverse order.
for (i = 0; i < 51; i++) sent[50-i] = b[i];
for (i = 0; i < 12; i++) sent[i+51] = dat1[11-i];
}
void rs_decode_(int *recd0, int *era0, int *numera0, int *decoded, int *nerr)
// Decode JT65 received data recd0[63], producing decoded[12].
// Erasures are indicated in era0[numera]. The number of corrected
// errors is *nerr. If the data are uncorrectable, *nerr=-1 is returned.
{
int numera;
int i;
int era_pos[50];
int recd[63];
if(first) {
rs=init_rs_int(6,0x43,3,1,51,0);
first=0;
}
numera=*numera0;
for(i=0; i<12; i++) recd[i]=recd0[62-i];
for(i=0; i<51; i++) recd[12+i]=recd0[50-i];
if(numera)
for(i=0; i<numera; i++) era_pos[i]=era0[i];
*nerr=decode_rs_int(rs,recd,era_pos,numera);
for(i=0; i<12; i++) decoded[i]=recd[11-i];
}
void rs_encode__(int *dgen, int *sent)
{
rs_encode_(dgen, sent);
}
void rs_decode__(int *recd0, int *era0, int *numera0, int *decoded, int *nerr)
{
rs_decode_(recd0, era0, numera0, decoded, nerr);
}

4
wsjtx.pro
View File

@ -54,7 +54,7 @@ RC_FILE = wsjtx.rc
unix {
INCLUDEPATH += $$quote(/usr/include/qwt-qt4)
LIBS += -lfftw3f /usr/lib/libgfortran.so.3
LIBS += ../wsjtx/libm65/libm65.a
LIBS += ../wsjtx/libm65/libjt9.a
LIBS += /usr/lib/libqwt-qt4.so
LIBS += -lportaudio
#LIBS +- -lusb
@ -62,7 +62,7 @@ LIBS += -lportaudio
win32 {
INCLUDEPATH += c:/qwt-6.0.1/include
LIBS += ../wsjtx/libm65/libm65.a
LIBS += ../wsjtx/libm65/libjt9.a
LIBS += ../wsjtx/libfftw3f_win.a
LIBS += ../QtSupport/palir-02.dll
LIBS += libwsock32