argonmcc_hbs.c

一维等离子体静电粒子模拟程序

	if(theRunWithXFlag)
	{
          s = j_pointer(ecolsp,nnp);
	  frac= (r_grid[s+1] - r[ecolsp][nnp])/dr_n[s];
	  temp = frac/vol[s];
	  mccrate[1][s]  += temp;
	  Pie[s]  += eexc*temp;
	  temp = (1.0 - frac)/vol[s+1]; 
	  mccrate[1][s+1]+= temp; 
	  Pie[s+1]+= eexc*temp; 
          nu_ex ++; 
	}
     #ifdef POW
        else nu_ex++;
     #endif 
      }  

      /*------------------------------------------------------------*/
      /* if the collision is ionization, add the released electron and ion first */
      
      else if(engy>eion && rand_num <= (sum_prob += asigma3(engy)/sigma_total))
      {
	/* first normalize the velocity components by the speed, vel */
	v_r[ecolsp][nnp] /= vel;
	v_theta[ecolsp][nnp] /= vel;
	v_z[ecolsp][nnp] /= vel;
	
        /*------------------------------------------------------------*/
	/* subtract the ionisn energy and partition the remaining energy */
	
	engy -= eion;
	rengy = 10.0*tan(frand()*atan(engy/20.0));
	engy -= rengy;
	
        /*------------------------------------------------------------*/
	/* scatter the created electron */
	
	vel = sqrt(rengy/Escale[ecolsp]);
	k = np[ecolsp];
	v_r[ecolsp][k] = v_r[ecolsp][nnp];
	v_theta[ecolsp][k] = v_theta[ecolsp][nnp];
	v_z[ecolsp][k] = v_z[ecolsp][nnp];
	newvel(rengy, vel, &v_r[ecolsp][k], &v_theta[ecolsp][k], &v_z[ecolsp][k], 0);
	r[ecolsp][k] = r[ecolsp][nnp];

        /* determine the charge density due to new electron */
        s = j_pointer(ecolsp,nnp);
	frac= (r_grid[s+1] - r[ecolsp][nnp])/dr_n[s];
        sp_n_mcc[ecolsp][s]  += frac;
        sp_n_mcc[ecolsp][s+1]+= (1.0 -  frac); 
	
        /*------------------------------------------------------------*/
	/* create ion and assign velocities if ion species exist */
	
        if (ionspecies)
        {
	  k = np[ionsp];
	  maxwellv(&v_r[ionsp][k], &v_theta[ionsp][k], &v_z[ionsp][k], vg_therm);
	  r[ionsp][k] = r[ecolsp][nnp];

          /* determine the charge density due to newly created ion */
          sp_n_mcc[ionsp][s]  += frac;
          sp_n_mcc[ionsp][s+1]+= (1.0 -  frac); 

        }

        /*------------------------------------------------------------*/
	/* finally scatter the incident electron */
	
	vel = sqrt(engy/Escale[ecolsp]);
	newvel(engy, vel, &v_r[ecolsp][nnp], &v_theta[ecolsp][nnp], &v_z[ecolsp][nnp], 0);
	
	if(++np[ionsp] >maxnp[ionsp] || ++np[ecolsp] >maxnp[ecolsp])
	{  
   #ifdef MPI_1D
      if (my_rank == ROOT) printf("mcc(Ionization): too many ions created ");
   #else
      printf("mcc(Ionization): too many ions created ");
   #endif

	 exit(1);
	}
	
        /*------------------------------------------------------------*/
	/* Determine the ionisation profile */
	if(theRunWithXFlag)
	{
	  temp = frac/vol[s];
	  mccrate[2][s]  += temp;
	  Pie[s]  += eion*temp;
	  temp = (1.0 - frac)/vol[s+1]; 
	  mccrate[2][s+1]+= temp; 
	  Pie[s+1]+= eion*temp; 
          nu_iz ++;
	 }
      #ifdef POW
         else nu_iz++;
      #endif
      }   /* end of test for type of collision */
    }     /* end of loop over number colliding electrons*/
  }       /* end of "if (ecollisional)" */

  /*------------------------------------------------------------*/

  if(icollisional && isp==icolsp)
  {
    if(theRunWithXFlag)
      for (j=0; j<ng; j++) mccrate[3][j] = mccrate[4][j] = Pcx[j] = 0.0;

 #ifdef POW
    cx_loss = 0.0;
 #endif
      
    icol_extra += np[icolsp]*col_prob_i;
    N = icol_extra;
    icol_extra -= N;
    nnp = np[icolsp];

   /*------------------------------------------------------------*/
    for(j=0; j< N; j++)
    {
      nnp--;
      index= nnp*frand();     /* Pick a random ion to collide */
      
/* Move colliding particle to end of array, so it isn't chosen twice */
      swap(r[icolsp],index,nnp);
      swap(v_r[icolsp],index,nnp);
      swap(v_theta[icolsp],index,nnp);
      swap(v_z[icolsp],index,nnp); 

   /*------------------------------------------------------------*/
   /* Check for type of collision */

      maxwellv(&vneutr, &vneutth, &vneutz, vg_therm);
      /* Determine relative velocity btwn neutral and ion */
      v_r[icolsp][nnp] -= vneutr;
      v_theta[icolsp][nnp] -= vneutth;
      v_z[icolsp][nnp] -= vneutz;
      v_sqrd = (v_r[icolsp][nnp]*v_r[icolsp][nnp]
          +v_theta[icolsp][nnp]*v_theta[icolsp][nnp] +v_z[icolsp][nnp]*v_z[icolsp][nnp]);
      engy= Escale[icolsp]*v_sqrd;
      vel = sqrt(v_sqrd);
      sigma_total = max_sigmav_i/(vel*vscale);
      rand_num= frand();
      
       /*------------------------------------------------------------*/
      /* if the collision is scattering */
      
      if (rand_num <= (sum_prob = asigma5(engy)/sigma_total))
      {
	float up1, up2, up3, mag;
	float r11, r12, r13, r21, r22, r23, r31, r32, r33;
	
	coschi= sqrt(frand());
	sinchi= sqrt(1. -coschi*coschi);
	
	phi1  = TWOPI*frand();
	cosphi= cos(phi1);
	sinphi= sin(phi1);
	
	r13 = v_r[icolsp][nnp]/vel;
	r23 = v_theta[icolsp][nnp]/vel;
	r33 = v_z[icolsp][nnp]/vel;
	
	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);
	
	r12/= mag;
	r22/= mag;
	r32/= mag;
	
	r11 = r22*r33 -r32*r23;
	r21 = r32*r13 -r12*r33;
	r31 = r12*r23 -r22*r13;
	
	v_r[icolsp][nnp]= vel*coschi*(r11*sinchi*cosphi 
              +r12*sinchi*sinphi +r13*coschi);
	v_theta[icolsp][nnp]= vel*coschi*(r21*sinchi*cosphi 
              +r22*sinchi*sinphi +r23*coschi); 
	v_z[icolsp][nnp]= vel*coschi*(r31*sinchi*cosphi 
              +r32*sinchi*sinphi +r33*coschi);

     /* Count number of elastic collisions */ 
        if(theRunWithXFlag)
	{
          s = j_pointer(icolsp,nnp);
	  frac= (r_grid[s+1] - r[icolsp][nnp])/dr_n[s];
	  mccrate[4][s]  += frac/vol[s]/sp_k[icolsp];
	  mccrate[4][s+1]+= (1.0 - frac)/vol[s+1]/sp_k[icolsp]; 
	}
      }   /* end of "if rand_num < sigma " loop */
      
       /*------------------------------------------------------------*/
      /* if the collision is charge exchange */
      
      else if (rand_num <= (sum_prob += asigma4(engy)/sigma_total))
      {
	v_r[icolsp][nnp] = v_theta[icolsp][nnp] = v_z[icolsp][nnp] = 0.0;
	
     /* Count number of charge exchange collisions */ 
        if(theRunWithXFlag)
	{
          s = j_pointer(icolsp,nnp);
	  frac= (r_grid[s+1] - r[icolsp][nnp])/dr_n[s];
	  temp = frac/vol[s]/sp_k[icolsp];
	  mccrate[3][s]  += temp;
	  Pcx[s] += engy*temp;
	  temp = (1.0 - frac)/vol[s+1]/sp_k[icolsp];
	  mccrate[3][s+1]+= temp; 
	  Pcx[s+1] += engy*temp;
       #ifdef POW
          cx_loss += engy;
       #endif
	}
     #ifdef POW
        else cx_loss += engy;
     #endif
      }  /* end of "else rand_num < sigma + dsig" loop */

    /* Add neutral velocity back to ion velocity */

      v_r[icolsp][nnp] += vneutr;
      v_theta[icolsp][nnp] += vneutth;
      v_z[icolsp][nnp] += vneutz;
    }    /* end of loop of colliding particles */
  }      /* end of if icollisional loop */

}       /* end of mcc subroutine */

/*************************************************************
**  Calculates new velocities for particles which have 
**  just made a collision
**************************************************************/
void newvel(float energy, float vel, float *vr, float *vth, float *vz, int e_flag)
{
  float phi1, cosphi, sinphi, coschi, sinchi, up1, up2, up3;
  float mag, r11, r12, r13, r21, r22, r23, r31, r32, r33;
  float tmp;

  if(energy < 1e-30)  coschi = 1.0;
  else  coschi = (energy +2.0 -2.0*pow(energy +1.0,frand()))/energy;
  sinchi= sqrt(1. - coschi*coschi);
  
  phi1  = TWOPI*frand();
  cosphi= cos(phi1);
  sinphi= sin(phi1);
  
  if (e_flag) {
     tmp = 2.0*mass[ecolsp]*(1-coschi)/mass[ionsp];
     vel *= sqrt(1.0 - tmp);
     en_loss = energy*tmp;
  }
  
  r13 = *vr;
  r23 = *vth;
  r33 = *vz;
  
  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);
  
  r12/= mag;
  r22/= mag;
  r32/= mag;
  
  r11 = r22*r33 -r32*r23;
  r21 = r32*r13 -r12*r33;
  r31 = r12*r23 -r22*r13;
  
  *vr= vel*(r11*sinchi*cosphi +r12*sinchi*sinphi +r13*coschi);
  *vth= vel*(r21*sinchi*cosphi +r22*sinchi*sinphi +r23*coschi); 
  *vz= vel*(r31*sinchi*cosphi +r32*sinchi*sinphi +r33*coschi);
}

/*---------------------------------------------------------------------------*/
/* XSections for Argon-Electron Collisions                    
Curve fittings from Helen Smith's thesis                   
Experimental data used for fitting:

Electron-Argon Neutral 
----------------------
Excitation: 
     deHeer, Jansen and Kaay, J.Phys.B, 12, 979 (1979)
Ionisation: 
     Krishnakumar & Srivastava J.Phys.B. 21, 1055 (1988)
     Vikor et al, FIZIKA 21, 345 (1989)
Elastic 0 - 20 : 
     Ferch, Granitza, Masche & Raith J.Phys.B, 18, 967 (1985)
       20 - 3000 :
     deHeer, Jansen and Kaay, J.Phys.B, 12, 979 (1979)

Ar+ - Ar
----------------------
Elastic Scattering: 
   McDaniel "Collision Phenomena in Ionised Gases" (1964)
Charge Exchange:
   McDaniel "Collision Phenomena in Ionised Gases" (1964)
   Rapp & Francis, J.Chem. Phys 37, 2631 (1962)
   Hegerberg, Stefansson & Elford, J.Phys.B 11, 133 (1978)    */

/*------------------------------------------------------------*/
/*  Total e- Scattering Cross-section                            */
/*-------------------------------------------------------------*/
float asigmaT(float energy)
{
  return(asigma1(energy) + asigma2(energy) + asigma3(energy));
}

/*------------------------------------------------------------*/
/*  e + Ar -> e + Ar  Elastic                                  */
/*-------------------------------------------------------------*/
float asigma1(float energy)
{
  int i;
  float els_sigma, alpha;
  static float *elastic1, *elastic2;
  static int init_flag=1;

  /******  Initialization   *********/
  if(init_flag) {
    float engy;
    float amp_1 = 8.32414e-22, pwr_1 = -1.14818;
    float amp_2 = 1.10478e-20, pwr_2 = 1.25583;
    float amp_3 = 1.17039e-18, pwr_3= -0.669049;

    elastic1 = (float *) malloc(101*sizeof(float));
    elastic2 = (float *) malloc(NEMAX*sizeof(float));

    /*****  With the Ramseur min.  *****/
    elastic1[0] = amp_1*pow(0.005,pwr_1);
    for(i=1; i<101; i++) {
      engy = .01*i;
      if(engy < 0.345) elastic1[i]= amp_1*pow(engy,pwr_1);
      else elastic1[i]= amp_2*pow(engy,pwr_2);
    }

    for(i=1; i<NEMAX; i++) {
      engy = i;
      if (engy < 12.0)
        elastic2[i]= amp_2*pow(engy,pwr_2);
      else elastic2[i] = amp_3*pow(engy,pwr_3);
    }
    init_flag =0;
  }
  /*--------------------------------------------------------------------*/

  if(energy < 1.0) {
    i= 100*energy;
    alpha= 100*energy -i;
    els_sigma = elastic1[i] + alpha*(elastic1[i+1] - elastic1[i]);
  }
  else {
    i= energy;
    if(i < NEMAX-1) {
      alpha= energy -i;
      els_sigma = elastic2[i] + alpha*(elastic2[i+1] - elastic2[i]);
    }
    else
      els_sigma = elastic2[NEMAX-1];
  }

  return(els_sigma);
}
/*------------------------------------------------------------------------*/
/* Momentum transfer cross-section - convert from elastic x-s assuming     
anisotropic distribution used by V.Vahedi et al                           */