program RWE2D_ZERO
  !! . . Riemannian Wavefield Extrapolation program for migration 
  !! . . in 2D generalized coordinates. 
  !!
  !!		The key concept is that wave-equation migration can be performed
  !!		on almost any type of grid called Riemannian coordintes
  !!		This program allows for user to input 
  !!		generalized source and receiver wavefields that are recorded on the 
  !!		first extrapolation step, and then propagate them into the user-specified
  !!		computational mesh.
  !!
  !! . . General interal processing steps
  !! 1) Input wavefields, coordinate system, Velocity model
  !! 2) Calculate geometrical coefficients for coordinate system and velocity model
  !! 3) Propagate wavefields and apply imaging condition
  !! 4) Output Image as: a) interpolated to Cartesian; and b) In original coordinates
  !!
  !!
  !! 	. . Sample call in Makefile
  !! 	Image.H:  rwf.H  swf.H
  !!		$(BINDIR)/RWE2D.x rwf="rwf.H" swf="swf.H" vel="vel.H"\
  !!			rays="rays.H" image="Image.H" Rimage="RayImage.H" \
  !!			nref=256 verbose=1  \
  !!			min_region_pct=1. del_dist=0.05 verb=1 \
  !!			niter_lloyd=25 ntap=10 nloop=4 migapp=1024 forward=1 \
  !!			norm=1 zmin=-2000 zmax=6000 xmin=0 xmax=20000 \
  !!			dxx=10 dzz=10 nsx=2 nsz=2 
  !!
  !!		where requisite files:
  !!
  !!		1) "rwf" tag is the input receiver wavefield in nx-nw-ns format
  !!		2) "swf" tag is the input source wavefield in nx-nw-ns format
  !! 		3) "vel" tag is the input velocity model in nz-nx format
  !!		4) "rays" tag is the input coordinate system in nz-nx-2 format
  !!		where nz-nx-1 is the x-coordinate field and nz-nx-2 the z-coordinate field
  !!		5) "image" tag is the Cartesian output imagein cnz-cnx-ns format
  !!		6) "Riamge" tag is the Ray based ouput image in nz-nx-ns formate
  !!
  !!		Required parameters:
  !!
  !!		1) nref=xx is the starting number of points for calculating reference velocities
  !!		2) verbose=1 is the level of verbosity for the RWE2D migration code
  !!		3) min_region=.1 controls how close a true coefficient is to a reference value 
  !!				before it is thrown away (in %)
  !!		4) del_dist=0.05 is the minimum distance between reference parameters (in %)
  !!		5) verb=1 is the verbage level of the Lloyd's code
  !!		6) niter_lloyd=25 the number of iterations of Lloyd's code
  !!		7) nloop=4 number of smoothing iterations for
  !!                    merging wavefields from different references
  !!		8) migapp=1024 migration aperature split evenly about shotpoint
  !!		9) forward=1 Forward/Backscattering switch (0/1)
  !!		10) norm=1 Normalize by Jacobian of matrix
  !!		11/12) zmin/zmax min/max output Cartesian image z-locations
  !!		13/14) xmin/xmax min/max output Cartesian image x-locations
  !!		15/16) dzz/dxx grid spacing in output Cartesian image
  !!
  !! 		Written by Jeff Shragge
  !!				May 15/2006
  !!				contact: jshragge@gmail.com
  !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  !! Required modules
  use sep                !! Basic SEPlib I/O 
  use util_              !! handles axis I/O and aux file I/O
  use RWE2D_ZERO_wemig  !! Migration driver module
  use C2R2D              !! Ray-to-Cartesian interpolation 
  implicit none
  type(myaxis) :: ar,at,aw,an,as,ax,az !! axes types 

  integer :: ix,iz,ii,iw,iis,mxl,mxr,oldn,ir,it,ishot,rrn
  integer :: nref,nsx,nsz,norm,cnx,cnz,migapp,i1,nr,nt,nz,nx
  real    :: forward,xmin,zmin,xmax,zmax,dxx,dzz,r0,rrd,rro,rscale
  real    :: dr,dt,or,ot,dx,dz,ox,oz,minang,maxang,pi,x,z,dx1,dx2,dd1,dd2
  real,    allocatable,dimension(:,:,:):: raysin    !! Input rays		
  real,    allocatable,dimension(:,:)  :: Velin,vel,gnorm, Velsave !! Velocity model, normalization matrix
  complex, allocatable,dimension(:,:)  :: swfsave,rwfsave
  complex, allocatable,dimension(:,:  ):: rwf, swf  ! receiver and source wavefiels
  real,    allocatable,dimension(:)    :: rline,tline
  real,    allocatable,dimension(:,:)  :: CCimg, Pimg, Image !! Cartesian and ray coordinate images
  logical                              :: verbose,kin,adj


  ! call sep_init(SOURCE)
  call sep_init()
  !! Routines to enable parallel processing with
  !! SEPlib "Parallel" command
  call getch_add_string("head=/dev/null")
  call getch_add_string("noheader=y")
  call set_no_putch()
  call sep_begin_prog()

  !! 		Read in input dimensions from the "rwf" tag
  !!		Determines the radial steps grid parameters
  call from_aux("rwf","n1",ar%n)
  call from_aux("rwf","d1",ar%d)
  call from_aux("rwf","o1",ar%o)
  nr=ar%n;or=ar%o;dr=ar%d

  !! 		Determine the frequencies parameters
  call from_aux("rwf","n2",aw%n)
  call from_aux("rwf","o2",aw%o)
  call from_aux("rwf","d2",aw%d)

  !!		Determine the number of shot parameters
  call from_aux("rwf","n3",as%n)
  call from_aux("rwf","d3",as%d)
  call from_aux("rwf","o3",as%o)

  !!    Read in velocity file
  call from_aux("vel","n1",az%n)
  call from_aux("vel","d1",az%d)
  call from_aux("vel","o1",az%o)
  call from_aux("vel","n2",ax%n)
  call from_aux("vel","d2",ax%d)
  call from_aux("vel","o2",ax%o)
  nx=ax%n;ox=ax%o;dx=ax%d
  nz=az%n;oz=az%o;dz=az%d

  !! 		Parameter input section.  
  !!	 	Call from command line, Makefile or a "par=pars.P" file
  call from_param("forward",forward,0.)      !! Forward scattering option
  call from_param("nref",nref,16)           !! Max # reference velocities
  call from_param("verbose",verbose,.false.) !! verbosity
  call from_param("kinematic",  kin,.true.)  !! Kinematic approximation 
  call from_param("norm",norm,1)             !! Whether (1) or not (0) to normalize by gnorm
  call from_param("nsx",nsx,2)
  call from_param("nsz",nsz,2)
  call from_param("dt",dt,0.003)
  call from_param("nt",nt,500)
  call from_param("ot",ot,0.) 
  call from_param("minang",minang, 30.)
  call from_param("maxang",maxang,150.)
  call from_param("rscale",rscale,1.  )
  at%n=nt;  at%o=ot;  at%d=dt;
  pi = acos(-1.)

  rrd = ar%d
  rro = ar%o
  rrn = ar%n

  !!r0=rscale*ar%o; 
  r0= rscale*( (ar%n-1.)*ar%d )/2.

  write(0,*) r0,ar%n,ar%d
  dr   = (maxang-minang)/180.*pi/dble(nr-1)
  ar%d = (maxang-minang)/180.*pi/dble(nr-1)
  ar%o = pi*minang/180.

  write(0,*) 'EXTRAPOLATION STEPS',nt,ot,dt
  write(0,*) 'OFFSETS' ,ar%n,ar%o,ar%d
  write(0,*) 'SHOTS' ,as%n,as%o,as%d
  write(0,*) 'FREQ ' ,aw%n,aw%o,aw%d

  !! Create output header file associated with "image" tag
  !! . . Output Image in Cartesian coordinates
  call auxputch("n1","i",az%n,"image")
  call auxputch("n2","i",ax%n,"image")
  call auxputch("d1","r",az%d,"image")
  call auxputch("d2","r",ax%d,"image")
  call auxputch("o1","r",az%o,"image")
  call auxputch("o2","r",ax%o,"image")
  call auxputch("esize","i",4,"image")

  !! Create output header file associated with "image" tag
  !! . . Output Image in RWE coordinates
 ! if (verbose) then
     call auxputch("n1","i",at%n,"Rimage")
     call auxputch("n2","i",ar%n,"Rimage")
     call auxputch("n3","i",as%n,"Rimage")
     call auxputch("d1","r",at%d,"Rimage")
     call auxputch("d2","r",ar%d,"Rimage")
     call auxputch("o1","r",at%o,"Rimage")
     call auxputch("o2","r",ar%o,"Rimage")
     call auxputch("esize","i",4,"Rimage")
  !end if

!  if (verbose) then
     call auxputch("n1","i",at%n,"rays.H")
     call auxputch("n2","i",ar%n,"rays.H")
     call auxputch("d1","r",at%d,"rays.H")
     call auxputch("d2","r",ar%d,"rays.H")
     call auxputch("o1","r",at%o,"rays.H")
     call auxputch("o2","r",ar%o,"rays.H")
     call auxputch("esize","i",8,"rays.H")

     call auxputch("n1","i",at%n,"AAA.H")
     call auxputch("n2","i",ar%n,"AAA.H")
     call auxputch("d1","r",at%d,"AAA.H")
     call auxputch("d2","r",ar%d,"AAA.H")
     call auxputch("o1","r",at%o,"AAA.H")
     call auxputch("o2","r",ar%o,"AAA.H")
     call auxputch("esize","i",4,"AAA.H")
!  end if


  !!	Allocate required arrays
  allocate( rwf(ar%n,aw%n), swf(ar%n,aw%n), swfsave(ar%n,aw%n),rwfsave(ar%n,aw%n)  )
  allocate( raysin(nt,ar%n,2),Velin(az%n,ax%n),vel(nt,ar%n),gnorm(nt,ar%n),Velsave(az%n,ax%n) )
  allocate( CCimg(az%n,ax%n) , Pimg(nt,ar%n), Image(az%n,ax%n) )
  allocate( tline(nt),rline(nr) )

  !! Zero all arrays
  gnorm=0.; swfsave=0.;  Image=0.; Velsave=0.

  !! Read in Velocity file
  call sep_read(Velsave,"vel ")
  Velin=Velsave 
  write(0,*) 'READ IN VELOCITY FILE'

  !! Make rays for interpolation
  rline=(/ (i1,i1=0,nr-1) /)* ar%d + ar%o 
  tline=(/ (i1,i1=0,nt-1) /)* dt   + ot
  write(0,*) 'MIN/MAX Angles: ',minval(rline)*180./pi,maxval(rline)*180./pi
  write(0,*) 'MIN/MAX Ex step:',minval(tline),maxval(tline)

  !! Generate Coordinate System
  do ir=1,nr
     do it=1,nt
        raysin(it,ir,1) = -abs(r0)*cosh(tline(it))*cos(rline(ir)) +(as%o+(r0+rro))
        raysin(it,ir,2) =  abs(r0)*sinh(tline(it))*sin(rline(ir))
        gnorm (it,ir  ) =r0*r0*(sinh(tline(it))**2+sin(rline(ir))**2) 
     end do
  end do
  call sep_write(cmplx(raysin(:,:,1),raysin(:,:,2)),"rays.H")

  call sep_write(gnorm,"AAA.H")
  write(0,*) "Rline min/max   :", minval(rline),maxval(rline)
  write(0,*) "X-coord min/max :", raysin(1,1,1),raysin(1,ar%n,1)
  write(0,*) "Z-coord min/max :", minval(raysin(:,:,2)),maxval(raysin(:,:,2))


  !!----------------------------------------------------------------------------------------
  !!
  !! MIGRATION SECTION
  !!
  !! . . Initialize module and image space
  !! . . Initialize the wave-equation migration driver
  write(0,*) 'INITALIZING WEMIG'
  call wemig_init(ar,at,aw,forward,kin,nref,verbose,r0,minang,maxang)

  do ishot=1,as%n   !! . . Loop over shot points
     write(0,*) 'SHOT POINT AT ', as%o+(ishot-1)*as%d
     write(0,*) 'First Receiver at ',raysin(1,1,1)
     write(0,*) 'Last Receiver at  ',raysin(1,ar%n,1)

     vel=0.; rwf=0.;  Velin=0.; CCimg=0.;  Pimg=0.

     Velin = Velsave

     !! Do Cart -> RWE velocity Interpolation
     adj=.false.
     call C2R_init(0,0,norm,dx,dz,nx,nz,nr,nt,ox,oz)
     call C2R_run(vel,raysin,Velin,gnorm,adj)  

     !! Read in next Receiver wavefield
     call sep_read(rwfsave ,"rwf ") !! Read in only current receiver wavefield

     !! Position wavefield on surface
     do ir = 1,nr
        x = raysin(1,ir,1) - as%o - (ishot-1)*as%d  
        ix = floor( ( x - rro ) / rrd) + 1  
        dx1 = x - (rro + dble(ix-1)*rrd )
        dx2 = (rro + dble(ix  )*rrd )-x
        dd1 = dx2/(dx1+dx2)
        dd2 = dx1/(dx1+dx2)
        if (ix .ge. nr) ix = nr-1
        if (ix .le. 0 ) ix = 1 
        rwf(ir,:) = dd1*rwfsave(ix,:)+dd2*rwfsave(ix+1,:)
     end do
     write(0,*) "MIN/MAX RWF: ", minval(real(rwfsave)),maxval(real(rwfsave))

     !! Do Imaging

     write(0,*) 'ishot in',ishot
     call wemig(rwf,Pimg,ishot,vel )   

     write(0,*) 'ishot out',ishot

     !! Do RWE -> Cart Image Interpolation
     call C2R_close()
     adj=.true.
     call C2R_init(nsx,nsz,norm,dx,dz,nx,nz,nr,nt,ox,oz)
     call C2R_run(Pimg,raysin,CCimg,gnorm,adj)
     call C2R_close()

     !! Add single image to global image
     do ix=1,ax%n
        do iz=1,az%n 
           Image(iz,ix) = Image(iz,ix) + CCimg(iz,ix) 
        end do
     end do

     !! Advance rays to next step
     do ir=1,ar%n
        do it=1,at%n
           raysin(it,ir,1)=raysin(it,ir,1)+as%d
        end do
     end do

     !! Write out RWE-based image
     call sep_write(Pimg,"Rimage")

     write(0,*) 'JUST DID SHOT ', ishot ,' ************************'
  end do !! . . End Shot-point loop

  !! Write out Cartesian-based image
  call sep_write(Image,"image")

  call sep_end_prog()
  call exit()
end program RWE2D_ZERO
