!! weighted smoothing
program WSmo_
  use sep
  use utilities
  implicit none

  type(paxis) :: a1,a2
  integer     :: i1,i2
  integer     :: j1,j2
  integer     :: k1,k2
  real, dimension(:,:), pointer :: din,dou  !! data
  real, dimension(:,:), pointer :: wgh      !! weight
  logical :: verb   !! verbose
  integer :: rect1, rect2, nrepeat,irepeat

  call sep_init()
  call from_param("verb",verb,.false.)
  call from_param("rect1",rect1,0)
  call from_param("rect2",rect2,0)
  call from_param("repeat",nrepeat,1)

  call hinaxis(a1,1); a1%l=" "; call houaxis(a1,1)
  call hinaxis(a2,2); a2%l=" "; call houaxis(a2,2)
  call to_history("n3",1)

  if(verb) write(0,*) "allocate"
  allocate(din(a1%n,a2%n))
  allocate(dou(a1%n,a2%n))
  allocate(wgh(a1%n,a2%n))

  if(verb) write(0,*) "read din"
  call sep_read(din)

  if(verb) write(0,*) "read wgh"
  call sep_read(wgh)
  wgh = wgh **2

  if(verb) write(0,*) "compute"

  do irepeat=1,nrepeat
     if(verb) write(0,*) "irepeat=",irepeat
     dou=0.
     do    i2=1+rect2,a2%n-rect2
!        if(verb) write(0,*) irepeat,i2,a2%n
        do i1=1+rect1,a1%n-rect1

           j1=i1-rect1
           k1=i1+rect1

           j2=i2-rect2
           k2=i2+rect2

           dou(i1,i2) = sum( wgh(j1:k1,j2:k2)*din(j1:k1,j2:k2) ) &
           /            sum( wgh(j1:k1,j2:k2)                  )

        end do
     end do
     if(verb) write(0,*) "update boundaries"
     do i1=1,rect1
        dou(i1,:) = dou(rect1+1,:)
     end do
     do i1=a1%n-rect1+1,a1%n
        dou(i1,:) = dou(a1%n-rect1,:) 
     end do

     do i2=1,rect2
        dou(:,i2) = dou(:,rect2+1)
     end do
     do i2=a2%n-rect2+1,a2%n
        dou(:,i2) = dou(:,a2%n-rect2) 
     end do

     din = dou
  end do


  if(verb) write(0,*) "write dou"
  call sep_write(dou)

  call exit(0)
end program WSmo_








