program pts_rays
use sep
use sep3d_struct_mod
use shootit_mod
use basic_util_mod
use verb_mod
implicit none
type(sep3d) :: pts_s,rays_s,vel
real, pointer :: vel2d(:,:),vel3d(:,:,:)
integer :: n(3),ndim,nh,ipts,i,ih
real :: o(3),d(3)
integer :: nang,jtrace,itheta,iphi
real    :: oang,dang,maxt,dt
real,allocatable :: rays(:,:,:),times(:)
logical, allocatable :: valid(:)
real, allocatable :: phi(:),theta(:),temp(:)
real, allocatable :: posx(:),posy(:),posz(:),var(:),dipx(:),dipy(:),gamma(:,:),vel_ar(:)
real, pointer :: x(:),y(:),z(:),v(:),dx(:),dy(:),g(:),th(:),ti(:),ph(:),v_ar(:)
integer, allocatable:: grid(:),coord(:,:)
integer :: ipt


call sep_init(SOURCE)
call init_3d()



call init_sep3d("in",pts_s,"INPUT")
call init_sep3d("vel",vel,"INPUT")
call set_n_o_d_3d(vel,n,o,d)
ndim=sep3d_ndims(vel)


call sep3d_grab_headers("in",pts_s,nh)
if(.not. init_shot_pars("vel",nh,rays_s))&
  call seperr("trouble initializing 2d shot parameters")

allocate(rays(3,par%nt*2,par%nphi*par%ntheta))
allocate(valid(((par%nphi*par%ntheta)-1)*2+1))

!PHI AND THETA FOR RAY GROUP
i=0
allocate(theta(par%nphi*par%ntheta))
allocate(phi(par%nphi*par%ntheta))
do iphi=1,par%nphi
  do itheta=1,par%ntheta
    i=i+1
    phi(i)=par%ophi+par%dphi*(iphi-1)
    theta(i)=par%otheta+par%dtheta*(itheta-1)
  end do
end do 

allocate(var(nh),dipx(nh),dipy(nh),posx(nh),posy(nh),posz(nh),gamma(1,nh),vel_ar(nh))
call sep3d_grab_key_vals(pts_s,"vel",vel_ar)
call sep3d_grab_key_vals(pts_s,"var",var)
call sep3d_grab_key_vals(pts_s,"pos_x",posx)
call sep3d_grab_key_vals(pts_s,"pos_y",posy)
call sep3d_grab_key_vals(pts_s,"pos_z",posz)
call sep3d_grab_key_vals(pts_s,"dip_x",dipx)
call sep3d_grab_key_vals(pts_s,"dip_y",dipy)
if(.not. sep3d_read_data("in",pts_s,gamma)) &
  call seperr("trouble reading data")

if(.not. init_verbose()) &
  call seperr("trouble initializing verbosity")

if(.not. init_print_progress(0,nh)) &
  call seperr("trouble initializing print progress")

allocate(times(size(theta)))
allocate(grid(size(theta)))
do ipt=1,nh
  call print_progress("tracing rays ")
  if(.not. make_rays(posx(ipt),posy(ipt),posz(ipt),dipx(ipt),dipy(ipt),valid,times,rays))&
   call seperr("trouble making rays")
  i=0
  do ih=1,size(valid)
    if(valid(ih)) then
      i=i+1
      rays(:,:,i)=rays(:,:,ih) 
    end if
  end do
    allocate(v_ar(i),x(i),y(i),z(i),v(i),dx(i),dy(i),g(i),ph(i),ti(i),th(i),coord(3,i))
    x=posx(ipt)
    y=posy(ipt)
    z=posz(ipt)
    v=var(ipt)
    dx=dipx(ipt)
    dy=dipy(ipt)
    v_ar=vel_ar(ipt)
    g=gamma(1,ipt)
    i=0
    grid=-1
    do ih=1,size(valid)
      if(valid(ih)) then
        i=i+1
        ti(i)=times(ih)
        ph(i)=phi(ih)
        th(i)=theta(ih)
        grid(ih)=i
        coord(1,i)=ih
        coord(2,i)=1
        coord(3,i)=ipt
      end if
    end do
    call sep3d_set_number_headers(rays_s,i)
    if(i>0) then
      call sep3d_set_key_vals(rays_s,"vel",v_ar)
      call sep3d_set_key_vals(rays_s,"var",v)
      call sep3d_set_key_vals(rays_s,"gamma",g)
      call sep3d_set_key_vals(rays_s,"pos_x",x)
      call sep3d_set_key_vals(rays_s,"pos_y",y)
      call sep3d_set_key_vals(rays_s,"pos_z",z)
      call sep3d_set_key_vals(rays_s,"dip_x",dx)
      call sep3d_set_key_vals(rays_s,"dip_y",dy)
      call sep3d_set_key_vals(rays_s,"phi",phi(1:i))
      call sep3d_set_key_vals(rays_s,"theta",theta(1:i))
      call sep3d_set_key_vals(rays_s,"time",ti)
      if(.not. sep3d_set_coords(rays_s,coord)) &
        call seperr("trouble setting coords")
    end if
    call sep3d_inorder(rays_s)
    if(.not. sep3d_write_data("out",rays_s,rays,nwind=(/rays_s%n(1),rays_s%n(2),&
      rays_s%n(3),1/),fwind=(/0,0,0,ipt-1/)))&
     call seperr("trouble writing out data")
    deallocate(x,y,z,v,dx,dy,g,ph,ti,th,coord)
end do
if(.not. sep3d_update_ntraces(rays_s)) &
  call seperr("trouble updating ntraces")
call sep3d_write_description("out",rays_s)
end program
