program tomo_prog
  use traceback_mod
  use sep3d_struct_mod
  use sep
  use verb_mod
  implicit none
  type(sep3d)          :: ray_s,slow_w,junk_s
  logical              :: tau
  integer              :: n(3),ngather,i1,i2,i3
  integer              :: nh,block,nb,iray,n123
  real                 :: o(3),d(3),p2a,fact,alpha,sc,tot1,tot2,tot,tot4,orient
  real,allocatable     :: rays(:,:,:),junk(:)
  real,allocatable     :: baz(:),err(:)
  integer, allocatable :: iloc(:,:,:),n1(:),ipos(:)
  real, allocatable    :: scale(:,:,:),fract(:,:,:),tr_scale(:),spos(:),var(:)
  integer :: i


  call sep_init()
  call init_3d()
  call init_sep3d("in",ray_s,"INPUT")
  tau=.false.
  p2a=atan(1.)/45.
  call from_param("orient",orient,0.)

  if(.not. init_verbose()) call seperr("trouble initializing verbose")

  if(.not. init_traceback_func("slow","sigma_x","sigma_y",ray_s%d(1))) &
    call seperr("trouble initializing traceback function")

  ngather=1
  allocate(rays(3,ray_s%n(1)/3,ngather))
  allocate(iloc(3,ray_s%n(1),ngather),scale(3,ray_s%n(1),ngather))
  allocate(n1(ngather),fract(3,ray_s%n(1),ngather),tr_scale(ngather))
  allocate(baz(ngather))
  allocate(err(ngather))
  allocate(var(ngather))
  allocate(ipos(ngather*8*ray_s%n(1)*ray_s%n(2)*ray_s%n(3)))
  allocate(spos(ngather*8*ray_s%n(1)*ray_s%n(2)*ray_s%n(3)))


n123=product(ray_s%n(1:sep3d_ndims(ray_s)))/ray_s%n(1)/ray_s%n(2)
 if(.not. init_print_progress(0,n123)) &
   call seperr("trouble initialzing print progress")


  tot2=0.
  tot1=0.
i=0
  do i3=1,ray_s%n(3)
   do i2=1,ray_s%n(2)
    call print_progress("rays->operator")
    call sep3d_grab_headers("in",ray_s,nh,nwind=(/1,1/),fwind=(/i2-1,i3-1/))
    if(nh>0) then 
      call sep3d_grab_key_vals(ray_s,"err",err)
      call sep3d_grab_key_vals(ray_s,"baz",baz);baz=(baz-orient)*p2a
      call sep3d_grab_key_vals(ray_s,"var",var)
      if(.not. sep3d_read_data('in',ray_s,rays)) &
        call seperr("trouble reading in data")
      do iray=1,nh
        if(.not. trace_back_depth(rays(:,:,iray),1.,scale(:,:,iray),&
          iloc(:,:,iray),fract(:,:,iray),n1(iray)))&
          call seperr("trouble tracing back depth")
        scale(:,:,iray)=scale(:,:,iray)/abs(cos(baz(iray)))
        i=i+1
        tr_scale(iray)=err(iray)
      end do
      call add_traceback(scale,iloc,fract,nh,n1,tr_scale,ipos,spos,var,.true.)
    end if
   end do
  end do
  call finish_tracebacks()
end program
