c--snm.f
c--program to compute power spectral density (PSD) and
c--seismic noise magnitude (SNM) modified from plotp.f and
c--plotps2.f (subroutine freqb) on October 1999
c-----------------------------------------------------------------------
      program SNM
      implicit real*8(a-h,o-y),integer*4(i-n),logical(z)
      parameter(ndim=132000)
      dimension data(ndim)
      character*60 title
      character*40 fname
      character*20 fmt,dunit,funit
c--common block statements
      common // data
c
c initialize
      write(*,'(/a)') 'Program SNM - computed from meanfft file'
      write(*,'(a)')  '----------------------------------------'
      write(*,'(a$)') 'input filename [.fft] --> '
      read(*,'(a)') fname
      call fnext(fname,'fft')
      len = lstr(fname,-40)
c
c read fft file
      open(1,file = fname)
      read(1,'(a)') title
      read(1,'(3i10,g12.5,2i10)') inorm,no,nyq,delf,iw,lcb
      read(1,'(3a20)') fmt,dunit,funit
      read(1,fmt) (data(i),phase,i=1,nyq)
      close(1)
c
c--call the frequency band subroutine
      call freqb(no,nyq,delf,sum2)
c
c--write the result
      snm = dlog10(sum2)+2.5d0
      write(*,'(2(a,f10.5))') 'mean PSD = ',sum2,' SNM = ',snm
      write(*,*) 'done'
c
      stop
      end
c--------------------------------------------------------------------------
      subroutine freqb(no,nyq,delf,sum2)
c
c for the line amplitudes, find the average values in the various
c frequency bands. k1, k2 are the integers for the frequency bands f1,f2
c
      implicit real*8(a-h,o-y),integer*4(i-n),logical(z)
      parameter(ndim=132000)
      dimension data(ndim)
      common // data
      data fpsf/1.d0/,ps2/3600.d0/
      data f1,f2/6.d0,18.d0/
c
c--initialise, note that because data files are in sec,
c--frequency units for .fft files are assumed to be cycles/sec
      np = 2*(nyq-1)
      rno = dfloat(no)
      rnp = dfloat(np)
c
c--fnyq is the Nyquist frequency
c--we use this to scale the frequency range for the various bands
      fnyq = delf*rnp/2.0
      delt = 1.0/(2.*fnyq)
c
c change these limits to current units [f1p,f2p] by dividing by the
c frequency scaling factor ps2
      f1p = f1/ps2
      f2p = f2/ps2
c
c--convert to integers
      k1 = nint(f1p/delf)+1
      k2 = nint(f2p/delf)+1
      rm = float(k2-k1+1)
c
c--zero the PSD sum
      sum2 = 0.0
c
c--compute the statistics
      do 20 k = k1,k2
        fact = 2.0
        if(k.eq.1.or.k.eq.nyq) fact = 1.0
        s = data(k)
c--unnormalised mean PSD, compute squared PSD
        s2 = s*s/rno
        if(k.eq.nyq) s2 = s2/2.0
        sum2 = sum2+s2
  20  continue
c
c--do the appropriate normalisation
      sum2 = delt*sum2/rm
c
c--multiply by a factor of 2 to include negative frequencies
c--for the mean PSD
      sum2 = 2.0*sum2
c
      return
      end
c---------------------------------------------------------------------------
      subroutine fnext(fname,ext)
c
c changes filename extension to '.ext'
c
      implicit real*8(a-h,o-y),integer*4(i-n),logical(z)
      character*40 fname,fnamec
      character*3 ext,extc
      character*1 dum1(40),dum3(3)
      equivalence (fnamec,dum1),(extc,dum3)
c
c original length of filename is 'len0'
c find length of extension 'len3' by searching last 4 characters
      fnamec = fname
      extc = ext
      len0 = lstr(fnamec,-40)
      len3 = 0
      ltmp = 0
      do 5 i = 1,4
        l = len0-i+1
        if(l.eq.0) goto 6
        if(dum1(l).eq.'.') then
          len3 = ltmp+1
          goto 6
        endif
        ltmp = ltmp+1
   5  continue
c
c find length of main segment of the filename, back to a '/'
c if it exists as part of the filename
   6  len2 = 0
      do 10 i = 1,8
        l = len0-len3-i+1
        if(l.eq.0.or.dum1(l).eq.'//') goto 15
        len2 = len2+1
  10  continue
c
c find length of first segment (before any '/')
c new overall length of filename is len1+(new)len2+'.ext'
  15  len1 = len0-len3-len2
      len = len1+len2+4
c
c now assign last 4 characters of new filename
      dum1(len-3) = '.'
      dum1(len-2) = dum3(1)
      dum1(len-1) = dum3(2)
      dum1(len) = dum3(3)
      fname = fnamec
c
      return
      end
c----------------------------------------------------------------------------
      function lstr(name,n)
c
c function to return the length of a string. On input the string
c 'name' has a nominal length n; if n is negative, blank characters
c at the beginning of 'name' are first truncated.
c
      implicit real*8(a-h,o-y),integer*4(i-n),logical(z)
      character*1 name(1)
c
c--zero length case first
      if(n.eq.0) then
        lstr = 0
        return
      endif
      k = 0
      ld = n
c--ld is negative, we truncate leading blanks
      if(ld.lt.0) then
        ld = -ld
        j = 0
  10    j = j+1
        if(name(j).eq.' ') goto 10
        nb = j-1
c--there are nb blanks to be erased in front of 'name'
        if(nb.gt.0) then
          do 20 i = nb+1,ld
            name(i-nb) = name(i)
            name(i) = ' '
  20      continue
c--reduce length of string by number of blanks
          ld = ld-nb
        endif
      endif
c--ld is non-zero, find length of string
c--first set true length equal to nominal length
      k = ld+1
c--and decrease k until a non-blank character appears or beginning
c--of string is reached
  40  k = k-1
      if(k.gt.0.and.name(k).eq.' ') goto 40
      lstr = k
c
      return
      end