# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
module resid_shared_mod
  use tomo_space_conv_mod
  use put_3d_mod
  use err_shared_mod
  implicit none
  contains
  logical function stolt_errors(pos,extra,time,dip,semb,axis,scale,zero&
    &,dt,bulk)
    real,dimension(:) :: pos,axis,scale,zero,dt,dip,time
    real :: semb,slow,fact,alpha,convert,gamma,extra
    real, pointer :: vel(:,:,:)
    integer :: i
    logical :: bulk
    stolt_errors=.false.
! tan Apparent dip = cos theta (tan true)
!
! Given two apparent dips 90 degrees apart
!
! (tan true)^2 = (tan ap1)^2 + (tan ap2)^2
!
! For angle:
!
! tan ap1/tan true = cos theta1
!
! If both same use same sine +x, +y
! If not same use sine of +x, but rotate 180
!
!
! Error
! From Biondi's submitted paper
!
!
! deltaz = (1-p) sin^2 gamma
! ----------------------------------------------- z0
! 1-p (1- cos alpha) (cos^2 alpha - sin^2 gamma)
!
!
!
! where p - is the stolt scaling factor
! gamma - incidence angle
! alpha - dip (true???)
! z0 - depth at zero incidence
!
    if (.not. vel_tomo_3d_pt(vel)) then
      write(0,*) "trouble grabbing velocity field"
      return
    end if
    call init_tomo_coord()
    slow=1./get_val_lin_pos(vel,pos(3),pos(1),pos(2))
    alpha= atan( sqrt(dip(1)**2+dip(2)**2))
    convert=atan(1.)/45.
! fact = (1. -p) (z0 *slow0)
    fact = pos(3) * (semb-1.) *slow/cos(alpha)
    convert=1.
! maps(idip,itheta,igam)=
! (1-gamma)/(1-gamma*(1-cos(mydip)))*(sin(theta)**2)/
! (cos(theta))**2-sin(mydip)**2
!fact=pos(3)*(semb-1.)*slow
    fact=pos(3)*(1.-semb)*slow
!write(0,*) "CHECK SLOW",pos(3),semb,slow,axis(1:2)
    if (.not. bulk) then
      do i=1,size(axis)
        gamma=convert*axis(i)
        dt(i)=fact*sin(gamma)**2/(cos(alpha)*(cos(alpha)**2-sin(gamma)&
          &**2))
        zero(i)=-2.*cos(axis(i))
      end do
    else
      do i=1,size(axis)
        zero(i)=-2.*cos(alpha)
      end do
      dt=semb**2-1. !MIGHT (semb-1.)
    end if
    stolt_errors=.true.
  end function
end module
