program RWE2D_WI
  use sep
  use util_
  use RWE2D_WI_routines

  implicit none
  !! ax is the Cartesian length dimension associated with Velocity 
  !! az is the Cartesian depth  dimension associated with Velocity 
  !! aw is the frequency axis
  !! as is the shot axis
  !! ag is the offset axis in Ray space
  !! at is the extrapolation direction in Ray space
  !! ar is the number of references
  type(myaxis) :: ax,az,aw,ar,as,at,ag

  integer :: ix,iz,ig,iw,it,is
  real    :: jscale,pi,w,OF
  real,    allocatable,dimension(:,:)  :: CVel
  real,    allocatable,dimension(:,:)  :: gradient
  real,    allocatable,dimension(:,:,:):: raysin
  complex, pointer    ,dimension(:,:,:):: data,swf
  character(len=128) :: name
  logical :: OF_req,GR_req

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

  !! Get axes information
  call Get_WI_axes(ag,aw,as,at,az,ax)

  !! Print out info
  call routine_info(at,ag,as,aw,ax,az)

  !! Input data files
  allocate( swf(ag%n,aw%n,as%n),data(at%n,aw%n,as%n) )
  allocate( CVel(at%n,ag%n),gradient(at%n,ag%n), raysin(at%n,ag%n,2) )

  !! . . Read in Raysss and Velocity
  call sep_read(raysin,"rays")
  call sep_read(CVel ,"vel ")
  write(0,*) 'VELOCITY MIN/MAX', minval(CVel),maxval(CVel)

  !! read in Source and data wavefields
  call sep_read(swf ,"swf")
  call sep_read(data,"data")

  !! . . data scaling
  call from_param("jscale",jscale,1.)
  pi = acos(-1.)

  !! Scale Data
  write(0,*) 'SCALING DATA'
  do ig=1,ag%n
     do iw=1,aw%n
        do is=1,as%n
           data(ig,iw,is)=data(ig,iw,is)/&
           maxval(abs(data(:,iw,is)))/sqrt( dble(iw-1)*aw%d+aw%d)*jscale
           swf (ig,iw,is)= swf(ig,iw,is)/&
           maxval(abs( swf(:,iw,is)))/sqrt( dble(iw-1)*aw%d+aw%d)*jscale
        end do
     end do
  end do


  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!
  !! . . MAIN PROGRAM
  !!
  !! Calculate the coefficients for the coordinate system for a 
  !! Constant velocity medium

  !! Initialize Gradient and Objective Function calculators
  write(0,*) 'INITIALIZING ROUTINES'
  call routines_init(ag,at,ar,as,ax,az,aw,raysin)

  !! Find location on grid
  write(0,*) 'GETTING GRID LOCATION'
  call ray_location(raysin)

  OF_req=.true.
  GR_req=.false.

  freq: do iw=1,aw%n !! Frequency Loop

     !! Frequency that we are doinggg
     w = 2.*pi*( (iw-1.)*aw%d+aw%o )

     !! If we need to calculate the Objective Function
     if (OF_req) write(0,*) 'Calculating Objective Function'
     if (OF_req) call Calc_ObjFun(CVel,raysin,data(:,iw,:),swf(:,iw,:),w,OF)

     !! If we need to calculate the Gradient
     if (GR_req) write(0,*)  'Calculating Gradient'
     if (GR_req) call Calc_Gradient(CVel,raysin,data(:,iw,:),swf(:,iw,:),w,gradient,OF)

     !call lpfs(  )

  end do freq !! End Frequency loop

  call sep_write(CVel,"image")

  call sep_end_prog()
  call exit()

end program RWE2D_WI
