move.c

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

/**************************************************************
**  This routine interpolates the field at the grids back to the
**  particles and then determines the new velocities, from the
**  Lorentz equations and conservation of angular momentum
**  dv_r/dt = q/m [ E_r + v_t . Bz ] + v_t^2/r
**  dv_t/dt = - q/m v_r . Bz
****************************************************************/
#include "pdc1.h"
#include <xgrafix.h>

void move(int isp)
{
  register int ii;
  int j, lcount, k;
  float frac, ef_part, vr_d;
  float x, y, cos_alpha, sin_alpha;
  float vz_incr;

  /* nb ef_part is the electric field at the particle, B_norm is
  normalised Bz, qm is q/m (species dependent), me_e is me/e */

  lcount = 0;      /* count number particles lost to walls */
  k = sp_k[isp];
  vz_incr = a_scale[isp]*E_z;

/*--------------- Loop over number of particles -------------------------*/
  for (ii=0; ii<np[isp]; ii++)
  {
/* First determine electric field at particle position */
      j = j_pointer(isp,ii);
/*  Linear weighting */
      frac = (r_grid[j+1] - r[isp][ii])/dr_n[j];
      ef_part = frac*efield[j] + (1.0 - frac)*efield[j+1];

/*   Calculate new velocities and positions - involves changing coordinate frames */
        v_r[isp][ii] += a_scale[isp]*ef_part + wc[isp]*v_theta[isp][ii]*k;
        x = r[isp][ii] + v_r[isp][ii]*k;
        y = v_theta[isp][ii]*k;
        r[isp][ii] = sqrt(x*x + y*y);
        cos_alpha = x/r[isp][ii];
        sin_alpha = y/r[isp][ii];

/*  Change coordinate frames */
        vr_d = v_r[isp][ii];
        v_r[isp][ii] =  vr_d*cos_alpha + v_theta[isp][ii]*sin_alpha;
        v_theta[isp][ii] =  -vr_d*sin_alpha + v_theta[isp][ii]*cos_alpha;

/* Determine v_z from E_z */
        v_z[isp][ii] += vz_incr;

/*  If particle has hit wall store index */
        if ((r[isp][ii]-r0)*(r1-r[isp][ii]) <= 0)
        {
          if (r[isp][ii] == 0)
          {
	    #ifdef MPI_1D
              if(my_rank == ROOT) {
	        printf("cos_alpha = %e, sin_alpha = %e \n",cos_alpha, sin_alpha);
                printf("v_r = %e, v_theta = %e \n",v_r[isp][ii], v_theta[isp][ii]);
	      }
	    #else
	      printf("cos_alpha = %e, sin_alpha = %e \n",cos_alpha, sin_alpha);
              printf("v_r = %e, v_theta = %e \n",v_r[isp][ii], v_theta[isp][ii]);
	    #endif
          }
          if (lcount > max_loss)
          {
	   #ifdef MPI_1D
             if(my_rank == ROOT) {
	       printf("cos_alpha = %e, sin_alpha = %e \n",cos_alpha, sin_alpha);
               printf("MOVE -error- too many particles lost \n");
               printf("lcount = %i \n", lcount);
	     }
	   #else
	     printf("cos_alpha = %e, sin_alpha = %e \n",cos_alpha, sin_alpha);
             printf("MOVE -error- too many particles lost \n");
             printf("lcount = %i \n", lcount);
	   #endif
             exit(1);
          }
          lost_index[isp][lcount] = ii;
          lcount ++;
        }

    }   /* end of loop over # particles */

    n_lost[isp] = lcount;
}       /* end of move routine */