c   .    1    .    2    .    3    .    4    .    5    .    6    .    7    .    8
c  spectral analysis
c  this version reads ASL format
c
      program unispek
      use nas_system
      parameter (nexp=20)
      parameter(ndim=2**nexp)
      real*4 x(ndim),y(ndim)
      character text1*72,filein1*24,runs*10,name*24

      narg=iargc()
      if(narg.eq.0) then
       write(6,'(/"call: unispek file [opt [freq1 [freq2 [resol]]]]"/
     & "Spectrum from freq1 to freq2 with resolution resol"/
     & "opt=1 decibel scale, opt=2 remove trend, opt=4 taper"/
     & "opt=8 fixed standard scale (-20 to +100 decibels)"/
     & "parameter opt is the sum of all desired options, default=15")')
      stop
      endif
      call getarg(1,filein1)
      runs='15'
      if(narg.ge.2) call getarg(2,runs)
      read(runs,*) iopt
      logar=mod(iopt,2)
      iopt=iopt/2
      itrend=mod(iopt,2)
      iopt=iopt/2
      itap=mod(iopt,2)
      iopt=iopt/2
      iscal=mod(iopt,2)
      if(iscal.eq.1) logar=1
      runs='0.'
      if(narg.ge.3) call getarg(3,runs)
      read(runs,*) freq1
      runs='1e25'
      if(narg.ge.4) call getarg(4,runs)
      read(runs,*) freq2
      if(freq2.eq.0.) freq2=1e12
      runs='1e-25'
      if(narg.ge.5) call getarg(5,runs)
      read(runs,*) resol

c  zeitserie einlesen
      n1=ndim
      call input(filein1,text1,n1,dt,tmin,tsec,x)
      name=trim(filein1)//'.spc'

      power1=0.
      do j=1,n1
       power1=power1+x(j)**2
      enddo
      power1=power1/n1

      freq2=min(freq2,0.5/dt)
      resol=max(resol,1./n1/dt)
      np=n1
      open(8,file=name)
      write(6,'(a,a)') ' file header: ',trim(text1)
      write(8,'(a,a)') ' file header: ',trim(text1)
      write(6,1) n1,dt
    1 format(1x,i6," samples, dt=",f10.3)
c  serienlaenge 2**n fuer die schnelle fouriertransformation
      do while (2**np1 .lt. np)
       np1=np1+1
      end do
      np1=min(np1,nexp)
      nmax=2**np1
      write(6,4) nmax,np1
      if(nband.gt.0) write(8,4) nmax,np1
    4 format(" Fourier transform uses ",i8,"  ( =2**",i2," ) samples")
      t=nmax*dt
c  beginn der analyse
      do i=1,np
      y(i)=0.
      enddo
      do i=np+1,nmax
       x(i)=0.
       y(i)=0.
      end do

      if(itrend.eq.1) then
       call trend(x,n1)
       write(*,*) 'the linear trend was removed'
      endif
      

      goto 100
c  alte Version des Tapers
      if(itap.eq.1) then
       fn1=float(n1+1)
       fn2=fn1/2
       do i=1,n1
         fi=float(i)
         x(i)=x(i)*fi/fn2*(fn1-fi)/fn2
       enddo
      endif
      goto 101
c  neues Taper
 100  if(itap.eq.1) then
       n14=nint(np/4.)
       fak=sqrt(4./3.)
       f1=2.*atan(1.)/n14
       do i=1,n14
         fakt=fak*sin(f1*i)
         x(i)=x(i)*fakt
         x(np+1-i)=x(np+1-i)*fakt  
       enddo
       do i=n14+1,np-n14
         x(i)=x(i)*fak 
       enddo
       write(*,*) 'A taper was applied'
      endif
      write(*,*) 'the power of the time series is', power1

  101 call sft(x,y,np1,-2,ndim)
      fa1=2.*float(nmax)/float(n1)
      do i=1,nmax/2+1
       x(i)=(x(i)**2+y(i)**2)*fa1
       y(i)=0.
      end do

      power=0.
      do j=1,nmax/2.+1
       power=power+x(j)
      enddo
c      write(*,*) 'total spectral power 1 is ', power
      call spektrum(x,y,nmax,dt,freq1,freq2,resol,logar,iscal,name)
      close(8)
      stop
      end
c ---------------------------------------------------------------------
      subroutine spektrum(x,y,nmax,dt,freq1,freq2,aufl,logar,iscal,name)
      dimension x(nmax),y(nmax),g(1024)
      character*20 name
      n=nmax/2
      t=nmax*dt
      bin=1./t
      index1=nint(freq1*t+1.)
      index2=min(nint(freq2*t)+1,n)
      pih=2.*atan(1.)

c  Ausschnitt
      ind1=index1-1
      do j=index1,index2
       y(j-ind1)=x(j)
      enddo
      n=index2-ind1

      power=0.
      do j=1,n
       power=power+y(j)
      enddo
c      write(*,*) 'total spectral power 2 is ', power

c  Dezimieren falls mehr als 4095 Punkte
      idez=min(max(1,n/2048),1024)
      if(idez.gt.1024) stop 1024
      if (idez.gt.1) then
       fidez=idez
       bin=bin*fidez
       gsum=-1.
       do 3 j=1,idez
         g(j)=(cos((j-idez)/fidez*pih))**2
    3     gsum=gsum+2.*g(j)
       n=n/idez
       do 5 k=1,n-1
         jk=k*idez+1
         y(k)=g(idez)*y(jk)
         do 6 l=1,idez-1
    6      y(k)=y(k)+g(l)*(y(jk-idez+l)+y(jk+idez-l))
    5   continue
         do 7 j=n,2,-1
    7       y(j)=y(j-1)
         y(1)=0.
      endif


      power=0.
      do j=1,n
       power=power+y(j)
      enddo
c      write(*,*) 'total spectral power 3 is ', power

      write(6,8) n,idez,bin
    8 format(" spectrum was decimated to ",i8," samples"/
     & " idez=",i3," spectral resolution is ",f10.7," Hz")

c  Glaetten
      id=min(nint(aufl*t/idez),1024)
      if (id.gt.1) then
       fid=id
       do 13 j=1,id
   13     g(j)=(cos((j-id)/fid*pih))**2
       gsum=g(id)
       do j=1,id-1
         gsum=gsum+2.*g(j)
       enddo
c       write(*,*)
c       write(*,*) id, gsum
c       write(*,*)
       do j=1,id
         g(j)=g(j)/gsum
       enddo
       do 15 k=id,n-id+1
         x(k)=g(id)*y(k)
         do 16 l=1,id-1
   16      x(k)=x(k)+g(l)*(y(k-id+l)+y(k+id-l))
   15   continue
       n1=n+1
       do 14 k=1,id-1
       x(k)=y(k)
       x(n1-k)=y(n1-k)
       do 18 l=1,k-1
         x(n1-k)=x(n1-k)+y(n1-k+l)+y(n1-k-l)
   18     x(k)=x(k)+y(k-l)+y(k+l)
         x(n1-k)=x(n1-k)/(2*k-1)
         x(k)=x(k)/(2*k-1)
   14 continue
      else
         do 17 j=1,n
   17       x(j)=y(j)
      endif

      power=0.
      do j=1,n
       power=power+x(j)
      enddo
      write(*,*) 'total  spectral  power  is', power
      write(*,*) 'over a positive ',freq2-freq1,' Hz bandwidth'

      write(*,*) 'split into ',n, ' bins of ',bin,' Hz each'

c  convert to spectral power density
      do j=1,n
       x(j)=x(j)/bin
      enddo

c  Logarithmieren
      if (logar.eq.1) then
       dbrange=120.
       xmax=-480.
       xmin=-xmax
       x(1)=x(2)
       do j=1,n
         x(j)=10.*alog10(max(1e-48,x(j)))
         xmax=max(xmax,x(j))
         xmin=min(xmin,x(j))
       enddo
       xmin=max(xmin,xmax-dbrange)
       xm=(xmax+xmin)/2.

c  fixed decibel scale
       if(iscal.eq.1) then
         xmin=-40.
         xm=20.
         xmax=80.
       endif

c       write(*,*) xmin,xm,xmax
       do j=1,n
         x(j)=max(x(j),xmin)-xm
       enddo
      else
       xmax=0.
       do j=1,n
         xmax=max(xmax,x(j))
       enddo
       do j=1,n
         x(j)=2.*x(j)-xmax
       enddo
      endif

c Ausgeben
      write(*,*) 'writing spectrum file ',name
      fnorm=1.0
      write(8,9) n,fnorm*idez/t,freq1,freq2,bin
    9 format(i10," (5f12.6)          ",4f10.6)
      write(8,'(5f12.6)') (x(j),j=1,n)

c generate plot-parameter file
      open(9,file='winplot.par')
      if(iscal.eq.1) then
       write(9,100)
  100   format("0, 1, 24., 16., 0, 0, -60.")
      else
       write(9,104)
  104   format("0, 1, 24., 16., 0, 0, 1.")
      endif
      write(9,'(a)') name
      if(logar.eq.1) write(9,102) xmax, xm, xmin
  102 format("log power spectrum, top=",f7.1," db, center=",f7.1,
     &" dB, bottom=",f7.1,"   dB re 1 count-squared per Hz")
      if(logar.eq.0) write(9,103) xmax,xmax/2.,0.
  103 format("power spectrum: top=",f10.3,"   center=",f10.3,
     &   "   bottom=",f10.3,"   counts-squared per Hz")
      close(9)
      write(*,*) 'writing plot-parameter file winplot.par'
      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(21)
c      write(6,'("Beginn SFT")')
      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
c      write(6,'("Ende   SFT")')
      return
      end
c ---------------------------------------------------------------------
      subroutine trend(x,n)
      dimension x(n)
      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

      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