c   .    1    .    2    .    3    .    4    .    5    .    6    .    7..
c
c     linear combination of six signals, to fit a seventh signal
c     May be used for the determination of a moment tensor
c     Oct. 2009 EW
c

      program linreg6
      use nas_system
      parameter(ndim=16500)
      dimension r(ndim),s(ndim),t(ndim),syn(ndim)
      dimension x(ndim),y(ndim),z(ndim),w(ndim),
     & wr(ndim),aik(6,6),rs(6),f(6)
      character*48 inf1,inf2,inf3,inf4,inf5,inf6,inf7,synf,resf,ddir
      character*20 form 
      character*72 text1,text2,text3,text4,text5,text6,text7


      write(6,2)
      write(6,2) ' LINREG6: solve w = a*r +b*s +c*t +d*x +e*y +f*z' 
      write(6,2) ' for a, b, c, d, e,f '
      write(6,2) ' r=file1, s = file2, etc.'
      write(6,2) ' the output file is the residue w -a*r -b*s -...'
      write(6,2) ' The parameterfile linreg.par must contain'
      write(6,2) ' name of the data directory'
      write(6,2) ' names of the r, s, t, x, y, z files, one line each'
      write(6,2) ' name of the w file'
      write(6,2) ' eps: stabilization weight, can normally be zero'
      write(6,2) ' mn: remove mean if =1, remove trend if =2'
    2 format(a)
      
      open(11,file='linreg.par',status='old')
      read(11,*) ddir
      read(11,*) inf1
      read(11,*) inf2
      read(11,*) inf3
      read(11,*) inf4
      read(11,*) inf5
      read(11,*) inf6
      read(11,*) inf7
      read(11,*) eps
      read(11,*) mn

      n=ndim
      inf1=trim(ddir)//'\'//trim(inf1)
      inf2=trim(ddir)//'\'//trim(inf2)
      inf3=trim(ddir)//'\'//trim(inf3)
      inf4=trim(ddir)//'\'//trim(inf4)
      inf5=trim(ddir)//'\'//trim(inf5)
      inf6=trim(ddir)//'\'//trim(inf6)
      inf7=trim(ddir)//'\'//trim(inf7)
      
c  read data
      call input(inf1,text1,n,form,dt,tmin,tsec,r)
      call input(inf2,text2,n,form,dt,tmin,tsec,s)
      call input(inf3,text4,n,form,dt,tmin,tsec,t)
      call input(inf4,text4,n,form,dt,tmin,tsec,x)
      call input(inf5,text5,n,form,dt,tmin,tsec,y)
      call input(inf6,text6,n,form,dt,tmin,tsec,z)
      call input(inf7,text7,n,form,dt,tmin,tsec,w)

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

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

      do i=1,n
        aik(1,1)=aik(1,1)+r(i)*r(i)
        aik(1,2)=aik(1,2)+r(i)*s(i)
        aik(1,3)=aik(1,3)+r(i)*t(i)
        aik(1,4)=aik(1,4)+r(i)*x(i)
        aik(1,5)=aik(1,5)+r(i)*y(i)
        aik(1,6)=aik(1,6)+r(i)*z(i)
        
        aik(2,2)=aik(2,2)+s(i)*s(i)
        aik(2,3)=aik(2,3)+s(i)*t(i)
        aik(2,4)=aik(2,4)+s(i)*x(i)
        aik(2,5)=aik(2,5)+s(i)*y(i)
        aik(2,6)=aik(2,6)+s(i)*z(i)
        
        aik(3,3)=aik(3,3)+t(i)*t(i)
        aik(3,4)=aik(3,4)+t(i)*x(i)
        aik(3,5)=aik(3,5)+t(i)*y(i)
        aik(3,6)=aik(3,6)+t(i)*z(i)
        
        aik(4,4)=aik(4,4)+x(i)*x(i)
        aik(4,5)=aik(4,5)+x(i)*y(i)
        aik(4,6)=aik(4,6)+x(i)*z(i)

        aik(5,5)=aik(5,5)+y(i)*y(i)
        aik(5,6)=aik(5,6)+y(i)*z(i)

        aik(6,6)=aik(6,6)+z(i)*z(i)
        
        rs(1)= rs(1)+r(i)*w(i)
        rs(2)= rs(2)+s(i)*w(i)
        rs(3)= rs(3)+t(i)*w(i)
        rs(4)= rs(4)+x(i)*w(i)
        rs(5)= rs(5)+y(i)*w(i)
        rs(6)= rs(6)+z(i)*w(i)
                 
      enddo
      
      do i=1,5
        do k=i+1,6
          aik(k,i)=aik(i,k)
        enddo
      enddo

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

      wq=0.
      wrq=0.
      do i=1,n
      syn(i)=f(1)*r(i)+f(2)*s(i)+f(3)*t(i)+f(4)*x(i)+f(5)*y(i)+f(6)*z(i)
        wr(i)=w(i)-syn(i)      
        wq = wq+ w(i)**2
        wrq=wrq+wr(i)**2
      enddo

      write(6,1) trim(inf7),trim(inf1),f(1),trim(inf2),f(2),
     & trim(inf3),f(3),trim(inf4),f(4),trim(inf5),f(5),trim(inf6),
     & f(6),wq,wrq,100.*sqrt(wrq/wq)
    1 format(/' Decomposition of signal w = ',a,':'/
     &        '   coefficient of series r = ',a,f13.6/
     &        '   coefficient of series s = ',a,f13.6/
     &        '   coefficient of series t = ',a,f13.6/ 
     &        '   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 residual       ',e10.3/
     & ' rms residual             ',f10.6,' %')

      synf=trim(inf7)//'.syn'
      call output(synf,text1,n,form,dt,tmin,tsec,syn)

      resf=trim(inf7)//'.res'
      call output(resf,text1,n,form,dt,tmin,tsec,wr)
      
c  compose winplot.par file
      open(12,file='winplot.par')      
      write(12,*) '0, 9, 24., 16., 1, 0, 0.8'
      write(12,'(a)') trim(inf7)
      write(12,'(a)') trim(synf)
      write(12,'(a)') trim(resf)      
      write(12,'(a)') trim(inf1)
      write(12,'(a)') trim(inf2)
      write(12,'(a)') trim(inf3)
      write(12,'(a)') trim(inf4)
      write(12,'(a)') trim(inf5)
      write(12,'(a)') trim(inf6)

      write(12,'(a)') 'linreg6 synthesis'
      close(12)
         
      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

