program RWE2D_source
  use sep
  use util_
  use RWE2D_source_wemig
  use C2R2D

  implicit none
  type(myaxis) :: ax,az,aw,ar,as

  integer :: ix,iz,ii,iw,ishot(1),iis,mxl,mxr,oldn
  integer :: nref,nsx,nsz,norm,cnx,cnz,migapp
  real    :: forward,xmin,zmin,xmax,zmax,dxx,dzz
  real,    allocatable,dimension(:,:,:):: raysin       !! Input Rayfield
  real,    allocatable,dimension(:,:)  :: Vel,gnorm    !! Velocity model  
  complex, pointer    ,dimension(:,:  ):: rwf
  real,    allocatable,dimension(:,:)  :: Pimg,CCimg   !! Single shot-image 
  logical                              :: verbose,kin
  character(len=128) :: name

  !! Geometry
  integer, allocatable :: mask(:,:),cutmask(:,:)
  real,    allocatable :: refs(:,:,:), fields(:,:,:),cutfields(:,:,:)

  call sep_init(SOURCE)
  call getch_add_string("head=/dev/null")
  call getch_add_string("noheader=y")
  call set_no_putch()
  call sep_begin_prog()

  call from_param("forward",forward,0.) !! Forward scattering option
  call from_param("nref",nref,1)        !! # reference velocities
  call from_param("verbose",verbose,.false.)
  call from_param("kinematic",  kin,.true.)
  call from_param("norm",norm,1)
  name="rays"

  !! Read in input dimensions
  call from_aux("rwf","n1",ax%n)
  call from_aux("rwf","n2",aw%n)
  call from_aux("rwf","n3",as%n)
  call from_aux("rwf","d1",ax%d)
  call from_aux("rwf","d2",aw%d)
  call from_aux("rwf","d3",as%d)
  call from_aux("rwf","o1",ax%o)
  call from_aux("rwf","o2",aw%o)
  call from_aux("rwf","o3",as%o)
  call from_aux("rays","n1",az%n)
  call from_aux("rays","d1",az%d)
  call from_aux("rays","o1",az%o)

  call from_param("migapp",migapp,ax%n)

  write(0,*) 'DEPTH',az%n,az%o,az%d
  write(0,*) 'X1   ',ax%n,ax%o,ax%d
  write(0,*) 'SHOTS',as%n,as%o,as%d
  write(0,*) 'FREQ ',aw%n,aw%o,aw%d

  !! Interpolation to Cartesian
  call from_param("xmin",xmin,ax%o)
  call from_param("zmin",zmin,az%o)
  call from_param("xmax",xmax,(ax%n-1)*ax%d+ax%o)
  call from_param("zmax",zmax,(az%n-1)*az%d+az%o)
  call from_param("dxx",dxx,ax%d)
  call from_param("dzz",dzz,az%d)
  call from_param("nsx",nsx,3)
  call from_param("nsz",nsz,3)
  cnx = floor((xmax - xmin)/dxx)+1
  cnz = floor((zmax - zmin)/dzz)+1

  !! . . Create offset axis
  call auxputch("n1","i",cnz ,"image")
  call auxputch("n2","i",cnx ,"image")
  call auxputch("n3","i",as%n,"image")
  call auxputch("d1","r",dzz ,"image")
  call auxputch("d2","r",dxx ,"image")
  call auxputch("d3","r",as%d,"image")
  call auxputch("o1","r",zmin,"image")
  call auxputch("o2","r",xmin,"image")
  call auxputch("o3","r",as%o,"image")
  call auxputch("esize","i",4,"image")

  !! . . Output Imaged WF in RWE coordinates
  call auxputch("n1","i",az%n,"Rimage")
  call auxputch("n2","i",ax%n,"Rimage")
  call auxputch("n3","i",as%n,"Rimage")
  call auxputch("d1","r",az%d,"Rimage")
  call auxputch("d2","r",ax%d,"Rimage")
  call auxputch("d3","r",as%d,"Rimage")
  call auxputch("o1","r",az%o,"Rimage")
  call auxputch("o2","r",ax%o,"Rimage")
  call auxputch("o3","r",as%o,"Rimage")
  call auxputch("esize","i",4,"Rimage")

  allocate( rwf(migapp,aw%n) )
  allocate( raysin(az%n,ax%n,2),Vel(az%n,ax%n),gnorm(az%n,ax%n) )
  allocate( mask(ax%n  ,az%n), refs(nref,4,az%n),fields(ax%n,4,az%n))
  allocate( cutmask(migapp,az%n), cutfields(migapp,4,az%n) )

  allocate( CCimg(cnz,cnx) , Pimg(az%n,migapp) )

  call sep_read(raysin,"rays")
  call sep_read(Vel   ,"vel ")


  write(0,*) 'VELOCITY MIN/MAX', minval(Vel),maxval(Vel)
write(0,*) ' GOT TO HERE'

  !! . . Calculate Geometry
  gnorm = 1.0
  call geometry_init(az%n,ax%n,az%d,ax%d,kin,nref,verbose)
  call geometry_calc(raysin,Vel,fields,refs,mask,gnorm)
  call geometry_close()
  deallocate( Vel )

  write(0,*) 'GNORM',minval(gnorm),maxval(gnorm)

  !! . . Initialize module and image space
  oldn=ax%n
  ax%n=migapp
  call wemig_init(ax,az,aw,ar,forward,kin,nref,verbose)
  write(0,*) 'Initialized WEMIG routine'

  call C2R_init(nsx,nsz,norm,dxx,dzz,cnx,cnz,migapp,az%n,xmin,zmin)
  write(0,*) 'Initialized C2R_INIT routine'

  do ii=1,as%n   !! . . Loop over shot points
     rwf   = 0.
     call sep_read(rwf ,"rwf ")
     Pimg = 0.; CCimg = 0.

     ishot = maxloc(cabs(rwf(:,1)));
     iis=ishot(1)
     mxl = max(   1,iis-migapp/2)
     mxr = min(oldn,iis+migapp/2)
     if (mxl .eq. 1   ) mxr=migapp
     if (mxr .eq. oldn) mxl=oldn-migapp+1
     
     write(0,*) 'ISHOT',iis,mxl,mxr

     write(0,*) 'RWF MIN/MAX', minval(abs(rwf)),maxval(abs(rwf))

     cutfields = fields(mxl:mxr,:,:)
     cutmask   = mask  (mxl:mxr,:  )
     call wemig(rwf,Pimg,ii,cutfields,refs,cutmask )
     call sep_write(Pimg,"Rimage")

write(0,*) 'Interpolating'
     call C2R_run(Pimg,raysin(:,mxl:mxr,:),CCimg,gnorm(:,mxl:mxr))
     call sep_write(CCimg,"image")

  end do !! . . End Shot-point loop

  call sep_end_prog()
  call exit()
end program RWE2D_source