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(:,:)
integer :: n(3),ndim,nh,ipts,i,ih,iin,iout
real :: o(3),d(3),orient
integer :: nang,jtrace,itheta,iphi,jt,i1,i2
real    :: oang,dang,maxt,dt,normal,z,timez
integer, allocatable:: grid(:),temp(:,:),coords(:,:)
integer :: ipt,nv(2),n_t,ineg,ipos
real :: ov(2),dv(2),d_t,max_t
logical :: found
real,allocatable :: rays(:,:,:),dip(:),baz(:),x(:),good(:),err(:,:),time(:),ray2(:,:,:)
real,allocatable :: elev(:),event(:),eout(:),stat(:)
real(kind=wp),allocatable  :: x_rays(:,:),z_rays(:,:),v_rays(:,:)




call sep_init(SOURCE)
call init_3d()



call init_sep3d("in",pts_s,"INPUT")
call init_sep3d("vel",vel,"INPUT")
call init_sep3d(pts_s,rays_s,"OUTPUT")
call from_param("orient",orient,25.)
call from_param("jt",jt,4)
call from_param("dt",d_t,.02)
call from_param("max_t",max_t,16.)
call from_param("nt",n_t,nint(max_t/d_t))
call init_shot(vel%n,vel%o,vel%d,n_t*jt,d_t/jt,1)

allocate(vel2d(vel%n(1),vel%n(2)))
rays_s%n(1)=n_t*3
rays_s%o(1)=0.
rays_s%d(1)=d_t
rays_s%keyname(8)="time"
rays_s%keyname(9)="err"
call sep3d_set_sep3d(rays_s)
z=(vel%n(1)-1)*vel%d(1)


call sreed("vel",vel2d,size(vel2d)*4)
call sep3d_grab_headers("in",pts_s,nh)
allocate(err(1,nh))
allocate(eout(nh))
if(.not. sep3d_read_data("in",pts_s,err)) &
  call seperr("trouble reading data")



allocate(x_rays(n_t*jt,1))
allocate(z_rays(n_t*jt,1))
allocate(v_rays(n_t*jt,1))

allocate(rays(3,n_t,nh))
allocate(good(nh),baz(nh),dip(nh),x(nh),grid(nh),time(nh))
allocate(temp(2,nh))
allocate(elev(nh),event(nh),stat(nh))
call sep3d_grab_key_vals(pts_s,"P",dip)
call sep3d_grab_key_vals(pts_s,"ELEV",elev)
call sep3d_grab_key_vals(pts_s,"EVENT_ID",event)
call sep3d_grab_key_vals(pts_s,"STATION_ID",stat)
call sep3d_grab_key_vals(pts_s,"BAZ",baz)
baz=baz-orient
where (baz <0) baz=baz+360.
where(baz <90. .or. baz > 270.) dip=-dip
call sep3d_grab_key_vals(pts_s,"good",good)
call sep3d_grab_key_vals(pts_s,"XSTA",x)
call sep3d_inorder(rays_s)

grid=-1
iin=0;iout=0
time=-1
  do i1=1,pts_s%n(2)
    iin=iin+1
    if(good(iin)>0.001) then
      iout=iout+1
      normal =45./atan(1.)*dip(iin)
      call shoot_rays_new(x_rays,z_rays,v_rays,vel2d,cmplx(x(iin),0.), &
        cmplx(normal,normal))
      dip(iout)=dip(iin)
      baz(iout)=baz(iin)
      good(iout)=good(iin)
      elev(iout)=elev(iin)
      event(iout)=event(iin)
      stat(iout)=stat(iin)
      x(iout)=x(iin)
      eout(iout)=err(1,iin)
      rays(1,:,iout)=x_rays(::jt,1)
      rays(2,:,iout)=0.
      rays(3,:,iout)=z_rays(::jt,1)
      temp(:,iout)=(/i1,1/)
      i=1; found=.false.
      do while (.not. found .and. i < size(z_rays,1))
        if(z_rays(i,1) > z) then
          found=.true.
        else
          i=i+1
        end if 
      end do
      if(found) then
        time(iout)=(z-z_rays(i-1,1))/(z_rays(i,1)-z_rays(i-1,1))+i-2
      end if
    end if 
end do
time=time*d_t/real(jt)
allocate(coords(2,iout))
coords=temp(:,:iout)
rays_s%n(2)=iout
call sep3d_set_sep3d(rays_s)
call sep3d_set_number_headers(rays_s,iout)
!if(.not. sep3d_set_coords(rays_s,coords))&
!  call seperr("trouble setting coords")
call sep3d_set_key_vals(rays_s,"P",dip)
call sep3d_set_key_vals(rays_s,"XSTA",x)
call sep3d_set_key_vals(rays_s,"ELEV",elev)
call sep3d_set_key_vals(rays_s,"EVENT_ID",event)
call sep3d_set_key_vals(rays_s,"STATION_ID",stat)
call sep3d_set_key_vals(rays_s,"BAZ",baz)
call sep3d_set_key_vals(rays_s,"err",eout)

call sep3d_set_key_vals(rays_s,"time",time)


allocate(ray2(size(rays,1),size(rays,2),iout))
ray2(:,:,:iout)=rays(:,:,:iout)
if(.not. sep3d_write_data("out",rays_s,ray2))&
     call seperr("trouble writing out data")
write(0,*) "CHECK WRITING",size(ray2,2),iout,iin

if(.not. sep3d_update_ntraces(rays_s)) &
  call seperr("trouble updating ntraces")
call sep3d_write_description("out",rays_s)

end program
