# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
module punish_change_mod
  use adj_mod
  use sep
  use sep3d_struct_mod
  use tomo_space_conv_mod
  use gsmooth_mod
  implicit none
  real, private, save :: myd,mye
  real, private, pointer :: scale(:)
  type(sep3d),private,save :: vel_s
  contains
  logical function init_punish_change(sc)
    real :: sc(:)
    integer :: n(3),i3,nw(3)
    real :: o(3),d(3)
    real, pointer :: junk(:,:,:)
    init_punish_change=.false.
    call set_n_o_d_3d(vel_s,n,o,d)
    allocate(scale(size(sc)))
    allocate(junk(n(1),n(2),n(3)))
    scale=mye
    nw=5
    do i3=1,3
      nw(i3)=min(5,n(i3))
    end do
    where(sc.ne.0)
      scale=0.
end where
    junk=reshape(scale,n)
    call gauss_smooth(junk,nw(1),nw(2),nw(3),1.)
    scale=reshape(junk,(/size(scale)/))
    deallocate(junk)
!call srite("mtest.H",scale,size(scale)*4)
    init_punish_change=.true.
  end function
  integer function punish_change_op(adj,add,model,data)
    logical, intent(in) :: adj,add
    real :: model(:),data(:)
    call adjnull(adj,add,model,data)
    if (.not. adj) then
      data= data+model*scale
    else
      model= model+data*scale
    end if
    punish_change_op=0
  end function
  real function return_change_sum()
    return_change_sum=myd
  end function
  logical function check_punish_change_params()
    if (.not. vel_tomo_structure(vel_s)) then
      write(0,*) "     ERRROR: Grabbing vel_tomo_structure "
      return
    end if
    check_punish_change_params=.false.
    call from_param("punish_eps",mye,.001)
! mye=mye/product(vel_s%n(1:sep3d_ndims(vel_s)))
    call from_param("punish_sum",myd,0.)
    check_punish_change_params=.true.
  end function
end module
