pds1.c

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

		vtpt[isp] /= vscale;

		Escale[isp]= .5*vscale*vscale*m[isp]/(1.602e-19);
		emin[isp] /= Escale[isp];
		de[isp]   /= Escale[isp];
		
		/* Conversions from physical to code units */
		qm[isp]= q[isp]/m[isp];
		q[isp] *= nc2p;
		
		/* Normalize  "enter" to = no. particles injected PER DT */
		if(enter[isp][0] = jj0[isp][0]*FOURPI*r0*r0*dt/fabs(q[isp]))
			inject[isp]= 1;
		if(enter[isp][1] = jj0[isp][1]*FOURPI*r1*r1*dt/fabs(q[isp]))
			inject[isp]= 1;
	
	}	/* end  "for isp=0 thru nsp-1 "  loop */

	if(!vrscale) vrscale = 1e4;
	vrscale *= 5;
	
	/*************************************************/
	/* Allocate space for the dist. functions arrays */
	
	fe      = (float far **)malloc(nsp*sizeof(float far *));
	e_array = (float far **)malloc(nsp*sizeof(float far *));
	ftheta  = (float far **)malloc(nsp*sizeof(float far *));
	th_array= (float far **)malloc(nsp*sizeof(float far *));
	
	for(isp=0; isp<nsp; isp++)
	{
		fe[isp]     = (float far *)_fmalloc(nbin[isp]*sizeof(float));
		ftheta[isp]  = (float far *)_fmalloc(nbin[isp]*sizeof(float));
		e_array[isp]= (float far *)_fmalloc(nbin[isp]*sizeof(float));
		th_array[isp]= (float far *)_fmalloc(nbin[isp]*sizeof(float));
	}

	if (nsp && (!th_array[nsp-1]))
	{
		puts("START: f(E) malloc failed");
		exit(1);
	}

	for(isp=0; isp<nsp; isp++)
	{
		for(i=0; i<nbin[isp]; i++)
		{
			fe[isp][i]      = 0.0;
			ftheta[isp][i]  = 0.0;
			e_array[isp][i] = emin[isp] +i*de[isp];
		   th_array[isp][i]= i*dtheta[isp];
		}
	}

	/*****************************************/
	/* Normalization of the input parameters */

	r0 /=dr;
	r1 /=dr;
	for ( i=0; i<= nc; i++)
	{
		rg[i] =r0 +i;
		rg3[i]= (rg[i]+.5)*(rg[i]+.5)*(rg[i]+.5);
	}
 
	/********************************************************/
	/* Load in all species' particles, properly distributed */
	
	if (iargc>2 && access(pargv[2],4))
	{
		/* Try to add default .DMP */
		for (i=0; pargv[2][i] != 0; i++) aachar[i] = pargv[2][i];
		aachar[i++] = '.';
		aachar[i++] = 'D';
		aachar[i++] = 'M';
		aachar[i++] = 'P';
		aachar[i++] = 0;
		pargv[2] = aachar;
	}

	/*********************************************************/
	/* Set up array of bit-reversed indices. "revers(num)"	*/
	/*	reverses the order of the base-i representation	of		*/
	/* the number "num",	yielding a fraction between 0 and 1.*/
	/* Index[] is also used in padjus.c.							*/
	
	for (i=0; i<1024; i++) index[i] = 1024*revers(i);

	/**************************************************/
	/* Start from a dump file if given otherwise load */
	/*	the system with the initial distribution		  */

	if (iargc>2 && !access(pargv[2],4))	Restore(pargv[2]);
	else load(initnp);

	/**************************************************************/
	/* Deciding which circuit solver to use  */

	if(extc < 1e-30)   /* When C tends to 0.0 */
	{
		if(fabs(dcbias)+fabs(ramp)+fabs(acbias) >0.0)
	 		printf("START:The active external source is ignored since C<1e-30\n");
	 
		(*bcptr) = bc1;
		(*circuitptr) = circuit1;
		ntri =1;  
		b0 = 1.;
	} 		
	else if(src== 'I' || src== 'i')	/* When current source is applied */
	{
		(*bcptr)= bc2;
		(*circuitptr)= circuit2;

		ntri =1;
		b0 = 1.;
	}
	else if(extl <1e-30 && extr <1e-30 &&
           (extc*(r1-r0)/(FOURPI*r0*r1*epsilon)) >= 1e5) 
	{	 /* When L=R=0. and C tends to infinite */
		(*bcptr) = bc3;
		(*circuitptr) = circuit3;
		ntri =2;
		b0 = -((rg[1]+0.5)*(rg[1]+0.5)+(r0+0.5)*(r0+0.5));
   }
	else  /* The general case with external voltage source */
	{ 
		(*bcptr) = bc4;	  
		(*circuitptr) = circuit4;
		ntri =1; 
		b0= 1/(2.25*extl/dt/dt +1.5*extr/dt +1/extc);
		b0 = 1.+r0*r0*dr*b0/((r0+0.5)*(r0+0.5)*epsilon*area);
   }
}	/* End START */

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

species(initnp, jj0, isp)
float jj0[][2];
int isp, initnp[][2];
{
	char aachar[80];
	float iinitn, j0l, j0r, v0l, v0r, vcl, vcr, vtl, vtr,
		fluxl, fluxr, temin, temax;

	/* Read in SPECIES parameters... */
	while (fscanf(InputDeck, "%f %f %f %f %f", &q[isp], &m[isp],
		&j0l, &j0r, &iinitn) <5)
		fscanf(InputDeck, "%s", aachar);

	while (fscanf(InputDeck, "%f %f %f %f %f %f %f", &v0l, &v0r, &vtl,
		&vtr, &vcl, &vcr, &vtpt[isp]) <7)
		fscanf(InputDeck, "%s", aachar);

	while (fscanf(InputDeck, "%d %f %f", &nbin[isp], &temin, &temax) <3)
		fscanf(InputDeck, "%s", aachar);

	/* Assign input parameters to arguments */
	jj0[isp][0]= fabs(j0l);
	jj0[isp][1]= fabs(j0r);
	v0[isp][0] = fabs(v0l);
	v0[isp][1] = -fabs(v0r);
	vc[isp][0] = fabs(vcl);
	vc[isp][1] = -fabs(vcr);
	vt[isp][0] = fabs(vtl);
	vt[isp][1] = fabs(vtr);
	vtpt[isp]  = fabs(vtpt[isp]);
	emin[isp]  = temin;
	de[isp]	  = (temax-temin)/(nbin[isp]-1);
	dtheta[isp]= PI/2/(nbin[isp]-1);

	if((vt[isp][0] || v0[isp][0]) && (vt[isp][1] ==0.0 && v0[isp][1] ==0.0)) 
	{
		initnp[isp][0] = iinitn*FOURPI*(r1*r1*r1-r0*r0*r0)/(3.*nc2p) +.5;
		initnp[isp][1] = 0.0;
	}
	else if((vt[isp][1] || v0[isp][1]) && (vt[isp][0] ==0.0 && v0[isp][0] ==0.0)) 
	{
		initnp[isp][0] = 0.0;
		initnp[isp][1] = iinitn*FOURPI*(r1*r1*r1-r0*r0*r0)/(3.*nc2p) +.5;
	}
	else if(!vt[isp][0] && !vt[isp][1] && !v0[isp][0] && !v0[isp][1])
	{
		initnp[isp][0] = iinitn*FOURPI*(r1*r1*r1-r0*r0*r0)/(3.*nc2p) +.5;
		initnp[isp][1] = 0.0;
	}
	else
	{
		initnp[isp][0] = .5*iinitn*FOURPI*(r1*r1*r1-r0*r0*r0)/(3.*nc2p) +.5;
		initnp[isp][1] = .5*iinitn*FOURPI*(r1*r1*r1-r0*r0*r0)/(3.*nc2p) +.5;
	}
}	/* End SPECIES */

/***************************************************************/
void load(initnp)
int initnp[][2];
{
	float fdist();
	int i, ix, j, k, isp, ksign, ii;
	float xx, v1, v2, ddv, vv, sf, svf, ddf, sumf, ss, ds, ff, 
	      Vpp, vmag, theta, vel[1024];

	/************************************************/
	/* Load one species and one direction at a time */

	for (isp=0; isp<nsp; isp++)
	{
		np[isp]= 0;

		for (k=0; k<2; k++)
		{
			
			if (!initnp[isp][k]) continue;
			
			if (np[isp] + initnp[isp][k] > MAXLEN)
			{
				printf("LOAD: too many particles, species %2d",isp);
				exit(1);
			}
			if (vt[isp][k]==0. )	/* Loader for a COLD beam   */
			{
				Vpp= (r1*r1*r1-r0*r0*r0)/initnp[isp][k];
				for(i=0; i< initnp[isp][k]; i++)
				{
					ix = np[isp] +i;
					if(!i)	r[isp][ix]= pow((Vpp +r0*r0*r0) ,1./3);
					else	r[isp][ix]= pow((Vpp+r[isp][ix-1]*r[isp][ix-1]*r[isp][ix-1]) ,1./3);
					vr[isp][ix] = v0[isp][k];
					vth[isp][ix] =0.;
					vph[isp][ix] =0.;
				}
				if(dde)
				{
					for(i=0; i< initnp[isp][k]; i++)
					{
						ix = np[isp] +i;
						xx = TWOPI*(r[isp][ix] -r0)/(r1 -r0);
						r[isp][ix] += dde*sin(xx);
					}
				}
			}	/*end V0 != 0. */
			else /* Loader for THERMAL distribution */
			{	
				ksign = 1 - 2*k;

				/* Integrate the [distsfcn] AND v*[dist fcn] from v1
					to v2, using Simpson's rule, to normalize them */

				if (k==0)
				{
					v1 = max(vc[isp][0], v0[isp][0]-6.*vt[isp][0] );
					v2 = max(v0[isp][0]+6.*vt[isp][0], v1+3.*vt[isp][0]);
				}
				else
				{
					v1 = min(vc[isp][1], v0[isp][1]+6.*vt[isp][1] );
					v2 = min(v0[isp][1]-6.*vt[isp][1], v1-3.*vt[isp][1]) ;
				}

				ddv = (v2-v1)/2048.;
				vv = v1;
				sf = fdist(vv,k,isp)/2.;
				svf= vv*sf;

				for (i=0; i<1024; i++)
				{
					vv += ddv;
					ddf = fdist(vv, k, isp);
					sf += ddf + ddf;
					svf += vv*(ddf + ddf);
					vv += ddv;
					ddf = fdist(vv, k, isp);
					sf += ddf;
					svf += vv*ddf;
				}

				ddf = fdist(v2, k, isp);
				sf -= ddf;
				svf -= v2*ddf;
				sf = (sf + sf)*fabs(ddv)/3.;
				svf = (svf + svf)*fabs(ddv)/3.;
				/* Integrate again, varying the stepsize to aim for the
					values of the integral at which particle velocities
					are desired */

				sumf = 0.;
				ii = 0;
				vv = v1;
				ds = sf/1024.;	 
				ss = ds/2.;

L12:			ff = fdist(vv, k, isp);
				ddv = ksign*vt[isp][k]/128.;
				if (ff > 0.) ddv = ksign*min(fabs(ddv), .45*ds/ff);
				do
				{
					ddf = fdist(vv, k, isp) + 4.*fdist(vv + ddv, k, isp)
						   + fdist(vv+ddv+ddv, k, isp);	
					ddf *= fabs(ddv)/3.;		
					if (ddf >= ds) ddv *= 0.5;
				} while (ddf >= ds);

				vv += ddv+ddv;
				sumf += ddf;
				if (sumf < ss) goto L12;

				/* When a desired velocity has been passed, interpolate
					to find it, and start looking for the next velocity */
				vel[ii++] = vv - (sumf - ss)*ddv/ddf;
				ss += ds;
				if (ii < 1024) goto L12;

				/* Compute the number of particles to be loaded, and load
					the 1024 velocities uniformly in space, in bit-rev order */

				Vpp = (r1*r1*r1-r0*r0*r0)/initnp[isp][k];
				for(i=0; i< initnp[isp][k]; i++)
				{
					ix = np[isp]+ i;
					if(!i)	r[isp][ix]= pow((Vpp +r0*r0*r0) ,1./3);
					else	r[isp][ix]=pow((Vpp +r[isp][ix-1]*r[isp][ix-1]*r[isp][ix-1]) ,1./3);
					vr[isp][ix]	 = vel[index[(i+k)%1024]];
					vmag = vtpt[isp]*sqrt(-2*log(frand()));
					theta = TWOPI*frand();
					vth[isp][ix] = vmag*cos(theta);
					vph[isp][ix] = vmag*sin(theta);
				}
				if(dde)
				{
					for(i=0; i< initnp[isp][k]; i++)
					{
						ix = np[isp] +i;
						xx = TWOPI*(r[isp][ix] -r0)/(r1 -r0);
						r[isp][ix] += dde*sin(xx);
					}
				}
			}	/* end of THERMAL loader  "if( VT>0.." */
		  
			np[isp] +=initnp[isp][k];

		}	/* end  "for( k=0,1...)" loop */
	}  /* end for(isp=0 thru nsp-1..)" loop */
}  /* end LOAD  */

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

display_title()
{
	printf("\n\nPDS1-Spherical Bounded Electrostatic 1 Dimensional Code\n");
	printf("version 2.1\n");
	printf("(c) Copyright 1989-1992 Regents of the University of California\n");
	printf("Plasma Theory and Simulation Group\n");
	printf("University of California - Berkeley\n");
}						   

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