wall.c

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

/**************************************************************
   These routines take particles out which have collided with   
   the walls. Indices of the escaping particles are obtained   
   from the move routine. Deleted particles are replaced by   
   those on the end of the array. The routines also include  
   secondary emission and reflection from the wall, due to  
   particle impact and as well as particle injection (specified
   as input)                                                  
    
    Parameters:                                                   
         n_lost  - # of particles lost to walls 
	 lost_index - index of lost particles                
****************************************************************/
#include "pdc1.h"
#include <xgrafix.h>

void wall(int isp)
{
  static int init_flag=1;
  static float extra[NSMAX][2];

  int ii, k, ind;
  int nnp, j, secountl=0, secountr=0;
  int ebin, abin;
  int ionsp, secsp;
  int nreflux[NSMAX];
  float del_t, energy, angle, add_part;
  float delta_e, delta_a; 
  
  /* INITIALIZE array for computing positions of injected particles */
  if (init_flag)
  {
     /* Set secondary and ion species to array indices */
     if (secondary) secsp = secondary-1;
     ionsp = ionspecies-1;

    /* "enter" is now the no. of particles injected each "dt" */
    for (ii=0; ii<nsp; ii++)
    {
      if (fabs(enter[ii][0]) > 0.) extra[ii][0] = 0.5123123;
      else extra[ii][0] = 0.;
      if (fabs(enter[ii][1]) > 0.) extra[ii][1] = 0.5123123;
      else extra[ii][1] = 0.;
    }
    init_flag = 0;
  }
  
  /* zero variables containing wall current, # secondaries emitted and reflected */
  iwall[isp] = 0.0;
  nreflux[isp]= 0;	    

#ifdef POW
  e_lost[isp] = 0.0;
#endif

/* remove particles impacting electrodes, allow for secondary electron emission, 
    accumulate charge at LH plate for calculating sigma */

  if (n_lost[isp] > 0.0) 
  {
      if (secondary) secountl = secountr = 0;
      nnp = np[isp] -1;     /* Find current no. of particles  */

/* Count down through indices of escaping particles backwards */
      for (ii=n_lost[isp]-1; ii>=0; ii--)
      {
        ind = lost_index[isp][ii];      

/* ---------- Look at particles lost at inner electrode ---------- */
        if (r[isp][ind] < r0) 
        {
             if (ind != nnp) 
             {
/* ------ Copy particle from end of array to escaped particle position ----- */
	       r[isp][ind] = r[isp][nnp];
	       v_r[isp][ind] = v_r[isp][nnp];
	       v_theta[isp][ind] = v_theta[isp][nnp];
	       v_z[isp][ind] = v_z[isp][nnp];
             }
	     nnp--;                      /* decrement # particles by 1 */

	     if(reflux) nreflux[isp]++;  /* Add 1 to no of reflected particles */
             iwall[isp] += intoi[isp]; /* increment wall current */
/*------- Check for secondary emission from electrode ------- */
	     if (secondary)
	     {
	       if (frand() < seec[isp]) /* if probability < random number, make secondary */
	       {
	         secountl++;
                 iwall[secsp] -= intoi[isp];   /* decrement wall charge */
	       }
	     }                 /* end of secondary loop */

	 }                   /* end of LH electrode loop */
	
/* ---------- Look at particles lost at RH electrode ---------- */
	 else if (r[isp][ind] >= r1)
	 {

           if (theRunWithXFlag) /* start diagnostics */
           {
/* Determine energy distribution of escaping particles */
             
 	     energy = v_r[isp][ind]*v_r[isp][ind] + 
	              v_theta[isp][ind]*v_theta[isp][ind] +
		      v_z[isp][ind]*v_z[isp][ind];
	   #ifdef POW
	     e_lost[isp] += Escale[isp]*energy; 
           #endif
		      
 	     energy = (energy - emin[isp])/de[isp];
	     ebin = (int)energy;
	     if (ebin < nbin[isp]-1 && energy >=0)
	     {
	       delta_e = energy - ebin;
	       fe[isp][ebin]  += 1 - delta_e;
	       fe[isp][ebin+1]+= delta_e;
	     }
	    
/* Determine angular distribution of escaping particles */
	     angle = atan(sqrt(v_theta[isp][ind]*v_theta[isp][ind]
			   +v_z[isp][ind]*v_z[isp][ind])/v_r[isp][ind])/dtheta[isp];
	     abin = angle;
	     if (abin < nbin[isp]-1)
	     {
	       delta_a = angle - abin;
	       ftheta[isp][abin]  += 1.0 - delta_a;
	       ftheta[isp][abin+1]+= delta_a;
	     }

           } /* end diagnostics */
	#ifdef POW
	   else {
	     e_lost[isp] += Escale[isp]*(v_r[isp][ind]*v_r[isp][ind] + 
	                                 v_theta[isp][ind]*v_theta[isp][ind] + 
		                         v_z[isp][ind]*v_z[isp][ind]);
           }
	#endif

           if (ind != nnp) 
           {
/* ------ Copy particle from end of array to escaped particle position ----- */
	     r[isp][ind] = r[isp][nnp];
	     v_r[isp][ind] = v_r[isp][nnp];
	     v_theta[isp][ind] = v_theta[isp][nnp];
	     v_z[isp][ind] = v_z[isp][nnp];
           }
	   nnp--;  /* decrement # particles */

/*------- Check for secondary emission from electrode ------- */
	   if (secondary)
	     if (frand() < seec[isp])	secountr++;
	 }
      }     /* end of loop over n_lost particles  */
      
      np[isp] = nnp +1;
    }	/* end of n_lost > 0 loop */
  
  /* ------- INJECT new particles at walls, one species at a time ----- */
  
    if (inject[isp] || reflux) /* if injection or reflection have occured */
    {
      for (k=0; k<2; k++)          /* k=0 LH electrode, k=1 RH electrode */
      {
        add_part = enter[isp][k]+nreflux[isp];
        if (add_part) 
        {
	  extra[isp][k] += add_part;   /* add injected & reflected particles */
	  while (extra[isp][k]>= 1.)
	  {
	    del_t= 1.0/add_part; /* dt for adding particles */
	    ii = np[isp];

	    np[isp]++;
	    if (ii > maxnp[isp]) 
	    {
	    #ifdef MPI_1D
	      if(my_rank == ROOT)
	        printf("WALL- inject- too many particles, species %d",isp);
	    #else
	      printf("WALL- inject- too many particles, species %d",isp);
	    #endif
	      exit(1);
	    }
	  
    /* Choose velocities for new particles  */
	    v_r[isp][ii] = distribution_flux(k,isp);
	    v_theta[isp][ii] = v0_th[isp] + vt_th[isp]*maxwellian(); /* thermal theta velocity */
	    v_z[isp][ii] = v0_z[isp] + vt_z[isp]*maxwellian();       /* thermal z velocity */
	  
    /* Adjust Vr for effect of electric field */
	    v_r[isp][ii] += qm[isp]*me_e*(del_t-0.5)*efield[k*nc];
	    r[isp][ii] = r0 +k*nc +del_t*v_r[isp][ii];

            if (!k) iwall[isp] -= intoi[isp];
	    extra[isp][k] -= 1;
	  }	/* end of while loop */
	}	/* end of add_part loop */
      }	/* end of loop over boundaries (k=0,1) */
    }	/* end of particle injection loop */
  
  /*------- If secondary electrons are turned on, create them here -------- */
  if(secountl+secountr)
  {
    if(np[secsp] > maxnp[secsp])
    {
    #ifdef MPI_1D
      if(my_rank == ROOT)
        printf("WALL(Secondaries): too many particles, MUST EXIT!");
    #else
      printf("WALL(Secondaries): too many particles, MUST EXIT!");
    #endif
      
      exit(1);
    }

/* Add secondaries at LH electrode */
    ii = np[secsp];
    for(j=ii; j< ii+secountl; j++)
    {
      r[secsp][j]= r0 +frand();
      v_r[secsp][j]= distribution_flux(0,secsp);
      v_theta[secsp][j] = v0_th[secsp] + vt_th[secsp]*maxwellian();
      v_z[secsp][j] = v0_z[secsp] + vt_z[secsp]*maxwellian();
    }
/* Add secondaries at RH electrode */
    ii += secountl;
    for(j=ii; j< ii+secountr; j++)
    {
      r[secsp][j]= r1 -frand();
      v_r[secsp][j]= distribution_flux(1,secsp);
      v_theta[secsp][j] = v0_th[secsp] + vt_th[secsp]*maxwellian();
      v_z[secsp][j] = v0_z[secsp] + vt_z[secsp]*maxwellian();
    }
    np[secsp] += secountl + secountr;    /* Add secondaries to particle number */
  }/* end of if(secondary) */
  
}/* end WALL */