c************************************************************* c This program is named MAX.FOR. Its purpose is to produce c an output table of local maxima from an input table. It c requests a threshold which may be set to zero. c************************************************************* implicit double precision (a-h,o-z) character*40 infile dimension freq(500000),fr(500000) c c name input file from tty c write(*,1) 1 format(1x,' input file is (N<500 000): ',$) read(*,2) infile 2 format(a40) open(9,file=infile,status='old') c write(*,120) 120 format(1x,' output file is: ',$) read(*,2) infile open(11, file=infile, status='unknown') c write(*,2000) 2000 format(1x,' Set peak detection threshold: ',$) read(*,*) thresh 2001 continue read(9,*)dummy do 2005 i=1,500000 read(9,*,end=2111,err=2111)freq(i),fr(i) 2005 continue nsft=i goto 2200 2111 nsft=i-1 c Find peaks. c 2200 nsft3=nsft-3 iflag=0 istart=1 200 do 14 i=istart,nsft3 igot=i if (fr(i) .gt. fr(i+1) .and. iflag .eq. 1) go to 15 if (fr(i) .lt. fr(i+1)) iflag=1 14 continue close(9) close(11) stop 15 a1=fr(i-1) a2=fr(i) a3=fr(i+1) iflag=0 istart=igot a1a2=log(a1/a2) a1a3=log(a1/a3) biga=a1a2/log(a2/a3) rat=(1.+biga)/(1.-biga) f1=freq(igot-1) f2=freq(igot) f3=freq(igot+1) dfreq=(f2-f1) p=-rat/2. f0=f2+p*dfreq temp=a1a3*rat/8. a0=a2*exp(temp) if(a0.lt.thresh)goto 200 ap0=1.0d0/f0 fm=f0*1e6 am=a0*1e3 write (11,1500) fm,ap0,am 1500 format(1x,f18.8,' muHz',1x,f16.8,' s',1x,1pe10.3,' mma') goto 200 end