module C2R2D

  implicit none

  integer, private :: nsx,nsz,norm,cnx,cnz,rnx,rnz,cnh
  real   , private :: xmin,xmax,zmin,zmax,dxx,dzz

contains

  !----------------------------------------------------------------  

  subroutine C2R2D_init(insx,insz,inorm,idxx,idzz,icnx,icnz,irnx,irnz,ixmin,izmin,icnh)
    integer :: insx,insz,inorm,icnx,icnz,irnx,irnz,icnh
    real    :: idxx,idzz,ixmin,izmin
    nsx = insx;  nsz = insz;   
    xmin=ixmin;  zmin= izmin
    cnx = icnx;  cnz = icnz
    rnx = irnx;  rnz = irnz
    dxx = idxx;  dzz = idzz
    norm= inorm; cnh = icnh

    !! Define max offsets
    xmax=dble(cnx-1.)*dxx+xmin
    zmax=dble(cnz-1.)*dzz+zmin

    write(0,*)                '        ************************************************'
    write(0,*)                '        * INTERP PARAMETERS'
    write(0,'(a,i6,i6)')         '         * nsx,nsz:',nsx,nsz
    write(0,'(a,i6,i6)')         '         * cnx,cnz:',cnx,cnz
    write(0,'(a,i6,i6)')         '         * rnx,rnz:',rnx,rnz
    write(0,'(a,2f11.2)')     '         * dxx,dzz:',dxx,dzz
    write(0,'(a,2f11.2)')     '         * xmn,zmn:',xmin,zmin
    write(0,*)                '        ************************************************'
    write(0,*) ' '
  end subroutine C2R2D_init

  !---------------------------------------------------------------- 

  subroutine C2R2D_IMAGE(RImg,rays,CImg,gnorm)
    integer:: ix1,ix3,ie1,ie3,ix,jx,iz,jz,ih
    real     :: CImg(:,:,:),RImg(:,:,:),rays(:,:,:),gnorm(:,:)
    real     :: x,z,ddx,ddz,r,s,dr,pt

    dr = sqrt( dxx**2 + dzz**2 )

    write(0,*)                          '        ************************************************'
    write(0,*)                          '        * FIELD PARAMETERS'
    write(0,'(a,2f11.2,2f11.2,2f11.2)') '         * Fields min/max:',minval(gnorm),maxval(gnorm)
    write(0,'(a,2f11.2,2f11.2,2f11.2)') '         * XLoc min/max:',minval(rays(:,1,1)),maxval(rays(:,1,1))
    write(0,*)                          '        ************************************************'
    write(0,*) ' '

    !$OMP PARALLEL DO SCHEDULE(DYNAMIC) PRIVATE(ih,ie3,ie1,x,z,ix,iz,pt,ddx,ddz,r,s)
    EHLoop: do ih=1,cnh
       E3Loop: do ie3=1,rnz
          E1Loop: do ie1=1,rnx

             x=rays(ie1,ie3,1)
             ix=floor((x-xmin)/dxx)+1
             if (ix < 1+nsx .or. ix > cnx-nsx ) cycle

             z=rays(ie1,ie3,2)
             iz=floor((z-zmin)/dzz)+1
             if (iz < 1+nsz .or. iz > cnz-nsz ) cycle

             X1Loop: do jx=ix-nsx,ix+nsx
                ddx=xmin+(jx-1.)*dxx-x

                X3Loop: do jz=iz-nsz,iz+nsz
                   ddz=zmin+(jz-1.)*dzz-z

                   r = sqrt( ddz**2 + ddx**2 ) / dr
                   s = min( sin(r)/r , 1.0 )
                   CImg(jz,jx,ih)=CImg(jz,jx,ih)+&
                               s*RImg(ie1,ie3,ih)*gnorm(ie1,ie3)

                end do X3Loop
             end do X1Loop

          end do E1Loop
       end do E3Loop
    end do EHLoop
    !$OMP END PARALLEL DO

  end subroutine C2R2D_IMAGE

!! -----------------------------------------------------------------

  subroutine C2R2D_VELOCITY(RC,rays,CC,gnorm)
    integer:: ix1,ix3,ie1,ie3,ix,jx,iz,jz,ih
    real   :: CC(:,:),RC(:,:),rays(:,:,:),gnorm(:,:)
    real   :: x,z,fx,fz
    real   :: w00,w01,w10,w11,sumw

    write(0,*)                          '        ************************************************'
    write(0,*)                          '        * FIELD PARAMETERS'
    write(0,'(a,2f11.2,2f11.2,2f11.2)') '         * Fields min/max:',minval(gnorm),maxval(gnorm)
    write(0,*)                          '        ************************************************'
    write(0,*) ' '

    !! . . Interpolating Cartesian velocity models -> RWE
    !$OMP PARALLEL DO PRIVATE(ie3,ie1,x,z,ix,iz,jx,jz,fx,fz,w00,w01,w10,w11,sumw) 
    do ie3=1,rnz
       do ie1=1,rnx
          x=rays(ie1,ie3,1)
          z=rays(ie1,ie3,2)
          x = min(x,xmax)
          x = max(x,xmin)
          z = min(z,zmax)
          z = max(z,zmin)

          ix=floor((x-xmin)/dxx)+1
          iz=floor((z-zmin)/dzz)+1

          ix = min(ix  ,cnx-1)
          ix = max(ix  , 1  )
          jx = min(ix+1,cnx)

          iz = min(iz  ,cnz-1)
          iz = max(iz  , 1  )
          jz = min(iz+1,cnz)

          fz= z - zmin-(iz-1)*dzz
          fx= x - xmin-(ix-1)*dxx

          fz = max(0., fz)
          fz = min(1., fz)

          fx = max(0., fx)
          fx = min(1., fx)

          !! . . Calculate weights
          w00 = (1.-fz)*(1.-fx)
          w01 = (1.-fz)*(   fx)
          w10 = (   fz)*(1.-fx)
          w11 = (   fz)*(   fx)
          sumw = w00+w01+w10+w11

          RC(ie1,ie3) = RC(ie1,ie3) &
          +( w00 * CC(iz,ix) + w10 * CC(jz,ix) &
          +  w01 * CC(iz,jx) + w11 * CC(jz,jx) )/sumw

       end do
    end do
    !$OMP END PARALLEL DO

  end subroutine C2R2D_VELOCITY

end module C2R2D
