pds1.c

美国berkeley大学开发的有界等离子体粒子1d3v计算机模拟程序,很实用

/* 	
**	PDS1 - ELECTROSTATIC 1 DIMENSIONAL BOUNDED PLASMA SIMULATION
**	
**	REVISION/PROGRAMMER/DATE
**	
**	1.0 modified version of PDP1 in cylindrical geometry
** /M. V. Alves and Vahid Vahedi/07-15-89
**
**	2.0 Addition of dumping and restoring a state of the system
**	/Vahid Vahedi and M. V. Alves/05-31-90
**
** 2.1 Fix the velocity scattering bug /Vahid Vahedi/03-01-92 
**
*/ 

#include <malloc.h>
#include "\win\wingraph.h"

#define GLOBALORIGIN
#include "pdefs.h"
#undef GLOBALORIGIN

FILE *InputDeck;
void start(), load();

main(argc, argv)
int argc;
char *argv[];
{
	int msg, flag=1;
	
	display_title();

	/*********************************************************/
	/* Read input file, and initialize params and vars. This */
	/*	function calls SPECIES and LOAD or Restore            */

	start(argc, argv);

	/*********************************************************/
	/* If collisions are included in the simulation load the */
	/* ion and electron-neutral collision cross-sections     */

	if(ecollisional) xsect();

	interval = 1;

	InputFile = argv[1];  			/* Passing the input deck name  */
											/* to the WinGraphics environ.  */
	InitWindows();	 					/* InitWindows()...				  */

	TimeStep = t;
	fields();							/* Compute electric potential	  */
											/*	and fields (includes CIRCUIT)*/
	history();							/* Initialize time histories    */

	while(msg != QUITMSG)
	{
		if(stepFlag) runFlag= TRUE;
	
      /* if(it== 12672 && flag)	runFlag= flag= 0;	 */
		if((np[0] >= MAXLEN-100 || np[1] >= MAXLEN-100) && flag)	runFlag= flag= 0; 
	
		if(runFlag)
		{
			t += dt;
			TimeStep = t;	 			/* Increment the time step */

			move();						/* Advance the particles one timestep  */

			adjust();					/* Adjust particle arrays to account
												for the several sources and boundary
												conditions.  This function uses others:
												VNEXT+FDIST+REVERS, REVERS, and GAUSIN. */
			fields();					/* Compute electric potential and fields 
												(includes CIRCUIT) */
			history();					/* accumulate and comb time histories */

			RefreshScreen();	   	/* Update diagnistics in the open windows */

			if(stepFlag)
			{
				stepFlag= FALSE;
				runFlag = FALSE;
			}
		}
		while(kbhit())					/* Window manager is invoked if a key */
		{									/*	is hit. It will then translate and */
			msg= getch();				/*	dispatch the given messages until  */
			msg= TranslateMsg(msg);	/*	the queue is empty */
			if(IsReadable(msg)) msg= DispatchMsg(msg);
		}
	}
}	/* End of MAIN */

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

void start(iargc, pargv)
int iargc;
char *pargv[];
{
	float revers();
	float jj0[NSMAX][2];
	char aachar[80];
	int i, isp, initnp[NSMAX][2];

	/***********************************************/
	/* Open files and read in "general" parameters */
	if (iargc < 2)
	{
		printf("no input deck? \n");
		exit(1);
	}
	InputDeck = fopen(pargv[1], "r");
	if (InputDeck == NULL)
	{
		/* Try to add default .INP */
		for (i=0; pargv[1][i] != 0; i++) aachar[i] = pargv[1][i];
		aachar[i++] = '.';
		aachar[i++] = 'i';
		aachar[i++] = 'n';
		aachar[i++] = 'p';
		aachar[i++] = 0;
		InputDeck = fopen(aachar, "r");
	}
	if (InputDeck==NULL)
	{
		printf("\nCan't find input file %s\n", pargv[1]);
		exit(1);
	}

	/* Read lines until we get to numbers */
	while (fscanf(InputDeck,"%d %d %f %f %f %f %f",
		&nsp, &nc, &nc2p, &dt, &r0, &r1, &epsilon) <7)
		fscanf(InputDeck, "%s", aachar);
		
	while (fscanf(InputDeck, "%f %f %f %f %f %f %f",
		&rhoback, &backj, &dde, &extr, &extl, &extc, &extq) <7)
		fscanf(InputDeck, "%s", aachar);

	while (fscanf(InputDeck, "%d %c %f %f %f %f %f",
		&dcramped, &src, &dcbias, &ramp, &acbias, &w0, &theta0) <7)
		fscanf(InputDeck, "%s", aachar);

	while (fscanf(InputDeck, "%d %d %d %d %d", &secsp, &ecolsp,
	   &icolsp, &reflux, &nfft) <5)
		fscanf(InputDeck, "%s", aachar);

	while (fscanf(InputDeck, "%f %f %d %f %f", &seec[0], &seec[1],
		&ionsp, &pressure, &gtemp) <5)
		fscanf(InputDeck, "%s", aachar);

	while (fscanf(InputDeck, "%f %f %f %f", &selsmax, &elsengy0,
		&elsengy1, &elsengy2) <4)
		fscanf(InputDeck, "%s", aachar);

	while (fscanf(InputDeck, "%f %f %f %f", &sextmax, &extengy0,
		&extengy1, &extengy2) <4)
		fscanf(InputDeck, "%s", aachar);

	while (fscanf(InputDeck, "%f %f %f %f", &sionmax, &ionengy0,
		&ionengy1, &ionengy2) <4)
		fscanf(InputDeck, "%s", aachar);

	while (fscanf(InputDeck, "%f %f %f %f", &achrgx, &bchrgx,
		&ascat, &bscat) <4)
		fscanf(InputDeck, "%s", aachar);

	/******************************************************/
	/* Check for errors in the "general" input parameters */
	
	if (extl < 0. || extr < 0. || extc < 0.)
	{
		printf("START: dont be cute with R,L,C < 0\n");
		exit(1);
	}
	if (!dt)
	{
		printf("START: dt == 0!!\n");
		exit(1);
	}
	if (nsp > NSMAX)
	{
		printf("START: nsp > NSMAX.\n");
		exit(1);
	}
	if (nc > NCMAX)
	{
		printf("START: nc > NCMAX.\n");
		exit(1);
	}
	
	/* Check to see if nfft is an integer power of 2 */
	if(nfft)
	{
		for(i=0; i<= 15; i++)
		{
			if(nfft == (1<<i)) 
			{
				i=0;
				break;
			}
		}
		if(i)
		{
			printf("START: nfft is not an integer power of 2.\n");
			exit(1);
		}
	}
	secondary = secsp;
	ecollisional= ecolsp;
	icollisional= icolsp;

	ng= nc+1;
	dr = (r1-r0)/nc; 
 												  
	/*************************************/
	/* Allocate space for field parameters */
	
	rg = (float far *)_fmalloc(ng*sizeof(float));
	rg3= (float far *)_fmalloc(ng*sizeof(float));
	r_array= (float far *)_fmalloc(ng*sizeof(float));
	for(i=0; i< ng; i++) r_array[i]= r0 +i*dr;
	srho= (float far **)malloc(nsp*sizeof(float far *));
	for(i=0; i<nsp; i++)	srho[i]= (float far *)_fmalloc(ng*sizeof(float));

	rho= (float far *)_fmalloc(ng*sizeof(float));
	phi= (float far *)_fmalloc(ng*sizeof(float));
	e  = (float far *)_fmalloc(ng*sizeof(float));

	/* Allocate space for xsections */
	if(ecollisional)
	{
		sels= (float far *)_fmalloc(NEMAX*sizeof(float));
		sext= (float far *)_fmalloc(NEMAX*sizeof(float));
		sion= (float far *)_fmalloc(NEMAX*sizeof(float));

		avene  = (float far *)_fmalloc(ng*sizeof(float));
		avenemp= (float far *)_fmalloc(ng*sizeof(float));
		ionrate= (float far *)_fmalloc(ng*sizeof(float));
		iontemp= (float far *)_fmalloc(ng*sizeof(float));
		for(i=0; i< ng; i++)	iontemp[i]= ionrate[i]= avene[i]= avenemp[i]= 0.;
	}
	
	if(icollisional)
	{
		chrgxrate= (float far *)_fmalloc(ng*sizeof(float));
		for(i=0; i< ng; i++) chrgxrate[i]= 0.;
	}

	/* Allocate space for time history diagnostics */

	t_array= (float far *)_fmalloc(HISTMAX*sizeof(float));
	ke_hist= (float far *)_fmalloc(HISTMAX*sizeof(float));
	te_hist= (float far *)_fmalloc(HISTMAX*sizeof(float));
	ese_hist= (float far *)_fmalloc(HISTMAX*sizeof(float));
	sigma_hist= (float far *)_fmalloc(HISTMAX*sizeof(float));
	
	com_cur_hist= (float far *)_fmalloc(HISTMAX*sizeof(float));
	com_phi_hist= (float far **)malloc(2*sizeof(float far *));
	com_phi_hist[0]= (float far *)_fmalloc(HISTMAX*sizeof(float));
	com_phi_hist[1]= (float far *)_fmalloc(HISTMAX*sizeof(float));
	
	Power_hist = (float far *)_fmalloc(HISTMAX*sizeof(float));
	
	np_hist= (float far **)malloc(nsp*sizeof(float far *));
	kes_hist= (float far **)malloc(nsp*sizeof(float far *));
	jwall_hist= (float far **)malloc(nsp*sizeof(float far *));
	
	for(i=0; i<nsp; i++)
	{
		np_hist[i]= (float far *)_fmalloc(HISTMAX*sizeof(float));
		kes_hist[i]= (float far *)_fmalloc(HISTMAX*sizeof(float));
		jwall_hist[i]= (float far *)_fmalloc(HISTMAX*sizeof(float));
	}
	
	/* Allocate space for frequency history diagnostics */

	if(nfft)
	{
		freq_hi= nfft/2;
		df= 1/(nfft*dt);
		f_array= (float far *)_fmalloc(nfft*sizeof(float));
		for(i=0; i< nfft; i++) f_array[i]= i*df;
		cur_hist= (float far *)_fmalloc(nfft*sizeof(float));
		phi_hist= (float far **)malloc(2*sizeof(float far *));
		phi_hist[0]= (float far *)_fmalloc(nfft*sizeof(float));
		phi_hist[1]= (float far *)_fmalloc(nfft*sizeof(float));

		Local_t_array = (float far *)_fmalloc(nfft*sizeof(float));

		z_fft= 	(float far *)_fmalloc(nfft*sizeof(float));
		cur_fft= (float far *)_fmalloc(nfft*sizeof(float));
		phi_fft= (float far *)_fmalloc(nfft*sizeof(float));
		mphi_fft=(float far *)_fmalloc(nfft*sizeof(float));
	}

	/* Allocate particle arrays */
	r  = (float far **)malloc(nsp*sizeof(float far *));
	vr = (float far **)malloc(nsp*sizeof(float far *));
	vth= (float far **)malloc(nsp*sizeof(float far *));
	vph= (float far **)malloc(nsp*sizeof(float far *));
	
	for(i=0; i<nsp; i++)
	{
		r[i]  = (float far *)_fmalloc(MAXLEN*sizeof(float));
		vr[i] = (float far *)_fmalloc(MAXLEN*sizeof(float));
		vth[i]= (float far *)_fmalloc(MAXLEN*sizeof(float));
		vph[i]= (float far *)_fmalloc(MAXLEN*sizeof(float));
	}
	
	if (nsp && (!r[nsp-1] || !vr[nsp-1] || !vth[nsp-1] || !vph[nsp-1]))
	{
		puts("START: r, v malloc failed");
		exit(1);
	}
	
	/**************************************/
	/* Conversions to dimensionless units */

	dttx = dt*dt/dr;
	vscale = dr/dt;

	gden = NperTORR*pressure/(gtemp+DBL_MIN);
	area = FOURPI*r0*r0;
	theta0 *= TWOPI/360.0;
   Trf=1/(w0+DBL_MIN);
	w0 *= TWOPI;
	epsilon*= EPS0;
	extq_1 = extq_2 = extq_3 = extq;

	if(dcramped)
	{
		if(ramp) risetime= fabs(dcbias)/ramp;
		else 		risetime= PI/w0;
	}
	
	ecolsp--;			 /* Fixing the indices into the array of species */
	ionsp--;
	secsp--;
	icolsp--;

	/* Set up each "species" parameters */
	for (isp=0; isp<nsp; ++isp)
	{
		species(initnp, jj0, isp);

		vrscale= max(vrscale, v0[isp][0]+vt[isp][0]);
		vrscale= max(vrscale, -v0[isp][1]+vt[isp][1]);

		/* Renormalize all velocities */
		v0[isp][0]/= vscale;
		v0[isp][1]/= vscale;
		vc[isp][0]/= vscale;
		vc[isp][1]/= vscale;
		vt[isp][0]/= vscale;
		vt[isp][1]/= vscale;