fieldemi.cpp

pic 模拟程序！面向对象

				}
		}
	else											// Alongx1()
		{
			int kemit = k1;
			int jtest = (j1 + j2)/2;  //sample electric field at center of emitter.

			Scalar 	EmitWhichWay = -(j2-j1)/(fabs(j2-j1));
			//EmitWhichWay = -1 for emission downward
			//EmitWhichWay = +1 for emission upward

			if (j1 < j2) {kemit--;} // faces downward if j1 < j2

			Vector2	Eposition1(jtest, kemit + 0.5);
			Etest = fields->E(Eposition1).e2();

			if (emitFlag == FALSE)  //ONLY emit when E-field is larger than threshold
				{
					if (  sign*Etest*EmitWhichWay  < threshold)
						return particleList;
					emitFlag = TRUE;
				}


			Scalar	vBeam;
			Scalar	beta;
			Scalar	uBeam;

			Vector3	u;




			Vector2	deltaY;
			Scalar 	Ecell;
			Scalar	cellsize;

			Scalar  CurrentDen;
			Scalar  DeltaI;
			int 	numberEmitted = 0;
			Scalar  current;
			Vector2	position;
			Vector3	dyMKS;

			Scalar jlow=j1;
			Scalar jhigh=j2;
			if (j2<j1)
				{
					jlow = j2;
					jhigh = j1;
				}

			for (int jcell = (int)jlow; jcell < jhigh; jcell++)
				{
					/*      metric elements for r only
									r1 = grid->getMKS(jemit+0.5, kcell).e2();
									r1s = r1*r1;
									r2 = grid->getMKS(jemit+0.5, kcell+1).e2();
									r2s = r2*r2;
									*/

					cellsize = 0.5*(grid->dl2(jcell, kemit)
													+ grid->dl2(jcell+1, kemit));

					//find normal electric field at cell center
					Ecell = fields->E(Vector2(jcell+0.5,kemit+0.5)).e2();
					//	cout << "Ecell is " << Ecell << endl;


					// calculate current density (A/m^2) needed at cell center:
					CurrentDen = EmitWhichWay*sign*
						sqrt(fabs(ELECTRON_CHARGE*Ecell*Ecell*Ecell)
								 /(6*0.5*cellsize*iEPS0*iEPS0*ELECTRON_MASS));

					//	vBeam NOT relativistic
					vBeam = sqrt(fabs(ELECTRON_CHARGE*Ecell*cellsize/ELECTRON_MASS));
					beta = vBeam*iSPEED_OF_LIGHT;
					uBeam = vBeam/sqrt(1 - beta*beta);
					uBeam *= EmitWhichWay;			// ensure proper sign
					u = Vector3(0, uBeam, 0);

					// get existing currents and average to cell center;
					current = 0.5*(fields->getIMesh(jcell,kemit).e2()
												 + fields->getIMesh(jcell+1,kemit).e2());

					//calculate additional total current and current density
					// at cell center which is needed
					// check this!
					DeltaI = CurrentDen*0.5*(grid->dS2(jcell,kemit)
																	 + grid->dS2(jcell,kemit+1)) - current;

					// find number of particles to be emitted and add nearly one to round-up
					// the value...
					numberEmitted = 1 + (int)fabs(DeltaI*dt/get_q());

					//	For electrons, add current when DeltaI < 0; for ions add
					//	current when DeltaI > 0, where DeltaI is the change in current
					//	required to make E_normal = 0 at the emitter.
					if (sign*DeltaI*EmitWhichWay > 0)
						{
							for (int ipart=0; ipart < numberEmitted; ipart++)
								{
									// below is position ON BOUNDARY where particle should be emitted
									position = Vector2(jcell+frand(), kemit + 0.5*(1.-EmitWhichWay) );
									Vector2 HalfCellOffBoundary = Vector2(position.e1(), kemit+0.5);

									// Now subtract to find emission direction vector
									Vector2 temp = fields->getMKS(HalfCellOffBoundary -position);
									dyMKS = Vector3(0, temp.e2(), 0);
									// below should check boundary intersection
									fields->translateAccumulate(position, dyMKS, get_q()/dt);
									particleList.push(new Particle(position, u, species, np2c));
								}
						}
				}
		}
	return particleList;
#endif
}

//--------------------------------------------------------------
// FieldEmitter2
//
// This bc injects charge sufficient to achieve E =0 normal to surface
// by Gauss's Law

FieldEmitter2::FieldEmitter2(oopicList <LineSegment> *segments,
			     Species* s, Scalar np2c, Scalar t, 
			     MaxwellianFlux *max) 
  : Emitter(segments, s, np2c)
{
  maxwellian = max;
  threshold = fabs(t);
//	t = 1E7 for most materials
  emitFlag = FALSE;
  emitDirection = normal*((species->get_q() > 0) ? 1 : -1);
  delta = normal*1e-4;
  grid = NULL;
  init = TRUE;
}

//---------------------------------------------------------------
// assume orthogonal boundary

ParticleList& FieldEmitter2::emit(Scalar t, Scalar dt, Species* species_to_push)
{
  Scalar N1, N2, R;

  if (species==species_to_push){
    if (init) initialize();
 
//    Scalar Dp, Dm;
    Boundary* bPtr;
    Scalar del_t;

    if (alongx1()){   // horizontal emitter
      Vector2 x(0, k1 + delta);
		Vector3 dxMKS(0,0,0);
      int kindex = (normal==1)?k1:k1-1;
      for (int j=j1; j<j2; j++){
		  N1 = normal*fields->get_epsi(j,kindex)*fields->getIntEdl()[j][kindex].e2()/grid->dl2(j, kindex)*grid->dSPrime()[j][kindex].e2()-rho[j][k1]*grid->get_halfCellVolumes()[j][k1];
		  N2 = normal*fields->get_epsi(j+1,kindex)*fields->getIntEdl()[j+1][kindex].e2()/grid->dl2(j+1, kindex)*grid->dSPrime()[j+1][kindex].e2()-rho[j+1][k1]*grid->get_halfCellVolumes()[j+1][k1];
		  N1 /= species->get_q()*np2c;
		  N2 /= species->get_q()*np2c;
		  Scalar N=0.5*(N1+N2);
		  // Scalar Nthres = fabs(fields->get_epsi(j,kindex)*threshold*
		  // grid->dS1(j,kindex)
		  // /(species->get_q()*np2c)); 
	
		  // N -= Nthres;
		  for (int i=0; i < (int)N; i++){
			 del_t = dt*(i+0.5)/((int) N);
			 R=frand();
			 // x.set_e1(j - (N1-sqrt(-(R-1)*sqr(N1)+R*sqr(N2)))/(N2-N1));
			 if (fabs(N2-N1)>1E-6)
				x.set_e1(j - (N1-sqrt(N1*N1*(1-R) + R*N2*N2))/(N2 - N1));
			 else 
				x.set_e1(j + R);
			 if (x.e1() < j || x.e1() > j+1) printf("FieldEmitter2: x = %G\n", x.e1());
			 Vector3 u=maxwellian->get_U(normal);
			 Particle* p = new Particle(x, u, species, np2c);
			 bPtr = initialPush(del_t, dt, *p);
			 if (!bPtr) particleList.add(p);
			 else bPtr->collect(*p, dxMKS);
		  }
      }
    }
    else { // vertical
      Vector2 x(j1 + delta, 0);
      int jindex = (normal==1)?j1:j1-1;
      for (int k=k1; k<k2; k++){
		  N1 = normal*fields->get_epsi(jindex,k)*fields->getIntEdl()[jindex][k].e1()/grid->dl1(jindex,k)* grid->dSPrime()[jindex][k].e1()-rho[j1][k]*grid->get_halfCellVolumes()[j1][k];
		  N2 =normal*fields->get_epsi(jindex,k+1)*fields->getIntEdl()[jindex][k+1].e1()/grid->dl1(jindex, k+1)*grid->dSPrime()[jindex][k+1].e1()-rho[j1][k+1]*grid->get_halfCellVolumes()[j1][k+1];
		  N1 /= species->get_q()*np2c;
		  N2 /= species->get_q()*np2c;
		  Scalar N=0.5*(N1+N2);
		  //Scalar Nthres = fabs(fields->get_epsi(jindex,k)*threshold*
		  // grid->dS1(jindex,k)
		  // /(species->get_q()*np2c)); 
	
		  // N-=Nthres;
		  //	cout << " number of particles " << N << " at j,k "<< jindex << " "<< k << endl;  
		  //	cout << " N1 " << N1 << " N2 " << N2 << endl;
		  for (int i=0; i<(int)N; i++){
			 del_t = dt*(i+0.5)/((int) N);
			 R=frand();
			 if (rweight) R=sqrt(k*k + R*(2*k + 1))-k;
			 x.set_e2(k - (N1-sqrt(-(R-1)*sqr(N1)+R*sqr(N2)))/(N2-N1));
			 Vector3 u=maxwellian->get_U(normal);
			 Particle* p = new Particle(x, u, species, np2c);
			 bPtr = initialPush(del_t, dt, *p);
			 if (!bPtr) particleList.add(p);
		  }
      }
    }
  }
  return particleList;
}

void FieldEmitter2::initialize()
{
  Emitter::initialize();
  grid = fields->get_grid();
  rho = fields->getrho();
}

#ifdef old_code
  if (species==species_to_push){
	 if (!grid) grid = fields->get_grid();
	 Vector3** ENode = fields->getENode();
	 if (!emitFlag){ // First, determine if it should begin emitting
		if (alongx1()){  // horizontal
		  for (int j=j1; j<j2; j++) 
			 if (ENode[j][k1].e2() > emitDirection*threshold)
				emitFlag = TRUE;
		}
		else {  // vertical
		  for (int k=k1; k<k2; k++)
			 if (ENode[j1][k].e1() > emitDirection*threshold)
				emitFlag = TRUE;
		}
	 }
	 
	 if (emitFlag){
		Scalar Dp, Dm;
		Boundary* bPtr;
		Scalar del_t;
		if (alongx1()){   // horizontal emitter
		  Vector2 x(0, k1 + delta);
		  for (int j=j1; j<j2; j++){
			 Scalar N = 0.5*fields->get_epsi(j,k1)*(ENode[j][k1].e2() + ENode[j+1][k1].e2())
				*grid->dS2(j,k1);
			 N /= species->get_q()*np2c;
			 for (int i=0; i < (int)N; i++){
				del_t = dt*(i+0.5)/((int) N);
				x.set_e1(j+frand());
				Vector3 u; // zero by default
				Particle* p = new Particle(x, u, species, np2c);
				bPtr = initialPush(del_t, dt, *p);
				if (!bPtr) particleList.add(p);
			 }
		  }
		}
		else { // vertical
		  Vector2 x(j1 + delta, 0);
		  for (int k=k1; k<k2; k++){
			 Scalar N = 0.5*fields->get_epsi(j1,k)*(ENode[j1][k].e1() + ENode[j1][k+1].e1())
				*grid->dS1(j1,k);
			 N /= species->get_q()*np2c;
			 for (int i=0; i<(int)N; i++){
				del_t = dt*(i+0.5)/((int) N);
				if (rweight) x.set_e2(sqrt(k*k + frand()*(2*k + 1)));
				else x.set_e2(k + frand());
				Vector3 u;
				Particle* p = new Particle(x, u, species, np2c);
				bPtr = initialPush(del_t, dt, *p);
				if (!bPtr) particleList.add(p);
			 }
		  }
		}
	 }
  }
  return particleList;
}
#endif