# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
 !
! Sun Jan 12 21:57:03 PST 2003
!
!
!
module back_project_mod
{
use basic_util_mod
logical, private, save :: verb
real, private, save :: beg1(50),end1(50),beg2(50),end2(50),beg3(50),end3(50)
real, private, save :: bege(50),ende(50),begt(50),endt(50)
integer, private, save :: nreg
contains
logical function back_project_standard(ang_tag,zero_tag){
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(:),zero(:)
real,allocatable :: s_x(:),s_y(:),s_z(:),time_err(:),extra(:),theta(:),scale(:)
integer :: ncz,nca,fwind(2),nwind(2),cang,czero
integer :: nbig,izero,bca,nerr_tot
logical :: reade1,reade2
integer, external :: sreed
character(len=128) :: err_method
back_project_standard=.false.
if (.not. vel_structure(vel_s))
{
  write(0,*) "trouble getting vel structure"
  return
}
call init_sep3d(ang_tag,ang_s,"INPUT")
call init_sep3d(zero_tag,zero_s,"INPUT")
call init_sep3d("err",err_s,"INPUT")
call set_n_o_d_3(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:12).eq."vertical_pts" .or. err_method(1:10).eq."stolt_refs".or. err_method(1:13).eq."vertical_refs")
{
}
else
{
  write(0,*) "back project standard doesn't understand error method",trim(err_method)
  return
}
if (ns(3).eq.1)
{
  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)))
}
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),theta(nbig),scale(nbig),zero(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){
  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=nwind)
  call sep3d_grab_grid_block(ang_tag,grid)
  if (czero > 0 .and. .and. cang >0)
{
    if (ns(3).eq.1)
{
      reade1=sep3d_read_data(zero_tag,zero_s,data2)
      dataz(1,:,nca)=real(data2(:,1))
      dataz(2,:,nca)=aimag(data2(:,1))
      reade2=sep3d_read_data(ang_tag,ang_s,data2)
      dataa(1,:,nca+1:nc+cang)=real(data2(:,1:cang))
      dataa(2,:,nca+1:nc+cang)=aimag(data2(:,1:cang))
      nca+=cang
}
    else
{
      reade1=sep3d_read_data(zero_tag,zero_s,data3)
      dataz(:,:,nca+1)=data3(:,:,1)
     reade2=sep3d_read_data(ang_tag,ang_s,data3)
     dataa(:,:,nca+cang)=data3(:,:,1:cang)
     nca+=cang
}
   bca=nca-cang+1
   if (.not. reade1 .or. .not. reade2)
{
     write(0,*) "trouble reading in rays at ",iloc
     return
}
   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,"theta",theta(bca:nca))
   call sep3d_grab_key_vals(ang_s,"description",extra(bca:nca))
   if (.not. parse_valid_back_pts(datac,bca,nca,dataz,s_x,s_y,s_z,theta,extra))
{
     write(0,*) "trouble parsing back project region"
     return
}
   if (nca >= bca)
{
      do i1=bca,nca
{
         dataz(:,:,:)=dataz(:,:,bca)
}
      if (.not. get_data_error(s_x(bca:nca),s_y(bca:nca),s_z(bca:nca),theta(bca:nca),extra(bca:nca),time_err(bca:nca),scale(bca:nca) zeros(bca:nca)))
{
           write(0,*) "trouble calculating data error"
           return
}
}
}
 if (iloc.eq.ang_s%n(3) .or. nca >= 1000)
{
ERROR BEFORE ERROR
               if (.not. trace_back(dataa(:,:,1:nca), dataz(:,:,1:nca),scale(1:nca),zeros(1:nca),time_err(1:nca)){ ERROR AFTER ERROR
       write(0,*) "trouble tracing back rays"
       return
}
    nca=0
}
}
if (ns(3).eq.1) deallocate(data2)
else deallocate(data3)
deallocate(grid,err,dataa,dataz,zero)
deallocate(s_x,s_y,s_z,extra,theta,scale)
call sep3df_delete(ang_s)
call sep3df_delete(zero_s)
call sep3df_delete(err_s)
back_project_standard=.true.
}
logical function parse_valid_back_pts(datac,bca,nca,dataz,s_x,s_y,s_z,theta,extra){
real :: datac(:,:,:),dataz(:,:,:)
integer :: nca,bca
logical :: kill
integer :: i,ncn,ikill
real :: s_z(:),s_x(:),s_y(:),theta(:),extra(:)
parse_valid_back_pts=.false.
if (nreg.eq.0){
  parse_valid_back_pts=.true.
  return
}
ncn=bca
do i=bca,nca
{
  kill=.false.
  do ikill=1,nreg
{
    if (
     s_x(i) >= beg2(ikill) .and. s_x(i) <= end2(ikill) .and.
     s_y(i) >= beg3(ikill) .and. s_y(i) <= end3(ikill) .and.
     s_z(i) >= beg1(ikill) .and. s_z(i) <= end1(ikill) .and.
     theta(i) >= begt(ikill) .and. theta(i) <= endt(ikill) .and.
     extra(i) >= bege(ikill) .and. extra(i) <= ende(ikill) ) kill=.true.
}
  if (.not. kill .and. ncn .ne. i){
    datac(:,:,ncn)=datac(:,:,i)
    s_x(ncn)=s_x(i)
    s_y(ncn)=s_y(i)
    s_z(ncn)=s_z(i)
    theta(ncn)=theta(i)
    extra(ncn)=extra(i)
    ncn+=1
}
}
ncn-=1
if (ncn < bca) nca=bca-1
parse_valid_back_pts=.true.
}
logical check_pars_back_project(file){
character(len=*) :: file
logical :: stat(20)
integer :: i
check_pars_back_project=.false.
beg1=default_real
beg2=default_real
beg3=default_real
end1=default_real
end2=default_real
end3=default_real
begt=default_real
endt=default_real
bege=default_real
ende=default_real
stat(1)= check_par_file("kill_back_beg1",beg1,tag,beg1)
stat(2)= check_par_file("kill_back_beg2",beg2,tag,beg2)
stat(3)= check_par_file("kill_back_beg3",beg3,tag,beg3)
stat(4)= check_par_file("kill_back_end1",end1,tag,end1)
stat(5)= check_par_file("kill_back_end2",end2,tag,end2)
stat(6)= check_par_file("kill_back_end3",end3,tag,end3)
stat(7)= check_par_file("kill_back_beg_extra",bege,tag,bege)
stat(8)= check_par_file("kill_back_end_extra",ende,tag,ende)
stat(9)= check_par_file("kill_back_beg_theta",begt,tag,begt)
stat(10)= check_par_file("kill_back_end_theta",endt,tag,endt)
nreg=1
do while (beg1(nreg).ne.default_real .or. beg2(nreg).ne.default_real .or. bege(nreg).ne.default_real .or. ende(nreg).ne.default_real .or. begt(nreg).ne.default_real .or. endt(nreg).ne.default_real .or. beg3(nreg).ne.default_real .or. end1(nreg).ne.default_real .or.
   end2(nreg).ne.default_real .or. end3(nreg).ne.default_real) nreg+=1
where(beg1.eq.default_real) beg1=-1000000.
where(beg2.eq.default_real) beg2=-1000000.
where(beg3.eq.default_real) beg3=-1000000.
where(end1.eq.default_real) end1=1000000.
where(end2.eq.default_real) end2=1000000.
where(end3.eq.default_real) end3=1000000.
where(bege.eq.default_real) bege=-1000000.
where(ende.eq.default_real) ende=1000000.
where(begt.eq.default_real) begt=-1000000.
where(endt.eq.default_real) endt=1000000.
stat=.true.
call from_param("verb",verb,.false.)
if (all(stat)) check_pars_back_project=.true.
}
}
