00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039 #include "G4NeutronField.hh"
00040 #include "G4PhysicalConstants.hh"
00041 #include "G4SystemOfUnits.hh"
00042 #include "G4VNuclearDensity.hh"
00043 #include "G4FermiMomentum.hh"
00044
00045 G4NeutronField::G4NeutronField(G4V3DNucleus * aNucleus) :
00046 G4VNuclearField(aNucleus), theDensity(theNucleus->GetNuclearDensity())
00047 {
00048 theA = theNucleus->GetMassNumber();
00049 theZ = theNucleus->GetCharge();
00050 theFermi.Init(theA, theZ);
00051 theR = 2.*theNucleus->GetOuterRadius();
00052 G4double aR=0;
00053 while(aR<theR)
00054 {
00055 G4ThreeVector aPosition(0,0,aR);
00056 G4double density = GetDensity(aPosition);
00057 G4double fermiMom = GetFermiMomentum(density);
00058 theFermiMomBuffer.push_back(fermiMom);
00059 aR+=0.3*fermi;
00060 }
00061 {
00062 G4ThreeVector aPosition(0,0,theR);
00063 G4double density = GetDensity(aPosition);
00064 G4double fermiMom = GetFermiMomentum(density);
00065 theFermiMomBuffer.push_back(fermiMom);
00066 }
00067 {
00068 G4ThreeVector aPosition(0,0,theR+0.001*fermi);
00069 theFermiMomBuffer.push_back(0);
00070 }
00071 {
00072 G4ThreeVector aPosition(0,0,1.*m);
00073 theFermiMomBuffer.push_back(0);
00074 }
00075 }
00076
00077 G4NeutronField::~G4NeutronField()
00078 { }
00079
00080 G4double G4NeutronField::GetField(const G4ThreeVector & aPosition)
00081 {
00082 G4double x = aPosition.mag();
00083 G4int index = static_cast<G4int>(x/(0.3*fermi) );
00084 if(index+2> static_cast<G4int>(theFermiMomBuffer.size())) return theFermiMomBuffer.back();
00085 G4double y1 = theFermiMomBuffer[index];
00086 G4double y2 = theFermiMomBuffer[index+1];
00087 G4double x1 = (0.3*fermi)*index;
00088 G4double x2 = (0.3*fermi)*(index+1);
00089 G4double fermiMom = y1 + (x-x1)*(y2-y1)/(x2-x1);
00090 return -1*(fermiMom*fermiMom)/(2*neutron_mass_c2);
00091 }
00092
00093 G4double G4NeutronField::GetBarrier()
00094 {
00095
00096
00097
00098
00099
00100
00101 return 0.;
00102 }
00103
00104
00105
00106
00107
00108