  !******************** start file nbline.for ; group nbline *************
  !
  !	this subroutine mocks up a neutral beam line.
  !	for each call it provides a randomly selected
  !	neutral path, characterized in terms of
  !	rsta, ysta, phia, thea
  !
  !	beam neutrals are started randomly across the injector
  !	face. they are given random angles, based on the beam
  !	divergeance, and tracked to the given aperture.
  !	if a particle is skimmed off at the aperture, a new
  !	trajectory is chosen, until a viable particle is found.
  !	then the neutral trajectory for tracking across the plasma
  !	is calculated.
  !
  !	rob goldston 6-15-82
  !
  !  mod dmc 12-19-84:
  !    generalized for beamlines whose ctrs do not lie in and/or are
  !  not parallel to the midplane.  beam grids and apertures are
  !  assumed to lie in planes perp to beam ctr line
  !
  !  mod dmc may 1985 - if icline=1, set up a sightline down the
  !   middle of the beam
  !
  !*  mod jec aug 2002
  !*      prevent infinite loop if aperture incorrectly set
  !*      and print aperture parameters..
  !*
  ! mg oct 2013
  ! considering just the z coordinate of a rectangular source.
  ! a point z is chosen over the beam width at the source.
  ! associated with this is an angle alpha=atan(z/f) which would direct the particle
  ! through the focal point with focal length f.
  ! to account for divergence, alpha is actually chosen as
  !  tan(alpha) = p(z/f, divza)
  ! where p(x,sigma) return a random value x* which has mean x and standard deviation sigma
  ! over a gaussian distribution. the particle must then make it through the aperture
  ! with this random alpha in order for it to be deposited in the plasma.

  subroutine nbline(ib0,icline)
    !
    !
    !  include transp common blocks
    use nbi_com
    use nbi_random_dep
    use nbdep_debug, only : znubeam_nbltes
    !
    !
    !  input information from trans.com
    !  --------------------------------
    !
    !	foclr
    !	foclz
    !	divr
    !	divz
    !	nbshap
    !	xlbapa - source to aperture distance, by beamline
    !	xlbtna - source to tangency radius distance, by beamline
    !	xybsca - dist (cm) ctr of source above/below midplane
    !  	xybapa - dist (cm) ctr of aperture above/below midplane
    !	bmwidr
    !	bmwidz
    !
    !  ** modif dmc oct 1984
    !     ib0 = beam line index passed to nbline
    !     (converted to ib locally; ib=ib0 for main beams; for trace element
    !     beams ib=abs(ntrace(ib0))...)
    !
    !	rtcena - tangency radius, by beamline
    !	cofrac
    !	xzedge -- superseded dmc may 1991
    !	redge  -- superseded dmc may 1991 - default maintained sr nbpdat
    !	degrad
    !
    !  ** mod dmc may 1991
    !     beam aperture description elaborated
    !  1st aperture for each beam:
    !       nbapsha(..)  =1 rectangular  =2 circular
    !         xlbapa(..),xybapa(..) as before
    !       xzpedga(..)
    !       rapedga(..)
    !  2nd aperture for each beam (optional):
    !       nbapsh2(..)  =0 none(default) =1 rectangular  =2 circular
    !       xlbapa2(..)
    !       xzpedg2(..)
    !       rapedg2(..)
    !
    !  mod dmc aug 2009 -- per user request: aperture offsets added;
    !    horizontal and vertical offsets for each of the two apertures;
    !    default values (0.0) means apertures centered around beam centerlines
    !    as before
    !
    !  output information to trans.com
    !  -------------------------------
    !
    !	rsta
    !	ysta
    !	phia
    !	thea
    !
    !  local common for debug routine nbltes
    !
    !
    !============
    implicit none
    integer :: icline,ib0,ib,idum
    !============
    real(fp) :: zzgrd,zrgrd,zdzdl,zdrdl,zdza,zzap,zdra,zrap,zdzdl0
    real(fp) :: zdrdl0,zradap,zdza2,zzap2,zdra2,zrap2,zradap2,zdymp
    real(fp) :: zsiny,zcosy,zthet0,zlbmp,zlbmpa,zhgt,zphil,zphi0,zbeta
    !============
    !
    integer ::   intry, intrylim
    parameter ( intrylim = 1000 )  !  give up after .. attempts
    !*                                   !  ( observed 8 ( max ) for jet run
    !*d    uncomment these lines to monitor number of tries
    !*d    ( & see end of icline block below )
    !*
    !*d    integer ::   intrymax
    !*d    save      intrymax
    !*d    data      intrymax / 1 /
    !
    !----------------------------------------------------------------
    !
    !  select from beam face
    !  ---------------------
    !
    !  return point for cases of beam scrape off
    !
    ib=ib0
    if(ntrace(ib0).ne.0) ib=abs(ntrace(ib0))
    !
    intry  =  0
    !
10  continue
    !----------------------------------
    intry  =  intry + 1
    if( intry .gt. intrylim )   then
       !..       no good trajectory found ..  problem with aperture defn ?
       write(nonlin,6000) ib, intry
       write(nonlin,6001) 1, nbapsha(ib),&
            xzpedga(ib), rapedga(ib), bmwidra(ib)
       if( xlbapa2(ib) .gt. 0._fp)  then
          write(nonlin,6001) 2, nbapsh2(ib), xzpedg2(ib), rapedg2(ib)
       end if
       !
       write(nonlin,*)&
            'nubeam/nbline : cant find good beam trajectory'
       
       rsta = -1.0_fp
       znubeam_nbltes(1:4)=zero

       return

    end if
6000 format(10x,'nubeam/nbline, beam ',i4,i6,' attempts ' )
6001 format(10x,'aperture',i2,&
          ' nbapsha,xzpedga,rapedga,bmwidra ',i4,3f12.3 )
    !
    if(icline.eq.1) then
       zzgrd=0.0_fp
       zrgrd=0.0_fp
       zdzdl=0.0_fp
       zdrdl=0.0_fp
       zdza=0.0_fp
       zzap=0.0_fp
       zdra=0.0_fp
       zrap=0.0_fp
    else
       if(nbshapa(ib).eq.1) goto 100
       !
       !  circular beam grids
50     continue
       zzgrd=bmwidra(ib)*(1.0_fp-2.0_fp*nbrandm_dep(idum))
       zrgrd=bmwidra(ib)*(1.0_fp-2.0_fp*nbrandm_dep(idum))
       !--     if(zzgrd**2+zrgrd**2.gt.bmwidra(ib)**2) goto 50
       if(zzgrd**2+zrgrd**2.gt.bmwidra(ib)**2) then
          !--        another loop trap..
          intry  =  intry + 1
          if( intry .gt. intrylim )  goto  10
          goto 50
       else
          goto 150
       end if
       !
       !  rectangular beam grids
100    continue
       zzgrd=bmwidza(ib)*(1.0_fp-2.0_fp*nbrandm_dep(idum))
       zrgrd=bmwidra(ib)*(1.0_fp-2.0_fp*nbrandm_dep(idum))
       !
       !  establish neutral trajectory within beam-line
       !  ---------------------------------------------
       !
       !	zdzdl and zdrdl are the tangents of the
       !	angles between the normal to the (assumed flat) grid
       !	surface and the particle trajectory, projected into
       !	the vertical and horizontal planes.
       !	i.e. dl in this definition is length along the grid-normal line
       !
150    continue
       if(foclra(ib).eq.0.0_fp) then
          write(nonlin,*)&
               'nubeam/nbline : input data error foclra(',ib,')= 0'

          rsta = -1.0_fp
          znubeam_nbltes(1:4)=zero
 
          return
       end if
       if(nbshapa(ib).eq.1) goto 1100
       !
       !  circular beam grids
       zdzdl0=-zzgrd/foclra(ib)
       zdrdl0=-zrgrd/foclra(ib)
       zdzdl=nbgauss_dep(divra(ib),zdzdl0)
       zdrdl=nbgauss_dep(divra(ib),zdrdl0)
       goto 1150
       !
       !  rectangular beam grids
1100   continue
       zdzdl0=-zzgrd/foclza(ib)
       zdrdl0=-zrgrd/foclra(ib)
       zdzdl=nbgauss_dep(divza(ib),zdzdl0)
       zdrdl=nbgauss_dep(divra(ib),zdrdl0)
       !
       !  scrape off beam at rectangular aperture
       !  ---------------------------------------
       !
1150   continue
       !
       !      	call expert(l1nbline,l2nbline,1)
       !       we no longer need expert calls jed 7/13/01
       !
       !
       zdza=zdzdl*xlbapa(ib)
       zzap=zzgrd+zdza
       zdra=zdrdl*xlbapa(ib)
       zrap=zrgrd+zdra
       if(nbapsha(ib).eq.1) then
          !  rectangular aperture
          if(abs(zzap-xzapoffa(ib)).gt.xzpedga(ib)) goto 10
          if(abs(zrap-xrapoffa(ib)).gt.rapedga(ib)) goto 10
       else
          !  circular aperture
          zradap=sqrt((zzap-xzapoffa(ib))**2+(zrap-xrapoffa(ib))**2)
          if(zradap.gt.rapedga(ib)) goto 10
       end if
       !
       !  optional 2nd aperture
       !
       if((nbapsh2(ib).gt.0).and.(xlbapa2(ib).gt.0.0_fp)) then
          zdza2=zdzdl*xlbapa2(ib)
          zzap2=zzgrd+zdza2
          zdra2=zdrdl*xlbapa2(ib)
          zrap2=zrgrd+zdra2
          if(nbapsh2(ib).eq.1) then
             !  rectangular aperture
             if(abs(zzap2-xzapoff2(ib)).gt.xzpedg2(ib)) goto 10
             if(abs(zrap2-xrapoff2(ib)).gt.rapedg2(ib)) goto 10
          else
             !  circular aperture
             zradap2= &
                  sqrt((zzap2-xzapoff2(ib))**2+(zrap2-xrapoff2(ib))**2)
             if(zradap2.gt.rapedg2(ib)) goto 10
          end if
          !
       end if
       !
       !  dmc - end of icline if block
    end if
    !----------------------------------
    !
    !      	call expert(l1nbline,l2nbline,2)
    !       expert no longer needed -jed 7/13/01
    !
    !d    if( intry .gt. intrymax )  then
    !d       intrymax = intry
    !d       print *,' nubeam/nbline intrymax = ',intrymax
    !d    end if
    !
    !  got a viable particle trajectory
    !  --------------------------------
    !
    !  dmc-- deal with beamline not pll to midplane
    !
    zdymp=xybapa(ib)-xybsca(ib)
    if(zdymp.eq.0.0_fp) then
       zsiny=0.0_fp
       zcosy=1.0_fp
       zthet0=0.0_fp
    else
       zsiny=zdymp/xlbapa(ib)
       zcosy=sqrt(1.0_fp-zsiny*zsiny)
       zthet0=atan(zsiny/zcosy)
    end if
    !
    !  following quantities include correction for deviation of selected
    !  sightline from beam centerline
    !
    !  distance, aperture to tangency radius, projected to midplane
    !  with correction for tilt of aperture plane
    !   note sign convention: zsiny.lt.0.0 if beamline is coming down
    !  from above (corresponding to a negative theta argument for
    !  trakrb).
    !
    zlbmp=(xlbtna(ib)-xlbapa(ib))*zcosy + zzap*zsiny
    !
    !  distance, sce to aperture, projected to midplane
    !
    zlbmpa=xlbapa(ib)*zcosy - zdza*zsiny
    !
    !  find ysta and rsta
    !
    if(nlsym2b) then
       zhgt=plhgt
    else
       zhgt=0.0_fp
    end if

    ysta=zzap*zcosy+xybapa(ib)-zhgt
    rsta=sqrt((rtcena(ib)-zrap)**2+zlbmp**2)
    !
    !  find phia and thea
    !
    !  positive zrap, zdra make beamlet more perp (==> smaller phia)
    !
    zphil = atan(rtcena(ib)/(zcosy*xlbtna(ib)))  ! angle from beam sightline to axis
    ! on equatorial plane
    !
    !  component, sce to aperture
    zphi0=atan(zdra/zlbmpa)
    !  component, from pt. crossing aperture plane on in
    phia=atan((rtcena(ib)-zrap)/zlbmp)
    !  add
    esta=phia-zphil   ! zeta angle relative to source
    phia=phia-zphi0
    !
    if(.not.nlco(ib)) then
       phia=-phia
       esta=-esta
    end if
    phia=phia/degrad
    esta=xbzeta(ib) + nsnccw*nsjdotb*esta/degrad  ! correct for ip direction
    ! make angle absolute in zeta
    !
    !	note that thea is **not** simply tan(zdzdl) since
    !	dl needs to be projected from the trajectory down to the
    !	horizontal midplane first, rather than being the grid-normal
    !
    !  dmc-- add theta contribution due to beamline itself
    !
    thea=atan(zdzdl/sqrt(1.0_fp+zdrdl**2))+zthet0
    thea=thea/degrad
    !
    !  dmc june 1992 -- support deposition in a box fast neutral density
    !  3d calculation for chers group
    !
    if(mdepbox.lt.0 .or. ndeptube.lt.0) then
       !
       !  special deposition call for neutral density calculation
       !
       xbgrid=zrgrd
       ybgrid=zzgrd
       dxdlbox=zdrdl
       dydlbox=zdzdl
       !
       !  compute length traveled along 3d box axis (= beam axis) per
       !  length traveled by *this* beam neutral
       zbeta=atan(sqrt(zdrdl**2+zdzdl**2))
       dlbxdtk=cos(zbeta)
       !
    end if
    !
    znubeam_nbltes(1)=zzgrd
    znubeam_nbltes(2)=zrgrd
    znubeam_nbltes(3)=zdzdl
    znubeam_nbltes(4)=zdrdl
    return
  end subroutine nbline
  !******************** end file nbline.for ; group nbline ***************
