e150posb.inp

pic 模拟程序！面向对象

e150posb
{

Modeling the E-150 plasma lens experiment at SLAC --

High-energy positron bunch enters a dense neutral gas in cylindrical 
geometry. The positron beam is Gaussian in z and r.

Boundary conditions:
a)  The simulation region must be bounded by either conductors or
    insulators, in order to capture lost particles.
b)  Conductors were chosen, to avoid any charge build up.
c)  The choice of conducting boundary conditions means that electric
    fields parallel to the boundaries are forced to zero;  however,
    fields near the boundaries of the simulation must be small anyway
    to accurately model a semi-infinite plasma, so this is OK.

Gas jet:
a)  The MCC algorithm is used to model a gas jet, which is
    confined to a subregion of the grid.
b)  The MCC model is turned off after the beam has fully traversed
    the jet.  Otherwise, the ionization cascade creates an exponentially
    growing number of ions and secondary electrons.

Moving window:
a)  Once the positron beam is beyond the gas jet, a moving window 
    algorithm is invoked so that the beam can then be modeled for 
    relatively long times.
b)  This has the added benefit of removing the now-irrelevant
    ions and secondary electrons from the simulation, which
    speeds things up quite a bit.
}

// Define variables that can be used throughout this input file.

Variables
{
// First, define some useful constants.
	pi = 3.14159265358979323846
	speedOfLight = 2.99792458e+08
	electronMass = 9.1093897e-31
	unitCharge = electronMass * 1.75881962e11
	electronCharge = -1. * unitCharge
	electronMassEV = electronMass * speedOfLight * speedOfLight / unitCharge
	ionCharge = unitCharge
	unitMassMKS = electronMass / 5.48579903e-04
	lithiumMassNum = 6.942
	lithiumMass = unitMassMKS * lithiumMassNum
	hydrogenMassNum = 1.00797
	hydrogenMass = unitMassMKS * hydrogenMassNum
	nitrogenMassNum = 14.007
	nitrogenMass = unitMassMKS * nitrogenMassNum

// Next, define the parameters of the high-energy electron beam.
	beamEnergyEV = 28.5e+09
  beamGammaMin1 = beamEnergyEV / electronMassEV
  beamGamma = 1 + beamGammaMin1
  beamBetaGamma = sqrt( beamGammaMin1 * (beamGammaMin1+2) )
  beamBeta = beamBetaGamma / beamGamma

	totalNumBeam = 1.5e+10
	totalBeamCharge = totalNumBeam * electronCharge
	rmsBeamRadius = 0.08e-04
	rmsBeamLength = 6.00e-04
	rmsBeamTime = rmsBeamLength / speedOfLight
	radialCutoffFac = 2.25
	axialCutoffFac = 2.6
	totalBeamRadius = radialCutoffFac * rmsBeamRadius
	totalBeamLength = 2 * axialCutoffFac * rmsBeamLength
	beamAspectRatio = totalBeamLength / totalBeamRadius
  totalBeamArea = pi * totalBeamRadius * totalBeamRadius
  rmsBeamVolume = pi * rmsBeamRadius * rmsBeamRadius * rmsBeamLength

  rmsEnergySpread = 0.001
	beamTempEV = rmsEnergySpread * beamEnergyEV
	thermalBeamSpeedEV = 0.5 * beamTempEV
  rmsNormalizedEmittance = 4.0e-06
  rmsBeamSize = rmsBeamRadius / sqrt(2)
  rmsThermalBeta = rmsNormalizedEmittance / rmsBeamSize
  rmsThermalGamma = 1. / sqrt(1.-rmsThermalBeta*rmsThermalBeta)
  rmsVelocityMKS = rmsThermalBeta * speedOfLight
  rmsVelocityEV = (rmsThermalGamma-1.)*electronMassEV
  rmsEfactor = 8.0e-04
  rmsVfactor = 1.0e-04

// Define the number of grids in R and Z
	simRadiusOverBeamRadius = 4
  totalSimulationLength = 0.0072
	simLengthOverBeamLength = totalSimulationLength / totalBeamLength
	numRgridsAcrossBeam = 6
	numZgridsAcrossBeam = 104
	numRgrids = numRgridsAcrossBeam * simRadiusOverBeamRadius
	numZgrids = numZgridsAcrossBeam * simLengthOverBeamLength
	numCells = numRgrids * numZgrids

// Number of beam particles
	numBeamPtclsPerCell = 1000
	numBeamCells = numRgridsAcrossBeam * numZgridsAcrossBeam
	numBeamPtcls = numBeamPtclsPerCell * numBeamCells
	beamNumRatio = totalNumBeam / numBeamPtcls

// Intermediate calculations for modeling Gaussian shape of the beam.
	invSigRsq = 1.0 / ( rmsBeamRadius * rmsBeamRadius )
	invSigZsq = 0.5 / ( rmsBeamLength * rmsBeamLength )
	invSigTsq = invSigZsq * speedOfLight * speedOfLight

// Calculate the size of the simulation region, grid spacings, time step.
// We are assuming the same grid size in both z and r	
  maxRadiusMKS = simRadiusOverBeamRadius * totalBeamRadius
  maxLengthMKS = simLengthOverBeamLength * totalBeamLength
	rGridSize = maxRadiusMKS / numRgrids
	zGridSize = maxLengthMKS / numZgrids
  dx = 1. / sqrt( 1./(rGridSize*rGridSize) + 1./(zGridSize*zGridSize) )
	timeStep = 0.7 * dx / speedOfLight

// Calculate peak currents for defining emission of the high-energy beam.
  peakCurrentDensity=totalBeamCharge*speedOfLight/rmsBeamVolume/sqrt(2.*pi)
	peakCurrent = peakCurrentDensity * totalBeamArea
	pulseLengthSec = totalBeamLength / speedOfLight
  oneHalfPulse = pulseLengthSec/2.
  oneEighthPulse = pulseLengthSec/8.
  threeEighthsPulse = 3.*oneEighthPulse
  sevenEighthsPulse = 7.*oneEighthPulse

// Define gas density, pressure and other MCC parameters
	gasTempEV       = 1.0e-06  // make gas cold (cannot set temperature to zero)
  nitrogenDensity = 2.2e+24  // density of neutral N2 (molecular nitrogen) in m^-3
                             // density of nitrogen atoms is twice N2 density
  gasDensityMKS   = 2.0 * nitrogenDensity
  gasPressureTorr = 1.20e-21 * gasDensityMKS * gasTempEV

  Min1MCC = 0.0006  // leading  edge of gas jet
  Max1MCC = 0.0036  // trailing edge of gas jet
  Avg1MCC = 0.5 * (Min1MCC + Max1MCC)  // center of gas jet
  RMS1MCC = (Max1MCC - Avg1MCC)/3.     // RMS width of gas jet
  Min2MCC = 0.
  Max2MCC = 2.*totalBeamRadius  // Gas jet does not go out to full radius.
                                // This reduces # of ions and secondary electrons

// Calculate when the MCC algorithm should turn on and off.
// Turn the moving window on after the MCC algorithm turns off.
  delayTimeMCC = Min1MCC / speedOfLight
  stopTimeMCC  = (Max1MCC+totalBeamLength) / speedOfLight
  windowDelay  = 0.98 * maxLengthMKS / speedOfLight

// Define parameters related to the adhoc photoionization by a "laser" pulse,
//   which models the copropagating optical field
	totalNumPhotoElectrons = 1.2e-04 * totalNumBeam
  laserPlasmaLength = 0.001
	laserPlasmaRadius = 2. * rmsBeamRadius
  laserPlasmaVolume = laserPlasmaLength*pi*laserPlasmaRadius*laserPlasmaRadius
  laserNumDensity     = totalNumPhotoElectrons / laserPlasmaVolume
  laserPlasmaNumRatio = 200.
  laserDelayTime = Min1MCC / speedOfLight
}

// This simulation has only one "region", which contains grid, all particles, etc.
Region
{

// Define the grid for this region.
Grid
{

// Define number of grids along Z-axis and physical coordinates.
	J = numZgrids
	x1s = 0.0
	x1f = maxLengthMKS
	n1 = 1.0

// Define number of grids along R-axis and physical coordinates.
	K = numRgrids
	x2s = 0.0
	x2f = maxRadiusMKS
	n2 = 1.0
}

// Specify "control" parameters for this region
Control
{
// Specify the time step.
	dt = timeStep

// Turn on the moving window algorithm.
  movingWindow   = 1
  shiftDelayTime = windowDelay

// Turn on damping for the high-frequency EM fields
	emdamping = 0.49
}

// Define the nitrogen ions.
Species
{
	name = ions
	m = nitrogenMass
	q = ionCharge
  subcycle = 100
  particleLimit = 8.0e+05 // prevents out-of-control growth in # of ptcls
}

// Define the secondary electrons.
Species
{
	name = electrons
	m = electronMass
	q = electronCharge
  collisionModel = 1
  particleLimit = 8.0e+05 // prevents out-of-control growth in # of ptcls
}

// Define the beam positrons.
Species
{
	name = positrons
	m = electronMass
	q = unitCharge
  collisionModel = 1
  rmsDiagnosticsFlag = 1
}

// Define the beam emitter, which introduces the high-energy beam into the
// simulation.

VarWeightBeamEmitter
{
  speciesName = positrons
  I = peakCurrent

// Define the 2-D function F(x,t) that specifies beam emission profile.
  xtFlag = 3
  F=exp(-invSigRsq*x*x)*exp(-invSigTsq*(t-oneHalfPulse)*(t-oneHalfPulse))*step(pulseLengthSec-t)

// Macroparticles are emitted from the left boundary, close to the axis of symmetry.
	j1 = 0
	j2 = 0
	k1 = 0
	k2 = numRgridsAcrossBeam
	normal = 1
	np2c = beamNumRatio

// Emit particles, directed along the Z-axis,  with specified energy and temperature.
	units = EV
	v1drift = beamEnergyEV
	v1thermal = rmsEfactor * rmsVelocityEV
	v2thermal = rmsVfactor * rmsVelocityEV
//	v3thermal = rmsVfactor * rmsVelocityEV
}

// Specify the Monte Carlo collision parameters for background gas
MCC
{
  gas = N
  relativisticMCC = 1
  pressure = gasPressureTorr
  temperature = gasTempEV
  analyticF = gasDensityMKS * exp(-(x1-Avg1MCC)*(x1-Avg1MCC)/RMS1MCC/RMS1MCC)
  eSpecies = electrons
  iSpecies = ions

  x1MinMKS = Min1MCC
  x1MaxMKS = Max1MCC
  x2MinMKS = Min2MCC
  x2MaxMKS = Max2MCC

  delayTime = delayTimeMCC  // delay time in seconds before MCC kicks in
  stopTime  = stopTimeMCC   // time after which MCC is no longer used
}

// Define an adhoc "laser plasma source", which simulates photoionization
//   by an assumed co-propagating optical field (due to synch. rad.)

LaserPlasmaSource
{
  name = plasmaSource
  A1 = Min1MCC
  B1 = Min1MCC + laserPlasmaLength
  A2 = 0
  B2 = laserPlasmaRadius
  units1 = EV
  speciesName1 = electrons	
  temperature1 = 0
  units2 = EV
  speciesName2 = ions
  temperature2 = 0	
  sourceRate = laserNumDensity/pulseLengthSec  // rate of plasma generation
  np2c = laserPlasmaNumRatio
  velocity = speedOfLight                      // velocity of ionization pulse
  x1Flag = 1                                   // 1 = horizontal motion
  pulseWidth = totalBeamLength                 // spatial width of pulse in m
  delayTime = laserDelayTime  // delay time in seconds before laser pulse starts
}

// Specify a perfect conductor along the left boundary.  This serves as a particle
//   boundary condition (catches particles that leave the simulation) and as a
//   field boundary condition (E_r is forced to vanish).

Conductor
{
	j1 = 0
	j2 = 0
	k1 = 0
	k2 = numRgrids
	normal = 1
}

// Specify a perfect conductor along the radial boundary.  This serves as a
//   particle boundary condition (catches particles that leave the simulation)
//   and as a field boundary condition (E_z is forced to vanish).

Conductor
{
	j1 = 0
	j2 = numZgrids
	k1 = numRgrids
	k2 = numRgrids
	normal = -1
}

// Specify a perfect conductor along the right boundary.  This serves as a
//   particle boundary condition (catches particles that leave the simulation)
//   and as a field boundary condition (E_r is forced to vanish).

Conductor
{
	j1 = numZgrids
	j2 = numZgrids
	k1 = numRgrids
	k2 = 0
	normal = -1
}

// Define the cylindrical symmetry axis.
CylindricalAxis
{
	j1 = 0
	j2 = numZgrids
	k1 = 0
	k2 = 0
	normal = 1
}

}