maxwellv.c

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

#include "pdc1.h"
#include "xgrafix.h"

float F(float);

/*****************************************************************
** Calculates a maxwellian velocity distribution for gas atoms 
** Used in loading new ions created by ionisation or in charge
** exchange collisions
****************************************************************/
void maxwellv(float *vx_local, float *vy_local, float *vz_local, float vth)
{
  static int nvel, init_flag= 1;
  static float *vsave_x, *vsave_y, *vsave_z;
  float nvts = 3.3;
          
  int i, n;
  float vmag, aphi, dv, rr, sintheta, costheta;
	          
  if (init_flag) 
  {
     nvel= 1./(1-F(nvts));
     if(nvel > 100000)
   #ifdef MPI_1D
     if (my_rank == ROOT)
       puts("Warning: Your choice of nvts has made nvel > 1e5"); 
   #else
     puts("Warning: Your choice of nvts has made nvel > 1e5"); 
   #endif
     dv = sqrt(PI)/(4.0*nvel);

     vsave_x= (float *) malloc(nvel*sizeof(float));
     vsave_y= (float *) malloc(nvel*sizeof(float));
     vsave_z= (float *) malloc(nvel*sizeof(float));
     init_flag= 0;
     i=n=0;
     for (n=0; n<nvel; n++) 
     {
        rr=(1.0*n)/nvel;
	while (F(i*dv)< rr) i++;
	vmag=i*dv;
	aphi=2*PI*frand();
	costheta = 1-2*frand();
	sintheta = sqrt(1-costheta*costheta);
	vsave_x[n] = vmag*sintheta*cos(aphi);
	vsave_y[n] = vmag*sintheta*sin(aphi);
	vsave_z[n] = vmag*costheta;
     }
  }
  n = (nvel-1)*frand();
  *vx_local = vth*vsave_x[n];
  *vy_local = vth*vsave_y[n];
  *vz_local = vth*vsave_z[n];
}
/**************************************************************/
float F(float v)
{
  return(-2*v*exp(-v*v)/sqrt(PI) +erf(v));
}
/*---------------------------------------------------------------------------*/
/* This routine swaps 2 array elements */
void swap (float v[], int i, int j)
{
      float temp;
 
      temp=v[i];  
      v[i]=v[j];  
      v[j]=temp;  
}
/*---------------------------------------------------------------------------*/