mcc.cpp

pic 模拟程序！面向对象

    int jNGD = static_cast<int>(xTemp.e2());
    Scalar cellVol = grid->cellVolume(iNGD, jNGD);
    
    Scalar numNeutralAtoms;
    try{
      numNeutralAtoms =  (ptrNGD->getNGD(iNGD,jNGD)) * cellVol; 
    }
    catch(Oops& oops2){
      oops2.prepend("MCCPackage::ionization: Error: \n");//OK
      throw oops2;
    }
    
    // Calculate number of physical ions created
    Scalar numIons = pg->get_np2c(index);
    // Check that numNeutral > numIons 
    if ( numNeutralAtoms > numIons ) { 
      // create the electron species computational particle with 
      // number of physical electrons equal to numIons
      pList->add(new Particle(pg->get_x(index), uEjected, eSpecies, 
                              numIons, pg->isVariableWeight() ) );
      // continue with the ionization and reset the 
      // Neutral Gas Density accordingly.
      Scalar stmp = (numNeutralAtoms - numIons)/cellVol; 
      ptrNGD->set(iNGD, jNGD,  stmp);
      
      // create n=ionzFraction ions from isotropic gas distribution
      for (int i=0; i<ionzFraction; i++) {
	uIon = gasDistribution->get_V();
	pList->add(new Particle(pg->get_x(index), uIon, iSpecies,
				pg->get_np2c(index)/ionzFraction, pg->isVariableWeight()));
      }
    } else { 
      // this is the case of numIons >= numNeutralAtoms
      // if so, set numIons = numNeutralAtoms
      //        set density to zero
      //        replace pg->getnp2c(index) with numNeutralAtoms and
      //        replace pg->isVariableWeight() with "true" when 
      //        creating the Particle(...).

      // create the electron species particle
      pList->add(new Particle(pg->get_x(index), uEjected, eSpecies, 
                              numNeutralAtoms, true) );
                              
      ptrNGD->set(iNGD, jNGD, 0.0);
      
      // create n=ionzFraction ions from isotropic gas distribution
      for (int i=0; i<ionzFraction; i++) {
      uIon = gasDistribution->get_V();
      pList->add(new Particle(pg->get_x(index), uIon, iSpecies,
			      numNeutralAtoms/ionzFraction, true));
      }
    }

    // Now generate the new normalized momentum vector for
    //   the scattered primary.
    Vector3 uScatter;
    uScatter = newMomentum(u, newEnergy, thetaP, phiP);

		//    cout << " initial gamma*v = " << u.magnitude()        << endl;
		//    cout << " final   gamma*v = " << uScatter.magnitude() << endl;
		//    cout << " difference is:    " << uScatter.magnitude()-u.magnitude() << endl;
		//    cout << endl << endl;

    u = uScatter;
		//    v = u.magnitude();
  }
}

void MCCPackage::ionElastic(CrossSection& cx)
{

  Scalar up1, up2, up3, mag;
  Scalar r11, r12, r13, r21, r22, r23, r31, r32, r33;
			
  Scalar coschi= sqrt(frand());
  Scalar sinchi = sqrt(1. -coschi*coschi);
	
  Scalar phi1  = TWOPI*frand();
  Scalar cosphi = cos(phi1);
  Scalar sinphi = sin(phi1);
	
  r13 = u.e1()/v;
  r23 = u.e2()/v;
  r33 = u.e3()/v;
	
  if(r33 == 1.0) { up1= 0;  up2= 1;  up3= 0; }
  else           { up1= 0;  up2= 0;  up3= 1; }
  
  r12 = r23*up3 -r33*up2;
  r22 = r33*up1 -r13*up3;
  r32 = r13*up2 -r23*up1;
  mag = sqrt(r12*r12 + r22*r22 + r32*r32) + 1E-30;
	
  r12/= mag;
  r22/= mag;
  r32/= mag;
	
  r11 = r22*r33 -r32*r23;
  r21 = r32*r13 -r12*r33;
  r31 = r12*r23 -r22*r13;
	
  u.set_e1(v*coschi*(r11*sinchi*cosphi +r12*sinchi*sinphi +r13*coschi));
  u.set_e2(v*coschi*(r21*sinchi*cosphi +r22*sinchi*sinphi +r23*coschi)); 
  u.set_e3(v*coschi*(r31*sinchi*cosphi +r32*sinchi*sinphi +r33*coschi));
}

void MCCPackage::chargeExchange(CrossSection& cx)
{
  u.set_e1(0);
  u.set_e2(0);
  u.set_e3(0);
}


/* Method for calculating energy imparted to a secondary (ejected)
 *   electron, using a parametric model that captures relativistic
 *   physics for large impact energy and agrees with data at low energy.
 *
 * The calculations underlying this method are contained in the
 *   Tech-X notes, "RNG Calculations for Scattering in XOOPIC," by
 *   David Bruhwiler (August 28, 2000).
 *
 * The parametric model is a relativistic generalization of the
 *   work by Rudd et al., Phys. Rev. A 44, 1644 (1991) and Phys.
 *   Rev. A 47, 1866 (1993).  The relativistic generalization is
 *   described in the Tech-X notes, "Rudd's Differential Cross
 *   Section for Electron Impact Ionization," by David Bruhwiler,
 *   (July 17, 2000).
 */
Scalar MCCPackage::ejectedEnergy(const CrossSection& cx,
                                 const Scalar& impactEnergy) {

	// First, we need to calculate a random # between 0 and A_>,
  //   where A_> = f1(T) * (T-I) / f2(T) / (T+I), with I the threshold
	//   ionization energy and T the kinetic energy of the impacting
	//   electron "impactEnergy".

  Scalar I = cx.get_threshold();
  if (impactEnergy <= I) return 0.;

  Scalar tPlusMC    = impactEnergy + ELECTRON_MASS_EV;
	Scalar twoTplusMC = impactEnergy + tPlusMC;
  Scalar tPlusI     = impactEnergy + I;
  Scalar tMinusI    = impactEnergy - I;

  Scalar invTplusMCsq = 1. / ( tPlusMC * tPlusMC );
  Scalar tPlusIsq     = tPlusI * tPlusI;
  Scalar iOverT       = I / impactEnergy;


  Scalar funcT1 = 14./3. + .25 * tPlusIsq * invTplusMCsq
                - ELECTRON_MASS_EV * twoTplusMC * iOverT * invTplusMCsq;

  Scalar funcT2 = 5./3. - iOverT - 2.*iOverT*iOverT/3.
                + .5 * I * tMinusI * invTplusMCsq
                + ELECTRON_MASS_EV*twoTplusMC*I*invTplusMCsq*log(iOverT)/tPlusI;

  Scalar aGreaterThan = funcT1 * tMinusI / funcT2 / tPlusI;

	// We may have to try several times before the "rejection method"
	//   lets us keep our ejected energy test value, wTest.
	//   Here we define a bool for the while construct and also
	//   declare some floating point variables needed inside.

  int needToTryAgain = 1;
  Scalar randomFW, wTest, wPlusI, wPlusIsq;
  Scalar invTminusW, invTminusWsq, invTminusW3;
  Scalar probabilityRatio;

  while (needToTryAgain==1) {

  	// The random number, called F(W) in the notes and randomFW here,
  	//   is the antiderivative of a phony_but_simple probability, which
  	//   is always >= the correct_but_messy probability.

		//    do {randomFW=frand();} while (randomFW<=0. || randomFW>=1.);
    randomFW = frand() * aGreaterThan;

  	// wTest is the inverse of F(W)
    wTest = I*funcT2*randomFW/(funcT1-funcT2*randomFW);

  	// Because we are not working directly with the distribution function,
    //   we must use the "rejection method".  This involves generating
    //   another random number and seeing whether it is > the ratio of
  	//   the true probability over the phony_but_simple probability.

    wPlusI       = wTest + I;
    wPlusIsq     = wPlusI * wPlusI;
    invTminusW   = 1./(impactEnergy-wTest);
    invTminusWsq = invTminusW * invTminusW;
    invTminusW3  = invTminusW * invTminusWsq;

    probabilityRatio = ( 1. + 4.*I/wPlusI/3. + wPlusIsq*invTminusWsq
                            + 4.*I*wPlusIsq*invTminusW3/3.
                         - ELECTRON_MASS_EV*twoTplusMC*wPlusI*invTminusW*invTplusMCsq
  											 + wPlusIsq*invTplusMCsq ) / funcT1;

    if ( probabilityRatio >= (Scalar)frand() ) needToTryAgain = 0;
		//    else cout << "** wTest was rejected..." << endl;
	}  // end of while construct
	//  cout << "** wTest was accepted!" << endl << endl;

  // wTest is the energy of the ejected particle;  the scattered primary
	//   has energy T-(W+I)

	//  cout << endl;
	//  cout << " impact  energy = " << impactEnergy << endl;
	//	cout << " ejected energy = " << wTest << endl;
  return wTest;
}


/* Method for calculating the scattering angles (theta, phi) with
 *   respect to the initial direction of the primary electron, for
 *   both the primary and secondary electrons in an impact-ionizaiton
 *   event.
 *
 * Unchanged imput variables are eImpact and eSecondary.
 * Empty variables passed in by reference to hold the desired results
 *   are:  thetaPrimary, thetaSecondary, phiPrimary and phiSecondary.
 *
 * The calculations underlying this method are contained in the
 *   Tech-X notes, "RNG Calculations for Scattering in XOOPIC," by
 *   David Bruhwiler (August 28, 2000).
 *
 * The parametric model is discussed in the Tech-X notes, "Rudd's
 *   Differential Cross Section for Electron Impact Ionization,"
 *   by David Bruhwiler, (July 17, 2000).
 */
void MCCPackage::primarySecondaryAngles(const CrossSection& cx,
            const Scalar& eImpact,  const Scalar& eSecondary,
                  Scalar& thetaPrimary,   Scalar& phiPrimary,
                  Scalar& thetaSecondary, Scalar& phiSecondary) {

    // Bail out, if impact energy < ionization threshold energy.

  Scalar I = cx.get_threshold();
  if (eImpact <= I) {
    thetaPrimary   = 0.;
    thetaSecondary = 0.;

    phiPrimary   = 0.;
    phiSecondary = 0.;

    return;
	}

    // First, we generate a Monte Carlo value for thetaPrimary --
    // We begin with the definition of some temporaries:

  Scalar ePrimary = eImpact - eSecondary - I;

  Scalar w         = ePrimary;
  Scalar wPlusI    = w       + I;
  Scalar wPlusI2MC = wPlusI  + 2.*ELECTRON_MASS_EV;
  Scalar tPlus2MC  = eImpact + 2.*ELECTRON_MASS_EV;

  Scalar alpha;
  alpha  = ELECTRON_MASS_EV / ( eImpact + ELECTRON_MASS_EV );
  alpha *= 0.6 * alpha;

  Scalar g2_T_W = sqrt( wPlusI * tPlus2MC / wPlusI2MC / eImpact );
  Scalar g3_T_W = alpha * sqrt( I * (1.-wPlusI/eImpact) / w );

	//  cout << "  g2-1 = " << g2_T_W - 1. << endl;
	//  cout << "  g3   = " << g3_T_W << endl;

  Scalar fTheta = frand();

  thetaPrimary = acos( g2_T_W + g3_T_W * tan(
                  (1.-fTheta) * atan( (1.-g2_T_W) / g3_T_W )
                     - fTheta * atan( (1.+g2_T_W) / g3_T_W ) ) );

	//  cout << " thetaPrimary = " << thetaPrimary << endl;

  	// Second, we generate a Monte Carlo value for phiPrimary --
  phiPrimary = TWOPI * frand();


    // Third, we generate a Monte Carlo value for thetaSecondary --
  w         = eSecondary;
  wPlusI    = w + I;
  wPlusI2MC = wPlusI + 2.*ELECTRON_MASS_EV;

  g2_T_W = sqrt( wPlusI * tPlus2MC / wPlusI2MC / eImpact );
  g3_T_W = alpha * sqrt( I * (1.-wPlusI/eImpact) / w );

  fTheta = frand();

  thetaSecondary = acos( g2_T_W + g3_T_W * tan(
                    (1.-fTheta) * atan2( (Scalar)(1.-g2_T_W) , (Scalar)g3_T_W )
                       - fTheta * atan2( (Scalar)(1.+g2_T_W) , (Scalar)g3_T_W ) ) );

  	// Fourth, we generate a Monte Carlo value for phiSecondary --
  	// In the low-energy limit, phiPrimary and phiSecondary are
  	//   uncorrelated.  In the high-energy limit, they are 100%
  	//   correlated.  Our approach is to implement a sort of sliding
  	//   scale, in which the correlation increases as ePrimary and
    //   eSecondary both increase.

  Scalar eCritical = 10. * I;
  Scalar a_T_W = 0.;
  if ( eSecondary>eCritical && ePrimary>eCritical ) {
    a_T_W = sqrt((eSecondary-eCritical)*(ePrimary-eCritical))/eCritical;
		//    cout << "  a_T_W = " << a_T_W << endl;
	}
  Scalar delta = TWOPI/(1.+a_T_W);
	//  cout << "  delta/(2*pi) = " << delta/TWOPI << endl;

  phiSecondary = phiPrimary + M_PI + (frand()-.5)*delta;

	//  cout << "  phiSecondary-(phiPrimary+PI) = "
	//       << phiSecondary-phiPrimary-M_PI << endl;

	return;
}


/* Method for calculating the scattering angles (theta, phi) with
 *   respect to the initial direction of the electron, for elastic
 *   electron-neutral scattering.
 *
 * Unchanged input variables are eImpact and the cross section cx.
 * Empty variables passed in by reference to hold the desired results
 *   are:  thetaScatter and phiScatter.
 *
 * The calculations underlying this method are contained in the
 *   Tech-X notes, "RNG Calculations for Scattering in XOOPIC," by
 *   David Bruhwiler (August 28, 2000).
 *
 * The parametric model is discussed in the Tech-X notes, "Elastic
 *   and Inelastic Scattering at all Energies," by David Bruhwiler,
 *   (July 26, 2000).
 */
void MCCPackage::elasticScatteringAngles(
                    const Scalar& eImpact, const CrossSection& cx,
                          Scalar& thetaScatter, Scalar& phiScatter) {

    // First, we generate a Monte Carlo value for thetaScatter --
    // We begin with the definition of some temporaries:

  Scalar atomicNum = cx.getAtomicNumber();

  Scalar gMinus1  = eImpact / ELECTRON_MASS_EV;
  Scalar gamma    = gMinus1 + 1.;
  Scalar beta     = sqrt((gamma+1.)*gMinus1)/gamma;
  Scalar betaMin  = 1.8e-3;

  Scalar thetaMin = pow( atomicNum, 1./3. ) * (1.+betaMin)
                  / ( 192. * gamma * (beta+betaMin) );

  Scalar thetaMinSq = thetaMin * thetaMin;
  Scalar fTheta = frand();
  thetaScatter = acos( (1.+thetaMinSq)
                         - thetaMinSq * (2.+thetaMinSq)
                         / (2.*(1-fTheta) + thetaMinSq) );

	//  if (eImpact>1e10) {
		//    cout << endl;
		//    cout << "impact energy = " << eImpact << endl;
		//    cout << "thetaMin      = " << thetaMin << endl;
		//    cout << "thetaElastic  = " << thetaScatter << endl;
	//	}

  	// Next, we generate a Monte Carlo value for phiScatter --
  phiScatter = TWOPI * frand();

	return;
}


/* Given an initial normalized momentum vel=p/(m*gamma) (m/s), a new
 *   energy, and the scattering angles (wrt to initial direction
 *   defined by u) theta and phi, this method calculates and
 *   returns a new momentum.
 *
 * Note that the magnitude of u is not important -- only the
 *   direction is used.  Thus, this method can be used for
 *   elastic scattering, for the primary in an ionization event,
 *   and also for the ejected secondary.
 *
 * All of the input variables are unchanged.
 * The new momentum is returned.