initwin.c

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

#include "pdc1.h"

#include <xgrafix.h>

/*********************************************************/
/* void XGWrite(float *, int, int, char *, char); */
/* void XGRead(float *, int, int, char *, char); */
void Dump(char *);
void Restore(char *);
void Quit(void);

void InitWindows(void)
{
  int ii, isp;
  char buffer[50];

  /*********************************************/
  
/*-------- Vr - r phase space --------*/
  for (isp=nsp-1; isp>=0; isp--) 
  {
     sprintf(buffer, "Vr-R Phase Space %d", isp+1);
     XGSet2D("linlin", "R", buffer, "closed", 300, 10, dr, dr/dt,
	  False, False, r0*dr, r1*dr, vscale*(v0[isp][1]+5*vt[isp][1]), 
          vscale*(v0[isp][0]+5*vt[isp][0]));
     XGScat2D(r[isp], v_r[isp], &np[isp], isp);
  }
/*-------- Vt - r phase space --------*/
  for (isp=nsp-1; isp>=0; isp--) 
  {
     sprintf(buffer, "Vtheta-R Phase Space %d", isp+1);
     XGSet2D("linlin", "R", buffer, "closed", 300, 10, dr, dr/dt,
	  False, False, r0*dr, r1*dr, vscale*(-v0_th[isp]-5*vt_th[isp]), 
          vscale*(v0_th[isp]+5*vt_th[isp]));
     XGScat2D(r[isp], v_theta[isp], &np[isp], isp);
  }
/*-------- Vz - r phase space --------*/
  for (isp=nsp-1; isp>=0; isp--) 
  {
     sprintf(buffer, "Vz-R Phase Space %d", isp+1);
     XGSet2D("linlin", "R", buffer, "closed", 300, 10, dr, dr/dt,
	  False, False, r0*dr, r1*dr, vscale*(-v0_z[isp]-5*vt_z[isp]), 
          vscale*(v0_z[isp]+5*vt_z[isp]));
     XGScat2D(r[isp], v_z[isp], &np[isp], isp);
  }
/*-------- Vtr - r phase space --------*/
  for (isp=nsp-1; isp>=0; isp--) 
  {
     sprintf(buffer, "rVt-R Phase Space %d", isp+1);
     XGSet2D("linlin", "R", buffer, "closed", 300, 10, dr, dr/dt,
	  False, True, r0*dr, r1*dr, vscale*(v0_th[isp]+5*vt_th[isp]), 
          vscale*(v0_th[isp]+5*vt_th[isp]));
     XGScat2D(r[isp], angular_mtm[isp], &np[isp], isp);
  }
/*-------- Vz - r phase space --------*/
  /*********************************************/

 XGSet2D("linlin", "Time", "Avg Density(t)", "closed", 100, 400, 1.0, 1.0, 
         True, True,1.0,1.0,1.0,1.0);

  for (isp=nsp-1; isp>=0; isp--) XGCurve(t_array, np_hist[isp], &hist_hi, isp);
  for (isp=0;isp<2;isp++) XGCurve(t_array, C_hist[isp], &hist_hi, isp+2);

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

  XGSet2D("linlin", "Time", "Computational Particle Number(t)", "closed", 100, 400, 1.0, 1.0, 
         True, True,1.0,1.0,1.0,1.0);

  for (isp=nsp-1; isp>=0; isp--) XGCurve(t_array, np_hist2[isp], &hist_hi, isp);

  /*********************************************/
  
  XGSet2D("linlin", "R", "E field(r)", "closed", 500, 10, 1.0, 1.0/(ptopn*dr), 
           False, True, r0*dr, r1*dr,1.0,1.0);
  
  /* field from Gauss's law integration */
  XGCurve(r_array, efield, &ng, 3);

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

  XGSet2D("linlin", "R", "Potential(r)", "closed", 400, 300, 1.0, 1.0/ptopn, 
           False, True, r0*dr,r1*dr,1.0,1.0);
  
  /* potential from Gauss's law calculation */
  XGCurve(r_array, phi, &ng, 4);
	
  /*********************************************/
  /* Time-dependance axial E field */
  if (src =='J' || src=='j' || src=='P'|| src=='p') 
  {
  XGSet2D("linlin", "Time", "E field(z)", "closed", 500, 10, 1.0, 1.0/(ptopn*dr),
           True, True, 1.0, 1.0,1.0,1.0);

  XGCurve(t_array, Ez_hist, &hist_hi, 3);
  }

  /*********************************************/
  /* Time-dependance of electron mobility */
  XGSet2D("linlin", "Time", "mu_e (t) ", "closed", 500, 10, 1.0, 1.0,
           True, True, 1.0, 1.0,1.0,1.0);

  XGCurve(t_array, mu_hist,  &hist_hi, 4);
  XGCurve(t_array, mu_hist2, &hist_hi, 3);
  /*XGCurve(t_array, mu_hist3, &hist_hi, 5);*/
  /*********************************************/
  XGSet2D("linlin", "Time", "I (z)", "closed", 500, 10, 1.0, 1.0,
           True, True, 1.0, 1.0,1.0,1.0);

     XGCurve(t_array, Iz_hist, &hist_hi, 4);
     XGCurve(t_array, Iz_hist2, &hist_hi, 6);
  /*********************************************/
/*
  XGSet2D("linlin", "Time", "Ionization rate", "closed", 500, 400,
           1.0, 1.0, True, True,1.0,1.0,1.0,1.0);

  XGCurve(t_array, nu_iz_hist, &hist_hi, 4);
*/
  /*********************************************/
/*
  XGSet2D("linlin", "Time", "Excitation rate", "closed", 500, 400,
           1.0, 1.0, True, True,1.0,1.0,1.0,1.0);

  XGCurve(t_array, nu_ex_hist, &hist_hi, 4);
*/
  /*********************************************/
/*
  XGSet2D("linlin", "Time", "Elastic scatt rate", "closed", 500, 400,
           1.0, 1.0, True, True,1.0,1.0,1.0,1.0);

  XGCurve(t_array, nu_el_hist, &hist_hi, 4);
*/
  /*********************************************/

  XGSet2D("linlin", "Time", "Mid Potential(t)", "closed", 500, 400, 
           1.0, 1.0/ptopn, True, True,1.0,1.0,1.0,1.0);
  
  XGCurve(t_array, com_phi_hist[1], &hist_hi, 4);

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

  XGSet2D("linlin", "Time", "LHS Potential(t)", "closed", 700, 400, 
           1.0, 1.0/ptopn, True, True,1.0,1.0,1.0,1.0);
  
  XGCurve(t_array, com_phi_hist[0], &hist_hi, 4);

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

if (src =='I' || src=='i' || src=='V'|| src=='v' || src=='K' || src=='W') {
  XGSet2D("linlin", "Time", "Current(t)", "closed", 600, 400, 1.0, 1.0, 
           True, True,1.0,1.0,1.0,1.0);
  
  XGCurve(t_array, com_cur_hist, &hist_hi, 4);

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

  XGSet2D("linlin", "Time", "Sigma(t)", "closed", 500, 500, 1.0, 1.0, 
           True, True,1.0,1.0,1.0,1.0);
  
  XGCurve(t_array, sigma_hist, &hist_hi, 3);
  
  /*********************************************/

  XGSet2D("linlin", "Time", "Power(t)", "close", 300, 550, 1.0, 
    	    1.0/ptopn, True, True, 0.0, 0.0, 0.0, 0.0);

  XGScat2D(t_array, com_power_hist, &hist_hi, 1);
  XGCurve(t_array, avg_pow_hist, &hist_hi, 2);
}
  /*********************************************/

 /* XGSet2D("linlin", "Time", "Flux to wall", "closed", 600, 200, 
           1.0, 1.0, True, True, 0.0, 0.0, 0.0, 0.0);
 
  for (isp=nsp-1; isp>=0; isp--) XGCurve(t_array, flux[isp], &hist_hi, isp);*/

  /*********************************************/
  /* Charge density profile */

  XGSet2D("linlin", "R", "rho(r) filt,unfilt", "closed", 400, 100, 1.0, 1.0, 
          False, True,r0*dr,r1*dr,1.0,1.0);

  XGCurve(r_array, rho, &ng, 2);
  XGCurve(r_array, rho_unfilt, &ng, 3);

  /*********************************************/
  /* Particle density profile */
 
  XGSet2D("linlin", "R", "Density n(r)", "closed", 600, 500, 1.0, 1.0,
	  False, True, r0*dr,r1*dr,1.0,1.0);

  for (isp=nsp-1; isp>=0; isp--) XGCurve(r_array, density[isp], &ng, isp);

  /*********************************************/
  /* J.E energy profile */
 
  XGSet2D("linlin", "R", "Ion J-E (r)", "closed", 600, 500, 1.0, 1.0,
	  False, True, r0*dr,r1*dr,1.0,1.0);

  for (isp=1; isp<nsp; isp++) XGCurve(r_array, jdote[isp], &ng, isp);

  /*********************************************/
  /* Jz.Ez energy profile */
 
  XGSet2D("linlin", "R", "Jz-Ez (r)", "closed", 600, 600, 1.0, 1.0,
	  False, True, r0*dr,r1*dr,1.0,1.0);

  for (isp=0; isp<nsp; isp++) XGCurve(r_array, jzdotez[isp], &ng, isp);

  /*********************************************/
  /* Time-averaged Current Density profile */

  XGSet2D("linlin", "R", "Ion J (r)", "closed", 600, 500, 1.0, 1.0,
          False, True, r0*dr,r1*dr,1.0,1.0);

  for (isp=1; isp<nsp; isp++) 
    XGCurve(r_array, jr_prof[isp], &ng, isp);
  
  /*********************************************/
  /* Electron & ion mobility profile */

  XGSet2D("linlin", "R", "mu_e (r) ", "closed", 400, 100, 1.0, 1.0,
          False, True,r0*dr,r1*dr,1.0,1.0);

   for (isp=0; isp<nsp; isp++) XGCurve(r_array, mu_prof[isp], &ng, isp);

  /*********************************************/
  /* Axial current density profile */
  XGSet2D("linlin", "R", "J_z (r)", "closed", 600, 500, 1.0, 1.0,
          False, True, r0*dr,r1*dr,1.0,1.0);

  for (isp=0; isp<nsp; isp++) XGCurve(r_array, jz_prof[isp], &ng, isp);
  /* for (isp=0; isp<nsp; isp++) XGCurve(r_array, jz_unfilt[isp], &ng, isp+2); */

  /*********************************************/
  /* Time-averaged KE profile */

  XGSet2D("linlin", "R", "KE (r)", "closed", 600, 500, 1.0, 1.0,
          False, True, r0*dr,r1*dr,1.0,1.0);

  for (isp=0; isp<nsp; isp++) XGCurve(r_array, ke_avg[isp], &ng, isp);

/*----------------------------------*/

  XGSet2D("linlin", "R", "KE r (r)", "closed", 600, 500, 1.0, 1.0,
          False, True, r0*dr,r1*dr,1.0,1.0);

  for (isp=0; isp<nsp; isp++) XGCurve(r_array, ker_avg[isp], &ng, isp);

/*----------------------------------*/

  XGSet2D("linlin", "R", "KE theta (r)", "closed", 600, 500, 1.0, 1.0,
          False, True, r0*dr,r1*dr,1.0,1.0);

  for (isp=0; isp<nsp; isp++) XGCurve(r_array, ket_avg[isp], &ng, isp);

/*----------------------------------*/

  XGSet2D("linlin", "R", "KE z (r)", "closed", 600, 500, 1.0, 1.0,
          False, True, r0*dr,r1*dr,1.0,1.0);

  for (isp=0; isp<nsp; isp++) XGCurve(r_array, kez_avg[isp], &ng, isp);

  /*********************************************/
  /*     Excited State Diagnostics             */
  /*********************************************/
if (RT_flag) {
  XGSet2D("linlin", "R", "Density_r (r)", "closed", 600, 500, Dr,
            dntod, False, True, 0.0, r1*dr, 0.0, 0.0);
  XGCurve(r_grid2, sp_ex, &nc_RT, 3);

  XGSet2D("linlin", "R", "Density_m (r)", "closed", 600, 500, Dr,
            dntod, False, True, 0.0, r1*dr, 0.0, 0.0);
  XGCurve(r_grid2, sp_exm, &nc_RT, 2);

  XGSet2D("linlin", "R", "Emission coeff (r)", "closed", 600, 500, Dr,
           F_factor, False, True, 0.0, r1*dr, 0.0, 0.0);
  XGCurve(r_grid2, emission_coeff, &nc_RT, 3);

  XGSet2D("linlin", "R", "Radiation Flux (r)", "closed", 600, 500, Dr,
           2*M_PI*Dr*F_factor, False, True, 0.0, r1*dr, 0.0, 0.0);
  XGCurve(r_grid2, radiation_flux, &nc_RT, 3);

  XGSet3D("linlinlin","i","j","A_ij", 45, -45, "closed", 600, 500, 1, 1, 1,
          False, False, True, 0, nc_RT, 0, nc_RT, 0, 0);
  XGSurf(r_grid, r_grid, A, &nc_RT, &nc_RT, 3);

/*XGSet3D("linlinlin", "Time", "R", "time ave N_ex(r,t)", 45.0, 225.0,
          "closed", 100, 150, 1.0, Dr, dntod, True, False, True,
          0.0, 0.0, 0.0, r1*dr, 0.0, 0.0);
  XGSurf(t_array2, r_grid2, st_nex, &hist_hi2, &nc_RT, 3);       */

  XGSet2D("linlin", "R", "Effective production (r)", "closed", 600, 500, Dr,
           dntod/(dt*sp_k_ex), False, True, 0.0, r1*dr, 0.0, 0.0);
  XGCurve(r_grid2, E_prod[0], &nc_RT, 3);

  XGSet2D("linlin", "R", "Effective production2 (r)", "closed", 600, 500, Dr,
           dntod/(dt*sp_k_ex), False, True, 0.0, r1*dr, 0.0, 0.0);
  for (isp=1;isp<4;isp++) XGCurve(r_grid2, E_prod[isp], &nc_RT, isp);

  XGSet2D("linlog", "Time", "Error_sum (t)", "closed", 100, 400, 1.0,
           1.0/nc_RT, True, False, 0.0, 0.0, 1e-7, 0.1);
  XGCurve(t_array3, error_sum, &hist_hi3, 2);
 
  XGSet2D("linlin", "R", "Error (r)", "closed", 700, 500, Dr,
            1, False, True, 0.0, r1*dr, 0.0, 0.0);
    XGCurve(r_grid2, error, &nc_RT, 3);

  XGSet2D("linlin", "Time", "Radiation Flux (t)", "closed", 600, 500, 1.0,
           2*M_PI*Dr*F_factor, True, True, 0.0, 0, 0.0, 0.0);
  XGCurve(t_array3, tot_intensity, &hist_hi3, 0);   

  XGSet2D("linlin", "normalized f", "Lineshape (f)", "closed", 600, 500, 1,
            1, False, True, xmin, xmax, 0.0, 0.0);
    XGCurve(x_array, L_array, &xbin, 3);

}
  /*********************************************/
  /* Distribution functions                    */
  /*********************************************/

  for(isp=0; isp<nsp; isp++)
  {
    sprintf(buffer, "f(E) %d", isp+1);

    XGSet2D("linlin", "E", buffer, "closed", 650, 350, 
         Escale[isp], 1.0, False, True, e_array[isp][0]*Escale[isp],
	 e_array[isp][nbin[isp]-1]*Escale[isp], 0.0, 0.0);
    XGCurve(e_array[isp], fe[isp], &nbin[isp], isp);

  }
  for(isp=0; isp<nsp; isp++)
  {
    sprintf(buffer, "f(E) mid %d", isp+1);

    XGSet2D("linlin", "E", buffer, "closed", 650, 350, 
         Escale[isp], 1.0, False, True, emid_array[isp][0]*Escale[isp],
         emid_array[isp][nbin_mid[isp]-1]*Escale[isp], 0.0, 0.0);
    XGCurve(emid_array[isp], fe_mid_show[isp], &nbin_mid[isp], isp);
  }

/*********************************************/
/* velocity distribution */
   if (vel_dist_accum)
   {
     for(isp=0; isp<nsp; isp++)
     {
      if (nvrbin[isp])
      {
          sprintf(buffer, "V_r dist species %d", isp+1);
          XGSet3D("linlinlin", "R", "Vr", buffer, -100.0, 40.0, "closed",
                700, 175, dr, 1, 1.0, True, True, True, 0, 1.0, 0, 1.0, 0, 1.0);
          XGSurf(r_grid, vr_array[isp], vr_dist[isp], &ng, &nvrbin[isp], 4);
      }
      if (nvtbin[isp])
      {