chris/js8call
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Fortran cleanup take two

Browse Source
This commit is contained in:
Jordan Sherer
2020-05-31 15:34:57 -04:00
parent df63f23c22
commit 365b9546de
78 changed files with 12 additions and 10439 deletions
Download Patch File Download Diff File Expand all files Collapse all files
lib/afc9.f90
-81
View File
@@ -1,81 +0,0 @@
subroutine afc9(c3a,npts,fsample,a,syncpk)
parameter (NZ2=1512)
complex c3a(0:NZ2-1)
complex c3(0:NZ2-1)
real a(3),deltaa(3)
a(1)=0. !f0
a(2)=0. !f1
a(3)=0. !f2
deltaa(1)=1.736
deltaa(2)=1.736
deltaa(3)=1.0
nterms=3
! Start the iteration
chisqr=0.
chisqr0=1.e6
c3=c3a
a3=a(3)
do iter=1,4
do j=1,nterms
if(a(3).ne.a3) call shft(c3a,a(3),a3,c3)
chisq1=fchisq(c3,npts,fsample,a)
fn=0.
delta=deltaa(j)
10 a(j)=a(j)+delta
if(a(3).ne.a3) call shft(c3a,a(3),a3,c3)
chisq2=fchisq(c3,npts,fsample,a)
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
if(a(3).ne.a3) call shft(c3a,a(3),a3,c3)
chisq3=fchisq(c3,npts,fsample,a)
if(chisq3.lt.chisq2) then
chisq1=chisq2
chisq2=chisq3
go to 20
endif
! Find minimum of parabola defined by last three points
delta=delta*(1./(1.+(chisq1-chisq2)/(chisq3-chisq2))+0.5)
a(j)=a(j)-delta
if(j.lt.3) deltaa(j)=deltaa(j)*fn/3.
! write(*,4000) iter,j,a,-chisq2
!4000 format(i1,i2,3f10.4,f11.3)
enddo
if(a(3).ne.a3) call shft(c3a,a(3),a3,c3)
chisqr=fchisq(c3,npts,fsample,a)
! write(*,4000) 0,0,a,-chisqr
if(chisqr/chisqr0.gt.0.99) exit
chisqr0=chisqr
enddo
syncpk=-chisqr
c3a=c3
! write(*,4001) a,syncpk
!4001 format(3x,3f10.4,f11.3)
return
end subroutine afc9
subroutine shft(c3a,a3a,a3,c3)
complex c3a(0:1359)
complex c3(0:1359)
a3=a3a
n=nint(a3)
c3=cshift(c3a,n)
if(n.gt.0) c3(1360-n:)=0.0
if(n.lt.0) c3(:n-1)=0.0
return
end subroutine shft
lib/ana932.f90
-21
View File
@@ -1,21 +0,0 @@
subroutine ana932(dat,npts0,cdat,npts)
real dat(npts0)
complex cdat(262145)
n=log(float(npts0))/log(2.0)
nfft1=2**(n+1)
nfft2=9*nfft1/32
df932=11025.0/nfft1
fac=2.0/nfft1
do i=1,npts0/2
cdat(i)=fac*cmplx(dat(2*i-1),dat(2*i))
enddo
cdat(npts0/2+1:nfft1/2)=0.
call four2a(cdat,nfft1,1,-1,0) !Forward r2c FFT
call four2a(cdat,nfft2,1,1,1) !Inverse c2c FFT
npts=npts0*9.0/32.0 !Downsampled data length
npts2=npts
return
end subroutine ana932
lib/analytic.f90
-75
View File
@@ -1,75 +0,0 @@
subroutine analytic(d,npts,nfft,c,pc,beq)
! Convert real data to analytic signal
parameter (NFFTMAX=1024*1024)
real d(npts) ! passband signal
real h(NFFTMAX/2) ! real BPF magnitude
real*8 pc(5),pclast(5) ! static phase coeffs
real*8 ac(5),aclast(5) ! amp coeffs
real*8 fp
complex corr(NFFTMAX/2) ! complex frequency-dependent correction
complex c(NFFTMAX) ! analytic signal
logical*1 beq ! boolean static equalizer flag
data nfft0/0/
data aclast/0.0,0.0,0.0,0.0,0.0/
data pclast/0.0,0.0,0.0,0.0,0.0/
! data ac/1.0,0.05532,0.11438,0.12918,0.09274/ ! amp coeffs for TS2000
data ac/1.0,0.0,0.0,0.0,0.0/
save corr,nfft0,h,ac,aclast,pclast,pi,t,beta
df=12000.0/nfft
nh=nfft/2
if( nfft.ne.nfft0 ) then
pi=4.0*atan(1.0)
t=1.0/2000.0
beta=0.1
do i=1,nh+1
ff=(i-1)*df
f=ff-1500.0
h(i)=1.0
if(abs(f).gt.(1-beta)/(2*t) .and. abs(f).le.(1+beta)/(2*t)) then
h(i)=h(i)*0.5*(1+cos((pi*t/beta )*(abs(f)-(1-beta)/(2*t))))
elseif( abs(f) .gt. (1+beta)/(2*t) ) then
h(i)=0.0
endif
enddo
nfft0=nfft
endif
if( any(aclast .ne. ac) .or. any(pclast .ne. pc) ) then
aclast=ac
pclast=pc
! write(*,3001) pc
!3001 format('Phase coeffs:',5f12.6)
do i=1,nh+1
ff=(i-1)*df
f=ff-1500.0
fp=f/1000.0
corr(i)=ac(1)+fp*(ac(2)+fp*(ac(3)+fp*(ac(4)+fp*ac(5))))
pd=fp*fp*(pc(3)+fp*(pc(4)+fp*pc(5))) ! ignore 1st two terms
corr(i)=corr(i)*cmplx(cos(pd),sin(pd))
enddo
endif
fac=2.0/nfft
c(1:npts)=fac*d(1:npts)
c(npts+1:nfft)=0.
call four2a(c,nfft,1,-1,1) !Forward c2c FFT
if( beq ) then
c(1:nh+1)=h(1:nh+1)*corr(1:nh+1)*c(1:nh+1)
else
c(1:nh+1)=h(1:nh+1)*c(1:nh+1)
endif
c(1)=0.5*c(1) !Half of DC term
c(nh+2:nfft)=0. !Zero the negative frequencies
call four2a(c,nfft,1,1,1) !Inverse c2c FFT
return
end subroutine analytic
lib/astro.f90
-107
View File
@@ -1,107 +0,0 @@
subroutine astro(nyear,month,nday,uth,freq8,Mygrid, &
NStation,MoonDX,AzSun,ElSun,AzMoon0,ElMoon0, &
ntsky,doppler00,doppler,dbMoon,RAMoon,DecMoon,HA,Dgrd,sd, &
poloffset,xnr,day,lon,lat,LST,techo)
! Computes astronomical quantities for display and tracking.
! NB: may want to smooth the Tsky map to 10 degrees or so.
character*6 MyGrid,HisGrid
real*8 freq8
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)
call MoonDopJPL(nyear,month,nday,uth,lon,lat,RAMoon,DecMoon, &
LST,HA,AzMoon,ElMoon,vr,techo)
RAMoon=rad*RAMoon
DecMoon=rad*DecMoon
dist=2.99792458d5*techo/2.d0
! 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)
doppler=-freq8*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.0d6/freq8)**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.0d6/freq8)**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 subroutine astro
lib/astro0.f90
-79
View File
@@ -1,79 +0,0 @@
subroutine astro0(nyear,month,nday,uth8,freq8,mygrid,hisgrid, &
AzSun8,ElSun8,AzMoon8,ElMoon8,AzMoonB8,ElMoonB8,ntsky,ndop,ndop00, &
dbMoon8,RAMoon8,DecMoon8,HA8,Dgrd8,sd8,poloffset8,xnr8,dfdt,dfdt0, &
width1,width2,xlst8,techo8)
parameter (DEGS=57.2957795130823d0)
character*6 mygrid,hisgrid
real*8 AzSun8,ElSun8,AzMoon8,ElMoon8,AzMoonB8,ElMoonB8
real*8 dbMoon8,RAMoon8,DecMoon8,HA8,Dgrd8,xnr8,dfdt,dfdt0,dt
real*8 sd8,poloffset8,width1,width2,xlst8
real*8 uth8,techo8,freq8
real*8 xl,b
common/librcom/xl(2),b(2)
data uth8z/0.d0/
save
uth=uth8
call astro(nyear,month,nday,uth,freq8,hisgrid,2,1, &
AzSun,ElSun,AzMoon,ElMoon,ntsky,doppler00,doppler, &
dbMoon,RAMoon,DecMoon,HA,Dgrd,sd,poloffset,xnr, &
day,xlon2,xlat2,xlst,techo)
AzMoonB8=AzMoon
ElMoonB8=ElMoon
xl2=xl(1)
xl2a=xl(2)
b2=b(1)
b2a=b(2)
call astro(nyear,month,nday,uth,freq8,mygrid,1,1, &
AzSun,ElSun,AzMoon,ElMoon,ntsky,doppler00,doppler, &
dbMoon,RAMoon,DecMoon,HA,Dgrd,sd,poloffset,xnr, &
day,xlon1,xlat1,xlst,techo)
xl1=xl(1)
xl1a=xl(2)
b1=b(1)
b1a=b(2)
techo8=techo
fghz=1.d-9*freq8
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
AzSun8=AzSun
ElSun8=ElSun
AzMoon8=AzMoon
ElMoon8=ElMoon
dbMoon8=dbMoon
RAMoon8=RAMoon/15.0
DecMoon8=DecMoon
HA8=HA
xlst8=xlst
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
lib/astrosub.f90
-55
View File
@@ -1,55 +0,0 @@
subroutine astrosub(nyear,month,nday,uth8,freq8,mygrid,hisgrid, &
AzSun8,ElSun8,AzMoon8,ElMoon8,AzMoonB8,ElMoonB8,ntsky,ndop,ndop00, &
RAMoon8,DecMoon8,Dgrd8,poloffset8,xnr8,techo8,width1,width2,bTx, &
AzElFileName,jpleph)
implicit real*8 (a-h,o-z)
character*6 mygrid,hisgrid,c1*1
character*6 AzElFileName*(*),jpleph*(*)
character*256 jpleph_file_name
logical*1 bTx
common/jplcom/jpleph_file_name
jpleph_file_name=jpleph
call astro0(nyear,month,nday,uth8,freq8,mygrid,hisgrid, &
AzSun8,ElSun8,AzMoon8,ElMoon8,AzMoonB8,ElMoonB8,ntsky,ndop,ndop00, &
dbMoon8,RAMoon8,DecMoon8,HA8,Dgrd8,sd8,poloffset8,xnr8,dfdt,dfdt0, &
width1,width2,xlst8,techo8)
if (len_trim(AzElFileName) .eq. 0) go to 999
imin=60*uth8
isec=3600*uth8
ih=uth8
im=mod(imin,60)
is=mod(isec,60)
open(15,file=AzElFileName,status='unknown',err=900)
c1='R'
nRx=1
if(bTx) then
c1='T'
nRx=0
endif
AzAux=0.
ElAux=0.
nfreq=freq8/1000000
doppler=ndop
doppler00=ndop00
write(15,1010,err=10) ih,im,is,AzMoon8,ElMoon8, &
ih,im,is,AzSun8,ElSun8, &
ih,im,is,AzAux,ElAux, &
nfreq,doppler,dfdt,doppler00,dfdt0,c1
! TXFirst,TRPeriod,poloffset,Dgrd,xnr,ave,rms,nRx
1010 format( &
i2.2,':',i2.2,':',i2.2,',',f5.1,',',f5.1,',Moon'/ &
i2.2,':',i2.2,':',i2.2,',',f5.1,',',f5.1,',Sun'/ &
i2.2,':',i2.2,':',i2.2,',',f5.1,',',f5.1,',Source'/ &
i5,',',f8.1,',',f8.2,',',f8.1,',',f8.2,',Doppler, ',a1)
! i1,',',i3,',',f8.1,','f8.1,',',f8.1,',',f12.3,',',f12.3,',',i1,',RPol')
10 close(15)
go to 999
900 print*,'Error opening azel.dat'
999 return
end subroutine astrosub
lib/averms.f90
-20
View File
@@ -1,20 +0,0 @@
subroutine averms(x,n,nskip,ave,rms)
real x(n)
integer ipk(1)
ns=0
s=0.
sq=0.
ipk=maxloc(x)
do i=1,n
if(abs(i-ipk(1)).gt.nskip) then
s=s + x(i)
sq=sq + x(i)**2
ns=ns+1
endif
enddo
ave=s/ns
rms=sqrt(sq/ns - ave*ave)
return
end subroutine averms
lib/badmsg.f90
-46
View File
@@ -1,46 +0,0 @@
subroutine badmsg(irc,dat,nc1,nc2,ng2)
! Get rid of a few QRA64 false decodes that cannot be correct messages.
integer dat(12) !Decoded message (as 12 integers)
ic1=ishft(dat(1),22) + ishft(dat(2),16) + ishft(dat(3),10)+ &
ishft(dat(4),4) + iand(ishft(dat(5),-2),15)
! Test for "......" or "CQ 000"
if(ic1.eq.262177560 .or. ic1.eq.262177563) then
irc=-1
return
endif
ic2=ishft(iand(dat(5),3),26) + ishft(dat(6),20) + &
ishft(dat(7),14) + ishft(dat(8),8) + ishft(dat(9),2) + &
iand(ishft(dat(10),-4),3)
ig=ishft(iand(dat(10),15),12) + ishft(dat(11),6) + dat(12)
! Test for blank, -01 to -30, R-01 to R-30, RO, RRR, 73
if(ig.ge.32401 .and. ig.le.32464) return
if(ig.ge.14220 .and. ig.le.14229) return !-41 to -50
if(ig.ge.14040 .and. ig.le.14049) return !-31 to -40
if(ig.ge.13320 .and. ig.le.13329) return !+00 to +09
if(ig.ge.13140 .and. ig.le.13149) return !+10 to +19
if(ig.ge.12960 .and. ig.le.12969) return !+20 to +29
if(ig.ge.12780 .and. ig.le.12789) return !+30 to +39
if(ig.ge.12600 .and. ig.le.12609) return !+40 to +49
if(ig.ge.12420 .and. ig.le.12429) return !R-41 to R-50
if(ig.ge.12240 .and. ig.le.12249) return !R-31 to R-40
if(ig.ge.11520 .and. ig.le.11529) return !R+00 to R+09
if(ig.ge.11340 .and. ig.le.11349) return !R+10 to R+19
if(ig.ge.11160 .and. ig.le.11169) return !R+20 to R+29
if(ig.ge.10980 .and. ig.le.10989) return !R+30 to R+39
if(ig.ge.10800 .and. ig.le.10809) return !R+40 to R+49
if(ic1.eq.nc1 .and. ic2.eq.nc2 .and. ng2.ne.32401 .and. ig.ne.ng2) irc=-1
return
end subroutine badmsg
lib/calibrate.f90
-96
View File
@@ -1,96 +0,0 @@
subroutine calibrate(data_dir,iz,a,b,rms,sigmaa,sigmab,irc)
! Average groups of frequency-calibration measurements, then fit a
! straight line for slope and intercept.
parameter (NZ=1000)
implicit real*8 (a-h,o-z)
character*(*) data_dir
character*256 infile,outfile
character*8 cutc,cutc1
character*1 c1
real*8 fd(NZ),deltaf(NZ),r(NZ),rmsd(NZ)
integer nn(NZ)
infile=trim(data_dir)//'fmt.all'
outfile=trim(data_dir)//'fcal2.out'
open(10,file=trim(infile),status='old',err=996)
open(12,file=trim(outfile),status='unknown',err=997)
nkhz0=0
sum=0.d0
sumsq=0.d0
n=0
j=0
do i=1,99999
read(10,*,end=10,err=995) cutc,nkHz,ncal,noffset,faudio,df,dblevel,snr
if((nkHz.ne.nkHz0) .and. i.ne.1) then
ave=sum/n
rms=0.d0
if(n.gt.1) then
rms=sqrt(abs(sumsq - sum*sum/n)/(n-1.d0))
endif
fMHz=0.001d0*nkHz0
j=j+1
fd(j)=fMHz
deltaf(j)=ave
r(j)=0.d0
rmsd(j)=rms
nn(j)=n
sum=0.d0
sumsq=0.d0
n=0
endif
dial_error=faudio-noffset
sum=sum + dial_error
sumsq=sumsq + dial_error**2
n=n+1
if(n.eq.1) then
cutc1=cutc
ncal0=ncal
endif
nkHz0=nkHz
enddo
10 ave=sum/n
rms=0.d0
if(n.gt.0) then
rms=sqrt((sumsq - sum*sum/n)/(n-1.d0))
endif
fMHz=0.001d0*nkHz
j=j+1
fd(j)=fMHz
deltaf(j)=ave
r(j)=0.d0
rmsd(j)=rms
nn(j)=n
iz=j
if(iz.lt.2) go to 998
call fitcal(fd,deltaf,r,iz,a,b,sigmaa,sigmab,rms)
write(12,1002)
1002 format(' Freq DF Meas Freq N rms Resid'/ &
' (MHz) (Hz) (MHz) (Hz) (Hz)'/ &
'----------------------------------------------------')
irc=0
do i=1,iz
fm=fd(i) + 1.d-6*deltaf(i)
c1=' '
if(rmsd(i).gt.1.0d0) c1='*'
write(12,1012) fd(i),deltaf(i),fm,nn(i),rmsd(i),r(i),c1
1012 format(f8.3,f9.3,f14.9,i4,f7.2,f9.3,1x,a1)
enddo
go to 999
995 irc=-4; iz=i; go to 999
996 irc=-1; go to 999
997 irc=-2; go to 999
998 irc=-3
999 continue
close(10)
close(12)
return
end subroutine calibrate
lib/ccf2.f90
-49
View File
@@ -1,49 +0,0 @@
subroutine ccf2(ss,nz,nflip,ccfbest,xlagpk)
! parameter (LAGMIN=-86,LAGMAX=258)
parameter (LAGMIN=-112,LAGMAX=258) ! Look for DT from -3.6s to +5.0s
real ss(nz)
real ccf(-LAGMAX:LAGMAX)
integer npr(126)
! 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=LAGMIN
lag2=LAGMAX
do lag=lag1,lag2
s0=0.
s1=0.
do i=1,126
j=16*(i-1)+1 + lag
if(j.ge.1 .and. j.le.nz-8) then
x=ss(j)
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
xlagpk=lagpk
endif
enddo
if( lagpk.gt.-LAGMAX .and. lagpk.lt.LAGMAX) then
call peakup(ccf(lagpk-1),ccf(lagpk),ccf(lagpk+1),dx)
xlagpk=lagpk+dx
endif
return
end subroutine ccf2
lib/chkcall.f90
-58
View File
@@ -1,58 +0,0 @@
subroutine chkcall(w,bc,cok)
! Check "w" to see if it could be a valid standard callsign or a valid
! compound callsign.
! Return base call "bc" and a logical "cok" indicator.
character w*13 !A putative callsign
character bc*6 !Base call (tentative)
character c*1
logical cok,isdigit,isletter
isdigit(c)=(ichar(c).ge.ichar('0')) .and. (ichar(c).le.ichar('9'))
isletter(c)=(ichar(c).ge.ichar('A')) .and. (ichar(c).le.ichar('Z'))
cok=.true.
bc=w(1:6)
n1=len_trim(w)
if(n1.gt.11) go to 100
if(index(w,'.').ge.1) go to 100
if(index(w,'+').ge.1) go to 100
if(index(w,'-').ge.1) go to 100
if(index(w,'?').ge.1) go to 100
if(n1.gt.6 .and. index(w,'/').le.0) go to 100
i0=index(w,'/')
if(max(i0-1,n1-i0).gt.6) go to 100 !Base call must be < 7 characters
if(i0.ge.2 .and. i0.le.n1-1) then !Extract base call from compound call
if(i0-1.le.n1-i0) bc=w(i0+1:n1)//' '
if(i0-1.gt.n1-i0) bc=w(1:i0-1)//' '
endif
nbc=len_trim(bc)
if(nbc.gt.6) go to 100 !Base call should have no more than 6 characters
! One of first two characters (c1 or c2) must be a letter
if((.not.isletter(bc(1:1))) .and. (.not.isletter(bc(2:2)))) go to 100
if(bc(1:1).eq.'Q') go to 100 !Calls don't start with Q
! Must have a digit in 2nd or 3rd position
i1=0
if(isdigit(bc(2:2))) i1=2
if(isdigit(bc(3:3))) i1=3
if(i1.eq.0) go to 100
! Callsign must have a suffix of 1-3 letters
if(i1.eq.nbc) go to 100
n=0
do i=i1+1,nbc
j=ichar(bc(i:i))
if(j.lt.ichar('A') .or. j.gt.ichar('Z')) go to 100
n=n+1
enddo
if(n.ge.1 .and. n.le.3) go to 200
100 cok=.false.
200 return
end subroutine chkcall
lib/chkhist.f90
-21
View File
@@ -1,21 +0,0 @@
subroutine chkhist(mrsym,nmax,ipk)
integer mrsym(63)
integer hist(0:63)
hist=0
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 subroutine chkhist
lib/code426.f90
-62
View File
@@ -1,62 +0,0 @@
program code426
parameter (MZ=26) !Number of 4-FSK symbols
parameter (JZMAX=64) !Desired number of codewords
integer ic(MZ,JZMAX),icsave(MZ)
real c(MZ)
character*12 arg
nargs=iargc()
if(nargs.ne.2) then
print*,'Usage: code426 <nmsgs> <iters>'
print*,'Example: code426 64 10000000'
go to 999
endif
call getarg(1,arg)
read(arg,*) nmsgs
call getarg(2,arg)
read(arg,*) iters
call init_random_seed()
open(13,file='code426.out',status='unknown')
write(*,1002) nmsgs,iters
write(13,1002) nmsgs,iters
1002 format('Nmsgs:',i4,' Iters:',i10/(66('-')))
do i=1,MZ !Create 4 mutually orthogonal codewords
ic(i,1)=mod(i-1,4)
ic(i,2)=mod(i,4)
ic(i,3)=mod(i+1,4)
ic(i,4)=mod(i+2,4)
enddo
do j=1,4 !Write them out
write(*,1000) j,MZ,ic(1:MZ,j)
write(13,1000) j,MZ,ic(1:MZ,j)
1000 format(2i5,3x,26i2)
enddo
do j=5,nmsgs !Find codewords up to j=nmsgs with maximum
npk=0 !distance from all the rest
do i=1,iters
call random_number(c) !Generate a random codeword candidate
ic(1:MZ,j)=int(4*c) !Convert real to integer
! nd=MZ
! do k=1,j-1 !Test candidate against all others in list
! n=count(ic(1:MZ,j).ne.ic(1:MZ,k))
! nd=min(n,nd)
! enddo
call dist426(ic,j,mind)
if(mind.gt.npk) then
npk=mind
icsave=ic(1:MZ,j) !Best candidate so far, save it
! if(npk.ge.19) exit !It won't get any better...
endif
enddo
write(*,1000) j,npk,ic(1:MZ,j)
write(13,1000) j,npk,ic(1:MZ,j)
enddo
999 end program code426
lib/conv232.f90
-38
View File
@@ -1,38 +0,0 @@
! Layland-Lushbaugh polynomials for a K=32, r=1/2 convolutional code,
! and 8-bit parity lookup table.
data npoly1/-221228207/,npoly2/-463389625/
integer*1 partab(0:255)
data partab/ &
0, 1, 1, 0, 1, 0, 0, 1, &
1, 0, 0, 1, 0, 1, 1, 0, &
1, 0, 0, 1, 0, 1, 1, 0, &
0, 1, 1, 0, 1, 0, 0, 1, &
1, 0, 0, 1, 0, 1, 1, 0, &
0, 1, 1, 0, 1, 0, 0, 1, &
0, 1, 1, 0, 1, 0, 0, 1, &
1, 0, 0, 1, 0, 1, 1, 0, &
1, 0, 0, 1, 0, 1, 1, 0, &
0, 1, 1, 0, 1, 0, 0, 1, &
0, 1, 1, 0, 1, 0, 0, 1, &
1, 0, 0, 1, 0, 1, 1, 0, &
0, 1, 1, 0, 1, 0, 0, 1, &
1, 0, 0, 1, 0, 1, 1, 0, &
1, 0, 0, 1, 0, 1, 1, 0, &
0, 1, 1, 0, 1, 0, 0, 1, &
1, 0, 0, 1, 0, 1, 1, 0, &
0, 1, 1, 0, 1, 0, 0, 1, &
0, 1, 1, 0, 1, 0, 0, 1, &
1, 0, 0, 1, 0, 1, 1, 0, &
0, 1, 1, 0, 1, 0, 0, 1, &
1, 0, 0, 1, 0, 1, 1, 0, &
1, 0, 0, 1, 0, 1, 1, 0, &
0, 1, 1, 0, 1, 0, 0, 1, &
0, 1, 1, 0, 1, 0, 0, 1, &
1, 0, 0, 1, 0, 1, 1, 0, &
1, 0, 0, 1, 0, 1, 1, 0, &
0, 1, 1, 0, 1, 0, 0, 1, &
1, 0, 0, 1, 0, 1, 1, 0, &
0, 1, 1, 0, 1, 0, 0, 1, &
0, 1, 1, 0, 1, 0, 0, 1, &
1, 0, 0, 1, 0, 1, 1, 0/
lib/count4.f90
-54
View File
@@ -1,54 +0,0 @@
program count4
parameter(NMAX=1000)
character*47 line
real snr(NMAX)
real dt(NMAX)
real f(NMAX)
open(10,file='/users/joe/appdata/local/wsjt-x/all.txt',status='old')
read(10,1000,end=10) line
1000 format(a47)
nsync1=0
nsync2=0
n1=0
n2=0
nerr=0
do i=1,99999
read(10,1000,end=10) line
if(line(47:47).ne.' ') cycle !Skip average decodes
if(line(20:20).eq.'*') nsync1=nsync1+1
if(line(20:20).eq.'#') nsync2=nsync2+1
if(line(22:34).eq.'CQ K1ABC FN42') then
n2=n2+1 !Correlation decode
read(line,1002) snr(n2),dt(n2),f(n2)
1002 format(4x,f4.0,f5.2,f5.0)
if(line(42:42).eq.'*') n1=n1+1 !Convolutional decode
else
if(line(22:34).ne.' ') nerr=nerr+1
endif
enddo
10 call stats(snr,n2,snrave,snrdev)
call stats(dt,n2,dtave,dtdev)
call stats(f,n2,fave,fdev)
write(*,1010) nsync1,nsync2,n1,n2,nerr,snrave,dtave,fave,snrdev,dtdev,fdev
1010 format(5i5,f7.1,f7.2,f7.0/25x,f7.1,f7.2,f7.0)
end program count4
subroutine stats(x,nz,ave,rms)
real x(nz)
ave=0.
rms=0.
if(nz.gt.0) ave=sum(x)/nz
x=x-ave
if(nz.gt.1) rms=sqrt(dot_product(x,x)/(nz-1))
return
end subroutine stats
lib/deep4.f90
-169
View File
@@ -1,169 +0,0 @@
subroutine deep4(sym0,neme,flip,mycall,hiscall,hisgrid,decoded,qual)
! Deep search routine for JT4
use prog_args
parameter (MAXCALLS=7000,MAXRPT=63)
real*4 sym0(206),sym(206)
character callsign*12,grid*4,message*22,hisgrid*6,ceme*3
character*12 mycall,hiscall
character mycall0*12,hiscall0*12,hisgrid0*6
character*22 decoded
character*22 testmsg(2*MAXCALLS + 2 + MAXRPT)
character*15 callgrid(MAXCALLS)
character*180 line
character*4 rpt(MAXRPT)
integer ncode(206)
real*4 code(206,2*MAXCALLS + 2 + MAXRPT)
real pp(2*MAXCALLS + 2 + MAXRPT)
data neme0/-99/
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 mycall0,hiscall0,hisgrid0,neme0,ntot,code,testmsg
sym=sym0
if(mycall.eq.mycall0 .and. hiscall.eq.hiscall0 .and. &
hisgrid.eq.hisgrid0 .and. neme.eq.neme0) go to 30
open(23,file=trim(data_dir)//'/CALL3.TXT',status='unknown')
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. 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 encode4(message,ncode)
code(1:206,k)=2*ncode(1:206)-1
if(n.ge.2) then
! Insert CQ message
if(j4.lt.1) callgrid(icall)=callsign(1:j2)//' '//grid
message='CQ '//callgrid(icall)
k=k+1
testmsg(k)=message
call encode4(message,ncode)
code(1:206,k)=2*ncode(1:206)-1
endif
enddo
10 continue
enddo
20 continue
close(23)
ntot=k
30 mycall0=mycall
hiscall0=hiscall
hisgrid0=hisgrid
neme0=neme
sq=0.
do j=1,206
sq=sq + sym(j)**2
enddo
rms=sqrt(sq/206.0)
sym=sym/rms
p1=-1.e30
p2=-1.e30
do k=1,ntot
pp(k)=0.
! Should re-instate the following:
! if(k.ge.2 .and. k.le.64 .and. flip.gt.0.0) cycle
! 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
p=0.
do j=1,206
p=p + code(j,k)*sym(j)
enddo
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. testmsg(i).ne.testmsg(ip1)) p2=pp(i)
enddo
qual=p1-max(1.15*p2,70.0)
! ### DO NOT REMOVE ###
rewind 77
if(ip1.ge.1 .and. ip1.le.2*MAXCALLS+2+MAXRPT) write(77,1001) p1,p2,ntot, &
rms,qual,ip1,testmsg(ip1)
1001 format(2f8.2,i8,2f8.2,i6,2x,a22)
call flush(77)
if(qual.gt.1.0) then
decoded=testmsg(ip1)
else
decoded=' '
qual=0.
endif
! 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
return
end subroutine deep4
lib/entail.f90
-30
View File
@@ -1,30 +0,0 @@
subroutine entail(dgen,data0)
! Move 72-bit packed data from 6-bit to 8-bit symbols and add a zero tail.
integer dgen(13)
integer*1 data0(13)
i4=0
k=0
m=0
do i=1,12
n=dgen(i)
do j=1,6
k=k+1
i4=i4+i4+iand(1,ishft(n,j-6))
i4=iand(i4,255)
if(k.eq.8) then
m=m+1
if(i4.gt.127) i4=i4-256
data0(m)=i4
k=0
endif
enddo
enddo
do m=10,13
data0(m)=0
enddo
return
end subroutine entail
lib/ephem.f90
-85
View File
@@ -1,85 +0,0 @@
subroutine ephem(mjd0,dut,east_long,geodetic_lat,height,nspecial, &
RA,Dec,Az,El,techo,dop,fspread_1GHz,vr)
implicit real*8 (a-h,o-z)
real*8 jd !Time of observationa as a Julian Date
real*8 mjd,mjd0 !Modified Julian Date
real*8 prec(3,3) !Precession matrix, J2000 to Date
real*8 rmatn(3,3) !Nutation matrix
real*8 rme2000(6) !Vector from Earth center to Moon, JD2000
real*8 rmeDate(6) !Vector from Earth center to Moon at Date
real*8 rmeTrue(6) !Include nutation
real*8 raeTrue(6) !Vector from Earth center to Obs at Date
real*8 rmaTrue(6) !Vector from Obs to Moon at Date
logical km,bary !Set km=.true. to get km, km/s from ephemeris
common/stcomx/km,bary,pvsun(6) !Common used in JPL subroutines
common/librcom/xl(2),b(2)
twopi=8.d0*atan(1.d0) !Define some constants
rad=360.d0/twopi
clight=2.99792458d5
au2km=0.1495978706910000d9
pi=0.5d0*twopi
pio2=0.5d0*pi
km=.true.
freq=1000.0d6
do jj=1,2
mjd=mjd0
if(jj.eq.1) mjd=mjd - 1.d0/1440.d0
djutc=mjd
jd=2400000.5d0 + mjd
djtt=mjd + sla_DTT(jd)/86400.d0
ttjd=jd + sla_DTT(jd)/86400.d0
if(nspecial.ne.8) then
call pleph(ttjd,10,3,rme2000) !RME (J2000) from JPL ephemeris
year=2000.d0 + (jd-2451545.d0)/365.25d0
call sla_PREC (2000.0d0, year, prec) !Get precession matrix
rmeDate(1:3)=matmul(prec,rme2000(1:3)) !Moon geocentric xyz at Date
rmeDate(4:6)=matmul(prec,rme2000(4:6)) !Moon geocentric vel at Date
else
call sla_DMOON(djtt,rmeDate) !No JPL ephemeris, use DMOON
rmeDate=rmeDate*au2km
endif
if(nspecial.eq.7) then
rmeTrue=rmeDate
else
!Nutation to true equinox of Date
call sla_NUT(djtt,rmatn)
call sla_DMXV(rmatn,rmeDate,rmeTrue)
call sla_DMXV(rmatn,rmeDate(4),rmeTrue(4))
endif
! Local Apparent Sidereal Time:
djut1=djutc + dut/86400.d0
if(nspecial.eq.6) djut1=djutc
xlast=sla_DRANRM(sla_GMST(djut1) + sla_EQEQX(djtt) + east_long)
call sla_PVOBS(geodetic_lat,height,xlast,raeTrue)
rmaTrue=rmeTrue - raeTrue*au2km
if(nspecial.ne.2) then
! Allow for planetary aberration
tl=499.004782D0*SQRT(rmaTrue(1)**2 + rmaTrue(2)**2 + rmaTrue(3)**2)
rmaTrue(1:3)=rmaTrue(1:3)-tl*rmaTrue(4:6)/au2km
endif
!Topocentric RA, Dec, dist, velocity
call sla_DC62S(rmaTrue,RA,Dec,dist,RAdot,DECdot,vr)
dop=-2.d0 * freq * vr/clight !EME doppler shift
techo=2.d0*dist/clight !Echo delay time (s)
call libration(jd,RA,Dec,xl(jj),b(jj))
enddo
fspread_1GHz=0.0d0
dldt=57.2957795131*(xl(2)-xl(1))
dbdt=57.2957795131*(b(2)-b(1))
rate=sqrt((2*dldt)**2 + (2*dbdt)**2)
fspread_1GHz=0.5*6741*rate
call sla_DE2H(xlast-RA,Dec,geodetic_lat,Az,El)
return
end subroutine ephem
lib/fano232.f90
-138
View File
@@ -1,138 +0,0 @@
subroutine fano232(symbol,nbits,mettab,ndelta,maxcycles,dat, &
ncycles,metric,ierr)
! Sequential decoder for K=32, r=1/2 convolutional code using
! the Fano algorithm. Translated from C routine for same purpose
! written by Phil Karn, KA9Q.
parameter (MAXBITS=103)
parameter (MAXBYTES=(MAXBITS+7)/8)
integer*1 symbol(0:2*MAXBITS-1) !Soft symbols (as unsigned i*1)
integer*1 dat(MAXBYTES) !Decoded user data, 8 bits per byte
integer mettab(-128:127,0:1) !Metric table
! These were the "node" structure in Karn's C code:
integer nstate(0:MAXBITS) !Encoder state of next node
integer gamma(0:MAXBITS) !Cumulative metric to this node
integer metrics(0:3,0:MAXBITS) !Metrics indexed by all possible Tx syms
integer tm(0:1,0:MAXBITS) !Sorted metrics for current hypotheses
integer ii(0:MAXBITS) !Current branch being tested
logical noback
include 'conv232.f90' !Polynomials defined here
ntail=nbits-31
! Compute all possible branch metrics for each symbol pair.
! This is the only place we actually look at the raw input symbols
i4a=0
i4b=0
do np=0,nbits-1
j=2*np
i4a=symbol(j)
i4b=symbol(j+1)
metrics(0,np) = mettab(i4a,0) + mettab(i4b,0)
metrics(1,np) = mettab(i4a,0) + mettab(i4b,1)
metrics(2,np) = mettab(i4a,1) + mettab(i4b,0)
metrics(3,np) = mettab(i4a,1) + mettab(i4b,1)
enddo
np=0
nstate(np)=0
n=iand(nstate(np),npoly1) !Compute and sort branch metrics
n=ieor(n,ishft(n,-16)) !from the root node
lsym=partab(iand(ieor(n,ishft(n,-8)),255))
n=iand(nstate(np),npoly2)
n=ieor(n,ishft(n,-16))
lsym=lsym+lsym+partab(iand(ieor(n,ishft(n,-8)),255))
m0=metrics(lsym,np)
m1=metrics(ieor(3,lsym),np)
if(m0.gt.m1) then
tm(0,np)=m0 !0-branch has better metric
tm(1,np)=m1
else
tm(0,np)=m1 !1-branch is better
tm(1,np)=m0
nstate(np)=nstate(np) + 1 !Set low bit
endif
ii(np)=0 !Start with best branch
gamma(np)=0
nt=0
do i=1,nbits*maxcycles !Start the Fano decoder
ngamma=gamma(np) + tm(ii(np),np) !Look forward
if(ngamma.ge.nt) then
! Node is acceptable. If first time visiting this node, tighten threshold:
if(gamma(np).lt.(nt+ndelta)) nt=nt + ndelta * ((ngamma-nt)/ndelta)
gamma(np+1)=ngamma !Move forward
nstate(np+1)=ishft(nstate(np),1)
np=np+1
if(np.eq.nbits) go to 100 !We're done!
n=iand(nstate(np),npoly1)
n=ieor(n,ishft(n,-16))
lsym=partab(iand(ieor(n,ishft(n,-8)),255))
n=iand(nstate(np),npoly2)
n=ieor(n,ishft(n,-16))
lsym=lsym+lsym+partab(iand(ieor(n,ishft(n,-8)),255))
if(np.ge.ntail) then
tm(0,np)=metrics(lsym,np) !We're in the tail, now all zeros
else
m0=metrics(lsym,np)
m1=metrics(ieor(3,lsym),np)
if(m0.gt.m1) then
tm(0,np)=m0 !0-branch has better metric
tm(1,np)=m1
else
tm(0,np)=m1 !1-branch is better
tm(1,np)=m0
nstate(np)=nstate(np) + 1 !Set low bit
endif
endif
ii(np)=0 !Start with best branch
else
do while(.true.)
noback=.false. !Threshold violated, can't go forward
if(np.eq.0) noback=.true.
if(np.gt.0) then
if(gamma(np-1).lt.nt) noback=.true.
endif
if(noback) then !Can't back up, either
nt=nt-ndelta !Relax threshold and look forward again
if(ii(np).ne.0) then
ii(np)=0
nstate(np)=ieor(nstate(np),1)
endif
exit
endif
np=np-1 !Back up
if(np.lt.ntail .and. ii(np).ne.1) then
ii(np)=ii(np)+1 !Search the next best branch
nstate(np)=ieor(nstate(np),1)
exit
endif
enddo
endif
enddo
i=nbits*maxcycles
100 metric=gamma(np) !Final path metric
nbytes=(nbits+7)/8 !Copy decoded data to user's buffer
np=7
do j=1,nbytes-1
i4a=nstate(np)
dat(j)=i4a
np=np+8
enddo
dat(nbytes)=0
ncycles=i+1
ierr=0
if(i.ge.maxcycles*nbits) ierr=-1
return
end subroutine fano232
lib/fcal.f90
-117
View File
@@ -1,117 +0,0 @@
program fcal
! Compute Intercept (A) and Slope (B) for a series of FreqCal measurements.
parameter(NZ=1000)
implicit real*8 (a-h,o-z)
real*8 fd(NZ),deltaf(NZ),r(NZ)
character infile*50
character line*80
character cutc*8
nargs=iargc()
if(nargs.ne.1) then
print*,'Usage: fcal <infile>'
print*,'Example: fcal fmtave.out'
go to 999
endif
call getarg(1,infile)
open(10,file=infile,status='old',err=997)
open(12,file='fcal.out',status='unknown')
open(13,file='fcal.plt',status='unknown')
i=0
do j=1,9999
read(10,1000,end=10) line
1000 format(a80)
i0=index(line,' 0 ')
i1=index(line,' 1 ')
if(i0.le.0 .and. i1.le.0) then
read(line,*,err=5) f,df
ncal=1
i=i+1
fd(i)=f
deltaf(i)=df
else if(i1.gt.0) then
i=i+1
read(line,*,err=5) f,df,ncal,nn,rr,cutc
fd(i)=f
deltaf(i)=df
r(i)=0.d0
endif
5 continue
enddo
10 iz=i
if(iz.lt.2) go to 998
call fit(fd,deltaf,r,iz,a,b,sigmaa,sigmab,rms)
write(*,1002)
1002 format(' Freq DF Meas Freq Resid'/ &
' (MHz) (Hz) (MHz) (Hz)'/ &
'-----------------------------------------')
do i=1,iz
fm=fd(i) + 1.d-6*deltaf(i)
calfac=1.d0 + 1.d-6*deltaf(i)/fd(i)
write(*,1010) fd(i),deltaf(i),fm,r(i)
write(13,1010) fd(i),deltaf(i),fm,r(i)
1010 format(f8.3,f9.3,f14.9,f9.3,2x,a6)
enddo
calfac=1.d0 + 1.d-6*b
err=1.d-6*sigmab
if(iz.ge.3) then
write(*,1100) a,b,rms
1100 format(/'A:',f8.2,' Hz B:',f9.4,' ppm StdDev:',f7.3,' Hz')
if(iz.gt.2) write(*,1110) sigmaa,sigmab
1110 format('err:',f6.2,9x,f9.4,23x,f13.9)
else
write(*,1120) a,b
1120 format(/'A:',f8.2,' Hz B:',f9.4)
endif
write(12,1130) a,b
1130 format(f10.4)
go to 999
997 print*,'Cannot open input file: ',infile
go to 999
998 print*,'Input file must contain at least 2 valid measurement pairs'
999 end program fcal
subroutine fit(x,y,r,iz,a,b,sigmaa,sigmab,rms)
implicit real*8 (a-h,o-z)
real*8 x(iz),y(iz),r(iz)
sx=0.d0
sy=0.d0
sxy=0.d0
sx2=0.d0
do i=1,iz
sx=sx + x(i)
sy=sy + y(i)
sxy=sxy + x(i)*y(i)
sx2=sx2 + x(i)*x(i)
enddo
delta=iz*sx2 - sx*sx
a=(sx2*sy - sx*sxy)/delta
b=(iz*sxy - sx*sy)/delta
sq=0.d0
do i=1,iz
r(i)=y(i) - (a + b*x(i))
sq=sq + r(i)**2
enddo
rms=0.
sigmaa=0.
sigmab=0.
if(iz.ge.3) then
rms=sqrt(sq/(iz-2))
sigmaa=sqrt(rms*rms*sx2/delta)
sigmab=sqrt(iz*rms*rms/delta)
endif
return
end subroutine fit
lib/fchisq0.f90
-23
View File
@@ -1,23 +0,0 @@
real function fchisq0(y,npts,a)
real y(npts),a(4)
! rewind 51
chisq = 0.
do i=1,npts
x=i
z=(x-a(3))/(0.5*a(4))
yfit=a(1)
if(abs(z).lt.3.0) then
d=1.0 + z*z
yfit=a(1) + a(2) * (1.0/d - 0.1)
endif
chisq=chisq + (y(i) - yfit)**2
! write(51,3001) i,y(i),yfit,y(i)-yfit
!3001 format(i5,3f10.4)
enddo
fchisq0=chisq
return
end function fchisq0
lib/fil3.f90
-159
View File
@@ -1,159 +0,0 @@
subroutine fil3(x1,n1,c2,n2)
! FIR real-to-complex filter designed using ScopeFIR
!
!-----------------------------------------------
! fsample (Hz) 12000 Input sample rate
! Ntaps 113 Number of filter taps
! fc (Hz) 500 Cutoff frequency
! fstop (Hz) 750 Lower limit of stopband
! Ripple (dB) 0.2 Ripple in passband
! Stop Atten (dB) 50 Stopband attenuation
! fmix (HZ) 1500 Mixing frequency
! fout (Hz) 1500 Output sample rate
! Resulting passband is 1000 - 2000 Hz
! Suggest calling with n1 = 8*n2 + 105, where n2 is the desired number
! of 1500 Hz output samples.
parameter (NTAPS=113)
parameter (NH=NTAPS/2)
parameter (NDOWN=8) !Downsample ratio = 1/8
real x1(n1)
complex z
complex c2(n1/NDOWN)
! Filter coefficients:
complex ca(-NH:NH)
data ca/ &
(-0.001818142144, 0.000000000000), &
(-0.000664066641,-0.000664066640), &
(-0.000000000000,-0.001044063550), &
( 0.000737290018,-0.000737290010), &
( 0.000908957610,-0.000000000000), &
( 0.000444156615, 0.000444156610), &
(-0.000000000000, 0.000202701460), &
( 0.000244876473,-0.000244876470), &
( 0.000978154552, 0.000000000000), &
( 0.001155650277, 0.001155650270), &
( 0.000000000000, 0.002243121590), &
(-0.001927618608, 0.001927618600), &
(-0.003006201675, 0.000000000000), &
(-0.002134087852,-0.002134087850), &
( 0.000000000000,-0.002717699570), &
( 0.001478946738,-0.001478946730), &
( 0.001162489032, 0.000000000000), &
(-0.000005589545,-0.000005589540), &
(-0.000000000000,-0.001321554800), &
( 0.001873767954,-0.001873767950), &
( 0.003843608784,-0.000000000000), &
( 0.003356874940, 0.003356874940), &
(-0.000000000000, 0.005218967040), &
(-0.003640348011, 0.003640348010), &
(-0.004470167307, 0.000000000000), &
(-0.002247131477,-0.002247131470), &
(-0.000000000000,-0.001335998900), &
(-0.000647656208, 0.000647656200), &
(-0.003386100636, 0.000000000000), &
(-0.004114456189,-0.004114456180), &
( 0.000000000000,-0.007939147960), &
( 0.006692816134,-0.006692816130), &
( 0.010145641899, 0.000000000000), &
( 0.006920770724, 0.006920770720), &
( 0.000000000000, 0.008285915750), &
(-0.003992321524, 0.003992321520), &
(-0.001995842303, 0.000000000000), &
( 0.001704388774, 0.001704388770), &
(-0.000000000000, 0.007202515550), &
(-0.008426458377, 0.008426458370), &
(-0.016028350845, 0.000000000000), &
(-0.013430355885,-0.013430355880), &
(-0.000000000000,-0.020297455950), &
( 0.013791263729,-0.013791263720), &
( 0.016298136197,-0.000000000000), &
( 0.007443596155, 0.007443596150), &
(-0.000000000000, 0.002223837360), &
( 0.005924356866,-0.005924356860), &
( 0.020854478160, 0.000000000000), &
( 0.024471928130, 0.024471928130), &
( 0.000000000000, 0.048909701460), &
(-0.044508219241, 0.044508219240), &
(-0.075874892030, 0.000000000000), &
(-0.061450241075,-0.061450241070), &
( 0.000000000000,-0.095332017640), &
( 0.071148679982,-0.071148679980), &
( 0.102420526192, 0.000000000000), &
( 0.071148679982, 0.071148679980), &
( 0.000000000000, 0.095332017640), &
(-0.061450241075, 0.061450241070), &
(-0.075874892030, 0.000000000000), &
(-0.044508219241,-0.044508219240), &
( 0.000000000000,-0.048909701460), &
( 0.024471928130,-0.024471928130), &
( 0.020854478160, 0.000000000000), &
( 0.005924356866, 0.005924356860), &
(-0.000000000000,-0.002223837360), &
( 0.007443596155,-0.007443596150), &
( 0.016298136197,-0.000000000000), &
( 0.013791263729, 0.013791263720), &
(-0.000000000000, 0.020297455950), &
(-0.013430355885, 0.013430355880), &
(-0.016028350845, 0.000000000000), &
(-0.008426458377,-0.008426458370), &
(-0.000000000000,-0.007202515550), &
( 0.001704388774,-0.001704388770), &
(-0.001995842303, 0.000000000000), &
(-0.003992321524,-0.003992321520), &
( 0.000000000000,-0.008285915750), &
( 0.006920770724,-0.006920770720), &
( 0.010145641899, 0.000000000000), &
( 0.006692816134, 0.006692816130), &
( 0.000000000000, 0.007939147960), &
(-0.004114456189, 0.004114456180), &
(-0.003386100636, 0.000000000000), &
(-0.000647656208,-0.000647656200), &
(-0.000000000000, 0.001335998900), &
(-0.002247131477, 0.002247131470), &
(-0.004470167307, 0.000000000000), &
(-0.003640348011,-0.003640348010), &
(-0.000000000000,-0.005218967040), &
( 0.003356874940,-0.003356874940), &
( 0.003843608784,-0.000000000000), &
( 0.001873767954, 0.001873767950), &
(-0.000000000000, 0.001321554800), &
(-0.000005589545, 0.000005589540), &
( 0.001162489032, 0.000000000000), &
( 0.001478946738, 0.001478946730), &
( 0.000000000000, 0.002717699570), &
(-0.002134087852, 0.002134087850), &
(-0.003006201675, 0.000000000000), &
(-0.001927618608,-0.001927618600), &
( 0.000000000000,-0.002243121590), &
( 0.001155650277,-0.001155650270), &
( 0.000978154552, 0.000000000000), &
( 0.000244876473, 0.000244876470), &
(-0.000000000000,-0.000202701460), &
( 0.000444156615,-0.000444156610), &
( 0.000908957610,-0.000000000000), &
( 0.000737290018, 0.000737290010), &
(-0.000000000000, 0.001044063550), &
(-0.000664066641, 0.000664066640), &
(-0.001818142144, 0.000000000000)/
save ca
n2=(n1-NTAPS+NDOWN)/NDOWN
k0=NH-NDOWN+1
! Loop over all output samples
do i=1,n2
z=0.
k=k0 + NDOWN*i
do j=-NH,NH
z=z + x1(j+k)*ca(j)
enddo
c2(i)=z
enddo
return
end subroutine fil3
lib/fil3c.f90
-72
View File
@@ -1,72 +0,0 @@
subroutine fil3c(c1,n1,c2,n2)
! FIR complex-to-complex low-pass filter designed with ScopeFIR
!
!-----------------------------------------------
! fsample (Hz) 12000 Input sample rate
! Ntaps 113 Number of filter taps
! fc (Hz) 500 Cutoff frequency
! fstop (Hz) 750 Lower limit of stopband
! Ripple (dB) 0.2 Ripple in passband
! Stop Atten (dB) 50 Stopband attenuation
! fout (Hz) 1500 Output sample rate
! Suggest calling with n1 = 8*n2 + 105, where n2 is the desired number
! of 1500 Hz output samples.
parameter (NTAPS=113)
parameter (NH=NTAPS/2)
parameter (NDOWN=8) !Downsample ratio = 1/8
complex c1(n1)
complex c2(n1/NDOWN)
complex z
! Filter coefficients:
real a(-NH:NH)
data a/ &
-0.001818142144,-0.000939132050,-0.001044063556,-0.001042685542, &
-0.000908957610,-0.000628132309,-0.000202701465, 0.000346307629, &
0.000978154552, 0.001634336295, 0.002243121592, 0.002726064379, &
0.003006201675, 0.003018055983, 0.002717699575, 0.002091546534, &
0.001162489032,-0.000007904811,-0.001321554806,-0.002649908053, &
-0.003843608784,-0.004747338068,-0.005218967042,-0.005148229529, &
-0.004470167307,-0.003177923811,-0.001335998901, 0.000915924193, &
0.003386100636, 0.005818719744, 0.007939147967, 0.009465071347, &
0.010145641899, 0.009787447819, 0.008285915754, 0.005645995244, &
0.001995842303,-0.002410369720,-0.007202515555,-0.011916811719, &
-0.016028350845,-0.018993391440,-0.020297455955,-0.019503792208, &
-0.016298136197,-0.010526834635,-0.002223837363, 0.008378305829, &
0.020854478160, 0.034608532659, 0.048909701463, 0.062944127288, &
0.075874892030, 0.086903764340, 0.095332017649, 0.100619428175, &
0.102420526192, 0.100619428175, 0.095332017649, 0.086903764340, &
0.075874892030, 0.062944127288, 0.048909701463, 0.034608532659, &
0.020854478160, 0.008378305829,-0.002223837363,-0.010526834635, &
-0.016298136197,-0.019503792208,-0.020297455955,-0.018993391440, &
-0.016028350845,-0.011916811719,-0.007202515555,-0.002410369720, &
0.001995842303, 0.005645995244, 0.008285915754, 0.009787447819, &
0.010145641899, 0.009465071347, 0.007939147967, 0.005818719744, &
0.003386100636, 0.000915924193,-0.001335998901,-0.003177923811, &
-0.004470167307,-0.005148229529,-0.005218967042,-0.004747338068, &
-0.003843608784,-0.002649908053,-0.001321554806,-0.000007904811, &
0.001162489032, 0.002091546534, 0.002717699575, 0.003018055983, &
0.003006201675, 0.002726064379, 0.002243121592, 0.001634336295, &
0.000978154552, 0.000346307629,-0.000202701465,-0.000628132309, &
-0.000908957610,-0.001042685542,-0.001044063556,-0.000939132050, &
-0.001818142144/
save a
n2=(n1-NTAPS+NDOWN)/NDOWN
k0=NH-NDOWN+1
! Loop over all output samples
do i=1,n2
z=0.
k=k0 + NDOWN*i
do j=-NH,NH
z=z + c1(j+k)*a(j)
enddo
c2(i)=z
enddo
return
end subroutine fil3c
lib/fil61.f90
-64
View File
@@ -1,64 +0,0 @@
12000 61 250 750 0.2 50, mix at 1500
-0.000000000000 0.001944450121
-0.000668730681 0.000668730681
-0.000974850191 -0.000000000000
-0.000581679123 -0.000581679123
0.000000000000 -0.000439648787
-0.000148911451 0.000148911451
-0.001140891736 -0.000000000000
-0.001653102965 -0.001653102965
0.000000000000 -0.003749915818
0.003740834397 -0.003740834397
0.006834087490 0.000000000000
0.005812808655 0.005812808655
-0.000000000000 0.009262713933
-0.006900370427 0.006900370427
-0.009503248519 -0.000000000000
-0.005874581677 -0.005874581677
0.000000000000 -0.006017530719
0.001785268072 -0.001785268072
-0.002214736448 -0.000000000000
-0.005777038427 -0.005777038427
0.000000000000 -0.015228682747
0.016402831440 -0.016402831440
0.031806920774 0.000000000000
0.028800401613 0.028800401613
-0.000000000000 0.049589395998
-0.041000303659 0.041000303659
-0.065514139214 -0.000000000000
-0.050781544715 -0.050781544715
0.000000000000 -0.076562341482
0.056225821996 -0.056225821996
0.080516569816 0.000000000000
0.056225821996 0.056225821996
-0.000000000000 0.076562341482
-0.050781544715 0.050781544715
-0.065514139214 -0.000000000000
-0.041000303659 -0.041000303659
0.000000000000 -0.049589395998
0.028800401613 -0.028800401613
0.031806920774 0.000000000000
0.016402831440 0.016402831440
-0.000000000000 0.015228682747
-0.005777038427 0.005777038427
-0.002214736448 -0.000000000000
0.001785268072 0.001785268072
-0.000000000000 0.006017530719
-0.005874581677 0.005874581677
-0.009503248519 -0.000000000000
-0.006900370427 -0.006900370427
0.000000000000 -0.009262713933
0.005812808655 -0.005812808655
0.006834087490 0.000000000000
0.003740834397 0.003740834397
-0.000000000000 0.003749915818
-0.001653102965 0.001653102965
-0.001140891736 -0.000000000000
-0.000148911451 -0.000148911451
-0.000000000000 0.000439648787
-0.000581679123 0.000581679123
-0.000974850191 -0.000000000000
-0.000668730681 -0.000668730681
0.000000000000 -0.001944450121
lib/fitcal.f90
-34
View File
@@ -1,34 +0,0 @@
subroutine fitcal(x,y,r,iz,a,b,sigmaa,sigmab,rms)
implicit real*8 (a-h,o-z)
real*8 x(iz),y(iz),r(iz)
sx=0.d0
sy=0.d0
sxy=0.d0
sx2=0.d0
do i=1,iz
sx=sx + x(i)
sy=sy + y(i)
sxy=sxy + x(i)*y(i)
sx2=sx2 + x(i)*x(i)
enddo
delta=iz*sx2 - sx*sx
a=(sx2*sy - sx*sxy)/delta
b=(iz*sxy - sx*sy)/delta
sq=0.d0
do i=1,iz
r(i)=y(i) - (a + b*x(i))
sq=sq + r(i)**2
enddo
rms=0.
sigmaa=0.
sigmab=0.
if(iz.ge.3) then
rms=sqrt(sq/(iz-2))
sigmaa=sqrt(rms*rms*sx2/delta)
sigmab=sqrt(iz*rms*rms/delta)
endif
return
end subroutine fitcal
lib/flat2.f90
-18
View File
@@ -1,18 +0,0 @@
subroutine flat2(s,nz,ref)
parameter (NSMAX=6827)
real s(NSMAX)
real ref(NSMAX)
nsmo=10
ia=nsmo+1
ib=nz-nsmo-1
do i=ia,ib
call pctile(s(i-nsmo),2*nsmo+1,5,ref(i))
enddo
ref(:ia-1)=ref(ia)
ref(ib+1:)=ref(ib)
return
end subroutine flat2
lib/fmeasure.f90
-75
View File
@@ -1,75 +0,0 @@
!-------------------------------------------------------------------------------
!
! This file is part of the WSPR application, Weak Signal Propagation Reporter
!
! File Name: fmeasure.f90
! Description:
!
! Copyright (C) 2001-2014 Joseph Taylor, K1JT
! License: GPL-3
!
! This program is free software; you can redistribute it and/or modify it under
! the terms of the GNU General Public License as published by the Free Software
! Foundation; either version 3 of the License, or (at your option) any later
! version.
!
! This program is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
! details.
!
! You should have received a copy of the GNU General Public License along with
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
! Street, Fifth Floor, Boston, MA 02110-1301, USA.
!
!-------------------------------------------------------------------------------
program fmeasure
parameter(NZ=1000)
implicit real*8 (a-h,o-z)
character infile*50
character line*80
nargs=iargc()
if(nargs.ne.1) then
print*,'Usage: fmeasure <infile>'
print*,'Example: fmeasure fmtave.out'
go to 999
endif
call getarg(1,infile)
open(10,file=infile,status='old',err=997)
open(11,file='fcal.out',status='old',err=998)
open(12,file='fmeasure.out',status='unknown')
read(11,*) a,b
write(*,1000)
write(12,1000)
1000 format(' Freq DF A+B*f Corrected Offset'/ &
' (MHz) (Hz) (Hz) (MHz) (Hz)'/ &
'-----------------------------------------------')
i=0
do j=1,9999
read(10,1010,end=999) line
1010 format(a80)
i0=index(line,' 0 ')
if(i0.gt.0) then
read(line,*,err=5) f,df
dial_error=a + b*f
fcor=f + 1.d-6*df - 1.d-6*dial_error
offset_hz=1.d6*(fcor-f)
write(*,1020) f,df,dial_error,fcor,offset_hz
write(12,1020) f,df,dial_error,fcor,offset_hz
1020 format(3f8.3,f15.9,f8.2)
endif
5 continue
enddo
go to 999
997 print*,'Cannot open input file: ',infile
go to 999
998 print*,'Cannot open fcal.out'
999 end program fmeasure
lib/fmtave.f90
-64
View File
@@ -1,64 +0,0 @@
program fmtave
! Average groups of frequency-calibration measurements.
implicit real*8 (a-h,o-z)
character infile*80
character*8 cutc,cutc1
nargs=iargc()
if(nargs.ne.1) then
print*,'Usage: fmtave <infile>'
print*,'Example: fmtave fmt.all'
go to 999
endif
call getarg(1,infile)
open(10,file=infile,status='old')
open(12,file='fmtave.out',status='unknown')
write(*,1000)
1000 format(' Freq DF CAL N rms UTC Call'/ &
' (kHz) (Hz) ? (Hz)'/ &
'----------------------------------------------------')
nkhz0=0
sum=0.d0
sumsq=0.d0
n=0
do i=1,99999
read(10,*,end=10) cutc,nkHz,ncal,noffset,faudio,df,dblevel,snr
if((nkHz.ne.nkHz0) .and. i.ne.1) then
ave=sum/n
rms=0.d0
if(n.gt.1) then
rms=sqrt(abs(sumsq - sum*sum/n)/(n-1.d0))
endif
fMHz=0.001d0*nkHz0
write(*,1010) fMHz,ave,ncal0,n,rms,cutc1
write(12,1010) fMHz,ave,ncal0,n,rms,cutc1
1010 format(f8.3,f9.3,i4,i5,f8.2,2x,a8,2x,a6)
sum=0.d0
sumsq=0.d0
n=0
endif
dial_error=faudio-noffset
sum=sum + dial_error
sumsq=sumsq + dial_error**2
n=n+1
if(n.eq.1) then
cutc1=cutc
ncal0=ncal
endif
nkHz0=nkHz
enddo
10 ave=sum/n
rms=0.d0
if(n.gt.0) then
rms=sqrt((sumsq - sum*sum/n)/(n-1.d0))
endif
fMHz=0.001d0*nkHz
write(*,1010) fMHz,ave,ncal,n,rms,cutc1
write(12,1010) fMHz,ave,ncal,n,rms,cutc1
999 end program fmtave
lib/fqso_first.f90
-32
View File
@@ -1,32 +0,0 @@
subroutine fqso_first(nfqso,ntol,ca,ncand)
! If a candidate was found within +/- ntol of nfqso, move it into ca(1).
type candidate
real freq
real dt
real sync
real flip
end type candidate
type(candidate) ca(300),cb
dmin=1.e30
i0=0
do i=1,ncand
d=abs(ca(i)%freq-nfqso)
if(d.lt.dmin) then
i0=i
dmin=d
endif
enddo
if(dmin.lt.float(ntol)) then
cb=ca(i0)
do i=i0,2,-1
ca(i)=ca(i-1)
enddo
ca(1)=cb
endif
return
end subroutine fqso_first
lib/freqcal.f90
-84
View File
@@ -1,84 +0,0 @@
subroutine freqcal(id2,k,nkhz,noffset,ntol,line)
parameter (NZ=30*12000,NFFT=55296,NH=NFFT/2)
integer*2 id2(0:NZ-1)
complex sp,sn
real x(0:NFFT-1)
real xi(0:NFFT-1)
real w(0:NFFT-1) !Window function
real s(0:NH)
character line*80,cflag*1
logical first
complex cx(0:NH)
equivalence (x,cx)
data n/0/,k0/9999999/,first/.true./
save n,k0,w,first,pi,fs,xi
if(first) then
pi=4.0*atan(1.0)
fs=12000.0
do i=0,NFFT-1
ww=sin(i*pi/NFFT)
w(i)=ww*ww/NFFT
xi(i)=2.0*pi*i
enddo
first=.false.
endif
if(k.lt.NFFT) go to 900
if(k.lt.k0) n=0
k0=k
x=w*id2(k-NFFT:k-1) !Apply window
call four2a(x,NFFT,1,-1,0) !Compute spectrum, r2c
df=fs/NFFT
if (ntol.gt.noffset) then
ia=0
ib=nint((noffset*2)/df)
else
ia=nint((noffset-ntol)/df)
ib=nint((noffset+ntol)/df)
endif
smax=0.
s=0.
do i=ia,ib
s(i)=real(cx(i))**2 + aimag(cx(i))**2
if(s(i).gt.smax) then
smax=s(i)
ipk=i
endif
enddo
call peakup(s(ipk-1),s(ipk),s(ipk+1),dx)
fpeak=df * (ipk+dx)
ap=(fpeak/fs+1.0/(2.0*NFFT))
an=(fpeak/fs-1.0/(2.0*NFFT))
sp=sum(id2((k-NFFT):k-1)*cmplx(cos(xi*ap),-sin(xi*ap)))
sn=sum(id2((k-NFFT):k-1)*cmplx(cos(xi*an),-sin(xi*an)))
fpeak=fpeak+fs*(abs(sp)-abs(sn))/(abs(sp)+abs(sn))/(2*NFFT)
xsum=0.
nsum=0
do i=ia,ib
if(abs(i-ipk).gt.10) then
xsum=xsum+s(i)
nsum=nsum+1
endif
enddo
ave=xsum/nsum
snr=db(smax/ave)
pave=db(ave) + 8.0
cflag=' '
if(snr.lt.20.0) cflag='*'
n=n+1
nsec=mod(time(),86400)
nhr=nsec/3600
nmin=mod(nsec/60,60)
nsec=mod(nsec,60)
ncal=1
ferr=fpeak-noffset
write(line,1100) nhr,nmin,nsec,nkhz,ncal,noffset,fpeak,ferr,pave, &
snr,cflag,char(0)
1100 format(i2.2,':',i2.2,':',i2.2,i7,i3,i6,2f10.3,2f7.1,2x,a1,a1)
900 return
end subroutine freqcal
lib/getlags.f90
-27
View File
@@ -1,27 +0,0 @@
subroutine getlags(nsps8,lag0,lag1,lag2)
if(nsps8.eq.864) then
lag1=39
lag2=291
lag0=123
else if(nsps8.eq.1920) then
lag1=70
lag2=184
lag0=108
else if(nsps8.eq.5120) then
lag1=84
lag2=129
lag0=99
else if(nsps8.eq.10368) then
lag1=91
lag2=112
lag0=98
else if(nsps8.eq.31500) then
lag1=93
lag2=102
lag0=96
else
stop 'Error in getlags'
endif
return
end subroutine getlags
lib/getmet4.f90
-56
View File
@@ -1,56 +0,0 @@
subroutine getmet4(mettab,ndelta)
! Return appropriate metric table for soft-decision convolutional decoder.
! Metric table (RxSymbol,TxSymbol)
! integer mettab(0:255,0:1)
integer mettab(-128:127,0:1)
real*4 xx0(0:255)
data xx0/ &
1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, &
1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, &
1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, &
1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, &
1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, &
1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, &
0.988, 1.000, 0.991, 0.993, 1.000, 0.995, 1.000, 0.991, &
1.000, 0.991, 0.992, 0.991, 0.990, 0.990, 0.992, 0.996, &
0.990, 0.994, 0.993, 0.991, 0.992, 0.989, 0.991, 0.987, &
0.985, 0.989, 0.984, 0.983, 0.979, 0.977, 0.971, 0.975, &
0.974, 0.970, 0.970, 0.970, 0.967, 0.962, 0.960, 0.957, &
0.956, 0.953, 0.942, 0.946, 0.937, 0.933, 0.929, 0.920, &
0.917, 0.911, 0.903, 0.895, 0.884, 0.877, 0.869, 0.858, &
0.846, 0.834, 0.821, 0.806, 0.790, 0.775, 0.755, 0.737, &
0.713, 0.691, 0.667, 0.640, 0.612, 0.581, 0.548, 0.510, &
0.472, 0.425, 0.378, 0.328, 0.274, 0.212, 0.146, 0.075, &
0.000,-0.079,-0.163,-0.249,-0.338,-0.425,-0.514,-0.606, &
-0.706,-0.796,-0.895,-0.987,-1.084,-1.181,-1.280,-1.376, &
-1.473,-1.587,-1.678,-1.790,-1.882,-1.992,-2.096,-2.201, &
-2.301,-2.411,-2.531,-2.608,-2.690,-2.829,-2.939,-3.058, &
-3.164,-3.212,-3.377,-3.463,-3.550,-3.768,-3.677,-3.975, &
-4.062,-4.098,-4.186,-4.261,-4.472,-4.621,-4.623,-4.608, &
-4.822,-4.870,-4.652,-4.954,-5.108,-5.377,-5.544,-5.995, &
-5.632,-5.826,-6.304,-6.002,-6.559,-6.369,-6.658,-7.016, &
-6.184,-7.332,-6.534,-6.152,-6.113,-6.288,-6.426,-6.313, &
-9.966,-6.371,-9.966,-7.055,-9.966,-6.629,-6.313,-9.966, &
-5.858,-9.966,-9.966,-9.966,-9.966,-9.966,-9.966,-9.966, &
-9.966,-9.966,-9.966,-9.966,-9.966,-9.966,-9.966,-9.966, &
-9.966,-9.966,-9.966,-9.966,-9.966,-9.966,-9.966,-9.966, &
-9.966,-9.966,-9.966,-9.966,-9.966,-9.966,-9.966,-9.966, &
-9.966,-9.966,-9.966,-9.966,-9.966,-9.966,-9.966,-9.966, &
-9.966,-9.966,-9.966,-9.966,-9.966,-9.966,-9.966,-9.966/
save
bias=0.5
scale=50
ndelta=nint(3.4*scale)
do i=0,255
xx=xx0(i)
if(i.ge.160) xx=xx0(160) - (i-160)*6.822/65.3
mettab(i-128,0)=nint(scale*(xx-bias))
if(i.ge.1) mettab(128-i,1)=mettab(i-128,0)
enddo
mettab(-128,1)=mettab(-127,1)
return
end subroutine getmet4
lib/go_chkfft.cmd
-1
View File
@@ -1 +0,0 @@
gfortran -o chkfft chkfft3.f90 four2a.f90 gran.c /jtsdk/fftw3f/libfftw3f-3.dll
lib/golay24_table.c
-28
View File
@@ -1,28 +0,0 @@
void golay24_table_(int cw[])
{
// Compute arnd return the table of 4096 codewords for the Golay (24,12) code.
// Array y contains the twelve rows (columns) of the parity-check matrix
int y[12] = { 0x7ff, 0xee2, 0xdc5, 0xb8b, 0xf16, 0xe2d,
0xc5b, 0x8b7, 0x96e, 0xadc, 0xdb8, 0xb71 };
int c[2]; /* Codeword composed of 12-bit info and 12-bit parity */
int i,j,k;
int aux;
int weight(int vector);
for(k=0; k<4096; k++) {
c[0] = k;
c[1] = 0;
for (i=0; i<12; i++) {
aux = 0;
for (j=0; j<12; j++) {
aux = aux ^ ((c[0] & y[i]) >> j & 1);
}
c[1] = (c[1] << 1) ^ aux;
}
cw[k]=4096*c[0] + c[1];
}
}
lib/graycode.f90
-9
View File
@@ -1,9 +0,0 @@
subroutine graycode(ia,n,idir,ib)
integer ia(n),ib(n)
do i=1,n
ib(i)=igray(ia(i),idir)
enddo
return
end subroutine graycode
lib/igray.c
-18
View File
@@ -1,18 +0,0 @@
int igray_(int *n0, int *idir)
{
int n;
unsigned long sh;
unsigned long nn;
n=*n0;
if(*idir>0) return (n ^ (n >> 1));
sh = 1;
nn = (n >> sh);
while (nn > 0) {
n ^= nn;
sh <<= 1;
nn = (n >> sh);
}
return (n);
}
lib/image.f90
-336
View File
@@ -1,336 +0,0 @@
subroutine imopen(plotfile)
character*(*) plotfile
common/imcom/ lu,npage
lu=80
open(lu,file=plotfile,status='unknown')
write(lu,1000)
1000 format('%!PS-Adobe-2.0'/ &
'/rightshow { dup stringwidth pop neg 0 rmoveto show } def'/ &
'/centershow { dup stringwidth pop neg 2 div ', &
'0 rmoveto show } def'/ &
'/lt { lineto } def'/'%%Page: 1 1')
npage=1
return
end subroutine imopen
subroutine impalette(palette)
character*(*) palette
integer r(0:8),g(0:8),b(0:8)
integer rr,gg,bb
common/imcom/ lu,npage
common/imcom2/rr(0:255),gg(0:255),bb(0:255)
if(palette.eq.'afmhot') then
do i=0,255
j=255-i
rr(i)=min(255,2*j)
gg(i)=max(0,min(255,2*j-128))
bb(i)=max(0,min(255,2*j-256))
enddo
else if(palette.eq.'hot') then
do i=0,255
j=255-i
rr(i)=min(255,3*j)
gg(i)=max(0,min(255,3*j-256))
bb(i)=max(0,min(255,3*j-512))
enddo
else
open(11,file="Palettes/"//palette,status="old")
do j=0,8
read(11,*) r(j),g(j),b(j)
enddo
close(11)
do i=0,255
j0=i/32
j1=j0+1
k=i-32*j0
rr(i)=r(j0) + int((k*(r(j1)-r(j0)))/31 + 0.5)
gg(i)=g(j0) + int((k*(g(j1)-g(j0)))/31 + 0.5)
bb(i)=b(j0) + int((k*(b(j1)-b(j0)))/31 + 0.5)
enddo
endif
return
end subroutine impalette
subroutine imclose
common/imcom/ lu,npage
write(lu,1000)
1000 format('showpage'/'%%Trailer')
close(lu)
return
end subroutine imclose
subroutine imnewpage
common/imcom/ lu,npage
npage=npage+1
write(lu,1000) npage,npage
1000 format('showpage'/'%%Page:',2i4)
return
end subroutine imnewpage
subroutine imxline(x,y,dx)
! Draw a line from (x,y) to (x+dx,y) integer r,g,b
common/imcom/ lu,npage
write(lu,1000) 72.0*x,72.0*y,72.0*dx
1000 format('newpath',2f7.1,' moveto',f7.1,' 0 rlineto stroke')
return
end subroutine imxline
subroutine imyline(x,y,dy)
! Draw a line from (x,y) to (x,y+dy)
common/imcom/ lu,npage
write(lu,1000) 72.0*x,72.0*y,72.0*dy
1000 format('newpath',2f7.1,' moveto 0',f7.1,' rlineto stroke')
return
end subroutine imyline
subroutine imwidth(width)
common/imcom/ lu,npage
write(lu,1000) width
1000 format(f7.1,' setlinewidth')
return
end subroutine imwidth
subroutine imfont(fontname,npoints)
character*(*) fontname
common/imcom/ lu,npage
write(lu,1000) fontname,npoints
1000 format('/',a,' findfont',i4,' scalefont setfont')
return
end subroutine imfont
subroutine imstring(string,x,y,just,ndeg)
character*(*) string
common/imcom/ lu,npage
write(lu,1000) 72.0*x,72.0*y,ndeg,string
1000 format(2f7.1,' moveto',i4,' rotate'/'(',a,')')
if(just.eq.1) write(lu,*) 'rightshow'
if(just.eq.2) write(lu,*) 'centershow'
if(just.eq.3) write(lu,*) 'show'
write(lu,1010) -ndeg
1010 format(i4,' rotate'/)
return
end subroutine imstring
subroutine imr4mat(z,IP,JP,imax,jmax,zz1,zz2,x,y,dx,dy,nbox)
real z(IP,JP)
integer idat(2048)
common/imcom/ lu,npage
z1=zz1
z2=zz2
if(z1.eq.0.0 .and. z2.eq.0.0) then
z1=z(1,1)
z2=z1
do i=1,imax
do j=1,jmax
z1=min(z(i,j),z1)
z2=max(z(i,j),z2)
enddo
enddo
endif
scale=255.99/(z2-z1)
write(lu,1002) 72.0*x,72.0*y,72.0*dx,72.0*dy
1002 format(2f7.1,' translate',2f7.1,' scale')
write(lu,*) imax,jmax,8,' [',imax,0,0,jmax,0,0,']'
write(lu,*) '{<'
do j=1,jmax
do i=1,imax
idat(i)=scale*(z(i,j)-z1)
idat(i)=max(idat(i),0)
idat(i)=min(idat(i),255)
idat(i)=255-idat(i)
enddo
write(lu,1004) (idat(i),i=1,imax)
1004 format(30z2.2)
enddo
write(lu,*) '>} image'
write(lu,1006) 1.0/(72.0*dx),1.0/(72.0*dy),-72.0*x,-72.0*y
1006 format(2f9.6,' scale',2f7.1,' translate')
if(nbox.ne.0) then
write(lu,1010) 72.0*x,72.0*y,72.0*dx,72.0*dy,-72*dx
1010 format('newpath',2f7.1,' moveto',f7.1,' 0 rlineto 0', &
f7.1,' rlineto',f7.1,' 0 rlineto closepath stroke')
endif
return
end subroutine imr4mat
subroutine imr4mat_color(z,IP,JP,imax,jmax,zz1,zz2,x,y,dx,dy,nbox)
real z(IP,JP)
integer idat(2048,3)
integer rr,gg,bb
common/imcom/ lu,npage
common/imcom2/rr(0:255),gg(0:255),bb(0:255)
z1=zz1
z2=zz2
if(z1.eq.0.0 .and. z2.eq.0.0) then
z1=z(1,1)
z2=z1
do i=1,imax
do j=1,jmax
z1=min(z(i,j),z1)
z2=max(z(i,j),z2)
enddo
enddo
endif
scale=255.99/(z2-z1)
write(lu,1002) 72.0*x,72.0*y,72.0*dx,72.0*dy
1002 format(2f7.1,' translate',2f7.1,' scale')
write(lu,1003) imax,jmax,8,imax,0,0,jmax,0,0
1003 format(3i5,' [',6i4,']')
write(lu,1004) imax
1004 format('{currentfile 3',i4,' mul string readhexstring pop} bind'/ &
'false 3 colorimage')
do j=1,jmax
do i=1,imax
n=scale*(z(i,j)-z1)
n=max(n,0)
n=min(n,255)
idat(i,1)=rr(n)
idat(i,2)=gg(n)
idat(i,3)=bb(n)
enddo
write(lu,1005) (idat(i,1),idat(i,2),idat(i,3),i=1,imax)
1005 format(30z2.2)
enddo
write(lu,1006) 1.0/(72.0*dx),1.0/(72.0*dy),-72.0*x,-72.0*y
1006 format(2f9.6,' scale',2f7.1,' translate')
if(nbox.ne.0) then
write(lu,1010) 72.0*x,72.0*y,72.0*dx,72.0*dy,-72*dx
1010 format('newpath',2f7.1,' moveto',f7.1,' 0 rlineto 0', &
f7.1,' rlineto',f7.1,' 0 rlineto closepath stroke')
endif
return
end subroutine imr4mat_color
subroutine imr4pro(p,imax,yy1,yy2,x,y,dx,dy,nbox)
real p(imax)
common/imcom/ lu,npage
y1=yy1
y2=yy2
if(y1.eq.0.0 .and. y2.eq.0.0) then
y1=p(1)
y2=y1
do i=1,imax
y1=min(p(i),y1)
y2=max(p(i),y2)
enddo
endif
xscale=72.0*dx/imax
xoff=72.0*x
yscale=72.0*dy
if(y1.ne.y2) yscale=yscale/(y2-y1)
yoff=72.0*y
write(lu,*) '1.416 setmiterlimit'
write(lu,1002) xoff+0.5*xscale,yoff+yscale*(p(1)-y1)
1002 format('newpath',2f7.1,' moveto')
do i=2,imax
write(lu,1004) xoff+(i-0.5)*xscale,yoff+yscale*(p(i)-y1)
1004 format(2f6.1,' lt')
enddo
write(lu,*) 'stroke'
if(nbox.ne.0) then
write(lu,1010) xoff,yoff,72.0*dx,72.0*dy,-72*dx
1010 format('newpath',2f7.1,' moveto',f7.1,' 0 rlineto 0', &
f7.1,' rlineto',f7.1,' 0 rlineto closepath stroke')
endif
return
end subroutine imr4pro
subroutine imline(x1,y1,x2,y2)
common/imcom/ lu,npage
write(lu,1000) 72*x1,72*y1,72*x2,72*y2
1000 format('newpath',2f7.1,' moveto',2f7.1,' lineto stroke')
return
end subroutine imline
subroutine imcircle(x,y,radius,shade)
common/imcom/ lu,npage
write(lu,1000) shade
1000 format(f7.1,' setgray')
write(lu,1002) 72*x,72*y,72*radius
1002 format('newpath',3f7.1,' 0 360 arc fill')
write(lu,1000) 0.0
write(lu,1004) 72*x,72*y,72*radius
1004 format('newpath',3f7.1,' 0 360 arc stroke')
return
end subroutine imcircle
subroutine imtriangle(x,y,rr,shade)
common/imcom/ lu,npage
write(lu,1000) shade
1000 format(f7.1,' setgray')
write(lu,1002) 72*x,72*(y+rr)
1002 format('newpath',2f7.1,' moveto ')
write(lu,1004) 72*(x-rr),72*(y-rr)
1004 format(2f7.1,' lineto ')
write(lu,1004) 72*(x+rr),72*(y-rr)
write(lu,*) 'closepath fill 0 setgray'
write(lu,1002) 72*x,72*(y+rr)
write(lu,1004) 72*(x-rr),72*(y-rr)
write(lu,1004) 72*(x+rr),72*(y-rr)
write(lu,*) 'closepath stroke'
return
end subroutine imtriangle
subroutine imr4prov(p,jmax,xx1,xx2,x,y,dx,dy,nbox)
real p(jmax)
common/imcom/ lu,npage
x1=xx1
x2=xx2
if(x1.eq.0.0 .and. x2.eq.0.0) then
x1=p(1)
x2=x1
do j=1,jmax
x1=min(p(j),x1)
x2=max(p(j),x2)
enddo
endif
xscale=72.0*dx
xoff=72.0*x
if(x1.ne.x2) xscale=xscale/(x2-x1)
yscale=72.0*dy/jmax
yoff=72.0*y
write(lu,*) '1.416 setmiterlimit'
write(lu,1002) xoff+xscale*(x2-p(1)),yoff+0.5*yscale
1002 format('newpath',2f7.1,' moveto')
do j=2,jmax
write(lu,1004) xoff+xscale*(x2-p(j)),yoff+(j-0.5)*yscale
1004 format(2f6.1,' lt')
enddo
write(lu,*) 'stroke'
if(nbox.ne.0) then
write(lu,1010) xoff,yoff,72.0*dx,72.0*dy,-72*dx
1010 format('newpath',2f7.1,' moveto',f7.1,' 0 rlineto 0', &
f7.1,' rlineto',f7.1,' 0 rlineto closepath stroke')
endif
return
end subroutine imr4prov
lib/jplsubs.f
-899
View File
@@ -1,899 +0,0 @@
C++++++++++++++++++++++++
C
SUBROUTINE FSIZER1(NRECL,KSIZE,NRFILE,NAMFIL)
C
C++++++++++++++++++++++++
C
C Version 1.0 uses the INQUIRE statement to find out the the record length
C of the direct access file before opening it. This procedure is non-standard,
C but seems to work for VAX machines.
C
C THE SUBROUTINE ALSO SETS THE VALUES OF NRECL, NRFILE, AND NAMFIL.
C *****************************************************************
C *****************************************************************
C
C THE PARAMETERS NAMFIL, NRECL, AND NRFILE ARE TO BE SET BY THE USER
C
C *****************************************************************
C NAMFIL IS THE EXTERNAL NAME OF THE BINARY EPHEMERIS FILE
CHARACTER*256 NAMFIL
c NAMFIL='JPLEPH'
C *****************************************************************
C NRECL=1 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN S.P. WORDS
C NRECL=4 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN BYTES
C (for a VAX, it is probably 1)
C
NRECL=4
C *****************************************************************
C NRFILE IS THE INTERNAL UNIT NUMBER USED FOR THE EPHEMERIS FILE
c NRFILE=12
C *****************************************************************
C FIND THE RECORD SIZE USING THE INQUIRE STATEMENT
c IRECSZ=0
INQUIRE(FILE=NAMFIL,RECL=IRECSZ)
C IF 'INQUIRE' DOES NOT WORK, USUALLY IRECSZ WILL BE LEFT AT 0
IF(IRECSZ .LE. 0) write(*,*)
. ' INQUIRE STATEMENT PROBABLY DID NOT WORK'
KSIZE=IRECSZ/NRECL
if(nrfile.eq.-99) stop !silence compiler warning
RETURN
END
C++++++++++++++++++++++++
C
SUBROUTINE FSIZER2(NRECL,KSIZE,NRFILE,NAMFIL)
C
C++++++++++++++++++++++++
C THIS SUBROUTINE OPENS THE FILE, 'NAMFIL', WITH A PHONY RECORD LENGTH, READS
C THE FIRST RECORD, AND USES THE INFO TO COMPUTE KSIZE, THE NUMBER OF SINGLE
C PRECISION WORDS IN A RECORD.
C
C THE SUBROUTINE ALSO SETS THE VALUES OF NRECL, NRFILE, AND NAMFIL.
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
SAVE
INTEGER OLDMAX
PARAMETER (OLDMAX = 400)
INTEGER NMAX
PARAMETER (NMAX = 1000)
CHARACTER*6 TTL(14,3),CNAM(NMAX)
CHARACTER*256 NAMFIL,jpleph_file_name
DIMENSION SS(3)
INTEGER IPT(3,13)
common/jplcom/jpleph_file_name
C *****************************************************************
C *****************************************************************
C
C THE PARAMETERS NRECL, NRFILE, AND NAMFIL ARE TO BE SET BY THE USER
C
C *****************************************************************
C NRECL=1 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN S.P. WORDS
C NRECL=4 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN BYTES
C (for UNIX, it is probably 4)
C
NRECL=4
C NRFILE IS THE INTERNAL UNIT NUMBER USED FOR THE EPHEMERIS FILE
NRFILE=12
C NAMFIL IS THE EXTERNAL NAME OF THE BINARY EPHEMERIS FILE
! NAMFIL='JPLEPH'
NAMFIL=jpleph_file_name
C *****************************************************************
C *****************************************************************
C ** OPEN THE DIRECT-ACCESS FILE AND GET THE POINTERS IN ORDER TO
C ** DETERMINE THE SIZE OF THE EPHEMERIS RECORD
MRECL=NRECL*1000
OPEN(NRFILE,
* FILE=NAMFIL,
* ACCESS='DIRECT',
* FORM='UNFORMATTED',
* RECL=MRECL,
* STATUS='OLD')
READ(NRFILE,REC=1)TTL,(CNAM(K),K=1,OLDMAX),SS,NCON,AU,EMRAT,
& ((IPT(I,J),I=1,3),J=1,12),NUMDE,(IPT(I,13),I=1,3)
CLOSE(NRFILE)
C FIND THE NUMBER OF EPHEMERIS COEFFICIENTS FROM THE POINTERS
KMX = 0
KHI = 0
DO I = 1,13
IF (IPT(1,I) .GT. KMX) THEN
KMX = IPT(1,I)
KHI = I
ENDIF
ENDDO
ND = 3
IF (KHI .EQ. 12) ND=2
KSIZE = 2*(IPT(1,KHI)+ND*IPT(2,KHI)*IPT(3,KHI)-1)
RETURN
END
C++++++++++++++++++++++++
C
SUBROUTINE FSIZER3(NRECL,KSIZE,NRFILE,NAMFIL)
C
C++++++++++++++++++++++++
C
C THE SUBROUTINE SETS THE VALUES OF NRECL, KSIZE, NRFILE, AND NAMFIL.
SAVE
CHARACTER*256 NAMFIL,jpleph_file_name
common/jplcom/jpleph_file_name
C *****************************************************************
C *****************************************************************
C
C THE PARAMETERS NRECL, NRFILE, AND NAMFIL ARE TO BE SET BY THE USER
C *****************************************************************
C NRECL=1 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN S.P. WORDS
C NRECL=4 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN BYTES
NRECL=4
C *****************************************************************
C NRFILE IS THE INTERNAL UNIT NUMBER USED FOR THE EPHEMERIS FILE (DEFAULT: 12)
NRFILE=12
C *****************************************************************
C NAMFIL IS THE EXTERNAL NAME OF THE BINARY EPHEMERIS FILE
! NAMFIL='JPLEPH'
NAMFIL=jpleph_file_name
C *****************************************************************
C KSIZE must be set by the user according to the ephemeris to be read
C For de200, set KSIZE to 1652
C For de405, set KSIZE to 2036
C For de406, set KSIZE to 1456
C For de414, set KSIZE to 2036
C For de418, set KSIZE to 2036
C For de421, set KSIZE to 2036
C For de422, set KSIZE to 2036
C For de423, set KSIZE to 2036
C For de424, set KSIZE to 2036
C For de430, set KSIZE to 2036
KSIZE = 2036
C *******************************************************************
RETURN
END
C++++++++++++++++++++++++++
C
SUBROUTINE PLEPH ( ET, NTARG, NCENT, RRD )
C
C++++++++++++++++++++++++++
C NOTE : Over the years, different versions of PLEPH have had a fifth argument:
C sometimes, an error return statement number; sometimes, a logical denoting
C whether or not the requested date is covered by the ephemeris. We apologize
C for this inconsistency; in this present version, we use only the four necessary
C arguments and do the testing outside of the subroutine.
C
C THIS SUBROUTINE READS THE JPL PLANETARY EPHEMERIS
C AND GIVES THE POSITION AND VELOCITY OF THE POINT 'NTARG'
C WITH RESPECT TO 'NCENT'.
C
C CALLING SEQUENCE PARAMETERS:
C
C ET = D.P. JULIAN EPHEMERIS DATE AT WHICH INTERPOLATION
C IS WANTED.
C
C ** NOTE THE ENTRY DPLEPH FOR A DOUBLY-DIMENSIONED TIME **
C THE REASON FOR THIS OPTION IS DISCUSSED IN THE
C SUBROUTINE STATE
C
C NTARG = INTEGER NUMBER OF 'TARGET' POINT.
C
C NCENT = INTEGER NUMBER OF CENTER POINT.
C
C THE NUMBERING CONVENTION FOR 'NTARG' AND 'NCENT' IS:
C
C 1 = MERCURY 8 = NEPTUNE
C 2 = VENUS 9 = PLUTO
C 3 = EARTH 10 = MOON
C 4 = MARS 11 = SUN
C 5 = JUPITER 12 = SOLAR-SYSTEM BARYCENTER
C 6 = SATURN 13 = EARTH-MOON BARYCENTER
C 7 = URANUS 14 = NUTATIONS (LONGITUDE AND OBLIQ)
C 15 = LIBRATIONS, IF ON EPH FILE
C
C (IF NUTATIONS ARE WANTED, SET NTARG = 14. FOR LIBRATIONS,
C SET NTARG = 15. SET NCENT=0.)
C
C RRD = OUTPUT 6-WORD D.P. ARRAY CONTAINING POSITION AND VELOCITY
C OF POINT 'NTARG' RELATIVE TO 'NCENT'. THE UNITS ARE AU AND
C AU/DAY. FOR LIBRATIONS THE UNITS ARE RADIANS AND RADIANS
C PER DAY. IN THE CASE OF NUTATIONS THE FIRST FOUR WORDS OF
C RRD WILL BE SET TO NUTATIONS AND RATES, HAVING UNITS OF
C RADIANS AND RADIANS/DAY.
C
C The option is available to have the units in km and km/sec.
C For this, set km=.true. in the STCOMX common block.
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
INTEGER NMAX
PARAMETER (NMAX = 1000)
DIMENSION RRD(6),ET2Z(2),ET2(2),PV(6,13)
DIMENSION PVST(6,11),PNUT(4)
DIMENSION SS(3),CVAL(NMAX),PVSUN(6),ZIPS(2)
DATA ZIPS/2*0.d0/
LOGICAL BSAVE,KM,BARY
LOGICAL FIRST
DATA FIRST/.TRUE./
INTEGER LIST(12),IPT(39),DENUM
COMMON/EPHHDR/CVAL,SS,AU,EMRAT,DENUM,NCON,IPT
COMMON/STCOMX/KM,BARY,PVSUN
C INITIALIZE ET2 FOR 'STATE' AND SET UP COMPONENT COUNT
C
ET2(1)=ET
ET2(2)=0.D0
GO TO 11
C ENTRY POINT 'DPLEPH' FOR DOUBLY-DIMENSIONED TIME ARGUMENT
C (SEE THE DISCUSSION IN THE SUBROUTINE STATE)
ENTRY DPLEPH(ET2Z,NTARG,NCENT,RRD)
ET2(1)=ET2Z(1)
ET2(2)=ET2Z(2)
11 IF(FIRST) CALL STATE(ZIPS,LIST,PVST,PNUT)
FIRST=.FALSE.
IF(NTARG .EQ. NCENT) RETURN
DO I=1,12
LIST(I)=0
ENDDO
C CHECK FOR NUTATION CALL
IF(NTARG.NE.14) GO TO 97
IF(IPT(35).GT.0) THEN
LIST(11)=2
CALL STATE(ET2,LIST,PVST,PNUT)
DO I=1,4
RRD(I)=PNUT(I)
ENDDO
RRD(5) = 0.d0
RRD(6) = 0.d0
RETURN
ELSE
DO I=1,4
RRD(I)=0.d0
ENDDO
WRITE(6,297)
297 FORMAT(' ***** NO NUTATIONS ON THE EPHEMERIS FILE *****')
STOP
ENDIF
C CHECK FOR LIBRATIONS
97 CONTINUE
DO I=1,6
RRD(I)=0.d0
ENDDO
IF(NTARG.NE.15) GO TO 98
IF(IPT(38).GT.0) THEN
LIST(12)=2
CALL STATE(ET2,LIST,PVST,PNUT)
DO I=1,6
RRD(I)=PVST(I,11)
ENDDO
RETURN
ELSE
WRITE(6,298)
298 FORMAT(' ***** NO LIBRATIONS ON THE EPHEMERIS FILE *****')
STOP
ENDIF
C FORCE BARYCENTRIC OUTPUT BY 'STATE'
98 BSAVE=BARY
BARY=.TRUE.
C SET UP PROPER ENTRIES IN 'LIST' ARRAY FOR STATE CALL
DO I=1,2
K=NTARG
IF(I .EQ. 2) K=NCENT
IF(K .LE. 10) LIST(K)=2
IF(K .EQ. 10) LIST(3)=2
IF(K .EQ. 3) LIST(10)=2
IF(K .EQ. 13) LIST(3)=2
ENDDO
C MAKE CALL TO STATE
CALL STATE(ET2,LIST,PVST,PNUT)
DO I=1,10
DO J = 1,6
PV(J,I) = PVST(J,I)
ENDDO
ENDDO
IF(NTARG .EQ. 11 .OR. NCENT .EQ. 11) THEN
DO I=1,6
PV(I,11)=PVSUN(I)
ENDDO
ENDIF
IF(NTARG .EQ. 12 .OR. NCENT .EQ. 12) THEN
DO I=1,6
PV(I,12)=0.D0
ENDDO
ENDIF
IF(NTARG .EQ. 13 .OR. NCENT .EQ. 13) THEN
DO I=1,6
PV(I,13) = PVST(I,3)
ENDDO
ENDIF
IF(NTARG*NCENT .EQ. 30 .AND. NTARG+NCENT .EQ. 13) THEN
DO I=1,6
PV(I,3)=0.D0
ENDDO
GO TO 99
ENDIF
IF(LIST(3) .EQ. 2) THEN
DO I=1,6
PV(I,3)=PVST(I,3)-PVST(I,10)/(1.D0+EMRAT)
ENDDO
ENDIF
IF(LIST(10) .EQ. 2) THEN
DO I=1,6
PV(I,10) = PV(I,3)+PVST(I,10)
ENDDO
ENDIF
99 DO I=1,6
RRD(I)=PV(I,NTARG)-PV(I,NCENT)
ENDDO
BARY=BSAVE
RETURN
END
C+++++++++++++++++++++++++++++++++
C
SUBROUTINE INTERP(BUF,T,NCF,NCM,NA,IFL,PV)
C
C+++++++++++++++++++++++++++++++++
C
C THIS SUBROUTINE DIFFERENTIATES AND INTERPOLATES A
C SET OF CHEBYSHEV COEFFICIENTS TO GIVE POSITION AND VELOCITY
C
C CALLING SEQUENCE PARAMETERS:
C
C INPUT:
C
C BUF 1ST LOCATION OF ARRAY OF D.P. CHEBYSHEV COEFFICIENTS OF POSITION
C
C T T(1) IS DP FRACTIONAL TIME IN INTERVAL COVERED BY
C COEFFICIENTS AT WHICH INTERPOLATION IS WANTED
C (0 .LE. T(1) .LE. 1). T(2) IS DP LENGTH OF WHOLE
C INTERVAL IN INPUT TIME UNITS.
C
C NCF # OF COEFFICIENTS PER COMPONENT
C
C NCM # OF COMPONENTS PER SET OF COEFFICIENTS
C
C NA # OF SETS OF COEFFICIENTS IN FULL ARRAY
C (I.E., # OF SUB-INTERVALS IN FULL INTERVAL)
C
C IFL INTEGER FLAG: =1 FOR POSITIONS ONLY
C =2 FOR POS AND VEL
C
C
C OUTPUT:
C
C PV INTERPOLATED QUANTITIES REQUESTED. DIMENSION
C EXPECTED IS PV(NCM,IFL), DP.
C
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
SAVE
C
DOUBLE PRECISION BUF(NCF,NCM,*),T(2),PV(NCM,*),PC(18),VC(18)
C
DATA NP/2/
DATA NV/3/
DATA TWOT/0.D0/
DATA PC(1),PC(2)/1.D0,0.D0/
DATA VC(2)/1.D0/
C
C ENTRY POINT. GET CORRECT SUB-INTERVAL NUMBER FOR THIS SET
C OF COEFFICIENTS AND THEN GET NORMALIZED CHEBYSHEV TIME
C WITHIN THAT SUBINTERVAL.
C
DNA=DBLE(NA)
DT1=DINT(T(1))
TEMP=DNA*T(1)
L=IDINT(TEMP-DT1)+1
C TC IS THE NORMALIZED CHEBYSHEV TIME (-1 .LE. TC .LE. 1)
TC=2.D0*(DMOD(TEMP,1.D0)+DT1)-1.D0
C CHECK TO SEE WHETHER CHEBYSHEV TIME HAS CHANGED,
C AND COMPUTE NEW POLYNOMIAL VALUES IF IT HAS.
C (THE ELEMENT PC(2) IS THE VALUE OF T1(TC) AND HENCE
C CONTAINS THE VALUE OF TC ON THE PREVIOUS CALL.)
IF(TC.NE.PC(2)) THEN
NP=2
NV=3
PC(2)=TC
TWOT=TC+TC
ENDIF
C
C BE SURE THAT AT LEAST 'NCF' POLYNOMIALS HAVE BEEN EVALUATED
C AND ARE STORED IN THE ARRAY 'PC'.
C
IF(NP.LT.NCF) THEN
DO 1 I=NP+1,NCF
PC(I)=TWOT*PC(I-1)-PC(I-2)
1 CONTINUE
NP=NCF
ENDIF
C
C INTERPOLATE TO GET POSITION FOR EACH COMPONENT
C
DO 2 I=1,NCM
PV(I,1)=0.D0
DO 3 J=NCF,1,-1
PV(I,1)=PV(I,1)+PC(J)*BUF(J,I,L)
3 CONTINUE
2 CONTINUE
IF(IFL.LE.1) RETURN
C
C IF VELOCITY INTERPOLATION IS WANTED, BE SURE ENOUGH
C DERIVATIVE POLYNOMIALS HAVE BEEN GENERATED AND STORED.
C
VFAC=(DNA+DNA)/T(2)
VC(3)=TWOT+TWOT
IF(NV.LT.NCF) THEN
DO 4 I=NV+1,NCF
VC(I)=TWOT*VC(I-1)+PC(I-1)+PC(I-1)-VC(I-2)
4 CONTINUE
NV=NCF
ENDIF
C
C INTERPOLATE TO GET VELOCITY FOR EACH COMPONENT
C
DO 5 I=1,NCM
PV(I,2)=0.D0
DO 6 J=NCF,2,-1
PV(I,2)=PV(I,2)+VC(J)*BUF(J,I,L)
6 CONTINUE
PV(I,2)=PV(I,2)*VFAC
5 CONTINUE
C
RETURN
C
END
C+++++++++++++++++++++++++
C
SUBROUTINE SPLIT(TT,FR)
C
C+++++++++++++++++++++++++
C
C THIS SUBROUTINE BREAKS A D.P. NUMBER INTO A D.P. INTEGER
C AND A D.P. FRACTIONAL PART.
C
C CALLING SEQUENCE PARAMETERS:
C
C TT = D.P. INPUT NUMBER
C
C FR = D.P. 2-WORD OUTPUT ARRAY.
C FR(1) CONTAINS INTEGER PART
C FR(2) CONTAINS FRACTIONAL PART
C
C FOR NEGATIVE INPUT NUMBERS, FR(1) CONTAINS THE NEXT
C MORE NEGATIVE INTEGER; FR(2) CONTAINS A POSITIVE FRACTION.
C
C CALLING SEQUENCE DECLARATIONS
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION FR(2)
C MAIN ENTRY -- GET INTEGER AND FRACTIONAL PARTS
FR(1)=DINT(TT)
FR(2)=TT-FR(1)
IF(TT.GE.0.D0 .OR. FR(2).EQ.0.D0) RETURN
C MAKE ADJUSTMENTS FOR NEGATIVE INPUT NUMBER
FR(1)=FR(1)-1.D0
FR(2)=FR(2)+1.D0
RETURN
END
C++++++++++++++++++++++++++++++++
C
SUBROUTINE STATE(ET2,LIST,PV,PNUT)
C
C++++++++++++++++++++++++++++++++
C
C THIS SUBROUTINE READS AND INTERPOLATES THE JPL PLANETARY EPHEMERIS FILE
C
C CALLING SEQUENCE PARAMETERS:
C
C INPUT:
C
C ET2 DP 2-WORD JULIAN EPHEMERIS EPOCH AT WHICH INTERPOLATION
C IS WANTED. ANY COMBINATION OF ET2(1)+ET2(2) WHICH FALLS
C WITHIN THE TIME SPAN ON THE FILE IS A PERMISSIBLE EPOCH.
C
C A. FOR EASE IN PROGRAMMING, THE USER MAY PUT THE
C ENTIRE EPOCH IN ET2(1) AND SET ET2(2)=0.
C
C B. FOR MAXIMUM INTERPOLATION ACCURACY, SET ET2(1) =
C THE MOST RECENT MIDNIGHT AT OR BEFORE INTERPOLATION
C EPOCH AND SET ET2(2) = FRACTIONAL PART OF A DAY
C ELAPSED BETWEEN ET2(1) AND EPOCH.
C
C C. AS AN ALTERNATIVE, IT MAY PROVE CONVENIENT TO SET
C ET2(1) = SOME FIXED EPOCH, SUCH AS START OF INTEGRATION,
C AND ET2(2) = ELAPSED INTERVAL BETWEEN THEN AND EPOCH.
C
C LIST 12-WORD INTEGER ARRAY SPECIFYING WHAT INTERPOLATION
C IS WANTED FOR EACH OF THE BODIES ON THE FILE.
C
C LIST(I)=0, NO INTERPOLATION FOR BODY I
C =1, POSITION ONLY
C =2, POSITION AND VELOCITY
C
C THE DESIGNATION OF THE ASTRONOMICAL BODIES BY I IS:
C
C I = 1: MERCURY
C = 2: VENUS
C = 3: EARTH-MOON BARYCENTER
C = 4: MARS
C = 5: JUPITER
C = 6: SATURN
C = 7: URANUS
C = 8: NEPTUNE
C = 9: PLUTO
C =10: GEOCENTRIC MOON
C =11: NUTATIONS IN LONGITUDE AND OBLIQUITY
C =12: LUNAR LIBRATIONS (IF ON FILE)
C
C OUTPUT:
C
C PV DP 6 X 11 ARRAY THAT WILL CONTAIN REQUESTED INTERPOLATED
C QUANTITIES (OTHER THAN NUTATION, STOERD IN PNUT).
C THE BODY SPECIFIED BY LIST(I) WILL HAVE ITS
C STATE IN THE ARRAY STARTING AT PV(1,I).
C (ON ANY GIVEN CALL, ONLY THOSE WORDS IN 'PV' WHICH ARE
C AFFECTED BY THE FIRST 10 'LIST' ENTRIES, AND BY LIST(12)
C IF LIBRATIONS ARE ON THE FILE, ARE SET.
C THE REST OF THE 'PV' ARRAYIS UNTOUCHED.)
C THE ORDER OF COMPONENTS STARTING IN PV(1,I) IS: X,Y,Z,DX,DY,DZ.
C
C ALL OUTPUT VECTORS ARE REFERENCED TO THE EARTH MEAN
C EQUATOR AND EQUINOX OF J2000 IF THE DE NUMBER IS 200 OR
C GREATER; OF B1950 IF THE DE NUMBER IS LESS THAN 200.
C
C THE MOON STATE IS ALWAYS GEOCENTRIC; THE OTHER NINE STATES
C ARE EITHER HELIOCENTRIC OR SOLAR-SYSTEM BARYCENTRIC,
C DEPENDING ON THE SETTING OF COMMON FLAGS (SEE BELOW).
C
C LUNAR LIBRATIONS, IF ON FILE, ARE PUT INTO PV(K,11) IF
C LIST(12) IS 1 OR 2.
C
C NUT DP 4-WORD ARRAY THAT WILL CONTAIN NUTATIONS AND RATES,
C DEPENDING ON THE SETTING OF LIST(11). THE ORDER OF
C QUANTITIES IN NUT IS:
C
C D PSI (NUTATION IN LONGITUDE)
C D EPSILON (NUTATION IN OBLIQUITY)
C D PSI DOT
C D EPSILON DOT
C
C * STATEMENT # FOR ERROR RETURN, IN CASE OF EPOCH OUT OF
C RANGE OR I/O ERRORS.
C
C COMMON AREA STCOMX:
C
C KM LOGICAL FLAG DEFINING PHYSICAL UNITS OF THE OUTPUT
C STATES. KM = .TRUE., KM AND KM/SEC
C = .FALSE., AU AND AU/DAY
C DEFAULT VALUE = .FALSE. (KM DETERMINES TIME UNIT
C FOR NUTATIONS AND LIBRATIONS. ANGLE UNIT IS ALWAYS RADIANS.)
C
C BARY LOGICAL FLAG DEFINING OUTPUT CENTER.
C ONLY THE 9 PLANETS ARE AFFECTED.
C BARY = .TRUE. =\ CENTER IS SOLAR-SYSTEM BARYCENTER
C = .FALSE. =\ CENTER IS SUN
C DEFAULT VALUE = .FALSE.
C
C PVSUN DP 6-WORD ARRAY CONTAINING THE BARYCENTRIC POSITION AND
C VELOCITY OF THE SUN.
C
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
SAVE
INTEGER OLDMAX
PARAMETER ( OLDMAX = 400)
INTEGER NMAX
PARAMETER ( NMAX = 1000)
DIMENSION ET2(2),PV(6,11),PNUT(4),T(2),PJD(4),BUF(1500),
. SS(3),CVAL(NMAX),PVSUN(6)
INTEGER LIST(12),IPT(3,13)
LOGICAL FIRST
DATA FIRST/.TRUE./
CHARACTER*6 TTL(14,3),CNAM(NMAX)
CHARACTER*256 NAMFIL
LOGICAL KM,BARY
COMMON/EPHHDR/CVAL,SS,AU,EMRAT,NUMDE,NCON,IPT
COMMON/CHRHDR/CNAM,TTL
COMMON/STCOMX/KM,BARY,PVSUN
C
C ENTRY POINT - 1ST TIME IN, GET POINTER DATA, ETC., FROM EPH FILE
C
IF(FIRST) THEN
FIRST=.FALSE.
C ************************************************************************
C ************************************************************************
C THE USER MUST SELECT ONE OF THE FOLLOWING BY DELETING THE 'C' IN COLUMN 1
C ************************************************************************
C CALL FSIZER1(NRECL,KSIZE,NRFILE,NAMFIL)
C CALL FSIZER2(NRECL,KSIZE,NRFILE,NAMFIL)
CALL FSIZER3(NRECL,KSIZE,NRFILE,NAMFIL)
IF(NRECL .EQ. 0) WRITE(*,*)' ***** FSIZER IS NOT WORKING *****'
C ************************************************************************
C ************************************************************************
IRECSZ=NRECL*KSIZE
NCOEFFS=KSIZE/2
OPEN(NRFILE,
* FILE=NAMFIL,
* ACCESS='DIRECT',
* FORM='UNFORMATTED',
* RECL=IRECSZ,
* STATUS='OLD')
READ(NRFILE,REC=1)TTL,(CNAM(K),K=1,OLDMAX),SS,NCON,AU,EMRAT,
& ((IPT(I,J),I=1,3),J=1,12),NUMDE,(IPT(I,13),I=1,3)
& ,(CNAM(L),L=OLDMAX+1,NCON)
IF(NCON .LE. OLDMAX)THEN
READ(NRFILE,REC=2)(CVAL(I),I=1,OLDMAX)
ELSE
READ(NRFILE,REC=2)(CVAL(I),I=1,NCON)
ENDIF
NRL=0
ENDIF
C ********** MAIN ENTRY POINT **********
IF(ET2(1) .EQ. 0.D0) RETURN
S=ET2(1)-.5D0
CALL SPLIT(S,PJD(1))
CALL SPLIT(ET2(2),PJD(3))
PJD(1)=PJD(1)+PJD(3)+.5D0
PJD(2)=PJD(2)+PJD(4)
CALL SPLIT(PJD(2),PJD(3))
PJD(1)=PJD(1)+PJD(3)
C ERROR RETURN FOR EPOCH OUT OF RANGE
IF(PJD(1)+PJD(4).LT.SS(1) .OR. PJD(1)+PJD(4).GT.SS(2)) GO TO 98
C CALCULATE RECORD # AND RELATIVE TIME IN INTERVAL
NR=IDINT((PJD(1)-SS(1))/SS(3))+3
IF(PJD(1).EQ.SS(2)) NR=NR-1
tmp1 = DBLE(NR-3)*SS(3) + SS(1)
tmp2 = PJD(1) - tmp1
T(1) = (tmp2 + PJD(4))/SS(3)
C READ CORRECT RECORD IF NOT IN CORE
IF(NR.NE.NRL) THEN
NRL=NR
READ(NRFILE,REC=NR,ERR=99)(BUF(K),K=1,NCOEFFS)
ENDIF
IF(KM) THEN
T(2)=SS(3)*86400.D0
AUFAC=1.D0
ELSE
T(2)=SS(3)
AUFAC=1.D0/AU
ENDIF
C INTERPOLATE SSBARY SUN
CALL INTERP(BUF(IPT(1,11)),T,IPT(2,11),3,IPT(3,11),2,PVSUN)
DO I=1,6
PVSUN(I)=PVSUN(I)*AUFAC
ENDDO
C CHECK AND INTERPOLATE WHICHEVER BODIES ARE REQUESTED
DO 4 I=1,10
IF(LIST(I).EQ.0) GO TO 4
CALL INTERP(BUF(IPT(1,I)),T,IPT(2,I),3,IPT(3,I),
& LIST(I),PV(1,I))
DO J=1,6
IF(I.LE.9 .AND. .NOT.BARY) THEN
PV(J,I)=PV(J,I)*AUFAC-PVSUN(J)
ELSE
PV(J,I)=PV(J,I)*AUFAC
ENDIF
ENDDO
4 CONTINUE
C DO NUTATIONS IF REQUESTED (AND IF ON FILE)
IF(LIST(11).GT.0 .AND. IPT(2,12).GT.0)
* CALL INTERP(BUF(IPT(1,12)),T,IPT(2,12),2,IPT(3,12),
* LIST(11),PNUT)
C GET LIBRATIONS IF REQUESTED (AND IF ON FILE)
IF(LIST(12).GT.0 .AND. IPT(2,13).GT.0)
* CALL INTERP(BUF(IPT(1,13)),T,IPT(2,13),3,IPT(3,13),
* LIST(12),PV(1,11))
RETURN
98 WRITE(*,198)ET2(1)+ET2(2),SS(1),SS(2)
198 FORMAT(' *** Requested JED,',f12.2,
* ' not within ephemeris limits,',2f12.2,' ***')
STOP
99 WRITE(*,'(2F12.2,A80)')ET2,'ERROR RETURN IN STATE'
STOP
END
C+++++++++++++++++++++++++++++
C
SUBROUTINE CONST(NAM,VAL,SSS,N)
C
C+++++++++++++++++++++++++++++
C
C THIS ENTRY OBTAINS THE CONSTANTS FROM THE EPHEMERIS FILE
C
C CALLING SEQEUNCE PARAMETERS (ALL OUTPUT):
C
C NAM = CHARACTER*6 ARRAY OF CONSTANT NAMES
C
C VAL = D.P. ARRAY OF VALUES OF CONSTANTS
C
C SSS = D.P. JD START, JD STOP, STEP OF EPHEMERIS
C
C N = INTEGER NUMBER OF ENTRIES IN 'NAM' AND 'VAL' ARRAYS
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
SAVE
INTEGER NMAX
PARAMETER (NMAX = 1000)
CHARACTER*6 NAM(*),TTL(14,3),CNAM(NMAX)
DOUBLE PRECISION VAL(*),SSS(3),SS(3),CVAL(NMAX),ZIPS(2)
DOUBLE PRECISION PVST(6,11),PNUT(4)
DATA ZIPS/2*0.d0/
INTEGER IPT(3,13),DENUM,LIST(12)
logical first
data first/.true./
COMMON/EPHHDR/CVAL,SS,AU,EMRAT,DENUM,NCON,IPT
COMMON/CHRHDR/CNAM,TTL
C CALL STATE TO INITIALIZE THE EPHEMERIS AND READ IN THE CONSTANTS
IF(FIRST) CALL STATE(ZIPS,LIST,PVST,PNUT)
first=.false.
N=NCON
DO I=1,3
SSS(I)=SS(I)
ENDDO
DO I=1,N
NAM(I)=CNAM(I)
VAL(I)=CVAL(I)
ENDDO
RETURN
END
lib/jt49sim.f90
-191
View File
@@ -1,191 +0,0 @@
program jt49sim
! Generate simulated data for testing JT4 and JT9
use wavhdr
use packjt
use jt4
parameter (NMAX=60*12000) ! = 648,000
parameter (NFFT=10*65536,NH=NFFT/2)
type(hdr) h !Header for .wav file
integer*2 iwave(NMAX) !Generated waveform
integer*4 itone(206) !Channel symbols (values 0-8)
real*4 xnoise(NMAX) !Generated random noise
real*4 dat(NMAX) !Generated real data
complex cdat(NMAX) !Generated complex waveform
complex cspread(0:NFFT-1) !Complex amplitude for Rayleigh fading
complex z
real*8 f0,dt,twopi,phi,dphi,baud,fsample,freq,dnsps
character message*22,fname*11,csubmode*2,arg*12
character msgsent*22
nargs=iargc()
if(nargs.ne. 7) then
print *, 'Usage: jt49sim "msg" nA-nE Nsigs fDop DT Nfiles SNR'
print *, 'Example jt49sim "K1ABC W9XYZ EN37" 4G 10 0.2 0.0 1 0'
print *, 'Example jt49sim "K1ABC W9XYZ EN37" 9A 1 0.0 0.0 1 -20'
print *, 'Use msg=@nnnn to generate a tone at nnnn Hz:'
print *, 'Example jt49sim "@1500" 9A 1 10.0 0.0 1 -20'
go to 999
endif
call getarg(1,message)
call fmtmsg(message, iz)
call getarg(2,csubmode)
imode=ichar(csubmode(1:1)) - ichar('0')
nsubmode=ichar(csubmode(2:2)) - ichar('A')
if(imode.ne.4 .and. imode.ne.9) go to 999
if(nsubmode.lt.0 .or. nsubmode.gt.7) go to 999
call getarg(3,arg)
read(arg,*) nsigs
call getarg(4,arg)
read(arg,*) fspread
call getarg(5,arg)
read(arg,*) xdt
call getarg(6,arg)
read(arg,*) nfiles
call getarg(7,arg)
read(arg,*) snrdb
rms=100.
fsample=12000.d0 !Sample rate (Hz)
dt=1.d0/fsample !Sample interval (s)
twopi=8.d0*atan(1.d0)
npts=60*12000 !Total samples in .wav file
h=default_header(12000,npts)
dfsig=2000.0/nsigs !Freq spacing between sigs in file (Hz)
ichk=0
if(imode.eq.4) then
nsym=206 !Number of channel symbols (JT4)
dnsps=12000.d0/4.375d0
baud=12000.d0/dnsps !Keying rate = 1.7361111111
else if(imode.eq.9) then
nsym=85 !Number of channel symbols (JT9)
dnsps=6912.d0 !Samples per symbol
baud=12000.d0/dnsps !Keying rate = 1.736...
endif
write(*,1000)
1000 format('File Sig Freq Mode S/N DT Dop Message'/60('-'))
do ifile=1,nfiles !Loop over requested number of files
write(fname,1002) ifile !Output filename
1002 format('000000_',i4.4)
open(10,file=fname//'.wav',access='stream',status='unknown')
xnoise=0.
cdat=0.
if(snrdb.lt.90) then
do i=1,npts
xnoise(i)=gran() !Generate gaussian noise
enddo
endif
do isig=1,nsigs !Generate requested number of sigs
if(mod(nsigs,2).eq.0) f0=1500.0 + dfsig*(isig-0.5-nsigs/2)
if(mod(nsigs,2).eq.1) f0=1500.0 + dfsig*(isig-(nsigs+1)/2)
if(nsigs.eq.1) f0=1000.0
xsnr=snrdb
if(snrdb.eq.0.0) xsnr=-20 - isig
if(imode.eq.4) call gen4(message,ichk,msgsent,itone,itype)
if(imode.eq.9) call gen9(message,ichk,msgsent,itone,itype)
bandwidth_ratio=2500.0/6000.0
sig=sqrt(2*bandwidth_ratio)*10.0**(0.05*xsnr)
if(xsnr.gt.90.0) sig=1.0
write(*,1020) ifile,isig,f0,csubmode,xsnr,xdt,fspread,message
1020 format(i4,i4,f10.3,2x,a2,2x,f5.1,f6.2,f6.1,1x,a22)
phi=0.d0
dphi=0.d0
k=(xdt+1.0)*12000 !Start audio at t = xdt + 1.0 s
isym0=-99
do i=1,npts !Add this signal into cdat()
isym=i/dnsps + 1
if(isym.gt.nsym) exit
if(isym.ne.isym0) then
if(message(1:1).eq.'@') then
read(message(2:),*) freq
else
if(imode.eq.4) freq=f0 + itone(isym)*baud*nch(1+nsubmode) !JT4
if(imode.eq.9) freq=f0 + itone(isym)*baud*(2**nsubmode) !JT9
endif
dphi=twopi*freq*dt
isym0=isym
endif
phi=phi + dphi
if(phi.gt.twopi) phi=phi-twopi
xphi=phi
z=cmplx(cos(xphi),sin(xphi))
k=k+1
if(k.ge.1) cdat(k)=cdat(k) + sig*z
enddo
enddo
if(fspread.ne.0) then !Apply specified Doppler spread
df=12000.0/nfft
twopi=8*atan(1.0)
cspread(0)=1.0
cspread(NH)=0.
b=6.0 !Lorenzian 3/28 onward
do i=1,NH
f=i*df
x=b*f/fspread
z=0.
a=0.
if(x.lt.3.0) then !Cutoff beyond x=3
a=sqrt(1.111/(1.0+x*x)-0.1) !Lorentzian
call random_number(r1)
phi1=twopi*r1
z=a*cmplx(cos(phi1),sin(phi1))
endif
cspread(i)=z
z=0.
if(x.lt.50.0) then
call random_number(r2)
phi2=twopi*r2
z=a*cmplx(cos(phi2),sin(phi2))
endif
cspread(NFFT-i)=z
enddo
do i=0,NFFT-1
f=i*df
if(i.gt.NH) f=(i-nfft)*df
s=real(cspread(i))**2 + aimag(cspread(i))**2
! write(13,3000) i,f,s,cspread(i)
!3000 format(i5,f10.3,3f12.6)
enddo
! s=real(cspread(0))**2 + aimag(cspread(0))**2
! write(13,3000) 1024,0.0,s,cspread(0)
call four2a(cspread,NFFT,1,1,1) !Transform to time domain
sum=0.
do i=0,NFFT-1
p=real(cspread(i))**2 + aimag(cspread(i))**2
sum=sum+p
enddo
avep=sum/NFFT
fac=sqrt(1.0/avep)
cspread=fac*cspread !Normalize to constant avg power
cdat(1:NFFT)=cspread*cdat(1:NFFT) !Apply Rayleigh fading
! do i=0,NFFT-1
! p=real(cspread(i))**2 + aimag(cspread(i))**2
! write(14,3010) i,p,cspread(i)
!3010 format(i8,3f12.6)
! enddo
endif
dat=aimag(cdat) + xnoise !Add the generated noise
fac=32767.0/nsigs
if(snrdb.ge.90.0) iwave(1:npts)=nint(fac*dat(1:npts))
if(snrdb.lt.90.0) iwave(1:npts)=nint(rms*dat(1:npts))
write(10) h,iwave(1:npts) !Save the .wav file
close(10)
enddo
999 end program jt49sim
lib/jt4a.f90
-2
View File
@@ -1,2 +0,0 @@
! The contents of this file have been migrated to lib/jt4_decode.f90
lib/jt4code.f90
-50
View File
@@ -1,50 +0,0 @@
program jt4code
! Provides examples of message packing, bit and symbol ordering,
! convolutional encoding, and other necessary details of the JT4
! protocol.
use jt4
use packjt
character*22 msg,decoded,bad*1,msgtype*13
integer i4tone(206)
include 'testmsg.f90'
nargs=iargc()
if(nargs.ne.1) then
print*,'Usage: jt4code "message"'
print*,' jt4code -t'
go to 999
endif
call getarg(1,msg)
nmsg=1
if(msg(1:2).eq."-t") nmsg=NTEST
write(*,1010)
1010 format(" Message Decoded Err? Type"/ &
74("-"))
do imsg=1,nmsg
if(nmsg.gt.1) msg=testmsg(imsg)
call fmtmsg(msg,iz) !To upper case, collapse multiple blanks
ichk=0
call gen4(msg,ichk,decoded,i4tone,itype)
msgtype=""
if(itype.eq.1) msgtype="Std Msg"
if(itype.eq.2) msgtype="Type 1 prefix"
if(itype.eq.3) msgtype="Type 1 suffix"
if(itype.eq.4) msgtype="Type 2 prefix"
if(itype.eq.5) msgtype="Type 2 suffix"
if(itype.eq.6) msgtype="Free text"
bad=" "
if(decoded.ne.msg) bad="*"
write(*,1020) imsg,msg,decoded,bad,itype,msgtype
1020 format(i2,'.',2x,a22,2x,a22,3x,a1,i3,": ",a13)
enddo
if(nmsg.eq.1) write(*,1030) i4tone
1030 format(/'Channel symbols'/(30i2))
999 end program jt4code
lib/jtmsg.f90
-131
View File
@@ -1,131 +0,0 @@
subroutine jtmsg(msg,iflag)
! Attempts to identify false decodes in JT-style messages.
! Returns iflag with sum of bits as follows:
! ------------------------------------------
! 1 Grid/Report invalid
! 2 Second callsign invalid
! 4 First callsign invalid
! 8 Very unlikely free text
! 16 Questionable free text
! 0 Message is probably OK
! ------------------------------------------
character*22 msg,t
character*13 w1,w2,w3,w
character*6 bc1,bc2,bc3
character*1 c
logical c1ok,c2ok,c3ok,isdigit,isletter,isgrid4
! Statement functions
isdigit(c)=(ichar(c).ge.ichar('0')) .and. (ichar(c).le.ichar('9'))
isletter(c)=(ichar(c).ge.ichar('A')) .and. (ichar(c).le.ichar('Z'))
isgrid4(w)=(len_trim(w).eq.4 .and. &
ichar(w(1:1)).ge.ichar('A') .and. ichar(w(1:1)).le.ichar('R') .and. &
ichar(w(2:2)).ge.ichar('A') .and. ichar(w(2:2)).le.ichar('R') .and. &
ichar(w(3:3)).ge.ichar('0') .and. ichar(w(3:3)).le.ichar('9') .and. &
ichar(w(4:4)).ge.ichar('0') .and. ichar(w(4:4)).le.ichar('9'))
t=trim(msg) !Temporary copy of msg
nt=len_trim(t)
! Check for standard messages
! Insert underscore in "CQ AA " to "CQ ZZ ", "CQ nnn " to make them one word.
if(t(1:3).eq.'CQ ' .and. isletter(t(4:4)) .and. &
isletter(t(5:5)) .and. t(6:6).eq.' ') t(3:3)='_'
if(t(1:3).eq.'CQ ' .and. isdigit(t(4:4)) .and. &
isdigit(t(5:5)) .and. isdigit(t(6:6)) .and. t(7:7).eq.' ') t(3:3)='_'
! Parse first three words
w1=' '
w2=' '
w3=' '
i1=index(t,' ')
if(i1.gt.0) w1(1:i1-1)=t(1:i1-1)
t=t(i1+1:)
i2=index(t,' ')
if(i2.gt.0) w2(1:i2-1)=t(1:i2-1)
t=t(i2+1:)
i3=index(t,' ')
if(i3.gt.0) w3(1:i3-1)=t(1:i3-1)
if(w1(1:3).eq.'CQ ' .or. w1(1:3).eq.'CQ_' .or. w1(1:3).eq.'DE ' .or. &
w1(1:4).eq.'QRZ ') then
! CQ/DE/QRZ: Should have one good callsign in w2 and maybe a grid/rpt in w3
call chkcall(w2,bc2,c2ok)
iflag=0
if(.not.c2ok) iflag=iflag+2
if(len_trim(w3).ne.0 .and. (.not.isgrid4(w3))) iflag=iflag+1
if(w1(1:3).eq.'DE ' .and. c2ok) iflag=0
if(iflag.eq.0) return
endif
! Check for two calls and maybe a grid, rpt, R+rpr, RRR, or 73
iflag=0
call chkcall(w1,bc1,c1ok)
call chkcall(w2,bc2,c2ok)
if(.not.c1ok) iflag=iflag+4
if(.not.c2ok) iflag=iflag+2
if(len_trim(w3).ne.0 .and. (.not.isgrid4(w3)) .and. &
w3(1:1).ne.'+' .and. w3(1:1).ne.'-' .and. &
w3(1:2).ne.'R+' .and. w3(1:2).ne.'R-' .and. &
w3(1:3).ne.'73 ' .and. w3(1:4).ne.'RRR ') iflag=iflag+1
call chkcall(w3,bc3,c3ok)
! Allow(?) non-standard messages of the form CQ AS OC K1JT
if(w1(1:3).eq.'CQ_'.and.isletter(w2(1:1)).and.isletter(w2(2:2)).and. &
w2(3:3).eq.' '.and.c3ok) iflag=0
if(iflag.eq.0 .or. nt.gt.13) return
! Check for plausible free text
nc=0
np=0
do i=1,13
c=msg(i:i)
if(c.ne.' ') nc=nc+1 !Number of non-blank characters
if(c.eq.'+') np=np+1 !Number of punctuation characters
if(c.eq.'-') np=np+1
if(c.eq.'.') np=np+1
if(c.eq.'/') np=np+1
if(c.eq.'?') np=np+1
enddo
nb=13-nc !Number of blanks
iflag=16 !Mark as potentially questionable
if(nc.ge.12 .or. (nc.ge.11 .and. np.gt.0)) then
iflag=8 !Unlikely free text, flag it
endif
! Save messages containing some common words
if(msg(1:3).eq.'CQ ') iflag=0
if(index(msg,'DE ').gt.0) iflag=0
if(index(msg,'TU ').gt.0) iflag=0
if(index(msg,' TU').gt.0) iflag=0
if(index(msg,'73 ').gt.0) iflag=0
if(index(msg,' 73').gt.0) iflag=0
if(index(msg,'TNX').gt.0) iflag=0
if(index(msg,'THX').gt.0) iflag=0
if(index(msg,'EQSL').gt.0) iflag=0
if(index(msg,'LOTW').gt.0) iflag=0
if(index(msg,'DECOD').gt.0) iflag=0
if(index(msg,'CHK').gt.0) iflag=0
if(index(msg,'CLK').gt.0) iflag=0
if(index(msg,'CLOCK').gt.0) iflag=0
if(index(msg,'LOG').gt.0) iflag=0
if(index(msg,'QRM').gt.0) iflag=0
if(index(msg,'QSY').gt.0) iflag=0
if(index(msg,'TEST').gt.0) iflag=0
if(index(msg,'CQDX').gt.0) iflag=0
if(index(msg,'CALL').gt.0) iflag=0
if(index(msg,'QRZ').gt.0) iflag=0
if(index(msg,'AUTO').gt.0) iflag=0
if(index(msg,'PHOTO').gt.0) iflag=0
if(index(msg,'HYBRID').gt.0) iflag=0
if(c1ok .and. w1(1:6).eq.bc1) iflag=0
if(c2ok .and. w2(1:6).eq.bc2) iflag=0
if(nb.ge.4) iflag=0
return
end subroutine jtmsg
lib/libration.f90
-38
View File
@@ -1,38 +0,0 @@
subroutine libration(jd,RA,Dec,xl,b)
! Compute optical libration of moon at jd: that is, the sub-observer
! point (xl,b) in selenographic coordinates. RA and Dec are
! topocentric values.
implicit real*8 (a-h,o-z)
parameter (RADS=0.0174532925199433d0)
parameter (TWOPI=6.28318530717959d0)
real*8 jd,j2000,mjd,lambda
j2000=2451545.0d0
RA2000=RA
Dec2000=Dec
year=2000.0d0 + (jd-j2000)/365.25d0
mjd=jd-2400000.d0
call sla_PRECES('FK5',year,2000.d0,RA2000,Dec2000)
call sla_EQECL(RA2000,Dec2000,mjd,lambda,beta)
day=jd - j2000
t = day / 36525.d0
xi = 1.54242 * RADS
ft = 93.2720993 + 483202.0175273 * t - .0034029 * t * t
b= ft / 360
a = 360 * (b - floor(b))
if (a.lt.0.d0) a = 360 + a;
f=a/57.2957795131d0
omega=sla_dranrm(2.182439196d0 - t*33.7570446085d0 + t*t*3.6236526d-5)
w = lambda - omega
y = sin(w) * cos(beta) * cos(xi) - sin(beta) * sin(xi)
x = cos(w) * cos(beta)
a = sla_dranrm(atan2(y, x))
xl = a - f
if(xl.lt.-0.25*TWOPI) xl=xl+TWOPI !Fix 'round the back' angles
if(xl.gt.0.25*TWOPI) xl=xl-TWOPI
b = asin(-sin(w) * cos(beta) * sin(xi) - sin(beta) * cos(xi))
return
end subroutine libration
lib/lorentzian.f90
-102
View File
@@ -1,102 +0,0 @@
subroutine lorentzian(y,npts,a)
! Input: y(npts); assume x(i)=i, i=1,npts
! Output: a(5)
! a(1) = baseline
! a(2) = amplitude
! a(3) = x0
! a(4) = width
! a(5) = chisqr
real y(npts)
real a(5)
real deltaa(4)
a=0.
df=12000.0/8192.0 !df = 1.465 Hz
width=0.
ipk=0
ymax=-1.e30
do i=1,npts
if(y(i).gt.ymax) then
ymax=y(i)
ipk=i
endif
! write(50,3001) i,i*df,y(i)
!3001 format(i6,2f12.3)
enddo
! base=(sum(y(ipk-149:ipk-50)) + sum(y(ipk+51:ipk+150)))/200.0
base=(sum(y(1:20)) + sum(y(npts-19:npts)))/40.0
stest=ymax - 0.5*(ymax-base)
ssum=y(ipk)
do i=1,50
if(ipk+i.gt.npts) exit
if(y(ipk+i).lt.stest) exit
ssum=ssum + y(ipk+i)
enddo
do i=1,50
if(ipk-i.lt.1) exit
if(y(ipk-i).lt.stest) exit
ssum=ssum + y(ipk-i)
enddo
ww=ssum/y(ipk)
width=2
t=ww*ww - 5.67
if(t.gt.0.0) width=sqrt(t)
a(1)=base
a(2)=ymax-base
a(3)=ipk
a(4)=width
! Now find Lorentzian parameters
deltaa(1)=0.1
deltaa(2)=0.1
deltaa(3)=1.0
deltaa(4)=1.0
nterms=4
! Start the iteration
chisqr=0.
chisqr0=1.e6
do iter=1,5
do j=1,nterms
chisq1=fchisq0(y,npts,a)
fn=0.
delta=deltaa(j)
10 a(j)=a(j)+delta
chisq2=fchisq0(y,npts,a)
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=fchisq0(y,npts,a)
if(chisq3.lt.chisq2) then
chisq1=chisq2
chisq2=chisq3
go to 20
endif
! 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.
! write(*,4000) iter,j,a,chisq2
!4000 format(i1,i2,4f10.4,f11.3)
enddo
chisqr=fchisq0(y,npts,a)
! write(*,4000) 0,0,a,chisqr
if(chisqr/chisqr0.gt.0.99) exit
chisqr0=chisqr
enddo
a(5)=chisqr
return
end subroutine lorentzian
lib/lpf1.f90
-29
View File
@@ -1,29 +0,0 @@
subroutine lpf1(dd,jz,dat,jz2)
parameter (NFFT1=64*11025,NFFT2=32*11025)
real dd(jz)
real dat(jz)
real x(NFFT1)
complex cx(0:NFFT1/2)
equivalence (x,cx)
save x,cx
fac=1.0/float(NFFT1)
x(1:jz)=fac*dd(1:jz)
x(jz+1:NFFT1)=0.0
call four2a(cx,NFFT1,1,-1,0) !Forwarxd FFT, r2c
cx(NFFT2/2:)=0.0
! df=11025.0/NFFT1
! do i=1,NFFT1/2
! sx=real(cx(i))**2 + aimag(cx(i))**2
! write(50,3000) i*df,sx
!3000 format(f15.6,e12.3)
! enddo
call four2a(cx,NFFT2,1,1,-1) !Inverse FFT, c2r
jz2=jz/2
dat(1:jz2)=x(1:jz2)
return
end subroutine lpf1
lib/makepings.f90
-26
View File
@@ -1,26 +0,0 @@
subroutine makepings(pings,npts,width,sig)
real pings(npts)
real*8 t
real t0(14)
iping0=-999
dt=1.0/12000.0
do i=1,14
t0(i)=i !Make pings at t=1, 2, ... 14 s.
enddo
amp=sig
do i=1,npts
iping=min(max(1,i/12000),14)
t=(i*dt-t0(iping))/width
if(t.lt.0.d0 .and. t.lt.10.0) then !????
fac=0.
else
fac=2.718*t*dexp(-t)
endif
pings(i)=fac*amp
enddo
return
end subroutine makepings
lib/mixlpf.f90
-25
View File
@@ -1,25 +0,0 @@
subroutine mixlpf(x1,nbfo,c0)
real*4 x1(512)
real*8 twopi,phi,dphi
complex c1(512),c2(105+512)
complex c0(64)
data phi/0.d0/
save phi,c2
twopi=8.d0*atan(1.d0)
dphi=twopi*nbfo/12000.d0
do i=1,512
phi=phi+dphi
if(phi.gt.twopi) phi=phi-twopi
xphi=phi
c1(i)=x1(i)*cmplx(cos(xphi),sin(xphi))
enddo
c2(106:105+512)=c1
call fil3c(c2,105+512,c0,n2)
c2(1:105)=c1(512-104:512) !Save 105 trailing samples
return
end subroutine mixlpf
lib/moondopjpl.f90
-43
View File
@@ -1,43 +0,0 @@
subroutine MoonDopJPL(nyear,month,nday,uth4,lon4,lat4,RAMoon4, &
DecMoon4,LST4,HA4,AzMoon4,ElMoon4,vr4,dist4)
implicit real*8 (a-h,o-z)
real*4 uth4 !UT in hours
real*4 lon4 !East longitude, degrees
real*4 lat4 !Latitude, degrees
real*4 RAMoon4 !Topocentric RA of moon, hours
real*4 DecMoon4 !Topocentric Dec of Moon, degrees
real*4 LST4 !Locat sidereal time, hours
real*4 HA4 !Local Hour angle, degrees
real*4 AzMoon4 !Topocentric Azimuth of moon, degrees
real*4 ElMoon4 !Topocentric Elevation of moon, degrees
real*4 vr4 !Radial velocity of moon wrt obs, km/s
real*4 dist4 !Echo time, seconds
twopi=8.d0*atan(1.d0) !Define some constants
rad=360.d0/twopi
clight=2.99792458d5
call sla_CLDJ(nyear,month,nday,djutc,j)
djutc=djutc + uth4/24.d0
dut=-0.460d0
east_long=lon4/rad
geodetic_lat=lat4/rad
height=40.
nspecial=0
call ephem(djutc,dut,east_long,geodetic_lat,height,nspecial, &
RA,Dec,Az,El,techo,dop,fspread_1GHz,vr)
RAMoon4=RA
DecMoon4=Dec
LST4=0. !These two variables not presently used
HA4=0.
AzMoon4=Az*rad
ElMoon4=El*rad
vr4=vr
dist4=techo
return
end subroutine MoonDopJPL
lib/noisegen.f90
-16
View File
@@ -1,16 +0,0 @@
subroutine noisegen(d4,nmax)
real*4 d4(4,nmax)
call init_random_seed() ! seed Fortran RANDOM_NUMBER generator
call sgran() ! see C rand generator (used in gran)
do i=1,nmax
d4(1,i)=gran()
d4(2,i)=gran()
d4(3,i)=gran()
d4(4,i)=gran()
enddo
return
end subroutine noisegen
lib/peakup.f90
-8
View File
@@ -1,8 +0,0 @@
subroutine peakup(ym,y0,yp,dx)
b=(yp-ym)/2.0
c=(yp+ym-2.0*y0)/2.0
dx=-b/(2.0*c)
return
end subroutine peakup
lib/ps4.f90
-27
View File
@@ -1,27 +0,0 @@
subroutine ps4(dat,nfft,s)
parameter (NMAX=2520+2)
parameter (NHMAX=NMAX/2-1)
real dat(nfft)
real dat2(NMAX)
real s(NHMAX)
complex c(0:NMAX)
equivalence(dat2,c)
nh=nfft/2
do i=1,nh
dat2(i)=dat(i)/128.0 !### Why 128 ??
enddo
do i=nh+1,nfft
dat2(i)=0.
enddo
call four2a(c,nfft,1,-1,0)
fac=1.0/nfft
do i=1,nh
s(i)=fac*(real(c(i))**2 + aimag(c(i))**2)
enddo
return
end subroutine ps4
lib/qratest.f90
-54
View File
@@ -1,54 +0,0 @@
program qratest
parameter (NMAX=60*12000)
real dd(NMAX)
character arg*8,mycall*12,hiscall*12,hisgrid*6,decoded*22
character c*1
logical loop
nargs=iargc()
if(nargs.lt.1 .or. nargs.gt.4) then
print*,'Usage: qratest nfile [sync f0 fTol]'
go to 999
endif
call getarg(1,arg)
read(arg,*) nfile
loop=arg(1:1).eq.'+'
minsync0=-1
nfqso0=-1
ntol0=-1
if(nargs.gt.1) then
call getarg(2,arg)
read(arg,*) minsync0
call getarg(3,arg)
read(arg,*) nfqso0
call getarg(4,arg)
read(arg,*) ntol0
endif
ndepth=3
nft=99
open(60,file='qra64_data.bin',access='stream')
do ifile=1,999
read(60,end=999) dd,npts,nutc,nf1,nf2,nfqso,ntol,mode64,minsync,ndepth, &
mycall,hiscall,hisgrid
if(ifile.lt.nfile) cycle
if(minsync0.ne.-1) minsync=minsync0
if(nfqso0.ne.-1) nfqso=nfqso0
if(ntol0.ne.-1) ntol=ntol0
call qra64a(dd,npts,nutc,nf1,nf2,nfqso,ntol,mode64,minsync,ndepth, &
mycall,hiscall,hisgrid,sync,nsnr,dtx,nfreq,decoded,nft)
c='a'
if(mode64.eq.2) c='b'
if(mode64.eq.4) c='c'
if(mode64.eq.8) c='d'
if(mode64.eq.16) c='e'
write(*,1000) ifile,c,nutc,nsnr,dtx,nfreq,decoded,nft-100,sync-3.4
1000 format(i4,1x,a1,1x,i4.4,i4,f6.2,i5,1x,a22,i3,f6.2)
if(ifile.eq.nfile .and. (.not.loop)) exit
enddo
999 end program qratest
lib/rectify_msk.f90
-64
View File
@@ -1,64 +0,0 @@
subroutine rectify_msk(c,msg0,imsg,freq2)
parameter (NSPM=1404)
complex c(0:NSPM-1) !Received data
complex cmsg(0:NSPM-1) !Message waveform
complex c1(0:NSPM-1) !Rectified signal
complex c2(0:NSPM-1) !Integral of rectified signal
complex c3(0:2*NSPM-1) !FFT of rectified signal
complex cfac
character*22 msg0,msg,msgsent
integer i4tone(234)
ichk=0
msg=msg0
nsym=234
if(imsg.ge.0) then
ichk=10000+imsg
msg="<C1ALL C2ALL> 73"
nsym=35
endif
call genmsk(msg,ichk,msgsent,i4tone,itype) !Get tone sequence for msg
twopi=8.0*atan(1.0)
dt=1.0/12000.0
f0=1000.0
f1=2000.0
phi=0.
dphi=0.
k=-1
c2=0.
do j=1,nsym !Generate Tx waveform for msg
if(i4tone(j).eq.0) dphi=twopi*f0*dt
if(i4tone(j).eq.1) dphi=twopi*f1*dt
do i=1,6
k=k+1
phi=phi+dphi
cmsg(k)=cmplx(cos(phi),sin(phi))
c1(k)=conjg(cmsg(k))*c(k)
if(k.ge.1) c2(k)=c2(k-1) + c1(k)
enddo
enddo
c2(0)=c2(1)
pha=atan2(aimag(c2(NSPM-1)),real(c2(NSPM-1)))
cfac=cmplx(cos(pha),-sin(pha))
c1=cfac*c1
c2=cfac*c2
nfft=2*NSPM
c3(0:NSPM-1)=c2
c3(NSPM:nfft-1)=0.
df=12000.0/nfft
call four2a(c3,nfft,1,-1,1)
smax=0.
do i=0,nfft-1
f=i*df
if(i.gt.nfft/2) f=f-12000.0
s=1.e-10*(real(c3(i))**2 + aimag(c3(i))**2)
if(s.gt.smax) then
smax=s
freq2=1500.0 + f
endif
enddo
return
end subroutine rectify_msk
lib/sfrsd.f90
-121
View File
@@ -1,121 +0,0 @@
subroutine sfrsd(mrsym,mrprob,mr2sym,mr2prob,ntrials,correct,indexes, &
param,ntry)
integer mrsym(0:62),mrprob(0:62),mr2sym(0:62),mr2prob(0:62)
integer correct(0:62),indexes(0:62),probs(0:62),thresh0(0:62)
integer rxdat(0:62),rxdat2(0:62),rxprob(0:62),rxprob2(0:62)
integer workdat(0:62),era_pos(0:50)
integer perr(0:7,0:7)
integer param(0:7)
real ratio0(0:62)
call init_rs_int()
do i=0,62
rxdat(i)=mrsym(62-i)
rxdat2(i)=mr2sym(62-i)
rxprob(i)=mrprob(62-i)
rxprob2(i)=mr2prob(62-i)
enddo
do i=0,62
indexes(i)=i
probs(i)=rxprob(i)
enddo
do ip=1,62
do k=0,63-ip
if(probs(k).lt.probs(k+1)) then
ntmp=probs(k)
probs(k)=probs(k+1)
probs(k+1)=ntmp
ntmp=indexes(k)
indexes(k)=indexes(k+1)
indexes(k+1)=ntmp
endif
enddo
enddo
era_pos=0
numera=0
workdat=rxdat
call decode_rs_int()
if(nerr.ge.0) then
correct=workdat
param=0
return
endif
call random_seed()
ncandidates=0
nsum=0
do i=0,62
nsum=nsum+rxprob(i)
j=indexes(62-i)
ratio0(i)=float(rxprob2(j))/(float(rxprob(j))+0.01)
ii=int(7.999*ratio0(i))
jj=(62-i)/8
thresh0(i)=nint(1.3*perr(jj,ii))
enddo
if(nsum.eq.0) return
do k=0,ntrials
era_pos=0
workdat=rxdat
numera=0
do i=0,62
j=indexes(62-i)
thresh=thresh0(i)
ir=rand()
if(...) then
era_pos(numera)=j
numera=numera+1
endif
enddo
call decode_rs_int()
if(nerr.ge.0) then
ncandidates=ncandidates+1
nhard=0
nsoft=0
nsofter=0
do i=0,62
if(workdat(i).ne.rxdat(i)) then
nhard=nhard+1
nsofter=nsofter+rxprob(i)
if(workdat(i).ne.rxdat2(i)) nsoft=nsoft+rxprob(i)
else
nsofter=nsofter-rxprob(i)
endif
enddo
nsoft=63*nsoft/nsum
nsofter=63*nsofter/nsum
ntotal=nsoft+nhard
if(ntotal.lt.ntotal_min) then
nsoft_min=nsoft
nhard_min=nhard
nsofter_min=nsofter
ntotal_min=ntotal
correct=workdat
nera_best=numera
ntry=k
endif
if(ntotal_min.lt.72 .and. nhard_min.lt.42) exit
endif
if(k.eq.ntrials-1) ntry=k+1
enddo
if(ntotal_min.ge.76 .or. nhard.ge.44) nhard_min=-1
param(0)=ncandidates
param(1)=nhard_min
param(2)=nsoft_min
param(3)=nera_best
param(4)=nsofter_min
if(param(0).eq.0) param(2)=-1
return
end subroutine sfrsd
lib/sgran.c
-7
View File
@@ -1,7 +0,0 @@
#include "init_random_seed.h"
/* Fortran wrapper to seed the C library rand */
void sgran_(void)
{
init_random_seed();
}
lib/slasubs.f
-3396
View File
File diff suppressed because it is too large Load Diff
lib/slope.f90
-40
View File
@@ -1,40 +0,0 @@
subroutine slope(y,npts,xpk)
! Remove best-fit slope from data in y(i). When fitting the straight line,
! ignore the peak around xpk +/- 4 bins
real y(npts)
sumw=0.
sumx=0.
sumy=0.
sumx2=0.
sumxy=0.
sumy2=0.
do i=1,npts
if(abs(i-xpk).gt.4.0) then
sumw=sumw + 1.0
x=i
sumx=sumx + x
sumy=sumy + y(i)
sumx2=sumx2 + x*x
sumxy=sumxy + x*y(i)
sumy2=sumy2 + y(i)**2
endif
enddo
delta=sumw*sumx2 - sumx**2
a=(sumx2*sumy - sumx*sumxy) / delta
b=(sumw*sumxy - sumx*sumy) / delta
sq=0.
do i=1,npts
y(i)=y(i)-(a + b*i)
if(abs(i-xpk).gt.4.0) sq=sq + y(i)**2
enddo
rms=sqrt(sq/(sumw-4.0))
y=y/rms
return
end subroutine slope
lib/ss.bat
-1
View File
@@ -1 +0,0 @@
svn status | grep -v "?"
lib/sun.f90
-88
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 !UT!in hours
real RA,Dec !RA and Dec of sun
! NB: Double caps here are single caps in the writeup.
! 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
! 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
!Ecliptic coords of sun (geocentric)
real xs,ys
!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/
! 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
! 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)
! Ecliptic coordinates of sun (rectangular):
xs = r * cos(lonsun/rad)
ys = r * sin(lonsun/rad)
! Equatorial coordinates of sun (rectangular):
xe = xs
ye = ys * cos(ecl/rad)
ze = ys * sin(ecl/rad)
! 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 subroutine sun
lib/t6.f90
-56
View File
@@ -1,56 +0,0 @@
program t6
parameter (MAXFFT=1404)
complex c(0:MAXFFT-1)
real s(0:MAXFFT-1)
m1=45
m2=67
m3=89
nsym=3*11 + m1 + m2 + m3
nfft=6*nsym
nh=nfft/2
best=9999.
! do m1=22,67
! do m2=37,97
do m1=30,67
do m2=26,100
m3=201-m2-m1
if(m3.lt.13) cycle
c=0.
n1=6*(11+m1)
n2=n1+6*(11+m2)
c(1:66)=1.
c(1+n1:66+n1)=1.
c(1+n2:66+n2)=1.
call four2a(c,nfft,1,-1,1) !c2c FFT
df=12000.0/nfft
smax=0.
do i=0,nfft-1
s(i)=real(c(i))**2 + aimag(c(i))**2
if(i.ne.0) smax=max(s(i),smax)
enddo
sidelobe=db(smax/s(0))
if(sidelobe.lt.best) then
write(*,1000) m1,m2,m3,sidelobe
1000 format(3i5,f8.2)
best=sidelobe
s=s/s(0)
rewind 13
do j=0,nfft-1
i=mod(j+nh,nfft)
f=i*df
if(i.gt.nh) f=f-12000.0
write(13,1020) f,s(i)
1020 format(2f12.4)
enddo
endif
enddo
enddo
end program t6
lib/testfast9.f90
-42
View File
@@ -1,42 +0,0 @@
program testfast9
parameter (NMAX=359424)
integer*2 id2(NMAX)
integer narg(0:11)
character*80 line(100)
character submode*1,infile*80
nargs=iargc()
if(nargs.ne.2) then
print*,'Usage: testfast9 submode infile'
print*,'Example: testfast9 E /data/VE1SKY/K1JT/JT9E/150806_123300.wav'
go to 999
endif
call getarg(1,submode)
call getarg(2,infile)
open(10,file=infile,access='stream',status='old')
read(10) id2(1:22) !Skip 44 header bytes
npts=NMAX
read(10,end=1) id2(1:npts) !Read the raw data
1 i1=index(infile,'.wav')
read(infile(i1-6:i1-1),*) narg(0)
narg(1)=NMAX
n=ichar(submode)
narg(2)=n-ichar('A')
if(n.ge.97 .and. n.le.104) narg(2)=n-ichar('a')
narg(3)=1
narg(4)=0
narg(5)=0
narg(6)=0
narg(7)=29951
narg(8)=1
narg(9)=102
narg(10)=700
narg(11)=500
call fast9(id2,narg,line)
print*,line(1)
999 end program testfast9
lib/tmoonsub.c
-518
View File
@@ -1,518 +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;
(void)day;
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;
(void)r;
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;
(void)r;
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);
}
lib/tweak1.f90
-23
View File
@@ -1,23 +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
complex w4
data twopi/0.d0/
save twopi
if(twopi.eq.0.d0) twopi=8.d0*atan(1.d0)
w=1.d0
dphi=twopi*f0/12000.d0
wstep=cmplx(cos(dphi),sin(dphi))
do i=1,jz
w=w*wstep
w4=w
cb(i)=w4*ca(i)
enddo
return
end subroutine tweak1
lib/update_hasharray.f90
-25
View File
@@ -1,25 +0,0 @@
subroutine update_hasharray(recent_calls,nrecent,nhasharray)
character*12 recent_calls(nrecent)
character*22 hashmsg
integer nhasharray(nrecent,nrecent)
nhasharray=-1
do i=1,nrecent
do j=i+1,nrecent
if( recent_calls(i)(1:1) .ne. ' ' .and. recent_calls(j)(1:1) .ne. ' ' ) then
hashmsg=trim(recent_calls(i))//' '//trim(recent_calls(j))
call fmtmsg(hashmsg,iz)
call hash(hashmsg,22,ihash)
ihash=iand(ihash,4095)
nhasharray(i,j)=ihash
hashmsg=trim(recent_calls(j))//' '//trim(recent_calls(i))
call fmtmsg(hashmsg,iz)
call hash(hashmsg,22,ihash)
ihash=iand(ihash,4095)
nhasharray(j,i)=ihash
endif
enddo
enddo
end subroutine update_hasharray
lib/update_recent_calls.f90
-19
View File
@@ -1,19 +0,0 @@
subroutine update_recent_calls(call,calls_hrd,nsize)
character*12 call,calls_hrd(nsize)
new=1
do ic=1,nsize
if( calls_hrd(ic).eq.call ) then
new=0
endif
enddo
if( new.eq.1 ) then
do ic=nsize-1,1,-1
calls_hrd(ic+1)(1:12)=calls_hrd(ic)(1:12)
enddo
calls_hrd(1)(1:12)=call(1:12)
endif
return
end subroutine update_recent_calls
lib/vit213.c
-221
View File
@@ -1,221 +0,0 @@
/* Viterbi decoder for arbitrary convolutional code
* viterbi27 and viterbi37 for the r=1/2 and r=1/3 K=7 codes are faster
* Copyright 1999 Phil Karn, KA9Q
* Modifications by Joe Taylor, K1JT
* May be used under the terms of the GNU Public License
*/
#include <stdlib.h>
/* Select code here */
#define V213
#ifdef V213
#define K 13 /* Constraint length */
#define N 2 /* Number of symbols per data bit */
#define Polys Poly213 /* Select polynomials here */
#endif
/* Rate 1/2 codes */
unsigned int Poly213[] = {012767,016461}; /* k = 13 */
#include <memory.h>
#define LONGBITS 32
#define LOGLONGBITS 5
#undef max
#define max(x,y) ((x) > (y) ? (x) : (y))
#define D (1 << max(0,K-LOGLONGBITS-1))
#define MAXNBITS 200 /* Maximum frame size (user bits) */
extern unsigned char Partab[]; /* Parity lookup table */
int Syms[1 << K];
int parity(int x)
{
x ^= (x >> 16);
x ^= (x >> 8);
return Partab[x & 0xff];
}
/* Convolutionally encode data into binary symbols */
int enc213(unsigned char symbols[], unsigned char data[],
unsigned int nbytes, unsigned int startstate,
unsigned int endstate)
{
unsigned int i,j,k;
int l,n=-1;
unsigned int encstate = startstate;
for(k=0; k<nbytes; k++) {
for(l=7;l>=0;l--){
encstate = (encstate + encstate) + ((data[k] >> l) & 1);
for(j=0;j<N;j++) {
n=n+1;
symbols[n] = parity(encstate & Polys[j]);
}
}
}
// Flush out with zero tail. (No need, if tail-biting code.)
for(i=0; i<K-1;i++){
encstate = (encstate << 1) | ((endstate >> i) & 1);
for(j=0;j<N;j++) {
n=n+1;
symbols[n] = parity(encstate & Polys[j]);
}
}
return 0;
}
/* Viterbi decoder */
int vit213(
int *metric, /* Final path metric (returned value) */
unsigned char *data, /* Decoded output data */
unsigned char *symbols, /* Raw deinterleaved input symbols */
unsigned int nbits, /* Number of output bits */
int mettab[2][256], /* Metric table, [sent sym][rx symbol] */
unsigned int startstate, /* Encoder starting state */
unsigned int endstate /* Encoder ending state */
){
int bitcnt = -(K-1);
int m0,m1;
int i,j,sym,ipp;
int mets[1 << N];
unsigned int paths[(MAXNBITS+K-1)*D];
unsigned int *pp,mask;
int cmetric[1 << (K-1)],nmetric[1 << (K-1)];
static int VDInit = 0;
memset(paths,0,sizeof(paths));
// Initialize on first time through:
if(!VDInit){
for(i=0;i<(1<<K);i++){
sym = 0;
for(j=0;j<N;j++)
sym = (sym << 1) + parity(i & Polys[j]);
Syms[i] = sym;
}
VDInit++;
}
// Keep only lower K-1 bits of specified startstate and endstate
startstate &= ~((1<<(K-1)) - 1);
endstate &= ~((1<<(K-1)) - 1);
/* Initialize starting metrics */
for(i=0;i< 1<<(K-1);i++)
cmetric[i] = -999999;
cmetric[startstate] = 0;
pp = paths;
ipp=0;
for(;;){ /* For each data bit */
/* Read input symbols and compute branch metrics */
for(i=0;i< 1<<N;i++){
mets[i] = 0;
for(j=0;j<N;j++){
mets[i] += mettab[(i >> (N-j-1)) & 1][symbols[j]];
}
}
symbols += N;
/* Run the add-compare-select operations */
mask = 1;
for(i=0;i< 1 << (K-1);i+=2){
int b1,b2;
b1 = mets[Syms[i]];
nmetric[i] = m0 = cmetric[i/2] + b1;
b2 = mets[Syms[i+1]];
b1 -= b2;
m1 = cmetric[(i/2) + (1<<(K-2))] + b2;
if(m1 > m0){
nmetric[i] = m1;
*pp |= mask;
}
m0 -= b1;
nmetric[i+1] = m0;
m1 += b1;
if(m1 > m0){
nmetric[i+1] = m1;
*pp |= mask << 1;
}
mask <<= 2;
if(mask == 0){
mask = 1;
pp++;
ipp++;
}
}
if(mask != 1){
pp++;
ipp++;
}
if(++bitcnt == (int)nbits){
*metric = nmetric[endstate];
break;
}
memcpy(cmetric,nmetric,sizeof(cmetric));
}
/* Chain back from terminal state to produce decoded data */
if(data == NULL)
return 0;/* Discard output */
memset(data,0,(nbits+7)/8); /* round up in case nbits % 8 != 0 */
for(i=nbits-1;i >= 0;i--){
// int a0,a1;
pp -= D;
ipp -= D;
m0=endstate >> LOGLONGBITS;
m1=1L << (endstate & (LONGBITS-1));
if(pp[m0] & m1) {
// a0=nmetric[endstate];
endstate |= (1 << (K-1));
// a1=nmetric[endstate];
data[i>>3] |= 0x80 >> (i&7);
// printf("B %d %d %d %d\n",*metric,i,a0,a1);
}
endstate >>= 1;
}
return 0;
}
// Wrapper for calling "encode" from Fortran:
void enc213_(
unsigned char data[], // User data, 8 bits per byte
int *nbits, // Number of user bits
unsigned char symbols[], // Encoded one-bit symbols, 8 per byte
int *nsymbols, // Number of symbols
int *kk, // K
int *nn) // N
{
int nbytes;
nbytes=(*nbits+7)/8; // Always encode multiple of 8 information bits
enc213(symbols,data,nbytes,0,0); // Do the encoding
*nsymbols=(*nbits+K-1)*N; // Return number of encoded symbols
*kk=K;
*nn=N;
}
// Wrapper for calling "viterbi" from Fortran:
void vit213_(
unsigned char symbols[], /* Raw deinterleaved input symbols */
unsigned int *Nbits, /* Number of decoded information bits */
int mettab[2][256], /* Metric table, [sent sym][rx symbol] */
unsigned char ddec[], /* Decoded output data */
int *Metric /* Final path metric (bigger is better) */
){
int metric;
vit213(&metric,ddec,symbols,*Nbits,mettab,0,0);
*Metric=metric;
}
lib/vit216.c
-219
View File
@@ -1,219 +0,0 @@
/* Viterbi decoder for arbitrary convolutional code
* viterbi27 and viterbi37 for the r=1/2 and r=1/3 K=7 codes are faster
* Copyright 1999 Phil Karn, KA9Q
* May be used under the terms of the GNU Public License
*/
/* Select code here */
#define V216
#ifdef V216
#define K 16 /* Constraint length */
#define N 2 /* Number of symbols per data bit */
#define Polys Poly216 /* Select polynomials here */
#endif
/* Rate 1/2 codes */
unsigned int Poly216[] = {0126723, 0152711}; /* k = 16 */
#include <stdlib.h>
#include <memory.h>
#define LONGBITS 32
#define LOGLONGBITS 5
#undef max
#define max(x,y) ((x) > (y) ? (x) : (y))
#define D (1 << max(0,K-LOGLONGBITS-1))
#define MAXNBITS 200 /* Maximum frame size (user bits) */
extern unsigned char Partab[]; /* Parity lookup table */
int Syms[1 << K];
int VDInit = 0;
int parity(int x)
{
x ^= (x >> 16);
x ^= (x >> 8);
return Partab[x & 0xff];
}
// Wrapper for calling "encode" from Fortran:
//void __stdcall ENCODE(
void enc216_(
unsigned char data[], // User data, 8 bits per byte
int *nbits, // Number of user bits
unsigned char symbols[], // Encoded one-bit symbols, 8 per byte
int *nsymbols, // Number of symbols
int *kk, // K
int *nn) // N
{
int nbytes;
nbytes=(*nbits+7)/8; // Always encode multiple of 8 information bits
enc216(symbols,data,nbytes,0,0); // Do the encoding
*nsymbols=(*nbits+K-1)*N; // Return number of encoded symbols
*kk=K;
*nn=N;
}
/* Convolutionally encode data into binary symbols */
enc216(unsigned char symbols[], unsigned char data[],
unsigned int nbytes, unsigned int startstate,
unsigned int endstate)
{
int i,j,k,n=-1;
unsigned int encstate = startstate;
for(k=0; k<nbytes; k++) {
for(i=7;i>=0;i--){
encstate = (encstate + encstate) + ((data[k] >> i) & 1);
for(j=0;j<N;j++) {
n=n+1;
symbols[n] = parity(encstate & Polys[j]);
}
}
}
// Flush out with zero tail. (No need, if tail-biting code.)
for(i=0; i<K-1;i++){
encstate = (encstate << 1) | ((endstate >> i) & 1);
for(j=0;j<N;j++) {
n=n+1;
symbols[n] = parity(encstate & Polys[j]);
}
}
return 0;
}
// Wrapper for calling "viterbi" from Fortran:
//void __stdcall VITERBI(
void vit216_(
unsigned char symbols[], /* Raw deinterleaved input symbols */
unsigned int *Nbits, /* Number of decoded information bits */
int mettab[2][256], /* Metric table, [sent sym][rx symbol] */
unsigned char ddec[], /* Decoded output data */
long *Metric /* Final path metric (bigger is better) */
){
long metric;
vit216(&metric,ddec,symbols,*Nbits,mettab,0,0);
*Metric=metric;
}
/* Viterbi decoder */
int vit216(
long *metric, /* Final path metric (returned value) */
unsigned char *data, /* Decoded output data */
unsigned char *symbols, /* Raw deinterleaved input symbols */
unsigned int nbits, /* Number of output bits */
int mettab[2][256], /* Metric table, [sent sym][rx symbol] */
unsigned int startstate, /* Encoder starting state */
unsigned int endstate /* Encoder ending state */
){
int bitcnt = -(K-1);
long m0,m1;
int i,j,sym,ipp;
int mets[1 << N];
unsigned long paths[(MAXNBITS+K-1)*D];
unsigned long *pp,mask;
long cmetric[1 << (K-1)],nmetric[1 << (K-1)];
memset(paths,0,sizeof(paths));
// Initialize on first time through:
if(!VDInit){
for(i=0;i<(1<<K);i++){
sym = 0;
for(j=0;j<N;j++)
sym = (sym << 1) + parity(i & Polys[j]);
Syms[i] = sym;
}
VDInit++;
}
// Keep only lower K-1 bits of specified startstate and endstate
startstate &= ~((1<<(K-1)) - 1);
endstate &= ~((1<<(K-1)) - 1);
/* Initialize starting metrics */
for(i=0;i< 1<<(K-1);i++)
cmetric[i] = -999999;
cmetric[startstate] = 0;
pp = paths;
ipp=0;
for(;;){ /* For each data bit */
/* Read input symbols and compute branch metrics */
for(i=0;i< 1<<N;i++){
mets[i] = 0;
for(j=0;j<N;j++){
mets[i] += mettab[(i >> (N-j-1)) & 1][symbols[j]];
}
}
symbols += N;
/* Run the add-compare-select operations */
mask = 1;
for(i=0;i< 1 << (K-1);i+=2){
int b1,b2;
b1 = mets[Syms[i]];
nmetric[i] = m0 = cmetric[i/2] + b1;
b2 = mets[Syms[i+1]];
b1 -= b2;
m1 = cmetric[(i/2) + (1<<(K-2))] + b2;
if(m1 > m0){
nmetric[i] = m1;
*pp |= mask;
}
m0 -= b1;
nmetric[i+1] = m0;
m1 += b1;
if(m1 > m0){
nmetric[i+1] = m1;
*pp |= mask << 1;
}
mask <<= 2;
if(mask == 0){
mask = 1;
pp++;
ipp++;
}
}
if(mask != 1){
pp++;
ipp++;
}
if(++bitcnt == nbits){
*metric = nmetric[endstate];
break;
}
memcpy(cmetric,nmetric,sizeof(cmetric));
}
/* Chain back from terminal state to produce decoded data */
if(data == NULL)
return 0;/* Discard output */
memset(data,0,(nbits+7)/8); /* round up in case nbits % 8 != 0 */
for(i=nbits-1;i >= 0;i--){
// int a0,a1;
pp -= D;
ipp -= D;
m0=endstate >> LOGLONGBITS;
m1=1L << (endstate & (LONGBITS-1));
if(pp[m0] & m1) {
// a0=nmetric[endstate];
endstate |= (1 << (K-1));
// a1=nmetric[endstate];
data[i>>3] |= 0x80 >> (i&7);
// printf("B %d %d %d %d\n",*metric,i,a0,a1);
}
endstate >>= 1;
}
return 0;
}
lib/wav11.f90
-27
View File
@@ -1,27 +0,0 @@
subroutine wav11(d2,npts,dd)
! Convert i*2 data sampled at 12000 Hz to r*4 sampled at 11025 Hz.
parameter (NZ11=60*11025,NZ12=60*12000)
parameter (NFFT1=64*12000,NFFT2=64*11025)
integer*2 d2(NZ12)
real*4 dd(NZ11)
real x(NFFT2)
complex cx(0:NFFT1/2)
equivalence (x,cx)
save x,cx
jz=min(NZ12,npts)
x(1:jz)=d2(1:jz)
x(jz+1:)=0.0
call four2a(cx,nfft1,1,-1,0) !Forward FFT, r2c
df=12000.0/NFFT1
ia=5000.0/df
cx(ia:)=0.0
call four2a(cx,nfft2,1,1,-1) !Inverse FFT, c2r
npts=jz*11025.0/12000.0
fac=1.e-6
dd(1:npts)=fac*x(1:npts)
return
end subroutine wav11
lib/wqdecode.f90
-338
View File
@@ -1,338 +0,0 @@
!-------------------------------------------------------------------------------
!
! This file is part of the WSPR application, Weak Signal Propagation Reporter
!
! File Name: wqdecode.f90
! Description:
!
! Copyright (C) 2001-2014 Joseph Taylor, K1JT
! License: GPL-3
!
! This program is free software; you can redistribute it and/or modify it under
! the terms of the GNU General Public License as published by the Free Software
! Foundation; either version 3 of the License, or (at your option) any later
! version.
!
! This program is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
! details.
!
! You should have received a copy of the GNU General Public License along with
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
! Street, Fifth Floor, Boston, MA 02110-1301, USA.
!
!-------------------------------------------------------------------------------
subroutine wqdecode(data0,message,ntype)
parameter (N15=32768)
integer*1 data0(11)
character*22 message
character*12 callsign
character*3 cdbm
character grid4*4,grid6*6
logical first
character*12 dcall(0:N15-1)
data first/.true./
save first,dcall
! May want to have a timeout (say, one hour?) on calls fetched
! from the hash table.
if(first) then
dcall=' '
first=.false.
endif
message=' '
call unpack50(data0,n1,n2)
! print*,data0,n1,n2
call unpackcall(n1,callsign)
i1=index(callsign,' ')
call unpackgrid(n2/128,grid4)
ntype=iand(n2,127) -64
! Standard WSPR message (types 0 3 7 10 13 17 ... 60)
if(ntype.ge.0 .and. ntype.le.62) then
nu=mod(ntype,10)
if(nu.eq.0 .or. nu.eq.3 .or. nu.eq.7) then
write(cdbm,'(i3)'),ntype
if(cdbm(1:1).eq.' ') cdbm=cdbm(2:)
if(cdbm(1:1).eq.' ') cdbm=cdbm(2:)
message=callsign(1:i1)//grid4//' '//cdbm
call hash(callsign,i1-1,ih)
dcall(ih)=callsign(:i1)
else
nadd=nu
if(nu.gt.3) nadd=nu-3
if(nu.gt.7) nadd=nu-7
ng=n2/128 + 32768*(nadd-1)
call unpackpfx(ng,callsign)
ndbm=ntype-nadd
write(cdbm,'(i3)'),ndbm
if(cdbm(1:1).eq.' ') cdbm=cdbm(2:)
if(cdbm(1:1).eq.' ') cdbm=cdbm(2:)
i2=index(callsign,' ')
message=callsign(:i2)//cdbm
call hash(callsign,i2-1,ih)
dcall(ih)=callsign(:i2)
endif
else if(ntype.lt.0) then
ndbm=-(ntype+1)
grid6=callsign(6:6)//callsign(1:5)
ih=(n2-ntype-64)/128
callsign=dcall(ih)
write(cdbm,'(i3)'),ndbm
if(cdbm(1:1).eq.' ') cdbm=cdbm(2:)
if(cdbm(1:1).eq.' ') cdbm=cdbm(2:)
i2=index(callsign,' ')
if(dcall(ih)(1:1).ne.' ') then
message='<'//callsign(:i2-1)//'> '//grid6//' '//cdbm
else
message='<...> '//grid6//' '//cdbm
endif
endif
return
end subroutine wqdecode
!-------------------------------------------------------------------------------
!
! This file is part of the WSPR application, Weak Signal Propagation Reporter
!
! File Name: unpack50.f90
! Description:
!
! Copyright (C) 2001-2014 Joseph Taylor, K1JT
! License: GPL-3
!
! This program is free software; you can redistribute it and/or modify it under
! the terms of the GNU General Public License as published by the Free Software
! Foundation; either version 3 of the License, or (at your option) any later
! version.
!
! This program is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
! details.
!
! You should have received a copy of the GNU General Public License along with
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
! Street, Fifth Floor, Boston, MA 02110-1301, USA.
!
!-------------------------------------------------------------------------------
subroutine unpack50(dat,n1,n2)
integer*1 dat(11)
i=dat(1)
i4=iand(i,255)
n1=ishft(i4,20)
i=dat(2)
i4=iand(i,255)
n1=n1 + ishft(i4,12)
i=dat(3)
i4=iand(i,255)
n1=n1 + ishft(i4,4)
i=dat(4)
i4=iand(i,255)
n1=n1 + iand(ishft(i4,-4),15)
n2=ishft(iand(i4,15),18)
i=dat(5)
i4=iand(i,255)
n2=n2 + ishft(i4,10)
i=dat(6)
i4=iand(i,255)
n2=n2 + ishft(i4,2)
i=dat(7)
i4=iand(i,255)
n2=n2 + iand(ishft(i4,-6),3)
return
end subroutine unpack50
!-------------------------------------------------------------------------------
!
! This file is part of the WSPR application, Weak Signal Propagation Reporter
!
! File Name: unpackcall.f90
! Description:
!
! Copyright (C) 2001-2014 Joseph Taylor, K1JT
! License: GPL-3
!
! This program is free software; you can redistribute it and/or modify it under
! the terms of the GNU General Public License as published by the Free Software
! Foundation; either version 3 of the License, or (at your option) any later
! version.
!
! This program is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
! details.
!
! You should have received a copy of the GNU General Public License along with
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
! Street, Fifth Floor, Boston, MA 02110-1301, USA.
!
!-------------------------------------------------------------------------------
subroutine unpackcall(ncall,word)
character word*12,c*37
data c/'0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ '/
n=ncall
word='......'
if(n.ge.262177560) go to 999 !Plain text message ...
i=mod(n,27)+11
word(6:6)=c(i:i)
n=n/27
i=mod(n,27)+11
word(5:5)=c(i:i)
n=n/27
i=mod(n,27)+11
word(4:4)=c(i:i)
n=n/27
i=mod(n,10)+1
word(3:3)=c(i:i)
n=n/10
i=mod(n,36)+1
word(2:2)=c(i:i)
n=n/36
i=n+1
word(1:1)=c(i:i)
do i=1,4
if(word(i:i).ne.' ') go to 10
enddo
go to 999
10 word=word(i:)
999 if(word(1:3).eq.'3D0') word='3DA0'//word(4:)
return
end subroutine unpackcall
!-------------------------------------------------------------------------------
!
! This file is part of the WSPR application, Weak Signal Propagation Reporter
!
! File Name: unpackgrid.f90
! Description:
!
! Copyright (C) 2001-2014 Joseph Taylor, K1JT
! License: GPL-3
!
! This program is free software; you can redistribute it and/or modify it under
! the terms of the GNU General Public License as published by the Free Software
! Foundation; either version 3 of the License, or (at your option) any later
! version.
!
! This program is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
! details.
!
! You should have received a copy of the GNU General Public License along with
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
! Street, Fifth Floor, Boston, MA 02110-1301, USA.
!
!-------------------------------------------------------------------------------
subroutine unpackgrid(ng,grid)
parameter (NGBASE=180*180)
character grid*4,grid6*6,digit*10
data digit/'0123456789'/
grid=' '
if(ng.ge.32400) go to 10
dlat=mod(ng,180)-90
dlong=(ng/180)*2 - 180 + 2
call deg2grid(dlong,dlat,grid6)
grid=grid6(1:4) !XXX explicitly truncate this -db
go to 100
10 n=ng-NGBASE-1
if(n.ge.1 .and.n.le.30) then
grid(1:1)='-'
grid(2:2)=char(48+n/10)
grid(3:3)=char(48+mod(n,10))
else if(n.ge.31 .and.n.le.60) then
n=n-30
grid(1:2)='R-'
grid(3:3)=char(48+n/10)
grid(4:4)=char(48+mod(n,10))
else if(n.eq.61) then
grid='RO'
else if(n.eq.62) then
grid='RRR'
else if(n.eq.63) then
grid='73'
endif
100 return
end subroutine unpackgrid
!-------------------------------------------------------------------------------
!
! This file is part of the WSPR application, Weak Signal Propagation Reporter
!
! File Name: unpackpfx.f90
! Description:
!
! Copyright (C) 2001-2014 Joseph Taylor, K1JT
! License: GPL-3
!
! This program is free software; you can redistribute it and/or modify it under
! the terms of the GNU General Public License as published by the Free Software
! Foundation; either version 3 of the License, or (at your option) any later
! version.
!
! This program is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
! details.
!
! You should have received a copy of the GNU General Public License along with
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
! Street, Fifth Floor, Boston, MA 02110-1301, USA.
!
!-------------------------------------------------------------------------------
subroutine unpackpfx(ng,call1)
character*12 call1
character*3 pfx
if(ng.lt.60000) then
! Add-on prefix of 1 to 3 characters
n=ng
do i=3,1,-1
nc=mod(n,37)
if(nc.ge.0 .and. nc.le.9) then
pfx(i:i)=char(nc+48)
else if(nc.ge.10 .and. nc.le.35) then
pfx(i:i)=char(nc+55)
else
pfx(i:i)=' '
endif
n=n/37
enddo
call1=pfx//'/'//call1
if(call1(1:1).eq.' ') call1=call1(2:)
if(call1(1:1).eq.' ') call1=call1(2:)
else
! Add-on suffix, one or teo characters
i1=index(call1,' ')
nc=ng-60000
if(nc.ge.0 .and. nc.le.9) then
call1=call1(:i1-1)//'/'//char(nc+48)
else if(nc.ge.10 .and. nc.le.35) then
call1=call1(:i1-1)//'/'//char(nc+55)
else if(nc.ge.36 .and. nc.le.125) then
nc1=(nc-26)/10
nc2=mod(nc-26,10)
call1=call1(:i1-1)//'/'//char(nc1+48)//char(nc2+48)
endif
endif
return
end subroutine unpackpfx
lib/wqencode.f90
-66
View File
@@ -1,66 +0,0 @@
subroutine wqencode(msg,ntype,data0)
! Parse and encode a WSPR message.
use packjt
parameter (MASK15=32767)
character*22 msg
character*12 call1,call2
character grid4*4
logical lbad1,lbad2
integer*1 data0(11)
integer nu(0:9)
data nu/0,-1,1,0,-1,2,1,0,-1,1/
! Standard WSPR message (types 0 3 7 10 13 17 ... 60)
i1=index(msg,' ')
i2=index(msg,'/')
i3=index(msg,'<')
call1=msg(:i1-1)
if(i1.lt.3 .or. i1.gt.7 .or. i2.gt.0 .or. i3.gt.0) go to 10
grid4=msg(i1+1:i1+4)
call packcall(call1,n1,lbad1)
call packgrid(grid4,ng,lbad2)
if(lbad1 .or. lbad2) go to 10
ndbm=0
read(msg(i1+5:),*) ndbm
if(ndbm.lt.0) ndbm=0
if(ndbm.gt.60) ndbm=60
ndbm=ndbm+nu(mod(ndbm,10))
n2=128*ng + (ndbm+64)
call pack50(n1,n2,data0)
ntype=ndbm
go to 900
10 if(i2.ge.2 .and. i3.lt.1) then
call packpfx(call1,n1,ng,nadd)
ndbm=0
read(msg(i1+1:),*) ndbm
if(ndbm.lt.0) ndbm=0
if(ndbm.gt.60) ndbm=60
ndbm=ndbm+nu(mod(ndbm,10))
ntype=ndbm + 1 + nadd
n2=128*ng + ntype + 64
call pack50(n1,n2,data0)
else if(i3.eq.1) then
i4=index(msg,'>')
call1=msg(2:i4-1)
call hash(call1,i4-2,ih)
i5=index(trim(msg(i1+1:)),' ')
! Convert grid to valid callsign format - first character moved to end
call2=msg(i1+2:i1+i5-1)//msg(i1+1:i1+1)//' '
call packcall(call2,n1,lbad1)
ndbm=0
read(msg(i1+i5+1:),*) ndbm
if(ndbm.lt.0) ndbm=0
if(ndbm.gt.60) ndbm=60
ndbm=ndbm+nu(mod(ndbm,10))
ntype=-(ndbm+1)
n2=128*ih + ntype + 64
call pack50(n1,n2,data0)
endif
go to 900
900 continue
return
end subroutine wqencode
lib/zplt.f90
-112
View File
@@ -1,112 +0,0 @@
subroutine zplt(z,iplt,sync,dtx,nfreq,flip,sync2,nplot,emedelay,dttol, &
nfqso,ntol)
real z(458,65)
real zz(458,65)
integer ij(2)
character*4 lab
call pctile(z,458*65,84,rms)
fac=0.05/rms
z=fac*z
dtq=0.114286
df=11025.0/(2.0*2520.0)
ia=nint((nfqso-ntol)/df) - 273
if(ia.lt.1) ia=1
ib=nint((nfqso+ntol)/df) - 273
if(ib.gt.458) ib=458
ja=(emedelay+0.8-dttol)/dtq
if(ja.lt.1) ja=1
jb=(emedelay+0.8+dttol)/dtq
if(jb.gt.65) jb=65
zz=0.
zz(ia:ib,ja:jb)=z(ia:ib,ja:jb)
zmin=minval(zz)
zmax=maxval(zz)
flip=1.0
if(abs(zmin).gt.abs(zmax)) flip=-1.0
ij=maxloc(zz)
if(flip.lt.0.0) ij=minloc(zz)
i0=ij(1)
j0=ij(2)
nfreq=nint((i0+273)*df)
dtx=j0*dtq-0.8
! write(69,3101) ia,ib,ja,jb,ij,dtx,nfreq
!3101 format(6i5,f8.2,i6)
ia=max(1,i0-72)
ib=min(458,i0+72)
sync=16.33*flip*(z(i0,j0) - 0.5*(z(ia,j0)+z(ib,j0)))
sync2=20.0*flip*z(i0,j0)
if(nplot.eq.0) go to 900
zmax=max(abs(zmin),abs(zmax),1.0)
zmin=-zmax
do j=1,65
write(61,1100) j*dtq-0.8,z(i0,j)
1100 format(2f10.3)
enddo
do i=1,458
write(62,1100) (i+273)*df,flip*z(i,j0)
enddo
xx=1.5
yy=7.5 - 3.0*iplt
width=6.0
height=2.0
IP=458
JP=65
imax=IP
jmax=JP
if(iplt.eq.0) then
call imopen("testjt4.ps")
call imfont("Helvetica",16)
call impalette("BlueRed.pal")
endif
call imr4mat_color(z,IP,JP,imax,jmax,zmin,zmax,xx,yy, &
width,height,1)
call imstring("Frequency (Hz)",xx+0.5*width,yy-0.5,2,0)
dy=0.1
do i=1,9
x=xx + 0.1*i*width
call imyline(x,yy,dy)
call imyline(x,yy+height,-dy)
enddo
do i=1,6
nf=(i-1)*200 + 600
write(lab,1020) nf
1020 format(i4)
x=xx + (i-1)*0.2*width
call imstring(lab,x,yy-0.25,2,0)
enddo
dx=0.1
do i=0,6
y=yy + height*(0.8+i)/(65.0*0.114286)
call imxline(xx,y,dx)
call imxline(xx+width,y,-dx)
enddo
do i=0,6,2
y=yy + height*(0.8+i)/(65.0*0.114286)
write(lab,1020) i
call imstring(lab(4:4),xx-0.15,y-0.08,2,0)
enddo
y=yy + height*(3.8)/(65.0*0.114286)
call imstring("DT", xx-0.5,y ,2,0)
call imstring("(s)",xx-0.5,y-0.25,2,0)
if(iplt.eq.2) call imclose
900 return
end subroutine zplt