 MODULE adaptsum_mod
 
 IMPLICIT NONE
 PUBLIC  ::  adapts, sto_sl

 contains 

!########################################################################

  subroutine adapts(n1,n2,data1,data2)
!************************************************ 
! Adaptive subtraction				*
!						*
! data1 = reference data 			*
! data2 = data before and after adaption	*
!************************************************

  implicit none

  integer					:: i1,i2,n1,n2
  real,dimension(:,:),intent(in)		:: data1
  real,dimension(:,:),intent(inout)		:: data2
  real,dimension(n1*n2)				:: d1,d2,d3,d4,AA,BB
  real,dimension(2*n1*n2-1)			:: CC
  real						:: a,b

!!
!! Generates vectors from the input seismograms
!!
  
  do i2=1,n2
	do i1=1,n1
		d1((i2-1)*n1+i1)=data1(i1,i2)
		d2((i2-1)*n1+i1)=data2(i1,i2)
	end do
  end do

!!
!! Performs correlations
!!
  do i1=1,n1*n2
	AA(i1)=dot_product(d2(1:n1*n2-(i1-1)),d2(i1:n1*n2))
	BB(i1)=dot_product(d1(1:n1*n2-(i1-1)),d2(i1:n1*n2))
  end do

!!
!! Solves for the filter
!! 
  CC=AA
  do i1=2,n1*n2
	CC(i1+n1*n2)=AA(i1)
  end do  
  call sto_sl ( n1*n2, CC, BB, BB, 0 )
!  write(0,*) "Filter ->", BB(:)
!!
!! Applies the filter
!!
  d3=d2
  do i1=1,n1*n2
	a=d2(i1)
	d2(n1*n2-(i1-1))=a
  end do
  do i1=1,n1*n2
	b=dot_product(BB(1:i1),d2(n1*n2-(i1-1):n1*n2))
	d4(i1)=d3(i1)-b
  end do

!!
!! Back to the original shape
!!
  data2=reshape(d4,(/n1, n2/))

  return
  end  subroutine adapts

!##################################################################

subroutine sto_sl ( n, a, b, x, job )

!***********************************************************************
!
!! STO_SL solves the real Toeplitz system A * X = B.
!
!  Modified:
!
!    11 March 2001
!
!  Author:
!
!    Oleg Arushanian, MK Samarin, Valentin Voevodin, Evgeny Tyrtyshnikov, 
!    Burton Garbow, James Boyle, Wayne Cowell, Kenneth Dritz
!
!    FORTRAN90 translation by John Burkardt
!
!  Reference:
!
!    Oleg Arushanian, MK Samarin, Valentin Voevodin, Evgeny Tyrtyshnikov, 
!    Burton Garbow, James Boyle, Wayne Cowell, Kenneth Dritz,
!    The TOEPLITZ Package Users Guide,
!    Argonne National Laboratory,
!    ANL-83-16, 1983.
!
!  Parameters:
!
!    Input, integer N, the order of the matrix.
!
!    Input, real A(2*N-1), the first row of the Toeplitz matrix, followed by
!    the first column of the Toeplitz matrix beginning with the second element.
!
!    Input, real B(N) the right hand side vector.
!
!    Output, real X(N), the solution vector.  X and B may share the
!    same storage.
!
!    Input, integer JOB,
!    0 to solve A*X=B,
!    nonzero to solve Transpose(A)*X=B.
!
  implicit none

  integer n
  real a(2*n-1)
  real b(n)
  real c1(n-1)
  real c2(n-1)
  integer i
  integer job
  integer nsub
  real r1
  real r2
  real r3
  real r5
  real r6
  real x(n)

  if ( n < 1 ) then
    return
  end if
!
!  Solve the system with the principal minor of order 1.
!
  r1 = a(1)
  x(1) = b(1) / r1

  if ( n == 1 ) then
    return
  end if
!
!  Recurrent process for solving the system with the Toeplitz matrix.
!
  do nsub = 2, n

!
!  Compute multiples of the first and last columns of the inverse of
!  the principal minor of order NSUB.
!
    if ( job == 0 ) then
      r5 = a(n+nsub-1)
      r6 = a(nsub)
    else
      r5 = a(nsub)
      r6 = a(n+nsub-1)
    end if

    if ( 2 < nsub ) then

      c1(nsub-1) = r2

      do i = 1, nsub-2
        if ( job == 0 ) then
          r5 = r5 + a(n+i) * c1(nsub-i)
          r6 = r6 + a(i+1) * c2(i)
        else
          r5 = r5 + a(i+1) * c1(nsub-i)
          r6 = r6 + a(n+i) * c2(i)
        end if
      end do

    end if

    r2 = - r5 / r1
    r3 = - r6 / r1
    r1 = r1 + r5 * r3

    if ( 2 < nsub ) then

      r6 = c2(1)
      c2(nsub-1) = 0.0E+00

      do i = 2, nsub-1
        r5 = c2(i)
        c2(i) = c1(i) * r3 + r6
        c1(i) = c1(i) + r6 * r2
        r6 = r5
      end do

    end if

    c2(1) = r3
!
!  Compute the solution of the system with the principal minor of order NSUB.
!
    r5 = 0.0E+00
    do i = 1, nsub-1
      if ( job == 0 ) then
        r5 = r5 + a(n+i) * x(nsub-i)
      else
        r5 = r5 + a(i+1) * x(nsub-i)
      end if
    end do

    r6 = ( b(nsub) - r5 ) / r1

    x(1:nsub-1) = x(1:nsub-1) + c2(1:nsub-1) * r6
    x(nsub) = r6

  end do

  return
end subroutine sto_sl

end module adaptsum_mod
