c   .    1    .    2    .    3    .    4    .    5    .    6    .    7    .    8
c
c  Bandspektrum
c  this version reads ASL format
c
c     max. ndim punkte je serie
c
c     Version Dez. 2005
      program bandnois
      use nas_system
      parameter (nexp=21)
      parameter(ndim=2**nexp)
      real*4 x(ndim),y(ndim)
      character text*72,infile*24

      narg=iargc()
      if (narg.eq.0) then
      write(6,*) 'Calculate noise in a fraction-of-decade bandwidth'
      write(6,*) 'Seismometer response: geophone or VBB type'
      write(6,'(" Filename? ")')
      read(5,'(a)') infile
      write(6,'(" how many microvolts per count ? ")')
      read(5,*) s
      write(6,'(" free period ? ")')
      read(5,*) ts
      write(6,'(" damping ? ")')
      read(5,*) hs
      write(6,'(" generator constant ? ")')
      read(5,*) gs
      write(6,'(" how many frequency bands per decade ? ")')
      read(5,*) nband
      else if (narg.eq.6) then
      call getarg(1,infile)
      call getarg(2,text)
      read(text,*) s
      call getarg(3,text)
      read(text,*) ts
      call getarg(4,text)
      read(text,*) hs
      call getarg(5,text)
      read(text,*) gs
      call getarg(6,text)
      read(text,*) nband
      else
      write(6,*) 'call BANDNOIS (interactive) or BANDNOIS filename micro
     &volt_per_count ts hs gs nband'
      stop
      endif
      zpi=8.*atan(1.)

      open(8,file='bandnois.out')
      write(6,*)
      write(6,3)
      write(8,3)
    3 format(" BANDNOIS:"/" rms amplitudes in frequency bands with a 
     &constant relative width")
c  read time series
      np=ndim
      call input(infile,text,np,dt,tmin,tsec,x)
      write(6,*) ' header line: ',trim(text)
      write(6,1) np,dt
    1 format(1x,i8," punkte, dt=",f10.5)
      write(6,*) "Seismometer: T = ",ts,"   h = ",hs,"   G = ",gs
      write(6,*) "Digitizer: 1 count = ",s," Mikrovolt"
      write(8,*) "Seismometer: T = ",ts,"   h = ",hs,"   G = ",gs
      write(8,*) "Digitizer: 1 count = ",s," Mikrovolt"

      np2=0
c  determine next higher power of 2 for FFT 
      do while (2**np2 .lt. np)
        np2=np2+1
      end do
      np2=min(np2,nexp)
      nfft=2**np2
      write(6,111) np2,nfft
  111 format(' FFT 2**',i2,' = ',i8,' samples')
      write(8,'(1x," Inputfile: ",a)') infile
      write(8,'(1x,a)') text
      write(8,1) np,dt
      fa1=sqrt(2.*float(nfft)/float(np))*s*zpi/ts/gs
      do i=1,np
        x(i)=x(i)*fa1
        y(i)=0.
      enddo
      do i=np+1,nfft
        x(i)=0.
        y(i)=0.
      end do
      call trend(x,np,ndim)
      call sft(x,y,np2,-2,ndim)
      call trenn2(x,y,x,y,nfft,ndim)
    2 call bandf(x,nfft,dt,nband,ndim,ts,hs)
      close(8)
      write(6,'(/"Results were saved as file bandnois.out")')
      stop
      end
c ---------------------------------------------------------------------
      subroutine sft(x,y,n,is,ndim)
c     schnelle fouriertransformation der 2**n komplexen werte (x,y)
c     is negativ - in den frequenzbereich / positiv - in den zeitbereich
c     normierung - abs(is) =1 - ohne / 2 - mit 1/ng / 3 - mit 1/sqrt(ng)
c                            4 - ohne normierung,ohne kontrollausdruck
      real*4 x(ndim),y(ndim)
      integer   zh(30)
      piz=8.*atan(1.)
c     tabelle der zweierpotenzen
      zh(1)=1
      do 1 l=1,n
    1 zh(l+1)=2*zh(l)
      ng=zh(n+1)
      gn=1./float(ng)
c     kernprogramm.dreifache schleife ueber schritt/index/teilserie
      do 2 m=1,n
      nar=zh(m)
      lar=ng/nar
      larh=lar/2
      alpha =  piz/float(isign(lar,is))
      do 3 jr=1,larh
      beta=alpha*float(jr-1)
      excos = cos(beta)
      exsin = sin(beta)
      ja=jr-lar
      do 4 nr=1,nar
      ja=ja+lar
      jb=ja+larh
      zx = x(ja)-x(jb)
      zy = y(ja)-y(jb)
      x(ja) = x(ja)+x(jb)
      y(ja) = y(ja)+y(jb)
      x(jb) = zx*excos-zy*exsin
      y(jb) = zx*exsin+zy*excos
    4 continue
    3 continue
    2 continue
c     normierung
      if(iabs(is).eq.1.or.iabs(is).eq.4) go to 10
      if(iabs(is).eq.3) gn=sqrt(gn)
    5 do 6  l=1,ng
      y(l) = y(l)*gn
    6 x(l)=x(l)*gn
c     umordnung nach "bitreversed" indizes
   10 do 7 j=1,ng
      js=j-1
      k=1
      nny=n+1
      do 8 ny=1,n
      nny=nny-1
      if(js.lt.zh(nny)) go to 8
      js=js-zh(nny)
      k=k+zh(ny)
    8 continue
      if(j-k) 9,7,7
    9 zx = x(j)
      zy = y(j)
      x(j) = x(k)
      y(j) = y(k)
      x(k) = zx
      y(k) = zy
    7 continue
      return
      end
c ---------------------------------------------------------------------
      subroutine trenn2(x,y,xt,yt,ng,ndim)
      real x(ndim),y(ndim),xt(ndim),yt(ndim)
      ngh = ng/2
      ngh1 = ngh+1
      xre = x(1)
      yt(1)=y(1)
      xt(1) = xre
      do 30 l=2,ngh
      m=ng-l+2
      xre=(x(l)+x(m))*0.5
      xim=(y(m)-y(l))*0.5
      yre=(y(l)+y(m))*0.5
      yim=(x(l)-x(m))*0.5
      xt(l)=xre
      xt(m)=xim
      yt(l)=yre
   30 yt(m)=yim
      xre = x(ngh1)
      yt(ngh1) = y(ngh1)
      xt(ngh1) = xre
      return
      end
c ---------------------------------------------------------------------
      subroutine trend(x,n,ndim)
      dimension x(ndim)
      gn = float(n)
      alpha = 0.5*gn*(gn+1.)
      beta = (2.*gn+1.)*(gn+1.)*gn/6.
      det = gn*beta-alpha*alpha
      sx = 0.
      sjx = 0.
      do 1001 j=1,n
      sx = sx+x(j)
 1001 sjx = sjx+x(j)*float(j)
      a = (sx*beta-sjx*alpha)/det
      b = (sjx*gn-sx*alpha)/det
      do 1002 j=1,n
 1002 x(j) = x(j)-a-b*float(j)
      return
      end
c----------------------------------------------------------------------
      subroutine bandf(x,nfft,dt,nband,ndim,ts,hs)
      character*72 form
      real x(ndim)
      zpi=8.*atan(1.)
      tsd=1000.
      w=360./zpi
      oms=zpi/ts
      n4=nfft+2
      t=nfft*dt
      domega=zpi/t
      xx=0.
      ncoef=0
      u=10.**(1./nband)
      ico=int(1024.+nband*alog10(2.*dt))-1023
      c0=ico/float(nband)
      co=10.**c0
      c1=u*co
      jo=t/co
      bbrel=10.**(0.5/nband)-10.**(-0.5/nband)
        write(6,25) bbrel
        write(6,26)
        write(8,25) bbrel
        write(8,26)
25    format(
     &  " the bandwidth is ",f5.3," times the center frequency"//
     &       "  bandctr  bandwidth     ampli   rms per       psd      NL
     &NM       dif")
26    format("    (m)Hz      (m)Hz    nm/s^2   root Hz        dB       
     &dB        dB"/)
      do i=2,jo
        j=jo+2-i
        jj=n4-j
        omega=(j-1)*domega
        omrel=omega/oms
        respq=(omrel-1./omrel)**2+4.*hs**2
        if(t/(j-1).ge.c1) then
          if(ncoef.ge.7) then
            bmitte=1./sqrt(co*c1)
            bmit=bmitte
	    if(bmitte.lt.0.01) bmitte=bmitte*1000.
            bbreit=1./co-1./c1
            breit=bbreit
	    if(bbreit.lt.0.01) bbreit=bbreit*1000.
          else
            return
          endif
      sxx=sqrt(xx)
      rld=sqrt(xx/breit)
      rldb=20.*alog10(max(rld,1e-24))-120.
      pdb=fnlnm(1./bmit)
      amax=sxx
      if(amax.lt.1000) then
         form='(4f10.3,3f10.1)'
      else if(amax.lt.10000) then
         form='(4f10.2,3f10.1)'
      else if(amax.lt.100000) then
         form='(4f10.1,3f10.1)'
      else
         form='(4f10.0,3f10.1)'
      endif
      write(6,form) bmitte,bbreit,tsd*sxx,tsd*rld,rldb,pdb,rldb-pdb
      write(8,form) bmitte,bbreit,tsd*sxx,tsd*rld,rldb,pdb,rldb-pdb
      if(ncoef.lt.7) return
        xx=0.
        co=c1
        c1=co*u
        ncoef=0
      end if
        xx=xx+(x(j)*x(j)+x(jj)*x(jj))*respq
        ncoef=ncoef+1
      end do
      return
      end
c ----------------------------------------------------------------------
c
      subroutine input(name,text,n,dt,tmin,tsec,x)
      dimension x(n)
      character iform*20,name*24,text*72,zeile*72,code1*12,code2*12
      logical seife

      write(6,*) 'Opening file ',trim(name)
        open(7,file=name,status='old')
        read(7,'(a)') zeile
        read(7,'(a)') zeile
        seife=index(zeile,'%').gt.0.or.index(zeile,'(').gt.0
        if(seife) then
          write(6,*) 'file ',trim(name),' assumed to be in SEIFE format'
        else
          write(6,*) 'file ',trim(name),' assumed to be in ASL format'
        endif
        close(7)
        
      if(seife) then
c read data in SEIFE format      
        open(7,file=name,status='old')
        read(7,'(a)') text
        write(6,*) ' header: ',trim(text)
   21   read(7,'(a)') zeile
        if(zeile(1:1).eq.'%') goto 21
   20   read(zeile,1) nn,iform,dt,tmin,tsec
    1   format(i10,a20,3f10.3)
        if(nn.gt.n) then
          write(*,*) 'sorry, too many data points. can handle only ', n
          write(*,*) 'for more, edit and recompile the source code'
        else
          n=nn
        endif
        write(6,22) n,trim(name)
   22   format(' reading',i7,' samples from file ',a)
        if(iform(1:3).eq.'fre'.or.index(iform,'i').gt.0) then
          read(7,*,err=25,end=23) (x(j),j=1,n)
        else
          read(7,iform,err=25,end=23) (x(j),j=1,n)
        endif
        close(7)
        return
      else
c read data in ASL format (such as written by Quanterra's Cimarron)
        open(7,file=name,status='old')
        read(7,'(a)') text
        write(6,*) 'header: ',trim(text)
        read(text,*) code1,code2,year,day,thr,tmin,tsec,t100,srate,nn
        if(nn.gt.n) then
          write(*,*) 'sorry, too many data points. can handle only ', n
          write(*,*) 'for more, edit and recompile the source code'
        else
          n=nn
        endif
        dt=1./srate
        tmin=tmin+60.*thr
        tsec=tsec+t100/100.
        read(7,*,err=25,end=23) (x(j),j=1,n)
        write(*,*) n,' samples read from file ',trim(name)
        close(7)
        return
      endif
      
   23 n=j-1
      write(6,24) n
      write(3,24) n
   24 format('end of file after ',i8,' samples')
      close(7)
      return
   25 write(6,26) j
      write(3,26) j
   26 format(' Input error at sample # ',i8)
      stop
      end

c ----------------------------------------------------------------------

c   evaluation of the Low Noise Model

      FUNCTION fnlnm(p)
      common /nlnm/ per(21),a(21),b(21)
c  
c   USGS New Low Noise Model (Peterson 1993)
c  
      data per(1),a(1),b(1)    / 0.1,-162.36,5.64 /
      data per(2),a(2),b(2)    / 0.17,-166.7,0 /
      data per(3),a(3),b(3)    / 0.4,-170,-8.3 /
      data per(4),a(4),b(4)    / 0.8,-166.4,28.9 /
      data per(5),a(5),b(5)    / 1.24,-168.6,52.48 /
      data per(6),a(6),b(6)    / 2.4,-159.98,29.81 /
      data per(7),a(7),b(7)    / 4.3,-141.1,0 /
      data per(8),a(8),b(8)    / 5,-71.36,-99.77 /
      data per(9),a(9),b(9)    / 6,-97.26,-66.49 /
      data per(10),a(10),b(10) / 10,-132.18,-31.57 /
      data per(11),a(11),b(11) / 12,-205.27,36.16 /
      data per(12),a(12),b(12) / 15.6,-37.65,-104.33 /
      data per(13),a(13),b(13) / 21.9,-114.37,-47.1 /
      data per(14),a(14),b(14) / 31.6,-160.58,-16.28 /
      data per(15),a(15),b(15) / 45,-187.5,0 /
      data per(16),a(16),b(16) / 70,-216.47,15.7 /
      data per(17),a(17),b(17) / 101,-185,0 /
      data per(18),a(18),b(18) / 154,-168.34,-7.61 /
      data per(19),a(19),b(19) / 328,-217.43,11.9 /
      data per(20),a(20),b(20) / 600,-258.28,26.6 /
      data per(21),a(21),b(21) / 10000,-346.88,48.75 /

      IF (p.lt.0.1) THEN
c        write(6,*) ' NLNM undefined, OLNM used'
        fnlnm=-168.
      ELSE IF (p.le.100000.) then
        do k=1,20
          IF (p.lt.per(k+1)) goto 20
        enddo
   20   fnlnm=a(k)+b(k)*LOG10(p)
      ELSE
c        write(6,*) ' NLNM undefined'
        fnlnm=0.
      ENDIF
      END