c   .    1    .    2    .    3    .    4    .    5    .    6    .    7..
c
c     linear combination of three signals, to fit a fourth signal
c     May be used for the orientation of borehole seismometers
c     April 1997 e.w. / updated May 2006
c

      program lincomb3
      use nas_system
      parameter(ndim=600000)
      dimension x(ndim),y(ndim),z(ndim),w(ndim),
     & wr(ndim),aik(3,3),rs(3),f(3)
c  wy: w nach Entfernung von x und y, usw.
      character inf1*48,inf2*48,inf3*48,inf4*48,form*20,par*12
      character ddir*48,outf*48,in4*48
      character*72 text1,text2,text3,text4


      write(6,*)
      write(6,2) ' LINCOMB3: solve w=a*x+b*y+c*z for a, b and c'
      write(6,2)  ' x=file1, y = file2, z=file3, w=file4'
      write(6,2)  ' the output file is the residue w-a*x-b*y-c*z.'
      narg=iargc()
      if(narg.lt.7) then
      write(6,*)
      write(6,2) ' call: lincomb datadir file1 file2 file3 file4 eps mt'
      write(6,2) ' eps can normally be set to zero. Small positive value
     &s'
      write(6,2) ' will stabilise the solution but deteriorate the fit.'
      write(6,2) ' remove mean if mt=1, remove trend if mt=2'
    2 format(a)
      stop
      endif

      call getarg(1,ddir)
      call getarg(2,inf1)
      call getarg(3,inf2)
      call getarg(4,inf3)
      call getarg(5,inf4)
      in4=inf4
      call getarg(6,par)
      read(par,*) eps
      call getarg(7,par)
      read(par,*) mn
      n=ndim
      inf1=trim(ddir)//'\'//trim(inf1)
      inf2=trim(ddir)//'\'//trim(inf2)
      inf3=trim(ddir)//'\'//trim(inf3)
      inf4=trim(ddir)//'\'//trim(inf4)
      
c  read data
      call input(inf1,text1,n,form,dt,tmin,tsec,x)
      call input(inf2,text2,n,form,dt,tmin,tsec,y)
      call input(inf3,text3,n,form,dt,tmin,tsec,z)
      call input(inf4,text4,n,form,dt,tmin,tsec,w)

      if(mn.eq.1) then
        call mean(x,n)
        call mean(y,n)
        call mean(z,n)
        call mean(w,n)
      else if (mn.eq.2) then
        call trend(x,n)
        call trend(y,n)
        call trend(z,n)
        call trend(w,n)
      endif

      do j=1,3
        rs(j)=0.
        do k=j,3
          aik(j,k)=0.
        enddo
      aik(j,j)=eps
      enddo

      do i=1,n
        aik(1,1)=aik(1,1)+x(i)**2
        aik(1,2)=aik(1,2)+x(i)*y(i)
        aik(1,3)=aik(1,3)+x(i)*z(i)
        aik(2,2)=aik(2,2)+y(i)**2
        aik(2,3)=aik(2,3)+y(i)*z(i)
        aik(3,3)=aik(3,3)+z(i)**2
         rs(1)= rs(1)+x(i)*w(i)
         rs(2)= rs(2)+y(i)*w(i)
         rs(3)= rs(3)+z(i)*w(i)
      enddo
      aik(2,1)=aik(1,2)
      aik(3,1)=aik(1,3)
      aik(3,2)=aik(2,3)

      call gauss(aik,3,3,rs,f)

      wq=0.
      wrq=0.
      do i=1,n
        wr(i)=w(i)-f(1)*x(i)-f(2)*y(i)-f(3)*z(i)
        wq = wq+ w(i)**2
        wrq=wrq+wr(i)**2
      enddo

      write(6,1) trim(inf4),trim(inf1),f(1),trim(inf2),f(2),
     & trim(inf3),f(3),wq,wrq,100.*sqrt(wrq/wq)
    1 format(/' Decomposition of signal w = ',a,':'/
     &        '   coefficient of series x = ',a,f13.6/
     &        '   coefficient of series y = ',a,f13.6/
     &        '   coefficient of series z = ',a,f13.6// 
     & ' Energy of series w       ',e10.3/
     & ' Energy of series w-xyz   ',e10.3/
     & ' rms residue              ',f10.6,' %')

      outf=trim(ddir)//'\'//trim(in4)//'.res'
      call output(outf,text1,n,form,dt,tmin,tsec,wr)
      stop
      end
c

      subroutine mean(x,n)
c  remove average
      dimension x(n)
      sum=0.
      do 1 j=1,n
    1 sum=sum+x(j)
      sum=sum/n
      do 2 j=1,n
    2 x(j)=x(j)-sum
      return
      end

c 
      subroutine trend(x,n)
c  remove linear trend
      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 1 j=1,n
      sx = sx+x(j)
    1 sjx = sjx+x(j)*float(j)
      a = (sx*beta-sjx*alpha)/det
      b = (sjx*gn-sx*alpha)/det
      do 2 j=1,n
    2 x(j) = x(j)-a-b*float(j)
      return
      end
c
      subroutine input(name,text,n,form,dt,tmin,tsec,x)
      dimension x(n)
      character name*(*),form*20,text*72
       call inpasc(name,text,n,form,dt,tmin,tsec,x)
      return
      end
c
      subroutine inpasc(name,text,n,form,dt,tmin,tsec,x)
      dimension x(n)
      character form*20,name*(*),text*72,pro*60
      open(7,file=name,status='old')
      read(7,'(a)') text
   21 read(7,'(a)') pro
      if(pro(1:1).eq.'%') goto 21
      read(pro,1) nn,form,dt,tmin,tsec
    1 format(i10,a20,3f10.3)
      write(6,22) trim(name),nn,n
   22 format(' file ',a,' wird gelesen:',i7,' /',i7,' punkte')
      if(n.lt.nn.or.nn.eq.0) then
c       write(6,'("punktzahl auf",i7," festgesetzt")') n
        nn=n
      endif
      n=nn
      if(form(1:4).eq.'frei') then
        read(7,*,err=25,end=23) (x(j),j=1,n)
      else
        read(7,form,end=23) (x(j),j=1,n)
      endif
      close(7)
      return
   23 n=j-1
      write(6,24) n
   24 format('Punktzahl ',i8,' durch Datenende bestimmt')
      close(7)
      return
   25 write(6,26) name,j
   26 format(' Input error reading file ',a,' at sample # ',i7)
      stop
      end
c
      subroutine output(name,text,n,form,dt,tmin,tsec,x)
      dimension x(n)
      character name*(*),text*72,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)
      i1=index(form,'(')
      i2=index(form,')')
      if(index(form,'frei').gt.0) then
        write(form,'(a,a,a)') '(5f13.',ndec,')'
      else if(i1.lt.1.or.i2.lt.i1+5) then
        write(6,'(a,a,a)') ' Ausgabeformat ',form,
     &' nicht interpretierbar'
        stop
      endif
      write(6,5) trim(name)
    5 format(/' writing residue into file ',a)
      open(8,file=name)
      write(8,'(a)') text
      write(8,1) n,form,dt,tmin,tsec
    1 format(i10,a20,f10.3,2f10.3)
      write(8,form) (x(j),j=1,n)
      close(8)
      return
      end

      subroutine gauss(aik,m,n,rs,f)
c  solve linear equations
      dimension aik(n,n),rs(n),f(n),h(14),imax(13)
      do 1401 j=1,m
      aikmax=0.
      do 1402 k=1,m
      h(k)=aik(j,k)
      if(abs(h(k)).le.aikmax) go to 1402
      aikmax=abs(h(k))
      index=k
 1402 continue
      h(m+1)=rs(j)
      do 1403 k=1,m
      q=aik(k,index)/h(index)
      do 1404 l=1,m
 1404 aik(k,l)=aik(k,l)-q*h(l)
 1403 rs(k)=rs(k)-q*h(m+1)
      do 1405 k=1,m
 1405 aik(j,k)=h(k)
      rs(j)=h(m+1)
 1401 imax(j)=index
      do 1406 j=1,m
      index=imax(j)
 1406 f(index)=rs(j)/aik(j,index)
      return
      end