program RWE3D_source
  use sep
  use util
  use RWE3D_source_wemig

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

  integer :: ix,iy,iz,ii,iw
  integer :: nref
  real    :: forward
  real,    allocatable,dimension(:,:,:,:,:):: raysin  !! Input Rayfield
  real,    allocatable,dimension(:,:,:,:)  :: rays    !! Single-shot rays
  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        !! Single shot-image 
  logical                              :: verbose,kin
  character(len=128) :: name

  call sep_init()
  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,.false.)
  name="rays"

  !! Read in input dimensions
  call hinaxis(ax,1)
  call hinaxis(ay,2)
  call hinaxis(aw,3) 
  call hinaxis(as,4)
  call ainaxis(az,1,name)

  if (verbose) then
     write(0,*) az%n,az%o,az%d
     write(0,*) ax%n,ax%o,ax%d
     write(0,*) ay%n,ay%o,ay%d
     write(0,*) as%n,as%o,as%d
     write(0,*) aw%n,aw%o,aw%d
  end if

  !! . . Create offset axis
  call houaxis(az,1)
  call houaxis(ax,2)
  call houaxis(ay,3) 
  call houaxis(as,4)
  call to_history("esize",4)

  !! . . Allocate input wavefields 
  allocate( rwfin(ax%n,ay%n,aw%n,as%n) , swfin(ax%n,ay%n,aw%n,as%n) )

  !! . . Allocate arrays for individual shots
  allocate(swf(ax%n,ay%n,aw%n) , rwf(ax%n,ay%n,aw%n) , Pimg(az%n,ax%n,ay%n) )

  !! . . Allocate coordinate system
  allocate(raysin(az%n,ax%n,ay%n,as%n,3), rays  (az%n,ax%n,ay%n,3))
  call sep_read(raysin,"rays")

  !! . . Allocate Velocity model (in ray coordinates)
  allocate( Vel(az%n,ax%n,ay%n) )
write(0,*) size(Vel,1),size(Vel,2),size(Vel,3)
  call sep_read(Vel,"vel")
  write(0,*) 'VELOCITY MIN/MAX', minval(Vel),maxval(Vel)

  !! . . Read in RWF and SWF
  call sep_read(rwfin)
  call sep_read(swfin,"swf")
  write(0,*) 'RWF MIN/MAX', minval(abs(rwfin)),maxval(abs(rwfin))
  write(0,*) 'SWF MIN/MAX', minval(abs(swfin)),maxval(abs(swfin))

  !! . . Initialize module and image space
  call wemig_init(ax,ay,az,aw,ar,forward,kin,nref,verbose)

  do ii=1,as%n   !! . . Loop over shot points
     Pimg=0.;rwf=0.;swf=0.;

     !! . . Individual Shots
     rwf  = rwfin(:,:,:,ii)
     swf  = swfin(:,:,:,ii)
     rays = raysin(:,:,:,ii,:)

     !! . . Image event
     call wemig(rwf,swf,Pimg,rays,Vel,ii)
     call sep_write(Pimg)

  end do !! Shot-point loop

  call sep_close()
  call exit()
end program RWE3D_source
