d.h

使用蒙特卡洛方法研究中子屏蔽以活求碰撞情况

bool Nucleus::return_type(Nreactiontype react)
{
	if(accompany) {mynuc.close(); accompany=false;}
	switch (nucleustype)
	{										//enum Nreactiontype{N2N,Elastic_Deflection,Nonelastic_Deflection,Capture};
		                                    //   enum Type{H,Be,B,C,N,O,F};
	case H:							
		switch (react)
		{
		case N2N:
			    return false;
			break;
		case Elastic_Deflection:
			mynuc.open("H-Elastic.txt",ios::in);
			if (mynuc.bad())
				return false;
			break;
		case Nonelastic_Deflection:
			    return false;
			break;
		case Capture:
			mynuc.open("H-Capture.txt",ios::in);
			if (mynuc.bad())	
			    return false;
			break;
		case Total_Reaction:
			mynuc.open("H-Total.txt",ios::in);
			if (mynuc.bad())
			{		    return false; }
            break;
		default:
			return false;
			break;
		}
		break;
	case O:
		switch (react)
		{
		case N2N:
			mynuc.open("O-N2N.txt",ios::in);
			if (mynuc.bad())	
			    return false;
			break;
		case Elastic_Deflection:
			mynuc.open("O-Elastic.txt",ios::in);
			if (mynuc.bad())	
			    return false;
			break;
		case Nonelastic_Deflection:
			mynuc.open("O-Nonelastic.txt",ios::in);
			if (mynuc.bad())	
			    return false;
			break;
		case Capture:
			mynuc.open("O-Capture.txt",ios::in);
			if (mynuc.bad())	
			    return false;
			break;
		case Total_Reaction:
			mynuc.open("O-Total.txt",ios::in);
			if (mynuc.bad())	
			    return false;
			break;
		default:
			return false;
			break;
		}
		break;
	case N:
		switch (react)
		{
		case N2N:
			mynuc.open("N-N2N.txt",ios::in);
			if (mynuc.bad())	
			    return false;
			break;
		case Elastic_Deflection:
			mynuc.open("N-Elastic.txt",ios::in);
			if (mynuc.bad())	
			    return false;
			break;
		case Nonelastic_Deflection:
			mynuc.open("N-Nonelastic.txt",ios::in);
			if (mynuc.bad())	
			    return false;
			break;
		case Capture:
			mynuc.open("N-Capture.txt",ios::in);
			if (mynuc.bad())	
			    return false;
			break;
		case Total_Reaction:
			mynuc.open("N-Total.txt",ios::in);
			if (mynuc.bad())	
			    return false;
			break;
		default:
			return false;
			break;
		}
		break;
	default:
		return false;
		break;
	}
    return true;
}





Neutron::Neutron(double px,double py,double pz,double cx,double cy,double cz,double E)
{
	PositionX=px;
	PositionY=py;
	PZ=pz;
	CosX=cx;
	CosY=cy;
	CosZ=cz;
	Energy=E;
	Mass=1;
	boardnode=NULL;
	previous=NULL;
	next=NULL;
	free=0;

}

void Neutron::Move(double TimeInterval)
{
	double Velocity=sqrt(2*Energy/Mass);
	PositionX+=Velocity*CosX*TimeInterval;
	PositionY+=Velocity*CosY*TimeInterval;
    PZ+=Velocity*CosZ*TimeInterval;
	free+=Velocity*TimeInterval;
}





void Neutron::Collision( Nucleus nuc)
{		
	double ccs=nuc.N_Capture_Cross_Section( *this);
	double edcs=nuc.N_Elastic_Deflection_Cross_Section(*this);
	double ndcs=nuc.N_Nonelastic_Deflection_Cross_Section( *this );
	double n2ncs=nuc.N2N_Cross_Section( *this );
	double total=ccs+edcs+ndcs+n2ncs;
	double k=(double)rand()/RAND_MAX;
	if (k<=ccs/total)            nuc.Reaction(Capture,this);                               
	else if (k<=(ccs+edcs)/total)       nuc.Reaction(Elastic_Deflection,this);
	else if (k<=(ccs+edcs+ndcs)/total)  nuc.Reaction(Nonelastic_Deflection,this);
	else if (k<=1)       nuc.Reaction(N2N, this);
}


MainSpace::MainSpace(Nucleus * ptr,double Size)
{
	head=end=unitcurrent=NULL;
	current=headneutron=endneutron=NULL;
	Pointer=NULL;
	for(int i=1;i<=Size/Interwidth;i++)
	    Addboard(ptr);
	headneutron=NULL;
	endneutron=NULL;
	count=false;
}

void MainSpace::Addboard(Nucleus * ptr)
{
	if (head==NULL) 
	{
		head=new Unitboard(NULL,0,Dist);
		end=new Unitboard(ptr,1,Interwidth);
		head->next=end;
		end->previous=head;
	}
	else
	{
		unitcurrent=new Unitboard(ptr,end->PlaceZ+1,Interwidth);
		end->next=unitcurrent;
		unitcurrent->previous=end;
		end=end->next;
	}
}


void MainSpace::Addneutron(double px,double py,double pz,double cx,double cy,double cz,double E,Unitboard * ptr)
{
	current=new Neutron(px,py,pz,cx,cy,cz,E);
	current->boardnode=ptr;
	current->mainnode=this;
//	Neutron a(1,1,1,1,-1,1,20);
//	Neutron b(1.2,1,1,1,-1,1,20);
	if (headneutron==NULL) 
	{
		headneutron=endneutron=current;
	}
	else
	{
		current->previous=endneutron;
		endneutron->next=current;
		endneutron=endneutron->next;
	}
	Total++;
}

bool MainSpace::Deleteneutron()
{
	count=true;
	if (Pointer==NULL) return 1;
	else if (Pointer==headneutron && Pointer==endneutron) {  headneutron=endneutron=NULL; delete Pointer; Pointer=NULL; return 1;}
	else if (Pointer==headneutron) { headneutron=Pointer->next;  headneutron->previous=NULL; delete Pointer; Pointer->next=headneutron; return 1;}
	else if (Pointer==endneutron) { endneutron=Pointer->previous; endneutron->next=NULL; delete Pointer; Pointer=NULL; return 1;}
	else 
	{
		Pointer->previous->next=Pointer->next;
		Pointer=Pointer->previous;
		delete Pointer->next->previous;
		Pointer->next->previous=Pointer;
		return 1;
	}
 }

 void MainSpace::Generateneutron(double emissivity)
 { 
	 //  add code here
	 //
	 if ((double)rand()/RAND_MAX<emissivity)
	 {
		 Addneutron((double)rand()/RAND_MAX,(double)rand()/RAND_MAX,0,0,0,1,30,head);
	 }
	 //

 }                                        // the emission of neutron

 void MainSpace::JudgeCollision()
 {
	 Nucleus * ptr=Pointer->boardnode->indicator;
	 double Total_Cross_Section=0;;
	 while(ptr!=NULL)
	 {
		 Total_Cross_Section+=ptr->NReaction_Cross_Section(*Pointer);
		 ptr=ptr->next;
	 }
	 if(Total_Cross_Section)
	 {
	    double mean_free_pass=1/Total_Cross_Section;
	    if (rand()/RAND_MAX<1-exp((-Pointer->free)/mean_free_pass))
		{
		    double k=(double)rand()/RAND_MAX;
		    double possibility=0;
		    Nucleus * p=Pointer->boardnode->indicator;
		    do
			{
			   possibility+=p->NReaction_Cross_Section(*Pointer);
			   if (k<=possibility/Total_Cross_Section) 
			   {
				   Pointer->free=0;
				   Pointer->Collision(*(p));
				   break;
			   }
			   else p=p->next;
			} while (p!=NULL);
		}
	 }
 }


 void MainSpace::JudgeEdge()
 {
	 if (Pointer->PZ>0 && Pointer->PZ<Pointer->boardnode->width) 
	 {}
	 else if (Pointer->PZ<0 && Pointer->boardnode==head->next)
	 {
		 Deleteneutron();
		 LostNumber++;
	 }
	 else if (Pointer->PZ>Pointer->boardnode->width && Pointer->boardnode==end)
	 {
		 Deleteneutron();
		 EscapeNumber++;
	 }
	 else if (Pointer->PZ>Pointer->boardnode->width && Pointer->boardnode!=end)
	 {
		 Pointer->boardnode->AddToBoard(Pointer);
		 Pointer->PZ-=Pointer->boardnode->width;
		 Pointer->boardnode=Pointer->boardnode->next;
	 }
	 else if (Pointer->PZ<Pointer->boardnode->width && Pointer->boardnode!=head)
	 {
		 Pointer->boardnode->AddToBoard(Pointer);
		 Pointer->boardnode=Pointer->boardnode->previous;
		 Pointer->PZ+=Pointer->boardnode->width;
	 }
 }



 void MainSpace::Action(double TimeInterval)
 {
	 Generateneutron(0.3);
	 Pointer=headneutron;
	 while (Pointer!=NULL)
	 {	 
		 JudgeCollision();
		 if (!count)
		 {
			 Pointer->Move(TimeInterval);
		     JudgeEdge();
		 }
		 if (Pointer) Pointer=Pointer->next;
		 count=false;
	 }
	 MainTime+=TimeInterval;
 }


 Unitboard * MainSpace::GetPlace(int i)
 {
	 int number=0;  Unitboard * temp=head;
	 do
	 {
		 if (number==i)	 { return temp;	  break;}
		 else { temp=temp->next;  number++;}
	 }while (temp!=NULL);
	 return NULL;
 }

#endif
 ///////////////////////////////////////////////////////////////////////////////////////////
 ///////////////////////////////////////////////////////////////////////////////////////////
 ///////////////////////////////////////////////////////////////////////////////////////////