# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
!
!
! Thu Jan 9 21:34:42 PST 2003
!
!
module tomo_full_mod
  use sep3d_struct_mod
  use param_check_mod
  use back_project_mod
  use ray_database_mod
  use traceback_mod
  use error_mod
  implicit none
  logical, private, save :: verb,tau
  character(len=128), private, save :: ray0_tag,ang_tag,err_method
  real, private, save :: oang,dang
  integer, private, save :: nang
  contains
  logical function calc_full_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
    calc_full_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)
    loc(1,:)=x
    loc(2,:)=y
    loc(3,:)=z
    dips(1,:)=dip1
    dips(2,:)=dip2
    descriptor(:)=descriptor
    if (verb) then
      write(0,*) "    Shooting  angle rays"
    end if
    if (.not. create_ray_database(ang_tag,oang,dang,nang,.true.,loc,dips&
      &,descriptor)) then
      write(0,*) "trouble creating ray database for tag ",trim(ang_tag&
        &)
      return
    end if
    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(ang_tag)) then
      write(0,*) "trouble creatting error map"
      return
    end if
    if (.not. new_traceback(1,iloc,fract,scale,nlen,data)) then
      write(0,*) "trouble create new back project"
      return
    end if
    if (verb) then
      write(0,*) "    Converting rays to operator"
    end if
    if (.not. back_project_full(ang_tag,ray0_tag,err_method,trace_back_full&
      &,.false.)) then
      write(0,*) "trouble back projecting ray"
      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("data",trim(data),"tomo_oper")
    call to_history("nlen",trim(nlen),"tomo_oper")
    call to_history("traceback_type","full","tomo_oper")
    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_full_tomo=.true.
  end function
  logical function load_full_tomo_oper(dat)
    real, pointer :: dat(:)
    load_full_tomo_oper=.false.
    if (.not. load_calc_traceback(1,"full",dat)) then
      write(0,*) "trouble loading traceback 0"
      return
    end if
    load_full_tomo_oper=.true.
  end function
  integer function tomo_full_op(adj,add,model,data)
    logical, intent (in) :: adj,add
    real :: model(:),data(:)
    call set_current(1)
    call adjnull(adj,add,model,data)
    tomo_full_op=0
    call tb_standard_op(adj,model,data)
    tomo_full_op=0
  end function
  logical function trace_back_full(nt,ang_ray,zero_ray,scale,zeros,err&
    &)
    real :: ang_ray(:,:,:),zero_ray(:,:,:)
    real :: scale(:),zeros(:),err(:)
    integer :: n2,i2,n1,j,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
    trace_back_full=.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
      ray_type=tau_type
    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. trace_back(zero_ray(:,:,i2),ray_type,zeros(i2),sc(:,n1a&
        &+1:),iloc(:,n1a+1:),fract(:,n1a+1:),n1z)) then
        write(0,*) "trouble tracing back zero ray"
        return
      end if
! sc(1,n1a+1:n1a+n1z)=0.
! sc(1,1:n1a)=0.
      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_full=.true.
  end function
  logical function trace_back_bulk(nt,ang_ray,zero_ray,scale,zeros,err&
    &)
    real :: ang_ray(:,:,:),zero_ray(:,:,:)
    real :: scale(:),zeros(:),err(:)
    integer :: n2,i2,n1,j,n1z,n1a,nt
    real :: fract(3,nt*2)
    integer :: iloc(3,nt*2)
    real :: sc(3,nt*2),a
    integer :: ray_type
    character(len=128) :: floc,ffract,fsc,flen,ferr
    trace_back_bulk=.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
      ray_type=tau_type
    else
      ray_type=z_type
    end if
    do i2=1,n2
      if (.not. trace_back_b(ang_ray(:,:,i2),ray_type,scale(i2),sc,iloc&
        &,fract,n1a,err(i2))) then
        write(0,*) "trouble tracing back ang ray"
        return
      end if
      if (.not. trace_back_b(zero_ray(:,:,i2),ray_type,zeros(i2),sc(:&
        &,n1a+1:),iloc(:,n1a+1:),fract(:,n1a+1:),n1z,a)) then
        write(0,*) "trouble tracing back zero 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_bulk=.true.
  end function
  logical function back_project_full(ang_tag,zero_tag,err_method,trace_b&
    &,bulk)
    interface
      logical function trace_b(nt,ang_ray,zero_ray,scale,zeros,err)
        real :: ang_ray(:,:,:),zero_ray(:,:,:)
        real :: scale(:),zeros(:),err(:)
        integer :: nt
      end function
    end interface
    character(len=*) :: ang_tag,zero_tag,err_method
    type(sep3d) :: ang_s,zero_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(:),x_loc(:),z_loc(:)
    real, pointer :: junk(:,:)
    integer :: ncz,nca,fwind(2),nwind(2),cang,czero
    integer :: nbig,izero,bca,nerr_tot,i1,nc,i2
    logical :: reade1,reade2,bulk
    integer, external :: sreed
    back_project_full=.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(zero_tag,zero_s,"INPUT")
    call init_sep3d("err",err_s,"INPUT")
    call sep3d_grab_headers("err",err_s,cang)
    allocate(err(err_s%n(2)),dip_xz(err_s%n(2)),dip_yz(err_s%n(2)))
    allocate(x_loc(err_s%n(2)))
    allocate(z_loc(err_s%n(2)))
    call sep3d_grab_key_vals(err_s,"s_x",x_loc)
    call sep3d_grab_key_vals(err_s,"s_z",z_loc)
    call sep3d_grab_key_vals(err_s,"dip_xz",dip_xz)
    call sep3d_grab_key_vals(err_s,"dip_yz",dip_yz)
    allocate(junk(vel_s%n(1),vel_s%n(2)))
    junk=1.
    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)
    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 init_coord(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:12).eq."vertical_pts"&
      & .or. err_method(1:4).eq."bulk" .or. err_method(1:10&
      &).eq."stolt_refs".or. err_method(1:13).eq."vertical_refs") then
    else
      write(0,*) "back project full doesn't understand error method="&
        &,trim(err_method)
      return
    end if
    if (ns(3).eq.1) then
      nt=ang_s%n(1)
      allocate(data2(ang_s%n(1),ang_s%n(2)))
    else
      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/2,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
      i2=(x_loc(iloc)-vel_s%o(2))/vel_s%d(2)+1.5
      i1=(z_loc(iloc)-vel_s%o(1))/vel_s%d(1)+1.5
      junk(i1,i2)=err(iloc)
      fwind(2)=iloc-1
      call sep3d_grab_headers(ang_tag,ang_s,cang,fwind=fwind,nwind=nwind&
        &)
      call sep3d_grab_headers(zero_tag,zero_s,czero,fwind=fwind,nwind=&
        &(/1,1/))
      call sep3d_grab_grid_block(ang_tag,grid)
      if (czero > 0 .and. cang >0) then
        if (ns(3).eq.1) then
          reade1=sep3d_read_data(zero_tag,zero_s,data2(:nt/2,:))
          dataz(1,:,nca+1)=real(data2(:,1))
          dataz(3,:,nca+1)=aimag(data2(:,1))
          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
          reade1=sep3d_read_data(zero_tag,zero_s,data3(:size(data3,1)/2&
            &,:))
          dataz(:,:,nca+1)=reshape(data3(:,1),(/3,size(dataz,2)/))
          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
            dataz(:,:,i1)=dataz(:,:,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),bulk)) 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_b(size(Dataa,2),dataa(:,:,1:nca), dataz(:,:,1:nca&
          &),scale(1:nca),zeros(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)
    call sep3df_delete(zero_s)
    back_project_full=.true.
  end function
  logical function check_par_tomo_full(struct,tag,ok)
    type(sep3d) :: struct
    character(len=*) :: tag
    logical :: stat_ok(20),stat_rerun(20),ok
    check_par_tomo_full=.false.
    stat_ok=.true.
    stat_rerun=.true.
    if (model_space(1:2).ne."ds") then
      write(0,*) "bulk tomo mode requires model_space=ds"
      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
    call from_param("verb",verb,.false.)
    call from_param("tau",tau,.false.)
    ray0_tag="ray0_tag"
    ang_tag="ang_tag"
    stat_rerun(1)=check_par_file("ray0_tag",ray0_tag,tag,ray0_tag)
    stat_rerun(2)=check_par_file("ang_tag",ang_tag,tag,ang_tag)
    stat_rerun(3)=check_par_file("oang",oang,tag,2.)
    stat_rerun(4)=check_par_file("dang",dang,tag,2.)
    stat_rerun(5)=check_par_file("nang",nang,tag,20)
    stat_rerun(7)=check_par_file("data_err_method",err_method,tag&
      &,"vertical_pts")
    stat_rerun(8)=data_errors_build_ok
    stat_ok(9)=check_pars_back_project(tag,stat_rerun(9))
    if (all(stat_ok)) then
      check_par_tomo_full=.true.
    end if
    if (all(stat_rerun)) then
      ok=.true.
    else
      ok=.false.
    end if
  end function
end module
