Try to avoid a crash in freqcal.f90 when data is all zeros.

This commit is contained in:
Steven Franke 2019-05-11 09:36:15 -05:00
parent c393740b0a
commit ab1454a24c

View File

@ -41,6 +41,7 @@ subroutine freqcal(id2,k,nkhz,noffset,ntol,line)
endif endif
smax=0. smax=0.
s=0. s=0.
ipk=-99
do i=ia,ib do i=ia,ib
s(i)=real(cx(i))**2 + aimag(cx(i))**2 s(i)=real(cx(i))**2 + aimag(cx(i))**2
if(s(i).gt.smax) then if(s(i).gt.smax) then
@ -49,24 +50,31 @@ subroutine freqcal(id2,k,nkhz,noffset,ntol,line)
endif endif
enddo enddo
call peakup(s(ipk-1),s(ipk),s(ipk+1),dx) if(ipk.ge.1) then
fpeak=df * (ipk+dx) call peakup(s(ipk-1),s(ipk),s(ipk+1),dx)
ap=(fpeak/fs+1.0/(2.0*NFFT)) fpeak=df * (ipk+dx)
an=(fpeak/fs-1.0/(2.0*NFFT)) ap=(fpeak/fs+1.0/(2.0*NFFT))
sp=sum(id2((k-NFFT):k-1)*cmplx(cos(xi*ap),-sin(xi*ap))) an=(fpeak/fs-1.0/(2.0*NFFT))
sn=sum(id2((k-NFFT):k-1)*cmplx(cos(xi*an),-sin(xi*an))) sp=sum(id2((k-NFFT):k-1)*cmplx(cos(xi*ap),-sin(xi*ap)))
fpeak=fpeak+fs*(abs(sp)-abs(sn))/(abs(sp)+abs(sn))/(2*NFFT) sn=sum(id2((k-NFFT):k-1)*cmplx(cos(xi*an),-sin(xi*an)))
xsum=0. fpeak=fpeak+fs*(abs(sp)-abs(sn))/(abs(sp)+abs(sn))/(2*NFFT)
nsum=0 xsum=0.
do i=ia,ib nsum=0
if(abs(i-ipk).gt.10) then do i=ia,ib
xsum=xsum+s(i) if(abs(i-ipk).gt.10) then
nsum=nsum+1 xsum=xsum+s(i)
endif nsum=nsum+1
enddo endif
ave=xsum/nsum enddo
snr=db(smax/ave) ave=xsum/nsum
pave=db(ave) + 8.0 snr=db(smax/ave)
pave=db(ave) + 8.0
else
snr=-99.9
pave=-99.9
fpeak=-99.9
ferr=-99.9
endif
cflag=' ' cflag=' '
if(snr.lt.20.0) cflag='*' if(snr.lt.20.0) cflag='*'
n=n+1 n=n+1
@ -77,7 +85,7 @@ subroutine freqcal(id2,k,nkhz,noffset,ntol,line)
ncal=1 ncal=1
ferr=fpeak-noffset ferr=fpeak-noffset
write(line,1100) nhr,nmin,nsec,nkhz,ncal,noffset,fpeak,ferr,pave, & write(line,1100) nhr,nmin,nsec,nkhz,ncal,noffset,fpeak,ferr,pave, &
snr,cflag,char(0) snr,cflag,char(0)
1100 format(i2.2,':',i2.2,':',i2.2,i7,i3,i6,2f10.3,2f7.1,2x,a1,a1) 1100 format(i2.2,':',i2.2,':',i2.2,i7,i3,i6,2f10.3,2f7.1,2x,a1,a1)
900 return 900 return