diagnostics.c

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

     d_skip *= 4;
   #ifdef MPI_1D
     if(my_rank == ROOT)
       printf("dskip = %d \n",d_skip);
   #else
     printf("dskip = %d \n",d_skip);
   #endif
  }

  /* --------- Accumulate Histories ------------------------ */

  ke_tot[hist_hi] = 0;
  for (isp=0; isp<nsp; isp++)  
  {
/* particle flux to outer wall */
/*    flux[isp][hist_hi] = n_lost[isp]*nc2p/(dt*sp_k[isp]*area1); */

    np_hist[isp][hist_hi] = np[isp]*nc2p/reactor_vol;      /* Simulation particle numbers */
    np_hist2[isp][hist_hi] = np[isp];

/* Average kinetic energy for each species */
    if (np[isp]) kes_hist[isp][hist_hi] = ke_sp[isp]*Escale[isp]/np[isp];

/* Total kinetic energy, summed over species */
    ke_tot[hist_hi] += ke_sp[isp]*Escale[isp];

  }
  t_array[hist_hi]= atime;
  com_cur_hist[hist_hi]= I_rf;                      /* Total Current */
  com_phi_hist[0][hist_hi] = phi[0];                /* LHS Potential */
  com_phi_hist[1][hist_hi] = phi[nc/2];             /* Mid Potential */
  com_power_hist[hist_hi] = power;                /* Instantaneous Power */
  avg_pow_hist[hist_hi] = pow_avg;                /* Instantaneous Power */
  sigma_hist[hist_hi] = sigma*sntos;              /* Charge density on LH electrode */
  ese_hist[hist_hi] = r0*sigma*phi[0]*fscale;
  for (i=0; i< ng; i++) 
    ese_hist[hist_hi] += vol[i]*rho[i]*phi[i]*fscale;
  Ez_hist[hist_hi] = E_z;                           /* Axial electric field */
  mu_hist[hist_hi] = mu_e;                                /* Electron mobility */
  
  if(fabs(E_z) > 1E-5)
    {
      mu_hist2[hist_hi] = ECHARGE*nc2p*np[0]*E_z*ef_norm;
      mu_hist2[hist_hi] = height*I_z/mu_hist2[hist_hi]; /* Alt electron mobility */      
    }
  else
    {
      mu_hist2[hist_hi] = 0;
    }
  /*mu_hist3[hist_hi] = mu_3;*/                                /* Electron mobility */
  Iz_hist[hist_hi] = I_z;                          /* Axial current - method 1 */
  Iz_hist2[hist_hi] = Iz_2;                        /* Axial current - method 2*/
/*  nu_iz_hist[hist_hi] = nu_iz/(sp_k[0]*dt)/np[0]; */ /* Ionization rate per electron */
/*   nu_ex_hist[hist_hi] = nu_ex/(sp_k[0]*dt)/np[0]; */ /* Excitation rate per electron */
/*   nu_el_hist[hist_hi] = nu_el/(sp_k[0]*dt)/np[0]; */ /* Elastic coll rate per electron */

  te_hist[hist_hi] = ke_tot[hist_hi]+ese_hist[hist_hi]; /* Total energy */
  /* ke_tot[hist_hi] = log10(fabs(ke_tot[hist_hi]+DBL_MIN));  */                 /* kinetic energy */
  /* ese_hist[hist_hi]= log10(fabs(ese_hist[hist_hi]+DBL_MIN));  */              /* field potential */

  C_hist[0][hist_hi]=C_hist[1][hist_hi]=0;                      /* LEEHJ */
  if (RT_flag) for (i=0;i<nc_RT;i++) 
  {
    temp=2.0*i+1;
    C_hist[0][hist_hi]+=sp_ex[i]*temp;
    C_hist[1][hist_hi]+=sp_exm[i]*temp;
  }
  temp=(float)(nc_RT); temp*=temp;
  C_hist[0][hist_hi]*=dntod/temp;
  C_hist[1][hist_hi]*=dntod/temp;

  /* ------ Calculate sheath width using ne/ni = 1/2                 */
 
  /* while ((sp_n[0][ii]/sp_n[1][ii] + sp_n[0][ii+1]/sp_n[1][ii+1] < 1.2) && (ii<nc/2)) ii++; */
  ii = sheath0();
  sheath[0][hist_hi] = ii;
  
  /* while ((sp_n[0][ii]/sp_n[1][ii]+ sp_n[0][ii-1]/sp_n[1][ii-1] < 1.2) && (ii>nc/2)) ii--; */

  ii = sheath1();
  sheath[1][hist_hi] = ii;

  hist_hi++;
  count = d_skip;
}

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

void freqanalysis(void)
{
  int i, j, init_freq_flag=1;
  static float *temp1, *temp2;
  
  if(init_freq_flag) 
  {
    if( !(temp1= (float  *)malloc((nfft/2)*sizeof(float)))
      || !(temp2= (float  *)malloc((nfft/2)*sizeof(float))))
       #ifdef MPI_1D
         if (my_rank == ROOT)
	   puts("Null ptr returned in freqanalysis()");
       #else
         puts("Null ptr returned in freqanalysis()");
       #endif
    init_freq_flag=0;
  }
  
   /******************************************************/
  for(i=0; i< nfft; i++)
    {
      cur_fft[i]= cur_hist[i];
      phi_fft[i]= phi_hist[0][i];
      mphi_fft[i]= phi_hist[1][i];
      pow_fft[i]= pow_hist[i];
    } 
  
  realft(phi_fft-1, freq_hi, 1);
  realft(mphi_fft-1,freq_hi, 1);
  realft(cur_fft-1, freq_hi, 1);
  realft(pow_fft-1, freq_hi, 1);
  
   /******************************************************/
  /**** Computing mag and phase of the current signal ***/

  temp2[0]= (cur_fft[0] > 0.0) ? 0.0 : 180.0;
  temp1[0]= fabs(cur_fft[0])*0.5/freq_hi;
/* Calculate amplitude = sqrt(A(2j)^2 + A(2j+1)^2) & phase = atan(A(2j)/A(2j+1))   */
  for(i=1, j=2; i< freq_hi; i++, j+=2)
    {
      temp1[i]= sqrt(cur_fft[j]*cur_fft[j]+ cur_fft[j+1]*cur_fft[j+1])/freq_hi;
      if(fabs(cur_fft[j+1]) < 1e-30 && fabs(cur_fft[j]) < 1e-30)
	temp2[i]= 0.0;
      else
	temp2[i]= (180./PI)*atan2(cur_fft[j], cur_fft[j+1]);
    }

/* Store amplitude and phase in curr_fft array */
  for(i=0; i< freq_hi; i++)
    {
      cur_fft[i]= temp1[i];
      cur_fft[freq_hi+i]= temp2[i];
    }
  
   /******************************************************/
  /**** Computing mag and phase of the voltage signal ***/
  temp2[0]= (phi_fft[0] > 0.0) ? 0.0 : 180.0;
  temp1[0]= fabs(phi_fft[0])/freq_hi/2;
  for(i=1, j=2; i< freq_hi; i++, j+=2)
    {
      temp1[i]= sqrt(phi_fft[j]*phi_fft[j]+ phi_fft[j+1]*phi_fft[j+1])/freq_hi;
      if(fabs(phi_fft[j+1]) < 1e-30 && fabs(phi_fft[j]) < 1e-30)
	temp2[i]= 0.0;
      else
	temp2[i]= (180./PI)*atan2(phi_fft[j], phi_fft[j+1]);
    }
  for(i=0; i< freq_hi; i++)
    {
      phi_fft[i]= temp1[i];
      phi_fft[freq_hi+i]= temp2[i];
    }
  
   /******************************************************/
  /**** Computing mag and phase of the mid-potential signal ***/
  temp2[0]= (mphi_fft[0] > 0.0) ? 0.0 : 180.0;
  temp1[0]= fabs(mphi_fft[0])/freq_hi/2;
  for(i=1, j=2; i< freq_hi; i++, j+=2)
    {
      temp1[i]= sqrt(mphi_fft[j]*mphi_fft[j]+ mphi_fft[j+1]*mphi_fft[j+1])/freq_hi;
      if(fabs(mphi_fft[j+1]) < 1e-30 && fabs(mphi_fft[j]) < 1e-30)
	temp2[i]= 0.0;
      else
	temp2[i]= (180./PI)*atan2(mphi_fft[j], mphi_fft[j+1]);
    }
  for(i=0; i< freq_hi; i++)
    {
      mphi_fft[i]= temp1[i] + DBL_MIN;
      mphi_fft[freq_hi+i]= temp2[i];
    }
  
   /******************************************************/
  /**** Computing mag and phase of the Plasma Impedance  ***/
  for(i=0; i< freq_hi; i++)
    {
      z_fft[i]= phi_fft[i] -cur_fft[i];
      z_fft[i+freq_hi]= phi_fft[i+freq_hi] -cur_fft[i+freq_hi];
      
      if(z_fft[i+freq_hi] > 180.0)
	z_fft[i+freq_hi] -= 360.0;
      
      else if(z_fft[i+freq_hi] < -180.0)
	z_fft[i+freq_hi] += 360.0;
    }
   /******************************************************/
  /**** Computing mag and phase of Plasma Power ***/
  temp2[0]= (pow_fft[0] > 0.0) ? 0.0 : 180.0;
  temp1[0]= fabs(pow_fft[0])/freq_hi/2;
  for(i=1, j=2; i< freq_hi; i++, j+=2)
    {
      temp1[i]= sqrt(pow_fft[j]*pow_fft[j]+ pow_fft[j+1]*pow_fft[j+1])/freq_hi;
      if(fabs(pow_fft[j+1]) < 1e-30 && fabs(pow_fft[j]) < 1e-30)
	temp2[i]= 0.0;
      else
	temp2[i]= (180./PI)*atan2(pow_fft[j], pow_fft[j+1]);
    }
  for(i=0; i< freq_hi; i++)
    {
      pow_fft[i]= temp1[i] + DBL_MIN;
      pow_fft[freq_hi+i]= temp2[i];
    }
  
}

/*-------------------------------------------------------------*/
void velocity(void)

/* This routine stores the velocity distribution at a specified position in the discharge */
{
   int i, isp, index, dum1;
   static int initflag=1;
   static float dv_r[NSMAX], roffset[NSMAX],dv_th[NSMAX], thoffset[NSMAX],
   dv_z[NSMAX], zoffset[NSMAX];

/********** Initialise velocity arrays **********/
   if(initflag) 
   {
      initflag =0;
      for(isp=0; isp<nsp; isp++) 
      {
         dv_r[isp] = (vr_u[isp] - vr_l[isp])/((nvrbin[isp]-1)*vscale);
         roffset[isp]= vr_l[isp]/vscale/dv_r[isp];
         dv_th[isp] = (vt_u[isp] - vt_l[isp])/((nvtbin[isp]-1)*vscale);
         thoffset[isp]= vt_l[isp]/vscale/dv_th[isp];
         dv_z[isp] = (vz_u[isp] - vz_l[isp])/((nvzbin[isp]-1)*vscale);
         zoffset[isp]= vz_l[isp]/vscale/dv_z[isp];
      }
      for(isp=0; isp<nsp;isp++)
      {
           for(i=0; i<nvrbin[isp]; i++) vr_array[isp][i] = vr_l[isp]+i*dv_r[isp]*vscale;
           for(i=0; i<nvtbin[isp]; i++) vth_array[isp][i] = vt_l[isp]+i*dv_th[isp]*vscale;
           for(i=0; i<nvzbin[isp]; i++) vz_array[isp][i] = vz_l[isp]+i*dv_z[isp]*vscale;
      }
   }

/********** Store velocities **********/

   for(isp=0;isp<nsp;isp++)
   {
      if (nvrbin[isp])
      for(i=0; i<np[isp];i++)
      { 
          index = (int)(v_r[isp][i]/dv_r[isp] - roffset[isp] +0.5);
          if((index>=0) && (index < nvrbin[isp]))
          {
/*              dum1 = (int) (r[isp][i] + 0.499); */
              dum1 = (int) (j_pointer(isp,i) + 0.499);
              vr_dist[isp][dum1][index] += 1;
          }
      }
      if (nvtbin[isp])
          for(i=0; i<np[isp];i++)
          {
              index = (int)(v_theta[isp][i]/dv_th[isp] - thoffset[isp] +0.5);
              if((index>=0) && (index < nvtbin[isp]))
              {
                /* dum1 = (int) (r[isp][i] + 0.499); */
                 dum1 = (int) (j_pointer(isp,i) + 0.499);
                 vth_dist[isp][dum1][index] += 1;
              }
          }
       if (nvzbin[isp])
           for(i=0; i<np[isp];i++)
           {
               index = (int)(v_z[isp][i]/dv_z[isp] - zoffset[isp] +0.5);
               if((index>=0) && (index < nvzbin[isp]))
               {
                  /* dum1 = (int) (r[isp][i] + 0.499); */
                   dum1 = (int) (j_pointer(isp,i) + 0.499);
                   vz_dist[isp][dum1][index] += 1;
               }
           }
   } /* end of for loop */
} /* end of function */

/*-------------------------------------------------------------*/
float mu_calc(void)
/* Calculate electron mobility using method 3 */
{
  int ii, ibin, ebin=300;
  static int init_flag=1;
  float ke_sum, energy, f_norm, df, mu;
  float static de_f, mu_norm, *eedf;

  if (init_flag)
  {
    eedf = (float *)malloc(ebin*sizeof(float *));
    mu_norm = -sqrt(2*ECHARGE/EMASS)/(3*gas_den);
    de_f = 30.0/ebin;
    init_flag = 0;
  }

/* zero eedf */
     for (ii=0; ii<ebin-1; ii++) eedf[ii] = 0;
     ke_sum =0;

/* Determine average eedf */
     for (ii=0; ii<np[0]; ii++)
     {
         energy = v_r[0][ii]*v_r[0][ii] + v_theta[0][ii]*v_theta[0][ii] + v_z[0][ii]*v_z[0][ii];
         ke_sum += energy;
         ibin = energy*Escale[0]/de_f;
         if (ibin < ebin) eedf[ibin] += 1;
         f_norm += 1;
     }
/* calc de for next time */
     de_f = ke_sum*Escale[0]/(30*np[0]);

/* normalise eedf */
     for (ii=0; ii<ebin; ii++) eedf[ii] /= sqrt((ii+0.5)*de_f);

     mu = 0.0;
     for (ii=0; ii<ebin-1; ii++)
     {
       energy = (ii+1)*de_f;
       df = eedf[ii+1] - eedf[ii];
       mu += energy*df/ (*sigmaMptr)(energy);
     }
     mu *= mu_norm/f_norm/de_f;  /* normalised mobility */

     return(mu);
}