module C2R2Dold

  implicit none

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

contains

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

 subroutine C2R_init(insx,insz,inorm,idxx,idzz,&
                     icnx,icnz,irnx,irnz,ixmin,izmin)
    integer :: insx,insz,inorm,icnx,icnz,irnx,irnz
    real    :: idxx,idzz,ixmin,izmin
    nsx = insx;  nsz = insz;   norm= inorm
    xmin=ixmin;  zmin= izmin
    cnx = icnx;  cnz = icnz
    rnx = irnx;  rnz = irnz
    dxx = idxx;  dzz = idzz
    
    xmax=(cnx-1)*dxx+xmin
    zmax=(cnz-1)*dzz+zmin
    call auxputch("n1","i",cnz,"Sline.H")
    call auxputch("o1","r",0.,"Sline.H")
    call auxputch("d1","r",1.,"Sline.H")

  end subroutine C2R_init

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

  subroutine C2R_run(RC,rays,CC,gnorm,adj)
    real    :: CC(:,:),RC(:,:),rays(:,:,:),gnorm(:,:)
    real,allocatable:: fold(:,:)
    integer :: ix,iz,it,ig,im,jx,jz
    real    :: x,z,zz,xx
    real    :: ddz,ddx,r,sincr,dr
    logical :: adj

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


    if (adj) then
       allocate(fold(cnz,cnx))
    else
       allocate(fold(rnz,rnx))
    end if

    fold=0.

!    if (adj) write(0,*) 'MIN/MAX Image to interp',minval(RC),maxval(RC)
!    write(0,*) ' RAY DIM',rnz,rnx
!    write(0,*) 'CART DIM',cnz,cnx
!    if (adj) write(0,*) 'MIN/MAX Gnorm',minval(gnorm),maxval(gnorm)
!    if (adj) write(0,*) 'MIN/MAX CC',minval(CC),maxval(CC)
    dr = sqrt( dxx**2 + dzz**2 )

    if (adj) then
       do ig=1,rnx
          do it=1,rnz   
             x=rays(it,ig,1)
             z=rays(it,ig,2)
             
             iz=floor((z-zmin)/dzz)+1
             ix=floor((x-xmin)/dxx)+1
             
             if(iz<   1+nsz) cycle
             if(iz> cnz-nsz) cycle
             if(ix<   1+nsx) cycle
             if(ix> cnx-nsx) cycle
             
             do jz=iz-nsz,iz+nsz
                zz=zmin+(jz-1)*dzz
                ddz=zz-z
                
                do jx=ix-nsx,ix+nsx
                   xx=xmin+(jx-1)*dxx
                   ddx=xx-x
                   
                   r=sqrt( ddz**2 + ddx**2 )
                   
                   if(abs(r)< 0.0001) then
                      sincr = 1.0
                   else
                      r = r / dr
                      sincr = sin(r)/r
                   end if
                   
                   CC(jz,jx)=CC(jz,jx)+RC(it,ig)*sincr!*sqrt(abs(gnorm(it,ig)))
                   fold(jz,jx) = fold(jz,jx) + 1.
                    
                end do
             end do
             
          end do !! End Tau
       end do !! End Mu

       where (fold .gt. 0.)
          CC=CC /fold
       end where

    else
       do ig=1,rnx
          do it=1,rnz

             x=rays(it,ig,1)
             z=rays(it,ig,2)

             z = min(z,zmax)
             z = max(z,zmin)
             x = min(x,xmax)
             x = max(x,xmin)
             
             iz=floor((z-zmin)/dzz)+1
             ix=floor((x-xmin)/dxx)+1
             
             iz = min(iz  ,cnz)
             iz = max(iz  , 1  )
             jz = min(iz+1,cnz)
             
             ix = min(ix  ,cnx)
             ix = max(ix  , 1  )
             jx = min(ix+1,cnx)
             
             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)
             
             !! 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

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

             w00=0.;w01=0.;w10=0.;w11=0.;sumw=0.
          end do
       end do


    end if

  deallocate ( fold )
  end subroutine C2R_run

  !----------------------------------------------------------------
  subroutine PixelBalance(velvel)
    integer   :: iz,ix
    real      :: velvel(:,:)

    do iz=2,cnz-1
       do ix=2,cnx-1
          if ( &
          abs( velvel(iz,ix  )- velvel(iz-1,ix  ))/velvel(iz,ix ) .gt. 0.10 .and. &
          abs( velvel(iz,ix  )- velvel(iz+1,ix  ))/velvel(iz,ix ) .gt. 0.10 .and. &
          abs( velvel(iz,ix  )- velvel(iz  ,ix-1))/velvel(iz,ix ) .gt. 0.10 .and. &
          abs( velvel(iz,ix  )- velvel(iz  ,ix+1))/velvel(iz,ix ) .gt. 0.10 ) then
             velvel(iz,ix) = 0.25*(velvel(iz+1,ix)+velvel(iz-1,ix)+velvel(iz,ix-1)+velvel(iz,ix+1))
  !           write(0,*) 'FIXED PIXEL'
         end if
       end do
    end do

  end subroutine PixelBalance

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

  subroutine CartSurface_filter(gtmp,www)
    integer :: iz,ix,ii
    real   :: www,gtmp(:,:),tmp
    real, allocatable:: Sline(:),Stt(:)

    allocate(Sline(cnz),Stt(cnz))

    Sline=1.
    Stt=0.
 
    tmp=1500./sqrt(sqrt(www))!!was 2100
    write(0,*) '-  Mininum Depth is: ',tmp

    do iz=1,cnz
 
      if (dble(iz-1.)*dzz+zmin .lt. -tmp ) then
          Stt(iz) = 1.
       else
          Stt(iz) = 0.
       end if
   !    write(0,*) iz,(dble(iz-1.)*dzz+zmin) ,tmp ,Stt(iz)
    end do

    do ii=1,11 !! Smoothing Iterations
       Sline(1 ) = Stt(1   )
       Sline(cnz) = Stt(cnz-1)
       do iz=2,cnz
          Sline(iz)=( Stt(iz-1)+Stt(iz) )/2.
       end do
       Stt=Sline
    end do

    !call srite("Sline.H",Sline,size(Sline)*4)

    do ix=1,cnx
       gtmp(:,ix) = gtmp(:,ix)*Sline
    end do
    
    

    deallocate( Sline,Stt )
  end subroutine CartSurface_filter

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




end module C2R2Dold
