mcc.cpp

pic 模拟程序！面向对象

  iCX[0] = new FuncCrossSection(sigmaIelastic, 2, 0.0, CrossSection::ION_ELASTIC);
  iCX[1] = new FuncCrossSection(sigmaCX, 2, 0.0, CrossSection::CHARGE_EXCHANGE);
}

Scalar HeMCC::sigmaElastic(Scalar energy)
{
  return (8.5e-19/(pow(energy+10.0, 1.1)));
}

Scalar HeMCC::sigmaExc(Scalar energy)
{
  const Scalar extengy0 = 19.8;
  if(energy < extengy0) return (0.0);
  else if(extengy0 <= energy && energy <27.0)
    return (2.08e-22*(energy -extengy0));
  return (3.4e-19/(energy +200));
}

Scalar HeMCC::sigmaIz(Scalar energy)
{
  const Scalar ionengy0 = 24.5;
  if(energy < ionengy0) return (0.0);
	return(1e-17*(energy -ionengy0)/((energy +50)*pow(energy+300.0, 1.2)));
}

Scalar HeMCC::sigmaIelastic(Scalar energy)
{
  Scalar crossx;
  const Scalar emin = 1.0;
  if(energy < emin) { energy = emin; }

  crossx = 3.6463e-19/sqrt(energy) - 7.9897e-21;
  //  fprintf(stderr, "\nenergy=%g, crossx=%g", energy, crossx);
  if(crossx < 0) { return 0.0; }
  return crossx;
}

Scalar HeMCC::sigmaCX(Scalar energy)
{
  Scalar crossx;
  const Scalar emin = 1.0;
  const Scalar cutoff = 377.8;

  if(energy < emin) { energy = emin; }
  if(energy < cutoff)
    {
      crossx = 1.2996e-19  - 7.8872e-23*energy + 1.9873e-19/sqrt(energy);
    }
  else
    {
      crossx = (1.5554e-18*log(energy)/energy) + 8.6407e-20;
    }
  return crossx;
}

//-------------------------------------------------------------------

//-------------------------------------------------------------
// LiMCC: Lithium MCC package    -- added 03/02/00 by Bruhwiler

LiMCC::LiMCC(Scalar p, Scalar temp, Species* eSpecies, Species* iSpecies, 
             SpeciesList& sList, Scalar dt, int ionzFraction,
             Scalar ecxFactor, Scalar icxFactor, int relativisticFlag,
             SpatialRegion* _SR, Scalar _Min1MKS, Scalar _Max1MKS,
             Scalar _Min2MKS, Scalar _Max2MKS,
             Scalar _delayTime, Scalar _stopTime, const ostring &analyticF,
             const int& discardDumpFileNGDDataFlag)
           : MCCPackage(Li, (ostring)"Li", p, temp, eSpecies, iSpecies, sList, dt,
                        ionzFraction, ecxFactor, icxFactor, _SR,
                        _Min1MKS, _Max1MKS, _Min2MKS, _Max2MKS,
                        _delayTime, _stopTime, analyticF, discardDumpFileNGDDataFlag) {

  relativisticMCC = relativisticFlag;
  staticRelMCC    = relativisticFlag;
  addCrossSections();
  init(sList, dt);
}


// Define static data members for the LiMCC class.
  Scalar LiMCC::threshold     = 5.392;
  Scalar LiMCC::youngerParamA =  1.35e-17;
  Scalar LiMCC::youngerParamB = -6.51e-18;
  Scalar LiMCC::youngerParamC =  5.99e-19;
  Scalar LiMCC::youngerParamD = -8.49e-18;
  Scalar LiMCC::fitParamA1    = 4.9785e-20;
  Scalar LiMCC::fitParamA2    = 2.0862e-20;
  Scalar LiMCC::fitParamA3    = 6.0173e-20;
  Scalar LiMCC::fitParamA4    = 1.462;
  Scalar LiMCC::fitParamN     = 0.439;
  int    LiMCC::staticRelMCC  = 1;

//
void LiMCC::addCrossSections() {
    // Specify number of electron collision models that are implemented --
  	//   at present, we have impact ionization and elastic scattering.
  neCX = 2;
    // create an array of CrossSection pointers
  eCX = new CrossSection*[neCX];
    // load up the array of CrossSection pointers
  eCX[1] = new FuncCrossSection(sigmaElastic,  3,  0., CrossSection::ELASTIC);
  eCX[0] = new FuncCrossSection(sigmaIz, 3, threshold, CrossSection::IONIZATION);
}

/* Method that calculates the elastic scattering cross-section
 * in m^2 as a function of the electron impact energy in eV.
 *
 * Cross section is a fit to LLNL's EEDL calculations at high
 *   energy and, for low energy, to calculations/data from
 *   Bray, Fursa and McCarthy, Phys. Rev. A 47, 1101 (1993).
 *
 * Units are eV for energy and m^2 for cross-sections.
 * NOTE: units are NOT cm^2 and NOT barns!
 */
Scalar LiMCC::sigmaElastic(Scalar energy)
{
  Scalar sigmaMin = 6.8e-24;
  Scalar a1 = 1.9e-6;
  Scalar a2 = 8.1e-3;
  Scalar a3 = 2.4e-4;

  Scalar gMin1 = energy / ELECTRON_MASS_EV;
  Scalar betaSq = (gMin1+2.)*gMin1/(gMin1+1.)/(gMin1+1.);

  Scalar crossSection;
  crossSection  = sigmaMin * (a1 + a2*betaSq + betaSq*betaSq);
  crossSection /= pow(a3+betaSq, 3.);

	//  cout << "energy sigmaElastic = " << energy << "; " << crossSection << endl;
  return crossSection;
}


/*  Method for calculating the e-impact ionization cross section
 *  for lithium.
 *
 *  The nonrelativistic model uses the fitting function from Eq. 4
 *  of the 1981 J.Q.S.R.T. paper by S.M. Younger.
 *
 *  The relativistic model uses a parametric function adapted from
 *  Younger, but modified to capture the relativistic rise.  The
 *  parameters were determined by least squares fit to a) Younger's
 *  low-energy values and b) the computational results of the LLNL
 *  report on the EEDL (electron evauated data library), by Perkins,
 *  Cullen and Seltzer (1991).
 *
 *  Units are eV for energy and m^2 for cross sections.
 *  NOTE:  cross section is NOT in cm^2 and NOT in barns!
 *
 *  First added March 2, 2000, by David Bruhwiler
 */
Scalar LiMCC::sigmaIz(Scalar energy) {

  	// This is the value that gets returned.
  Scalar crossSection;

    // normalized impact energy
    // -- assumed that both energy and threshold are in eV
  Scalar U = energy / threshold;
  if ( U <= 1. ) return 0.;

  	// This is needed for both relativistic and NR models
  Scalar oneMinusUinv = 1. - 1. / U;

  // nonrelativistic model from Younger's paper
  if (staticRelMCC == 0) {
    Scalar logU = log(U);
    crossSection = youngerParamA * oneMinusUinv
                 - youngerParamB * oneMinusUinv * oneMinusUinv
                 + youngerParamC * logU
                 - youngerParamD * logU / U;
    crossSection /= (energy * threshold);
	}

    // relativistic model -- fit to Younger and also EEDL calc.'s
  else {
    Scalar gMin1 = energy / ELECTRON_MASS_EV;
    Scalar gamma = gMin1 + 1.;
    Scalar beta  = sqrt((gamma+1.)*gMin1)/gamma;

    Scalar oneMinusUinvSq = 1.-1./(U*U);

    // "nonrelativistic" energy, scaled to threshold
    Scalar Unr = 0.5 * ELECTRON_MASS_EV * beta * beta / threshold;
    Scalar logUnr = log(Unr);

    // Now we apply the empirical function to actually calculate
    //   the cross section:
    crossSection = fitParamA1*oneMinusUinv/U
                 + fitParamA2*logUnr/(U*U)
                 + fitParamA3*logUnr/Unr
                 + fitParamA3*oneMinusUinvSq
                   * ( log( fitParamA4*pow(gamma,fitParamN) ) / Unr
 									 - 2.*threshold/ELECTRON_MASS_EV );

		//    cout << "energy sigmaIz = " << energy << "; " << crossSection << endl;
  }

	// Return the calculated cross section (units = m^2)
  return crossSection;
}



//-------------------------------------------------------------
// NMCC: Nitrogen MCC package  -- added 09/03/2000 by Bruhwiler

NMCC::NMCC(Scalar p, Scalar temp, Species* eSpecies, Species* iSpecies, 
           SpeciesList& sList, Scalar dt, int ionzFraction,
           Scalar ecxFactor, Scalar icxFactor, int relativisticFlag,
           SpatialRegion* _SR, Scalar _Min1MKS, Scalar _Max1MKS,
           Scalar _Min2MKS, Scalar _Max2MKS,
           Scalar _delayTime, Scalar _stopTime, const ostring &analyticF,
           const int& discardDumpFileNGDDataFlag)
          : MCCPackage(N, (ostring)"N", p, temp, eSpecies, iSpecies, sList, dt,
                       ionzFraction, ecxFactor, icxFactor, _SR,
                       _Min1MKS, _Max1MKS, _Min2MKS, _Max2MKS,
                       _delayTime, _stopTime, analyticF, discardDumpFileNGDDataFlag) {

  relativisticMCC = relativisticFlag;
  staticRelMCC    = relativisticFlag;
  addCrossSections();
  init(sList, dt);
}


// Define static data members for the LiMCC class.
  Scalar NMCC::threshold     = 11.48;
  Scalar NMCC::fitParamA1    = -1.6112e-19;
  Scalar NMCC::fitParamA2    = -3.6835e-19;
  Scalar NMCC::fitParamA3    =  8.7133e-20;
  Scalar NMCC::fitParamA4    = 15.;
  Scalar NMCC::fitParamN     = 0.28;
  int    NMCC::staticRelMCC  = 1;

//
void NMCC::addCrossSections() {
    // Specify number of electron collision models that are implemented --
  	//   at present, we have impact ionization and elastic scattering.
  neCX = 2;
    // create an array of CrossSection pointers
  eCX = new CrossSection*[neCX];
    // load up the array of CrossSection pointers
  eCX[1] = new FuncCrossSection(sigmaElastic,  7,  0., CrossSection::ELASTIC);
  eCX[0] = new FuncCrossSection(sigmaIz, 7, threshold, CrossSection::IONIZATION);
}

/* Method that calculates the elastic scattering cross-section
 * in m^2 as a function of the electron impact energy in eV.
 *
 * Cross section is a fit to LLNL's EEDL calculations at high
 *   energy and, for low energy, to calculations/data from
 *   Bray, Fursa and McCarthy, Phys. Rev. A 47, 1101 (1993).
 *
 * Units are eV for energy and m^2 for cross-sections.
 * NOTE: units are NOT cm^2 and NOT barns!
 */
Scalar NMCC::sigmaElastic(Scalar energy)
{
  Scalar sigmaMin = 2.035e-23;
  Scalar u1 = 1.759e-8;
  Scalar u2 = 6.460e-5;
  Scalar l1 = 8.021e-5;

  Scalar gMin1 = energy / ELECTRON_MASS_EV;
  Scalar betaSq = (gMin1+2.)*gMin1/(gMin1+1.)/(gMin1+1.);

  Scalar crossSection;
  crossSection  = sigmaMin * (u1 + u2*betaSq + betaSq*betaSq);
  crossSection /= pow(l1+betaSq, 3.);

	//  cout << "energy sigmaElastic = " << energy << "; " << crossSection << endl;
  return crossSection;
}


/*  Method for calculating the e-impact ionization cross section
 *  for nitrogen.
 *
 *  This relativistic model uses a parametric function adapted from
 *  Younger, but modified to capture the relativistic rise.  The
 *  parameters were determined by least squares fit to the computational
 *  results of the LLNL report on the EEDL (electron evauated data
 *  library), by Perkins, Cullen and Seltzer (1991).
 *
 *  Units are eV for energy and m^2 for cross sections.
 *  NOTE:  cross section is NOT in cm^2 and NOT in barns!
 *
 *  First added September 3, 2000, by David Bruhwiler
 */
Scalar NMCC::sigmaIz(Scalar energy) {

  	// This is the value that gets returned.
  Scalar crossSection;

    // normalized impact energy
    // -- assumed that both energy and threshold are in eV
  Scalar U = energy / threshold;
  if ( U <= 1. ) return 0.;

  Scalar gMin1 = energy / ELECTRON_MASS_EV;
  Scalar gamma = gMin1 + 1.;
  Scalar beta  = sqrt((gamma+1.)*gMin1)/gamma;

  Scalar oneMinusUinv = 1. - 1. / U;
  Scalar oneMinusUinvSq = 1.-1./(U*U);

  // "nonrelativistic" energy, scaled to threshold
  Scalar Unr = 0.5 * ELECTRON_MASS_EV * beta * beta / threshold;
  Scalar logUnr = log(Unr);

  // Now we apply the empirical function to actually calculate
  //   the cross section:
  crossSection = fitParamA1*oneMinusUinv/U
               + fitParamA2*logUnr/(U*U)
               + fitParamA3*logUnr/Unr
               + fitParamA3*oneMinusUinvSq
                 * ( log( fitParamA4*pow(gamma,fitParamN) ) / Unr
								 - 2.*threshold/ELECTRON_MASS_EV );

	//  cout << "energy sigmaIz = " << energy << "; " << crossSection << endl;

	// Return the calculated cross section (units = m^2)
  return crossSection;
}


//
void MCCPackage::BoltzmannInit()
{
  
  const int NumberBins=100;
  //	const int IS; //IntegrationSteps
  const Scalar Tmin = .05;  // this needs to match min temp in boltzman.cpp
  const Scalar Tmax = 10;  // This is hopefull the hightest temperture.
  const Scalar EnergyStep=.5;
  
  Scalar T;
  G = new Scalar[NumberBins];
  E = new Scalar[NumberBins];
  for (int n=0; n<NumberBins; n++){
    G[n] = 0;
    E[n] = 0;
  }
  
  Scalar temp = 0;
  
  for (int i=0; i<neCX; i++)
    switch (eCX[i]->get_type())
      {
      case CrossSection::IONIZATION: {
        for (int n=0; n<NumberBins; n++){
          T = n*(Tmax-Tmin)/(NumberBins-1)+Tmin;
          for (Scalar energy=EnergyStep/2; energy<MAX_ENERGY; energy+=EnergyStep){
            if (energy/T<50)
              temp += energy*exp(-energy/T)*eCX[i]->sigma(energy);
          }
          temp *= gasDensity*sqrt(2/(T*M_PI*eSpecies->get_m()))*2/T;
          G[n] += temp;
          E[n] -= temp*eCX[i]->get_threshold();
        }
        break;
      }
      case CrossSection::ELASTIC:
        {
          for (int n=0; n<NumberBins; n++){
            T = n*(Tmax-Tmin)/(NumberBins-1)+Tmin;
            for (Scalar energy=EnergyStep/2; energy<MAX_ENERGY; energy+=EnergyStep){
              if (energy/T<50)
                temp += energy*exp(-energy/T)*eCX[i]->sigma(energy);
            }
            temp *= gasDensity*sqrt(2/(T*M_PI*eSpecies->get_m()))*2/T;
            E[n] -= temp*eCX[i]->get_threshold();
          }
          break;
        }
      case CrossSection::EXCITATION:
        {
          for (int n=0; n<NumberBins; n++){
            T = n*(Tmax-Tmin)/(NumberBins-1)+Tmin;
            for (Scalar energy=EnergyStep/2; energy<MAX_ENERGY; energy+=EnergyStep){
              if (energy/T<50)