argonmcc_hbs.c

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

/*------------------------------------------------------------------------*/
float asigma1a(float energy)
{
  float temp, convert, q_mtm;
  temp = log(1 + energy);
  convert = 2.0*(energy - temp)/(energy*temp);
  q_mtm = convert*asigma1(energy);
  return(q_mtm);
}
/*---------------------------------------------------------------------------*/
/*  e + Ar -> e + Ar  Momentum Transfer elastic scattering */
/* Data taken from 
       M. Suzuki, T.Taniguchi, H.Tagashira, J.Phys.D.: Appl. Phys. 23, 842 (1990)
*/
/*---------------------------------------------------------------------------*/
float asigma1b(float energy)
{
  int i;
  static int init_flag=1;
  float mtm_sigma, alpha;
  static float *mtm_1, *mtm_2;

  /******  Initialization   *********/
  if(init_flag)
  {
    int N=22, i, j=1;
    float engy;
    float energy[23] = { 0.0, 0.05, 0.11, 0.14, 0.2, 0.25, 0.32, 0.4,
         0.5, 0.65, 0.8, 1.0, 1.5, 2.0, 3.0, 6.0, 8.0, 10.0, 12.0,
         15.0, 20.0, 26.0, 30.0};

    float sigma[23] = { 8.0, 1.84, 0.342, 0.235, 0.177, 0.156, 0.151,
         0.182, 0.283, 0.47, 0.68, 1.05, 1.74, 2.48, 4.07, 8.7, 11.7, 13.8,
         14.5, 13.2, 10.4, 8.3, 7.2};

    for (i=0; i<N+1; i++) sigma[i] *= 1.0e-20; /* convert x-sect from angstroms to m^2 */

    mtm_1 = (float *) malloc(101*sizeof(float));
    mtm_2 = (float *) malloc(NEMAX*sizeof(float));

    for(i=0; i<101; i++)
    {
        engy = 0.01*i;
        while (engy > energy[j] && j<N) j++;
        while (engy < energy[j-1]) j--;
        mtm_1[i]=sigma[j]+(sigma[j]-sigma[j-1])*(engy-energy[j])/(energy[j]-energy[j-1]);
    }
    for(i=1; i<NEMAX; i++)
    {
       engy = i;
       while (engy > energy[j] && j<N) j++;
       while (engy < energy[j-1]) j--;
       mtm_2[i]=sigma[j]+(sigma[j]-sigma[j-1])*(engy-energy[j])/(energy[j]-energy[j-1]);
       if (mtm_2[i] < 0) mtm_2[i] = 0;
    }
    init_flag =0;
  }
  /*--------------------------------------------------------------------*/
  if(energy < 1.0)
  { 
    i= 100*energy;
    alpha= 100*energy -i;
    mtm_sigma = mtm_1[i] + alpha*(mtm_1[i+1] - mtm_1[i]);
  }
  else
  { 
    i= energy;
    if(i < NEMAX-1)
    {
      alpha= energy -i;
      mtm_sigma = mtm_2[i] + alpha*(mtm_2[i+1] - mtm_2[i]);
    }
    else
      mtm_sigma = mtm_2[NEMAX-1];
  }
  
  return(mtm_sigma);
}
/*---------------------------------------------------------------------------*/
/*  e + Ar -> e + Ar  Excitation    */
/*---------------------------------------------------------------------------*/
float asigma2(float energy)
{
  int i;
  float exc_sigma, alpha;
  static int init_flag=1;
  static float *excit;

  /*-----  Initialization   --------*/
  if(init_flag)
  {
      float engy, x;
      float a = 0.18, b = 0.85, sig_exc0=1.1e-17;

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

      for(i=0; i<NEMAX; i++)
      {
          engy = i;
          x = engy/eexc;
          if(x < 1.0)
              excit[i] = 0;
          else
             excit[i] = sig_exc0/(pow(eexc,2.0)*pow(x,b))*(a*x-a)/(a*x+1);
      }
      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 asigma3(float energy)
{
  int i;
  float ion_sigma, alpha;
  static int init_flag=1;
  static float *ioniz;

  /******  Initialization   *********/
  if(init_flag)
  {
    float engy, x;
    float a = 0.5, b = 0.87, c=1.4, sig_ion0=6.5e-17;
    ioniz = (float *) malloc(NEMAX*sizeof(float));

    for(i=0; i<NEMAX; i++)
    {
      engy = i;
      x = engy/eion;
      if(x < 1.0) ioniz[i] = 0;
      else
        ioniz[i] = sig_ion0/(eion*eion)/pow(x,b)*pow((a*x - a)/(a*x + 1),c);
    }
    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
*/
/*---------------------------------------------------------------------------*/

float asigma4(float energy)
{
  int i;
  float cx_sigma, alpha;
  static int init_flag=1;
  static float *cx1, *cx2;

  /******  Initialization   *********/
  if(init_flag)
  {
    float engy;

    cx1 = (float *) malloc(1001*sizeof(float));
    cx2 = (float *) malloc(NIMAX*sizeof(float));
    cx1[0] = 11.02e-19; /* cx[en = 1e-4] */
    cx2[0] = 11.02e-19;

    for(i=1; i<1001; i++) {
      engy = 0.001*i;
      cx1[i] = 1.0e-20*(7.0 - 0.38*log(engy))*(7.0 - 0.38*log(engy));
    }

    for(i=1; i<NIMAX; i++)
    {
      engy = i;
      cx2[i] = 1.0e-20*(7.0 - 0.38*log(engy))*(7.0 - 0.38*log(engy));
    }
    init_flag =0;
  }
  /*--------------------------------------------------------------------*/

  if(energy < 1.0) {
    i= 1000*energy;
    alpha= 1000*energy -i;
    cx_sigma = cx1[i] + alpha*(cx1[i+1] - cx1[i]);
  }
  else {
    i= energy;
    if(i < NIMAX-1)
    {
      alpha= energy -i;
      cx_sigma = cx2[i] + alpha*(cx2[i+1] - cx2[i]);
    }
    else
      cx_sigma = cx2[NIMAX-1];
  }

  return(cx_sigma);
}

/*---------------------------------------------------------------------------*/
/*  Ar + Ar+ -> Ar + Ar+   Elastic  Scattering
*/
/*---------------------------------------------------------------------------*/

float asigma5(float energy)
{
  int i;
  float sc_sigma, alpha;
  static int init_flag=1;
  static float *sc1, *sc2;

  /******  Initialization   *********/
  if(init_flag)
  {
    float engy;

    sc1 = (float *) malloc(1001*sizeof(float));
    sc2 = (float *) malloc(NIMAX*sizeof(float));
    sc1[0] = 9.63e-19; /* sc[en = 1e-4] */
    sc2[0] = 9.63e-19;

    for(i=1; i<1001; i++)
    {
      engy = 0.001*i;
      sc1[i] = 1.0e-20*(6.45 - 0.365*log(engy))*(6.45 - 0.365*log(engy));
    }
    
    for(i=1; i<NIMAX; i++)
    {
      engy = i;
      sc2[i] = 1.0e-20*(6.45 - 0.365*log(engy))*(6.45 - 0.365*log(engy));
    }

    init_flag =0;
  }
  /*--------------------------------------------------------------------*/
   
  if(energy < 1.0) {
    i= 1000*energy;
    alpha= 1000*energy -i;
    sc_sigma = sc1[i] + alpha*(sc1[i+1] - sc1[i]);
  }
  else {
    i= energy;
    if(i < NIMAX-1)
    {
      alpha= energy -i;
      sc_sigma = sc2[i] + alpha*(sc2[i+1] - sc2[i]);
    }
    else
    sc_sigma = sc2[NIMAX-1];
  }

  return(sc_sigma);
}

/*---------------------------------------------------------------------------*/
/*  Ar* + e -> Ar+ + 2e  Ionization of excited states by electron impact
*/
/*---------------------------------------------------------------------------*/

float asigma6(float energy)
{
    int i, j, N;
    float engy, alpha, q;
    static int init_flag=1;
    static float *ioniz2;
/*
    static float eng[13] = { 4.16, 5.0, 10., 15., 20., 30., 40., 50., 60.,
                             70., 80., 100, 200};
    static float sig[13] = { 0.0, 4.6, 8.1, 8.4, 7.8, 7.0, 5.5, 5.1, 4.5,
                             4.2, 3.8, 3.2, 2.0};
*/
    static float eng[13] = { 0.0, 3.4, 4.21, 4.76, 5.44, 6.8, 10.2, 17.0, 27.2,
                             40.8, 54.4, 81.6, 108.8 };
    static float sig[13] = { 0.0, 0.0, 7.18, 8.05, 8.37, 8.33, 7.17, 5.07,
                             3.38, 2.40, 1.84, 1.25, 0.951};
    if (init_flag) {
      ioniz2 = (float *) malloc(NEMAX*sizeof(float));
      N=12;
      for(i=0; i<NEMAX; i++) {
        engy = i;
        if (engy>eng[N]) {
          ioniz2[i]=sig[N]+(sig[N]-sig[N-1])*(engy-eng[N])/(eng[N]-eng[N-1]);
          if (ioniz2[i]<0) ioniz2[i]=0;
        }
        else if (engy<4.16) ioniz2[i]=0;
        else
          for (j=1;j<=N;j++) if (engy>eng[j-1] && engy<=eng[j])
            ioniz2[i]=sig[j]+(sig[j]-sig[j-1])*(engy-eng[j])/(eng[j]-eng[j-1]);
        ioniz2[i]*=1e-20;
      }
      init_flag=0;
    }

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

    return(q);
}

/*---------------------------------------------------------------------------*/
/*  Ar + e -> Ar^m + 2e  Excitation to metastable states by electron impact
*/
/*------------------------------------------------------------------LEEHJ----*/

float asigma7(float energy)
{
    int i, j, N;
    float engy, alpha, q;
    static int init_flag=1;
    static float *meta;
    static float en[18] = { 11.55, 12.0, 13.0,  15.0, 17.0, 18.0, 20.0, 21.0,
                             22.0, 27.0, 35.0, 40.0, 55.0, 60.0, 82.0, 102.0,
                             122.0, 142.0};
    static float si[18] = { 0.0, 0.015, 0.07, 0.113, 0.26, 0.3, 0.345, 0.35,
                             0.355, 0.3, 0.188, 0.139, 0.06, 0.05, 0.039,
                             0.035, 0.031, 0.028};

    if (init_flag) {
      meta = (float *) malloc(NEMAX*sizeof(float));
      N=17;
      for(i=0; i<NEMAX; i++) {
        engy = i;
        if (engy>en[N]) {
          meta[i]=si[N]+(si[N]-si[N-1])*(engy-en[N])/(en[N]-en[N-1]);
          if (meta[i]<0) meta[i]=0;
        }
        else if (engy<11.55) meta[i]=0;
        else
          for (j=1;j<=N;j++)
            if (engy>en[j-1] && engy<=en[j])
              meta[i]=si[j]+(si[j]-si[j-1])*(engy-en[j])/(en[j]-en[j-1]);
        meta[i]*=1e-20;
      }
      init_flag=0;
    }

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

    return(q);
}
/*---------------------------------------------------------------------------*/