/* **	ES1 - ELECTROSTATIC 1 DIMENSIONAL PLASMA SIMULATION******	REVISION/PROGRAMMER/DATE****	1.0 Conversion from FORTRAN to C/T. Lasinski****	2.0 Conversion to MS C 5.0, improvement of speed and graphics resolution (to**	EGA), addition of diagnostics (E,P,KE,TE)/John P. Verboncoeur/12-12-87****	2.1 Dynamic array allocation/John P. Verboncoeur/02-23-88**	2.101 CGA code/John P. Verboncoeur/03-27-88**	2.102 Fix Maxwellian thermal loaders/John P. Verboncoeur/04-01-88**	2.103 Fix Magnetized velocity rotator/John P. Verboncoeur/04-18-88**	2.104 Register optimizations/John P. Verboncoeur/04-21-88**	2.105 Add kperp (charge disk)/John P. Verboncoeur/05-07-88**	2.106 Add Field Energy histories for Fourier modes 1-3/John P.**			Verboncoeur/05-08-88**	2.107 Enhance display paging scheme/John P. Verboncoeur/06-14-88**	2.108 Make '.INP' the default input deck extension/John P.**			Verboncoeur/06-18-88**	2.109 Enable arbitrary number of text lines in input decks/John P.**			Verboncoeur/06-20-88**	2.11  Semi-log plots/John P. Verboncoeur/07-29-88**	**	3.0 WinGraphics interface/John P. Verboncoeur and Vahid Vahedi/02-28-89****      4.0 XGRAFIX interface/Vahid Vahedi and John Verboncoeur/01-04-91****      4.1 Velocity distribution diagnostics, quadratic and cubic spline**          weighting for move and accel, with energy conserving algorithms,**	    Vx-Vy diagnostic.		**          Keith Cartwright and Peter Mardahl***/ #define GLOBALORIGIN#include "es1.h"#undef GLOBALORIGIN#include <xgrafix.h>main(argc, argv)int argc;char *argv[];{  int i;  char outfilename[100];  int J;  float Vmax;  FILE * outfile;  display_title();  /*********************************************************/  /* Read input file, and initialize params and vars. This */  /* function calls SPECIES and LOAD or Restore            */  XGInit(argc, argv, &t);  Start(argc,argv);  InitWindows(argc, argv);    /*  history();*/  XGStart();}void XGMainLoop(){  int i;  char outfilename[100];  int J;  float Vmax;  FILE * outfile;  /*********************************************************/  /* If collisions are included in the simulation load the */  /* ion and electron-neutral collision cros-sections     */    p = 0.0;    for (i=1; i <= nsp; i++)     {      accel(ins[i], ins[i+1]-1, qs[i], ms[i], ts[i], &pxs_hist[i], &kes_hist[i]);      p += pxs_hist[i];    }    for (i=1;i <= ng;i++) rho[i] = rho0;      rho[ng1] = 0.0;    rho[ng+2]=rho[0]=0.0;    /* ch:zeroing all extra points */        for (i=1; i <= nsp; i++) move(ins[i],ins[i+1]-1,qs[i]);    rho[1]+=rho[ng1];    rho[2]+=rho[ng+2];    rho[ng]+=rho[0];    rho[ng1]=rho[1];    rho[0]=rho[ng];    rho[ng+2]=rho[2];    t= it*dt;    history();    velocity();    fields(ith);        ith = ++it -ithl;    }/****************************************************************/display_title(){  printf("\n\nES1 - Electrostatic 1 Dimensional Code\n");  printf("version 4.1\n");  printf("(c) Copyright 1987-92 Regents of the University of California\n");  printf("Plasma Theory and Simulation Group\n");  printf("University of California - Berkeley\n");}/****************************************************************/Start(argc,argv) int argc;char *argv[];{  char a_char[80];  int i, j;    /**********************************************************/  /* Setting the global parameters to their default values. */    l=6.283185307;  dt=0.2;  nsp=1;    epsi=1.0;  ng=32;  iw=2;  ec=0;    ins[1] = 1;  vbins[1] = 1;  /*  sets up the start of vbin.  */  interval=1;  /***********************************************************/    if (!WasInputFileGiven) InputDeck = fopen("es1data","r");  else InputDeck = fopen(theInputFile,"r");  /*  if (argc >1 && !InputDeck)  {    for (i=0; argv[1][i] != 0; i++) a_char[i] = argv[1][i];    a_char[i++] = '.';		    a_char[i++] = 'i';    a_char[i++] = 'n';    a_char[i++] = 'p';    a_char[i++] = 0;    InputDeck = fopen(a_char,"r");  } */  if (!InputDeck)  {    printf("\nCan't find input file %s\n",argv[1]);    printf("\nCorrect syntax is: ES1 -i file.inp\n");    exit(-1);  }    /* read lines until we get to numbers */    while (fscanf(InputDeck,"%d %g %g %d %d %g %d",		&nsp, &l, &dt, &nt, &mmax, &la, &accum) <7)    fscanf(InputDeck, "%s", a_char);  /* note: la is l/a */    while (fscanf(InputDeck," %d %d %d %g %g %g %g %g",		&ng, &iw, &ec, &epsi, &a1, &a2, &e0, &w0) < 8)    fscanf(InputDeck, "%s", a_char);    if (nsp > NSPM) {    printf("Number of species nsp cannot exceed NSPM\n");    exit(-1);  }  if(accum<0) { printf("\nError:  accum can't be negative! \n"); exit(1);	      }  if(!(ec==1||ec==0)) {    printf("\n Error:  What are you thinking?  There are only two possible values of ec\n");    printf("0 and 1.  %d is not 0 or 1.",ec);    exit(1);  }  if(!iw&&ec) {    printf("\nError:  There IS no energy-conserving algorithm for NGP\n");    exit(1);  }  ecconst=0.0;  if (ec) {           ecconst=0.5;  }  if(iw>3 || iw<0)    {      printf("\nError:  bad iw flag!  Please check your input deck!\n");      exit(1);    }  if (ng > NGMAX) {    printf("\nNumber of grids ng cannot exceed NGMAX\n");    exit(-1);  }  dx = l/ng;  ng1= ng+1;  k_hi= ng/2;     /*******************************/  /* Allocating space for arrays */    nms= (float *)malloc((nsp+1)*sizeof(float));  ms = (float *)malloc((nsp+1)*sizeof(float));  qs = (float *)malloc((nsp+1)*sizeof(float));  ts = (float *)malloc((nsp+1)*sizeof(float));    x_array= (float *)malloc((ng+1)*sizeof(float));  for (i=0; i<= ng; i++) x_array[i] = i*dx;    rho= (float *)malloc((ng+3)*sizeof(float));   phi= (float *)malloc((ng+2)*sizeof(float));   phik=(float *)malloc((ng+2)*sizeof(float));    k_array= (float *)malloc(ng*sizeof(float));  for(i=0; i< k_hi; i++) k_array[i]= i*2*PI/l;    e  = (float *)malloc((ng+2)*sizeof(float));  acc= (float *)malloc((ng+3)*sizeof(float));   t_array= (float *)malloc(HISTMAX*sizeof(float));  ese= (float *)malloc(HISTMAX*sizeof(float));  ke = (float *)malloc(HISTMAX*sizeof(float));  te = (float *)malloc(HISTMAX*sizeof(float));    kes_hist= (float *)malloc((nsp+1)*sizeof(float));  pxs_hist= (float *)malloc((nsp+1)*sizeof(float));  esem_hist=(float *)malloc((mmax+1)*sizeof(float));    kes = (float **)malloc(nsp*sizeof(float *));  for (i=0; i< nsp; i++)     kes[i] = (float *)malloc(HISTMAX*sizeof(float));    pxs = (float **)malloc(nsp*sizeof(float *));  for (i=0; i< nsp; i++)     pxs[i] = (float *)malloc(HISTMAX*sizeof(float));    esem = (float **)malloc(mmax*sizeof(float *));  for (i=0; i< mmax; i++)     esem[i] = (float *)malloc(HISTMAX*sizeof(float));    x  = (float *)malloc(MAXPARTICLES*sizeof(float));  vx = (float *)malloc(MAXPARTICLES*sizeof(float));  vy = (float *)malloc(MAXPARTICLES*sizeof(float));    vbint= (float *)malloc(NVBINMAX*sizeof(float));  vbin = (float *)malloc(nsp*NVBINMAX*sizeof(float));  vbin_inst= (float *)malloc(nsp*NVBINMAX*sizeof(float));  dvbin = (float *)malloc((nsp+1)*sizeof(float));  vbinstart = (float *)malloc((nsp+1)*sizeof(float));  v_array = (float *)malloc( NVBINMAX*nsp * sizeof(float));  for(i=0; i<nsp*NVBINMAX; i++) vbin_inst[i]= 0.0;  if (!x || !vx || !vy || !vbint || !vbin_inst || !dvbin || !v_array )  {    printf("START: Could not get enough memory for x or v's.\n");    exit(1);  }    printf(" nsp = %2d     l = %8.5f \n",nsp,l);  printf(" dt = %4.2f    nt = %4d \n",dt,nt);  printf(" ng = %5d   iw = %2d   ec = %2d  accum = %4d\n",ng,iw,ec,accum);  printf(" epsi = %4.2f  a1 = %4.2f  a2 = %4.2f \n",epsi,a1,a2);    for (i=1;i <= nsp;i++)     init(&ins[i],&ins[i+1],&ms[i],&qs[i],&ts[i],&nms[i],&vbins[i],&vbins[i+1],&dvbin[i],&vbinstart[i],&nvbin[i]);  /* added vbins to param list */  fclose(InputDeck);    for (i=1; i<= nsp; i++) np[i]= ins[i+1] -ins[i];    for (i=0; i < nsp; i++)  {    for (j=0; j < HISTMAX; j++)     {      kes[i][j] = 0.0;      pxs[i][j] = 0.0;    }  }  rho[1]+=rho[ng1];   /*  These resolve the periodic boundary conditions.  */  rho[2]+=rho[ng+2];  rho[ng]+=rho[0];  rho[0]=rho[ng];  rho[ng+2]=rho[2];  rho[ng1] = rho[1];  fields(0);    for (i=1; i <= nsp; i++) {    setv(ins[i], ins[i+1]-1, qs[i], ms[i], ts[i], &pxs_hist[i], &kes_hist[i]);    /*  scale all the velocities  properly */    dvbin[i] *= dt/dx;    vbinstart[i] *= dt/dx;    }  startvel();}/**************************************************************/history(){  int i, j, k;  static int count=1;    if (--count) return;				/* only accum every interval steps */  if (hist_hi >= HISTMAX)			/* comb time histories */  {    for (j=0; j< nsp; j++)    {      for (i=1, k=4; i<HISTMAX/4; i++, k+=4)      {	kes[j][i] = kes[j][k];	pxs[j][i] = pxs[j][k];      }    }    for (j=0; j< mmax; j++)    {      for (i=1, k=4; i<HISTMAX/4; i++, k+=4)	esem[j][i] = esem[j][k];    }    for (i=1, k=4; i<HISTMAX/4; i++, k+=4)    {      ese[i] = ese[k];      ke[i] = ke[k];       te[i] = te[k];      t_array[i] = t_array[k];    }    hist_hi = i;    interval *= 4;  }  t_array[hist_hi]= t;  for (i=0; i< mmax; i++) esem[i][hist_hi] = fabs(esem_hist[i+1]) +1e-30;    ke[hist_hi]= 1e-30;  for (j=0; j< nsp; j++) {    kes[j][hist_hi] = kes_hist[j+1];    ke[hist_hi] += kes_hist[j+1];    pxs[j][hist_hi] = pxs_hist[j+1];  }  te[hist_hi] = ke[hist_hi] + ese_hist;  ese[hist_hi]= ese_hist+ 1e-30;  hist_hi++;  count = interval;}/**********************************************************************  velocity()  This code maintains the velocity distribution code.    It works thusly:  I wanted to accumulate velocities in bins for  ten iterations, then slap up a velocity distribution to be seen  by the users of xes1.  To do that I had to play the following tricks:   I do all the work in vbin_inst, and I'm always displaying vbin.  When I've accumulated velocities for accum timesteps  I copy the current accumulated vbin_inst into vbin.**************************************************************/velocity(){  int  s,nbinmaxi;  register int i,j;  static int count=5;  float pval;  float *vsps; 	  /*  vsps is a pointer to the part of the		      velocity bin array where the species 		      starts  */    float vst,dvt;		/*  this is the vstart[isp],dv[isp]  */    if(!accum) return; /* if the accum is set to zero, don't bother making			a current velocity distribution profile.  */  pval=1.0/accum;  /*  This code does the velocity-distribution stuff.  */  count--;  if(!count ) {    for (j=0; j<=nvbin[1];j++) vbint[j]=0.0;    for(i=1;i<=nsp;i++) {      for (j=vbins[i];j<=vbins[i]+nvbin[i];j++) {	vbin[j]=ms[i]*vbin_inst[j];	vbin_inst[j]=0.0;      }      for (j=0; j<=nvbin[1];j++) vbint[j]+= vbin[j+vbins[i]];    }  }  if( !(count)) count=accum;  for(i=1;i<=nsp;i++) {    vsps = & (vbin_inst[vbins[i]] );    vst=vbinstart[i];    dvt=dvbin[i];    nbinmaxi=nvbin[i];        if(nbinmaxi>0)      for (j=ins[i];j<ins[i+1];j++) {	s= (vx[j]-vst)/dvt;	if(s>=0 && s<nbinmaxi) vsps[s]+=pval;      }	  }}  /************************************************************************//* frand() returns values 0 through 1.                                  *//* From "Random number generators: good ones are hard to find", S. Park *//* and K. Miller, Communications of ACM, October 1988, pp 1192-1201.    *//* This is from page 1195, and is to work on any system for which       *//* maxint is 2**31-1 or larger. Due earlier to Schrage, as cited by P&M.*//*                                                                      *//* Note: OK values for iseed are 1 through 2**31-2. Give it 0 or 2*31-1 *//* and it will return the same values thereafter!                       *//*                                                                      *//* C version 6/91, Bruce Langdon.                                       *//*                                                                      *//* Algorithm replaces seed by mod(a*seed,m). First represent            *//* seed = q*hi + lo.  Then                                              *//* a*seed = a*q*hi + lo = (m - r)*hi + a*lo = (a*lo - r*hi) + m*hi,     *//* and new seed = a*lo - r*hi unless negative; if so, then add m.       */float frand(){  long a = 16807, m = 2147483647, q = 127773, r = 2836;  long hi, lo;  static long seed=31207321;  float fnumb;  hi = seed/q;  lo = seed - q*hi;  seed = a*lo - r*hi;  /* "seed" will always be a legal integer of 32 bits (including sign). */  if(seed <= 0) seed = seed + m;  fnumb = seed/2147483646.0;  return(fnumb);}/***********************************************************************/