# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
program ab
  use sep3d_struct_mod
  use sep
  use calc_ray_mod
  implicit none
  type(sep3d) :: vel_s,ray_s
  real ::a1
  real, allocatable :: vel(:,:),ray(:,:),pray(:,:)
  complex, allocatable :: cray(:)
  integer :: i2,n,i
  call init_3d()
  call init_sep3d("in",ray_s,"INPUT")
  call init_sep3d("vel",vel_s,"INPUT")
  allocate(vel(vel_s%n(1),vel_s%n(2)),cray(ray_s%n(1)),ray(3,ray_s%n(1&
    &)),pray(3,ray_s%n(1)))
  call to_history("esize",4)
  call sreed("vel",vel,size(vel)*4)
  write(0,*) "CHECK MIN MAX",minval(vel),maxval(vel),size(cray),size&
    &(pray,2),ray_s%n
  do i2=1,ray_s%ntraces
    call sreed("in",cray,size(Cray)*8)
    if (.not. convert_3d(cray,ray)) then
      write(0,*) "trouble converting ray 3d ",i2
      return
    end if
    n=tb_length(ray)
    if (.not. get_vel(vel,vel_s,ray(:,1:n),pray(1,1:n))) then
      write(0,*) "trouble getting ray parameter"
      return
    end if
    pray(1,n+1:)=0.
    do i=1,n
      if (pray(1,i) >.001) then
        pray(1,i)=1./pray(1,i)
      end if
    end do
    a1=sum(pray(1,:))/n
    write(0,*) "a1",a1
    call srite("out",pray(1,:),size(pray,2)*4)
  end do
end program
