c   .    1    .    2    .    3    .    4    .    5    .    6    .    7    .    8
c
c  kreuzkorrelation zum vergleich zweier seismischer signale
c
c     max. ndim punkte je serie
c
c     created: <840118.1411>
c     Einfach-Version fuer SUN 1991
c     an NAS-Fortran angepasst Nov. 2005
      program twocrosp
      use nas_system
      parameter (nexp=21)
      parameter(ndim=2**nexp)
      real*4 x(ndim),y(ndim)
      character text1*72,text2*72,filein1*24,filein2*24
      logical repeat
      
      open(8,file='twocrosp.out')
      write(6,*) 'standard dual-signal coherency analysis'
      write(8,*) 'standard dual-signal coherency analysis'
      narg=iargc()
      if (narg.eq.0) then
      write(6,'(" File 1? ")')
      read(5,'(a)') filein1
      write(6,'(" File 2 ? ")')
      read(5,'(a)') filein2
      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
      nband=6
      else if (narg.eq.2) then
      call getarg(1,filein1)
      call getarg(2,filein2)
      nband=6
      else
      call getarg(1,filein1)
      call getarg(2,filein2)
      call getarg(3,text1)
      read(text1,*) nband
      endif
      repeat=(nband.gt.0)
      nband=abs(nband)

c  zwei zeitserien einlesen
      n1=ndim
      n2=ndim
      call input(filein1,text1,n1,dt1,tmin,tsec,x)
      write(6,*) ' header line: ',trim(text1)
      write(6,1) n1,dt1
    1 format(1x,i8," samples,   dt=",f10.5)
      if(filein2.eq.filein1) then
        text2=text1
        n2=n1
        dt2=dt1
        do j=1,ndim
          y(j)=x(j)
        enddo
      else
        call input(filein2,text2,n2,dt2,tmin,tsec,y)
        write(6,*) ' header line: ',trim(text2)
        write(6,1) n2,dt2
        if(dt1.ne.dt2) stop 12
      endif
      np=n1
      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)
      nmax=2**np2
      t=nmax*dt1
      write(6,111) np2,nmax
  111 format(' FFT 2**',i2,' = ',i8,' samples')
c  vorbereitung
      write(8,'(1x," 1. inputfile: ",a)') filein1
      write(8,'(1x,a)') text1
      write(8,1) n1,dt1
      write(8,'(1x," 2. inputfile: ",a)') filein2
      write(8,'(1x,a)') text2
      write(8,1) n2,dt2
c  beginn der analyse
      do i=np+1,nmax
        x(i)=0.
        y(i)=0.
      end do
      call trend(x,n1,ndim)
      call trend(y,n2,ndim)
      call sft(x,y,np2,-2,ndim)
      call trenn2(x,y,x,y,nmax,ndim)
      fa1=sqrt(2.*float(nmax)/float(n1))
      fa2=sqrt(2.*float(nmax)/float(n2))
      x(1)=0.
      y(1)=0.
      do i=2,nmax
        x(i)=x(i)*fa1
        y(i)=y(i)*fa2
      end do
      tshift=0.
      fak=1.
    2 call bandf(x,y,nmax,dt1,tshift,fak,nband,ndim)
      if(repeat) then
      write(6,'(/"gain and delay of signal 2? [0,0=end]: ")')
      read(5,*) fak,tshift
      write(8,3) fak,tshift
    3 format(//"signal 2 was divided by ",f10.3," and delayed by",f10.3)
      if(fak.ne.0.) goto 2
      endif
      close(8)
      write(6,'("results were stored in file twocrosp.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(24)
      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,n2,dt,tshift,fak,nband,ndim)
      character*72 form
      real x(ndim),y(ndim)
      complex cy
      w=45./atan2(1.,1.)
      n4=n2+2
      t=n2*dt
      domega=360./w/t
      xx=0.
      xy=0.
      yx=0.
      yy=0.
      ncoef=0
      xxtot=0.
      xytot=0.
      yxtot=0.
      yytot=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(/"  twocrosp: rms amplitudes in counts"/
     & "  use noise1, noise2 if both systems have the same response"/
     & "  use nois21 if system 1 has substantially lower noise than 2"/
     & "  equalize the responses before the analysis"/
     & "  use noise21 when coherency is > 0.5, otherwise use noise2"/
     & "  the bandwidth is ",f6.3," times the center frequency"//
     & " ctr.frq  amp 1  noise1   amp 2  noise2  nois21  phase2   gain2
     &     coh")
26    format("  (m)Hz   ------------   counts   -------------    deg "/)
      do i=2,jo
        j=jo+2-i
        jj=n4-j
        if(t/(j-1).ge.c1) then
          if(ncoef.ge.7) 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.
            xxtot=xxtot+xx
            xytot=xytot+xy
            yxtot=yxtot+yx
            yytot=yytot+yy
          else
            xx=xxtot
            xy=xytot
            yx=yxtot
            yy=yytot
            bmitte=1e12
            bbreit=1e12
          endif
          xyxy=xy*xy+yx*yx
          wxy=sqrt(xyxy)
          rn1=dim(xx,wxy)
          rn2=dim(yy,wxy)
          phi=-atan2(yx,xy)*w
          renz=min(wxy/sqrt(xx*yy),1.)
c          srld=max(wxy/bbreit,1e-11)
c          srld=10.*log10(srld)
           r21=sqrt(xx*(1.-renz**2)/renz**2)
      sxx=sqrt(xx)
      srn1=sqrt(rn1)
      syy=sqrt(yy)
      srn2=sqrt(rn2)

      amax=max(sxx,srn1,syy,srn2)
      if(amax.lt.1000) then
         form='(f7.3,1x,5(f7.3,1x),f7.1,1x,2(f7.4,1x))'
      else if(amax.lt.10000) then
         form='(f7.3,1x,5(f7.2,1x),f7.1,1x,2(f7.4,1x))'
      else if(amax.lt.100000) then
         form='(f7.3,1x,5(f7.1,1x),f7.1,1x,2(f7.4,1x))'
      else
         form='(f7.3,1x,5(f7.0,1x),f7.1,1x,2(f7.4,1x))'
      endif

            write(6,form) bmitte,sxx,srn1,
     &      syy,srn2,r21,phi,wxy/xx,renz
            write(8,form) bmitte,sxx,srn1,
     &      syy,srn2,r21,phi,wxy/xx,renz
          if(ncoef.lt.7) return
          xx=0.
          xy=0.
          yx=0.
          yy=0.
          co=c1
          c1=co*u
          ncoef=0
        end if
        omega=(j-1)*domega
        cy=cmplx(y(j),-y(jj))*cexp(cmplx(0.,-omega*tshift))
        yre=real(cy)/fak
        yim=-aimag(cy)/fak
        xx=xx+x(j)*x(j)+x(jj)*x(jj)
        xy=xy+x(j)*yre +x(jj)*yim
        yx=yx+x(j)*yim -x(jj)*yre
        yy=yy+yre *yre +yim  *yim
        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