!
!
! Thu Jan  9 21:34:42 PST 2003
!
!
module tomo_normal_mod 
  use sep3d_struct_mod
  use back_project_mod
  use traceback_mod
  use param_check_mod
  use error_mod
  use handle_rays_3d_mod
  use handle_rays_2d_mod
  use ray_database_mod
  logical, private, save :: verb,tau
  integer, private, save :: nt
  character(len=128), private, save :: ray0_tag
  contains
  logical function calc_normal_tomo()
    type(sep3d)    :: err_s
    integer        :: npts
    real, pointer  :: x(:),y(:),z(:),dip1(:),dip2(:),loc(:,:),dips(:,:),descriptor(:)
    character(len=128) :: iloc,fract,scale,nlen,data
    logical :: tau
    calc_normal_tomo=.false.
    call init_sep3d("err",err_s,"INPUT")
    call sep3d_grab_headers("err",err_s,npts)
    allocate(x(npts),y(npts),z(npts),dip1(npts),dip2(npts))
    allocate(loc(3,npts),dips(2,npts),descriptor(npts))
    call sep3d_grab_key_vals(err_s,"s_x",x)
    call sep3d_grab_key_vals(err_s,"s_y",y)
    call sep3d_grab_key_vals(err_s,"s_z",z)
    call sep3d_grab_key_vals(err_s,"dip_xz",dip1)
    call sep3d_grab_key_vals(err_s,"dip_yz",dip2)
    call sep3d_grab_key_vals(err_s,"descriptor",descriptor)
    call sep3df_delete(err_s,"err")
    loc(1,:)=x
    loc(2,:)=y
    loc(3,:)=z
    dips(1,:)=dip1
    dips(2,:)=dip2
    descriptor(:)=descriptor
    if (verb) then
      write(0,*) "    Shooting  zero angle rays"
    end if
    if (.not. create_ray_database(ray0_tag,0.,1.,1,.false.,loc,dips,descriptor)) then
      write(0,*) "trouble creating ray database for tag ",trim(ray0_tag)
      return
    end if
    if (.not. init_err_map(ray0_tag)) then
      write(0,*) "trouble initializing error mapping"
    end if
    return 
  end function 
  if (.not. new_traceback(1,iloc,fract,scale,nlen,data)) then
    write(0,*) "trouble create new back project"
    return
  end if
  call to_history("iloc",trim(iloc),"tomo_oper")
  call to_history("fract",trim(fract),"tomo_oper")
  call to_history("scale",trim(scale),"tomo_oper")
  call to_history("nlen",trim(nlen),"tomo_oper")
  call to_history("data",trim(data),"tomo_oper")
  call to_history("traceback_type","normal","tomo_oper")
  if (verb) then
    write(0,*) "    Back projecting normal rays"
  end if
  if (.not. back_project_normal("err",ray0_tag)) then
    write(0,*) "trouble creating standard back projection"
    return
  end if
  if (.not. finish_cur_trace_back()) then
    write(0,*) "trouble finishing back project"
    return
  end if
  deallocate(x,y,z,dip1,dip2,loc,dips)
  calc_normal_tomo=.true.
end module
logical function load_normal_tomo_oper(dat)
  real, pointer :: dat(:)
  load_normal_tomo_oper=.false.
  if (.not. load_calc_traceback(1,"normal",dat)) then
    write(0,*) "trouble loading traceback 0"
    return
  end if
  load_normal_tomo_oper=.true.
end function 
integer function tomo_normal_op(adj,add,model,data)
  logical, intent (in) :: adj,add
  real :: model(:),data(:)
  call set_current(1)
  call adjnull(adj,add,model,data)
  tomo_normal_op=0
  call tb_standard_op(adj,model,data)
end function 
logical function trace_back_normal(nt,ang_ray,scale,err)
  real :: ang_ray(:,:,:)
  real :: scale(:),err(:)
  integer :: n2,i2,n1,n1z,n1a,nt
  real  ::   fract(3,nt*2)
  integer :: iloc(3,nt*2)
  real    :: sc(3,nt*2)
  integer :: ray_type
  character(len=128) :: floc,ffract,fsc,flen,ferr,fdat
  trace_back_normal=.false.
  if (.not. create_tb_names(1,floc,ffract,fsc,flen,ferr)) then
    write(0,*) "trouble creating traceback names"
    return
  end if
  n2=size(ang_ray,3)
  if (tau) then
    write(0,*) "Don't know how to do normal ray tomography in tau"
    return
  else
    ray_type=z_type
  end if 
  do i2=1,n2 
    if (.not. trace_back(ang_ray(:,:,i2),ray_type,scale(i2),sc,iloc,fract,n1a)) then
      write(0,*) "trouble tracing back ang ray"
      return
    end if
    if (.not. write_trace_back(floc,ffract,fsc,flen,ferr,iloc,fract,sc,n1a+n1z,err(i2))) then
      write(0,*) "trouble writing trace back"
      return
    end if
  end do 
  trace_back_normal=.true.
end function 
logical function back_project_normal(err_tag,ang_tag)
  character(len=*) :: err_tag,ang_tag
  type(sep3d)      :: ang_s,err_s,vel_s
  integer          :: iloc,iang,ns(3),ndim,nt
  real             :: os(3),ds(3)
  real,allocatable :: dataa(:,:,:),dataz(:,:,:)
  complex,allocatable :: data2(:,:)
  real,   allocatable :: data3(:,:),err(:)
  integer, allocatable :: grid(:)
  real,allocatable :: fact(:),zeros(:),dip_xz(:),dip_yz(:)
  real,allocatable :: s_x(:),s_y(:),s_z(:),time_err(:),extra(:),theta(:),scale(:),rayt(:)
  integer          :: ncz,nca,fwind(2),nwind(2),cang,czero
  integer          :: nbig,izero,bca,nerr_tot,i1
  logical          :: reade1,reade2
  integer, external :: sreed
  character(len=128) :: err_method
  back_project_normal=.false.
  if (.not. vel_structure(vel_s)) then
    write(0,*) "trouble getting vel structure"
    return
  end if
  call init_sep3d(ang_tag,ang_s,"INPUT")
  call init_sep3d(err_tag,err_s,"INPUT")
  allocate(err(err_s%n(2)),dip_xz(err_s%n(2)),dip_yz(err_s%n(2)))
  call sep3d_grab_key_vals(err_s,"dip_xz",dip_xz)
  call sep3d_grab_key_vals(err_s,"dip_yz",dip_yz)
  if (.not. sreed_field("err",err)) then
    write(0,*) "trouble reading in error file"
    return
  end if
  call sep3d_grab_headers("err",err_s,cang,fwind=fwind,nwind=nwind)
  call sep3df_delete(err_s,"err")
  if (size(err).ne.ang_s%n(3)) then
    write(0,*) "number of errors not the same as number of errors"
    return
  end if
  call set_n_o_d_3d(vel_s,ns,os,ds)
!call from_aux("aux","data_err_method",err_method,"undefined")
  if (err_method(1:9).eq."stolt_pts" .or.  err_method(1:10).eq."stolt_refs") then
  else
    write(0,*) "back project normal doesn't understand error method",trim(err_method)
    return
  end if 
  if (ns(3).eq.1) then
    ndim=2
    nt=ang_s%n(1)
    allocate(data2(ang_s%n(1),ang_s%n(2)))
  else
    ndim=3
    nt=ang_s%n(1)/3
    allocate(data3(ang_s%n(1),ang_s%n(2)))
  end if 
  nbig=1000*ang_s%n(2)
  allocate(dataa(3,nt,nbig))
  allocate(dataz(3,nt,nbig))
  allocate(grid(ang_s%n(2)))
  allocate(s_x(nbig),s_y(nbig),s_z(nbig),extra(nbig),rayt(nbig),theta(nbig),scale(nbig),zeros(nbig),time_err(nbig))
  fwind=0
  nwind=(/ang_s%n(2),1/)
  dataa(2,:,:)=os(3)+ds(3)/2.
  dataz(2,:,:)=os(3)+ds(3)/2.
  scale=1.
  nca=0
  nerr_tot=0
  do iloc=1,ang_s%n(3)
    if (verb .and. 1.eq.mod(iloc, nint(ang_s%n(3)/10.))) then
      write(0,*) "      working on ",iloc," of ",ang_s%n(3)
    end if
    call sep3d_grab_headers(ang_tag,ang_s,cang,fwind=fwind,nwind=nwind)
    call sep3d_grab_grid_block(ang_tag,grid)
    if (cang >0) then
      if (ns(3).eq.1) then
        reade2=sep3d_read_data(ang_tag,ang_s,data2)
        dataa(1,:,nca+1:nca+cang)=real(data2(:,1:cang))
        dataa(3,:,nca+1:nca+cang)=aimag(data2(:,1:cang))
        nca=      nca+cang
      else
        reade2=sep3d_read_data(ang_tag,ang_s,data3)
        dataa(:,:,nca+1:nca+cang)=reshape(data3(:,1:cang),(/3,size(dataa,2),cang/))
        nca=      nca+cang
      end if 
      bca=nca-cang+1
      if (.not. reade1 .or. .not. reade2) then
        write(0,*) "trouble reading in rays at ",iloc
        return
      end if
      call sep3d_grab_key_vals(ang_s,"s_x",s_x(bca:nca))
      call sep3d_grab_key_vals(ang_s,"s_y",s_y(bca:nca))
      call sep3d_grab_key_vals(ang_s,"s_z",s_z(bca:nca))
      call sep3d_grab_key_vals(ang_s,"time",rayt(bca:nca))
      call sep3d_grab_key_vals(ang_s,"theta",theta(bca:nca))
      call sep3d_grab_key_vals(ang_s,"description",extra(bca:nca))
      if (.not. parse_valid_back_pts(dataa,bca,nca,dataz,s_x,s_y,s_z,theta,rayt,extra)) then
        write(0,*) "trouble parsing back project region"
        return
      end if
      if (nca >= bca) then
        do i1=bca,nca 
          dataa(:,:,i1)=dataa(:,:,bca)
        end do 
        if (.not. get_data_errors((/s_x(bca),s_y(bca),s_z(bca)/),extra(bca),rayt(bca:nca),(/dip_xz(iloc),dip_yz(iloc)/),err(iloc),theta(bca:nca),time_err(bca:nca),scale(bca:nca),zeros(bca:nca),.false.)) then
          write(0,*) "trouble calculating data error"
          return
        end if
      end if
    end if
    if (iloc.eq.ang_s%n(3) .or. nca >= 1000) then
      if (.not. trace_back_normal(size(Dataa,2),dataa(:,:,1:nca),scale(1:nca),time_err(1:nca))) then
        write(0,*) "trouble tracing back rays"
        return
      end if
      nca=0
    end if
  end do 
  if (ns(3).eq.1) then
    deallocate(data2)
  else
    deallocate(data3)
  end if 
  deallocate(grid,err,dataa,dataz,zeros)
  deallocate(dip_xz,dip_yz)
  deallocate(s_x,s_y,s_z,extra,theta,scale)
  call sep3df_delete(ang_s,ang_tag)
  call sep3df_delete(err_s,"err")
  back_project_normal=.true.
end function 
logical function check_par_tomo_normal(struct,tag,ok)
  type(sep3d)        :: struct
  character(len=*)   :: tag
  logical        :: stat_ok(20),stat_rerun(20),ok
  character(len=256)   :: data_err_method
  check_par_tomo_normal=.false.
  stat_ok=.true.
  stat_rerun=.true.
  if (model_space(1:2).ne."v2") then
    write(0,*) "bulk tomo mode requires model_space=v2"
    return
  end if
  call from_param('tau',tau,.false.)
  if (tau) then
    write(0,*) "normal ray tomography won't work in tau domain"
    return
  end if
  if (sep3d_ndims(struct) .eq.2) then
    stat_ok(1)=read_check_ray_params_2d(tag,stat_rerun(1))
  else
    stat_ok(1)=read_check_ray_params_3d(tag,stat_rerun(1))
  end if 
  stat_rerun(2)=check_par_file("ray0_tag",ray0_tag,tag,ray0_tag)
!stat_rerun(3)=check_par_file("data_err_method",data_err_method,tag)
  stat_rerun(4)=data_errors_build_ok
  stat_ok(5)=check_pars_back_project(tag,stat_rerun(5))
  if (all(stat_rerun)) then
    check_par_tomo_normal=.true.
  end if
  if (all(stat_ok)) then
    ok=.true.
  else
    ok=.false.
  end if 
end function 
ERROR BEFORE ERROR
}
ERROR AFTER