padjuss.c

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

			if(random <prob) 
			{
				/* determine the type of collision and calculate new velocities, etc */
				random /= prob;
				sig1 /= stot +DBL_MIN;
				
				/**********************************/
				/* if the collision is scattering */
            if (random <= sig1)
				{
					float up1, up2, up3, mag;
					float r11, r12, r13, r21, r22, r23, r31, r32, r33;
	
					coschi = sqrt(1. - frand());
					sinchi = sqrt(1 - coschi*coschi);

					phi1  = TWOPI*frand();
					cosphi= cos(phi1);
					sinphi= sin(phi1);
	
					r13 = vr[icolsp][j]/vel;
					r23 =vth[icolsp][j]/vel;
					r33 = vph[icolsp][j]/vel;
	
					if(r33 == 1.0) { up1= 0;  up2= 1;  up3= 0; }
					else           { up1= 0;  up2= 0;  up3= 1; }
	
					r12 = r23*up3 -r33*up2;
					r22 = r33*up1 -r13*up3;
					r32 = r13*up2 -r23*up1;
					mag = sqrt(r12*r12 + r22*r22 + r32*r32);
	
					r12/= mag;
					r22/= mag;
					r32/= mag;
	
					r11 = r22*r33 -r32*r23;
					r21 = r32*r13 -r12*r33;
					r31 = r12*r23 -r22*r13;
	
					vr[icolsp][j] = vel*coschi*(r11*sinchi*cosphi +r12*sinchi*sinphi +r13*coschi);
					vth[icolsp][j]= vel*coschi*(r21*sinchi*cosphi +r22*sinchi*sinphi +r23*coschi); 
					vph[icolsp][j] = vel*coschi*(r31*sinchi*cosphi +r32*sinchi*sinphi +r33*coschi);
				}

				/************************************/
				/* if the collision is charge exch. */
				else
				{
				   engy = -gtemp*log(frand()+DBL_MIN);
               vel = sqrt(fabs(engy)/Escale[icolsp]);

				   vr[icolsp][j] = 1 -2*frand();
				   vth[icolsp][j]= 1 -2*frand();
				   vph[icolsp][j]= 1 -2*frand();
         	   vel /= sqrt(vr[icolsp][j]*vr[icolsp][j]
							+vth[icolsp][j]*vth[icolsp][j]
							+vph[icolsp][j]*vph[icolsp][j]);
				
					vr[icolsp][j]  *= vel;
					vth[icolsp][j] *= vel;
					vph[icolsp][j] *= vel;

					if ( r[icolsp][j] >= r0 && r[icolsp][j] <= r1)
					   {
						  s= r[icolsp][j] -r0;

						  if (r[icolsp][j] < r0+.5)
					     del = ((r[icolsp][j]+.5)*(r[icolsp][j]+.5)*(r[icolsp][j]+.5)-rg3[0])/
						        ((r[icolsp][j]+.5)*(r[icolsp][j]+.5)*(r[icolsp][j]+.5)-r0*r0*r0);
					     else
						  {						   
				           del =((r[icolsp][j]+.5)*(r[icolsp][j]+.5)*(r[icolsp][j]+.5)-rg3[s])
						          /(3.*r[icolsp][j]*r[icolsp][j]+0.25);
						  }

						  if (r[icolsp][j] > r1-.5)
					        del =(r1*r1*r1-r[icolsp][j]*r[icolsp][j]*r[icolsp][j])/
						  		  (r1*r1*r1-(r[icolsp][j]-.5)*(r[icolsp][j]-.5)*(r[icolsp][j]-.5));

						}
					chrgxrate[s] += 1 - del;
					chrgxrate[s+1] += del;

				}
			}
	  	}
	}

	if(secondary)
	{
		i = np[secsp];
		np[secsp] += secountl + secountr;
		
		if(np[secsp] > MAXLEN)
			ExitMsg("ADJUST(Secondaries): too many particles, MUST EXIT!");
		
	   for(j=i; j< i+secountl; j++)
		{
			r[secsp][j]= r0 +frand();
			vr[secsp][j]= vnext(0,secsp);
			vmag= vtpt[secsp]*sqrt(2*fabs(log(frand()+DBL_MIN)));
			theta= TWOPI*frand();
			vth[secsp][j] = vmag*sin(theta);
			vph[secsp][j] = vmag*cos(theta);
		}
		i += secountl;
		for(j=i; j< i+secountr; j++)
		{
			r[secsp][j]= r1 -frand();
			vr[secsp][j]= vnext(1,secsp);
			vmag= vtpt[secsp]*sqrt(2*fabs(log(frand()+DBL_MIN)));
			theta= TWOPI*frand();
			vth[secsp][j] = vmag*sin(theta);
			vph[secsp][j] = vmag*cos(theta);
		}
	}/* end of if(secondary) */

	for (isp=0; isp<nsp; isp++) jwall[isp] /= dt;

}/* end ADJUST */

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

float vnext(k, isp)
/* Function to give bit-reversed velocities to the injection
	scheme of ADJUST, one at a time */
int k, isp;
{
	char locbuf[80];
	extern int index[];
	extern float v0[][2], vt[][2], vc[][2];
	static int iv[NSMAX][2], qinit=1;
	static float vel[1024][NSMAX][2], svf[NSMAX][2];
	int iiv, ksign, ii, i;
	float v1, v2, ddv, vv, vf, df, sumvf, ss, ds, fdist(), revers();

	/* At t=0., set up iv arrays */
	if (qinit)
	{
		for (i=0; i<NSMAX; i++)	iv[i][0] = iv[i][1] = 0;
		qinit = 0;
	}

	/* For a COLD beam, injection is simple: */
	if (vt[isp][k]==0.) return (v0[isp][k]);

	if (iv[isp][k]==0)
	{
		/* Integrate velocity-distribution function
			using Simpsons's rule (to normalize it) */
		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;
		svf[isp][k] = vv*fdist(vv,k,isp)/2.;

		for (i=0; i<1023; i++)
		{
			vv += ddv;
			df = fdist(vv,k,isp);
			svf[isp][k] += vv*(df+df);
			vv += ddv;
			df = fdist(vv,k,isp);
			svf[isp][k] += vv*df;
		}	/* end "for(i...)" loop */

		df = fdist(v2,k,isp)/2.;
		svf[isp][k] = (svf[isp][k] + v2*df)*2.*ddv/3.;
	}	/* end "if(iv[...)" loop */

	iiv = iv[isp][k]%1024;
	if (iiv==0)
	{
		ksign =  1 - 2*k;
		/* ksign = +1(k=0)  or -1(k=1) */
		sumvf = 0.;
		ii = 0;
		if (k==0) vv = max(vc[isp][0], v0[isp][0]-6.*vt[isp][0]);
		else vv = min(vc[isp][1], v0[isp][1]+6.*vt[isp][1]);
		ds = svf[isp][k]/1024.;
		iiv = iv[isp][k]/1024.;
		ss = revers(iiv)*ds;

L11:	ddv = ksign*vt[isp][k]/128.;
		vf = vv*fdist(vv,k,isp);
		
		if (vf!=0.) ddv=ksign*min(fabs(ddv),.45*ds/fabs(vf));
		/* Vary the integration step to aim for the next
			desired value for the integral */

		do
		{
			df = ddv*(vf + 4.*(vv+ddv)*fdist(vv+ddv,k,isp) +
				(vv+ddv+ddv)*fdist(vv+ddv+ddv,k,isp))/3.;
			ddv *= 0.5;
		} while (df >= ds);
		/* If integral increases too much, go back w. dv/2 step */
		ddv *= 2.;

		vv += 2.*ddv;
		sumvf += df;
		if (sumvf < ss) goto L11;

		vel[ii][isp][k] = vv - (sumvf-ss)*ddv/df;
		/* Interpolate to find the desired velocity */

		ii += 1;
		ss += ds;
		if (ii < 1024) goto L11;

		if (iv[isp][k]==0) iv[isp][k] = 1;
		/* Skip the very first velocity */
	}  /* end "if(iiv..)"  loop */

	ii = iv[isp][k]%1024;
	iiv = index[ii];

	iv[isp][k] += 1;
	if (iv[isp][k] < 0) iv[isp][k] = 1;
	/* Keep track of the total number of particles injected */

	return (vel[iiv][isp][k]);
}	/* end of VNEXT */

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

newvel(energy, vel, vr, vth, vph, e_flag)
int e_flag;
float energy, vel;
float far *vr, far *vth, far *vph;
{
  float phi1, cosphi, sinphi, coschi, sinchi, up1, up2, up3;
  float mag, r11, r12, r13, r21, r22, r23, r31, r32, r33;

  if(energy < 1e-30)  coschi = 1;
  else  coschi = (energy +2 -2*pow(energy +1,frand()))/energy;
  sinchi= sqrt(1. - coschi*coschi);
  
  phi1  = TWOPI*frand();
  cosphi= cos(phi1);
  sinphi= sin(phi1);
  
  if(e_flag)  vel *= sqrt(1 - 2*m[ecolsp]*(1-coschi)/m[ionsp]);
  
  r13 = *vr;
  r23 = *vth;
  r33 = *vph;
  
  if(r33 == 1.0) { up1= 0;  up2= 1;  up3= 0; }
  else           { up1= 0;  up2= 0;  up3= 1; }
  
  r12 = r23*up3 -r33*up2;
  r22 = r33*up1 -r13*up3;
  r32 = r13*up2 -r23*up1;
  mag = sqrt(r12*r12 + r22*r22 + r32*r32);
  
  r12/= mag;
  r22/= mag;
  r32/= mag;
  
  r11 = r22*r33 -r32*r23;
  r21 = r32*r13 -r12*r33;
  r31 = r12*r23 -r22*r13;
  
  *vr= vel*(r11*sinchi*cosphi +r12*sinchi*sinphi +r13*coschi);
  *vth= vel*(r21*sinchi*cosphi +r22*sinchi*sinphi +r23*coschi); 
  *vph= vel*(r31*sinchi*cosphi +r32*sinchi*sinphi +r33*coschi);
}

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