argonmcc_lk.c

XPDC1是针对静电模型中的等离子体模拟程序！通过INP文件输入参数有较好的通用性

	/* 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];
	newvelLK(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]);
	newvelLK(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 charge exchange */
      
      if (rand_num <= 4.0e-19/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 newvelLK(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;

  coschi = 1.0 - 2.0*frand();  /* isotropic scattering angle */
  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);
}

/*------------------------------------------------------------*/
/*  Total e- Scattering Cross-section                            */
/*-------------------------------------------------------------*/
float LK_asigmaT(float energy)
{
  return(LK_asigma1(energy) + LK_asigma2(energy) + LK_asigma3(energy));
 }
  

/*---------------------------------------------------------------------------*/
/* XSections for Argon-Electron Collisions                      */
/*------------------------------------------------------------*/
/*  e + Ar -> e + Ar  Elastic                                  */
/*-------------------------------------------------------------*/
float LK_asigma1(float energy)
{
  int i;
  float els_sigma, alpha;
  static int init_flag=1;
  static float *elastic;

  /*-----  Initialization   --------*/
  if(init_flag) {
    float engy;
    float amp_1 = 1.59e-19;

    elastic = (float *) malloc(NEMAX*sizeof(float));
  
    for(i=0; i<NEMAX; i++)
    {
      engy = i;
      if (engy <= eexc) elastic[i]= amp_1*engy/eexc;
      else elastic[i] = amp_1*sqrt(eexc/engy);
    }

    init_flag =0;
  }
/*------------------------------------------------------------------------*/
  i= energy;
  if(i < NEMAX-1)
  {
    alpha= energy -i;
      els_sigma = elastic[i] + alpha*(elastic[i+1] - elastic[i]);
    }
    else
      els_sigma = elastic[NEMAX-1];

  return(els_sigma);
}
/*------------------------------------------------------------------------*/
/* Momentum transfer cross-section - calculate from the elastic x-s       */
/*------------------------------------------------------------------------*/
float LK_asigma1a(float energy)
{
  float q_mtm;  
  q_mtm = LK_asigma1(energy);
  return(q_mtm);
}
/*---------------------------------------------------------------------------*/
/*  e + Ar -> e + Ar  Excitation    */
/*---------------------------------------------------------------------------*/
float LK_asigma2(float energy)
{
  int i;
  float exc_sigma, alpha;
  static int init_flag=1;
  static float *excit;

  /*-----  Initialization   --------*/
  if(init_flag)
  {
       float x;
       float amp = 1.56e-20;

       excit = (float *) malloc(NEMAX*sizeof(float));

       for(i=0; i<NEMAX; i++)
       {
            x = (float)i/eexc;     /* x = energy/eexc */
            if(x < 1.0)
                excit[i] = 0;
            else
                excit[i] = amp*log(x)/x;
       }
       init_flag =0;
  }
  /*--------------------------------------------------------------------*/

  i= energy;
  if(i < NEMAX-1)
  {
    alpha= energy -i;
    exc_sigma = excit[i] + alpha*(excit[i+1] - excit[i]);
  }
  else
    exc_sigma = excit[NEMAX-1];

  return(exc_sigma);

}

/*---------------------------------------------------------------------------*/
/*  e + Ar -> e + e + Ar+  Ion.     */
/*---------------------------------------------------------------------------*/

float LK_asigma3(float energy)
{
  int i;
  float ion_sigma, alpha;
  static int init_flag=1;
  static float *ioniz;

  /*-----  Initialization   --------*/
  if(init_flag)
  {
    float amp=3.18e-20, x;

    ioniz = (float *) malloc(NEMAX*sizeof(float));

    for(i=0; i<NEMAX; i++)
    {
      x = (float)i/eion;    /* x = energy/eion */
      if (x < 1.0)
          ioniz[i] = 0;
      else
          ioniz[i] = amp*log(x)/x;
    }
    init_flag =0;
  }
  /*--------------------------------------------------------------------*/
  i= energy;
  if(i < NEMAX-1)
  {
    alpha= energy -i;
    ion_sigma = ioniz[i] + alpha*(ioniz[i+1] - ioniz[i]);
  }
  else
    ion_sigma = ioniz[NEMAX-1];

  return(ion_sigma);
}
/*---------------------------------------------------------------------------*/
/* Argon Ion-neutral Collision Cross-sections                                */
/*---------------------------------------------------------------------------*/
/*  Ar + Ar+ -> Ar+ + Ar  Charge Exchange       LK_asigma4 = 4.0e-19 */
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
/*  Ar + Ar+ -> Ar + Ar+   No Elastic Scattering in K&L Model               */
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
/*  Ar* + e -> Ar+ + 2e  Ionization of excited states by electron impact     */
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
/*  Ar + e -> Ar^m + 2e  Excitation to metastable states by electron impact */
/* Not included in K&L cross-sections                                       */
/*------------------------------------------------------------------LEEHJ----*/