mcc.cpp

pic 模拟程序！面向对象

 *
 * The calculations underlying this method are contained in the
 *   Tech-X notes, "Scattered Velocity Calculation for XOOPIC,"
 *   by David Bruhwiler (August 31, 2000).
 */
Vector3 MCCPackage::newMomentum(const Vector3 U_initial,
                                const Scalar& newEnergy,
             										const Scalar& theta,
                                const Scalar& phi) {

  Scalar u10 = U_initial.e1();
  Scalar u20 = U_initial.e2();
  Scalar u30 = U_initial.e3();
  Scalar u0  = sqrt( u10*u10 + u20*u20 + u30*u30 );

  Scalar cosT0 = u30 / u0;
  Scalar sinT0 = sqrt( 1. - cosT0*cosT0 );

  Scalar cosPhi0;
  Scalar root = sqrt(u0*u0-u30*u30);
  if ( root <= fabs(u10) ) cosPhi0 = 1;
  else cosPhi0 = u10 / root;

  Scalar sinPhi0 = sin( acos( cosPhi0 ) );

  Scalar cosT = cos(theta);
  Scalar sinT = sin(theta);

  Scalar cosPhi = cos(phi);
  Scalar sinPhi = sin(phi);

  Scalar gMin1 = newEnergy / ELECTRON_MASS_EV;
  Scalar u0new = SPEED_OF_LIGHT * sqrt( (gMin1+2) * gMin1 );  // =beta*gamma*c

  Scalar u10new = u0new * ( cosPhi*sinT*cosPhi0*cosT0
                          - sinPhi*sinT*sinPhi0
													+        cosT*cosPhi0*sinT0 );

  Scalar u20new = u0new * ( cosPhi*sinT*sinPhi0*cosT0
                          + sinPhi*sinT*cosPhi0
													+        cosT*sinPhi0*sinT0 );

  Scalar u30new = u0new * (-cosPhi*sinT*sinT0 + cosT*cosT0);

  Vector3 uNew(u10new,u20new,u30new);
  return uNew;
}



void MCCPackage::newVelocity(Scalar energy, Scalar v, Vector3& u, int eFlag)
{
	Scalar	up1, up2, up3;
	Scalar	coschi;
	Scalar	sinchi;
	if (energy < 1e-30) coschi = 1.0;
	else coschi = 1 + 2.0*(1.0 - pow(energy + 1.0, frand()))/energy;
	sinchi = sqrt(1 - coschi*coschi);
	Scalar	phi = TWOPI*frand();
	Scalar	cosphi = cos(phi);
	Scalar	sinphi = sin(phi);

	// transfer momentum:

//	if (eFlag) v *= sqrt(1 - 2*eSpecies->get_m()/iSpecies->get_m()*(1 - coschi));
	// doesn't work when m/M>.5

	Scalar me = eSpecies->get_m();
	Scalar mi = iSpecies->get_m();
	v*= sqrt(1 - 2*me*mi/(me+mi)/(me+mi)*(1-coschi));

	Scalar	r13 = u.e1();
	Scalar	r23 = u.e2();
	Scalar	r33 = u.e3();

	if (r33 == 1)	{up1 = 0; up2 = 1; up3 = 0;}
	else				{up1 = 0; up2 = 0; up3 = 1;}

	Scalar	r12 = r23*up3 - r33*up2;
	Scalar	r22 = r33*up1 - r13*up3;
	Scalar	r32 = r13*up2 -r23*up1;
	Scalar	mag = sqrt(r12*r12 + r22*r22 + r32*r32);
  
	r12 /= mag;
	r22 /= mag;
	r32 /= mag;
  
	Scalar	r11 = r22*r33 - r32*r23;
	Scalar	r21 = r32*r13 - r12*r33;
	Scalar	r31 = r12*r23 - r22*r13;

	u.set_e1(v*(r11*sinchi*cosphi + r12*sinchi*sinphi + r13*coschi));
	u.set_e2(v*(r21*sinchi*cosphi + r22*sinchi*sinphi + r23*coschi)); 
	u.set_e3(v*(r31*sinchi*cosphi + r32*sinchi*sinphi + r33*coschi));
}


//-------------------------------------------------------------
// ArMCC: argon MCC package

ArMCC::ArMCC(Scalar p, Scalar temp, Species* eSpecies, Species* iSpecies,
      SpeciesList& sList, Scalar dt, int ionzFraction, Scalar ecxFactor,
      Scalar icxFactor, SpatialRegion* _SR,
      Scalar _Min1MKS, Scalar _Max1MKS, Scalar _Min2MKS, Scalar _Max2MKS,
      Scalar _delayTime, Scalar _stopTime, const ostring &analyticF, 
      const int& discardDumpFileNGDDataFlag)
 : MCCPackage(Ar, (ostring) "Ar", p, temp, eSpecies, iSpecies, sList, dt, ionzFraction,
      ecxFactor, icxFactor, _SR, _Min1MKS, _Max1MKS, _Min2MKS, _Max2MKS,
      _delayTime, _stopTime, analyticF, discardDumpFileNGDDataFlag)
{
  addCrossSections();
  init(sList, dt);
}

void ArMCC::addCrossSections()
{
  neCX = 3; // electron model
  eCX = new CrossSection*[neCX];
  eCX[0] = new FuncCrossSection(sigmaElastic, 18, 0.0, CrossSection::ELASTIC);
  eCX[1] = new FuncCrossSection(sigmaExc, 18, 12.0, CrossSection::EXCITATION);
  eCX[2] = new FuncCrossSection(sigmaIz, 18, 15.76, CrossSection::IONIZATION);
  // eCX[0] = new eCrossSection(0.3, 15, 20, 1.2e-19, CrossSection::ELASTIC);
  //  eCX[1] = new eCrossSection(11.55, 30, 100, 7e-21, CrossSection::EXCITATION);
  //	ionizationCX = new eCrossSection(15.76, 30, 100, 3e-20, IONIZATION);


  niCX = 2; // ion model
  iCX = new CrossSection*[niCX];
  
  // The original cross-sections.  Cause problems at lower energies. 
  // Commented out by: Jeff Hammel, 11/4/02
  //  iCX[0] = new iCrossSection(1.8e-19, 4e-19, CrossSection::ION_ELASTIC);
  //  iCX[1] = new iCrossSection(2e-19, 5.5e-19, CrossSection::CHARGE_EXCHANGE);
  
  iCX[0] = new FuncCrossSection(sigmaIelastic, 18, 0.0, CrossSection::ION_ELASTIC);
  iCX[1] = new FuncCrossSection(sigmaCX, 18, 0.0, CrossSection::CHARGE_EXCHANGE);
}

Scalar ArMCC::sigmaElastic(Scalar energy)
{
    
  if(energy < 1.0) 
    {
      if(energy < 0.2) { return 1./pow(10.0, 19.0 +energy/.11); }
      else 
	{ 
	  return 9.07e-19*pow(energy, 1.55)*pow(energy+70.0, 1.10)/pow(14.+energy, 3.25); 
	}
    }
  else 
    {
      return 9.07e-19*pow(energy, 1.55)*pow(energy+70.0, 1.10)/pow(14.+energy, 3.25);
    }

}

Scalar ArMCC::sigmaExc(Scalar energy)
{
 
  if(energy < 12.0) { return 0.; }  
  return (3.85116e-19*log(energy/3.4015) -4.85227e-19)/energy;
}

Scalar ArMCC::sigmaCX(Scalar energy)
{
  if(energy > 4.0) return(2.0e-19 +5.5e-19/sqrt(energy));
  return(-2.95e-19*sqrt(energy) +10.65e-19);
}
Scalar ArMCC::sigmaIelastic(Scalar energy)
{
  if(energy > 4.0) return(1.8e-19 +4.0e-19/sqrt(energy));
  return(-2.0e-19*sqrt(energy) +7.8e-19);
}


//-------------------------------------------------------------------
// Argon ionization cross section (Rapp & Golden) - energies in eV

Scalar ArMCC::sigmaIz(Scalar energy)
{
   if (energy < 15.76) return 0;
	else if (energy < 79) 
	  return 1.7155e-18*(energy-15.76)/(energy*energy)*log(0.0634*energy);
	else return 2.648e-18*(energy-15.76)/(energy*energy)*log(0.0344*energy);
}

//-------------------------------------------------------------
// NeMCC: neon MCC package

NeMCC::NeMCC(Scalar p, Scalar temp, Species* eSpecies, Species* iSpecies, 
     SpeciesList& sList, Scalar dt, int ionzFraction, Scalar ecxFactor,
     Scalar icxFactor, SpatialRegion* _SR,
     Scalar _Min1MKS, Scalar _Max1MKS, Scalar _Min2MKS, Scalar _Max2MKS,
     Scalar _delayTime, Scalar _stopTime, const ostring &analyticF, 
     const int& discardDumpFileNGDDataFlag)
 : MCCPackage(Ne, (ostring)"Ne", p, temp, eSpecies, iSpecies, sList, dt, ionzFraction,
     ecxFactor, icxFactor, _SR, _Min1MKS, _Max1MKS, _Min2MKS, _Max2MKS,
     _delayTime, _stopTime, analyticF, discardDumpFileNGDDataFlag)
{
  addCrossSections();
  init(sList, dt);
}

void NeMCC::addCrossSections()
{
  neCX = 3; // electron model
  eCX = new CrossSection*[neCX];
  eCX[0] = new FuncCrossSection(sigmaElastic, 10, 0.,  CrossSection::ELASTIC);
  eCX[1] = new FuncCrossSection(sigmaExc, 10, 15.8489, CrossSection::EXCITATION);
  eCX[2] = new FuncCrossSection(sigmaIz, 10,  20.893,  CrossSection::IONIZATION);

  niCX = 1; // ion model
  iCX = new CrossSection*[niCX];
  iCX[0] = new iCrossSection(0, 5.98e-20, CrossSection::CHARGE_EXCHANGE);
}

Scalar NeMCC::sigmaElastic(Scalar energy)
{
  if (energy < 2.63027) return 1.27808e-20 + energy*(4.22329e-21 - 7.31854e-22*energy);
  else if (energy < 19.4984) return 1.66333e-20 + energy*(8.66972e-22 - 1.14153e-23*energy);
  else if (energy < 33.8844) return 2.2297e-20 + energy*(5.8066e-22 - 1.18023e-23*energy);
  else if (energy < 190.546) return 3.862e-20 + energy*(-3.19062e-22 + 8.00046e-25*energy);
  else if (energy < 501.187) return 1.37254e-20 + energy*(-4.66403e-23 + 4.6479e-26*energy);
  else if (energy < 1000) return 4.38798e-21 + energy*(-6.19168e-24 + 2.5747e-27*energy);
  else return 0;
}

Scalar NeMCC::sigmaExc(Scalar energy) // summed excitations
{
  if (energy <= 15.8489) return 0;
  else if (energy < 23.9883) return -2.37959e-22 + energy*(-5.87101e-23 + 4.6417e-24*energy);
  else if (energy < 36.3078) return -3.49005e-21 + energy*(2.76934e-22 - 3.69893e-24*energy);
  else if (energy < 54.9541) return 2.75471e-22 + energy*(6.0912e-23 - 5.99808e-25*energy);
  else if (energy < 251.189) return 2.271e-21 + energy*(-8.78552e-24 + 1.41951e-26*energy);
  else if (energy < 616.595) return 1.53078e-21 + energy*(-2.7892e-24 + 1.90737e-27*energy);
  else if (energy < 1000) return 1.05134e-21 + energy*(-1.07925e-24 + 3.83108e-28*energy);
  else return 0;
}

Scalar NeMCC::sigmaIz(Scalar energy)
{
  if (energy <= 20.893) return 0;
  else if (energy < 54.9541) return -2.34018e-21 + energy*1.10755e-22;
  else if (energy < 125.893) return -2.4556e-21 + energy*(1.42218e-22 - 5.37266e-25*energy);
  else if (energy < 190.546) return 2.62391e-21 + energy*(5.39463e-23 - 1.54319e-25*energy);
  else if (energy < 436.516) return 8.04001e-21 + energy*(-2.68853e-24 - 6.11372e-27*energy);
  else if (energy < 1000) return 9.16045e-21 + energy*(-9.59021e-24 + 3.85096e-27*energy);
  else return 0;
}

//-------------------------------------------------------------
// XeMCC: xenon MCC package

XeMCC::XeMCC(Scalar p, Scalar temp, Species* eSpecies, Species* iSpecies, 
     SpeciesList& sList, Scalar dt, int ionzFraction, Scalar ecxFactor,
     Scalar icxFactor, SpatialRegion* _SR,
     Scalar _Min1MKS, Scalar _Max1MKS, Scalar _Min2MKS, Scalar _Max2MKS,
     Scalar _delayTime, Scalar _stopTime, const ostring &analyticF,
     const int& discardDumpFileNGDDataFlag)
 : MCCPackage(Xe, (ostring)"Xe", p, temp, eSpecies, iSpecies, sList, dt, ionzFraction,
     ecxFactor, icxFactor, _SR, _Min1MKS, _Max1MKS, _Min2MKS, _Max2MKS,
     _delayTime, _stopTime, analyticF, discardDumpFileNGDDataFlag)
{
  addCrossSections();
  init(sList, dt);
}

void XeMCC::addCrossSections()
{
  neCX = 3; // electron model
  eCX = new CrossSection*[neCX];
  eCX[0] = new FuncCrossSection(sigmaElastic, 54, 0.,  CrossSection::ELASTIC);
  eCX[1] = new FuncCrossSection(sigmaExc, 54, 7.94328, CrossSection::EXCITATION);
  eCX[2] = new FuncCrossSection(sigmaIz,  54, 11.2202, CrossSection::IONIZATION);

  niCX = 1; // ion model
  iCX = new CrossSection*[niCX];
  iCX[0] = new iCrossSection(0, 4.541e-20, CrossSection::CHARGE_EXCHANGE);
}

Scalar XeMCC::sigmaElastic(Scalar energy)
{
  if (energy <= 1) return 0;
  else if (energy < 3.46737) return -2.56426e-20 + energy*(3.05131e-20 + 1.01401e-20*energy);
  else if (energy < 4.2658) return -7.02412e-19 + energy*(4.15012e-19 - 4.46826e-20*energy);
  else if (energy < 5.62341) return -3.14744e-19 + energy*(2.19764e-19 - 2.027176e-20*energy);
  else if (energy < 12.0226) return 2.23972e-19 + energy*(2.61272e-20 - 2.81077e-21*energy);
  else if (energy < 25.704) return 3.39014e-19 + energy*(-2.28038e-20 + 4.64864e-22*energy);
  else if (energy < 72.4436) return 9.50404e-20 + energy*(-1.80585e-21 + 1.24365e-23*energy);
  else if (energy < 190.546) return 4.51567e-20 + energy*(-2.78976e-22 + 6.93189e-25*energy);
  else if (energy < 537.032) return 2.59891e-20 + energy*(-5.72321e-23 + 4.79915e-26*energy);
  else if (energy < 1000) return 1.46698e-20 + energy*(-1.28135e-23 + 3.98769e-27*energy);
  else return 0;
}

Scalar XeMCC::sigmaExc(Scalar energy) // summed excitations
{
  if (energy <= 7.94328) return 0;
  else if (energy < 11.2202) return 7.71935e-20 + energy*(-1.8566e-20 + 1.117e-21*energy);
  else if (energy < 20.893) return -5.72625e-20 + energy*(8.29355e-21 - 2.02235e-22*energy);
  else if (energy < 33.8844) return 9.4784e-21 + energy*(1.36322e-21 - 2.20758e-23*energy);
  else if (energy < 190.546) return 3.66375e-20 + energy*(-1.87228e-22 + 3.94988e-25*energy);
  else if (energy < 467.735) return 2.51242e-20 + energy*(-6.45713e-23 + 6.14246e-26*energy);
  else if (energy < 1000) return 1.4849e-20 + energy*(-1.78916e-23 + 7.74732e-27*energy);
  else return 0;
}

Scalar XeMCC::sigmaIz(Scalar energy)
{
  if (energy <= 11.2202) return 0;
  else if (energy < 44.6684) return -4.78153e-20 + energy*(4.75275e-21 - 5.83101e-23*energy);
  else if (energy < 67.6083) return 4.90931e-20 + energy*(-6.99912e-23 + 1.11671e-24*energy);
  else if (energy < 102.329) return 3.73779e-20 + energy*(3.25673e-22 - 2.15897e-24*energy);
  else if (energy < 204.174) return 6.20189e-20 + energy*(-1.62295e-22 + 2.70937e-25*energy);
  else if (energy < 575.44) return 4.1e-20 - 3.49995e-24*energy;
  else if (energy < 812.831) return 2.19858e-20 + energy*(9.27581e-23 - 1.09139e-25*energy);
  else if (energy < 1000) return 8.37604e-20 + energy*(-1.01201e-22 + 3.62937e-26*energy);
  else return 0;
}

//-------------------------------------------------------------
// HMCC:  Hydrogen MCC package

HMCC::HMCC(Scalar p, Scalar temp, Species* eSpecies, Species* iSpecies, 
     SpeciesList& sList, Scalar dt, int ionzFraction, Scalar ecxFactor,
     Scalar icxFactor, SpatialRegion* _SR,
     Scalar _Min1MKS, Scalar _Max1MKS, Scalar _Min2MKS, Scalar _Max2MKS,
     Scalar _delayTime, Scalar _stopTime, const ostring &analyticF, 
     const int& discardDumpFileNGDDataFlag)
 : MCCPackage(H, (ostring)"H", p, temp, eSpecies, iSpecies, sList, dt, ionzFraction,
     ecxFactor, icxFactor, _SR, _Min1MKS, _Max1MKS, _Min2MKS, _Max2MKS,
     _delayTime, _stopTime, analyticF, discardDumpFileNGDDataFlag)
{
  addCrossSections();
  init(sList, dt);
}

void HMCC::addCrossSections()
{
  neCX = 1; // electron model
  eCX = new CrossSection*[neCX];
  eCX[0] =  new FuncCrossSection(sigmaIz, 1, 13.6, CrossSection::IONIZATION);

}

Scalar HMCC::sigmaIz(Scalar energy)
{
  if (energy <= 13.6) return 0;
  else if (energy <= 17.789) return 8.208063e-21*log(energy)-2.071614e-20;
  else if (energy <= 53.6999) return 3.751578e-21*log(energy)-7.887766e-21;
  else return (2.996231e-19*log(energy)-8.145988e-19)/energy;
}
//-------------------------------------------------------------------

//-------------------------------------------------------------
// HeMCC: Helium MCC package    -- added 10/24/02 by Hammel;
// electron cross-sections from pdp1, 
// ion cross sections from McDaniel et al, 1993  (curve fit)

HeMCC::HeMCC(Scalar p, Scalar temp, Species* eSpecies, Species* iSpecies, 
     SpeciesList& sList, Scalar dt, int ionzFraction, Scalar ecxFactor,
     Scalar icxFactor, SpatialRegion* _SR,
     Scalar _Min1MKS, Scalar _Max1MKS, Scalar _Min2MKS, Scalar _Max2MKS,
     Scalar _delayTime, Scalar _stopTime, const ostring &analyticF, 
     const int& discardDumpFileNGDDataFlag)
 : MCCPackage(H, (ostring)"He", p, temp, eSpecies, iSpecies, sList, dt, ionzFraction,
     ecxFactor, icxFactor, _SR, _Min1MKS, _Max1MKS, _Min2MKS, _Max2MKS,
     _delayTime, _stopTime, analyticF, discardDumpFileNGDDataFlag)
{
  addCrossSections();
  init(sList, dt);
}

void HeMCC::addCrossSections()
{
  neCX = 3; // electron model
  eCX = new CrossSection*[neCX];
  eCX[0] =  new FuncCrossSection(sigmaElastic, 2, 0.0, CrossSection::ELASTIC);
  eCX[1] =  new FuncCrossSection(sigmaExc, 2, 19.8, CrossSection::EXCITATION);
  eCX[2] =  new FuncCrossSection(sigmaIz, 2, 24.5, CrossSection::IONIZATION);

  niCX = 2; // ion model
  iCX = new CrossSection*[niCX];