c   .    1    .    2    .    3    .    4    .    5    .    6    .    7    .    8
c
c  calculates a time series with NLNM noise
c
c     max. ndim samples per series
c
c     Version Dez. 2005
      program lownoise
      use nas_system
      parameter (nexp=19)
      parameter(ndim=2**nexp)
      real*4 x(ndim),y(ndim)
      character text*50,outfile*12
      outfile='lownoise.sfe'
      text='synthetic NLNM noise from Nyqvist freq. to 3600 s '


      open(7,file='lownoise.par',status='old',err=201)
      goto 202
  201 write(6,*) 'File lownoise.par not found'
      write(6,*) 'This file must contain in one line each:'
      write(6,*) 'number of samples'
      write(6,*) 'sampling interval'
      write(6,*) 'free period of the seismometer'
      write(6,*) 'fraction of critical damping'
      write(6,*) 'generator constant'
      write(6,*) 'microvolts per digitizer count'
      write(6,*) 'number of bands per decade for the BANDNOIS analysis'
      write(6,*) '===================================================='
      write(6,*) 'The output is a time series in SEIFE format'
      write(6,*) 'representing noise according to the NLNM up to 7200 s'
      write(6,*) 'as it would be recorded by the seismometer and'
      write(6,*) 'digitizer specified above'
      stop
  202 continue
      read(7,*) np
      np=min(np,ndim)
      read(7,*) dt
      read(7,*) ts
      read(7,*) hs
      read(7,*) gs
      read(7,*) s
      read(7,*) nband
      close(7)
      zpi=8.*atan(1.)
      pih=zpi/4.

      np2=0
c  serienlaenge 2**n fuer die schnelle fouriertransformation
      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,' Punkte')

      nffth=nfft/2
      dauer=nfft*dt
      amp0=sqrt(0.5/dt)*gs/s*ts/zpi
      oms=zpi/ts
      om7rel=zpi/7200./oms
      x(1)=0.
      y(1)=0.
      do i=2, nffth
        ii=nfft+2-i
        per=dauer/(i-1)
        if(per.lt.7200.) then
          omrel=zpi/per/oms
          respq=(omrel-1./omrel)**2+4.*hs**2
          amp=amp0/sqrt(respq)*10**(fnlnm(per)/20.+5.85)
          if(per.gt.3600.) then
            amp=amp*sin(pih*(omrel/om7rel-1.))**2
          endif
          x(i)=amp*gaussrand(dum)
          y(i)=amp*gaussrand(dum)
        else
          x(i)=0.
          y(i)=0.
        endif
        x(ii)=x(i)
        y(ii)=-y(i)
      enddo
      x(nffth+1)=0.
      y(nffth+1)=0.

c      call output('lownoise1.re',text,nffth,dt,x)
c      call output('lownoise1.im',text,nffth,dt,y)

      call sft(x,y,np2,3,ndim)

      write(6,*)
      write(6,*) text
      call output(outfile,text,np,dt,x)
      open(8,file='lownoise.out')
      write(6,*)
      write(6,3)
      write(8,3)
    3 format(" BANDNOIS:"/" rms amplitudes in frequency bands of constant 
     &t width in a logarithmic scale")
      write(6,*) 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," microvolt"
      write(8,*) "seismometer: T = ",ts,"   h = ",hs,"   G = ",gs
      write(8,*) "digitizer: 1 count = ",s," microvolt"

c  vorbereitung
      write(8,'(1x," Outputfile: ",a)') outfile
      write(8,'(1x,a)') text
      write(8,1) np,dt
c  beginn der analyse
      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,'(/" BANDNOIS type analysis written to lownoise.out")')
c produce plot-parameter file
      open(9,file='winplot.par')
      write(9,7) 1,25.,16.,np
  7   format("0,",i3,",",f5.1,",",f5.1,",","  1,",i8,", 0.8")
      write(9,'(a)') trim(outfile)
      close(9)
      write(6,*) 'plot-parameter file winplot.par was generated'
      stop
      end
c ---------------------------------------------------------------------
      function gaussrand(dum)
      sum=0.
      do i=1,12
        call random_number(r)
        sum=sum+r
      enddo
      gaussrand=sum-6.
      return
      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 central frequency"//
     &       "   bandctr     bandw       amp     rtPSD       PSD      NL
     &NM       dif")
26    format("    (m)Hz      (m)Hz    nm/s^2   .../rtHz       dB       
     &dB        dB"/)
      bmit=1e12
      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
c          if(ncoef.ge.7.and.bmit.gt.0.00018) 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   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
c ========================================================
      subroutine output(name,text,n,dt,x)
      dimension x(n)
      character name*12,text*50,form*20
      xmax=0.
      do 2 j=1,n
    2 xmax=max(xmax,abs(x(j)))
      nvor=int(log10(max(xmax,2.))+1.)
      ndec=max(0,10-nvor)
      write(form,'(a,i1,a)') '(5f13.',ndec,')'
      write(6,5) trim(name)
    5 format(' writing synthetic NLNM noise into file ',a)
      open(8,file=name)
      write(8,'(a)') text
      write(8,1) n,form,dt,0.,0.
    1 format(i10,a20,f10.3,2f10.3)
      write(8,form) (x(j),j=1,n)
      close(8)
      return
      end