mcc.cpp

pic 模拟程序！面向对象

  Vector2 xTemp;
  
  // loop over random selection of particles and see if collision occurs:
  while (nCollisions--){
    index = (int)(n*frand());   //	get extended particle index
    // find it
    for (pgIter.restart(); !pgIter.Done(); pgIter++) {
      if (index < pgIter()->get_n()) break;
      else index -= (int)pgIter()->get_n();
    }
    /** 
     * the initialization of the class member data:
     * ParticleGroup* pg; (declared in mcc.h) follows.
     * Note that this pointer is not initialized in the 
     * constructor of MCCPackage.
     */
    pg = pgIter(); 
    // This is new code to determine if particle lies within the region
    //   of the grid that contains neutral gas. (Bruhwiler 09/26/2000)
    //   If the particle is outside the gas region, then we simply
    //   "continue" back to the top of the while construct.
    xTemp = pg->get_x(index);
    if ( xTemp.e1() < p1Grid.e1() || xTemp.e1() > p2Grid.e1() ||
         xTemp.e2() < p1Grid.e2() || xTemp.e2() > p2Grid.e2() ) continue;
    /**
     * indices for the cell on the grid, the particle is on:
     */
    int iNGD = static_cast<int>(xTemp.e1());
    int jNGD = static_cast<int>(xTemp.e2());
    
    // check that the NeutralGasDensity in the cell is 
    // greater than zero so ionization events can occur, if not
    // return at this point.
    try{
      if ( ptrNGD->getNGD(iNGD, jNGD) > 0.0 ) {
        // do nothing and continue the rest of the loop's body
      } else { 
        // start over at the beginning of the loop by attempting 
        // a collision with another particle
        continue; 
      }
    }
    catch(Oops& oops2){
      oops2.prepend("MCCPackage::collide: Error: \n"); //MCC::collide:
      throw oops2;
    }

    // Particle lies within gas region, so proceed with previous logic
    // mindgame: start with gasVelocity stuff
    Vector3 gasVelocity;
    Vector3 *uPtr = pg->get_u();
    // play with the velocity, not gamma*velocity
    u = uPtr[index]/pg->gamma(index);
    if (collisionModel == ion) {
       // velocity should be relative velocity to the neutral gas.
	gasVelocity = gasDistribution->get_V();
	u -= gasVelocity;
    }

    energy = pg->energy2_eV(u, index);
    v = u.magnitude() + 1E-30;
    
    // mindgame: make it simpler and self-consistent.
#if 0
    // This is part of the original, nonrelativistic algorithm.
    // For this algorithm, in the "relativistic" case, the gamma
    //  factor is extracted from the normalized momentum, leaving
    //  v to be the velocity in m/s.
    BOOL relativistic = FALSE;
    if ( relativisticMCC == 0 ) {
      if (v > 1.e8) {
        relativistic = TRUE;
        u /= pg->gamma(index);
        v = u.magnitude() + 1E-30;
      }
      energy = pg->energy_eV(index); // in eV
    }
    // This is for the new, truly relativistic algorithm.
    // This logic is a mess.  It was done this way to
    //   introduce new capabilities with minimal impact
    //   on existing code.  At some point, we need to
    //   clean this up.
    else {
      Scalar gamma = sqrt(1.+v*v*iSPEED_OF_LIGHT_SQ);
      energy = ELECTRON_MASS_EV*(gamma-1.);
      
      //      cout << "^^^^^^^^" << endl;
      //      cout << " Initially, v=gamma*beta*c = " << v      << endl;
      //      cout << " gamma  is calculated to be: " << gamma  << endl;
      //      cout << " energy is calculated to be: " << energy << endl;
      
      v /= gamma;     // v must be the velocity in m/s (no gamma)
      
      //      cout << " Finally, v=beta*c is found: " << v      << endl;
      //      cout << endl << endl;
    }
#endif
    
    Scalar	random = frand()/v;
    Scalar	sumSigma = 0.;
    if ((collisionModel == electron)||(collisionModel == test)) {
      random *= maxSigmaVe;

      // loop over all possible collisions with this gas
      for (int i=0; i<neCX; i++) {
	
	/**
	 * instead of the old line of code:
	 * sumSigma += ecxFactor*eCX[i]->sigma(energy);
	 * which worked for the uniform gas density we use 
	 * the scaled neutral gas density: dad
	 */
        try{
	        sumSigma += ecxFactor*eCX[i]->sigma(energy)*(ptrNGD->getNGD(iNGD,jNGD))/maxNGD;
        }
        catch(Oops& oops2){
          oops2.prepend("MCCPackage::collide: Error: \n");//OK
          throw oops2;
        }
	///  sumSigma += ecxFactor*eCX[i]->sigma(energy);
      
	if (random <= sumSigma) {
          dynamic(*eCX[i]);
	  break; // leave for loop if collision done
	}
      }
    }
    else if (collisionModel == ion) {
      random *= maxSigmaVi;
      // loop over all possible collisions with this gas
      for (int i=0; i<niCX; i++) {
        try{
	        sumSigma += icxFactor*iCX[i]->sigma(energy)*(ptrNGD->getNGD(iNGD,jNGD))/maxNGD;
        }
        catch(Oops& oops2){
          oops2.prepend("MCCPackage::collide: Error: \n");//OK
          throw oops2;
        }
        ///sumSigma += icxFactor*iCX[i]->sigma(energy);	
	      if (random <= sumSigma) {
	        dynamic(*iCX[i]);
	        break; // leave for loop if collision done
	      }
      }
    }
    
#if 0
    // This isn't needed any more
    if (relativisticMCC == 0) {
      if (relativistic) u *= pg->gamma(u);
    }
#endif
    if (collisionModel == ion) {
      u += gasVelocity;
    }
    u *= pg->gamma2(u);
    // mindgame: save gamma*velocity to the array.
    uPtr[index] = u;
    
  } // end of "while(nCollisions--)" construct
  
  return;
}

// compute null (max) cross section for e-neutral collisions
Scalar MCCPackage::sumSigmaVe()
{
  maxSigmaVe = 0;
  for (Scalar energy=0; energy<MAX_ENERGY; energy+=0.5) 
    {
      Scalar temp = 0.;
      for (int i=0; i<neCX; i++) 
	{ 
	  temp += eCX[i]->sigma(energy);
	}
      maxSigmaVe = MAX(maxSigmaVe, sqrt(energy)*temp);
    }
  maxSigmaVe *= ecxFactor/sqrt(eEnergyScale);
  return maxSigmaVe;
}

// compute null (max) cross section for ion-neutral collisions
Scalar MCCPackage::sumSigmaVi()
{
  maxSigmaVi = 0;
  for (Scalar energy=0; energy<MAX_ENERGY; energy+=0.5) {
	 Scalar temp = 0;
	 for (int i=0; i<niCX; i++) temp += iCX[i]->sigma(energy);
	 // fprintf(stderr, "\nenergy: %g, sigmav=%g, max=%g", 
	 //	 energy, sqrt(energy)*temp, maxSigmaVi);
	 maxSigmaVi = MAX(maxSigmaVi, sqrt(energy)*temp);
  }
  maxSigmaVi *= icxFactor/sqrt(iEnergyScale);
  return maxSigmaVi;
}

void MCCPackage::dynamic(CrossSection& cx)
throw(Oops){
  switch(cx.get_type()) 
	 {
	 case CrossSection::ELASTIC:
		elastic(cx);
		break;
	 case CrossSection::EXCITATION:
		excitation(cx);
		break;
	 case CrossSection::IONIZATION:
     try{
		   ionization(cx);
     }
     catch(Oops& oops3){
       oops3.prepend("MCCPackage::dynamic: Error: \n");
       throw oops3;
     }
		 break;
	 case CrossSection::ION_ELASTIC:
		ionElastic(cx);
		break;
	 case CrossSection::CHARGE_EXCHANGE:
		chargeExchange(cx);
		break;
	 default:
		printf("unrecognized collision type!\n");
	 }
}

void MCCPackage::elastic(CrossSection& cx)
{
  	// This is the original, nonrelativistic algorithm
  if (relativisticMCC == 0) {
    u /= v; // normalize
    newVelocity(energy, v, u, TRUE); // scatter electron
	}

  	// This is the new relativistic algorithm
  else {

		// Calculate the scattering angles.
    Scalar thetaS, phiS;
    elasticScatteringAngles(energy, cx, thetaS, phiS);

		// Now, calculate the new normalized momentum vector
    //   for the scattered electron.
    Vector3 uScatter;
    uScatter = newMomentum(u, energy, thetaS, phiS);
    u = uScatter;
		//    v = u.magnitude();
	}
}

void MCCPackage::excitation(CrossSection& cx)
{
  u /= v;   // is this in the right place? Should be after energy loss?
  energy -= cx.get_threshold();   // subtract inelastic energy loss
  v = sqrt(energy/eEnergyScale);  // recompute magnitude of velocity
  newVelocity(energy, v, u, TRUE);
}

void MCCPackage::ionization(CrossSection& cx)
  throw(Oops){
  	// This is the original, nonrelativistic algorithm
  if (relativisticMCC == 0) {
    u /= v;
    energy -= cx.get_threshold(); // subtract ionization energy
    Scalar rEnergy = 10*tan(frand()*atan(energy/20));
    energy -= rEnergy;			//	partition remaining energy
				
    //	scatter created electron:
    v = sqrt(rEnergy/eEnergyScale);
    Vector3 uNew = u;
    newVelocity(rEnergy, v, uNew, FALSE);
    
    // holds position vector of particle that undergos collision
    Vector2 xTemp;
    xTemp = pg->get_x(index);
    
    /**
     * indices for the cell on the grid, the particle is on:
     */
    int iNGD = static_cast<int>(xTemp.e1());
    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");//done
      throw oops2;
    }
    // Calculate number of physical ions created
    Scalar numIons = pg->get_np2c(index);
    // Check that numNeutral > numIons 
    if ( numNeutralAtoms > numIons ) { 
      // create a computational particle with electron species with
      // number of physical electrons equal to the number of physical
      // ions: numIons
      pList->add(new Particle(pg->get_x(index), uNew*pg->gamma(uNew), eSpecies, 
                              numIons, pg->isVariableWeight() ) );
      
      // continue with creation of the ions 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++) {
	uNew = gasDistribution->get_V();
	pList->add(new Particle(pg->get_x(index), uNew, iSpecies,
				numIons/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 first the electrons:
      pList->add(new Particle(pg->get_x(index), uNew*pg->gamma(uNew), eSpecies, 
                              numNeutralAtoms, true ) );
      
      ptrNGD->set(iNGD, jNGD, 0.0 );
      
      // create n=ionzFraction ions from isotropic gas distribution
      for (int i=0; i<ionzFraction; i++) {
        uNew = gasDistribution->get_V();
        pList->add(new Particle(pg->get_x(index), uNew, iSpecies,
                                numNeutralAtoms/ionzFraction, true));
      }
    }

    // scatter incident electron
    v = sqrt(energy/eEnergyScale);
    newVelocity(energy, v, u, FALSE);
  }
  
  // This is the new relativistic algorithm
  
  /// All new ionization/gas code implemented above is repeated below.
  /// Eventually, the nonrelativistic code above will go away.
  /// --- Let's put the relevant code into helper functions. ---

  else {

		// First, calculate the energy of the ejected electron.
    Scalar eEjected = ejectedEnergy(cx, energy);
    Scalar newEnergy = energy - cx.get_threshold() - eEjected;
		//    cout << "*******************************" << endl;
		//    cout << " impact     energy = " << energy             << endl;
		//    cout << " ejected    energy = " << eEjected           << endl;
		//    cout << " ionization energy = " << cx.get_threshold() << endl;
		//    cout << " remaining  energy = " << newEnergy          << endl;
		//    cout << " energy difference = " << newEnergy-energy   << endl;
		//    cout << endl;

		// Next, calculate the scattering angles of both the
		//   primary electron and the ejected electron (these
		//   angles are wrt initial direction of primary).
    Scalar thetaP, phiP, thetaS, phiS;
    primarySecondaryAngles(cx, energy, eEjected,
                           thetaP, phiP, thetaS, phiS);

    // Now, calculate the new normalized momentum vector
    // for the ejected secondary electron, and add it
    // to the appropriate particle list when the 
    // electron species particle is created.
    Vector3 uEjected = newMomentum(u, eEjected, thetaS, phiS);

    // create n=ionzFraction ions from isotropic gas distribution

    // repeat the NeutralGasDensity reset and creat the ions
    Vector3 uIon;
    // holds position vector of particle that undergos collision
    Vector2 xTemp;
    xTemp = pg->get_x(index);
    
    /**
     * indices for the cell on the grid, the particle is on:
     */
    int iNGD = static_cast<int>(xTemp.e1());