program RWE3D
  use sep
  use util_
  use RWE3D_wemig
  use C2R

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

  integer :: ix,iy,iz,ii,iw
  integer :: nref,nsx,nsy,nsz,norm,cnx,cny,cnz
  real    :: forward,xmin,ymin,zmin,xmax,ymax,zmax,dxx,dyy,dzz
  real,    allocatable,dimension(:,:,:,:):: raysin      !! Input Rayfield
  real,    allocatable,dimension(:,:,:)  :: Vel         !! Velocity model  
  complex, allocatable,dimension(:,:,:,:):: rwfin,swfin !! input S&R wflds
  complex, pointer,    dimension(:,:,:)  :: rwf,swf     !! S&R wfld
  real,    allocatable,dimension(:,:,:)  :: Pimg,CCimg  !! Single shot-image 
  logical                                :: verbose,kin
  character(len=128) :: name

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

  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",ay%n)
  call from_aux("rwf","n3",aw%n)
  call from_aux("rwf","n4",as%n)
  call from_aux("rwf","d1",ax%d)
  call from_aux("rwf","d2",ay%d)
  call from_aux("rwf","d3",aw%d)
  call from_aux("rwf","d4",as%d)
  call from_aux("rwf","o1",ax%o)
  call from_aux("rwf","o2",ay%o)
  call from_aux("rwf","o3",aw%o)
  call from_aux("rwf","o4",as%o)
  call from_aux("rays","n1",az%n)
  call from_aux("rays","d1",az%d)
  call from_aux("rays","o1",az%o)

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


  call from_param("xmin",xmin,ax%o)
  call from_param("ymin",ymin,ay%o)
  call from_param("zmin",zmin,az%o)
  call from_param("xmax",xmax,(ax%n-1)*ax%d+ax%o)
  call from_param("ymax",ymax,(ay%n-1)*ay%d+ay%o)
  call from_param("zmax",zmax,(az%n-1)*az%d+az%o)
  call from_param("dxx",dxx,ax%d)
  call from_param("dyy",dyy,ay%d)
  call from_param("dzz",dzz,az%d)
  call from_param("nsx",nsx,3)
  call from_param("nsy",nsy,3)
  call from_param("nsz",nsz,3)

  cnx =       floor((xmax - xmin)/dxx)+1
  cny = max(1,floor((ymax - ymin)/dyy))
  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",cny,"image")
  call auxputch("n4","i",as%n,"image")
  call auxputch("d1","r",dzz,"image")
  call auxputch("d2","r",dxx,"image")
  call auxputch("d3","r",dyy,"image")
  call auxputch("d4","r",as%d,"image")
  call auxputch("o1","r",zmin,"image")
  call auxputch("o2","r",xmin,"image")
  call auxputch("o3","r",ymin,"image")
  call auxputch("o4","r",as%o,"image")
  call auxputch("esize","i",4,"image")

  allocate( rwfin(ax%n,ay%n,aw%n,as%n) , swfin(ax%n,ay%n,aw%n,as%n) )
  allocate(swf(ax%n,ay%n,aw%n),rwf(ax%n,ay%n,aw%n),Pimg(az%n,ax%n,ay%n) )
  allocate(raysin(az%n,ax%n,ay%n,3) )
  allocate( Vel(az%n,ax%n,ay%n),CCimg(cnz,cnx,cny) )
  allocate(mask(ax%n,ay%n,az%n),refs(nref,10,az%n),fields(ax%n,ay%n,10,az%n) )

  call sep_read(raysin,"rays")
  call sep_read(Vel,"vel")
  call sep_read(rwfin,"rwf")
  call sep_read(swfin,"swf")

  write(0,*) 'VELOCITY MIN/MAX', minval(Vel),maxval(Vel)
  write(0,*) 'RWF MIN/MAX', minval(abs(rwfin)),maxval(abs(rwfin))
  write(0,*) 'SWF MIN/MAX', minval(abs(swfin)),maxval(abs(swfin))

  !! . . Calculate Geometry

  call geometry_init(az%n,ax%n,ay%n,az%d,ax%d,ay%d,kin,nref,verbose)
  call geometry_calc(raysin,Vel,fields,refs,mask)
  call geometry_close()

  deallocate( Vel )

  !! . . Initialize module and image space
  call wemig_init(ax,ay,az,aw,ar,forward,kin,nref,verbose)
  call C2R_init(nsx,nsy,nsz,norm,dxx,dyy,dzz,cnx,cny,cnz,ax%n,az%n,ay%n,xmin,ymin,zmin)

  do ii=1,as%n   !! . . Loop over shot points

     Pimg=0.;rwf=0.;swf=0.;
     rwf  = rwfin(:,:,:,ii)
     swf  = swfin(:,:,:,ii)

     call wemig(rwf,swf,Pimg,ii,fields,refs,mask)

     call C2R_run(Pimg,raysin,CCimg)

     call sep_write(CCimg,"image")

  end do !! . . End Shot-point loop

  call sep_end_prog()
  call exit()
end program RWE3D
