!
! ---------------------------- nbline_tube -------------------------
! this subroutine tests the intersection of a beam with an infinite tube. the
! beam is considered to have a rectangular source and aperture.  it is better
! to accept a possible intersection then to fail to find an intersection.  therefore
! no attempt is made to consider whether the intersection happens outside
! the region of tube interest, outside the vacuum vessel, behind a limiter ...
!
! this subroutine should stay consistent with nbline.for.  the parameters
! are entered through arguments to allow its use from outside nubeam before
! nubeam initialization.  their interpretation should be the same as the
! corresponding nbispec.dat namelist variables.
!
subroutine nbline_tube(rtcena, xlbtna, xlbapa, xybsca, xybapa, xbzeta, nlco, &
     bmwidra, bmwidza, rapedga, xzpedga, &
     rtube, ytube, xzetatube, phitube, thetatube, rhotube, &
     ipccw, ib, &
     isintersect, cinter,i_cx)
 
 
  use iso_c_binding, only: fp => c_double
  use physconst_mod, only: z0=>zero,z1=>one,rad2deg,pi
  implicit none
  real(fp),  intent(in)  :: rtcena      ! tangential radius of beam sightline
  real(fp),  intent(in)  :: xlbtna      ! distance from source to tangent
  real(fp),  intent(in)  :: xlbapa      ! distance from source to aperture
  real(fp),  intent(in)  :: xybsca      ! height of source
  real(fp),  intent(in)  :: xybapa      ! height of aperture
  real(fp),  intent(in)  :: xbzeta      ! toroidal angle of source
  logical, intent(in)  :: nlco        ! true if beam points co-rotating with plasma
  real(fp),  intent(in)  :: bmwidra     ! half width of source
  real(fp),  intent(in)  :: bmwidza     ! half height of source
  real(fp),  intent(in)  :: rapedga     ! half width of aperture
  real(fp),  intent(in)  :: xzpedga     ! half height of aperture
  real(fp),  intent(in)  :: rtube       ! radius of tube start
  real(fp),  intent(in)  :: ytube       ! height of tube start
  real(fp),  intent(in)  :: xzetatube   ! toroidal angle of start of tube
  real(fp),  intent(in)  :: phitube     ! azimuthal sightline of tube in +ip direction (deg)
  real(fp),  intent(in)  :: thetatube   ! elevation sightline of tube (deg)
  real(fp),  intent(in)  :: rhotube     ! radius of tube
  integer, intent(in)  :: ipccw       ! +1 if ip is ccw
  integer, intent(in)  :: ib          ! beam line index
  integer, intent(in)  :: i_cx          ! cx sightline, added for better diagnostic (mg)
  logical, intent(out) :: isintersect ! .true. if tube intersects beam
  real(fp),  intent(out) :: cinter      ! interception distance or zero
 
  logical, parameter :: debug = .false.
  real(fp), parameter :: zsmall  = 1.e-8_fp  ! a very small number
 
  integer :: nnco                    ! +1 if beam points with the current
  real(fp)  :: dy                      ! height of aperture above source
  real(fp)  :: sin_thetab, cos_thetab  ! trig of thetab
  real(fp)  :: xproj                   ! projection of beam sightline on midplane
  real(fp)  :: rb                      ! radius of beam source
  real(fp)  :: yb                      ! height of beam source
  real(fp)  :: zetab_rad               ! toroidal angle of source (rad)
  real(fp)  :: phib_rad                ! azimuthal beam sightline in +toroidal (rad)
  real(fp)  :: s(3)                    ! source position
  real(fp)  :: xl                      ! distance from source to aperture
  real(fp)  :: dsrc_r, dsrc_z          ! source half height and width
  real(fp)  :: dap_r, dap_z            ! aperture half height and width
  real(fp)  :: rt                      ! radius of tube start
  real(fp)  :: yt                      ! height of tube start
  real(fp)  :: zetat_rad               ! toroidal angle of tube start (rad)
  real(fp)  :: thetat_rad              ! elevation angle of tube (rad)
  real(fp)  :: sin_thetat, cos_thetat  ! trig of thetat
  real(fp)  :: phit_rad                ! azimuth angle of tube
  real(fp)  :: a(3)                    ! location of tube start, later in beam frame
  real(fp)  :: b(3)                    ! unit vector along tube, later in beam frame
  real(fp)  :: as(3)                   ! tube start relative to source
  real(fp)  :: delta                   ! toroidal rotation of beam sightline
  real(fp)  :: sin_delta,cos_delta     ! trig of delta
  real(fp)  :: m(3,3)                  ! rotation matrix from lab frame to beam sightline frame
  real(fp)  :: rsum                    ! maximal r displacement from source to aperture
  real(fp)  :: zsum                    ! maximal z displacement from source to aperture
  real(fp)  :: xcross                  ! common term, looks like a cross product
  integer :: ihv                     ! loop on horizontal and vertical fans
  integer :: ipm                     ! loop on up/down or left/right fans
  real(fp)  :: e1(3), e2(3)            ! unit vectors along edge of fans
  real(fp)  :: p0(3)                   ! location of apex of fans
  real(fp)  :: e3(3)                   ! unit vector orthoganol to e1,e2
  real(fp)  :: ap(3)                   ! tube start relative to fan apex
  real(fp)  :: b3,left,right           ! temps
  real(fp)  :: drp(3)                  ! projection of ap onto fan plane
  real(fp)  :: y(3)                    ! projection of fan apex into tube orthoganol plane
  real(fp)  :: v1(3),v2(3)             ! e1,e2 projected into tube orthoganol plane
  real(fp)  :: y1,y2                   ! closest approach of fan edges to tube
  real(fp)  :: zz,zm                   ! temp
  real(fp)  :: ctan                    ! distance to tangency
  real(fp)  :: sinter(4)               ! distance to intercepts
  integer :: iinter                  ! number of intercepts
 
  integer :: io           ! temporary i/o unit number

  ! --- start ---
  isintersect = .false.  ! assume null hypothesis
  cinter=z0
 
  ! ------ establish beam data -------
  if (nlco) then
     nnco = 1
  else
     nnco = -1
  end if
 
  dy = xybapa-xybsca
 
  sin_thetab = max(-z1, min(z1,dy/xlbapa))
  cos_thetab = sqrt(z1 - sin_thetab**2)
 
  xproj    = cos_thetab*xlbtna
  rb       = sqrt(rtcena**2 + xproj**2)
  yb       = xybsca
  zetab_rad= xbzeta/rad2deg
 
  !sin_phib = (rtcena*ipccw*nnco)/rb
  !tan_phib = (rtcena*ipccw*nnco)/xproj
  !cos_phib = xproj/rb
  phib_rad  = btan2(rtcena*ipccw*nnco, xproj)
 
  s  = (/rb*cos(zetab_rad), rb*sin(zetab_rad), yb/)
 
  xl = xlbapa
 
  dsrc_r = bmwidra
  dsrc_z = bmwidza
 
  dap_r  = rapedga
  dap_z  = xzpedga
 
  ! ------- establish tube data ------
  rt        = rtube
  yt        = ytube
  zetat_rad = xzetatube/rad2deg
 
  thetat_rad = thetatube/rad2deg
  sin_thetat = sin(thetat_rad)
  cos_thetat = cos(thetat_rad)
 
  phit_rad   = ipccw*phitube/rad2deg
 
  a = (/rt*cos(zetat_rad), rt*sin(zetat_rad), yt/)
  b = (/-cos_thetat*cos(phit_rad-zetat_rad), cos_thetat*sin(phit_rad-zetat_rad), sin_thetat/)
 
  ctan = rt*cos(phit_rad)/cos_thetat
 
  if (debug) call gnuplot(ib)
 
  !
  ! ---------- rotation of tube into beam sightline coordinates ---------
  ! the beam sightline frame has the beam pointing along the -x axis
  ! with the source at the origin.  transformation to this frame consists
  ! of a rotation around z by (phib-zetab) then a rotation around y by
  ! -thetab (down)
  !
  as = a-s  ! tube start relative to source
 
  delta = phib_rad-zetab_rad
  cos_delta = cos(delta)
  sin_delta = sin(delta)
  m(1,:) = (/cos_delta*cos_thetab, -sin_delta*cos_thetab, -sin_thetab/)
  m(2,:) = (/sin_delta,             cos_delta,             z0        /)
  m(3,:) = (/cos_delta*sin_thetab, -sin_delta*sin_thetab,  cos_thetab/)
 
  a = matmul(m,as)   ! a now becomes tube start in beam sightline frame
  b = matmul(m,b)    ! b now becomes tube direction in beam sightline frame
 
  !
  ! -------------------- test each fan -----------------
  ! the maixmal extent of the beam is formed by connecting opposite corners
  ! of the source rectangle with the aperture rectangle with lines.  a "fan" is the
  ! two dimensional cone formed by a pair of lines.  the lines have unit vectors
  ! e1 and e2 and intersect at the beam sightline coordinate p0.
  !
  rsum   = max(zsmall,dap_r+dsrc_r)
  zsum   = max(zsmall,dap_z+dsrc_z)
  xcross = (dsrc_z*dap_r-dsrc_r*dap_z)
  iinter = 0
  do ihv=1,2           ! loop horizontal,vertical
     do ipm=-1,1,2     ! loop +/-
        if (ihv==1) then
           ! ihv=1 -> horizontal
           e1 = normalize((/-xl,  rsum, -ipm*zsum/))
           e2 = normalize((/-xl, -rsum, -ipm*zsum/))
           p0 = (/-xl*dsrc_r/rsum, z0, ipm*xcross/rsum/)
        else
           ! ihv=2 -> vertical
           e1 = normalize((/-xl, -ipm*rsum, -zsum/))
           e2 = normalize((/-xl, -ipm*rsum,  zsum/))
           p0 = (/-xl*dsrc_z/zsum, -ipm*xcross/zsum, z0/)
        end if
 
        e3 = normalize(cross(e1,e2))
        ap = a-p0
 
        ! --- test for tube axis intersection with fan ---
        b3 = sum(b*e3)
        if (abs(b3)>zsmall) then
           drp = ap - b*(sum(ap*e3)/b3)
 
           left = sum(e1*drp)*sum(e2*drp)
           right= sum(e1*e2) *sum(drp*drp)
 
           if (left>=right) then
              isintersect=.true.
              iinter=iinter+1
              sinter(iinter) = -sum(ap*e3)/b3
              if (debug) print *, 'side intersection with beam ', ib, &
                   ' dist=',sinter(iinter)
           end if
        end if
 
        ! --- test tube within radius of edges of fan ---
        if (ihv==1) then                    ! redundant to test vertical, edges are the same
           y  = -ap + b*sum(ap*b)           ! projections orthogonal to tube
           v1 = normalize(e1-b*sum(e1*b))
           v2 = normalize(e2-b*sum(e2*b))
 
           if (norm(v1)>zsmall) then
              y1 = norm(cross(y,v1))
           else
              y1 = norm(y)
           end if
 
           if (norm(v2)>zsmall) then
              y2 = norm(cross(y,v2))
           else
              y2 = norm(y)
           end if
 
           zz = min(y1,y2)
           if (debug .and. zz<=4._fp*rhotube) print *, ' min dist beam ', ib, ' = ', zz
           if (zz<=rhotube) then
              isintersect=.true.
              if (debug) print *, 'radial intersection with beam ', ib
           end if
        end if
     end do    ! end +/- loop
  end do       ! end horizontal,vertical loop
 
  if (isintersect) then
     cinter = 2.0_fp*rhotube   ! minimum interception distance
     if (iinter>0) then
        zz = min(2.0_fp*ctan, maxval(sinter(1:iinter)))  ! furthest intercept but not too far
        zm = minval(sinter(1:iinter))
        if (zm<0._fp) then
           print *, '!nbline_tube: cx sightline #',i_cx,' instersects beam #',ib,&
                ' behind start of cx sightline -- check geometry or ignore'
           zm=0._fp
        end if
        cinter = max(cinter, zz-zm)
     end if
     if (debug) print *, 'estimated intercept dist with beam ', ib, ' = ', cinter
  end if
contains
  !
  ! return a unit vector
  !
  function normalize(vector) result(vout)
    real(fp), intent(in)              :: vector(:)  ! input vector to normalize
    real(fp), dimension(size(vector)) :: vout       ! normalized
    real(fp)                          :: vmag       ! magnitude
 
    vmag = sqrt(sum(vector*vector))
    if (vmag>z0) then
       vout = vector/vmag
    else
       vout = z0
    end if
  end function normalize
 
  !
  ! find the cross product of two vectors
  !
  function cross(va,vb) result(vout)
    real(fp), intent(in)  :: va(:), vb(:)    ! input vectors
    real(fp)              :: vout(3)         ! only makes sense for 3d
 
    vout(1) = va(2)*vb(3)-va(3)*vb(2)
    vout(2) = va(3)*vb(1)-va(1)*vb(3)
    vout(3) = va(1)*vb(2)-va(2)*vb(1)
  end function cross
 
  !
  ! return norm
  !
  function norm(vector) result (vmag)
    real(fp), intent(in)              :: vector(:)  ! input vector to normalize
    real(fp)                          :: vmag       ! magnitude
 
    vmag = sqrt(sum(vector*vector))
  end function norm
 
  !
  ! stupid fortran
  !
  function btan2(y,x) result(z)
    real(fp), intent(in) :: x,y
    real(fp)             :: z
    if (x==0._fp.and. y==0._fp) then
       z=0._fp
    else
       z=atan2(y,x)
    end if
  end function btan2
 
  !
  ! --------------- for testing --------------------
  ! writes out gnuplotable files
  !
  subroutine gnuplot(ig)
 
    integer,intent(in) :: ig
    integer :: i,j,k
    real(fp), dimension(3)    :: that,zhat,hhat,vhat,sap,au,bu,ba,c
    real(fp), dimension(3,50) :: ak,bk
    character*16 :: filenum
    real(fp)   ncross(3)
 
    that = (/ -cos_thetab*cos(phib_rad-zetab_rad), &
         cos_thetab*sin(phib_rad-zetab_rad), &
         sin_thetab /)
 
    zhat = (/0._fp, 0._fp, 1._fp/)
 
    ncross = cross(that,zhat)
    hhat = normalize(ncross)
    ncross = cross(hhat,that)
    vhat = normalize(ncross)
 
    sap = s + that*xl
 
    call find_io_unit(io)
    write(filenum, '(i5)') ig
    open(unit=io, file="gnuplot"//trim(adjustl(filenum))//".dat", status="replace")
 
 
    k=0
    do i=-1,1,2
       do j=-i*1,i*1,i*2
          au = s   + (i*dsrc_z)*vhat + (j*dsrc_r)*hhat
          bu = sap - (i*dap_z) *vhat - (j*dap_r) *hhat
          ba = normalize(bu-au)
 
          c = au + 1.5_fp*xproj*ba
 
          k=k+1
          ak(:,k) = au
          bk(:,k) = bu
 
          write(io,150)
          write(io,140) au,c
140       format(3(f15.6,2x)/3(f15.6,2x))
       end do
    end do
150 format(//)
 
    ! beam line
    c = s + 1.5_fp*xproj*that
    write(io,150)
    write(io,140) s,c
 
    ! write out source box
    ak(:,5)=ak(:,1)
    write(io,150)
    do k=1,5
       write(io,'(3(f15.6,2x))') ak(:,k)
    end do
 
    ! write out aperture box
    bk(:,5)=bk(:,1)
    write(io,150)
    do k=1,5
       write(io,'(3(f15.6,2x))') bk(:,k)
    end do
 
    ! tube
    ncross = cross(b,zhat)
    hhat = normalize(ncross)
    ncross = cross(hhat,b)
    vhat = normalize(ncross)
 
    write(io,150)
    c = a + xproj*b
    do k=1,17
       au = cos(2._fp*pi*(k-1)/16._fp)*vhat + sin(2._fp*pi*(k-1)/16._fp)*hhat
       ak(:,k) = a + rhotube*au
       bk(:,k) = c + rhotube*au
       write(io,140) ak(:,k),bk(:,k)
    end do
 
    write(io,150)
    do k=1,17
       write(io,'(3(f15.6,2x))') ak(:,k)
    end do
 
    write(io,150)
    do k=1,17
       write(io,'(3(f15.6,2x))') bk(:,k)
    end do
 
    write(io,150)
    write(io,140) a,c
 
    au = (/0._fp, 0._fp, -.5_fp*xproj/)
    bu = (/0._fp, 0._fp,  .5_fp*xproj/)
 
    ! axis
    write(io,150)
    write(io,140) au,bu
 
    close(unit=io)
  end subroutine gnuplot
end subroutine nbline_tube
! 11sep2004 fgtok -s r8_precision.sub nubeam_test_r8.sub "r8con.csh conversion"
