load.c

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

#include "pdc1.h"
/***********************************************************/
/* This routine loads particles into a new simulation.     */
/* There are 3 density profiles available for loading, the */
/* one to use is specified in the species input:           */
/* profile =  - uniform profile                            */
/* 0 - uniform profile                                     */
/* 1 - cosine profile                                      */ 
/* 2 - J0 profile                                          */
/* The velocity distribution function                      */ 
/* is determined from the input parameters: if temperature */
/* is 0, a cold beam is loaded, otherwise a thermal        */
/* (maxwellian) distribution is calculated and loaded      */
/***********************************************************/
void unif_load(int initnp[NSMAX][2])
{
  int i, ix, k, isp, thermal=1;
  float inv_den, nend;
  
  /************************************************/
  /* Load one species and one direction at a time */
  
  for (isp=0; isp<nsp; isp++) 
  {
    np[isp]= 0;

    for (k=0; k<2; k++) 
    {
      nend = initnp[isp][k];
      if (!nend) continue;
      
/* Check that number of particles is within bounds */
      if (np[isp] + nend > maxnp[isp]) 
      {
       #ifdef MPI_1D
	 if (my_rank == ROOT)
	   printf("LOAD: too many particles, species %d ",isp);
       #else
         printf("LOAD: too many particles, species %d ",isp);
       #endif
         exit(1);
      }

      if (vt[isp][k]==0. ) thermal = 0;   /* determine whether thermal or cold beam */

/*----- Load particle positions  & velocities  -------------- */

/*---- Different loaders for r0>0 ----------- */
      if (r0) 
      {
        inv_den = (r1-r0)*(r1+r0)/nend;
        ix = np[isp];
        r[isp][ix]= sqrt(inv_den + r0*r0);
        if (thermal) 
           v_r[isp][ix] =  distribution(k,isp);
        else 
           v_r[isp][ix] = v0[isp][k];

        for(i=1; i<nend; i++) 
        {
	 ix = i +np[isp];
	 r[isp][ix] = sqrt(inv_den + r[isp][ix-1]*r[isp][ix-1]);
         if (thermal) 
	     v_r[isp][ix] =  distribution(k,isp); 
         else 
	     v_r[isp][ix] = v0[isp][k];
         v_theta[isp][ix] = v0_th[isp] + vt_th[isp]*maxwellian();
         v_z[isp][ix] = v0_z[isp] + vt_z[isp]*maxwellian();
        }
/*------------------ and r0=0 ----------- */
      } else
      {
        for(i=0; i<nend; i++) 
        {
	  ix = i + np[isp];
	  r[isp][ix] = r1*sqrt((i+0.5)/nend);
          if (thermal) 
	     v_r[isp][ix] =  distribution(k,isp); 
          else 
	     v_r[isp][ix] = v0[isp][k];
          v_theta[isp][ix] = v0_th[isp] + vt_th[isp]*maxwellian();
          v_z[isp][ix] = v0_z[isp] + vt_z[isp]*maxwellian();
        }
      }

 /*   Put a sinusoidal perturbation on uniform distribution */
      if(dde) 
      {
      	 for(i=0; i<nend; i++) 
	 {
	   ix = np[isp] +i;
	   r[isp][ix] += dde*sin(TWOPI*(r[isp][ix] -r0)/(r1 -r0));
	 }
      }

      np[isp] +=nend;
    }	/* end  "for( k=0,1...)" loop */
      #ifdef MPI_1D
        if (my_rank == ROOT)
	  printf("proc 0: Species %d: %d particles loaded \n",isp, np[isp]);
      #else
        printf("Species %d: %d particles loaded \n",isp, np[isp]);
      #endif
  }	/* end for(isp=0 thru nsp-1..)" loop */
}	/* end LOAD  */

/***************************************************************/
void cos_load(int initnp[NSMAX][2], float iinitn)
{
  int i, ix, k, isp, thermal=1;
  int nend, np_tot;
  float d0, frac, length;
  
  d0 = nc2p/(PI*height*iinitn*dr*dr);
  length = r1 - r0;

  /************************************************/
  /* Load one species and one direction at a time */
  
  for (isp=0; isp<nsp; isp++) 
  {
    np_tot = initnp[isp][0] + initnp[isp][1];
    np[isp]= 0;

    for (k=0; k<2; k++) 
    {
      nend = initnp[isp][k];
      frac = (float)np_tot/nend;
      if (!nend) continue;
      
/* Check that number of particles is within bounds */
      if (np[isp] + nend > maxnp[isp]) 
      {
       #ifdef MPI_1D
         if (my_rank == ROOT)
	   printf("LOAD: too many particles, species %d ",isp);
       #else
         printf("LOAD: too many particles, species %d ",isp);
       #endif
         exit(1);
      }

      if (vt[isp][k]==0. ) thermal = 0;   /* determine whether thermal or cold beam */

/*----- Load particle positions  & velocities  -------------- */
      ix = np[isp];
      r[isp][ix]= sqrt(r0*r0 + d0*frac);
      if (thermal) 
           v_r[isp][ix] =  distribution(k,isp);
      else 
           v_r[isp][ix] = v0[isp][k];

      for(i=1; i<nend; i++) 
      {
	 ix = i +np[isp];
	 r[isp][ix] = sqrt(r[isp][ix-1]*r[isp][ix-1] + frac*d0/cos(PI*0.5*(r[isp][ix-1]-r0)/length));
         if (r[isp][ix] >= nc) 
         {
           nend=i-1;
           i = nend;
         }
         if (thermal) 
	     v_r[isp][ix] =  distribution(k,isp); 
         else 
	     v_r[isp][ix] = v0[isp][k];
         v_theta[isp][ix] = v0_th[isp] + vt_th[isp]*maxwellian();
         v_z[isp][ix] = v0_z[isp] + vt_z[isp]*maxwellian();
      }

 /*   Put a sinusoidal perturbation on uniform distribution */
      if(dde) 
      {
      	 for(i=0; i<nend; i++) 
	 {
	   ix = np[isp] +i;
	   r[isp][ix] += dde*sin(TWOPI*(r[isp][ix] -r0)/(r1 -r0));
	 }
      }

      np[isp] +=nend;
    }	/* end  "for( k=0,1...)" loop */
      #ifdef MPI_1D
        if (my_rank == ROOT)
	  printf("proc 0: Species %d: %d particles loaded \n",isp, np[isp]);
      #else
        printf("Species %d: %d particles loaded \n",isp, np[isp]);
      #endif
    
  }	/* end for(isp=0 thru nsp-1..)" loop */
}	/* end LOAD  */

/***************************************************************/
void bes_load(int initnp[NSMAX][2], float iinitn)
{
  int i, ix, k, isp, thermal=1;
  int nend, np_tot;
  float d0, frac, length;
  
  d0 = nc2p/(PI*height*iinitn*dr*dr);
  length = r1 - r0;

  /************************************************/
  /* Load one species and one direction at a time */
  
  for (isp=0; isp<nsp; isp++) 
  {
    np_tot = initnp[isp][0] + initnp[isp][1];
    np[isp]= 0;

    for (k=0; k<2; k++) 
    {
      nend = initnp[isp][k];
      frac = (float)np_tot/nend;
      if (!nend) continue;
      
/* Check that number of particles is within bounds */
      if (np[isp] + nend > maxnp[isp]) 
      {
       #ifdef MPI_1D
	 if(my_rank == ROOT)
           printf("LOAD: too many particles, species %d ",isp);
       #else
         printf("LOAD: too many particles, species %d ",isp);
       #endif
         exit(1);
      }

      if (vt[isp][k]==0. ) thermal = 0;   /* determine whether thermal or cold beam */

/*----- Load particle positions  & velocities  -------------- */
      ix = np[isp];
      r[isp][ix]= sqrt(r0*r0 + d0*frac);
      if (thermal) 
           v_r[isp][ix] =  distribution(k,isp);
      else 
           v_r[isp][ix] = v0[isp][k];

      for(i=1; i<nend; i++) 
      {
         ix = i +np[isp];
         r[isp][ix] = sqrt(r[isp][ix-1]*r[isp][ix-1] + frac*d0/J_0((r[isp][ix-1]-r0)/length));
         if (r[isp][ix] >= nc) 
         {
           nend=i-1;
           i = nend;
         }
         if (thermal) 
             v_r[isp][ix] =  distribution(k,isp); 
         else 
             v_r[isp][ix] = v0[isp][k];
         v_theta[isp][ix] = v0_th[isp] + vt_th[isp]*maxwellian();
         v_z[isp][ix] = v0_z[isp] + vt_z[isp]*maxwellian();
      }

 /*   Put a sinusoidal perturbation on uniform distribution */
      if(dde) 
      {
         for(i=0; i<nend; i++) 
         {
           ix = np[isp] +i;
           r[isp][ix] += dde*sin(TWOPI*(r[isp][ix] -r0)/(r1 -r0));
         }
      }

      np[isp] +=nend;
    }   /* end  "for( k=0,1...)" loop */
      #ifdef MPI_1D
        if (my_rank == ROOT)
	  printf("proc 0: Species %d: %d particles loaded \n",isp, np[isp]);
      #else
        printf("Species %d: %d particles loaded \n",isp, np[isp]);
      #endif
  }     /* end for(isp=0 thru nsp-1..)" loop */
}       /* end LOAD  */

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