# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
!
!
! 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
  implicit none
  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"
      return
    end if
    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 function
  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
end module
