module codmap_
  use util_
  implicit none

  integer, pointer, private :: lmap(:)
  real,    pointer, private :: lwgh(:)
  type(myaxis),     private :: la1,la2,la3
  real,             private :: lrho

contains
  subroutine codmap_init(a1,a2,a3,rho)
    type(myaxis), intent(in) :: a1,a2,a3
    real,         intent(in) :: rho
    integer :: ia1
    real    :: zr
    integer :: zn
    real    :: zw

    la1=a1
    la2=a2
    la3=a3
    lrho=rho

    allocate(lmap(la1%n),lwgh(la1%n))

    !! loop over stretched axis
    do ia1=1,la1%n
       !! the stretch is relative to zero
       zr=(la1%o+(ia1-1)*la1%d)*lrho

       !! find index relative to the origin of the fast axis 
       zn = 1 +floor((zr- la1%o)               / la1%d)
       zw = 1.-     ((zr-(la1%o+(zn-1)*la1%d)) / la1%d)

       !! I am outside the grid
       if(zn<1) then
          zn=1
          zw=-1.
       end if

       if(zn>la1%n-1) then
          zn=la1%n-1
          zw=-1.
       end if

       lmap(ia1) = zn
       lwgh(ia1) = zw
    end do
  end subroutine codmap_init

  subroutine codmap_execute(din,dou)
    real    :: din(:,:,:),dou(:,:,:)
    integer :: ia1,ia2,ia3

    a3:do ia3=1,la3%n
       a2:do ia2=1,la2%n
          a1:do ia1=1,la1%n
             if(lwgh(ia1)>0) then
                dou(     ia1,   ia2,ia3) = &
                din(lmap(ia1)  ,ia2,ia3) *    lwgh(ia1)  + & 
                din(lmap(ia1)+1,ia2,ia3) * (1-lwgh(ia1))
             else
                dou( ia1,ia2,ia3) =0.
             end if

          end do a1
       end do a2
    end do a3

  end subroutine codmap_execute
end module codmap_




