program plot_spectrum
c
c	standalone program to plot amplitude spectrum contained in
c	file produced by FIND (dunc@naic.edu - July 18, 2000)
c       N.B. Needs to be linked with local PGPLOT libraries
c
      implicit none
      integer mp,mc
      parameter(mp=2**22,mc=20)
      real fbin(mp),samp(mp),smin,smax,dm,ac
      real*8 tsamp,freq
      integer iargc,fold,npf,i
      character*80 sfile
      
      if (iargc().ne.1) stop 'usage: plot_spectrum <filename>'
      call getarg(1,sfile) 
      write(*,*) 'reading file...'
      call readspec(sfile,fold,samp,dm,ac,tsamp,npf)

      call pgbegin(0,'?',1,1)
      call pgscf(2)
      call pgsch(1.5)
      call pgvport(0.15,0.85,0.15,0.85)
      
      call pgadvance
      
      fold=1
      write(*,*) 'calculating frequencies...'
      do i=1,npf
        fbin(i)=real(freq(tsamp,npf,fold,i))
      enddo
      call minmax(samp,npf,smin,smax)

      write(*,*) 'plotting spectrum, wait a moment...'
      call pgwindow(fbin(1),fbin(npf),smin-abs(smin*0.1),smax*1.1)
      call pgbox('bcnst',0.0,0,'bcnst',0.0,0)
      call pgline(npf,fbin,samp)
      call pglabel('Frequency (Hz)','Amplitude',' ')
      
      write(*,*) 'all done' 
      call pgend
      end
c=============================================================================
      real*8 function freq(tsamp,npf,fold,k)
c=============================================================================
c
c     Returns the fluctuation frequency (Hz) of bin k in the Fourier
c     spectrum having npf points. Tsamp is the sampling interval
c     of the corresponding time domain data (seconds), whilst fold
c     refers to 1 plus the number of harmonic sums that have produced
c     the present spectrum. e.g fold=1 -> refers to the raw spectrum
c     fold=5 refers to 4 harmonic sums (16 harmonics).
c
      real*8 tsamp
      integer npf,fold,k
      freq=real(k)/(2.0*tsamp*npf*2.0**(fold-1))
      end
c=============================================================================
      integer function fbin(tsamp,npf,fold,freq)
c=============================================================================
c
c     Just the mathematical inverse of freq. Returns the bin number for
c     a given frequency.
c      
      real*8 tsamp
      integer npf,fold
      real freq
      fbin=real(freq)*(2.0*real(tsamp)*npf*2.0**(fold-1))
      end
c=============================================================================
      subroutine minmax(data,n,dmin,dmax)
      integer i,n
      real data(n),dmin,dmax
      dmin=+1.0e32
      dmax=-1.0e32
      do i=1,n
      dmin=min(dmin,data(i))
      dmax=max(dmax,data(i))
      enddo
      end
c=============================================================================
      subroutine readspec(sfile,fold,ssnr,dm,ac,tsamp,npf)
c=============================================================================
c
c     Reads in a power spectrum 
c 
c     99/07/12 - dunc@naic.edu -- added dm and ac to header info
c      
c=============================================================================
      implicit none
      character*80 sfile,fname
      real dm, ac
      integer fold,npf
      real*8 tsamp
      real ssnr(*)
      integer i,lun
C      byte bdat(2**23)
C      real scale,offset
      logical filex

      if (sfile.ne.' ') then
         fname=sfile
      else
         write(fname,'(''fold'',i1,''.spc'')')fold
      endif

      inquire(file=fname,exist=filex)
      if (.not.filex) stop 'File does not exist!'
      call glun(lun)
      open(unit=lun,file=fname,status='unknown',form='unformatted')
      read(lun) dm, ac ! new 
      read(lun) tsamp,npf,fold
      read(lun) (ssnr(i),i=1,npf)
C      read(lun) scale,offset
C      read(lun) (bdat(i),i=1,npf)
C      call reatfbin(ssnr,npf,bdat,scale,offset,-1)
      close(unit=lun)
      end
      
c==============================================================================
C nicked from pgplot
C*GRGLUN -- get a Fortran logical unit number (Sun/Convex-UNIX)
C+
      SUBROUTINE GLUN(LUN)
      INTEGER LUN
C
C Get an unused Fortran logical unit number. 
C Returns a Logical Unit Number that is not currently opened.
C After GRGLUN is called, the unit should be opened to reserve
C the unit number for future calls.  Once a unit is closed, it
C becomes free and another call to GRGLUN could return the same
C number.  Also, GRGLUN will not return a number in the range 1-9
C as older software will often use these units without warning.
C
C Arguments:
C  LUN    : receives the logical unit number, or -1 on error.
C--
C 12-Feb-1989 [AFT/TJP].
C DRL adapted to subroutine GLUN for use with stand-alone software
C 16-Jul-1993 @ JB
c==============================================================================
      INTEGER I
      LOGICAL QOPEN
C---
      DO 10 I=99,10,-1
          INQUIRE (UNIT=I,  OPENED=QOPEN)
          IF (.NOT.QOPEN) THEN
              LUN = I
              RETURN
          END IF
   10 CONTINUE
C none left
      STOP 'RAN OUT OF LUNs!' 
      END