module C2R3D

  implicit none

  integer, private :: nsx,nsy,nsz,norm,cnx,cny,cnz,rnx,rny,rnz
  real   , private :: xmin,xmax,ymin,ymax,zmin,zmax,dxx,dyy,dzz

contains

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

  subroutine C2R3D_init(insx,insy,insz,inorm,idxx,idyy,idzz,&
         icnx,icny,icnz,irnx,irny,irnz,ixmin,iymin,izmin)
    integer :: insx,insy,insz,inorm,icnx,icny,icnz,irnx,irny,irnz
    real    :: idxx,idyy,idzz,ixmin,iymin,izmin
    nsx = insx;  nsy=insy;  nsz = insz;   
    norm= inorm
    xmin=ixmin;  ymin = iymin; zmin= izmin
    cnx = icnx;  cny=icny;  cnz = icnz
    rnx = irnx;  rny = irny; rnz = irnz
    dxx = idxx;  dyy=idyy; dzz = idzz

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

    write(0,*)                          '        ************************************************'
    write(0,*)                          '        * INTERP PARAMETERS'
    write(0,'(a,i4,i4,i4)')             '         * nsx,nsy,nsz:',nsx,nsy,nsz
    write(0,'(a,i4,i4,i4)')             '         * cnx,cny,cnz:',cnx,cny,cnz
    write(0,'(a,i4,i4,i4)')             '         * rnx,rny,rnz:',rnx,rny,rnz
    write(0,'(a,2f11.2,2f11.2,2f11.2)') '         * dxx,dyy,dzz:',dxx,dyy,dzz
    write(0,'(a,2f11.2,2f11.2,2f11.2)') '         * xmn,ymn,zmn:',xmin,ymin,zmin
    write(0,*)                          '        ************************************************'
    write(0,*) ' '
  end subroutine C2R3D_init

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

  subroutine C2R3D_CART(RC,rays,CC,gnorm,adj)
    integer:: ix1,ix2,ix3,ie1,ie2,ie3,ix,jx,iz,jz,jy,iy
    real     :: CC(:,:,:),RC(:,:,:),rays(:,:,:,:),gnorm(:,:,:)
    real     :: x,y,z,xx,yy,zz,ddx,ddy,ddz,r,s,dr,fx,fz,fy,pt
    real     :: w000,w001,w010,w011,w100,w101,w110,w111,sumw
    logical :: adj

    !! For OMP
    integer, external :: omp_get_num_threads,omp_get_thread_num
    integer  :: nth,id

    real, allocatable, dimension(:,:,:,:) :: CCtmp

    !$OMP PARALLEL
    nth = omp_get_num_threads()
    !$OMP END PARALLEL

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

    if (adj) then !! . . Interpolating RWE wavefields -> Cart

       allocate ( CCtmp(cnz,cnx,cny,nth) )
       !$OMP PARALLEL DO PRIVATE(ie3,ie2,ie1,x,y,z,ix,iy,iz,pt,id,xx,ddx,yy,ddy,zz,ddz,r,s)
       E3Loop: do ie3=1,rnz
          E2Loop: do ie2=1,rny
             E1Loop: do ie1=1,rnx

                !! Check X
                x=rays(ie1,ie2,ie3,1)
                ix=floor((x-xmin)/dxx)+1

                !! Check y
                y=rays(ie1,ie2,ie3,2)
                iy=floor((y-ymin)/dyy)+1

                !! Check z
                z=rays(ie1,ie2,ie3,3)
                iz=floor((z-zmin)/dzz)+1
 
                if(ix .lt. 1+nsx .or. ix .gt. cnx-nsx) cycle
                if(iy .lt. 1+nsy .or. iy .gt. cny-nsy) cycle
                if(iz .lt. 1+nsz .or. iz .gt. cnz-nsz) cycle
                
                pt = RC(ie1,ie2,ie3)*gnorm(ie1,ie2,ie3)**2
                id = omp_get_thread_num()+1

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

                   X2Loop: do jy=iy-nsy,iy+nsy
                      yy=ymin+(jy-1.)*dyy
                      ddy=yy-y

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

                         r = sqrt( ddz**2 + ddy**2 + ddx**2 ) / dr
                         s = 1.! min(sin(r)/r,1.0)
                         CCtmp(jz,jx,jy,id)=CCtmp(jz,jx,jy,id)+pt*s +10.

                      end do X3Loop
                   end do X2Loop
                end do X1Loop

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

       do id=1,nth
          !$OMP PARALLEL DO PRIVATE(ix,iy,iz)
          do iy=1,cny
             do ix=1,cnx
                do iz=1,cnz
                   CC(iz,ix,iy) = CC(iz,ix,iy) + CCtmp(iz,ix,iy,id)
                end do
             end do
          end do
          !$OMP END PARALLEL DO
       end do

       deallocate ( CCtmp )

    else  !! . . Interpolating Cartesian velocity models -> RWE
       !$OMP PARALLEL DO PRIVATE(ie3,ie2,ie1,x,y,z,ix,iy,iz,jx,jy,jz,fx,fy,fz,&
                  w000,w001,w010,w011,w100,w101,w110,w111,sumw) 

       do ie3=1,rnz
          do ie2=1,rny
             do ie1=1,rnx
                x=rays(ie1,ie2,ie3,1)
                y=rays(ie1,ie2,ie3,2)
                z=rays(ie1,ie2,ie3,3)
                x = min(x,xmax)
                x = max(x,xmin)
                y = min(y,ymax)
                y = max(y,ymin)
                z = min(z,zmax)
                z = max(z,zmin)

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

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

                iy = min(iy  ,cny)
                iy = max(iy  , 1  )
                jy = min(iy+1,cny)

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

                fz= z - zmin-(iz-1)*dzz
                fx= x - xmin-(ix-1)*dxx
                fy= y - ymin-(iy-1)*dyy

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

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

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

                !! . . Calculate weights
                w000 = (1.-fz)*(1.-fx)*(1.-fy)
                w001 = (1.-fz)*(1.-fx)*(   fy)
                w010 = (1.-fz)*(   fx)*(1.-fy)
                w100 = (   fz)*(1.-fx)*(1.-fy)
                w101 = (   fz)*(1.-fx)*(   fy)
                w110 = (   fz)*(   fx)*(1.-fy)
                w011 = (1.-fz)*(   fx)*(   fy)
                w111 = (   fz)*(   fx)*(   fy)
                sumw = w000+w001+w010+w011+w100+w101+w110+w111

                RC(ie1,ie2,ie3) = RC(ie1,ie2,ie3) &
                +( w000 * CC(iz,ix,iy) + w100 * CC(jz,ix,iy) &
                +  w010 * CC(iz,jx,iy) + w110 * CC(jz,jx,iy) &
                +  w001 * CC(iz,ix,jy) + w101 * CC(jz,ix,jy) &
                +  w011 * CC(iz,jx,jy) + w111 * CC(jz,jx,jy) ) /sumw 
             end do
          end do
       end do
       !$OMP END PARALLEL DO
    end if
  end subroutine C2R3D_CART


end module C2R3D
