program RWE3D_SI
  use sep
  use util_
  use RWE3D_SI_wemig
  use geometry

  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(:,:,:)  :: Vel     !! Velocity model  
  complex, allocatable,dimension(:,:,:,:):: rwfin   !! input S&R wflds
  complex, pointer    ,dimension(:,:,:)  :: rwf     !! S&R wfld
  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.)
  name="rays"
write(0,*) 'GOT TO HERE'
  !! 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",1.)

  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)
  as%n=1.
  as%d=1.
  as%o=1
  if (verbose) then
     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
  end if

  !! . . Create offset axis
!  call auxputch("n1","i",ax%n,"SP")
!  call auxputch("n2","i",ay%n,"SP")
!  call auxputch("n3","i",aw%n,"SP")
!  call auxputch("n4","i",as%n,"SP")
!  call auxputch("d1","r",ax%d,"SP")
!  call auxputch("d2","r",ay%d,"SP")
!  call auxputch("d3","r",aw%d,"SP")
!  call auxputch("d4","r",as%d,"SP")
!  call auxputch("o1","r",ax%o,"SP")
!  call auxputch("o2","r",ay%o,"SP")
!  call auxputch("o3","r",aw%o,"SP")
!  call auxputch("o4","r",as%o,"SP")

  call auxputch("n1","i",ax%n,"SP")
  call auxputch("n2","i",ay%n,"SP")
  call auxputch("n3","i",az%n,"SP")
  call auxputch("n4","i",aw%n,"SP")

  call auxputch("d1","r",ax%d,"SP")
  call auxputch("d2","r",ay%d,"SP")
  call auxputch("d3","r",az%d,"SP")
  call auxputch("d4","r",aw%d,"SP")

  call auxputch("o1","r",ax%o,"SP")
  call auxputch("o2","r",ay%o,"SP")
  call auxputch("o3","r",az%o,"SP")
  call auxputch("o4","r",aw%o,"SP")

  call auxputch("esize","i",8,"SP")

  allocate(rwfin(ax%n,ay%n,aw%n,as%n),rwf(ax%n,ay%n,aw%n) )
  allocate(raysin(az%n,ax%n,ay%n,3) , Vel(az%n,ax%n,ay%n) )
  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 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( raysin, Vel )

  call wemig_init(ax,ay,az,aw,ar,forward,kin,nref,verbose)

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

     rwf=0.;
     rwf  = rwfin(:,:,:,ii)
     call wemig(rwf,ii,fields,refs,mask)
  end do !! . . End Shot-point loop
  call sep_end_prog()
  call exit()
  write(0,*) 'FINISHED PROGRAM'
end program RWE3D_SI
