m_a3_interp.f90

module m_a3_interp
  use m_a3_types

  implicit none
  private

  public :: a3_interp0
  public :: a3_interp1
  public :: a3_interp0_to_grid

contains

  function a3_interp0(tree, rr, ivs, n_var) result(vals)
    use m_a3_utils, only: a3_get_loc
    type(a3_t), intent(in) :: tree
    real(dp), intent(in)    :: rr(3)
    integer, intent(in)     :: n_var
    integer, intent(in)     :: ivs(n_var)
    real(dp)                :: vals(n_var)
    type(a3_loc_t)         :: loc

    loc = a3_get_loc(tree, rr)

    if (loc%id == -1) then
       print *, "a3_interp0 error, no box at ", rr
       stop
    end if

    vals = tree%boxes(loc%id)%cc(loc%ix(1), loc%ix(2), loc%ix(3), ivs)
  end function a3_interp0

  function a3_interp1(tree, r, ivs, n_var, id_guess) result(vals)
    use m_a3_utils, only: a3_get_id_at
    type(a3_t), intent(in)       :: tree
    real(dp), intent(in)          :: r(3)
    integer, intent(in)           :: n_var
    integer, intent(in)           :: ivs(n_var)
    integer, intent(in), optional :: id_guess
    real(dp)                      :: vals(n_var)
    integer                       :: i, iv, id, ix(3)
    real(dp)                      :: r_loc(3), dvec(3), ovec(3), w(2, 2, 2)

    id = a3_get_id_at(tree, r, id_guess)

    if (id <= af_no_box) then
       print *, "a3_interp1: point outside domain", r
       stop
    end if

    ! Compute ix such that r lies between cell centers at ix and ix + 1
    ix = nint((r - tree%boxes(id)%r_min) / tree%boxes(id)%dr)
    r_loc = a3_r_cc(tree%boxes(id), ix)
    dvec  = r - r_loc

    ! Normalize dvec to a value [0, 1]
    dvec = dvec / tree%boxes(id)%dr
    ovec = 1 - dvec

    ! Compute weights of linear interpolation
    w(1, 1, 1) = ovec(1) * ovec(2) * ovec(3)
    w(2, 1, 1) = dvec(1) * ovec(2) * ovec(3)
    w(1, 2, 1) = ovec(1) * dvec(2) * ovec(3)
    w(2, 2, 1) = dvec(1) * dvec(2) * ovec(3)
    w(1, 1, 2) = ovec(1) * ovec(2) * dvec(3)
    w(2, 1, 2) = dvec(1) * ovec(2) * dvec(3)
    w(1, 2, 2) = ovec(1) * dvec(2) * dvec(3)
    w(2, 2, 2) = dvec(1) * dvec(2) * dvec(3)

    do i = 1, size(ivs)
       iv = ivs(i)
       vals(i) = sum(w * tree%boxes(id)%cc(ix(1):ix(1)+1, &
            ix(2):ix(2)+1, ix(3):ix(3)+1, iv))
    end do
  end function a3_interp1

  subroutine a3_interp0_to_grid(tree, rr, iv, amount, to_density)
    use m_a3_utils, only: a3_get_loc
    type(a3_t), intent(inout) :: tree
    integer, intent(in)        :: iv
    real(dp), intent(in)       :: rr(3)
    real(dp), intent(in)       :: amount
    logical, intent(in)        :: to_density
    real(dp)                   :: actual_amount
    type(a3_loc_t)            :: loc
    integer                    :: id, ix(3)

    loc = a3_get_loc(tree, rr)

    if (loc%id == -1) then
       print *, "a3_interp0_to_grid error, no box at ", rr
       stop
    end if

    id = loc%id
    ix = loc%ix

    if (to_density) then
       actual_amount = amount / tree%boxes(id)%dr**3
    else
       actual_amount = amount
    end if

    tree%boxes(id)%cc(ix(1), ix(2), ix(3), iv) = &
         tree%boxes(id)%cc(ix(1), ix(2), ix(3), iv) + &
         actual_amount
  end subroutine a3_interp0_to_grid

end module m_a3_interp