mcc.h

pic 模拟程序！面向对象

/*
=====================================================================
mcc.h

Describes the MCC (Monte Carlo collision) class.  This class provides 
collisions by operating on ParticleGroups.

0.9	 (JohnV 02-06-95) Initial draft.
1.1  (JohnV 12-12-96) Added Rapp&Golden Argon ionization cross section
2.5  (Bruhwiler 03-01-00) Added electron impact ionization cross section for Li
2.51 (Bruhwiler 08-29-00) added relativisticMCC flag and modified the Li
        ionization cross section to include relativistic effects
2.52 (Bruhwiler 09-03-00) added relativistic elastic scattering for
        Li -- this required changes in CrossSection and FuncCrossSection
2.53 (Bruhwiler 09-26-00) added capability to limit MCC region (i.e.
        limit the neutral gas) to a subregion of the grid.
2.54 (Bruhwiler 10-04-00) added capability to bound MCC algorithm in time.

=====================================================================
*/

#ifndef	__MCC_H
#define	__MCC_H

#include "ovector.h"
#include "misc.h"
#include "particle.h"
#include "sptlrgn.h"
#include "ngd.h"

#include <oops.h>

 
class Maxwellian;
class	Species;

class CrossSection
{
  private:

  void newVelocity(Scalar energy, Scalar v, Vector3& u, int eFlag);

  public:

    enum Type { ELASTIC, EXCITATION, IONIZATION, ION_ELASTIC,
                CHARGE_EXCHANGE, TEST_TO_BOLTZMANN } type;

	/** Return the collision cross section.
	 *  Pure virtual -- must be implemented in derived classes.
	 */
  virtual Scalar sigma(Scalar energy)=0;

	/** return inelastic threshold energy
	 */
  virtual Scalar get_threshold() const {return 0.;}

	/** return atomic number
	 */
	virtual int getAtomicNumber() const {return 1;}

	/** return collision type
	 */
  Type get_type() {return type;}
};

class eCrossSection : public CrossSection
{
	Scalar	energy0;						//	(eV) threshold
	Scalar	energy1;						//	(eV) start of plateau
	Scalar	energy2;						//	(eV) begin logE/E drop 
	Scalar	sigmaMax;						//	peak cross section m^2
	Scalar	const1;							//	internal constant for energy range 1
	Scalar	const3;							//	internal constant for energy range 3
public:
	eCrossSection(Scalar e0, Scalar e1, Scalar e2, Scalar sMax, Type type);
	virtual Scalar sigma(Scalar energy);
	//Scalar	get_energy0() {return energy0;}
	virtual Scalar get_threshold() const {return energy0;}
};

class FuncCrossSection : public CrossSection
{
protected:
   Scalar threshold;
   int atomicNumber;
   Scalar (*sigmaPtr)(Scalar energy);
public:
	 FuncCrossSection(Scalar (*sigma)(Scalar), int _Z,
                    Scalar _threshold, Type _type) {
	   sigmaPtr     = sigma;
     atomicNumber = _Z;
     threshold    = _threshold;
     type         = _type;
}
	virtual Scalar sigma(Scalar energy) {return (*sigmaPtr)(energy);}
	virtual Scalar get_threshold() const {return threshold;}
	virtual int getAtomicNumber() const  {return atomicNumber;}
};


class argonIoniz : public eCrossSection
{
 public:
   argonIoniz() : eCrossSection(15.76f, 30, 100, 1, IONIZATION) {};
   virtual Scalar sigma(Scalar energy);
};

class iCrossSection : public CrossSection
{
	Scalar	a;									//	cross section coeff.
	Scalar	b;									// ditto
public:
	iCrossSection(Scalar a, Scalar b, Type type);
	virtual Scalar	sigma(Scalar energy);
};

/**
 * The BaseMCPackage (base Monte Carlo package) will contain data and methods 
 * common for the derived classes MCCPackage and LIMCPackage. 
 * dad 01/18/2001
 */
class BaseMCPackage {
 public:
  enum GAS_TYPE {Ar, Ne, Xe, H, He, Li, N, Cs} type;
  BaseMCPackage(GAS_TYPE _type, const ostring& strGasType, 
                const Scalar& temperature, 
                Species* _eSpecies, Species* _iSpecies, 
                const Scalar& pressure, SpatialRegion* _SR,
                Scalar _Min1MKS, Scalar _Max1MKS, 
                Scalar _Min2MKS, Scalar _Max2MKS,
                const ostring& _analyticF, 
                const int& discardDumpFileNGDDataFlag) throw(Oops);
  ~BaseMCPackage() {};
 protected: 
  // 
  // variable which are used by both the MCCPackage and the LIMCPackage
  // 
  // the following two lines of variables 
  // are used for confinement of gas density to subregion of grid
  // 
  Scalar Min1MKS, Max1MKS, Min2MKS, Max2MKS;
  Vector2 deltaP, p1, p2, p1Grid, p2Grid;
  
  SpatialRegion* region;  // spatial region that owns this collision model
  Grid* grid;             // the grid associated with the spatial region
  int rz;                 // =0 by default;  =1 for cylindrical geometry
  int J, K;               // maximum grid values

  /**
   * the following two pointers are used in both the collision Monte Carlo
   * and the tunneling ionization. 
   */
  ParticleList* pList;
  ParticleGroup* pg;

  /**
   * electron and ion species will be created both in the 
   * Monte Carlo collisions and the tunneling ionization
   */
  Species* eSpecies;   
  Species* iSpecies;   
  
  /**
   * data members for gas type and gas temperature 
   */
  ostring gasTypeString;  // character representation of the gas type
  Scalar temp;            // gas temp in eV

  /**
   * When the Neutral Gas Density varies, 
   * the collision probability should depend on the maximum value of
   * this density. We'll calcualte it at initialization. 
   * The maximum Neutral Gas Density "maxNGD" is equal to
   * the "gasDensity" in the uniform case, and in the case of varying
   * neutral gas density, the "gasDensity" will be initialized to the value
   * of "maxNGD". dad 12/04/2000. 
   */
 Scalar gasDensity;   // gas number density m^-3 UNIFORM GAS DENSITY CASE
 Scalar maxNGD; 
 NGD* ptrNGD; // pointer to an NGD class (NeutralGasDensity) 
 
 private:
 BaseMCPackage(const BaseMCPackage&);
 BaseMCPackage& operator=(const BaseMCPackage&);
}; 

class MCCPackage : public BaseMCPackage
{
 private:
  Maxwellian* gasDistribution;
  //	Scalar MinPot;
  Scalar* collProb;
  Scalar ecxFactor;    // factor to multiply by electron cx
  Scalar icxFactor;    // factor to multiply by ion cx
  Scalar maxSigmaVe;   // MAX(sigma*v)
  Scalar maxSigmaVi;
  Scalar* extra;       // fractional collisions from previous steps
  int ionzFraction;    // fraction of particles to create in ionization event
  Scalar	eEnergyScale;		//	m/2e (converts (m/s)^2 -> eV)
  Scalar	iEnergyScale;		//	M/2e
  void newVelocity(Scalar energy, Scalar v, Vector3& u, int eFlag);
  Scalar sumSigmaVe();    // returns MAX(sigma*ve)
  Scalar sumSigmaVi();    // returns MAX(sigma*vi)
  
  // following are used in both collision() and elastic(), excitation(), etc.
  Vector3 u;   // 3-D velocity (gamma factor stripped out)
  Scalar v;    // v is the magnitude of u
  Scalar energy;
  int index;
  virtual void addCrossSections();
  
  // used to restrict use of MCC algorithm to specified time interval
  Scalar delayTime;  // don't start MCC until this much time has elapsed
  Scalar stopTime;   // stop use of MCC after this time (if value is > 0.)
  
 protected:
  
  /** Flag indicating whether to use the original nonrelativistic
   *  MCC algorithms (= 0) or the new relativistic models (= 1).
   *
   *  The default value is 0.  Can be set to 1 via the input file.
   *
   *  At present, only the Li (N coming soon) model is relativistic.
   *  The Li MCC model can be used in nonrelativistic mode, but the
   *  default is fully relativistic.
   *
   *  -- Bruhwiler 08/29/2000
   */
  int relativisticMCC;
  
  /** 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 ejectedEnergy(const CrossSection& cx, const Scalar& impactEnergy);
  
  /** 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 input variables are impactEnergy and ejectedEnergy.
   *  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 primarySecondaryAngles(const CrossSection& cx,
                              const Scalar& impactEnergy,
                              const Scalar& ejectedEnergy,
                              Scalar& thetaPrimary,   Scalar& phiPrimary,
                              Scalar& thetaSecondary, Scalar& phiSecondary);
  
  /** 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 elasticScatteringAngles(
                               const Scalar& eImpact, const CrossSection& cx,
                               Scalar& thetaScatter, Scalar& phiScatter);

  /** Given an initial normalized momentum u=p/m=gamma*v (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.
   *
   *  The calculations underlying this method are contained in the
   *    Tech-X notes, "Scattered Velocity Calculation for XOOPIC,"
   *    by David Bruhwiler (August 31, 2000).
   */
  Vector3 newMomentum(const Vector3 U_initial, const Scalar& energy,
                      const Scalar& theta, const Scalar& phi);
  
 public:
  
  MCCPackage( GAS_TYPE type, const ostring& strGasType, 
              Scalar pressure, 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);
  virtual ~MCCPackage();
  virtual void init(SpeciesList& sList, Scalar dt);
  virtual void BoltzmannInit();
  Scalar *G,*E;
  int neCX; // number of cross sections
  int niCX; // number of cross sections
  int ntCX; 
  int nbCX;
  CrossSection** eCX;
  CrossSection** iCX;
  CrossSection** tCX;
  CrossSection** bCX;
  void collide(ParticleGroupList &pgList, ParticleList& pList)throw(Oops);
  virtual void dynamic(CrossSection& cx)throw(Oops);
  virtual void elastic(CrossSection& cx);
  virtual void excitation(CrossSection& cx);
  virtual void ionization(CrossSection& cx) throw(Oops);
  virtual void ionElastic(CrossSection& cx);
  virtual void chargeExchange(CrossSection& cx);
  
  virtual void setRelativisticMCC(int flag) {relativisticMCC = flag;}
  virtual int  getRelativisticMCC()  {return relativisticMCC;}
};


class ArMCC : public MCCPackage
{
  static Scalar sigmaElastic(Scalar energy);
  static Scalar sigmaExc(Scalar energy);
  static Scalar sigmaIz(Scalar energy);
  static Scalar sigmaCX(Scalar energy);
  static Scalar sigmaIelastic(Scalar energy);
public:
  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);
  virtual void addCrossSections();
};

class NeMCC : public MCCPackage
{
  static Scalar sigmaElastic(Scalar energy);
  static Scalar sigmaExc(Scalar energy);
  static Scalar sigmaIz(Scalar energy);
public:
  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);
  virtual void addCrossSections(); 
};

class XeMCC : public MCCPackage
{
  static Scalar sigmaElastic(Scalar energy);
  static Scalar sigmaExc(Scalar energy);
  static Scalar sigmaIz(Scalar energy);
public:
  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);

  virtual void addCrossSections();
};


class HMCC : public MCCPackage
{
	static Scalar sigmaIz(Scalar energy);
public:
  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);
  virtual void addCrossSections();
};

/** Monte Carlo collision (MCC) for Helium -- Jeff Hammel **/
class HeMCC : public MCCPackage
{
  static Scalar sigmaElastic(Scalar energy);
  static Scalar sigmaExc(Scalar energy);
  static Scalar sigmaIz(Scalar energy);
  static Scalar sigmaIelastic(Scalar energy);
  static Scalar sigmaCX(Scalar energy);