!----------------------------------------------------------------------
!  nbltes -- test nbline; generate beam line scatter at target plane
!
subroutine nbltes
  use nbi_com
  use nbi_random
  use sglib_mod
  use nbdep_debug, only : znubeam_nbltes
 !
  implicit none
  LOGICAL, EXTERNAL :: LCHAR,LXOUT,LYES,LYESC                !
  !
  real(fp) :: zdrrms(mib),zdzrms(mib),zdrav(mib),zdzav(mib)
  !
  integer, parameter :: imaxpl=250
  real(fp) :: zrpl(imaxpl),zzpl(imaxpl)
  !
  real(fp) :: zr2a(5),zz2a(5)
  !
  character(len=80) :: zlbl1,zlbl2,zlbl3

  integer :: instep                    !   new stuff  sds 10/91
  parameter(instep=400)             !
  real(fp) :: z_zgrid(instep)              !
  real(fp) :: z_rgrid(instep)              !
  real(fp) :: zstep                        !
  real(fp) :: z_zsum(instep)               !
  real(fp) :: z_rsum(instep)               !
  real(fp) :: z_totsum                     !
  real(fp) :: z_max                        !
  !
  integer, parameter :: maxpath=250
  real(fp) :: xpath(maxpath),ypath(maxpath),zpath(maxpath)
  real(fp) :: zdist(maxpath),zwk1(maxpath),zwk2(maxpath),zwk3(maxpath)
  integer :: iwk(maxpath)
  !
  integer, parameter :: maxtrack=300000
  real(fp) :: zx0(maxtrack),zy0(maxtrack),zz0(maxtrack)
  real(fp) :: zx1(maxtrack),zy1(maxtrack),zz1(maxtrack)
  !
  integer :: lnout
  integer :: idcod
  !
  integer :: itrks,itest,ib,ilunt,j,izb,iab,iert,ians,itrk,ipl
  integer :: ifrac,ierr,i,isb,itype,ith,icount,imode
  real(fp) :: zrtcen,zrtarg,zl,zlhor,zdydl,zcosy,zdz,zdz2,zdr,zdr2,zz,zr
  real(fp) :: zindr,zindz,zdyl,zsiny,zwgtot,ztheta,zrho,zdum1,zdum2
  real(fp) :: zphilen,zlen,zlp
  !
  real(fp) :: zxmin,zxmax,zymin,zymax,zdiff
  real(fp) :: zrmin,zrmax,zzmin,zzmax
  real(fp) :: zrminp,zrmaxp,zzminp,zzmaxp
  !
  logical :: ipol,iready
  integer :: iseed,iscreen,ind,idum
  !-----------------------------------------------
  !  from nbline
  !
  real(fp) :: zzgrd,zrgrd,zdzdl,zdrdl
  !
  !-----------------------------------------------
  !
  zzgrd=znubeam_nbltes(1)
  zrgrd=znubeam_nbltes(2)
  zdzdl=znubeam_nbltes(3)
  zdrdl=znubeam_nbltes(4)

  iseed=abs(123454321*nbrandm(idum))
  iscreen=1620
  ind=21
  call nbi_random_init_dep(iseed,ind,iscreen)

  itest=1
  !
  zstep = 1.0_fp  ! cm
  z_max = (instep/2)*zstep          ! maximum elevation in grid

  do ib=1,mib
     zdrrms(ib)=0.0_fp
     zdzrms(ib)=0.0_fp
     zdrav(ib)=0.0_fp
     zdzav(ib)=0.0_fp

  end do
  !
  ilunt=lnout(0)
  !
20 continue

  z_totsum = 0._fp

  do j=1,instep
     z_zsum(j) = 0._fp
     z_rsum(j) = 0._fp

  end do

  call neupag
  if(lxout(0)) then
     write(ilunt,1001) nbeam
1001 format(/' number of beams = ',i2/ &
          ' #  ab zb co? rtcena      <dr**2>**.5 <dz**2>**.5 <dr>', &
          '           <dz>')
     do ib=1,nbeam
        izb=xzbeama(ib)+.1_fp
        iab=abeama(ib)+.5_fp
        write(ilunt,1002) ib,iab,izb,nlco(ib),rtcena(ib), &
             zdrrms(ib),zdzrms(ib),zdrav(ib),zdzav(ib)
1002    format(1x,i2,2(1x,i2),2x,l1,2x,5(1pe11.4,1x))

     end do
  end if
  !
  if(lxout(0)) then
     call out_test
     write(ilunt,*) ' enter "t" to change test'
     if(itest.ne.3) &
          write(ilunt,*) ' enter "s" to preset plot scale.'
     write(ilunt,*) ' enter "q" to quit.'
  end if
  call uread( &
       'nbltes:  enter beam # for testing, or "t" or "s" or "q":$')
  !
  if(lchar('q')) return
  !
  if(lchar('s')) then
     if(itest.eq.3) then
        write(6,*) ' test no. 3: "s" option ignored.'
     else
        call gscdef
     end if
     goto 20
  end if
  !
  if(lchar('t')) then
     call ureadi('enter # of desired test (1 2 3 or 4):$',itest)
     itest=max(1,min(4,itest))
     goto 20
  end if
  !
  ians=idcod(iert)
  if(iert.ne.0) goto 20
  if((ians.le.0).or.(ians.gt.nbeam)) then
     call uermsi('%nbltes, beam no. out of range',ians,1)
     goto 20
  end if
  !
  ib=ians
  !
  if(itest.eq.1) then
     itrks=1000
     call ureadi( &
          'nbltes:  enter # of beamlets to track (statistical control):$', &
          itrks)
     !
     if(itrks.le.0) goto 20
     !
     zrtcen=rtcena(ib)
     !
     call nbline(ib,1)
     !
     if(lxout(0)) write(ilunt,1003) rsta,zrtcen
1003 format(/' aperture radius = ',1pe11.4,' cm'/ &
          ' tangency radius = ',1pe11.4,' cm'// &
          ' % radius of target plane ctr must lie between these two...')
     !
     zrtarg=zrtcen
     call ureadr8('nbltes:  enter target plane center radius:$', &
          zrtarg)
     !
     if((zrtarg.lt.zrtcen).or.(zrtarg.gt.rsta)) then
        call uermsr8('%nbltes:  target plane radius out of range', &
             zrtarg,1)
        goto 20
     end if
     !
     if(zrtarg.eq.zrtcen) then
        zl=xlbtna(ib)
     else
        zlhor=sqrt(rsta**2-zrtcen**2)-sqrt(rsta**2-zrtarg**2)
        zdyl=xybapa(ib)-xybsca(ib)
        zsiny=zdyl/xlbapa(ib)
        zcosy=sqrt(1.0_fp-zsiny**2)
        zl=zlhor/zcosy+xlbapa(ib)
     end if
     !
     zdz=0.0_fp
     zdz2=0.0_fp
     zdr=0.0_fp
     zdr2=0.0_fp
     !
     !  distance zl from source pt to test plane assumed not to change; actual
     !  variation is indeed small since the beamlets are all nearly parallel
     !
     do itrk=1,itrks
        call nbline(ib,0)

        zz=zdzdl*zl+zzgrd
        zr=zdrdl*zl+zrgrd

        if(itrk.le.imaxpl) then

           zzpl(itrk) = zz
           zrpl(itrk) = zr

        end if

        zindz = (zz+z_max)/zstep
        zindr = (zr+z_max)/zstep
        if(((zindz.gt.0.0_fp).and.(zindz.lt.2.0_fp*z_max)) .and. &
             ((zindr.gt.0.0_fp).and.(zindr.lt.2.0_fp*z_max))) then

           z_totsum =  z_totsum + 1

           j         =  1 + int(zindz) ! increment z bin
           z_zsum(j) =  z_zsum(j) + 1._fp

           j         =  1 + int(zindr) ! increment r bin
           z_rsum(j) =  z_rsum(j) + 1._fp

        end if
        !
        zdz=zdz+zz
        zdz2=zdz2+zz**2
        zdr=zdr+zr
        zdr2 = zdr2 + zr**2
        !

     end do
     !
     do j=1,instep

        z_zsum(j)  = z_zsum(j)/z_totsum
        z_rsum(j)  = z_rsum(j)/z_totsum

        z_zgrid(j) = zstep * (-instep/2 + (j-0.5_fp) )
        z_rgrid(j) = zstep * (-instep/2 + (j-0.5_fp) )


     end do

     zdzav(ib)=zdz/itrks
     zdrav(ib)=zdr/itrks
     zdzrms(ib)=sqrt(zdz2/itrks)
     zdrrms(ib)=sqrt(zdr2/itrks)
     !
     zr2a(1)=zdrrms(ib)
     zr2a(2)=zdrrms(ib)
     zr2a(3)=-zdrrms(ib)
     zr2a(4)=-zdrrms(ib)
     zr2a(5)=zdrrms(ib)
     !
     zz2a(1)=zdzrms(ib)
     zz2a(2)=-zdzrms(ib)
     zz2a(3)=-zdzrms(ib)
     zz2a(4)=zdzrms(ib)
     zz2a(5)=zdzrms(ib)
     !
     ipl=min(imaxpl,itrks)
     !
     zlbl1=' '
     zlbl2=' '
     zlbl3=' '
     write(zlbl1,2001) ib,nlco(ib),zrtcen,zrtarg
     write(zlbl2,2002) zdrav(ib),zdzav(ib),itrks,runid(1:lrunid)
     write(zlbl3,2003) zdrrms(ib),zdzrms(ib)
     !
2001 format( &
          '(0,0)=ctrline ib=',i2,' co=',l1, &
          ' rtan=',1pe11.4,' rtarg=',1pe11.4)
     !
2002 format('  <zdr>=',1pe11.4,'  <zdz>=',1pe11.4,' n=',i8,1x,a)
2003 format(' sqrt(<zdr**2>)=',1pe11.4,' sqrt(<zdz**2>)=',1pe11.4)

     !
     call r8_grafxf(zr2a,zz2a,5,zrpl,zzpl,ipl,'dr(cm)','dz(cm)', &
          zlbl1,zlbl2,zlbl3)
     !
     !  new plots of distribution vs x and y   ... sds 10/91
     !
     call r8_grafx1(z_zgrid,z_zsum,instep,'z (cm)','p(z)', &
          zlbl1,zlbl2,zlbl3)

     call r8_grafx1(z_rgrid,z_rsum,instep,'r (cm)','p(r)', &
          zlbl1,zlbl2,zlbl3)

  else if(itest.eq.2) then
     call set_lsbeam

     call nbdep_debug_activate
     call nbdep_debug_center
     call nbdep_debug_setbeam(ib)

     if(abs(xzbeami-1.0_fp).lt.0.1_fp) then
        ifrac=1
        call ureadi('nbltes: enter 1 2 or 3 for full/half or'// &
             '1/3 energy beam:$',ifrac)
        ifrac=max(1,min(3,ifrac))
        call nbdep_debug_set_ien(ifrac)
     end if

     zwgtot=0
     call nbidep(zwgtot,ierr)
     if(ierr.ne.0) &
          write(6,*) ' ?nbltes: dep call returned ierr=',ierr

     call nbdep_debug_plot

     call nbdep_debug_set_ien(0)
     call nbdep_debug_setbeam(0)
     call nbdep_debug_decenter
     call nbdep_debug_deactivate

  else if(itest.eq.3) then
     !
     !  dmc -- projection plots added oct 2003
     !
     itrks=min(100,maxtrack)
     ipol=.true.
201  continue                       ! re-enter her for reconfiguration
     if(lxout(0)) then
        write(6,*) ' %max no. of tracks for projection test: ', &
             maxtrack
        write(6,2011) itrks, ipol, .not.ipol
2011    format(/' options:'/ &
             ' n -- set number of tracks to plot (currently: ',i6,')'/ &
             ' p -- select poloidal plane projection (',l1,')'/ &
             ' h -- select horizontal plane projection (',l1,')'/ &
             ' a -- accept settings, proceed to make plot.'/ &
             ' q -- quit this menu.'/ &
             )
     end if
     call uread('nbltes: choose projection plot option:$')
     iready=.false.
     if(lchar('a')) then
        iready=.true.
     else if(lchar('n')) then
        call ureadi('nbltes: number of tracks to generate:$',itrks)
     else if(lchar('p')) then
        ipol=.true.
     else if(lchar('h')) then
        ipol=.false.
     else if(lchar('q')) then
        goto 20
     else
        call uermsg('?nbltes: unrecognized option code.',' ',0)
     end if
     if(.not.iready) goto 201

     if(.not.ipol) then
        !     centerline calculation for plot scale only, horizontal plane
        call nbline(ib,1)
     end if

     call eq_bdlims(itype,zrmin,zrmax,zzmin,zzmax,ierr)
     if(ierr.ne.0) then
        write(6,*) ' ?nbltes: eq_bdlims ierr=',ierr
        call uclear
        goto 20
     end if
     !
     ! *** careful ***
     !   xplasma uses meters, but nubeam uses cm, for length units!
     !
     if(.not.ipol) then
        zxmax=1.05_fp*max(rsta,100.0_fp*zrmax)
        zxmin=-zxmax
        zymax=zxmax
        zymin=zxmin
     else
        zxmin=zrmin*100
        zxmax=zrmax*100
        zymin=zzmin*100
        zymax=zzmax*100
        call make_square
     end if

     !  compute tracks now...
     do i=1,itrks
        call nbline(ib,0)
        zx0(i)=rsta*cos(esta*degrad)
        zy0(i)=rsta*sin(esta*degrad)
        zz0(i)=ysta

        !  aiming angle, looking back in from starting point, with phi
        !  adjustment (sensitive to direction of current).

        zphilen=(esta*degrad+xpi)-nsnccw*nsjdotb*phia*degrad

        !  trial length
        zlen=2*max(100.0_fp*zrmax,rsta)*0.01_fp

        !  compute trial paths
        xpath(1)=zx0(i)*0.01_fp
        ypath(1)=zy0(i)*0.01_fp
        zpath(1)=zz0(i)*0.01_fp
        do j=2,maxpath
           zl=(j-1)*zlen/(maxpath-1)
           xpath(j)=xpath(1)+zl*cos(thea*degrad)*cos(zphilen)
           ypath(j)=ypath(1)+zl*cos(thea*degrad)*sin(zphilen)
           zpath(j)=zpath(1)+zl*sin(thea*degrad)
           zdist(j)=0.0_fp
        end do
        !
        !  compute relationship to limiters
        !
        call eqx_dist(maxpath,xpath,ypath,zpath,zdist,0, &
             zwk1,zwk2,zwk3,ierr)
        if(ierr.ne.0) write(6,*) ' ?nbltes: eqx_dist call error.'
        !
        !  scan for entry to and then exit from vv

        imode=0
        zl=0.0_fp
        zx1(1)=zx0(1)/100.0_fp
        zy1(1)=zy0(1)/100.0_fp
        zz1(1)=zz0(1)/100.0_fp
        do j=1,maxpath
           zlp=zl
           zl=(j-1)*zlen/(maxpath-1)
           if(imode.eq.0) then
              if(zdist(j).lt.0.0_fp) imode=1
           end if
           if(imode.eq.1) then
              if(zdist(j).gt.0.0_fp) then
                 imode=2
                 zlen=(zlp*zdist(j)-zl*zdist(j-1))/ &
                      (zdist(j)-zdist(j-1))
                 zx1(i)=xpath(1)+zlen*cos(thea*degrad)*cos(zphilen)
                 zy1(i)=ypath(1)+zlen*cos(thea*degrad)*sin(zphilen)
                 zz1(i)=zpath(1)+zlen*sin(thea*degrad)
                 exit
              end if
           end if
        end do
        if(imode.eq.0) then
           write(6,*) ' track ',i,': entry to plasma not found.'
        else if(imode.eq.1) then
           write(6,*) ' track ',i,': exit from plasma not found.'
        end if
     end do
     zx1(1:itrks)=100.0_fp*zx1(1:itrks)
     zy1(1:itrks)=100.0_fp*zy1(1:itrks)
     zz1(1:itrks)=100.0_fp*zz1(1:itrks)

202  continue                       ! re-enter here for rescaled plot

     call initt(0)
     call binitt
     call erase

     call slimx(120,720)
     call slimy(100,700)

     call r8_dlimx(zxmin,zxmax)
     call r8_dlimy(zymin,zymax)

     call xneat(1)
     call yneat(1)

     call xfrm(2)
     call yfrm(2)

     icount=0
     if(ipol) then
        call plot1(4)               ! vv -- long dash
        zlbl2='poloidal plane projection'
     else
        call plot2(100.0_fp*zrmin,100.0_fp*zrmax,4) ! vv -- long dash
        zlbl2='horizontal plane projection'
     end if

     do i=1,5
        zrho=min(0.99999999_fp,i*0.2_fp)
        if(ipol) then
           call plot1(1)            ! plasma -- short dash
        else
           call eq_glimrz(zrho,zrminp,zrmaxp,zzminp,zzmaxp,ierr)
           if(ierr.ne.0) then
              write(6,*) ' ?nbltes: eq_glimrz ierr=',ierr
              call uclear
           end if
           call plot2(100.0_fp*zrminp,100.0_fp*zrmaxp,1) ! plasma -- short dash
        end if
     end do

     do i=1,itrks
        do j=1,maxpath
           xpath(j)=((maxpath-j)*zx0(i)+(j-1)*zx1(i))/(maxpath-1)
           ypath(j)=((maxpath-j)*zy0(i)+(j-1)*zy1(i))/(maxpath-1)
           zpath(j)=((maxpath-j)*zz0(i)+(j-1)*zz1(i))/(maxpath-1)
           if(ipol) then
              xpath(j)=sqrt(xpath(j)**2+ypath(j)**2)
              ypath(j)=zpath(j)
           end if
        end do
        call plot4(0)
     end do

     zlbl1=' '
     if(nlco(ib)) then
        write(zlbl1,2021) ib,rtcena(ib)
2021    format('nbline test, beam#',i2,' coinjected; rtan=', &
             1pe11.4,' cm')
     else
        write(zlbl1,2022) ib,rtcena(ib)
2022    format('nbline test, beam#',i2,' ctrinjected; rtan=', &
             1pe11.4,' cm')
     end if

     call agtcgt(0,zlbl1)
     call agtchs(zlbl2)
     call agtcvs('cm')

     call agtmhs(735,650,'long dash:')
     call agtmhs(735,628,'  vacuum vessel')

     call agtmhs(735,584,'short dash:')
     call agtmhs(735,562,'  plasma')

     call agtmhs(735,518,'solid lines:')
     call agtmhs(735,496,'  beam tracks')

     call agtmhs(735,222,'beam tracks start')
     if(ipol) then
        call agtmhs(735,200,'at aperture.')
     else
        call agtmhs(735,200,'at source.')
     end if

     call agtmhs(735,156,'viewed from above:')
     if(nsnccw*nsjdotb.eq.1) then
        call agtmhs(735,134,'j is ctr-clockwise')
     else
        call agtmhs(735,134,'j is clockwise')
     end if

     call upause(0)
     call finitt(0,0)

     call uread('nbltes: change scale of plot? (y/n):$')
     if(lyes(0)) then
        do
           zdum1=zxmin
           zdum2=zxmax
           call ureadr8('nbltes: new xmin:$',zdum1)
           call ureadr8('nbltes: new xmax:$',zdum2)
           if(zdum1.eq.zdum2) then
              call uermsg(' x limits must be distinct.',' ',0)
              cycle
           end if
           zxmin=min(zdum1,zdum2)
           zxmax=max(zdum1,zdum2)
           exit
        end do
        do
           zdum1=zymin
           zdum2=zymax
           call ureadr8('nbltes: new ymin:$',zdum1)
           call ureadr8('nbltes: new ymax:$',zdum2)
           if(zdum1.eq.zdum2) then
              call uermsg(' y limits must be distinct.',' ',0)
              cycle
           end if
           zymin=min(zdum1,zdum2)
           zymax=max(zdum1,zdum2)
           exit
        end do
        !
        call make_square
        !
        goto 202
        !
     end if

  end if

  goto 20
  !------------------------
contains
  subroutine out_test
    use iso_c_binding, only: fp => c_double

    character(len=5) :: arrow
    character(len=60) :: ztests(4)

    ztests(1)='normal plane scatter test'
    ztests(2)='beam centerline test'
    ztests(3)='beamlet projection tests'

    write(6,*) ' current test indicated with arrow:'
    do i=1,3
       if (i.eq.itest) then
          arrow='====>'
       else
          arrow=' '
       end if
       write(ilunt,8881) arrow,i,ztests(i)
8881   format(1x,a,1x,i1,'.  ',a)
    end do

  end subroutine out_test

  subroutine make_square
    use iso_c_binding, only: fp => c_double
    !
    !  enforce squareness
    !
    zdiff=(zxmax-zxmin)-(zymax-zymin)
    if(zdiff.gt.0) then
       zymax=zymax+zdiff/2
       zymin=zymin-zdiff/2
    else
       zxmax=zxmax-zdiff/2
       zxmin=zxmin+zdiff/2
    end if

  end subroutine make_square

  subroutine set_lsbeam
    use iso_c_binding, only: fp => c_double

    ! last minute setup for dep call...

    do isb=1,nsbeam
       if(abs(abeama(ib)-abeams(isb)).lt.0.5_fp) then
          if(abs(xzbeama(ib)-xzbeams(isb)).lt.0.5_fp) then
             lsbeam=isb
             call inispc
          end if
       end if
    end do

  end subroutine set_lsbeam

  subroutine plot1(iline)
    use iso_c_binding, only: fp => c_double
    integer, intent(in) :: iline    ! line style

    integer :: inum
    real(fp) :: ztol

    icount = icount + 1

    ztol=1.0e-6_fp

    if(icount.eq.1) then
       call eq_limcon(maxpath,inum,xpath,ypath,ztol,ierr)
       if(ierr.ne.0) then
          write(6,*) ' ?nbltes: eq_limcon error.'
          inum=10
          xpath=0
          ypath=0
       else
          xpath(1:inum)=100.0_fp*xpath(1:inum)
          ypath(1:inum)=100.0_fp*ypath(1:inum)
       end if
       call npts(inum)
       call line(iline)
       call r8_check(xpath,ypath)
       call r8_dsplay(xpath,ypath)
       call frame

    else

       do ith=1,maxpath
          ztheta=(ith-1)*twopi/(maxpath-1)
          zwk1(ith)=zrho
          zwk2(ith)=ztheta
          zwk3(ith)=0.0_fp
       end do
       call eq_rzget(maxpath,zwk1,zwk2,zwk3,xpath,ypath, &
            iwk,ierr)
       if(ierr.ne.0) then
          write(6,*) ' ?nbltes: eq_rzget error.'
          xpath=0
          ypath=0
       else
          xpath=100.0_fp*xpath
          ypath=100.0_fp*ypath
       end if
       call npts(maxpath)
       call line(iline)
       call r8_cplot(xpath,ypath)
    end if

  end subroutine plot1

  subroutine plot2(zrmin,zrmax,iline)
    use iso_c_binding, only: fp => c_double
    real(fp),intent(in) :: zrmin,zrmax  ! radii of circles to plot
    integer, intent(in) :: iline    ! line style

    integer :: ir
    real(fp) :: zr(2)

    zr(1)=zrmax
    zr(2)=zrmin
    do ir=1,2
       do ith=1,maxpath
          ztheta=(ith-1)*twopi/(maxpath-1)
          xpath(ith)=zr(ir)*cos(ztheta)
          ypath(ith)=zr(ir)*sin(ztheta)
       end do
       icount = icount + 1
       if(icount.eq.1) then
          call npts(maxpath)
          call line(iline)
          call r8_check(xpath,ypath)
          call r8_dsplay(xpath,ypath)
          call frame
       else
          call npts(maxpath)
          call line(iline)
          call r8_cplot(xpath,ypath)
       end if
    end do

  end subroutine plot2

  subroutine plot4(iline)
    use iso_c_binding, only: fp => c_double
    integer, intent(in) :: iline    ! line style

    call npts(maxpath)
    call line(iline)
    call r8_cplot(xpath,ypath)

  end subroutine plot4

end subroutine nbltes
