d.h

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

#if !defined(D_H_H__5A8BF52B_D92C_4512_8DDC_79FA24173FF3__INCLUDED_)
#define D_H__5A8BF52B_D92C_4512_8DDC_79FA24173FF3__INCLUDED_

#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000


#define	 Pi     3.14159265357
#include<iomanip>
#include"d.h"
#include<cmath>
#include<complex>
#include<fstream>
#include<ctime>
#include<cstdlib>
#include<vector>
using namespace std;
enum Nreactiontype{N2N,Elastic_Deflection,Nonelastic_Deflection,Capture,Total_Reaction};
enum Type{H,Be,B,C,N,O,F};
ifstream  mynuc;
bool accompany;
int k;


///////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////



class Nucleus;
class Neutron;
class MainSpace;
Nucleus * prime;

class Unitboard
{
public:
	int EnergyDistribution[200];
	std::vector<complex<double> > Points;
	int PlaceZ;
	double width;
	Nucleus * indicator;
	Unitboard * next;
	Unitboard * previous;

public:
	Unitboard(Nucleus *,int,double);
	void AddToBoard( Neutron *);

};

class Nucleus
{
public:
	double Mass;
	double Ep;
	Nucleus * previous;
	Nucleus * next;
	Type  nucleustype;

public:
	Nucleus(Type);
	double N_Capture_Cross_Section( Neutron );
	double N_Elastic_Deflection_Cross_Section( Neutron );
	double N_Nonelastic_Deflection_Cross_Section( Neutron );
	double N2N_Cross_Section( Neutron );
	double NReaction_Cross_Section( Neutron);
    void   Reaction(Nreactiontype,Neutron *);
	void   Deflection(double,Neutron *);
	bool   return_type(Nreactiontype);
	int    findline(double,double,int,int);
	double return_energy(Neutron);
};


class Neutron
{
public:
	double CosX;
	double CosY;
	double CosZ;
	double Energy;
	double PositionX;
	double PositionY;
	double PZ;      // note this is a length within width
	double Mass;
	Unitboard * boardnode;
	Neutron * next;
	Neutron * previous;
	MainSpace * mainnode;
	double free;



public:
	Neutron(double,double,double,double,double,double,double);
	void Move(double TimeInterval);
	void Collision(Nucleus);
};


class MainSpace
{
public:
	Unitboard * head;
	Unitboard * end;
	Unitboard * unitcurrent;
	Neutron * headneutron;
	Neutron * endneutron;
	Neutron * current;
	Neutron * Pointer;
	bool count;

public:
	MainSpace(Nucleus * ptr,double);
	void Addboard(Nucleus *);
	void Addneutron(double,double,double,double,double,double,double,Unitboard *);
	bool Deleteneutron();
	void Generateneutron(double);
	void JudgeCollision();
	void JudgeEdge();
	void Action(double);
	Unitboard * GetPlace(int);
};
/////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////
MainSpace ms(prime,10);
int CaptureNumber=0;
int EscapeNumber=0;
int LostNumber=0;
int Total=0;
int CollisionNum=0; int N2NCN=0; int ECN=0; int NCN=0;
double MainTime=0;
double Dist=2; double Thickness=0.3;
double Interwidth=0.1;

///////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////


Unitboard::Unitboard(Nucleus * nuc,int pz,double wid)
{
	width=wid;
	PlaceZ=pz;
	indicator=nuc;
	next=NULL;
	previous=NULL;
	for(int j=0;j<200;j++)
		EnergyDistribution[j]=0;

}


void Unitboard::AddToBoard(Neutron * nut)
{
	EnergyDistribution[(int)(nut->Energy/32767*200)]++;
	complex<double> a(nut->PositionX,nut->PositionY);
	Points.push_back(a);
}

Nucleus::Nucleus(Type p)
{
	Mass=10;
	Ep=2;
	next=NULL;
	previous=NULL;
	nucleustype=p;
}

int Nucleus::findline(double pt,double qt,int i, int j)
{
	int p;
	mynuc.seekg(i);  int a1; mynuc>>a1;	 mynuc.seekg(i-14);	 double b1; mynuc>>b1;
	mynuc.seekg(j);  int a2; mynuc>>a2;	 mynuc.seekg(j-14);  double b2; mynuc>>b2;
	if (pt>a2 || pt<a1 || (pt==a2 && qt>b2) || (pt==a1 && qt<b1)) return 0;
	else 
	{
		mynuc.seekg(j-27); int a3; mynuc>>a3;  mynuc.seekg(j-27-14);  double b3; mynuc>>b3;
		mynuc.seekg(i+27); int a4; mynuc>>a4;  mynuc.seekg(i+27-14);  double b4; mynuc>>b4;
		if ((pt>a3 || (pt==a3 && qt>b3)) && (pt<a2 ||(pt==a2 && qt<b2) )) 
			return j-27;
		if ((pt>a1 || (pt==a1 && qt>b1)) && (pt<a4 ||(pt==a4 && qt<b4) )) 
			return i;
		int q=i+27*(((j-i)/27+1)/2-1);
		mynuc.seekg(q); int a5; mynuc>>a5;	mynuc.seekg(q-14); double b5; mynuc>>b5;
		if (pt>a5 || (pt==a5 && qt>b5)) 
		{
			p=findline(pt,qt,q,j);
			return p;
		}			
		else if (pt<=a5) 
		{
			p=findline(pt,qt,i,q);
			return p;
		}
	}
	return 0;
}

double Nucleus::return_energy(Neutron n)
{
	mynuc.seekg(119);
	char a;  mynuc>>a;
	int totalnumber;  mynuc>>totalnumber; // totalnumber for line number
	double temp=log(n.Energy)/log(10);	
	int exponent=(temp>=0)?(int)temp:(int)temp-1;   
	double base=n.Energy/pow(10,exponent);
	int k=findline(exponent,base,241,241+(totalnumber-1)*27);
	if (k)
	{
	   mynuc.seekg(k+11); mynuc>>exponent;
	   mynuc.seekg(k+2);  mynuc>>base;
	   return (base*pow(10,exponent));
	}
	return 0;
}


double Nucleus::N_Capture_Cross_Section( Neutron n)	
{
	if(return_type(Capture))
	{	accompany=true;
	return return_energy(n);}
	else  return 0;

 
}

double Nucleus::N_Elastic_Deflection_Cross_Section( Neutron n)
{   
	if(return_type(Elastic_Deflection))
	{	accompany=true;
	return return_energy(n);}
	else  return 0;

}

double Nucleus::N_Nonelastic_Deflection_Cross_Section( Neutron n)
{   
	if(return_type(Nonelastic_Deflection))
	{	accompany=true;
	return return_energy(n);}
	else  return 0;

}		 

double Nucleus::N2N_Cross_Section( Neutron n)
{
	if(return_type(N2N))
	{	accompany=true;
	return return_energy(n);}
	else  return 0;

}

double Nucleus::NReaction_Cross_Section(Neutron n)
{
	if(return_type(Total_Reaction))
	{	accompany=true;
	return return_energy(n);}
	else  return 0;

}

void Nucleus::Deflection(double Ep,Neutron * nut)
{
	double v=sqrt(2*(nut->Energy)/nut->Mass);
	double v1=Mass*v/(Mass+1);  v1*=v1; v1-=2*Mass*Ep/(Mass+1); v1=sqrt(v1);
	double cosfi=2*(double)rand()/RAND_MAX-1;
	double sinfi=sqrt(1-cosfi*cosfi);
	double angle=(double)rand()/RAND_MAX*2*Pi;
	double SinZ=sqrt(1-(nut->CosZ)*(nut->CosZ));
	double vx=v1*cosfi*(nut->CosX);
	double vy=v1*cosfi*(nut->CosY);
    if (SinZ!=0)
	{
		   vx-=v1*sinfi*cos(angle)*(nut->CosZ)*(nut->CosX)/SinZ;
		   vx-=v1*sinfi*sin(angle)*(nut->CosY)/SinZ;
		   vx+=v*(nut->CosX)/(Mass+1);
	       vy-=v1*sinfi*cos(angle)*(nut->CosZ)*(nut->CosY)/SinZ;
		   vy+=v1*sinfi*sin(angle)*(nut->CosX)/SinZ;
		   vy+=v*(nut->CosY)/(Mass+1);
	}
	else
	{
		vx=vx-v1*sinfi*cos(angle);
		vy=vy+v1*sinfi*sin(angle);
	}
	double vz=v1*cosfi*(nut->CosZ)+v1*sinfi*cos(angle)*SinZ+v*(nut->CosZ)/(Mass+1);
	nut->Energy=0.5*(vx*vx+vy*vy+vz*vz);
	nut->CosX=vx/sqrt(2*(nut->Energy));
	nut->CosY=vy/sqrt(2*(nut->Energy));
	nut->CosZ=vz/sqrt(2*(nut->Energy));
}


void Nucleus::Reaction(Nreactiontype react,Neutron * nut)
{
	CollisionNum++;
	switch (react)  
	{
	case N2N:
		//
		nut->CosX=0.3;
		nut->CosY=0.4;
		nut->CosZ=sqrt(1-0.09-0.16);
		nut->mainnode->Addneutron(nut->PositionX,nut->PositionY,nut->PZ,0.5,0.6,sqrt(1-0.25-0.36),nut->Energy*0.13,nut->boardnode);
		nut->Energy*=0.87;
		N2NCN++;
		// add reaction information
		break;
	case Elastic_Deflection:
		//
//		nut->CosX=0.2;
//		nut->CosY=0.5;
//		nut->CosZ=sqrt(1-0.04-0.25);
		Deflection(0,nut);
		ECN++;
		// add reaction information
		break;
	case Nonelastic_Deflection:
//		nut->CosX=0.7;
//		nut->CosY=0.3;
//		nut->CosZ=sqrt(1-0.49-0.09);
//		nut->Energy*=0.56;
		Deflection(Ep,nut);
		NCN++;
		// add reaction information
		break;
	case Capture:
		nut->mainnode->Deleteneutron();
		CaptureNumber++;
		// add reaction information
	default:
		// do something
		break;
	}
}