program RCslow_
  use sep
  use util_

  implicit none
  character(len=128) :: name

  real,    pointer :: velCC(:,:), velRC(:,:)
  complex, pointer :: rays(:,:)

  type(myaxis) :: ax,az,at,ag
  integer      :: ix,iz,it,ig, jx,jz
  real         :: x,z

  real         :: xmin,xmax,zmin,zmax
  real         :: fx,fz
  real         :: w00,w01,w10,w11, sumw

  call sep_init()

  name="rays"
  call ainaxis(at,1,name)
  call ainaxis(ag,2,name)

  call hinaxis(az,1)
  call hinaxis(ax,2)

  call houaxis(at,1)
  call houaxis(ag,2)

  call sep_close()

  allocate(velCC(az%n,ax%n)); velCC=0.
  allocate(velRC(at%n,ag%n)); velRC=0.
  allocate(rays(at%n,ag%n))

  call sreed("rays",rays,size( rays)*8)
  call sreed("in" ,velCC,size(velCC)*4)

  xmin = ax%o
  xmax = ax%o + (ax%n-1)*ax%d
  zmin = az%o
  zmax = az%o + (az%n-1)*az%d

  do ig=1,ag%n
     do it=1,at%n

        z=aimag(rays(it,ig))
        x= real(rays(it,ig))

        z = min(z,zmax)
        z = max(z,zmin)
        x = min(x,xmax)
        x = max(x,xmin)

        iz=floor((z-az%o)/az%d)+1
        ix=floor((x-ax%o)/ax%d)+1

        iz = min(iz  ,az%n)
        iz = max(iz  , 1  )
        jz = min(iz+1,az%n)

        ix = min(ix  ,ax%n)
        ix = max(ix  , 1  )
        jx = min(ix+1,ax%n)

        fz= (z-az%o-(iz-1)*az%d)
        fx= (x-ax%o-(ix-1)*ax%d)
        fz = max(0., fz)
        fz = min(1., fz)
        fx = max(0., fx)
        fx = min(1., fx)

        !! 00   01       --> x
        !!    .         |
        !! 10   11      v
        !!               z
        w00=(1.-fx)*(1.-fz)
        w10=(   fx)*(1.-fz)
        w01=(1.-fx)*(   fz)
        w11=(   fx)*(   fz)

        !! normalize the sum of the weights
        sumw = w00+w10+w01+w11
        w00 = w00 / sumw
        w10 = w10 / sumw
        w01 = w01 / sumw
        w11 = w11 / sumw

        velRC(it,ig) = velRC(it,ig) + w00 * velCC(iz,ix)
        velRC(it,ig) = velRC(it,ig) + w10 * velCC(jz,ix)
        velRC(it,ig) = velRC(it,ig) + w01 * velCC(iz,jx)
        velRC(it,ig) = velRC(it,ig) + w11 * velCC(jz,jx)

     end do
  end do

  call srite("out",velRC,size(velRC)*4)

  call exit(0)
end program RCslow_






