module smooth_mod
                                                                       
 implicit none
                                                                       
contains

subroutine goodie( n1, n2, id, nth,data, sum, rect1, rect2, tridiag, scaleup, diff1, diff2, absval, scale, dist1, dist2)
integer i1,n1,i2,n2,id,nth
real	data(n1,n2,nth), sum(n1,n2,nth)
integer rect1, rect2, tridiag, scaleup, diff1, diff2, absval
real  scale, dist1, dist2
!do i1=1,n1
!  do i2=1,n2
!    sum(i1,i2) = 0.
!  end do
!end do 

!  if ( absval > 0 ) then
!    do i2=1,n2
!      do i1=1,n1
!        data(i1,i2) = abs(data(i1,i2))
!      end do
!    end do
!  end if

  if ( rect1 > 1 .and. n1 > 1) then
    do i2=1,n2  
      if ( tridiag .eq. 0 ) then
        call triangle( rect1, 1, n1, data(1,i2,id), data(1,i2,id))
      else
        call leaky(   dist1, 1, n1, data(1,i2,id), data(1,i2,id))
      end if
    end do
  end if

!  if ( rect2 > 1 .and. n2 > 1) then
!   do i1=1,n1  
!      if ( tridiag .eq. 0 ) then
!        call triangle( rect2, n1,n2, data(i1,1,id), data(i1,1,id))
!      else
!        call leaky(   dist2, n1,n2, data(i1,1,id), data(i1,1,id))
!      end if
!    end do
!  end if

!  if ( scaleup > 0 ) then
!    scale = 1./( rect1**2 * rect2**2 )
!    do i2=1,n2
!      do i1=1,n1
!        data(i1,i2,id) = data(i1,i2,id) * scale
!      end do
!    end do
!  end if

!  if ( diff1 > 0 ) then
!    do i1= n1, 2, -1
!      do i2= 1, n2
!        data(i1,i2) = data(i1,i2) - data(i1-1,i2)
!      end do
!    end do 
!    i1 = 1
!    do i2= 1, n2
!      data(i1,i2) = 0.
!    end do
!  end if
!  if ( diff2 > 0 ) then
!    do i2= n2, 2, -1
!      do i1= 1, n1
!        data(i1,i2) = data(i1,i2) - data(i1,i2-1)
!      end do
!    end do 
!    i2 = 1
!    do i1= 1, n1
!      data(i1,i2) = 0.
!    end do
!  end if

end subroutine goodie

! keyword: tridiagonal smoothing on 1-axis or 2-axis
subroutine leaky( distance, m1, n12, uu, vv )
integer i,                  m1, n12
real    distance        ! input:  1. < distance < infinity
real    uu(m1,n12)      ! data in  is the vector (uu( 1, i), i=1,n12)
real    vv(m1,n12)      ! data out is the vector (vv( 1, i), i=1,n12)
real a, b, dc, side
real vecin( n12), vecout( n12)
a  = - (1.-1./distance)
b = 1.+a*a
dc = b+a+a
a = a/dc
b = b/dc
side = a + b
do i= 1,n12  
  vecin(i) = uu(1,i)
end do 
if ( distance<=1..or. n12.eq.1) then
  call copy( n12, vecin, vecout)
else
  call tris( n12, side, a, b, a, side, vecin, vecout)
end if 
do i= 1,n12   
  vv(1,i) = vecout(i)
end do 

end subroutine leaky

! Convolve with triangle
!
subroutine triangle( nr, m1, n12, uu, vv)
! input:        nr      rectangle width (points) (Triangle base twice as wide.)
! input:        uu(m1,i2),i2=1,n12      is a vector of data.
! output:       vv(m1,i2),i2=1,n12      may be on top of uu
integer nr,m1,n12, i,np,nq
real uu( m1, n12),  vv( m1, n12)        
real pp(n12+nr-1), qq(n12+nr+nr-2), tt(n12)
do i=1,n12  
  qq(i) = uu(1,i)
end do 
if ( n12 .eq. 1 ) then
  call copy( n12, qq, tt)
else
  call boxconv( nr, n12, qq, pp)
  np = nr+n12-1
  call boxconv( nr, np , pp, qq)
  nq = nr+np-1
  do i= 1, n12
    tt(i) = qq(i+nr-1)
  end do 
  do i= 1, nr-1                                   
  ! fold back near end
    tt(i) = tt(i) + qq(nr-i)
  end do 
  do i= 1, nr-1                                   ! fold back far end
    tt(n12-i+1) = tt(n12-i+1) + qq(n12+(nr-1)+i)
  end do
end if 
do i=1,n12  
  vv(1,i) = tt(i)
end do 

end subroutine triangle

subroutine copy( n, xx, yy)
integer i, n
real xx(n), yy(n)
do i= 1, n
  yy(i) = xx(i)
end do 

end subroutine copy 

subroutine boxconv( nb, nx, xx, yy)
! inputs:       nx,  xx(i), i=1,nx      the data
!               nb                      the box length
! output:       yy(i),i=1,nx+nb-1       smoothed data
integer nx, ny, nb, i
real xx(nx), yy(1)
real bb(nx+nb)
if ( nb < 1 .or. nb > nx) then
  call erexit('boxconv') 
end if
ny = nx+nb-1
do i= 1, ny
  bb(i) = 0.
end do 
bb(1) = xx(1)
do i= 2, nx
  bb(i) = bb(i-1) + xx(i)         ! make B(Z) = X(Z)/(1-Z)
end do 
do i= nx+1, ny
  bb(i) = bb(i-1)
end do 
do i= 1, nb
  yy(i) = bb(i)
end do 
do i= nb+1, ny
  yy(i) = bb(i) - bb(i-nb)        ! make Y(Z) = B(Z)*(1-Z**nb)
end do 
do i= 1, ny
  yy(i) = yy(i) / nb
end do 

end subroutine boxconv

! tridiagonal simultaneous equations as in FGDP and IEI
!
subroutine tris( n, endl, a, b, c, endr, d, t )
integer i,       n
real endl, a, b, c, endr, d(n), t(n)
real e(n), f(n), deni(n)
if ( n .eq. 1) then
  t(1) = d(1) / b
  return
end if
e(1) = - a / endl
do i= 2, n-1  
  deni(i) = 1. / ( b + c * e(i-1) )
  e(i) = - a * deni(i)
end do 
f(1) = d(1) / endl
do i= 2, n-1
  f(i) = (d(i) - c * f(i-1)) * deni(i)
end do 
t(n) = ( d(n) - c * f(n-1) )  /  ( endr + c * e(n-1) )
do i= n-1, 1, -1
  t(i) = e(i) * t(i+1) + f(i)
end do 

end subroutine tris 

end module smooth_mod
