diff --git a/CMakeLists.txt b/CMakeLists.txt index 00a1fc325..ea9055557 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1077,6 +1077,12 @@ target_link_libraries (msk144sim wsjt_fort wsjt_cxx) add_executable (msk144d2 lib/msk144d2.f90 wsjtx.rc) target_link_libraries (msk144d2 wsjt_fort wsjt_cxx) +add_executable (fmtave lib/fmtave.f90 wsjtx.rc) + +add_executable (fcal lib/fcal.f90 wsjtx.rc) + +add_executable (fmeasure lib/fmeasure.f90 wsjtx.rc) + add_executable (jt9 ${jt9_FSRCS} ${jt9_CXXSRCS} wsjtx.rc) if (${OPENMP_FOUND} OR APPLE) if (APPLE) @@ -1229,7 +1235,7 @@ install (TARGETS udp_daemon message_aggregator ) install (TARGETS jt9 jt65code qra64code qra64sim jt9code jt4code - msk144code wsprd + msk144code wsprd fmtave fcal fmeasure RUNTIME DESTINATION ${WSJT_BIN_DESTINATION} COMPONENT runtime BUNDLE DESTINATION ${WSJT_BIN_DESTINATION} COMPONENT runtime ) diff --git a/lib/fcal.f90 b/lib/fcal.f90 new file mode 100644 index 000000000..88ff99e50 --- /dev/null +++ b/lib/fcal.f90 @@ -0,0 +1,117 @@ +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 ' + 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 Call'/ & + ' (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 diff --git a/lib/fmeasure.f90 b/lib/fmeasure.f90 new file mode 100644 index 000000000..6613f301b --- /dev/null +++ b/lib/fmeasure.f90 @@ -0,0 +1,76 @@ +!------------------------------------------------------------------------------- +! +! 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) + real*8 fd(NZ),deltaf(NZ),r(NZ) + character infile*50 + character line*80 + + nargs=iargc() + if(nargs.ne.1) then + print*,'Usage: fmeasure ' + 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 diff --git a/lib/fmtave.f90 b/lib/fmtave.f90 new file mode 100644 index 000000000..e35dbe565 --- /dev/null +++ b/lib/fmtave.f90 @@ -0,0 +1,64 @@ +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 ' + 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 diff --git a/lib/freqcal.f90 b/lib/freqcal.f90 index ad3a09db6..4e2fcaef9 100644 --- a/lib/freqcal.f90 +++ b/lib/freqcal.f90 @@ -4,13 +4,12 @@ subroutine freqcal(id2,k,nkhz,noffset,ntol,line) integer*2 id2(0:NZ-1) real x(0:NFFT-1) real s(NH) - character line*80,cflag*1,callsign*6 + character line*80,cflag*1 complex cx(0:NH) equivalence (x,cx) data n/0/,k0/9999999/ save n,k0 - callsign=' ' if(k.lt.NFFT) go to 900 if(k.lt.k0) n=0 k0=k @@ -53,10 +52,8 @@ subroutine freqcal(id2,k,nkhz,noffset,ntol,line) ncal=1 ferr=fpeak-noffset write(line,1100) nhr,nmin,nsec,nkhz,ncal,noffset,fpeak,ferr,pave, & - snr,callsign,cflag,char(0) - write(61,1100) nhr,nmin,nsec,nkhz,ncal,noffset,fpeak,ferr,pave, & - snr,callsign,cflag,char(0) -1100 format(i2.2,':',i2.2,':',i2.2,i7,i3,i6,2f10.3,2f7.1,2x,a6,2x,a1,a1) + 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