# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
module smooth_diff_type_mod
    implicit none
  type sm_type
    integer,pointer :: iclose(:)
    real,pointer :: wt(:)
  end type
end module
module smooth_diff_mod
  use smooth_diff_type_mod
  use adj_mod
  use velocity_mod
  use chain_mod
  use error_mod
  implicit none
  type(sm_type), pointer,private,save :: vals(:)
  real, pointer ,private,save:: var1(:)
  integer, pointer,private,save :: imap(:)
  real, pointer,private,save :: scale(:),temp(:)
  contains
  logical function smooth_diff_init()
    type(sep3d) :: err_s,vel_s,semb
    integer :: nmax
    real, allocatable :: dist(:),var(:,:,:),x(:),y(:),z(:),loc(:,:)
    integer, allocatable :: iloc(:)
    real :: big,pos(3),maxd,myd
    integer :: i1,i2,ipos(1),nots,ix,iy,iz,npts,nclose
    integer :: ns(3),ii
    real :: os(3),ds(3)
    smooth_diff_init=.false.
    call from_param("nclose",nclose,5)
    if (.not. vel_structure(vel_s)) then
      write(0,*) "trouble grabbing velocity structure"
      return
    end if
    if (.not. get_error_map(imap)) then
      write(0,*) "trouble getting error map"
      return
    end if
    allocate(scale(maxval(imap)))
    allocate(temp(maxval(imap)))
    scale=0
    do i1=1,size(imap)
      scale(imap(i1))= scale(imap(i1))+1
    end do
    big=10.*(vel_s%o(2)+vel_s%d(2)*vel_s%n(2)-1)
    call from_param("maxd",maxd,vel_s%d(2)*20.)
    call init_sep3d("err",err_s,"INPUT")
    call sep3d_grab_headers("err",err_s,npts)
    allocate(x(npts),y(npts),z(npts))
    allocate(loc(3,npts))
    allocate(vals(npts))
    call init_sep3d("semb_coher",semb,"INPUT")
    call set_n_o_d_3d(semb,ns,os,ds)
    call sep3d_grab_headers("err",err_s,npts)
    call sep3d_grab_key_vals(err_s,"s_x",x)
    call sep3d_grab_key_vals(err_s,"s_y",y)
    call sep3d_grab_key_vals(err_s,"s_z",z)
    call sep3df_delete(err_s)
    loc(1,:)=x
    loc(2,:)=y
    loc(3,:)=z
    deallocate(x,y,z)
    allocate(dist(nclose),iloc(nclose),var1(npts))
    allocate(var(ns(1),ns(2),ns(3)))
    maxd=maxd*maxd
    do i2=1,size(loc,2)
      pos=loc(:,i2)
      dist=big
      do i1=1,size(loc,2)
        myd=(pos(1)-loc(1,i1))**2+ (pos(2)-loc(2,i1))**2+(pos(3)-loc(3&
          &,i1))**2
        if (myd< maxd .and. i1.ne.i2) then
          ipos=maxloc(dist)
          if (dist(ipos(1)) > myd) then
            iloc(ipos(1))=i1
            dist(ipos(1))=myd
          end if
        end if
      end do
      ii=0
      do i1=1,size(dist)
        if (dist(i1) < maxd) then
          ii= ii+1
        end if
      end do
      allocate(vals(i2)%iclose(ii))
      allocate(vals(i2)%wt(ii))
      ii=0
      do i1=1,size(dist)
        if (dist(i1) < maxd) then
          ii= ii+1
          vals(i2)%iclose(ii)=iloc(i1)
          vals(i2)%wt(ii)=1./dist(i1)
! write(0,*) "DYING 1",i2,i1,vals(i2)%iclose(ii),vals(i2)%wt(ii)
        end if
      end do
      iz=max(1,min(ns(3),nint((loc(3,i2)-os(1))/ds(1)+1.)))
      ix=max(1,min(ns(1),nint((loc(1,i2)-os(2))/ds(2)+1.)))
      iy=max(1,min(ns(2),nint((loc(2,i2)-os(3))/ds(3)+1.)))
      var1(i2)=var(iz,ix,iy)
      vals(i2)%wt=vals(i2)%wt/sum(vals(i2)%wt)
    end do
    scale=1./scale
    deallocate(dist,iloc,loc,var,scale)
    var1=var1/(sum(var1)/size(var1))
    smooth_diff_init=.true.
  end function
  integer function compress_gamma(adj, add, model,data)
    logical, intent(in) :: adj,add
    real :: model(:),data(:)
    integer :: i
    call adjnull(adj,add,model,data)
    if (.not. adj) then
      do i=1,size(model)
        data(imap(i))= data(imap(i))+model(i)*scale(imap(i))
      end do
    else
      do i=1,size(model)
        model(i)= model(i)+data(imap(i))*scale(imap(i))
      end do
    end if
    compress_gamma=0
  end function
  integer function close_op(adj,add,model,data)
    logical, intent(in) :: adj,add
    real :: model(:),data(:)
    call chain(smooth_op,compress_gamma,adj,add,model,data,temp)
    close_op=0
  end function
  integer function smooth_op(adj,add,model,data)
    logical, intent(in) :: adj,add
    real :: model(:),data(:)
    integer :: i1,i2,idat
    call adjnull(adj,add,model,data)
    if (.not. adj) then
      do i2=1,size(model)
        if (associated(vals(2)%iclose)) then
          data(idat)= data(idat)+model(i2)*var1(i2)
          do i1=1,size(vals(i2)%iclose)
            data(i2)= data(i2)-model(vals(i2)%iclose(i1))*vals&
              &(i2)%wt(i1)*var1(i2)
          end do
        end if
      end do
    else
      do i2=1,size(model)
        if (associated(vals(i2)%iclose)) then
          model(i2)= model(i2)+data(i2)*var1(i2)
          do i1=1,size(vals(i2)%iclose)
            model(vals(i2)%iclose(i1))= model(vals(i2)%iclose(i1&
              &))-data(i2)*vals(i2)%wt(i1)*var1(i2)
          end do
        end if
      end do
    end if
    smooth_op=0
  end function
  integer function ndomain_smooth()
    ndomain_smooth=maxval(imap)
  end function
end module
