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
  integer              :: nh,f,block,nb,iray
  real                 :: o(3),d(3),p2a,fact,alpha,sc,tot1,tot2,tot,tot4
  real,allocatable     :: rays(:,:,:),junk(:)
  real,allocatable     :: gamma(:),theta(:),posz(:),slow(:),dipx(:),dipy(:),posx(:)
  integer, allocatable :: iloc(:,:,:),n1(:),ipos(:)
  real, allocatable    :: scale(:,:,:),fract(:,:,:),tr_scale(:),spos(:),var(:)
  real, allocatable    :: scalez(:,:),scalea(:,:),fractz(:,:),fracta(:,:)
  integer,allocatable  :: ilocz(:,:),iloca(:,:)
  integer              :: nz,na,nout


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

  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")

 if(abs(ray_s%o(2)) >.0001) call seperr("must include zero")

  ngather=ray_s%n(2)*ray_s%n(3)
  allocate(rays(3,ray_s%n(1)/3,ngather))
  allocate(gamma(ngather),theta(ngather))

  allocate(iloca(3,ray_s%n(1)*2),scalea(3,ray_s%n(1)*2))
  allocate(fracta(3,ray_s%n(1)*2))

  allocate(ilocz(3,ray_s%n(1)*2),scalez(3,ray_s%n(1)*2))
  allocate(fractz(3,ray_s%n(1)*2))

  allocate(iloc(3,ray_s%n(1)*2,ngather),scale(3,ray_s%n(1)*2,ngather))
  allocate(fract(3,ray_s%n(1)*2,ngather))

  allocate(n1(ngather),tr_scale(ngather))
  allocate(slow(ngather),posx(ngather),posz(ngather),dipx(ngather),dipy(ngather),var(ngather))
  allocate(ipos(ngather*16*ray_s%n(1)/3))
  allocate(spos(ngather*16*ray_s%n(1)/3))

 if(.not. init_print_progress(0,ray_s%n(4))) &
   call seperr("trouble initialzing print progress")


  tot2=0.
  tot1=0.
  nout=0
  do f=1,ray_s%n(4)
    call print_progress("rays->operator")
    call sep3d_grab_headers("in",ray_s,nh,nwind=(/ray_s%n(2),ray_s%n(3),1/),&
      fwind=(/0,0,f-1/))
    
    if(nh>0) then
      call sep3d_grab_key_vals(ray_s,"gamma",gamma)
      call sep3d_grab_key_vals(ray_s,"vel",slow);slow=1./slow
      call sep3d_grab_key_vals(ray_s,"pos_x",posx)
      call sep3d_grab_key_vals(ray_s,"pos_z",posz)
      call sep3d_grab_key_vals(ray_s,"theta",theta);theta=theta*p2a
      call sep3d_grab_key_vals(ray_s,"dip_x",dipx)
      call sep3d_grab_key_vals(ray_s,"dip_y",dipy)
      call sep3d_grab_key_vals(ray_s,"var",var)
      fact=posz(1)*(1.-gamma(1))*slow(1)
      if(.not. sep3d_read_data('in',ray_s,rays)) &
        call seperr("trouble reading in data")
      alpha=atan(dipx(1))
      if(abs(theta(1)) < .001) then !if we have the zero angle
        if(tau) then
          if(.not. trace_back_tau(rays(:,:,1),1.,scalez,&
            ilocz,fractz,nz))&
            call seperr("trouble tracing back tau")
        else
          if(.not. trace_back_depth(rays(:,:,1),1.,scalez,&
            ilocz,fractz,nz))&
            call seperr("trouble tracing back depth")
        end if
       
  nout=nout+nh-1
        do iray=2,nh
           call srite("pos",posz,4)
           call srite("pos",posx,4)
           call srite("pos",theta(iray)/p2a,4)
           call srite("pos",gamma,4)
          !for now the 1-D dip
          if(tau) then
            if(.not. trace_back_tau(rays(:,:,iray),1.,scalea,&
              iloca,fracta,na))&
              call seperr("trouble tracing back tau")
          else
            if(.not. trace_back_depth(rays(:,:,iray),1.,scalea,&
              iloca,fracta,na))&
              call seperr("trouble tracing back depth")
          end if
          scale(:,:na,iray-1)         =scalea(:,:na)
          iloc(:,:na,iray-1)          =iloca(:,:na)
          fract(:,:na,iray-1)         =fracta(:,:na)
          scale(:,na+1:na+nz,iray-1)  =-cos(theta(iray))*scalez(:,:nz)
       !   scale(:,:na,iray-1)         =cos(theta(iray))*scalea(:,:na)
        !  scale(:,na+1:na+nz,iray-1)  =scalez(:,:nz)
          
          
          
          iloc(:,na+1:na+nz,iray-1)   =ilocz(:,:nz)
          fract(:,na+1:na+nz,iray-1)  =fractz(:,:nz)
          n1(iray-1)=nz+na
          !n1(iray-1)=na

          tr_scale(iray-1)=fact*sin(theta(iray))**2/&
            (cos(alpha)*(cos(alpha)**2-sin(theta(iray))**2))
          tr_scale(iray-1)=fact*sin(theta(iray))**2/&
            (cos(alpha)*(cos(alpha)**2-sin(theta(iray))**2))
          tr_scale(iray-1)=-fact*tan(theta(iray))**2
     
        end do
        call add_traceback(scale,iloc,fract,nh-1,n1,tr_scale,ipos,spos,var,.false.)
      end if
    end if
  end do
  call to_history("n1",4,"pos")
  call to_history("n2",nout,"pos")
  call finish_tracebacks()
end program
