! smooth by convolving with triangle in two dimensions.
!
module triang
implicit none
contains
!====================================================
subroutine triangle3( rect1, rect2, n1, n2, uu, vv)
integer       :: i1,i2, rect1, rect2, n1, n2,minz,maxz
real          :: uu(n1,n2), vv(n1,n2),c2(rect2),c1(rect1),temp1(n1),temp2(n2)
real          :: ss(n1,n2)

do i2=1,rect2/2
   c2(i2)=(2.*(i2-1)/float(rect2-1))/int(rect2/2)
end do

do i2=rect2/2+1,rect2
   c2(i2)=(1-2.*(i2-rect2/2-1)/float(rect2-1))/int(rect2/2)
end do

do i1=1,rect1/2
   c1(i1)=(2.*(i1-1)/float(rect1-1))/int(rect1/2)
end do

do i1=rect1/2+1,rect1
   c1(i1)=(1-2.*(i1-rect1/2-1)/float(rect1-1))/int(rect1/2)
end do

do i1=1,n1
   do i2=1,n2
      minz=max(i2-rect2/2,1)
      maxz=min(i2+rect2/2,n2)
      temp2(i2)=sum(uu(i1,minz:maxz)*c2(max(1,rect2/2-i2+2):min(rect2,n2-(i2-rect2/2)+1)))
   end do
   uu(i1,:)=temp2
end do

do i2=1,n2
   do i1=1,n1
      minz=max(i1-rect1/2,1)
      maxz=min(i1+rect1/2,n1)
      temp1(i1)=sum(uu(minz:maxz,i2)*c1(max(1,rect1/2-i1+2):min(rect1,n1-(i1-rect1/2)+1)))
   end do
   vv(:,i2)=temp1
end do

return
end subroutine
!======================================================

!====================================================
subroutine triangle2( rect1, rect2, n1, n2, uu, vv)
integer       :: i1,i2, rect1, rect2, n1, n2
real          :: uu(n1,n2), vv(n1,n2)
real          :: ss(n1,n2)

write(0,*) "from triangle2",rect1,rect2,n1,n2

do i1= 1, n1
  call triangle( rect2, n1, n2, uu(i1,1), ss(i1,1))
  do i2= 1, n2
    call triangle( rect1, 1, n1, ss(1,i2), vv(1,i2))
  end do
end do
return
end subroutine
!======================================================
! 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)

qq(1:n12) = uu(1,1:n12)
if( n12 == 1 ) then
  qq=tt
else 
  call boxconv( nr, n12, qq, pp)
  np = nr+n12-1
  call boxconv( nr, np , pp, qq)
  nq = nr+np-1
end if
tt = qq(nr:n12+nr-1)
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
vv(1,1:n12) = tt
end subroutine
!========================================
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) return
ny = nx+nb-1
bb(1:ny) = 0.
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

yy(1:nb) = bb(1:nb)

do i= nb+1, ny
  yy(i) = bb(i) - bb(i-nb) ! make Y(Z) = B(Z)*(1-Z**nb)
end do

yy(1:ny) = yy(1:ny) / nb
return
end subroutine
!========================================
end module
