c   .    1    .    2    .    3    .    4    .    5    .    6    .    7    .    8
c
c  Kohaerenz von drei Signalen, Methode SLT (Sleeman Wettum Trampert)
c
c  Albuquerque, Mai 2006 EW
c
c     max. ndim punkte je serie
c
      program tricrosp
      use nas_system
      parameter (nexp=21)
      parameter(ndim=2**nexp)
      real*4 x(ndim),y(ndim),z(ndim),u(ndim)
      character filein1*24,filein2*24,filein3*24,
     &  text1*72,text2*72,text3*72
      logical repeat

      narg=iargc()

      write(*,*) 'Triple-signal coherency analysis according to Sleeman,
     & Wettum, Trampert'
      if (narg.eq.0) then
        write(6,'(" File 1? ")')
        read(5,'(a)') filein1
        write(6,'(" File 2 ? ")')
        read(5,'(a)') filein2
        write(6,'(" File 3 ? ")')
        read(5,'(a)') filein3
        write(6,5)
    5   format(' number of frequency bands per decade ?')
        read(5,*) nband
      else if (narg.eq.1) then
        call getarg(1,filein1)
        filein2=filein1
        filein3=filein1
        nband=6
      else if (narg.eq.3) then
        call getarg(1,filein1)
        call getarg(2,filein2)
        call getarg(3,filein3)
        nband=6
      else if (narg.eq.4) then
        call getarg(1,filein1)
        call getarg(2,filein2)
        call getarg(3,filein3)
        call getarg(4,text1)
        read(text1,*) nband
      else
        write(*,*) ' improper number of runstring arguments - stop'
        stop
      endif
      open(8,file='tricrosp.out')
      write(8,*) 'Triple-signal coherency analysis according to Sleeman,
     & Wettum, Trampert'

c  drei zeitserien einlesen
      n1=ndim
      n2=ndim
      n3=ndim
      call input(filein1,text1,n1,dt1,tmin,tsec,x)
      write(6,'(1x," read input file 1: ",a)') filein1
      write(8,'(1x," read input file 1: ",a)') filein1
      write(6,*) ' header line: ',trim(text1)
      write(8,*) ' header line: ',trim(text1)
      write(6,1) n1,dt1
      write(8,1) n1,dt1
    1 format(1x,i8," samples, dt=",f10.5)
      if(narg.eq.1) then
        text2=text1
        text3=text1
        n2=n1
        n3=n1
        dt2=dt1
        dt3=dt1
        do j=1,ndim
          y(j)=x(j)
          z(j)=x(j)
        enddo
      else
        call input(filein2,text2,n2,dt2,tmin,tsec,y)
        write(6,'(1x," read input file 2: ",a)') filein2
        write(8,'(1x," read input file 2: ",a)') filein2
        write(6,*) ' header line: ',trim(text2)
        write(8,*) ' header line: ',trim(text2)
        write(6,1) n2,dt2
        write(8,1) n2,dt2
        if(dt1.ne.dt2) stop 'inconsistent sampling interval - stop'
        call input(filein3,text3,n3,dt3,tmin,tsec,z)
        write(6,'(1x," read input file 3: ",a)') filein3
        write(8,'(1x," read input file 3: ",a)') filein3
        write(6,*) ' header line: ',trim(text2)
        write(8,*) ' header line: ',trim(text2)
        write(6,1) n3,dt3
        write(8,1) n3,dt3
        if(dt1.ne.dt3) stop 'inconsistent sampling interval - stop'
      endif
      np=min(n1,n2,n3)
      do i=1,np
        u(i)=0.
      enddo

c  Taper
      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
        y(i)=y(i)*fakt
        z(i)=z(i)*fakt
        x(np+1-i)=x(np+1-i)*fakt
        y(np+1-i)=y(np+1-i)*fakt
        z(np+1-i)=z(np+1-i)*fakt
      enddo
      do i=n14+1,np-n14
        x(i)=x(i)*fak
        y(i)=y(i)*fak
        z(i)=z(i)*fak
      enddo
      write(6,*) 'The signals were tapered'
      write(8,*) 'The signals were tapered'

c  serienlaenge 2**n fuer die schnelle fouriertransformation
      np2=0
      do while (2**np2 .lt. np)
        np2=np2+1
      end do
      np2=min(np2,nexp)
      nmax=2**np2
      t=nmax*dt1
      write(6,111) np2,nmax,np
      write(8,111) np2,nmax,np
  111 format(" FFT 2**",i2," = ",i8," points of which ",i8," are data")

c  beginn der analyse
      do i=np+1,nmax
        x(i)=0.
        y(i)=0.
        z(i)=0.
        u(i)=0.
      end do
      call trend(x,n1,ndim)
      call trend(y,n2,ndim)
      call trend(z,n3,ndim)
      call sft(x,y,np2,-2,ndim)
      call trenn2(x,y,x,y,nmax,ndim)
      call sft(z,u,np2,-2,ndim)
      call trenn2(z,u,z,u,nmax,ndim)
      fa=sqrt(2.*float(nmax)/float(np))
      x(1)=0.
      y(1)=0.
      z(1)=0.
      do i=2,nmax
        x(i)=x(i)*fa
        y(i)=y(i)*fa
        z(i)=z(i)*fa
      end do
    2 call bandf(x,y,z,nmax,dt1,nband,ndim)
      close(8)
      write(6,*)
      write(6,'("results were stored in file tricrosp.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,y,z,n2,dt,nband,ndim)
      character*72 form
      real x(ndim),y(ndim),z(ndim)
      complex c(3),cc(3,3),cnull,t2,t3,ee1,ee2
      cnull=cmplx(0.,0.)
      w=45./atan2(1.,1.)
      n4=n2+2
      t=n2*dt
      do 30 k=1,3
        do 30 kk=1,3
   30     cc(k,kk)=cnull
      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)
        write(6,26)
        write(8,25)
        write(8,26)
25      format(/"  tricrosp: rms amplitudes in counts. amp1 is the coher
     &ent signal in ch.1"/
     &  "  bandm    amp1      t2      t3      n1      n2      n3    coh21
     &1    coh31")
26    format(
     &  "  (m)Hz     cnt     rel     rel     cnt     cnt     cnt      sqd
     &d      sqd"/)

c  Testausdruck
c          write(6,*) 'ccc'
c          write(6,'(9(2x,2i3))') ((k,kk,kk=1,3),k=1,3)

      do i=2,jo
        j=jo+2-i
        jj=n4-j        
        if(t/(j-1).ge.c1) then
          bmitte=1./sqrt(co*c1)
	  if(bmitte.lt.0.01) bmitte=bmitte*1000.
          bbreit=1./co-1./c1
	  if(bbreit.lt.0.01) bbreit=bbreit*1000.

coherency squared
          coh21=real(cc(2,1)*cc(1,2)/(cc(1,1)*cc(2,2)))
          coh31=real(cc(3,1)*cc(1,3)/(cc(1,1)*cc(3,3)))

c  calculate complex gain (relative to the x channel)
          t2=cc(2,3)/cc(1,3)*coh21
          t2q=real(t2*conjg(t2))
          t3=cc(3,2)/cc(1,2)*coh31
          t3q=real(t3*conjg(t3))

c  calculate signal and noise power from the cross-spectral coefficients
c  statistical model: signal = transfer fct * input + noise
          ee1=cc(2,1)*cc(1,3)/cc(2,3) ! first  estimate of coh. power
          ee2=cc(3,1)*cc(1,2)/cc(3,2) ! second estimate of coh. power
          eer=real((ee1+ee2)/2.)      ! avg.   estimate of coh. power
          eei=aimag((ee1-ee2)/2.)     ! should be zero
          qn1=real(cc(1,1))-eer       ! noise power 1
          qn2=real(cc(2,2))-t2q*eer   ! noise power 2
          qn3=real(cc(3,3))-t3q*eer   ! noise power 3

          amax=max(eer,qn1,qn2,qn3)
          if(amax.lt.1000) then
            form='(f7.3,1x,f7.3,1x,2(f7.3,1x),3(f7.3,1x),f8.4,1x,f8.4)'
          else if(amax.lt.10000) then
            form='(f7.3,1x,f7.2,1x,2(f7.3,1x),3(f7.2,1x),f8.4,1x,f8.4)'
          else if(amax.lt.100000) then
            form='(f7.3,1x,f7.1,1x,2(f7.3,1x),3(f7.1,1x),f8.4,1x,f8.4)'
          else
            form='(f7.3,1x,f7.0,1x,2(f7.3,1x),3(f7.0,1x),f8.4,1x,f8.4)'
          endif

c  convert to amplitude
          erms=sqrt(max(0.,eer))
          rn1=sqrt(max(0.,qn1))
          rn2=sqrt(max(0.,qn2))
          rn3=sqrt(max(0.,qn3))

          write(6,form) bmitte,erms,cabs(t2),cabs(t3),rn1,rn2,rn3
     &   ,coh21,coh31
          write(8,form) bmitte,erms,cabs(t2),cabs(t3),rn1,rn2,rn3
     &   ,coh21,coh31

c  Testausdruck
c          write(6,'(9f8.3)') ((cabs(cc(k,kk)),kk=1,3),k=1,3)

          if(ncoef.lt.7) return
          
          do 31 k=1,3
            do 31 kk=1,3
   31         cc(k,kk)=cnull
          co=c1
          c1=co*u
          ncoef=0
        end if

c  calculate cross-spectral coefficients
        c(1)=cmplx(x(j),x(jj))
        c(2)=cmplx(y(j),y(jj))
        c(3)=cmplx(z(j),z(jj))
        
        do 32 k=1,3
          do 32 kk=1,3
   32       cc(k,kk)=cc(k,kk)+c(k)*conjg(c(kk))
        
        ncoef=ncoef+1
      end do
      return
      end
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