#<
#shoot.x
#
#Usage:
#shoot.x  <in.H >out.H 
#
#Input Parameters:
# 
# stdin   -   sepfile   Velocity file
# nangle  -   integer   Number of angles to shoot
# verb    -   integer   Verbosity level
# output_tau- integer   [0] wheter (1) or not (0) to output rays,etc in tau
# sources  -  sepfile   file containg source locations (optional)
# oxsrc    -  real      begining source location x
# dxsrc    -  real      sources sampling x
# ozsrc    -  real      begining source location z
# dzsrc    -  real      sources sampling z
# nsrc     -  integer   number of sources
# angles   -  sepfile   file containing angle range (complex (first,last))
# oang     -  real      first angle to shoot
# dang     -  real      sampling in angle
# nang     -  integer   number of angles
# nt       -  integer   number of time samples
# ot       -  real      first time sample
# dt       -  real      tme sampling
# zero_offset-integer   [0] no, [1] down, -1 up
# j1      -   integer   [1]  subsampling of the first axis to write out
#
#
#Output Parameters:
# stdout  -   sepfile   Sepfile contain the rays
# raypar  -   sepfile   (optional) Ray parameter
# raymod  -   sepfile   (optional) Ray mod
#
#
#
#Description:
#
#
#
#>
#%end of self-documentaiton
#-------------------------------------------------
#
# THIS IS PAUL'S CODE, F90 version
#
#Author: Robert Clapp, ESMB 463, 7230253
#
#Date Created:Fri Aug 21 13:40:16 PDT 1998
#
#Purpose: 
#

program shoot{ 
	use shoot_mod
  use tau_z_conv_mod
  integer :: i1,i2,i3,verb
	real     :: sx,sz
	type(axis) :: vel_axis(2),ray_axis(2)
  real,dimension(:,:), allocatable :: velocity,map
  real(kind=wp),dimension(:,:), allocatable :: x_ray,z_ray,v_ray
  complex,dimension(:,:), allocatable :: boundary,raypar
	complex,dimension(:) ,pointer :: sources
	complex, dimension(:) ,pointer :: angles
	integer                         ::  nloc,scale,j1
	logical                         :: rite_raypar,rite_raymod,output_tau


	if(verb>1) write(0,*) "before get_params"
	call get_shot_params(vel_axis,verb,rite_raypar,rite_raymod,output_tau,j1)

	if(verb>1) write(0,*) "after get_params"
	call init_shot_params(vel_axis,nloc,verb,sources,angles,j1)

	if(verb>1) write(0,*) "after init_params"
	call rite_shot_params(ray_axis(2)%n,nloc,rite_raypar,rite_raymod)

	if(verb>1) write(0,*) "after rite_params"
	allocate(velocity(vel_axis(1)%n,vel_axis(2)%n))

	call read_vel(velocity)

	if(verb>0) write(0,*) "before allocates",nt,nangle
	allocate(x_ray(nt,nangle))
	allocate(z_ray(nt,nangle))
	allocate(v_ray(nt,nangle))

	if(rite_raypar)
		allocate(raypar(nt,nangle))

	if(output_tau){
		if(verb>0) write(0,*) "Output rays, etc in tau"
		allocate(map(vel_axis(1)%n,vel_axis(2)%n))
		call create_tau_map(velocity,map)
	}
		

	do i1=1,nloc{

		
		sx=real(sources(i1)); sz=aimag(sources(i1))
		if(verb>1) write(0,*) "before shooting",sx,sz
		call shoot_rays_new(x_ray,z_ray,v_ray,velocity,sources(i1),angles(i1))

		call srite("out",cmplx(x_ray(::j1,:),z_ray(::j1,:)),nangle*ntout*8)

		if(rite_raymod)
			call srite("raymod",cmplx(real(v_ray(::j1,:)),0.),ntout*nangle*8)

		if(rite_raypar){
			call calc_raypar(x_ray,z_ray,v_ray,nangle,raypar)
			call srite("raypar",raypar(::j1,:),8*nangle*ntout)
		}
			

		if(output_tau)
			call convert_z_tau(x_ray,z_ray)
		

	}
}