c   .    1    .    2    .    3    .    4    .    5    .    6    .    7..
c
c     linear combination of two signals, to fit a third one
c     July 1996 e.w.
c

      program linreg2
      use nas_system

      parameter(ndim=60000)
      dimension x(ndim),y(ndim),z(ndim),zz(ndim)
      character inf1*24,inf2*24,inf3*24,outf*28,form*20
      character*72 text1,text2,text3

      write(6,*)
      write(6,'(a,a)') 'linreg2 file1 file2 file3 [n1 n2]:', 
     & '   solve z = a*x + b*y + residue'
      write(6,'(a,a)') 'where x=signal1, y=signal2, z=signal3, and the',
     & ' constants a and b are '
      write(6,'(a,a)') 'determined for best fit. A window may be speci',
     & 'fied with n1 and n2.'
      write(6,'(a)') 'The residue is written into file3.res'
      narg=iargc()
      if(narg.lt.3) stop

      call getarg(1,inf1)
      call getarg(2,inf2)
      call getarg(3,inf3)
      
      if (narg.eq.5) then
        call getarg(4,form)
        read(form,*) nstart
        call getarg(5,form)
        read(form,*) nstop
      else
        nstart=1
        nstop=60000
      endif

c  read data
      n1=ndim
      call input(inf1,text1,n1,form,dt,tmin,tsec,x)
      n2=ndim
      call input(inf2,text2,n2,form,dt,tmin,tsec,y)
      n3=ndim
      call input(inf3,text3,n3,form,dt,tmin,tsec,z)
      n=min(n1,n2,n3)
      nstart=max(nstart,1)
      nstop=min(nstop,n)

      call mean(x,n)
      call mean(y,n)
      call mean(z,n)

      a1=0.
      a2=0.
      a3=0.
      b1=0.
      b2=0.

      do i=nstart,nstop
        a1=a1+x(i)**2
        a2=a2+y(i)**2
        a3=a3+x(i)*y(i)
        b1=b1+x(i)*z(i)
        b2=b2+y(i)*z(i)
      enddo
      det=a1*a2-a3*a3
      a=(b1*a2-b2*a3)/det
      b=(b2*a1-b1*a3)/det

      do i=1,n
        zz(i)=z(i)-a*x(i)-b*y(i)
      enddo

      call mean(zz,n)

      xq=0.
      yq=0.
      zq=0.
      zzq=0.
      do i=nstart,nstop
        xq=xq+x(i)**2
        yq=yq+y(i)**2
        zq=zq+z(i)**2
        zzq=zzq+zz(i)**2
      enddo
      nf=nstop-nstart+1
      xq=sqrt(xq/nf)
      yq=sqrt(yq/nf)
      zq=sqrt(zq/nf)
      zzq=sqrt(zzq/nf)

      write(6,1) a,b,xq,yq,zq,zzq,100.*zzq/zq,nstart,nstop
    1 format(' Koeff.   series x = ',f15.6,'  series y = ',f15.6/
     &       ' rms Amp. series x = ',f15.6,'  series y = ',f15.6/
     &       '          series z = ',f15.6,'  residue = ',f15.6/
     &       ' the residue is ',f6.1,' % of the rms amplitude'/
     &       ' data were fitted from samples ',i7,' to ',i7,
     &       ' (mean removed)')

c  Restsignal in File inf3.res schreiben

      text3='Rest'
      outf=trim(inf3)//'.res'
      call output(outf,text3,n,form,dt,tmin,tsec,zz)

      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 form*20,name*24,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('reading file ',a,', ',i7,' /',i7,' samples')
      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(' end of file encountered after ',i7,' samples')
      close(7)
      return
   25 write(6,26) trim(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*28,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)') ' illegal format specification ',form
        stop
      endif
      write(6,5) trim(name),n,form
    5 format('writing file ',a,i10,' samples in format ',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