# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
module map_pts_reg_mod
  use sep3d_struct_mod
  use basic_util_mod
  implicit none
  contains
  logical function get_selected_points(space,field,loc,val)
    type(sep3d) :: space
    real :: field(:,:,:),loc(:,:)
    real, pointer :: val(:),temp(:)
    integer :: n(3)
    real :: o(3),d(3)
    integer :: i1,i2,i3,i,iout
    get_selected_points=.false.
    allocate(temp(size(loc,2)))
    call set_n_o_d_3d(space,n,o,d)
    iout=0
    do i=1,size(loc,2)
      i1=min(n(1),max(1,nint((loc(3,i)-o(1))/d(1)+1.5)))
      i2=min(n(2),max(1,nint((loc(1,i)-o(2))/d(2)+1.5)))
      i3=min(n(3),max(1,nint((loc(2,i)-o(3))/d(3)+1.5)))
      if (i3 >0 .and. i3 <= n(3) .and. i1>0 .and. i1 <=n(1) .and. i2 >0&
        & .and. i2 <=n(2)) then
        iout= iout+1
        temp(iout)=field(i1,i2,i3)
      else
        write(0,*) "      ERROR:  Trouble grabbing point at ",loc(:3,i&
          &)
        write(0,*) "       INDICIES", i1,i2,i3
        write(0,*) "        LIMITS  ",n
        return
      end if
    end do
    if (iout.eq.0) then
      write(0,*) "       ERROR: found zero points "
      return
    end if
    allocate(val(iout))
    val(1:iout)=temp(1:iout)
    deallocate(temp)
    get_selected_points=.true.
  end function
  subroutine interp_ireg_to_reg(struct,pts,vals,field,sub)
    type(sep3d) :: struct
    real :: pts(:,:),vals(:)
    real,pointer :: field(:,:,:),sumf(:,:,:)
    real, pointer :: filter(:,:,:)
    integer :: n(3)
    integer :: bm1,bm2,bm3
    integer :: em1,em2,em3
    integer :: bf1,bf2,bf3
    integer :: ef1,ef2,ef3
    integer :: nf1,nf2,nf3
    integer :: c1,c2,c3
    real :: o(3),d(3),dist2,sub(3)
    integer :: i
    call set_n_o_d_3d(struct,n,o,d)
    allocate(sumf(n(1),n(2),n(3)))
    allocate(field(n(1),n(2),n(3)))
    nf1=floor(n(1)/sub(1))
    nf2=floor(n(2)/sub(2))
    nf3=floor(n(3)/sub(3))
    allocate(filter(nf1*2+1,nf2*2+1,nf3*2+1))
    do c3=1,nf3*2+1
      do c2=1,nf2*2+1
        do c1=1,nf1*2+1
          dist2=(c1-nf1-1)**2+(c2-nf2-1)**2+(c3-nf3-1)**2+.05
          !1/(dist^2+eps)
          filter(c1,c2,c3)=1./dist2
        end do
      end do
    end do
    sumf=0
    field=0
    do i=1,size(pts,2)
      c1=(pts(3,i)-o(1))/d(1)+1.5
      c2=(pts(1,i)-o(2))/d(2)+1.5
      c3=(pts(2,i)-o(3))/d(3)+1.5
      if (c1>0 .and. c1 <= size(field,1) .and. c2>0 .and. c2 <= size&
        &(field,2) .and.c3>0 .and. c3 <= size(field,3)) then
        bm1=max(1,c1-nf1)
        em1=min(n(1),c1+nf1)
        bm2=max(1,c2-nf2)
        em2=min(n(2),c2+nf2)
        bm3=max(1,c3-nf3)
        em3=min(n(3),c3+nf3)
        bf1=nf1+1+(bm1-c1)
        ef1=nf1+1+(em1-c1)
        bf2=nf2+1+(bm2-c2)
        ef2=nf2+1+(em2-c2)
        bf3=nf3+1+(bm3-c3)
        ef3=nf3+1+(em3-c3)
        sumf(bm1:em1,bm2:em2,bm3:em3) = sumf(bm1:em1,bm2:em2,bm3:em3&
          &) +filter(bf1:ef1,bf2:ef2,bf3:ef3)
        field(bm1:em1,bm2:em2,bm3:em3)= field(bm1:em1,bm2:em2,bm3:em3&
          &)+filter(bf1:ef1,bf2:ef2,bf3:ef3)*vals(i)
      end if
    end do
    where(sumf>0.)
      field=field/sumf
end where
    deallocate(filter)
  end subroutine
  subroutine map_ireg_to_reg(struct,pts,vals,field)
    type(sep3d) :: struct
    real :: pts(:,:),vals(:)
    real,pointer :: field(:,:,:)
    integer :: n(3),i1,i2,i3,i
    real :: o(3),d(3)
    call set_n_o_d_3d(struct,n,o,d)
    allocate(field(n(1),n(2),n(3)))
    field=0.
    do i=1,size(pts,2)
      i1=(pts(3,i)-o(1))/d(1)+1.5
      i2=(pts(1,i)-o(2))/d(2)+1.5
      i3=(pts(2,i)-o(3))/d(3)+1.5
      if (i1>0 .and. i1 <= size(field,1) .and. i2>0 .and. i2 <= size&
        &(field,2) .and.i3>0 .and. i3 <= size(field,3)) then
        field(i1,i2,i3)=vals(i)
      end if
    end do
  end subroutine
end module
