iondrv Subroutine

  subroutine iondrv(aaion,bmax,dpp,dtheta,edrive,emitt,etai,lf, &
       nbeams,qion,sigma,sigma0,tauf,theta,vi,gain,etadrv)

    !! Routine to calculate parameters of a heavy ion driver
    !! suitable for inertial fusion energy
    !! author: P J Knight, CCFE, Culham Science Centre
    !! aaion  : input real : Ion mass (amu)
    !! bmax   : input real : Maximum field at the superconductor (T)
    !! dpp    : input real : Beam momentum spread
    !! dtheta : input real : Pointing error (rad)
    !! edrive : input real : Driver energy (J)
    !! emitt  : input real : Normalised beam emittance (metre-rad)
    !! etai   : input real : Axial quadrupole packing fraction at injection
    !! lf     : input real : Distance from final focussing quad to target (m)
    !! nbeams : input integer : Number of beams
    !! qion   : input real : Ion charge state
    !! sigma  : input real : Depressed tune (incl. space charge effects) (deg)
    !! sigma0 : input real : Phase advance per lattice period (tune) (deg)
    !! tauf   : input real : Post-acceleration pulse duration (s)
    !! theta  : input real : Final focussing half angle (rad)
    !! vi     : input real : Injection voltage (V)
    !! gain   : output real : Target gain
    !! etadrv : output real : Driver efficiency
    !! This routine calculates the parameters of a heavy ion driver
    !! suitable for inertial fusion energy.
    !! <P>Currently, the complicated model taken from the reference is not
    !! complete, so it is recommended that the simple model is used
    !! (set <CODE>ISIMP=1</CODE>)
    !! Heavy-ion Driver Design and Scaling, R. Bieri et al.,
    !! Fusion Technology, vol.21 (1992) 1583
    !! Meier and Bieri, Fusion Technology, vol.21 (1992) 1547
    !!
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    use constants, only: electron_charge, proton_mass, degrad, pi
    implicit none

    !  Arguments

    real(dp), intent(in) :: aaion,bmax,dpp,dtheta,edrive,emitt,etai,lf,qion, &
         sigma,sigma0,tauf,theta,vi
    integer, intent(in) :: nbeams
    real(dp), intent(out) :: gain,etadrv

    !  Local variables

    real(dp), parameter :: c2 = 8.98755178737D16
    integer, parameter :: isimp = 1  !  Switch for simple model (1=yes)

    real(dp) :: ci,de,dgap,dlcore,drcore,e,eomc2,fins,floss, &
         ibf,ibfo,ibi,ibpc,lfocus,lpf,lpfo,lpi,lppc,lq,phif,phifo, &
         phii,phipc,rion,rs,rs1,rs2,rs3,rs4,sig,sig0,taufo,taui, &
         taupc,tbrad,vf,vfo,vpc,vs,vscore,vshe,vsle,vspc,xhe,xle,xpc
    real(dp), dimension(10) :: eve,gve
    integer :: i,ie,j,nche,ncle,ncpc,nqarrs,nqle,nqpche,nquads

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    if (isimp /= 1) then

       !  Complex IONDRV model (incomplete)...
       !  ====================

       !  Electron charge / (proton mass * c**2)

       eomc2 = electron_charge / (proton_mass*c2)

       !  Degrees to radians

       sig = sigma * degrad
       sig0 = sigma0 * degrad

       !  Beam current coefficient

       ci = cbeam(etai,bmax,aaion,qion,sig,sig0,emitt,vi,eomc2)

       !  Final voltage (V)

       vf = (edrive/(dble(nbeams)*tauf*ci))**(0.666666D0)

       !  Final beam current (A)

       ibf = ci*sqrt(vf)

       !  Ion range (g/cm^2)

       rion = (3.04D-5 * aaion + 349.0D0*aaion**(-2.21D0)) * &
            (qion*vf/1.0D9)**(1.72D0-0.00275D0*aaion)

       !  Convert to kg/m^2

       rion = rion*10.0D0

       !  Total spot size (rms of four terms):
       !  (1) Emittance term

       rs1 = emitt / (betgam(aaion,qion,vf)*theta)

       !  (2) Dispersion (momentum spread) term

       rs2 = 8.0D0*dpp*lf*theta

       !  (3) Alignment and jitter term

       rs3 = dtheta*lf

       !  (4) Space charge term (inaccurate at present)

       rs4 = 1.0D-3

       rs = sqrt(rs1**2 + rs2**2 + rs3**2 + rs4**2)

       !  Target gain (inaccurate at present)

       gain = 62.6D0 - 142.3D0*rs**1.25D0 * rion**0.633D0 + &
            (32.8D0+82.1D0*rs)**0.783D0 * rion**0.647D0 * &
            log(1.0D-6*edrive)

       !  Source/Injector characteristics:
       !  Transport beam radius (m)

       tbrad = ( (emitt**2 * sig0 * c2)/(sig**2 * bmax/1.5D0 * etai * &
            sqrt(vi)) )**(0.333333D0) * sqrt(2.0D0 * aaion * proton_mass / &
            (qion * electron_charge))

       !  Extractor voltage (V)

       vs = 4.0D0/3.0D0 * vi

       !  Source extraction gap width (m)

       dgap = 0.8D0

       !  Injected current per beam (A)

       ibi = 4.0D0 * pi * tbrad**2 * 5.46D-8 * &
            sqrt( qion*vs**3/(aaion*dgap**4) )

       !  Beam loss fraction between injector and target

       floss = 0.08D0

       !  Pulse duration at injection (s)

       taui = 1.0D0/(1.0D0-floss) * tauf * ibf / ibi

       !  Pulse length at injection (m)

       lpi = taui * sqrt( 2.0D0 * qion * electron_charge * vi / (aaion * proton_mass) )

       !  Initial voltage gradient (V/m)

       phii = 0.3D0 * vi / lpi

       !  End of Low Energy Transport Stage:
       !  Insulated fraction of driver length

       fins = 0.85D0

       !  Maximum voltage gradient, limited by insulator flash-over (V/m)

       phifo = 1.0D6 * fins

       !  Voltage at which this limit is reached (V)
       !  The end of the low energy stage is set by the point at which
       !  the limit is reached.

       vfo = (phifo/phii)**(2.0D0/3.0D0) * vi

       !  Beam current (A)

       ibfo = ibi * (vfo/vi)

       !  Pulse duration (s)

       taufo = taui * (vi/vfo)

       !  Pulse length (m)

       lpfo = taui * vi / vfo * sqrt( 2.0D0 * qion * electron_charge * vfo / &
            (aaion * proton_mass) )

       !  Length of the low energy transport stage (m)
       !  (rearrangement of integral of PHI(V).dl)

       xle = 2.0D0*vi**1.5D0 / phii * (1.0D0/sqrt(vi) - 1.0D0/sqrt(vfo))

       !  End of Pulse Compression Stage:
       !  Pulse length = final pulse length (m)

       lpf = tauf * sqrt( 2.0D0 * qion * electron_charge * vf / (aaion * proton_mass) )
       lppc = lpf

       !  Length of the pulse compression region (m)

       xpc = (lpfo - lpf)/0.3D0

       !  Voltage gradient (V/m) - limit already reached

       phipc = phifo

       !  Voltage (V)

       vpc = vfo + (xpc*phipc)

       !  Beam current (A)

       ibpc = ibfo * (lpfo/lppc) * sqrt(vpc/vfo)

       !  Pulse duration (s)

       taupc = taufo * (lppc/lpfo) * sqrt(vfo/vpc)

       !  End of High Energy Transport Stage:
       !  Voltage gradient (V/m) - limit already reached

       phif = phifo

       !  Length of the high energy transport stage (m)

       xhe = (vf - vfo)/phifo

       write(*,*) '                 Injector:'
       write(*,*) '   vi = ',vi
       write(*,*) '  ibi = ',ibi
       write(*,*) ' taui = ',taui
       write(*,*) ' phii = ',phii
       write(*,*) '  lpi = ',lpi
       write(*,*) '                 End of LET stage:'
       write(*,*) '  vfo = ',vfo
       write(*,*) ' ibfo = ',ibfo
       write(*,*) 'taufo = ',taufo
       write(*,*) 'phifo = ',phifo
       write(*,*) ' lpfo = ',lpfo
       write(*,*) '  xle = ',xle
       write(*,*) '                 End of PC stage:'
       write(*,*) '  vpc = ',vpc
       write(*,*) ' ibpc = ',ibpc
       write(*,*) 'taupc = ',taupc
       write(*,*) 'phipc = ',phipc
       write(*,*) ' lppc = ',lppc
       write(*,*) '  xpc = ',xpc
       write(*,*) '                 End of HET stage:'
       write(*,*) '   vf = ',vf
       write(*,*) '  ibf = ',ibf
       write(*,*) ' tauf = ',tauf
       write(*,*) ' phif = ',phif
       write(*,*) '  lpf = ',lpf
       write(*,*) '  xhe = ',xhe

       !  Volt-second requirements and number of cores

       !  Volt-seconds per core (2.4T is field swing of the core
       !  material metglas)

       drcore = 0.8D0
       dlcore = 0.2D0
       vscore = 2.4D0 * drcore * dlcore

       !  Low energy stage

       vsle = taui * vi * log(vfo/vi)
       ncle = 1 + int(vsle/vscore)

       !  Pulse compression stage

       vspc = 2.0D0 * tauf * sqrt(vf) * (sqrt(vpc)-sqrt(vfo))
       ncpc = 1 + int(vspc/vscore)

       !  High energy stage

       drcore = 0.4D0
       dlcore = 0.1D0
       vscore = 2.4D0 * DRCORE * DLCORE

       vshe = 2.0D0 * tauf * vf * (1.0D0 - sqrt(vpc/vf))
       nche = 1 + int(vshe/vscore)

       write(*,*) '                 Number of cores: '
       write(*,*) ' ncle = ',ncle
       write(*,*) ' ncpc = ',ncpc
       write(*,*) ' nche = ',nche

       !  Quadrupole requirements

       !  Effective quadrupole length (m)

       lq = ( emitt * etai * sig0**2 * sqrt(c2*vi) / &
            (sig*(bmax/1.5D0)**2) )**(0.333333D0) * &
            sqrt(2.0D0*aaion*proton_mass/(qion*electron_charge))

       write(*,*) '   lq = ',lq

       !  Actual quadrupole length (m)

       write(*,*) '  alq = ',lq/etai

       !  Length of quadrupole focussing fields: low energy stage

       lfocus = etai * vi/phii * (1.0D0 - vi/vfo)
       nqle = 1 + int(lfocus/lq)

       !  Pulse compression + high energy stages

       lfocus = 2.0D0 * etai * sqrt(vi)/phifo * (sqrt(vf) - sqrt(vfo))

       nqpche = 1 + int(lfocus/lq)

       !  Number of quadrupole arrays

       nqarrs = nqle + nqpche

       !  Number of quadrupoles

       nquads = nqarrs * nbeams

       write(*,*) 'nqarrs = ',nqarrs
       write(*,*) 'nquads = ',nquads

       etadrv = 0.28D0

    else

       !  Simple model
       !  ============

       !  gve(k): target gain at edrive = k MegaJoules

       gve(1)  = 25.0D0
       gve(2)  = 44.0D0
       gve(3)  = 62.0D0
       gve(4)  = 76.0D0
       gve(5)  = 87.0D0
       gve(6)  = 97.0D0
       gve(7)  = 107.0D0
       gve(8)  = 115.0D0
       gve(9)  = 125.0D0
       gve(10) = 132.0D0

       !  eve(k): driver efficiency at edrive = k MegaJoules

       eve(1)  = 0.232D0
       eve(2)  = 0.256D0
       eve(3)  = 0.269D0
       eve(4)  = 0.276D0
       eve(5)  = 0.282D0
       eve(6)  = 0.286D0
       eve(7)  = 0.290D0
       eve(8)  = 0.292D0
       eve(9)  = 0.294D0
       eve(10) = 0.296D0

       e = 1.0D-6 * edrive
       ie = int(e)
       de = e - dble(IE)

       !  Assume linear interpolations and extrapolations
       !  Would be better to prevent extrapolation

       if (ie <= 1) then

          gain   = gve(2) - 1.0D-6*(edrive-2.0D6)*(gve(1)-gve(2))
          etadrv = eve(2) - 1.0D-6*(edrive-2.0d6)*(eve(1)-eve(2))

       else if (ie >= 9) then

          gain   = gve(9) + 1.0D-6*(edrive-9.0D6)*(gve(10)-gve(9))
          etadrv = eve(9) + 1.0D-6*(edrive-9.0D6)*(eve(10)-eve(9))

       else

          gain   = gve(ie) + de*(gve(ie+1)-gve(ie))
          etadrv = eve(ie) + de*(eve(ie+1)-eve(ie))

       end if

       !  Ensure sensible values

       gain = max(0.01D0,gain)
       etadrv = max(0.01D0,etadrv)

    end if

  contains

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    function cbeam(etai,bmax,aaion,qion,sigma,sigma0,emitt,vi,eomc2)

      !! Routine to evaluate the beam current coefficient, CI
      !! author: P J Knight, CCFE, Culham Science Centre
      !! etai   : input real : Axial quadrupole packing fraction at injection
      !! bmax   : input real : Maximum field at the superconductor (T)
      !! aaion  : input real : Ion mass (amu)
      !! qion   : input real : Ion charge state
      !! sigma  : input real : Depressed tune (incl. space charge effects) (rad)
      !! sigma0 : input real : Phase advance per lattice period (tune) (rad)
      !! emitt  : input real : Normalised beam emittance (metre-rad)
      !! vi     : input real : Injection voltage (V)
      !! eomc2  : input real : Electron charge / (proton mass * light speed**2)
      !! This routine calculates the beam current coefficient
      !! of a heavy ion driver suitable for inertial fusion energy.
      !! Heavy-ion Driver Design and Scaling, R. Bieri et al.,
      !! Fusion Technology, vol.21 (1992) 1583
      !!       !
      ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      implicit none

      real(dp) :: cbeam

      !  Arguments

      real(dp), intent(in) :: etai,bmax,aaion,qion,sigma,sigma0, &
           emitt,vi,eomc2

      !  Local variables

      real(dp) :: bbe

      ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      !  Maximum field at the beam edge

      bbe = bmax/1.5D0

      !  Perform calculation

      cbeam = 2.89D6*(1.0D0-(sigma/sigma0)**2) * &
           ( sigma0**4 * etai**2 * vi * aaion/qion * (emitt/sigma)**2 &
           * bbe**2 )**0.333333D0 &
           * (2.0D0 * qion * eomc2 / aaion)**(5.0D0/6.0D0)

    end function cbeam

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    function betgam(aaion,qion,v)

      !! Routine to calculate the relativistic factor (beta*gamma) for the
      !! heavy ions in the beam driver
      !! author: P J Knight, CCFE, Culham Science Centre
      !! aaion  : input real : Ion mass (amu)
      !! qion   : input real : Ion charge state
      !! v      : input real : Acceleration voltage (V)
      !! This routine calculates the relativistic factor (beta*gamma) for the
      !! heavy ions in a beam driver suitable for inertial fusion energy.
      !! Heavy-ion Driver Design and Scaling, R. Bieri et al.,
      !! Fusion Technology, vol.21 (1992) 1583
      !!       !
      ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      implicit none

      real(dp) :: betgam

      !  Arguments

      real(dp), intent(in) :: aaion,qion,v

      !  Local variables

      real(dp), parameter :: c2 = 8.98755178737D16

      ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      betgam = sqrt( 2.0D0*qion*electron_charge*v / (aaion*proton_mass*c2) )

    end function betgam

  end subroutine iondrv