f_thermal.c

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

/************************************************************
** These routines are used for calculating a thermal (maxwellian) 
** velocity distribution for initialising the particle velocities
** and also a thermal flux from the walls, used in wall for injection,
** reflection and secondary emission of particles.
**
*************************************************************/
#include "pdc1.h"

/**************************************************************/
/* Calculates a maxwellian distribution */
float maxwellian()
{
    static int iset=0;
    static float gset;
    float fac,r2,v1,v2;
  
    if(iset==0)
    {
      do
      {
	 v1=2.0*frand()-1.0;
	 v2=2.0*frand()-1.0;
	 r2 = v1*v1+v2*v2;
      } while (r2>=1.0);
      fac=sqrt(-log(r2)/r2);
      gset=v1*fac;
      iset = 1;
      return v2*fac;
    }
    else
    {
       iset=0;
       return gset;
    }
}     /* end of maxwellian */

/*******************************************************************
**  Calculates a maxwellian distribution for initialising particle 
**  particle velocities at startup
*********************************************************************/

float distribution(int k, int isp)
{
    static int init=1;
    static int cutoff=1;
    static int init_sp[NSMAX][2];
    static int nvel[NSMAX][2];
    static float *v[NSMAX][2];
    float erfvl, diff_erf, vlower, vupper, flower, fupper, f, fmid;
    float dv, dvi;
    float nvtc=3.3;    /* number of thermal velocity for upper cutoff */
    int i, n, ksign;
    float vnew;


 /* Set tag to initialise distribution for each new species / direction of motion */

    if(init)
    {
      for(i=0; i<NSMAX; i++)
      for(ksign=0; ksign<2; ksign++) init_sp[i][ksign] = 1;
      init = 0;
    }

/*----- For a cold beam return the beam velocity -----*/
          if (vt[isp][k]==0) return vnew=v0[isp][k];

/*----- For a warm distribution :
             Initialise the velocity distribution -----*/

    if (init_sp[isp][k]) 
    {

/*  If there is a cut-off velocity calculate modified distribution */
       if (vc[isp][k]!=0)  
       {
         erfvl = erf(vc[isp][k]);
	      diff_erf = erf(nvtc)-erfvl;
         nvel[isp][k] = 1.0/(1-diff_erf/(1-erfvl));
	      v[isp][k] = (float *) malloc(nvel[isp][k]*sizeof(float));
	      v[isp][k][0] = vc[isp][k];
	      v[isp][k][nvel[isp][k]-1] = nvtc;
	      dvi = (nvtc-vc[isp][k])/((float)nvel[isp][k]);
	      dv = dvi;
	      fmid = 0;
	      vlower = vc[isp][k];
	      f = 0;
	      flower = erf(vlower)-f*(diff_erf)-erfvl;
	      vupper = vlower + dv;
	      for (n=1; n<(nvel[isp][k]-1); n++)
	      {
	         f = (float)n/(float)(nvel[isp][k]-1);
	         flower = erf(vlower) - f*diff_erf - erfvl;
	         fupper = erf(vupper) - f*diff_erf - erfvl;
	         while ((flower*fupper)>0.0)
	         {
	          dv*=2;
		       vupper = vlower + dv;
		       fupper = erf(vupper) - f*diff_erf - erfvl;
	         }
            while (dv>1e-8)
	         {
	          dv /=2;
		       fmid = erf(vlower+dv) - f*diff_erf - erfvl; 
		       if (fmid<=0) 
		         vlower+=dv;
		       else
		         vupper-=dv;
	         }
	         v[isp][k][n] = v0[isp][k]+vt[isp][k]*vlower;
	      }
       }
       else
           cutoff = 0;
           init_sp[k][isp] = 0;
    }

    if (cutoff)
      vnew = v[isp][k][(int)(nvel[isp][k]*frand())];
    else
      vnew = v0[isp][k]+vt[isp][k]*fabs(maxwellian());

    /* #ifdef MPI_1D
          if (my_rank == ROOT) printf("DISTRIBUTION: vnew = %e \n",vnew);
       #else
          printf("DISTRIBUTION: vnew = %e \n",vnew);
       #endif
    */

    return vnew;                    
}

/**********************************************************************/

float maxwellian_flux(void)
{
  return sqrt(-log(frand()));
}

/********************************************************************
** Determines a maxellian flux f.v Used for determining velocities for
** injection, reflection or secondary emission from the electrodes */

float distribution_flux(int k, int isp)
{
    int ii, jj;
    static int init=1;
    static float vc2[NSMAX][2];
    float v;

    if (init)
    {
       for (ii=0;ii<nsp;ii++) 
         for (jj=0;jj<2;jj++) vc2[ii][jj] = vc[ii][jj]*vc[ii][jj];
       init=0;
    }

    v = v0[isp][k]+vt[isp][k]*sqrt(vc2[isp][k]-log(frand()));
              
    return v;
}
/**************************************************************/