# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
!this is from Paul's hwt3d code
module shoot_types
  use double_mod
    implicit none
  type point3d
    real :: x,y,z,v,t,l
  end type
  type v3d
    real(kind=wp) :: dx,dy,dz
  end type
end module
module shoot3d_mod
  use shoot_types
  implicit none
  real, pointer,private,save :: vel(:,:,:)
  integer,private,save :: nv(3)
  real,private,save :: ov(3),dv(3),dt,min1,min2,max1,max2,min3,max3
  contains
  subroutine init_shot3d(velocity,n_in,o_in,d_in,d1)
    real, pointer :: velocity(:,:,:)
    real :: o_in(3),d_in(3),d1
    integer :: n_in(3)
    nv=n_in
    dv=d_in
    ov=o_in
    vel=>velocity
    min1=ov(1)
    min2=ov(2)
    min3=ov(3)
    max1=ov(1)+dv(1)*(nv(1)-1)
    max2=ov(2)+dv(2)*(nv(2)-1)
    max3=ov(3)+dv(3)*(nv(3)-1)
    dt=d1
  end subroutine
!it=is the time
!ii=is the inclination
!ia=is the azimuth
!pVV=is the velocity
!pRR=is the ray
!
!
  logical function shoot3d_rays(rays,sources,angx,angy)
    real rays(:,:), sources(:),angx,angy
    real xn,yn,zn,rn,xm,ym,zm,rm
    real(kind=wp) denom,x1,y1,z1,x2,y2,z2
    real(kind=wp) ax,bx,ay,by,a,b,cs,delta
    type(point3D) T,C,S1,S2,S,N1,N2,M1,M2
    real a1,a2,LL,LL0
    type(v3D) v1,v2,vT,vx,vz
    real aa,bb
    real deps
    real wico,dist,tempv
    logical ::end
    integer :: nt,i1,i2,i3,i,ia
    shoot3d_rays=.false.
    wico=2.
    nt=size(rays,2)
    if (size(rays,1).ne.6) then
      call erexit("expected dimension 1 of rays to be 6 ")
    end if
!we will paralize of this one
!first do the initialization (not here)
    C%x=sources(1)
    C%y=sources(2)
    C%z=sources(3)
    tempv=get_vel_lin(C)
    rays(1:4,1)=(/c%x,c%y,c%z,tempv/)
    dist=tempv*dt
    c%x=dist*sin(pi/180.*angx)
    c%y=dist*sin(pi/180.*angy)
    c%z=-dist*sqrt(cos(pi/180.*angx)**2+cos(pi/180*angy**2))
! c%z=c%z*cos(pi/180.*angx)/abs(cos(pi/180.*angx))
    rays(1:3,2)=(/c%x,c%y,c%z/)+rays(1:3,1)
    rays(4,2)=tempv
    i1=2
    end=.false.
    do while (i1 <= nt-1 .and. .not. end)
      C%x=rays(1,i1)
      T%x=rays(1,i1-1)
      C%y=rays(2,i1)
      T%y=rays(2,i1-1)
      C%z=rays(3,i1)
      T%z=rays(3,i1-1)
      C%v=rays(4,i1)
      T%v=rays(4,i1-1)
      C%t=rays(5,i1)
      T%t=rays(5,i1-1)
      C%l=rays(6,i1)
      T%l=rays(6,i1-1)
!unitary vectors
      vz%dx=0
      vz%dy=0
      vz%dz=1
      vx%dx=1
      vx%dy=0
      vx%dz=0
      vT%dx=C%x-T%x
      vT%dy=C%y-T%y
      vT%dz=C%z-T%z
      LL =wico*v_l(vT)
! the vector T->C
      aa=v_a(vT,vz)
! angle between TC and the vertical axis
      vT%dz=0
      LL0=wico*v_l(vT)
! horizontal projection of TC
      if (LL0.ne.0) then
        bb=signum(v_mp(vT,vx,vz))*v_a(vT,vx)
      else
!angle between the horizontal projection of TC and the x axis
        bb=0
      end if
      aa=PI*aa/180.
      bb=PI*bb/180.
!points in the vertical plane containing the ray
      N1%x=C%x+LL*cos(aa)*cos(bb)
      N2%x=C%x-LL*cos(aa)*cos(bb)
      N1%y=C%y-LL*cos(aa)*sin(bb)
      N2%y=C%y+LL*cos(aa)*sin(bb)
      N1%z=C%z-LL*sin(aa)
      N2%z=C%z+LL*sin(aa)
      N1%v=get_vel_lin(N1)
      N2%v=get_vel_lin(N2)
!points in a horizontal plane
      M1%x=C%x-LL*sin(bb)
      M2%x=C%x+LL*sin(bb)
      M1%y=C%y-LL*cos(bb)
      M2%y=C%y+LL*cos(bb)
      M1%z=C%z
      M2%z=C%z
      M1%v=get_vel_lin(M1)
      M2%v=get_vel_lin(M2)
      xn= N2%x-N1%x
      xm= M2%x-M1%x
      yn= N2%y-N1%y
      ym= M2%y-M1%y
      zn= N2%z-N1%z
      zm= M2%z-M1%z
      rn=(N1%v-N2%v)*dt
      rm=(M1%v-M2%v)*dt
      deps=0.01*v_l(vT)
      denom=xn*ym-xm*yn
      if (denom.eq.0) then
        denom=deps
      end if
      ax=(rn*ym-rm*yn)/denom
      bx=(zn*ym-zm*yn)/denom
      ay=(xn*rm-xm*rn)/denom
      by=(xn*zm-xm*zn)/denom
      a=1+bx*bx+by*by
      b=ax*bx+ay*by
      cs=ax*ax+ay*ay-1
      delta=b*b-a*cs
      if (delta<0) then
        delta=0
        write(0,*) "Trouble! Imaginary solutions...\n"
        return
      end if
! line intersects sphere = 1
      z1=(b+sqrt(delta))/a
      S1%z=(C%z)+(C%v)*dt*z1
      x1=ax-bx*z1
      S1%x=(C%x)+(C%v)*dt*x1
      y1=ay-by*z1
      S1%y=(C%y)+(C%v)*dt*y1
! line intersects sphere = 2
      z2=(b-sqrt(delta))/a
      S2%z=(C%z)+(C%v)*dt*z2
      x2=ax-bx*z2
      S2%x=(C%x)+(C%v)*dt*x2
      y2=ay-by*z2
      S2%y=(C%y)+(C%v)*dt*y2
      vT%dx= T%x-C%x
      v1%dx=S1%x-C%x
      v2%dx=S2%x-C%x
      vT%dy= T%y-C%y
      v1%dy=S1%y-C%y
      v2%dy=S2%y-C%y
      vT%dz= T%z-C%z
      v1%dz=S1%z-C%z
      v2%dz=S2%z-C%z
      a1=v_a(v1,vT)
      a2=v_a(v2,vT)
      if (a1>a2) then
        S=S1
      else
        S=S2
      end if
      S%v=get_vel_lin(S)
      rays(1,i1+1)=s%x
      rays(2,i1+1)=s%y
      rays(3,i1+1)=s%z
      rays(4,i1+1)=s%v
      rays(5,i1+1)=s%t
      rays(6,i1+1)=s%l
      if (rays(3,i1+1) < 0. .and. i1+1 < nt) then
        do i=i1+1,nt
          rays(:,i)=rays(:,i1+1)
        end do
        end=.true.
      end if
      i1= i1+1
    end do
    shoot3d_rays=.true.
  end function
  real(kind=wp) function v_l(pv)
    type(v3d) :: pv
!returns the length of a vector
    v_l=sqrt ( (pV%dx)*(pV%dx) + (pV%dy)*(pV%dy) + (pV%dz)*(pV%dz) )
  end function
  integer function signum(n)
    real(kind=wp) :: n
!returns the sign of a double (positive if zero)
    if (n>=0) then
      signum= 1
    else
      signum= -1
    end if
  end function
  real(kind=wp) function v_a(pV1,pV2)
    type(v3D) :: pV1, pV2
    integer sdp
! returns the measure of the angle between two vectors
    sdp=signum(v_cosa(pV1,pV2))
    v_a= 180*((1-sdp)*PI/2+sdp*asin(real(v_sina(pV1,pV2))))/PI
  end function
  real (kind=wp) function v_cosa(pV1,pV2)
    type(v3D) :: pV1, pV2
! returns the cosine of the angle between two vectors
    v_cosa=v_dp(pV1,pV2)/(v_l(pV1)*v_l(pV2))
  end function
  real(kind=wp) function v_dp(pV1,pV2)
    type(v3D) :: pV1,pV2
    v_dp=(pV1%dx)*(pV2%dx) + (pV1%dy)*(pV2%dy) + (pV1%dz)*(pV2%dz)
  end function
  real(kind=wp) function v_sina(pV1,pV2)
    type(v3D) pV1,pV2
!returns the sine of the angle between two vectors
    type(v3D) v
    v=v_vp(pV1,pV2)
    v_sina=v_l(v)/(v_l(pV1)*v_l(pV2))
  end function
  type(v3d) function v_vp(pv1,pv2)
    type(v3d) :: pv1,pv2
!returns the vector product of two vectors
    v_vp%dx=(pV1%dy*pV2%dz) - (pV2%dy*pV1%dz)
    v_vp%dy=(pV2%dx*pV1%dz) - (pV1%dx*pV2%dz)
    v_vp%dz=(pV1%dx*pV2%dy) - (pV2%dx*pV1%dy)
  end function
  real(kind=wp) function v_mp(pV1,pV2,pV3)
    type(v3D) pV1, pV2, pV3
! returns the signed mixt product of three vectors (v1 x v2) * v3
    v_mp= (pV3%dx)* (pV1%dy)*(pV2%dz) -(pV3%dx)*(pV2%dy)*(pV1%dz) +&
      &(pV2%dx)*(pV3%dy)*(pV1%dz) -(pV2%dx)*(pV1%dy)*(pV3%dz) +&
      &(pV1%dx)*(pV2%dy)*(pV3%dz) -(pV1%dx)*(pV3%dy)*(pV2%dz)
  end function
  real(kind=wp) function get_vel_lin(p)
    type(point3d) :: p
    real(kind=wp) x,y,z
    integer i1,i2,i3
    real(kind=wp) v
    real(kind=wp) v000,v0x0,v00y,v0xy,vz00,vzx0,vz0y,vzxy
    real(kind=wp) v00,v0y,vz0,vzy,v0,vz
    real(kind=wp) m,n,q
    if (p%z<=min1) then
      z=min1
    end if
    if (p%z>=max1-dv(1)) then
      z=max1-dt
    end if
    if (p%z>min1 .and. p%z<max1-dv(1)) then
      z=p%z
    end if
    if (p%x<=min2) then
      x=min2
    end if
    if (p%x>=max2-dv(2)) then
      x=max2-dv(2)
    end if
    if (p%x>min2 .and. p%x<max2-dv(2)) then
      x=p%x
    end if
    if (p%y<=min3) then
      y=min3
    end if
    if (p%y>=max3-dv(3)) then
      y=max3-dv(3)
    end if
    if (p%y>min3 .and. p%y<max3-dv(3)) then
      y=p%y
    end if
    i1=floor((z-ov(1))/dv(1))+1
    i2=floor((x-ov(2))/dv(2))+1
    i3=floor((y-ov(3))/dv(3))+1
    v000=1/VEL(i1,i2,i3)
    v0x0=1/VEL(i1,i2+1,i3)
    v00y=1/VEL(i1,i2,i3+1)
    v0xy=1/VEL(i1,i2+1,i3+1)
    vz00=1/VEL(i1+1,i2,i3)
    vzx0=1/VEL(i1+1,i2+1,i3)
    vz0y=1/VEL(i1+1,i2,i3+1)
    vzxy=1/VEL(i1+1,i2+1,i3+1)
    m=(x-(ov(2)+(i2-1)*dv(2)))/dv(2)
    n=(y-(ov(3)+(i3-1)*dv(3)))/dv(3)
    q=(z-(ov(1)+(i1-1)*dv(1)))/dv(1)
    v00=(1-m)*v000+m*v0x0
    v0y=(1-m)*v00y+m*v0xy
    vz0=(1-m)*vz00+m*vzx0
    vzy=(1-m)*vz0y+m*vzxy
    v0 =(1-n)*v00 +n*v0y
    vz =(1-n)*vz0 +n*vzy
    v=(1-q)*v0 +q*vz
    get_vel_lin=1./v
  end function
end module
