diff --git a/libm65/CMakeLists.txt b/libm65/CMakeLists.txt index 77cd2454e..12b33e697 100644 --- a/libm65/CMakeLists.txt +++ b/libm65/CMakeLists.txt @@ -55,28 +55,46 @@ FortranCInterface_HEADER (FC.h MACRO_NAMESPACE "FC_" SYMBOL_NAMESPACE "FC_" set (FSRCS + afc65b.f90 + astro.f90 astro0.f90 astrosub.f90 + ccf2.f90 ccf65.f90 cgen65.f90 + chkhist.f90 + chkmsg.f90 + coord.f90 + dcoord.f90 decode0.f90 decode1a.f90 + decode65b.f90 deep65.f90 + deg2grid.f90 demod64a.f90 display.f90 + dot.f90 dpol.f90 + encode65.f90 extract.f90 + fchisq.f90 + fil6521.f90 + filbig.f90 four2a.f90 ftninit.f90 ftnquit.f90 gen65.f90 + geocentric.f90 getdphi.f90 + getpfx1.f90 + getpfx2.f90 graycode65.f90 iqcal.f90 iqfix.f90 jt65code.f90 map65a.f90 noisegen.f90 + pfxdump.f90 recvpkt.f90 rfile3a.f90 s3avg.f90 @@ -88,24 +106,7 @@ set (FSRCS tm2.f90 zplot.f90 - afc65b.f - astro.f - ccf2.f - chkhist.f - chkmsg.f - coord.f - dcoord.f - decode65b.f - deg2grid.f - dot.f - encode65.f f77_wisdom.f - fchisq.f - fil6521.f - filbig.f - geocentric.f - getpfx1.f - getpfx2.f graycode.f grid2deg.f grid2k.f @@ -121,7 +122,6 @@ set (FSRCS packmsg.f packtext.f pctile.f - pfxdump.f set.f setup65.f sort.f diff --git a/libm65/afc65b.f b/libm65/afc65b.f deleted file mode 100644 index 1b294a3cf..000000000 --- a/libm65/afc65b.f +++ /dev/null @@ -1,71 +0,0 @@ - subroutine afc65b(cx,cy,npts,nfast,fsample,nflip,ipol,xpol, - + ndphi,iloop,a,ccfbest,dtbest) - - logical xpol - complex cx(npts) - complex cy(npts) - real a(5),deltaa(5) - - a(1)=0. - a(2)=0. - a(3)=0. - a(4)=45.0*(ipol-1.0) - deltaa(1)=2.0 - deltaa(2)=2.0 - deltaa(3)=2.0 - deltaa(4)=22.5 - deltaa(5)=0.05 - nterms=3 - if(xpol) nterms=4 - -! Don't fit polarization when solving for dphi - if(ndphi.ne.0) nterms=3 - -! Start the iteration - chisqr=0. - chisqr0=1.e6 - do iter=1,3 !One iteration is enough? - do j=1,nterms - chisq1=fchisq(cx,cy,npts,nfast,fsample,nflip,a,ccfmax,dtmax) - fn=0. - delta=deltaa(j) - 10 a(j)=a(j)+delta - chisq2=fchisq(cx,cy,npts,nfast,fsample,nflip,a,ccfmax,dtmax) - if(chisq2.eq.chisq1) go to 10 - if(chisq2.gt.chisq1) then - delta=-delta !Reverse direction - a(j)=a(j)+delta - tmp=chisq1 - chisq1=chisq2 - chisq2=tmp - endif - 20 fn=fn+1.0 - a(j)=a(j)+delta - chisq3=fchisq(cx,cy,npts,nfast,fsample,nflip,a,ccfmax,dtmax) - if(chisq3.lt.chisq2) then - chisq1=chisq2 - chisq2=chisq3 - go to 20 - endif - -C Find minimum of parabola defined by last three points - delta=delta*(1./(1.+(chisq1-chisq2)/(chisq3-chisq2))+0.5) - a(j)=a(j)-delta - deltaa(j)=deltaa(j)*fn/3. - enddo - chisqr=fchisq(cx,cy,npts,nfast,fsample,nflip,a,ccfmax,dtmax) - if(chisqr/chisqr0.gt.0.9999) go to 30 - chisqr0=chisqr - enddo - - 30 ccfbest=ccfmax * (1378.125/fsample)**2 - dtbest=dtmax - - if(a(4).lt.0.0) a(4)=a(4)+180.0 - if(a(4).ge.180.0) a(4)=a(4)-180.0 - if(nint(a(4)).eq.180) a(4)=0. - ipol=nint(a(4)/45.0) + 1 - if(ipol.gt.4) ipol=ipol-4 - - return - end diff --git a/libm65/afc65b.f90 b/libm65/afc65b.f90 new file mode 100644 index 000000000..eb7e74d51 --- /dev/null +++ b/libm65/afc65b.f90 @@ -0,0 +1,71 @@ +subroutine afc65b(cx,cy,npts,nfast,fsample,nflip,ipol,xpol,ndphi,iloop, & + a,ccfbest,dtbest) + + logical xpol + complex cx(npts) + complex cy(npts) + real a(5),deltaa(5) + + a(1)=0. + a(2)=0. + a(3)=0. + a(4)=45.0*(ipol-1.0) + deltaa(1)=2.0 + deltaa(2)=2.0 + deltaa(3)=2.0 + deltaa(4)=22.5 + deltaa(5)=0.05 + nterms=3 + if(xpol) nterms=4 + +! Don't fit polarization when solving for dphi + if(ndphi.ne.0) nterms=3 + +! Start the iteration + chisqr=0. + chisqr0=1.e6 + do iter=1,3 !One iteration is enough? + do j=1,nterms + chisq1=fchisq(cx,cy,npts,nfast,fsample,nflip,a,ccfmax,dtmax) + fn=0. + delta=deltaa(j) +10 a(j)=a(j)+delta + chisq2=fchisq(cx,cy,npts,nfast,fsample,nflip,a,ccfmax,dtmax) + if(chisq2.eq.chisq1) go to 10 + if(chisq2.gt.chisq1) then + delta=-delta !Reverse direction + a(j)=a(j)+delta + tmp=chisq1 + chisq1=chisq2 + chisq2=tmp + endif +20 fn=fn+1.0 + a(j)=a(j)+delta + chisq3=fchisq(cx,cy,npts,nfast,fsample,nflip,a,ccfmax,dtmax) + 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. + enddo + chisqr=fchisq(cx,cy,npts,nfast,fsample,nflip,a,ccfmax,dtmax) + if(chisqr/chisqr0.gt.0.9999) go to 30 + chisqr0=chisqr + enddo + +30 ccfbest=ccfmax * (1378.125/fsample)**2 + dtbest=dtmax + + if(a(4).lt.0.0) a(4)=a(4)+180.0 + if(a(4).ge.180.0) a(4)=a(4)-180.0 + if(nint(a(4)).eq.180) a(4)=0. + ipol=nint(a(4)/45.0) + 1 + if(ipol.gt.4) ipol=ipol-4 + + return +end subroutine afc65b diff --git a/libm65/astro.f b/libm65/astro.f deleted file mode 100644 index a92b6ad12..000000000 --- a/libm65/astro.f +++ /dev/null @@ -1,109 +0,0 @@ - subroutine astro(nyear,month,nday,uth,nfreq,Mygrid, - + NStation,MoonDX,AzSun,ElSun,AzMoon0,ElMoon0, - + ntsky,doppler00,doppler,dbMoon,RAMoon,DecMoon,HA,Dgrd,sd, - + poloffset,xnr,day,lon,lat,LST) - -C Computes astronomical quantities for display and tracking. -C NB: may want to smooth the Tsky map to 10 degrees or so. - - character*6 MyGrid,HisGrid - real LST - real lat,lon - integer*2 nt144(180) - -! common/echo/xdop(2),techo,AzMoon,ElMoon,mjd - real xdop(2) - - data rad/57.2957795/ - data nt144/ - + 234, 246, 257, 267, 275, 280, 283, 286, 291, 298, - + 305, 313, 322, 331, 341, 351, 361, 369, 376, 381, - + 383, 382, 379, 374, 370, 366, 363, 361, 363, 368, - + 376, 388, 401, 415, 428, 440, 453, 467, 487, 512, - + 544, 579, 607, 618, 609, 588, 563, 539, 512, 482, - + 450, 422, 398, 379, 363, 349, 334, 319, 302, 282, - + 262, 242, 226, 213, 205, 200, 198, 197, 196, 197, - + 200, 202, 204, 205, 204, 203, 202, 201, 203, 206, - + 212, 218, 223, 227, 231, 236, 240, 243, 247, 257, - + 276, 301, 324, 339, 346, 344, 339, 331, 323, 316, - + 312, 310, 312, 317, 327, 341, 358, 375, 392, 407, - + 422, 437, 451, 466, 480, 494, 511, 530, 552, 579, - + 612, 653, 702, 768, 863,1008,1232,1557,1966,2385, - + 2719,2924,3018,3038,2986,2836,2570,2213,1823,1461, - + 1163, 939, 783, 677, 602, 543, 494, 452, 419, 392, - + 373, 360, 353, 350, 350, 350, 350, 350, 350, 348, - + 344, 337, 329, 319, 307, 295, 284, 276, 272, 272, - + 273, 274, 274, 271, 266, 260, 252, 245, 238, 231/ - save - - call grid2deg(MyGrid,elon,lat) - lon=-elon - call sun(nyear,month,nday,uth,lon,lat,RASun,DecSun,LST, - + AzSun,ElSun,mjd,day) - - freq=nfreq*1.e6 - if(nfreq.eq.2) freq=1.8e6 - if(nfreq.eq.4) freq=3.5e6 - - call MoonDop(nyear,month,nday,uth,lon,lat,RAMoon,DecMoon, - + LST,HA,AzMoon,ElMoon,vr,dist) - -C Compute spatial polarization offset - xx=sin(lat/rad)*cos(ElMoon/rad) - cos(lat/rad)* - + cos(AzMoon/rad)*sin(ElMoon/rad) - yy=cos(lat/rad)*sin(AzMoon/rad) - if(NStation.eq.1) poloffset1=rad*atan2(yy,xx) - if(NStation.eq.2) poloffset2=rad*atan2(yy,xx) - - techo=2.0 * dist/2.99792458e5 !Echo delay time - doppler=-freq*vr/2.99792458e5 !One-way Doppler - - call coord(0.,0.,-1.570796,1.161639,RAMoon/rad,DecMoon/rad,el,eb) - longecl_half=nint(rad*el/2.0) - if(longecl_half.lt.1 .or. longecl_half.gt.180) longecl_half=180 - t144=nt144(longecl_half) - tsky=(t144-2.7)*(144.0/nfreq)**2.6 + 2.7 !Tsky for obs freq - - xdop(NStation)=doppler - if(NStation.eq.2) then - HisGrid=MyGrid - go to 900 - endif - - doppler00=2.0*xdop(1) - doppler=xdop(1)+xdop(2) -! if(mode.eq.3) doppler=2.0*xdop(1) - dBMoon=-40.0*log10(dist/356903.) - sd=16.23*370152.0/dist - -! if(NStation.eq.1 .and. MoonDX.ne.0 .and. -! + (mode.eq.2 .or. mode.eq.5)) then - if(NStation.eq.1 .and. MoonDX.ne.0) then - poloffset=mod(poloffset2-poloffset1+720.0,180.0) - if(poloffset.gt.90.0) poloffset=poloffset-180.0 - x1=abs(cos(2*poloffset/rad)) - if(x1.lt.0.056234) x1=0.056234 - xnr=-20.0*log10(x1) - if(HisGrid(1:1).lt.'A' .or. HisGrid(1:1).gt.'R') xnr=0 - endif - - tr=80.0 !Good preamp - tskymin=13.0*(408.0/nfreq)**2.6 !Cold sky temperature - tsysmin=tskymin+tr - tsys=tsky+tr - dgrd=-10.0*log10(tsys/tsysmin) + dbMoon - 900 AzMoon0=Azmoon - ElMoon0=Elmoon - ntsky=nint(tsky) - -! auxHA = 15.0*(LST-auxra) !HA in degrees -! pi=3.14159265 -! pio2=0.5*pi -! call coord(pi,pio2-lat/rad,0.0,lat/rad,auxha*pi/180.0, -! + auxdec/rad,azaux,elaux) -! AzAux=azaux*rad -! ElAux=ElAux*rad - - return - - end diff --git a/libm65/astro.f90 b/libm65/astro.f90 new file mode 100644 index 000000000..d34819cc6 --- /dev/null +++ b/libm65/astro.f90 @@ -0,0 +1,105 @@ +subroutine astro(nyear,month,nday,uth,nfreq,Mygrid,NStation,MoonDX, & + AzSun,ElSun,AzMoon0,ElMoon0,ntsky,doppler00,doppler,dbMoon,RAMoon, & + DecMoon,HA,Dgrd,sd,poloffset,xnr,day,lon,lat,LST) + +! 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 LST + real lat,lon + integer*2 nt144(180) + +! common/echo/xdop(2),techo,AzMoon,ElMoon,mjd + real xdop(2) + + data rad/57.2957795/ + data nt144/ & + 234, 246, 257, 267, 275, 280, 283, 286, 291, 298, & + 305, 313, 322, 331, 341, 351, 361, 369, 376, 381, & + 383, 382, 379, 374, 370, 366, 363, 361, 363, 368, & + 376, 388, 401, 415, 428, 440, 453, 467, 487, 512, & + 544, 579, 607, 618, 609, 588, 563, 539, 512, 482, & + 450, 422, 398, 379, 363, 349, 334, 319, 302, 282, & + 262, 242, 226, 213, 205, 200, 198, 197, 196, 197, & + 200, 202, 204, 205, 204, 203, 202, 201, 203, 206, & + 212, 218, 223, 227, 231, 236, 240, 243, 247, 257, & + 276, 301, 324, 339, 346, 344, 339, 331, 323, 316, & + 312, 310, 312, 317, 327, 341, 358, 375, 392, 407, & + 422, 437, 451, 466, 480, 494, 511, 530, 552, 579, & + 612, 653, 702, 768, 863,1008,1232,1557,1966,2385, & + 2719,2924,3018,3038,2986,2836,2570,2213,1823,1461, & + 1163, 939, 783, 677, 602, 543, 494, 452, 419, 392, & + 373, 360, 353, 350, 350, 350, 350, 350, 350, 348, & + 344, 337, 329, 319, 307, 295, 284, 276, 272, 272, & + 273, 274, 274, 271, 266, 260, 252, 245, 238, 231/ + save + + call grid2deg(MyGrid,elon,lat) + lon=-elon + call sun(nyear,month,nday,uth,lon,lat,RASun,DecSun,LST,AzSun,ElSun,mjd,day) + + freq=nfreq*1.e6 + if(nfreq.eq.2) freq=1.8e6 + if(nfreq.eq.4) freq=3.5e6 + + call MoonDop(nyear,month,nday,uth,lon,lat,RAMoon,DecMoon,LST,HA, & + AzMoon,ElMoon,vr,dist) + +! Compute spatial polarization offset + xx=sin(lat/rad)*cos(ElMoon/rad) - cos(lat/rad)*cos(AzMoon/rad)*sin(ElMoon/rad) + yy=cos(lat/rad)*sin(AzMoon/rad) + if(NStation.eq.1) poloffset1=rad*atan2(yy,xx) + if(NStation.eq.2) poloffset2=rad*atan2(yy,xx) + + techo=2.0 * dist/2.99792458e5 !Echo delay time + doppler=-freq*vr/2.99792458e5 !One-way Doppler + + call coord(0.,0.,-1.570796,1.161639,RAMoon/rad,DecMoon/rad,el,eb) + longecl_half=nint(rad*el/2.0) + if(longecl_half.lt.1 .or. longecl_half.gt.180) longecl_half=180 + t144=nt144(longecl_half) + tsky=(t144-2.7)*(144.0/nfreq)**2.6 + 2.7 !Tsky for obs freq + + xdop(NStation)=doppler + if(NStation.eq.2) then + HisGrid=MyGrid + go to 900 + endif + + doppler00=2.0*xdop(1) + doppler=xdop(1)+xdop(2) +! if(mode.eq.3) doppler=2.0*xdop(1) + dBMoon=-40.0*log10(dist/356903.) + sd=16.23*370152.0/dist + +! if(NStation.eq.1 .and. MoonDX.ne.0 .and. +! + (mode.eq.2 .or. mode.eq.5)) then + if(NStation.eq.1 .and. MoonDX.ne.0) then + poloffset=mod(poloffset2-poloffset1+720.0,180.0) + if(poloffset.gt.90.0) poloffset=poloffset-180.0 + x1=abs(cos(2*poloffset/rad)) + if(x1.lt.0.056234) x1=0.056234 + xnr=-20.0*log10(x1) + if(HisGrid(1:1).lt.'A' .or. HisGrid(1:1).gt.'R') xnr=0 + endif + + tr=80.0 !Good preamp + tskymin=13.0*(408.0/nfreq)**2.6 !Cold sky temperature + tsysmin=tskymin+tr + tsys=tsky+tr + dgrd=-10.0*log10(tsys/tsysmin) + dbMoon +900 AzMoon0=Azmoon + ElMoon0=Elmoon + ntsky=nint(tsky) + +! auxHA = 15.0*(LST-auxra) !HA in degrees +! pi=3.14159265 +! pio2=0.5*pi +! call coord(pi,pio2-lat/rad,0.0,lat/rad,auxha*pi/180.0, +! + auxdec/rad,azaux,elaux) +! AzAux=azaux*rad +! ElAux=ElAux*rad + + return +end subroutine astro diff --git a/libm65/ccf2.f b/libm65/ccf2.f deleted file mode 100644 index c2dfadbce..000000000 --- a/libm65/ccf2.f +++ /dev/null @@ -1,45 +0,0 @@ - subroutine ccf2(ss,nz,nflip,ccfbest,lagpk) - -! parameter (LAGMAX=60) - parameter (LAGMAX=200) - real ss(nz) - real ccf(-LAGMAX:LAGMAX) - integer npr(126) - -C The JT65 pseudo-random sync pattern: - data npr/ - + 1,0,0,1,1,0,0,0,1,1,1,1,1,1,0,1,0,1,0,0, - + 0,1,0,1,1,0,0,1,0,0,0,1,1,1,0,0,1,1,1,1, - + 0,1,1,0,1,1,1,1,0,0,0,1,1,0,1,0,1,0,1,1, - + 0,0,1,1,0,1,0,1,0,1,0,0,1,0,0,0,0,0,0,1, - + 1,0,0,0,0,0,0,0,1,1,0,1,0,0,1,0,1,1,0,1, - + 0,1,0,1,0,0,1,1,0,0,1,0,0,1,0,0,0,0,1,1, - + 1,1,1,1,1,1/ - save - - ccfbest=0. - lag1=-LAGMAX - lag2=LAGMAX - do lag=lag1,lag2 - s0=0. - s1=0. - do i=1,126 - j=2*(8*i + 43) + lag - if(j.ge.1 .and. j.le.nz-8) then - x=ss(j)+ss(j+8) !Add two half-symbol contributions - if(npr(i).eq.0) then - s0=s0 + x - else - s1=s1 + x - endif - endif - enddo - ccf(lag)=nflip*(s1-s0) - if(ccf(lag).gt.ccfbest) then - ccfbest=ccf(lag) - lagpk=lag - endif - enddo - - return - end diff --git a/libm65/ccf2.f90 b/libm65/ccf2.f90 new file mode 100644 index 000000000..287e70ffb --- /dev/null +++ b/libm65/ccf2.f90 @@ -0,0 +1,45 @@ +subroutine ccf2(ss,nz,nflip,ccfbest,lagpk) + +! parameter (LAGMAX=60) + parameter (LAGMAX=200) + 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=-LAGMAX + lag2=LAGMAX + do lag=lag1,lag2 + s0=0. + s1=0. + do i=1,126 + j=2*(8*i + 43) + lag + if(j.ge.1 .and. j.le.nz-8) then + x=ss(j)+ss(j+8) !Add two half-symbol contributions + if(npr(i).eq.0) then + s0=s0 + x + else + s1=s1 + x + endif + endif + enddo + ccf(lag)=nflip*(s1-s0) + if(ccf(lag).gt.ccfbest) then + ccfbest=ccf(lag) + lagpk=lag + endif + enddo + + return +end subroutine ccf2 diff --git a/libm65/chkhist.f b/libm65/chkhist.f deleted file mode 100644 index 1b1a694b1..000000000 --- a/libm65/chkhist.f +++ /dev/null @@ -1,23 +0,0 @@ - subroutine chkhist(mrsym,nmax,ipk) - - integer mrsym(63) - integer hist(0:63) - - do i=0,63 - hist(i)=0 - enddo - do j=1,63 - i=mrsym(j) - hist(i)=hist(i)+1 - enddo - - nmax=0 - do i=0,63 - if(hist(i).gt.nmax) then - nmax=hist(i) - ipk=i+1 - endif - enddo - - return - end diff --git a/libm65/chkhist.f90 b/libm65/chkhist.f90 new file mode 100644 index 000000000..c814c3ca3 --- /dev/null +++ b/libm65/chkhist.f90 @@ -0,0 +1,23 @@ +subroutine chkhist(mrsym,nmax,ipk) + + integer mrsym(63) + integer hist(0:63) + + do i=0,63 + hist(i)=0 + enddo + do j=1,63 + i=mrsym(j) + hist(i)=hist(i)+1 + enddo + + nmax=0 + do i=0,63 + if(hist(i).gt.nmax) then + nmax=hist(i) + ipk=i+1 + endif + enddo + + return +end subroutine chkhist diff --git a/libm65/chkmsg.f b/libm65/chkmsg.f deleted file mode 100644 index 9dd7ef577..000000000 --- a/libm65/chkmsg.f +++ /dev/null @@ -1,32 +0,0 @@ - subroutine chkmsg(message,cok,nspecial,flip) - - character message*22,cok*3 - - nspecial=0 - flip=1.0 - cok=" " - - do i=22,1,-1 - if(message(i:i).ne.' ') go to 10 - enddo - i=22 - - 10 if(i.ge.11) then - if ((message(i-3:i).eq.' OOO') .or. - + (message(20:22).eq.' OO')) then - cok='OOO' - flip=-1.0 - if(message(20:22).eq.' OO') then - message=message(1:19) - else - message=message(1:i-4) - endif - endif - endif - - if(message(1:3).eq.'RO ') nspecial=2 - if(message(1:4).eq.'RRR ') nspecial=3 - if(message(1:3).eq.'73 ') nspecial=4 - - return - end diff --git a/libm65/chkmsg.f90 b/libm65/chkmsg.f90 new file mode 100644 index 000000000..457a80924 --- /dev/null +++ b/libm65/chkmsg.f90 @@ -0,0 +1,31 @@ +subroutine chkmsg(message,cok,nspecial,flip) + + character message*22,cok*3 + + nspecial=0 + flip=1.0 + cok=" " + + do i=22,1,-1 + if(message(i:i).ne.' ') go to 10 + enddo + i=22 + +10 if(i.ge.11) then + if ((message(i-3:i).eq.' OOO') .or. (message(20:22).eq.' OO')) then + cok='OOO' + flip=-1.0 + if(message(20:22).eq.' OO') then + message=message(1:19) + else + message=message(1:i-4) + endif + endif + endif + + if(message(1:3).eq.'RO ') nspecial=2 + if(message(1:4).eq.'RRR ') nspecial=3 + if(message(1:3).eq.'73 ') nspecial=4 + + return +end subroutine chkmsg diff --git a/libm65/coord.f b/libm65/coord.f deleted file mode 100644 index 3fb7c647d..000000000 --- a/libm65/coord.f +++ /dev/null @@ -1,41 +0,0 @@ - SUBROUTINE COORD(A0,B0,AP,BP,A1,B1,A2,B2) - -C Examples: -C 1. From ha,dec to az,el: -C call coord(pi,pio2-lat,0.,lat,ha,dec,az,el) -C 2. From az,el to ha,dec: -C call coord(pi,pio2-lat,0.,lat,az,el,ha,dec) -C 3. From ra,dec to l,b -C call coord(4.635594495,-0.504691042,3.355395488,0.478220215, -C ra,dec,l,b) -C 4. From l,b to ra,dec -C call coord(1.705981071d0,-1.050357016d0,2.146800277d0, -C 0.478220215d0,l,b,ra,dec) -C 5. From ra,dec to ecliptic latitude (eb) and longitude (el): -C call coord(0.d0,0.d0,-pio2,pio2-23.443*pi/180,ra,dec,el,eb) -C 6. From ecliptic latitude (eb) and longitude (el) to ra,dec: -C call coord(0.d0,0.d0,-pio2,pio2-23.443*pi/180,el,eb,ra,dec) - - - SB0=sin(B0) - CB0=cos(B0) - SBP=sin(BP) - CBP=cos(BP) - SB1=sin(B1) - CB1=cos(B1) - SB2=SBP*SB1 + CBP*CB1*cos(AP-A1) - CB2=SQRT(1.e0-SB2**2) - B2=atan(SB2/CB2) - SAA=sin(AP-A1)*CB1/CB2 - CAA=(SB1-SB2*SBP)/(CB2*CBP) - CBB=SB0/CBP - SBB=sin(AP-A0)*CB0 - SA2=SAA*CBB-CAA*SBB - CA2=CAA*CBB+SAA*SBB - TA2O2=0.0 !Shut up compiler warnings. -db - IF(CA2.LE.0.e0) TA2O2=(1.e0-CA2)/SA2 - IF(CA2.GT.0.e0) TA2O2=SA2/(1.e0+CA2) - A2=2.e0*atan(TA2O2) - IF(A2.LT.0.e0) A2=A2+6.2831853 - RETURN - END diff --git a/libm65/coord.f90 b/libm65/coord.f90 new file mode 100644 index 000000000..8934b612b --- /dev/null +++ b/libm65/coord.f90 @@ -0,0 +1,40 @@ +SUBROUTINE COORD(A0,B0,AP,BP,A1,B1,A2,B2) + +! Examples: +! 1. From ha,dec to az,el: +! call coord(pi,pio2-lat,0.,lat,ha,dec,az,el) +! 2. From az,el to ha,dec: +! call coord(pi,pio2-lat,0.,lat,az,el,ha,dec) +! 3. From ra,dec to l,b +! call coord(4.635594495,-0.504691042,3.355395488,0.478220215, +! ra,dec,l,b) +! 4. From l,b to ra,dec +! call coord(1.705981071d0,-1.050357016d0,2.146800277d0, +! 0.478220215d0,l,b,ra,dec) +! 5. From ra,dec to ecliptic latitude (eb) and longitude (el): +! call coord(0.d0,0.d0,-pio2,pio2-23.443*pi/180,ra,dec,el,eb) +! 6. From ecliptic latitude (eb) and longitude (el) to ra,dec: +! call coord(0.d0,0.d0,-pio2,pio2-23.443*pi/180,el,eb,ra,dec) + + SB0=sin(B0) + CB0=cos(B0) + SBP=sin(BP) + CBP=cos(BP) + SB1=sin(B1) + CB1=cos(B1) + SB2=SBP*SB1 + CBP*CB1*cos(AP-A1) + CB2=SQRT(1.e0-SB2**2) + B2=atan(SB2/CB2) + SAA=sin(AP-A1)*CB1/CB2 + CAA=(SB1-SB2*SBP)/(CB2*CBP) + CBB=SB0/CBP + SBB=sin(AP-A0)*CB0 + SA2=SAA*CBB-CAA*SBB + CA2=CAA*CBB+SAA*SBB + TA2O2=0.0 !Shut up compiler warnings. -db + IF(CA2.LE.0.e0) TA2O2=(1.e0-CA2)/SA2 + IF(CA2.GT.0.e0) TA2O2=SA2/(1.e0+CA2) + A2=2.e0*atan(TA2O2) + IF(A2.LT.0.e0) A2=A2+6.2831853 + RETURN +END SUBROUTINE COORD diff --git a/libm65/dcoord.f b/libm65/dcoord.f deleted file mode 100644 index b3d12d846..000000000 --- a/libm65/dcoord.f +++ /dev/null @@ -1,40 +0,0 @@ - SUBROUTINE DCOORD(A0,B0,AP,BP,A1,B1,A2,B2) - - implicit real*8 (a-h,o-z) -C Examples: -C 1. From ha,dec to az,el: -C call coord(pi,pio2-lat,0.,lat,ha,dec,az,el) -C 2. From az,el to ha,dec: -C call coord(pi,pio2-lat,0.,lat,az,el,ha,dec) -C 3. From ra,dec to l,b -C call coord(4.635594495,-0.504691042,3.355395488,0.478220215, -C ra,dec,l,b) -C 4. From l,b to ra,dec -C call coord(1.705981071d0,-1.050357016d0,2.146800277d0, -C 0.478220215d0,l,b,ra,dec) -C 5. From ecliptic latitude (eb) and longitude (el) to ra, dec: -C call coord(0.d0,0.d0,-pio2,pio2-23.443*pi/180,ra,dec,el,eb) - - SB0=sin(B0) - CB0=cos(B0) - SBP=sin(BP) - CBP=cos(BP) - SB1=sin(B1) - CB1=cos(B1) - SB2=SBP*SB1 + CBP*CB1*cos(AP-A1) - CB2=SQRT(1.D0-SB2**2) - B2=atan(SB2/CB2) - SAA=sin(AP-A1)*CB1/CB2 - CAA=(SB1-SB2*SBP)/(CB2*CBP) - CBB=SB0/CBP - SBB=sin(AP-A0)*CB0 - SA2=SAA*CBB-CAA*SBB - CA2=CAA*CBB+SAA*SBB - TA2O2=0.0 !Shut up compiler warnings. -db - IF(CA2.LE.0.D0) TA2O2=(1.D0-CA2)/SA2 - IF(CA2.GT.0.D0) TA2O2=SA2/(1.D0+CA2) - A2=2.D0*atan(TA2O2) - IF(A2.LT.0.D0) A2=A2+6.2831853071795864D0 - - RETURN - END diff --git a/libm65/dcoord.f90 b/libm65/dcoord.f90 new file mode 100644 index 000000000..5ef6877aa --- /dev/null +++ b/libm65/dcoord.f90 @@ -0,0 +1,40 @@ +SUBROUTINE DCOORD(A0,B0,AP,BP,A1,B1,A2,B2) + + implicit real*8 (a-h,o-z) +! Examples: +! 1. From ha,dec to az,el: +! call coord(pi,pio2-lat,0.,lat,ha,dec,az,el) +! 2. From az,el to ha,dec: +! call coord(pi,pio2-lat,0.,lat,az,el,ha,dec) +! 3. From ra,dec to l,b +! call coord(4.635594495,-0.504691042,3.355395488,0.478220215, +! ra,dec,l,b) +! 4. From l,b to ra,dec +! call coord(1.705981071d0,-1.050357016d0,2.146800277d0, +! 0.478220215d0,l,b,ra,dec) +! 5. From ecliptic latitude (eb) and longitude (el) to ra, dec: +! call coord(0.d0,0.d0,-pio2,pio2-23.443*pi/180,ra,dec,el,eb) + + SB0=sin(B0) + CB0=cos(B0) + SBP=sin(BP) + CBP=cos(BP) + SB1=sin(B1) + CB1=cos(B1) + SB2=SBP*SB1 + CBP*CB1*cos(AP-A1) + CB2=SQRT(1.D0-SB2**2) + B2=atan(SB2/CB2) + SAA=sin(AP-A1)*CB1/CB2 + CAA=(SB1-SB2*SBP)/(CB2*CBP) + CBB=SB0/CBP + SBB=sin(AP-A0)*CB0 + SA2=SAA*CBB-CAA*SBB + CA2=CAA*CBB+SAA*SBB + TA2O2=0.0 !Shut up compiler warnings. -db + IF(CA2.LE.0.D0) TA2O2=(1.D0-CA2)/SA2 + IF(CA2.GT.0.D0) TA2O2=SA2/(1.D0+CA2) + A2=2.D0*atan(TA2O2) + IF(A2.LT.0.D0) A2=A2+6.2831853071795864D0 + + RETURN +END SUBROUTINE DCOORD diff --git a/libm65/decode65b.f b/libm65/decode65b.f deleted file mode 100644 index b94c05534..000000000 --- a/libm65/decode65b.f +++ /dev/null @@ -1,49 +0,0 @@ - subroutine decode65b(s2,flip,mycall,hiscall,hisgrid,mode65, - + neme,ndepth,nqd,nkv,nhist,qual,decoded,s3,sy) - - real s2(66,126) - real s3(64,63),sy(63) - logical first,ltext - character decoded*22,deepmsg*22 - character mycall*12,hiscall*12,hisgrid*6 - common/prcom/pr(126),mdat(126),mref(126,2),mdat2(126),mref2(126,2) - data first/.true./ - save - - if(first) call setup65 - first=.false. - - do j=1,63 - k=mdat(j) !Points to data symbol - if(flip.lt.0.0) k=mdat2(j) - do i=1,64 - s3(i,j)=s2(i+2,k) - enddo - k=mdat2(j) !Points to data symbol - if(flip.lt.0.0) k=mdat(j) - sy(j)=s2(1,k) - enddo - - nadd=mode65 - call extract(s3,nadd,ncount,nhist,decoded,ltext) !Extract the message -C Suppress "birdie messages" and other garbage decodes: - if(decoded(1:7).eq.'000AAA ') ncount=-1 - if(decoded(1:7).eq.'0L6MWK ') ncount=-1 - if(flip.lt.0.0 .and. ltext) ncount=-1 - nkv=1 - if(ncount.lt.0) then - nkv=0 - decoded=' ' - endif - - qual=0. - if(ndepth.ge.1 .and. (nqd.eq.1 .or. flip.eq.1.0)) then - call deep65(s3,mode65,neme,flip,mycall,hiscall, - + hisgrid,deepmsg,qual) - if(nqd.ne.1 .and. qual.lt.10.0) qual=0.0 - if(ndepth.lt.2 .and. qual.lt.6.0) qual=0.0 - endif - if(nkv.eq.0 .and. qual.ge.1.0) decoded=deepmsg - - return - end diff --git a/libm65/decode65b.f90 b/libm65/decode65b.f90 new file mode 100644 index 000000000..9b3ebccda --- /dev/null +++ b/libm65/decode65b.f90 @@ -0,0 +1,48 @@ +subroutine decode65b(s2,flip,mycall,hiscall,hisgrid,mode65,neme,ndepth, & + nqd,nkv,nhist,qual,decoded,s3,sy) + + real s2(66,126) + real s3(64,63),sy(63) + logical first,ltext + character decoded*22,deepmsg*22 + character mycall*12,hiscall*12,hisgrid*6 + common/prcom/pr(126),mdat(126),mref(126,2),mdat2(126),mref2(126,2) + data first/.true./ + save + + if(first) call setup65 + first=.false. + + do j=1,63 + k=mdat(j) !Points to data symbol + if(flip.lt.0.0) k=mdat2(j) + do i=1,64 + s3(i,j)=s2(i+2,k) + enddo + k=mdat2(j) !Points to data symbol + if(flip.lt.0.0) k=mdat(j) + sy(j)=s2(1,k) + enddo + + nadd=mode65 + call extract(s3,nadd,ncount,nhist,decoded,ltext) !Extract the message +! Suppress "birdie messages" and other garbage decodes: + if(decoded(1:7).eq.'000AAA ') ncount=-1 + if(decoded(1:7).eq.'0L6MWK ') ncount=-1 + if(flip.lt.0.0 .and. ltext) ncount=-1 + nkv=1 + if(ncount.lt.0) then + nkv=0 + decoded=' ' + endif + + qual=0. + if(ndepth.ge.1 .and. (nqd.eq.1 .or. flip.eq.1.0)) then + call deep65(s3,mode65,neme,flip,mycall,hiscall,hisgrid,deepmsg,qual) + if(nqd.ne.1 .and. qual.lt.10.0) qual=0.0 + if(ndepth.lt.2 .and. qual.lt.6.0) qual=0.0 + endif + if(nkv.eq.0 .and. qual.ge.1.0) decoded=deepmsg + + return +end subroutine decode65b diff --git a/libm65/deg2grid.f b/libm65/deg2grid.f deleted file mode 100644 index 8c64028a8..000000000 --- a/libm65/deg2grid.f +++ /dev/null @@ -1,30 +0,0 @@ - subroutine deg2grid(dlong0,dlat,grid) - - real dlong !West longitude (deg) - real dlat !Latitude (deg) - character grid*6 - - dlong=dlong0 - if(dlong.lt.-180.0) dlong=dlong+360.0 - if(dlong.gt.180.0) dlong=dlong-360.0 - -C Convert to units of 5 min of longitude, working east from 180 deg. - nlong=60.0*(180.0-dlong)/5.0 - n1=nlong/240 !20-degree field - n2=(nlong-240*n1)/24 !2 degree square - n3=nlong-240*n1-24*n2 !5 minute subsquare - grid(1:1)=char(ichar('A')+n1) - grid(3:3)=char(ichar('0')+n2) - grid(5:5)=char(ichar('a')+n3) - -C Convert to units of 2.5 min of latitude, working north from -90 deg. - nlat=60.0*(dlat+90)/2.5 - n1=nlat/240 !10-degree field - n2=(nlat-240*n1)/24 !1 degree square - n3=nlat-240*n1-24*n2 !2.5 minuts subsquare - grid(2:2)=char(ichar('A')+n1) - grid(4:4)=char(ichar('0')+n2) - grid(6:6)=char(ichar('a')+n3) - - return - end diff --git a/libm65/deg2grid.f90 b/libm65/deg2grid.f90 new file mode 100644 index 000000000..9ca3602f8 --- /dev/null +++ b/libm65/deg2grid.f90 @@ -0,0 +1,30 @@ +subroutine deg2grid(dlong0,dlat,grid) + + real dlong !West longitude (deg) + real dlat !Latitude (deg) + character grid*6 + + dlong=dlong0 + if(dlong.lt.-180.0) dlong=dlong+360.0 + if(dlong.gt.180.0) dlong=dlong-360.0 + +! Convert to units of 5 min of longitude, working east from 180 deg. + nlong=60.0*(180.0-dlong)/5.0 + n1=nlong/240 !20-degree field + n2=(nlong-240*n1)/24 !2 degree square + n3=nlong-240*n1-24*n2 !5 minute subsquare + grid(1:1)=char(ichar('A')+n1) + grid(3:3)=char(ichar('0')+n2) + grid(5:5)=char(ichar('a')+n3) + +! Convert to units of 2.5 min of latitude, working north from -90 deg. + nlat=60.0*(dlat+90)/2.5 + n1=nlat/240 !10-degree field + n2=(nlat-240*n1)/24 !1 degree square + n3=nlat-240*n1-24*n2 !2.5 minuts subsquare + grid(2:2)=char(ichar('A')+n1) + grid(4:4)=char(ichar('0')+n2) + grid(6:6)=char(ichar('a')+n3) + + return +end subroutine deg2grid diff --git a/libm65/dot.f b/libm65/dot.f deleted file mode 100644 index 8e2d826bf..000000000 --- a/libm65/dot.f +++ /dev/null @@ -1,11 +0,0 @@ - real*8 function dot(x,y) - - real*8 x(3),y(3) - - dot=0.d0 - do i=1,3 - dot=dot+x(i)*y(i) - enddo - - return - end diff --git a/libm65/dot.f90 b/libm65/dot.f90 new file mode 100644 index 000000000..5829e8787 --- /dev/null +++ b/libm65/dot.f90 @@ -0,0 +1,11 @@ +real*8 function dot(x,y) + + real*8 x(3),y(3) + + dot=0.d0 + do i=1,3 + dot=dot+x(i)*y(i) + enddo + + return +end function dot diff --git a/libm65/encode65.f b/libm65/encode65.f deleted file mode 100644 index 670c2e583..000000000 --- a/libm65/encode65.f +++ /dev/null @@ -1,13 +0,0 @@ - subroutine encode65(message,sent) - - character message*22 - integer dgen(12) - integer sent(63) - - call packmsg(message,dgen) - call rs_encode(dgen,sent) - call interleave63(sent,1) - call graycode(sent,63,1) - - return - end diff --git a/libm65/encode65.f90 b/libm65/encode65.f90 new file mode 100644 index 000000000..e7567bd3e --- /dev/null +++ b/libm65/encode65.f90 @@ -0,0 +1,13 @@ +subroutine encode65(message,sent) + + character message*22 + integer dgen(12) + integer sent(63) + + call packmsg(message,dgen) + call rs_encode(dgen,sent) + call interleave63(sent,1) + call graycode(sent,63,1) + + return +end subroutine encode65 diff --git a/libm65/fchisq.f b/libm65/fchisq.f deleted file mode 100644 index fa311d95f..000000000 --- a/libm65/fchisq.f +++ /dev/null @@ -1,77 +0,0 @@ - real function fchisq(cx,cy,npts,nfast,fsample,nflip,a, - + ccfmax,dtmax) - - parameter (NMAX=60*96000) !Samples per 60 s - complex cx(npts),cy(npts) - real a(5) - complex w,wstep,za,zb,z - real ss(3000) - complex csx(0:NMAX/64),csy(0:NMAX/64) - data twopi/6.283185307/a1,a2,a3/99.,99.,99./ - save - - call timer('fchisq ',0) - baud=nfast*11025.0/4096.0 - nsps=nint(fsample/baud) !Samples per symbol - nsph=nsps/2 !Samples per half-symbol - ndiv=16 !Output ss() steps per symbol - nout=ndiv*npts/nsps - dtstep=1.0/(ndiv*baud) !Time per output step - - if(a(1).ne.a1 .or. a(2).ne.a2 .or. a(3).ne.a3) then - a1=a(1) - a2=a(2) - a3=a(3) - -C Mix and integrate the complex X and Y signals - csx(0)=0. - csy(0)=0. - w=1.0 - x0=0.5*(npts+1) - s=2.0/npts - do i=1,npts - x=s*(i-x0) - if(mod(i,100).eq.1) then - p2=1.5*x*x - 0.5 -! p3=2.5*(x**3) - 1.5*x -! p4=4.375*(x**4) - 3.75*(x**2) + 0.375 - dphi=(a(1) + x*a(2) + p2*a(3)) * (twopi/fsample) - wstep=cmplx(cos(dphi),sin(dphi)) - endif - w=w*wstep - csx(i)=csx(i-1) + w*cx(i) - csy(i)=csy(i-1) + w*cy(i) - enddo - endif - -C Compute 1/2-symbol powers at 1/16-symbol steps. - fac=1.e-4 - pol=a(4)/57.2957795 - aa=cos(pol) - bb=sin(pol) - - do i=1,nout - j=i*nsps/ndiv - k=j-nsph - ss(i)=0. - if(k.ge.1) then - za=csx(j)-csx(k) - zb=csy(j)-csy(k) - z=aa*za + bb*zb - ss(i)=fac*(real(z)**2 + aimag(z)**2) - endif - enddo - - ccfmax=0. - call timer('ccf2 ',0) - call ccf2(ss,nout,nflip,ccf,lagpk) - call timer('ccf2 ',1) - if(ccf.gt.ccfmax) then - ccfmax=ccf - dtmax=lagpk*dtstep - endif - fchisq=-ccfmax - - call timer('fchisq ',1) - return - end diff --git a/libm65/fchisq.f90 b/libm65/fchisq.f90 new file mode 100644 index 000000000..ced1a92e7 --- /dev/null +++ b/libm65/fchisq.f90 @@ -0,0 +1,76 @@ +real function fchisq(cx,cy,npts,nfast,fsample,nflip,a,ccfmax,dtmax) + + parameter (NMAX=60*96000) !Samples per 60 s + complex cx(npts),cy(npts) + real a(5) + complex w,wstep,za,zb,z + real ss(3000) + complex csx(0:NMAX/64),csy(0:NMAX/64) + data twopi/6.283185307/a1,a2,a3/99.,99.,99./ + save + + call timer('fchisq ',0) + baud=nfast*11025.0/4096.0 + nsps=nint(fsample/baud) !Samples per symbol + nsph=nsps/2 !Samples per half-symbol + ndiv=16 !Output ss() steps per symbol + nout=ndiv*npts/nsps + dtstep=1.0/(ndiv*baud) !Time per output step + + if(a(1).ne.a1 .or. a(2).ne.a2 .or. a(3).ne.a3) then + a1=a(1) + a2=a(2) + a3=a(3) + +! Mix and integrate the complex X and Y signals + csx(0)=0. + csy(0)=0. + w=1.0 + x0=0.5*(npts+1) + s=2.0/npts + do i=1,npts + x=s*(i-x0) + if(mod(i,100).eq.1) then + p2=1.5*x*x - 0.5 +! p3=2.5*(x**3) - 1.5*x +! p4=4.375*(x**4) - 3.75*(x**2) + 0.375 + dphi=(a(1) + x*a(2) + p2*a(3)) * (twopi/fsample) + wstep=cmplx(cos(dphi),sin(dphi)) + endif + w=w*wstep + csx(i)=csx(i-1) + w*cx(i) + csy(i)=csy(i-1) + w*cy(i) + enddo + endif + +! Compute 1/2-symbol powers at 1/16-symbol steps. + fac=1.e-4 + pol=a(4)/57.2957795 + aa=cos(pol) + bb=sin(pol) + + do i=1,nout + j=i*nsps/ndiv + k=j-nsph + ss(i)=0. + if(k.ge.1) then + za=csx(j)-csx(k) + zb=csy(j)-csy(k) + z=aa*za + bb*zb + ss(i)=fac*(real(z)**2 + aimag(z)**2) + endif + enddo + + ccfmax=0. + call timer('ccf2 ',0) + call ccf2(ss,nout,nflip,ccf,lagpk) + call timer('ccf2 ',1) + if(ccf.gt.ccfmax) then + ccfmax=ccf + dtmax=lagpk*dtstep + endif + fchisq=-ccfmax + call timer('fchisq ',1) + + return +end function fchisq diff --git a/libm65/fil6521.f b/libm65/fil6521.f deleted file mode 100644 index 32feb8f44..000000000 --- a/libm65/fil6521.f +++ /dev/null @@ -1,44 +0,0 @@ - subroutine fil6521(c1,n1,c2,n2) - -C FIR lowpass filter designed using ScopeFIR - -C Pass #1 Pass #2 -C----------------------------------------------- -C fsample (Hz) 1378.125 Input sample rate -C Ntaps 21 Number of filter taps -C fc (Hz) 40 Cutoff frequency -C fstop (Hz) 172.266 Lower limit of stopband -C Ripple (dB) 0.1 Ripple in passband -C Stop Atten (dB) 38 Stopband attenuation -C fout (Hz) 344.531 Output sample rate - - parameter (NTAPS=21) - parameter (NH=NTAPS/2) - parameter (NDOWN=4) !Downsample ratio = 1/4 - complex c1(n1) - complex c2(n1/NDOWN) - -C Filter coefficients: - real a(-NH:NH) - data a/ - + -0.011958606980,-0.013888627387,-0.015601306443,-0.010602249570, - + 0.003804023436, 0.028320058273, 0.060903935217, 0.096841904411, - + 0.129639871228, 0.152644580853, 0.160917511283, 0.152644580853, - + 0.129639871228, 0.096841904411, 0.060903935217, 0.028320058273, - + 0.003804023436,-0.010602249570,-0.015601306443,-0.013888627387, - + -0.011958606980/ - - n2=(n1-NTAPS+NDOWN)/NDOWN - k0=NH-NDOWN+1 - -C Loop over all output samples - do i=1,n2 - c2(i)=0. - k=k0 + NDOWN*i - do j=-NH,NH - c2(i)=c2(i) + c1(j+k)*a(j) - enddo - enddo - - return - end diff --git a/libm65/fil6521.f90 b/libm65/fil6521.f90 new file mode 100644 index 000000000..b4c808c10 --- /dev/null +++ b/libm65/fil6521.f90 @@ -0,0 +1,44 @@ +subroutine fil6521(c1,n1,c2,n2) + +! FIR lowpass filter designed using ScopeFIR + +! Pass #1 Pass #2 +!----------------------------------------------- +! fsample (Hz) 1378.125 Input sample rate +! Ntaps 21 Number of filter taps +! fc (Hz) 40 Cutoff frequency +! fstop (Hz) 172.266 Lower limit of stopband +! Ripple (dB) 0.1 Ripple in passband +! Stop Atten (dB) 38 Stopband attenuation +! fout (Hz) 344.531 Output sample rate + + parameter (NTAPS=21) + parameter (NH=NTAPS/2) + parameter (NDOWN=4) !Downsample ratio = 1/4 + complex c1(n1) + complex c2(n1/NDOWN) + +! Filter coefficients: + real a(-NH:NH) + data a/ & + -0.011958606980,-0.013888627387,-0.015601306443,-0.010602249570, & + 0.003804023436, 0.028320058273, 0.060903935217, 0.096841904411, & + 0.129639871228, 0.152644580853, 0.160917511283, 0.152644580853, & + 0.129639871228, 0.096841904411, 0.060903935217, 0.028320058273, & + 0.003804023436,-0.010602249570,-0.015601306443,-0.013888627387, & + -0.011958606980/ + + n2=(n1-NTAPS+NDOWN)/NDOWN + k0=NH-NDOWN+1 + +! Loop over all output samples + do i=1,n2 + c2(i)=0. + k=k0 + NDOWN*i + do j=-NH,NH + c2(i)=c2(i) + c1(j+k)*a(j) + enddo + enddo + + return +end subroutine fil6521 diff --git a/libm65/filbig.f b/libm65/filbig.f deleted file mode 100644 index e4c062130..000000000 --- a/libm65/filbig.f +++ /dev/null @@ -1,157 +0,0 @@ - subroutine filbig(dd,nmax,nfast,f0,newdat,nfsample,xpol, - + c4a,c4b,n4) - -C Filter and downsample complex data stored in array dd(4,nmax). -C Output is downsampled from 96000 Hz to 1375.125 Hz. - - parameter (MAXFFT1=5376000,MAXFFT2=77175) - real*4 dd(4,nmax) !Input data - complex ca(MAXFFT1),cb(MAXFFT1) !FFTs of input - complex c4a(MAXFFT2),c4b(MAXFFT2) !Output data - real*8 df - real halfpulse(8) !Impulse response of filter (one sided) - complex cfilt(MAXFFT2) !Filter (complex; imag = 0) - real rfilt(MAXFFT2) !Filter (real) - integer*8 plan1,plan2,plan3,plan4,plan5 - logical first,xpol - include 'fftw3.f' - equivalence (rfilt,cfilt) - data first/.true./,npatience/1/,nfast0/0/ - data halfpulse/114.97547150,36.57879257,-20.93789101, - + 5.89886379,1.59355187,-2.49138308,0.60910773,-0.04248129/ - save - - if(nmax.lt.0) go to 900 - - if(nfast.eq.1) then - nfft1=MAXFFT1 - nfft2=MAXFFT2 - if(nfsample.eq.95238) then - nfft1=5120000 - nfft2=74088 - endif - else - nfft1=2621440 - nfft2=37632 - if(nfsample.eq.95238) then - nfft1=2560000 - nfft2=37044 - endif - endif - - if(nfast.ne.nfast0) then - if(nfast0.ne.0) then - call sfftw_destroy_plan(plan1) - call sfftw_destroy_plan(plan2) - call sfftw_destroy_plan(plan3) - call sfftw_destroy_plan(plan4) - call sfftw_destroy_plan(plan5) - endif - endif - - if(first .or. nfast.ne.nfast0) then - nflags=FFTW_ESTIMATE - if(npatience.eq.1) nflags=FFTW_ESTIMATE_PATIENT - if(npatience.eq.2) nflags=FFTW_MEASURE - if(npatience.eq.3) nflags=FFTW_PATIENT - if(npatience.eq.4) nflags=FFTW_EXHAUSTIVE -C Plan the FFTs just once - call timer('FFTplans ',0) - call sfftw_plan_dft_1d(plan1,nfft1,ca,ca, - + FFTW_BACKWARD,nflags) - call sfftw_plan_dft_1d(plan2,nfft1,cb,cb, - + FFTW_BACKWARD,nflags) - call sfftw_plan_dft_1d(plan3,nfft2,c4a,c4a, - + FFTW_FORWARD,nflags) - call sfftw_plan_dft_1d(plan4,nfft2,c4b,c4b, - + FFTW_FORWARD,nflags) - call sfftw_plan_dft_1d(plan5,nfft2,cfilt,cfilt, - + FFTW_BACKWARD,nflags) - call timer('FFTplans ',1) - -C Convert impulse response to filter function - do i=1,nfft2 - cfilt(i)=0. - enddo - fac=0.00625/nfft1 - cfilt(1)=fac*halfpulse(1) - do i=2,8 - cfilt(i)=fac*halfpulse(i) - cfilt(nfft2+2-i)=fac*halfpulse(i) - enddo - call timer('FFTfilt ',0) - call sfftw_execute(plan5) - call timer('FFTfilt ',1) - - base=cfilt(nfft2/2+1) - do i=1,nfft2 - rfilt(i)=real(cfilt(i))-base - enddo - - df=96000.d0/nfft1 - if(nfsample.eq.95238) df=95238.1d0/nfft1 - first=.false. - endif - nfast0=nfast - -C When new data comes along, we need to compute a new "big FFT" -C If we just have a new f0, continue with the existing ca and cb. - - if(newdat.ne.0) then - nz=min(nmax,nfft1) - do i=1,nz - ca(i)=cmplx(dd(1,i),dd(2,i)) - if(xpol) cb(i)=cmplx(dd(3,i),dd(4,i)) - enddo - - if(nmax.lt.nfft1) then - do i=nmax+1,nfft1 - ca(i)=0. - if(xpol) cb(i)=0. - enddo - endif - call timer('FFTbig ',0) - call sfftw_execute(plan1) - if(xpol) call sfftw_execute(plan2) - call timer('FFTbig ',1) - newdat=0 - endif - -C NB: f0 is the frequency at which we want our filter centered. -C i0 is the bin number in ca and cb closest to f0. - - i0=nint(f0/df) + 1 - nh=nfft2/2 - do i=1,nh !Copy data into c4a and c4b, - j=i0+i-1 !and apply the filter function - if(j.ge.1 .and. j.le.nfft1) then - c4a(i)=rfilt(i)*ca(j) - if(xpol) c4b(i)=rfilt(i)*cb(j) - else - c4a(i)=0. - if(xpol) c4b(i)=0. - endif - enddo - do i=nh+1,nfft2 - j=i0+i-1-nfft2 - if(j.lt.1) j=j+nfft1 !nfft1 was nfft2 - c4a(i)=rfilt(i)*ca(j) - if(xpol) c4b(i)=rfilt(i)*cb(j) - enddo - -C Do the short reverse transform, to go back to time domain. - call timer('FFTsmall',0) - call sfftw_execute(plan3) - if(xpol) call sfftw_execute(plan4) - call timer('FFTsmall',1) - n4=min(nmax/64,nfft2) - go to 999 - - 900 call sfftw_destroy_plan(plan1) - call sfftw_destroy_plan(plan2) - call sfftw_destroy_plan(plan3) - call sfftw_destroy_plan(plan4) - call sfftw_destroy_plan(plan5) - - 999 return - end diff --git a/libm65/filbig.f90 b/libm65/filbig.f90 new file mode 100644 index 000000000..7ca1a92ed --- /dev/null +++ b/libm65/filbig.f90 @@ -0,0 +1,152 @@ +subroutine filbig(dd,nmax,nfast,f0,newdat,nfsample,xpol,c4a,c4b,n4) + +! Filter and downsample complex data stored in array dd(4,nmax). +! Output is downsampled from 96000 Hz to 1375.125 Hz. + + parameter (MAXFFT1=5376000,MAXFFT2=77175) + real*4 dd(4,nmax) !Input data + complex ca(MAXFFT1),cb(MAXFFT1) !FFTs of input + complex c4a(MAXFFT2),c4b(MAXFFT2) !Output data + real*8 df + real halfpulse(8) !Impulse response of filter (one sided) + complex cfilt(MAXFFT2) !Filter (complex; imag = 0) + real rfilt(MAXFFT2) !Filter (real) + integer*8 plan1,plan2,plan3,plan4,plan5 + logical first,xpol + include 'fftw3.f' + equivalence (rfilt,cfilt) + data first/.true./,npatience/1/,nfast0/0/ + data halfpulse/114.97547150,36.57879257,-20.93789101, & + 5.89886379,1.59355187,-2.49138308,0.60910773,-0.04248129/ + save + + if(nmax.lt.0) go to 900 + + if(nfast.eq.1) then + nfft1=MAXFFT1 + nfft2=MAXFFT2 + if(nfsample.eq.95238) then + nfft1=5120000 + nfft2=74088 + endif + else + nfft1=2621440 + nfft2=37632 + if(nfsample.eq.95238) then + nfft1=2560000 + nfft2=37044 + endif + endif + + if(nfast.ne.nfast0) then + if(nfast0.ne.0) then + call sfftw_destroy_plan(plan1) + call sfftw_destroy_plan(plan2) + call sfftw_destroy_plan(plan3) + call sfftw_destroy_plan(plan4) + call sfftw_destroy_plan(plan5) + endif + endif + + if(first .or. nfast.ne.nfast0) then + nflags=FFTW_ESTIMATE + if(npatience.eq.1) nflags=FFTW_ESTIMATE_PATIENT + if(npatience.eq.2) nflags=FFTW_MEASURE + if(npatience.eq.3) nflags=FFTW_PATIENT + if(npatience.eq.4) nflags=FFTW_EXHAUSTIVE + +! Plan the FFTs just once + call timer('FFTplans ',0) + call sfftw_plan_dft_1d(plan1,nfft1,ca,ca,FFTW_BACKWARD,nflags) + call sfftw_plan_dft_1d(plan2,nfft1,cb,cb,FFTW_BACKWARD,nflags) + call sfftw_plan_dft_1d(plan3,nfft2,c4a,c4a,FFTW_FORWARD,nflags) + call sfftw_plan_dft_1d(plan4,nfft2,c4b,c4b,FFTW_FORWARD,nflags) + call sfftw_plan_dft_1d(plan5,nfft2,cfilt,cfilt,FFTW_BACKWARD,nflags) + call timer('FFTplans ',1) + +! Convert impulse response to filter function + do i=1,nfft2 + cfilt(i)=0. + enddo + fac=0.00625/nfft1 + cfilt(1)=fac*halfpulse(1) + do i=2,8 + cfilt(i)=fac*halfpulse(i) + cfilt(nfft2+2-i)=fac*halfpulse(i) + enddo + call timer('FFTfilt ',0) + call sfftw_execute(plan5) + call timer('FFTfilt ',1) + + base=cfilt(nfft2/2+1) + do i=1,nfft2 + rfilt(i)=real(cfilt(i))-base + enddo + + df=96000.d0/nfft1 + if(nfsample.eq.95238) df=95238.1d0/nfft1 + first=.false. + endif + nfast0=nfast + +! When new data comes along, we need to compute a new "big FFT" +! If we just have a new f0, continue with the existing ca and cb. + + if(newdat.ne.0) then + nz=min(nmax,nfft1) + do i=1,nz + ca(i)=cmplx(dd(1,i),dd(2,i)) + if(xpol) cb(i)=cmplx(dd(3,i),dd(4,i)) + enddo + + if(nmax.lt.nfft1) then + do i=nmax+1,nfft1 + ca(i)=0. + if(xpol) cb(i)=0. + enddo + endif + call timer('FFTbig ',0) + call sfftw_execute(plan1) + if(xpol) call sfftw_execute(plan2) + call timer('FFTbig ',1) + newdat=0 + endif + +! NB: f0 is the frequency at which we want our filter centered. +! i0 is the bin number in ca and cb closest to f0. + + i0=nint(f0/df) + 1 + nh=nfft2/2 + do i=1,nh !Copy data into c4a and c4b, + j=i0+i-1 !and apply the filter function + if(j.ge.1 .and. j.le.nfft1) then + c4a(i)=rfilt(i)*ca(j) + if(xpol) c4b(i)=rfilt(i)*cb(j) + else + c4a(i)=0. + if(xpol) c4b(i)=0. + endif + enddo + do i=nh+1,nfft2 + j=i0+i-1-nfft2 + if(j.lt.1) j=j+nfft1 !nfft1 was nfft2 + c4a(i)=rfilt(i)*ca(j) + if(xpol) c4b(i)=rfilt(i)*cb(j) + enddo + +! Do the short reverse transform, to go back to time domain. + call timer('FFTsmall',0) + call sfftw_execute(plan3) + if(xpol) call sfftw_execute(plan4) + call timer('FFTsmall',1) + n4=min(nmax/64,nfft2) + go to 999 + +900 call sfftw_destroy_plan(plan1) + call sfftw_destroy_plan(plan2) + call sfftw_destroy_plan(plan3) + call sfftw_destroy_plan(plan4) + call sfftw_destroy_plan(plan5) + +999 return +end subroutine filbig diff --git a/libm65/geocentric.f b/libm65/geocentric.f deleted file mode 100644 index 1767f9cf8..000000000 --- a/libm65/geocentric.f +++ /dev/null @@ -1,17 +0,0 @@ - subroutine geocentric(alat,elev,hlt,erad) - - implicit real*8 (a-h,o-z) - -C IAU 1976 flattening f, equatorial radius a - f = 1.d0/298.257d0 - a = 6378140.d0 - c = 1.d0/sqrt(1.d0 + (-2.d0 + f)*f*sin(alat)*sin(alat)) - arcf = (a*c + elev)*cos(alat) - arsf = (a*(1.d0 - f)*(1.d0 - f)*c + elev)*sin(alat) - hlt = datan2(arsf,arcf) - erad = sqrt(arcf*arcf + arsf*arsf) - erad = 0.001d0*erad - - return - end - diff --git a/libm65/geocentric.f90 b/libm65/geocentric.f90 new file mode 100644 index 000000000..59ed90b51 --- /dev/null +++ b/libm65/geocentric.f90 @@ -0,0 +1,17 @@ +subroutine geocentric(alat,elev,hlt,erad) + + implicit real*8 (a-h,o-z) + +! IAU 1976 flattening f, equatorial radius a + f = 1.d0/298.257d0 + a = 6378140.d0 + c = 1.d0/sqrt(1.d0 + (-2.d0 + f)*f*sin(alat)*sin(alat)) + arcf = (a*c + elev)*cos(alat) + arsf = (a*(1.d0 - f)*(1.d0 - f)*c + elev)*sin(alat) + hlt = datan2(arsf,arcf) + erad = sqrt(arcf*arcf + arsf*arsf) + erad = 0.001d0*erad + + return +end subroutine geocentric + diff --git a/libm65/getpfx1.f b/libm65/getpfx1.f deleted file mode 100644 index 21489f522..000000000 --- a/libm65/getpfx1.f +++ /dev/null @@ -1,96 +0,0 @@ - subroutine getpfx1(callsign,k,nv2) - - character*12 callsign0,callsign,lof,rof - character*8 c - character addpfx*8,tpfx*4,tsfx*3 - logical ispfx,issfx,invalid - common/pfxcom/addpfx - include 'pfx.f' - - callsign0=callsign - nv2=0 - iz=index(callsign,' ') - 1 - if(iz.lt.0) iz=12 - islash=index(callsign(1:iz),'/') - k=0 - c=' ' - if(islash.gt.0 .and. islash.le.(iz-4)) then -! Add-on prefix - c=callsign(1:islash-1) - callsign=callsign(islash+1:iz) - do i=1,NZ - if(pfx(i)(1:4).eq.c) then - k=i - go to 10 - endif - enddo - if(addpfx.eq.c) then - k=449 - go to 10 - endif - - else if(islash.eq.(iz-1)) then -! Add-on suffix - c=callsign(islash+1:iz) - callsign=callsign(1:islash-1) - do i=1,NZ2 - if(sfx(i).eq.c(1:1)) then - k=400+i - go to 10 - endif - enddo - endif - - 10 if(islash.ne.0 .and.k.eq.0) then -! Original JT65 would force this compound callsign to be treated as -! plain text. In JT65v2, we will encode the prefix or suffix into nc1. -! The task here is to compute the proper value of k. - lof=callsign0(:islash-1) - rof=callsign0(islash+1:) - llof=len_trim(lof) - lrof=len_trim(rof) - ispfx=(llof.gt.0 .and. llof.le.4) - issfx=(lrof.gt.0 .and. lrof.le.3) - invalid=.not.(ispfx.or.issfx) - if(ispfx.and.issfx) then - if(llof.lt.3) issfx=.false. - if(lrof.lt.3) ispfx=.false. - if(ispfx.and.issfx) then - i=ichar(callsign0(islash-1:islash-1)) - if(i.ge.ichar('0') .and. i.le.ichar('9')) then - issfx=.false. - else - ispfx=.false. - endif - endif - endif - - if(invalid) then - k=-1 - else - if(ispfx) then - tpfx=lof - k=nchar(tpfx(1:1)) - k=37*k + nchar(tpfx(2:2)) - k=37*k + nchar(tpfx(3:3)) - k=37*k + nchar(tpfx(4:4)) - nv2=1 - i=index(callsign0,'/') - callsign=callsign0(:i-1) - callsign=callsign0(i+1:) - endif - if(issfx) then - tsfx=rof - k=nchar(tsfx(1:1)) - k=37*k + nchar(tsfx(2:2)) - k=37*k + nchar(tsfx(3:3)) - nv2=2 - i=index(callsign0,'/') - callsign=callsign0(:i-1) - endif - endif - endif - - return - end - diff --git a/libm65/getpfx1.f90 b/libm65/getpfx1.f90 new file mode 100644 index 000000000..b9e0522f5 --- /dev/null +++ b/libm65/getpfx1.f90 @@ -0,0 +1,96 @@ +subroutine getpfx1(callsign,k,nv2) + + character*12 callsign0,callsign,lof,rof + character*8 c + character addpfx*8,tpfx*4,tsfx*3 + logical ispfx,issfx,invalid + common/pfxcom/addpfx + include 'pfx.f90' + + callsign0=callsign + nv2=0 + iz=index(callsign,' ') - 1 + if(iz.lt.0) iz=12 + islash=index(callsign(1:iz),'/') + k=0 + c=' ' + if(islash.gt.0 .and. islash.le.(iz-4)) then +! Add-on prefix + c=callsign(1:islash-1) + callsign=callsign(islash+1:iz) + do i=1,NZ + if(pfx(i)(1:4).eq.c) then + k=i + go to 10 + endif + enddo + if(addpfx.eq.c) then + k=449 + go to 10 + endif + + else if(islash.eq.(iz-1)) then +! Add-on suffix + c=callsign(islash+1:iz) + callsign=callsign(1:islash-1) + do i=1,NZ2 + if(sfx(i).eq.c(1:1)) then + k=400+i + go to 10 + endif + enddo + endif + +10 if(islash.ne.0 .and.k.eq.0) then +! Original JT65 would force this compound callsign to be treated as +! plain text. In JT65v2, we will encode the prefix or suffix into nc1. +! The task here is to compute the proper value of k. + lof=callsign0(:islash-1) + rof=callsign0(islash+1:) + llof=len_trim(lof) + lrof=len_trim(rof) + ispfx=(llof.gt.0 .and. llof.le.4) + issfx=(lrof.gt.0 .and. lrof.le.3) + invalid=.not.(ispfx.or.issfx) + if(ispfx.and.issfx) then + if(llof.lt.3) issfx=.false. + if(lrof.lt.3) ispfx=.false. + if(ispfx.and.issfx) then + i=ichar(callsign0(islash-1:islash-1)) + if(i.ge.ichar('0') .and. i.le.ichar('9')) then + issfx=.false. + else + ispfx=.false. + endif + endif + endif + + if(invalid) then + k=-1 + else + if(ispfx) then + tpfx=lof + k=nchar(tpfx(1:1)) + k=37*k + nchar(tpfx(2:2)) + k=37*k + nchar(tpfx(3:3)) + k=37*k + nchar(tpfx(4:4)) + nv2=1 + i=index(callsign0,'/') + callsign=callsign0(:i-1) + callsign=callsign0(i+1:) + endif + if(issfx) then + tsfx=rof + k=nchar(tsfx(1:1)) + k=37*k + nchar(tsfx(2:2)) + k=37*k + nchar(tsfx(3:3)) + nv2=2 + i=index(callsign0,'/') + callsign=callsign0(:i-1) + endif + endif + endif + + return +end subroutine getpfx1 + diff --git a/libm65/getpfx2.f b/libm65/getpfx2.f deleted file mode 100644 index 26a46cce5..000000000 --- a/libm65/getpfx2.f +++ /dev/null @@ -1,24 +0,0 @@ - subroutine getpfx2(k0,callsign) - - character callsign*12 - include 'pfx.f' - character addpfx*8 - common/pfxcom/addpfx - - k=k0 - if(k.gt.450) k=k-450 - if(k.ge.1 .and. k.le.NZ) then - iz=index(pfx(k),' ') - 1 - callsign=pfx(k)(1:iz)//'/'//callsign - else if(k.ge.401 .and. k.le.400+NZ2) then - iz=index(callsign,' ') - 1 - callsign=callsign(1:iz)//'/'//sfx(k-400) - else if(k.eq.449) then - iz=index(addpfx,' ') - 1 - if(iz.lt.1) iz=8 - callsign=addpfx(1:iz)//'/'//callsign - endif - - return - end - diff --git a/libm65/getpfx2.f90 b/libm65/getpfx2.f90 new file mode 100644 index 000000000..d747e7f29 --- /dev/null +++ b/libm65/getpfx2.f90 @@ -0,0 +1,24 @@ +subroutine getpfx2(k0,callsign) + + character callsign*12 + include 'pfx.f90' + character addpfx*8 + common/pfxcom/addpfx + + k=k0 + if(k.gt.450) k=k-450 + if(k.ge.1 .and. k.le.NZ) then + iz=index(pfx(k),' ') - 1 + callsign=pfx(k)(1:iz)//'/'//callsign + else if(k.ge.401 .and. k.le.400+NZ2) then + iz=index(callsign,' ') - 1 + callsign=callsign(1:iz)//'/'//sfx(k-400) + else if(k.eq.449) then + iz=index(addpfx,' ') - 1 + if(iz.lt.1) iz=8 + callsign=addpfx(1:iz)//'/'//callsign + endif + + return +end subroutine getpfx2 + diff --git a/libm65/pfx.f b/libm65/pfx.f deleted file mode 100644 index 6a83179e1..000000000 --- a/libm65/pfx.f +++ /dev/null @@ -1,50 +0,0 @@ - parameter (NZ=339) !Total number of prefixes - parameter (NZ2=12) !Total number of suffixes - character*1 sfx(NZ2) - character*5 pfx(NZ) - - data sfx/'P','0','1','2','3','4','5','6','7','8','9','A'/ - data pfx/ - + '1A ','1S ','3A ','3B6 ','3B8 ','3B9 ','3C ','3C0 ', - + '3D2 ','3D2C ','3D2R ','3DA ','3V ','3W ','3X ','3Y ', - + '3YB ','3YP ','4J ','4L ','4S ','4U1I ','4U1U ','4W ', - + '4X ','5A ','5B ','5H ','5N ','5R ','5T ','5U ', - + '5V ','5W ','5X ','5Z ','6W ','6Y ','7O ','7P ', - + '7Q ','7X ','8P ','8Q ','8R ','9A ','9G ','9H ', - + '9J ','9K ','9L ','9M2 ','9M6 ','9N ','9Q ','9U ', - + '9V ','9X ','9Y ','A2 ','A3 ','A4 ','A5 ','A6 ', - + 'A7 ','A9 ','AP ','BS7 ','BV ','BV9 ','BY ','C2 ', - + 'C3 ','C5 ','C6 ','C9 ','CE ','CE0X ','CE0Y ','CE0Z ', - + 'CE9 ','CM ','CN ','CP ','CT ','CT3 ','CU ','CX ', - + 'CY0 ','CY9 ','D2 ','D4 ','D6 ','DL ','DU ','E3 ', - + 'E4 ','EA ','EA6 ','EA8 ','EA9 ','EI ','EK ','EL ', - + 'EP ','ER ','ES ','ET ','EU ','EX ','EY ','EZ ', - + 'F ','FG ','FH ','FJ ','FK ','FKC ','FM ','FO ', - + 'FOA ','FOC ','FOM ','FP ','FR ','FRG ','FRJ ','FRT ', - + 'FT5W ','FT5X ','FT5Z ','FW ','FY ','M ','MD ','MI ', - + 'MJ ','MM ', 'MU ','MW ','H4 ','H40 ','HA ', - + 'HB ','HB0 ','HC ','HC8 ','HH ','HI ','HK ','HK0 ', - + 'HK0M ','HL ','HM ','HP ','HR ','HS ','HV ','HZ ', - + 'I ','IS ','IS0 ', 'J2 ','J3 ','J5 ','J6 ', - + 'J7 ','J8 ','JA ','JDM ','JDO ','JT ','JW ', - + 'JX ','JY ','K ','KG4 ','KH0 ','KH1 ','KH2 ','KH3 ', - + 'KH4 ','KH5 ','KH5K ','KH6 ','KH7 ','KH8 ','KH9 ','KL ', - + 'KP1 ','KP2 ','KP4 ','KP5 ','LA ','LU ','LX ','LY ', - + 'LZ ','OA ','OD ','OE ','OH ','OH0 ','OJ0 ','OK ', - + 'OM ','ON ','OX ','OY ','OZ ','P2 ','P4 ','PA ', - + 'PJ2 ','PJ7 ','PY ','PY0F ','PT0S ','PY0T ','PZ ','R1F ', - + 'R1M ','S0 ','S2 ','S5 ','S7 ','S9 ','SM ','SP ', - + 'ST ','SU ','SV ','SVA ','SV5 ','SV9 ','T2 ','T30 ', - + 'T31 ','T32 ','T33 ','T5 ','T7 ','T8 ','T9 ','TA ', - + 'TF ','TG ','TI ','TI9 ','TJ ','TK ','TL ', - + 'TN ','TR ','TT ','TU ','TY ','TZ ','UA ','UA2 ', - + 'UA9 ','UK ','UN ','UR ','V2 ','V3 ','V4 ','V5 ', - + 'V6 ','V7 ','V8 ','VE ','VK ','VK0H ','VK0M ','VK9C ', - + 'VK9L ','VK9M ','VK9N ','VK9W ','VK9X ','VP2E ','VP2M ','VP2V ', - + 'VP5 ','VP6 ','VP6D ','VP8 ','VP8G ','VP8H ','VP8O ','VP8S ', - + 'VP9 ','VQ9 ','VR ','VU ','VU4 ','VU7 ','XE ','XF4 ', - + 'XT ','XU ','XW ','XX9 ','XZ ','YA ','YB ','YI ', - + 'YJ ','YK ','YL ','YN ','YO ','YS ','YU ','YV ', - + 'YV0 ','Z2 ','Z3 ','ZA ','ZB ','ZC4 ','ZD7 ','ZD8 ', - + 'ZD9 ','ZF ','ZK1N ','ZK1S ','ZK2 ','ZK3 ','ZL ','ZL7 ', - + 'ZL8 ','ZL9 ','ZP ','ZS ','ZS8 ','KC4 ','E5 '/ diff --git a/libm65/pfx.f90 b/libm65/pfx.f90 new file mode 100644 index 000000000..724b0f8a3 --- /dev/null +++ b/libm65/pfx.f90 @@ -0,0 +1,50 @@ + parameter (NZ=339) !Total number of prefixes + parameter (NZ2=12) !Total number of suffixes + character*1 sfx(NZ2) + character*5 pfx(NZ) + + data sfx/'P','0','1','2','3','4','5','6','7','8','9','A'/ + data pfx/ & + '1A ','1S ','3A ','3B6 ','3B8 ','3B9 ','3C ','3C0 ', & + '3D2 ','3D2C ','3D2R ','3DA ','3V ','3W ','3X ','3Y ', & + '3YB ','3YP ','4J ','4L ','4S ','4U1I ','4U1U ','4W ', & + '4X ','5A ','5B ','5H ','5N ','5R ','5T ','5U ', & + '5V ','5W ','5X ','5Z ','6W ','6Y ','7O ','7P ', & + '7Q ','7X ','8P ','8Q ','8R ','9A ','9G ','9H ', & + '9J ','9K ','9L ','9M2 ','9M6 ','9N ','9Q ','9U ', & + '9V ','9X ','9Y ','A2 ','A3 ','A4 ','A5 ','A6 ', & + 'A7 ','A9 ','AP ','BS7 ','BV ','BV9 ','BY ','C2 ', & + 'C3 ','C5 ','C6 ','C9 ','CE ','CE0X ','CE0Y ','CE0Z ', & + 'CE9 ','CM ','CN ','CP ','CT ','CT3 ','CU ','CX ', & + 'CY0 ','CY9 ','D2 ','D4 ','D6 ','DL ','DU ','E3 ', & + 'E4 ','EA ','EA6 ','EA8 ','EA9 ','EI ','EK ','EL ', & + 'EP ','ER ','ES ','ET ','EU ','EX ','EY ','EZ ', & + 'F ','FG ','FH ','FJ ','FK ','FKC ','FM ','FO ', & + 'FOA ','FOC ','FOM ','FP ','FR ','FRG ','FRJ ','FRT ', & + 'FT5W ','FT5X ','FT5Z ','FW ','FY ','M ','MD ','MI ', & + 'MJ ','MM ', 'MU ','MW ','H4 ','H40 ','HA ', & + 'HB ','HB0 ','HC ','HC8 ','HH ','HI ','HK ','HK0 ', & + 'HK0M ','HL ','HM ','HP ','HR ','HS ','HV ','HZ ', & + 'I ','IS ','IS0 ', 'J2 ','J3 ','J5 ','J6 ', & + 'J7 ','J8 ','JA ','JDM ','JDO ','JT ','JW ', & + 'JX ','JY ','K ','KG4 ','KH0 ','KH1 ','KH2 ','KH3 ', & + 'KH4 ','KH5 ','KH5K ','KH6 ','KH7 ','KH8 ','KH9 ','KL ', & + 'KP1 ','KP2 ','KP4 ','KP5 ','LA ','LU ','LX ','LY ', & + 'LZ ','OA ','OD ','OE ','OH ','OH0 ','OJ0 ','OK ', & + 'OM ','ON ','OX ','OY ','OZ ','P2 ','P4 ','PA ', & + 'PJ2 ','PJ7 ','PY ','PY0F ','PT0S ','PY0T ','PZ ','R1F ', & + 'R1M ','S0 ','S2 ','S5 ','S7 ','S9 ','SM ','SP ', & + 'ST ','SU ','SV ','SVA ','SV5 ','SV9 ','T2 ','T30 ', & + 'T31 ','T32 ','T33 ','T5 ','T7 ','T8 ','T9 ','TA ', & + 'TF ','TG ','TI ','TI9 ','TJ ','TK ','TL ', & + 'TN ','TR ','TT ','TU ','TY ','TZ ','UA ','UA2 ', & + 'UA9 ','UK ','UN ','UR ','V2 ','V3 ','V4 ','V5 ', & + 'V6 ','V7 ','V8 ','VE ','VK ','VK0H ','VK0M ','VK9C ', & + 'VK9L ','VK9M ','VK9N ','VK9W ','VK9X ','VP2E ','VP2M ','VP2V ', & + 'VP5 ','VP6 ','VP6D ','VP8 ','VP8G ','VP8H ','VP8O ','VP8S ', & + 'VP9 ','VQ9 ','VR ','VU ','VU4 ','VU7 ','XE ','XF4 ', & + 'XT ','XU ','XW ','XX9 ','XZ ','YA ','YB ','YI ', & + 'YJ ','YK ','YL ','YN ','YO ','YS ','YU ','YV ', & + 'YV0 ','Z2 ','Z3 ','ZA ','ZB ','ZC4 ','ZD7 ','ZD8 ', & + 'ZD9 ','ZF ','ZK1N ','ZK1S ','ZK2 ','ZK3 ','ZL ','ZL7 ', & + 'ZL8 ','ZL9 ','ZP ','ZS ','ZS8 ','KC4 ','E5 '/ diff --git a/libm65/pfxdump.f b/libm65/pfxdump.f deleted file mode 100644 index 1fbddc392..000000000 --- a/libm65/pfxdump.f +++ /dev/null @@ -1,13 +0,0 @@ - subroutine pfxdump(fname) - character*(*) fname - include 'pfx.f' - - open(11,file=fname,status='unknown') - write(11,1001) sfx - 1001 format('Supported Suffixes:'/(11('/',a1,2x))) - write(11,1002) pfx - 1002 format(/'Supported Add-On DXCC Prefixes:'/(15(a5,1x))) - close(11) - - return - end diff --git a/libm65/pfxdump.f90 b/libm65/pfxdump.f90 new file mode 100644 index 000000000..7587dbf72 --- /dev/null +++ b/libm65/pfxdump.f90 @@ -0,0 +1,13 @@ +subroutine pfxdump(fname) + character*(*) fname + include 'pfx.f90' + + open(11,file=fname,status='unknown') + write(11,1001) sfx +1001 format('Supported Suffixes:'/(11('/',a1,2x))) + write(11,1002) pfx +1002 format(/'Supported Add-On DXCC Prefixes:'/(15(a5,1x))) + close(11) + + return +end subroutine pfxdump