# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
module calc_ray_mod
  use sep3d_struct_mod
  implicit none
  contains
  logical function get_ray_parameter(ray_in,ray_out,vel,vel_s)
    real :: vel(:,:)
    type(sep3d) :: vel_s
    real :: ray_in(:,:)
    real :: ray_out(:,:)
    real :: delta(4,size(ray_in,2))
    logical :: junk
    integer :: i
    junk= calc_delta(ray_in,delta)
    junk= get_vel(vel,vel_s,ray_in,ray_out(1,:))
    ray_out(2,:)=delta(1,:)/delta(4,:)
    do i =1,size(ray_out,2)
      if (delta(3,i)<.000) then
        ray_out(2,i)=-ray_out(2,i)
      end if
    end do
    ray_out(3,:)=ray_out(2,:)/ray_out(1,:)
    get_ray_parameter=.true.
  end function
  logical function get_vel(vel,vel_s,ray_in,ray_out)
    real :: vel(:,:)
    type(sep3d) :: vel_s
    real :: ray_in(:,:),ray_out(:)
    integer :: i
    integer :: ix,iz
    real :: fx,fz
    do i=1,size(ray_out)
      fx=(ray_in(1,i)-vel_s%o(2))/vel_s%d(2)+1.
      fz=(ray_in(3,i)-vel_s%o(1))/vel_s%d(1)+1.
      ix=fx
      iz=fz
      fx= fx-ix
      fz= fz-iz
      if (ix < 1) then
        ix=1
        fx=0.
      end if
      if (iz < 1) then
        iz=1
        fz=0.
      end if
      if (ix >= vel_s%n(2)) then
        ix=vel_s%n(2)-1
        fx=1.
      end if
      if (iz >= vel_s%n(1)) then
        iz=vel_s%n(1)-1
        fz=1.
      end if
      ray_out(i)=(1.-fx)*(1.-fz)*vel(iz,ix) +(fx)*(1.-fz)*vel(iz,ix+1&
        &) +(1.-fx)*(fz)*vel(iz,ix) +(fx)*(fz)*vel(iz+1,ix+1)
    end do
    get_vel=.true.
  end function
  integer function tb_length(ray)
    real :: ray(:,:)
    real :: last
    if (any(ray(:,size(ray,2)-1) .ne. ray(:,size(ray,2)))) then
      tb_length=size(ray,2)
    else
      tb_length=1
      do while (any(ray(:,tb_length).ne.ray(:,tb_length+1)))
        tb_length= tb_length+1
      end do
    end if
  end function
  logical function convert_3d(ray_in,ray_out)
    complex :: ray_in(:)
    real :: ray_out(:,:)
    integer :: i
    ray_out(1,:)=real(ray_in)
    ray_out(3,:)=aimag(ray_in)
    ray_out(2,:)=0
    do i=1,size(ray_out,2)
      if (ray_out(3,i) < 0.) then
        ray_out(:,i)=0.
      end if
    end do
    convert_3d=.true.
  end function
  logical function calc_center(ray_in,center)
    real ray_in(:,:),center(:,:)
    calc_center=.false.
    center(:,1:size(ray_in)-1)=0.5*(ray_in(:,2:size(ray_in)-1)+ray_in(:&
      &,1:size(ray_in)-1))
    calc_center=.true.
  end function
  logical function calc_delta(ray_in,delta)
    real ray_in(:,:),delta(:,:)
    calc_delta=.false.
    delta(1:3,1:size(ray_in,2)-1)=(ray_in(:,2:size(ray_in,2)-1)-ray_in&
      &(:,1:size(ray_in,2)-1))
    delta(4,:)=sqrt(sum(delta(1:3,:)**2,1))
    delta(:,size(ray_in,2))=delta(:,size(ray_in,2)-1)
    calc_delta=.true.
  end function
end module
