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
