NIC - NEUTRON INELASTIC COOLING TO UCN ENERGIES

a tool for studying UCN production.
Roger Hill.  August, 1996


Needs["Graphics`Graphics`"]

-I. In thermal equilibrium with a moderator:

Clear["Global`*"];
T = 25;     (* Moderator temperature in K *)
WeightT = 1.; (*Temperature weighting factor for creating composite spectra *) 
Vmax = 8.0;   (*Definition of UCN = Vmax for UCN in m/s *)
rho0 = 1.835 10^9; (* incident neutron density, neutrons/cm^3  *)

kbzm = 1.3807 10^-23; (* Boltzman Constant in MKS *)
mass = 1.6749 10^-27;  (* neutron mass in kg *)
lVConstant = 3957; (* V Lamda = 3957 for V in m/s and Lamda in Angstroms for neutrons*)
ThetaDebye = 110;
alpha = Sqrt[2 kbzm T/mass]; (* most probable (peak) velocity of the Maxwell Boltzman Distribution *)
F = WeightT*rho0*2*alpha/Sqrt[N[Pi]]*100; (* neutron flux in neutrons/cm^2/s *)
Vmax = 3*alpha; (*TEMPORARY. Use this to calculate the full Maxwellian spectrum rather than just the UCN spectrum
                 by using Vmax = 8 m/s.*)
Emax = 0.5 mass Vmax^2/kbzm;  (*Emax in K*)
lMin = lVConstant/Vmax;
SetOptions[Plot, DefaultFont->{"Courier",12}];
labelV1 := NumberForm[StringForm["Equilibrium N(v) for a Moderator at `` K",T],5];
label2 := NumberForm[StringForm["Incident neutron flux = `` neutron/cm^2/s", F],5];
labelV2 := NumberForm[StringForm[ " `` neutrons/cm^3 with V < `` m/s", DensityV, Vmax],5];
labelE1 := NumberForm[StringForm["Equilibrium N(E) for a Moderator at `` K",T],5];
labelE2 := NumberForm[StringForm["Equilibrium Flux(E) for a Moderator at `` K",T],5];
labelE3 := NumberForm[StringForm[" `` neutrons/cm^3 with E < `` K", DensityE, Emax],5];
labelE4 := NumberForm[StringForm["`` neutron/cm^2/s with E < `` K",  FluxE, Emax],5];
labelE5 := NumberForm[StringForm["Vbar = `` m/s",Vbar],5];
labelL1 := NumberForm[StringForm["Equilibrium N(lamda) for a Moderator at `` K",T],5];
labelL2 := NumberForm[StringForm[ " `` neutrons/cm^3 with lamda > `` Angstroms", DensityL, lMin],5];
blank := "  ";

Nv[v_] := 2*F 10^-2 /(alpha^2)*(v/alpha)^2*Exp[-(v/alpha)^2];
(* units are #/cm^3 when v is in m/s and F is #/cm^2/s (this is why there is a 10^-2 in the formula) *)
Nvt[vupper_] := NIntegrate[Nv[v],{v,0,vupper}];
DensityV := Nvt[Vmax];
SaveTheNv[fileName_,numberOfPoints_] := Do[Write[fileName,CForm[SequenceForm[j*Vmax/numberOfPoints,blank,Nv[j*Vmax/numberOfPoints]]]],{j,1,numberOfPoints}];  
PlotNv := Plot[Nv[x], {x,0,Vmax}, PlotLabel->ColumnForm[{labelV1,label2,labelV2,blank},Center],AxesLabel->{"Velocity (m/s)", "Density (#/cm^3/(m/s))"}];

Nl[la_] := 2*F 10^-2  lVConstant^3 1/(la alpha)^4 Exp[-(lVConstant/la/alpha)^2];
(* units are #/cm^3 when v is in m/s and F is #/cm^2/s (this is why there is a 10^-2 in the formula) *)
BrightL[l_] := Nl[l] lVConstant/l 100/8/N[Pi]; (* Brightness in n/cm^2/s/Å/sr at l Å*)
Nlt[lower_] := NIntegrate[Nl[l],{l,lower,Infinity}];
DensityL := Nlt[lMin];
PlotNl := Plot[Nl[x], {x,lMin,10*lMin}, PlotLabel->ColumnForm[{labelL1,label2,labelL2,blank},Center],AxesLabel->{"Wave Length (A) ", "Density (#/cm^3/A)"}];

Ne[e_] :=  N[Sqrt[mass/2/kbzm] F/T^2 10^-2  Sqrt[e]*Exp[-e/T]];
(* units are #/cm^3/K when v is in m/s and F is #/cm^2/s (this is why there is a 10^-2 in the formula) *)
Net[eupper_] := NIntegrate[Ne[x],{x,0,eupper}];
DensityE := Net[Emax];
PlotNE := Plot[Ne[x], {x,0,Emax}, PlotLabel->ColumnForm[{labelE1,label2,labelE3,blank},Center],AxesLabel->{"Energy (K)", "Density (#/cm^3/K)"}];
SaveTheNe[fileName_,numberOfPoints_] := Do[Write[fileName,CForm[SequenceForm[j*Emax/numberOfPoints,blank,Ne[j*Emax/numberOfPoints]]]],{j,1,numberOfPoints}];

PhiE[E_] := F*E/(T^2)*Exp[-E/T]; (* E in units of K *)
PhiEt[eupper_] := NIntegrate[PhiE[x],{x,0,0,eupper}];
FluxE := PhiEt[Emax];
Vbar := FluxE/DensityE/100;
PlotPhiE := Plot[PhiE[x], {x,0,Emax},  PlotLabel->ColumnForm[{labelE2,label2,labelE4,labelE5},Center],AxesLabel->{"Energy (K)", "Flux (#/cm^2/s/K)"}];

(* enter the commands here*)
(*PlotNv
PlotNl*)
(*PlotNE*)
(*PlotPhiE*)
(*SaveTheNe["UniMax_150K",100]
Close["UniMax_150K"]
DensityE*)
(*Close["UniMax_150K"]*)
(*BrightL[8.9]*)
(*SaveTheNv["MaxV_300_UCN",10]*)
(*DensityE
Vbar*)
(*Close["MaxV_20_UCN"]*)


(*SaveTheNv["MaxV_20_UCN",10]*)
(*Close["Max_50%_20"]*)
(*Ne[12]*)
(*SaveTheNe["Max_10%_80",200]*)
(*Close["Max_10%_80"]*)
(*DensityE*)
(*SaveTheNe["Max_100%_100",200]*)
(*Ne[72]
PhiE[72]*)
(*PlotPhiE*)
(*Net[0.001]*)
(*DensityV*)
DensityE
PlotPhiE

-II. By phonon scattering in a Debye Solid:

Clear["Global`*"];
(*The following 3 quantitities are required when you are not reading in a flux spectrum.*)
Tn =25.; (* Maxwell temperature of the incident neutron flux in K. *);
rho0= 1.835 10^9; (* incident neutron density, neutrons/cm^3  *)
WeightN = 1; (* spectral weight factor (0 -> 1.) for Tn  - useful when doing a composite of Maxwellians*)

ThetaDebye = 110; (* Debye temperature of the solid *)
Vucn= 8.; (*Definition of UCN = Vmax for UCN in m/s *)

FM = WeightN*rho0*alpha[Tn]*2/Sqrt[N[Pi]]*100; (* Maxwell distribution at Tn for incident neutrons *) 

F = FM;

ReferenceFlux = 12.28; (* UCN Flux from a gas converter at 300 K in n/cm^2/s *)

kbzm = 1.3807 10^-23; (* Maxwell- Boltzman constant in MKS units *)
mass = 1.6749 10^-27; (* Neutron mass in kg *)
lEConstant = 30.84 ; (*  E^1/2 Lamda = 30.84 for E in K and Lamda in Angstroms for neutrons *)
lVConstant = 3957 ; (*  V Lamda = 3957 for V in m/s and Lamda in Angstroms for neutrons *)

(* for D2*)
mu = 1.9968; (* M/(neutron mass) for the gas "atoms" *)
NucDen = 0.06; (* Nuclear density in units of 10^24 nuclei/cm^3 *)
Sigma0C = 0.907 ; (* cross-section in barns for (incoherent) n scattering from a free deuteron at rest *)
Sigma0T = 3.396 ; (* total cross-section in barns for n scattering from a free deuteron at rest *)
SigmaAv = N[0.0005 2224.]; (* absorption cross section * v. Set =0.0005 at 300 K (v = 2224 m/s) *)
SigmaAe = N[0.0005 17.321]; (* absorption cross section * Sqrt[E]. Set =0.0005 at 300 K (v = 2224 m/s) *)
SigmaAl = N[0.0005 0.562]; (* absorption cross section * 1/lamda. Set =0.0005 at 300 K (lamda = 1.779 Angstroms) *)


(* for H2
mu = 0.99917; (* M/(neutron mass) for the gas "atoms" *)
NucDen = 0.066; (* Nuclear density in units of 10^24 nuclei/cm^3 *)
Sigma0C = 20.05 ; (* cross-section in barns for (incoherent) n scattering from a free deuteron at rest *)
Sigma0T = 20.49 ; (* total cross-section in barns for n scattering from a free deuteron at rest *)
SigmaAv = N[0.333 2224.]; (* absorption cross section * v. Set =0.0005 at 300 K (v = 2224 m/s) *)
SigmaAe = N[0.333 17.321]; (* absorption cross section * Sqrt[E]. Set =0.0005 at 300 K (v = 2224 m/s) *)
SigmaAl = N[0.333 0.562]; (* absorption cross section * 1/lamda. Set =0.0005 at 300 K (lamda = 1.779 Angstroms) *)
*)

Eucn =  N[0.5 mass Vucn^2/kbzm]; (* Emax for UCN in K*)
DeToDv = N[mass/kbzm];

lMin = N[lVConstant/Vucn];
lMax = 10*lMin;

vf[Ei_] := N[Sqrt[2 kbzm Ei/mass]];
En[v_] := N[mass v^2/2/kbzm];
El[L_] := N[(30.84/L)^2];
Lam[Ei_] := N[30.84/Sqrt[Ei]]; 
alpha[T_] := Sqrt[2 kbzm T/mass]; (* most probable velocity of neutron in a Maxwell Dbn at T *)
X[Ex_,T_] := N[Sqrt[Ex/T]];

label2 := StringForm["Incident flux = `` neutron/cm^2/s at Tn = `` K", F,Tn];
labelV1[T_] := StringForm["N(v) for a Converter (with mu = ``) at `` K",mu,T];
labelV3[T_] := StringForm[" `` neutrons/cm^3 with V < `` m/s", DensityV[T], Vucn];
LlabelV3 := StringForm[" `` neutrons/cm^3 with V < `` m/s", ListDensityV, Vucn];
labelE1[T_] := StringForm["N(E) for a Converter (with mu = ``) at `` K",mu,T];
labelE2[T_] := StringForm["Flux(E) for a Converter (mu = `` ) at `` K",mu,T];
labelE3[T_] := StringForm[" `` neutrons/cm^3 with E < `` K", DensityE[T], Eucn];
LlabelE3 := StringForm[" `` neutrons/cm^3 with E < `` K", ListDensityE, Eucn];
labelE4[T_] := StringForm["`` neutron/cm^2/s with E < `` K",  FluxE[T], Eucn];
LlabelE4 := StringForm["`` neutron/cm^2/s with E < `` K",  ListFluxE, Eucn];
labelL1[T_] := StringForm["N(lamda) for a Converter (with mu = ``) at `` K ",mu,T];
labelL2[T_] := StringForm[ " ``  neutrons/cm^3 with lamda > `` Angstroms", lamda[T], lMin];
LlabelL2 := StringForm[ " ``  neutrons/cm^3 with lamda > `` Angstroms", ListLamda, lMin];
DBlabelV1[T_] := StringForm["Equilibrium N(v) for a Converter (with mu = ``) at `` K",mu,T];
DBlabelV3[T_] := StringForm[" `` neutrons/cm^3 with V < `` m/s", DBDensityV[T], Vucn];
DBLlabelV3 := StringForm[" `` neutrons/cm^3 with V < `` m/s", ListDBDensityV, Vucn];
DBlabelE1[T_] := StringForm["Equilibrium N(E) for a Converter (with mu = ``) at `` K",mu,T];
DBlabelE2[T_] := StringForm["Equilibrium Flux(E) for a Converter (mu = `` ) at `` K",mu,T];
DBlabelE3[T_] := StringForm[" `` neutrons/cm^3 with E < `` K", DBDensityE[T], Eucn];
DBLlabelE3 := StringForm[" `` neutrons/cm^3 with E < `` K", ListDBDensityE, Eucn];
DBlabelE4[T_] := StringForm["`` neutron/cm^2/s with E < `` K",  DBFluxE[T], Eucn];
DBLlabelE4 := StringForm["`` neutron/cm^2/s with E < `` K",  ListDBFluxE, Eucn];
DBlabelL1[T_] := StringForm["Equilibrium N(lamda) for a Converter (with mu = ``) at `` K",mu,T];
DBlabelL2[T_] := StringForm[ " ``  neutrons/cm^3 with lamda > `` Angstroms", DBlamda[T], lMin];
DBLlabelL2 := StringForm[ " ``  neutrons/cm^3 with lamda > `` Angstroms", ListDBLamda, lMin];
blank = "  ";
com = " , ";

EMax[T_,Ei_] := Min[7*Tn,ThetaDebye];
LamMin[T_,l_] := Min[lMin,Lam[EMax[T,El[l]]]];
EMin = Eucn;

Phonons[eps_]:= If[Abs[eps]<=ThetaDebye,3./ThetaDebye^3*eps^2,0]; 

G[eps_,T_]:= N[Phonons[Abs[eps]]/\
			If[Abs[eps]<10^-8,1,eps*(1-Exp[-eps/T])]]; 

(*tau[T_] := N[1/(NIntegrate[1/e Coth[e/2/T] Phonons[e], {e,10^-6,ThetaDebye},AccuracyGoal->5])];*)

Debye ={{0,0.6667},{0.1,0.6267},{0.2,0.5367},{0.3,0.4437},
		{0.4,0.3697},{0.5,0.3011},{0.6,0.2583},{0.7,0.2255},
		{0.8,0.1998},{0.9,0.1791},{1.0,0.1622},{1.2,0.1362},
		{1.4,0.1174},{1.6,0.1040},{1.8,0.09252},{2.0,0.08326},
		{2.5,0.06663},{3.0,0.05555},{3.5,0.04762},{4.0,0.04166}}; (* Turchin,"Slow Neutrons",1965,
		                                                             Appendix 3. page 299 *)
fnDebye = Interpolation[Debye];
tau[T_]:= fnDebye[T/ThetaDebye]*ThetaDebye;



SigmaE[T_,Ei_,Ef_]:= Return[argH = 2*Sqrt[Ei*Ef]/mu/tau[T]; 
								Sig = Sigma0T*(1+mu)^2/mu^2*tau[T]/2/Ei*G[Ei-Ef,T]*Exp[-(Ei+Ef)/mu/tau[T]]\
								       *((Ei+Ef+mu*tau[T])*Sinh[argH] - 2*Sqrt[Ei*Ef]*Cosh[argH]); (*the use of Sigma0T here is the incoherent approximation *)
						  Sig ]; (*this equation is from Yu, Malik and Golub, 1986 *)

SigmaTot[T_,Ei_] := Return[DSigmaCurve[T,Ei];
							DSigF = Interpolation[myList];
							ans = NIntegrate[DSigF[x],{x,EMin,EMax[T,Ei]},AccuracyGoal->5]; 
							ans];
							
SigmaTotCurve[T_] := Return[xmin = EMin;
						xmax = EMax[T,Ei];
						delta = .5;
						myList2 = {{xmin,SigmaTot[T,xmin]}};
						For[z2=xmin+delta*10^MantissaExponent[xmin/10][[2]],z2<xmax||(z2-xmax)/xmax<.001,z2=z2+delta*10^MantissaExponent[z2/10][[2]],
							myList2 = Append[myList2,{z2,SigmaTot[T,z2]}];	
							];
							Do[Write["SigTotDebye",CForm[SequenceForm[myList2[[kk,1]],blank,myList2[[kk,2]]]]],{kk,1,Length[myList2]}];
	      Close["SigTotDebye"];
						myList2];
						
PlotTheTotSigmaCurve[T_] := {SigmaTotCurve[T];
							LogLogListPlot[myList2, PlotJoined->True,PlotRange->{{.001,200.},{1. 10^-6,1. 10^-3.}},\
							FrameLabel->{" Energy (K)"," Total Inelastic Sigma (barns)",StringForm["Debye Solid (Single phonon) at `` K. ",T]," "},\
							AxesOrigin->{.001,1 10^-6},Frame->True]};

DSigmaCurve[T_,Ei_] := Return[xmin = EMin;
						xmax = EMax[T,Ei];
						delta = 0.5;
						myList = {{xmin,SigmaE[T,Ei,xmin]}};
						For[z=xmin+delta*10^MantissaExponent[xmin/10][[2]],z<xmax||(z-xmax)/xmax<.001,z=z+delta*10^MantissaExponent[z/10][[2]],
							If[z!=Ei,myList = Append[myList,{z,SigmaE[T,Ei,z]}]];	
							];
						myList];
					
DSigmaCurveF[T_,Ef_] := Return[xmin = EMin;
						xmax = EMax[T,Ef];
						delta = 0.5;
						myListF = {{xmin,SigmaE[T,xmin,Ef]}};
						For[z=xmin+delta*10^MantissaExponent[xmin/10][[2]],z<xmax||(z-xmax)/xmax<.001,z=z+delta*10^MantissaExponent[z/10][[2]],
							If[z!=Ef,myListF = Append[myListF,{z,SigmaE[T,z,Ef]}]];	
							];
						myListF];
						
						
PlotTheDSigmaCurve[T_,Ei_] :={DSigmaCurve[T,Ei];
							LogLogListPlot[myList, PlotJoined->True,PlotRange->{{.001,2000},{10^-6,0.1}},\
			FrameLabel->{" Energy (K)"," d(Sigma)/dE (barns/K)",StringForm["Solid D2 at `` K. Incident Energy = `` K ",T,Ei]," "},\
			AxesOrigin->{.001,1 10^-6},Frame->True]};

PlotTheProductionCurve[T_,Ef_] :={DSigmaCurveF[T,Ef];
							LogLogListPlot[myListF, PlotJoined->True,PlotRange->{{.1,1000},{10^-4,1.}},\
			FrameLabel->{" Energy (K)"," d(Sigma)/dE (barns/K)",StringForm["Solid D2 at `` K. Final Energy = `` K ",T,Ef]," "},\
			AxesOrigin->{.1,1 10^-4},Frame->True]};

MFPCurve[T_,Ef_] := Return[xmin = EMin;
						xmax = EMax[T,Ef];
						delta = 0.5;
						myList = {{xmin,MFPScatt[T,xmin,Ef]}};
						For[z=xmin+delta*10^MantissaExponent[xmin/10][[2]],z<xmax||(z-xmax)/xmax<.001,z=z+delta*10^MantissaExponent[z/10][[2]],
							myList = Append[myList,{z,MFPScatt[T,z,Ef]}];	
							];
						myList];

PlotTheMFPCurve[T_,Ef_] :={MFPCurve[T,Ef];
							LogLogListPlot[myList, PlotJoined->True,PlotRange->{{.001,200},{10^4,10^8}},\
			FrameLabel->{" Incident Energy (K)"," E-MFP (cm-K)",StringForm["Solid D2 at `` K. Final Energy = `` K ",T,Ef]," "},\
			AxesOrigin->{.001,1 10^3},Frame->True]};


IntDownE[T_,e_] := NIntegrate[ee*Exp[-ee/Tn]*SigmaE[T,ee,e], {ee,e, EMax[T,e]},AccuracyGoal->5];
IntDownEF[T_,e_] := NIntegrate[FluxEDbn[ee]*SigmaE[T,ee,e], {ee,e, EMax[T,e]},AccuracyGoal->5];
MFPDown[T_,e_] := 1/(NucDen NIntegrate[SigmaE[T,ee,e], {ee,e, EMax[T,e]},AccuracyGoal->5]);
(*MFPScatt[T_,ei_,ef_] := 1/(NucDen NIntegrate[SigmaE[T,ei,ee], {ee,ef*0.99, ef*1.01},AccuracyGoal->5]);*)
MFPScatt[T_,ei_,ef_] := 1/(NucDen SigmaE[T,ei,ef]);

IntDownEDB[T_,e_] := NIntegrate[ee*Exp[-ee/T]*SigmaE[T, ee, e]  , {ee,e,EMax [T,e]},AccuracyGoal->5]; 
IntUpE[T_,e_] := NIntegrate[SigmaE[T, e, eee]  , {eee,1.1*e,EMax[T,e]},AccuracyGoal->5];
MFPA[T_,e_] := 1/(NucDen*(SigmaAv/vf[e]));
MFPUp[T_,e_] := 1/(NucDen*(IntUpE[T,e]));
MFPLoss[T_,e_] := 1/(NucDen*(SigmaAv/vf[e]+IntUpE[T,e]));



IntDownVE[T_,v_] := NIntegrate[vv^3*Exp[-En[vv]/Tn]*SigmaE[T,En[vv],En[v]], {vv,v, vf[EMax[T,En[v]]]},AccuracyGoal->5];
IntDownVEDB[T_,v_] :=  NIntegrate[vv^3*Exp[-En[vv]/T]*SigmaE[T, En[vv],En[v]]  , {vv,v, vf[EMax[T,En[v]]]},AccuracyGoal->5]; 
IntUpVE[T_,v_] := NIntegrate[SigmaE[T, En[v], En[vv]]*vv  , {vv,v, vf[EMax[T,En[v]]]},AccuracyGoal->5];

IntDownLE[T_,lu_] :=   NIntegrate[Exp[-El[lp]/Tn]*SigmaE[T,El[lp],El[lu]]/lp^5, {lp, LamMin[T,lu], lu},AccuracyGoal->5];
IntProdLE[T_,lu_] :=   NIntegrate[Exp[-El[lp]/Tn]*SigmaE[T,El[lp],El[lu]]/lp^5, {lp, LamMin[T,lu], lu-0.17 },AccuracyGoal->5]+NIntegrate[Exp[-El[lp]/Tn]*SigmaE[T,El[lp],El[lu]]/lp^5, {lp,  lu+0.17, lMin },AccuracyGoal->5];
IntDownLDBE[T_,lu_] := NIntegrate[Exp[-El[lpd]/T]*SigmaE[T,El[lpd],El[lu]]/lpd^5, {lpd,LamMin[T,lu], lu},AccuracyGoal->5];
IntUpLE[T_,lu_] :=  NIntegrate[SigmaE[T, El[lu], El[lx]]/lx^3  , {lx,LamMin[T,lu],lu},AccuracyGoal->5];
IntScatLE[T_,lu_] :=  NIntegrate[SigmaE[T, El[lu], El[lx]]/lx^3  , {lx,LamMin[T,lu],lu-0.17},AccuracyGoal->5]+NIntegrate[SigmaE[T, El[lu], El[lx]]/lx^3  , {lx,lu+0.17,lMin},AccuracyGoal->5];


(*  Using Detailed Balancing (Equilibrium), Rigorous Formulation. *)
 DBNv[T_,v_] := N[2 F 10^-2 1/alpha[Tn]^4 DeToDv v IntDownVE[T,v]/(SigmaAv + Exp[En[v]/T]/v*DeToDv*IntDownVEDB[T,v])]; 
 DBDensityV[T_] := NIntegrate[DBNv[T,v] ,{v,0,Vucn}, AccuracyGoal->5];
 PlotTheDBDensityVDistribution[T_] := Plot[DBNv[T,v], {v,Vucn/100,Vucn},PlotLabel->ColumnForm[{DBlabelV1[T],label2,DBlabelV3[T],blank},Center],AxesLabel->{"Velocity (m/s)", "Density (#/cm^3/(m/s))"}];
 SaveTheDBNv[fileName_,T_,numberOfPoints_] := Do[Write[fileName,CForm[SequenceForm[j*Vucn/numberOfPoints,blank,DBNv[T,j*Vucn/numberOfPoints]]]],{j,1,numberOfPoints}];  
 ListPlotDBDensityV[T_,numberOfPoints_]:=
  	{ n = numberOfPoints;
      ListDBNv = Table[{j*Vucn/numberOfPoints,DBNv[T,j*Vucn/numberOfPoints]},{j,1,n}];
      ListDBNvFit = Interpolation[ListDBNv];
      ListDBDensityV = NIntegrate[ListDBNvFit[x],{x,0,Vucn}];
      ListPlot[ListDBNv,PlotJoined->True,PlotLabel->ColumnForm[{DBlabelV1[T],label2,DBLlabelV3,blank},Center],AxesLabel->{"Velocity (m/s)", "Density (#/cm^3/(m/s))"}];
     };
 
 DBNe[T_,e_] := N[Sqrt[mass/2/kbzm] F/Tn^2 10^-2 IntDownE[T,e]/(SigmaAe + Exp[e/T]/Sqrt[e]*IntDownEDB[T,e])];
 DBDensityE[T_] := NIntegrate[DBNe[T,ee] ,{ee,0,Eucn}, AccuracyGoal->5];
 PlotTheDBDensityEDistribution[T_] := Plot[DBNe[T,e], {e,Eucn/100,Eucn},PlotLabel->ColumnForm[{DBlabelE1[T],label2,DBlabelE3[T],blank},Center],AxesLabel->{"Energy (K)", "Density (#/cm^3/K)"}];
 SaveTheDBNe[fileName_,T_,numberOfPoints_] := Do[Write[fileName,CForm[SequenceForm[j*Eucn/numberOfPoints,blank,DBNe[T,j*Eucn/numberOfPoints]]]],{j,1,numberOfPoints}];  
 ListPlotDBDensityE[T_,numberOfPoints_]:=
  	{ n = numberOfPoints;
      ListDBNe = Table[{j*Eucn/numberOfPoints,DBNe[T,j*Eucn/numberOfPoints]},{j,1,n}];
      ListDBNeFit = Interpolation[ListDBNe];
      ListDBDensityE = NIntegrate[ListDBNeFit[x],{x,0,Eucn}];
      ListPlot[ListDBNe,PlotJoined->True,PlotLabel->ColumnForm[{DBlabelE1[T],label2,DBLlabelE3,blank},Center],AxesLabel->{"Energy (K)", "Density (#/cm^3/K)"}];
     };
 
 DBNl[T_,lu_] := N[4 F 1/alpha[Tn]^4 10^-2 lVConstant^3 lEConstant^2/lu^3 IntDownLE[T,lu]/(SigmaAl + IntDownLDBE[T,lu]*Exp[El[lu]/T]*lu*2*lEConstant^2)]; 
 DBlamda[T_] := NIntegrate[DBNl[T,lz], {lz,lMin,lMax},AccuracyGoal->5];
 PlotTheDBLamdaDistribution[T_] := Plot[DBNl[T,x], {x,lMin,lMax},PlotLabel->ColumnForm[{DBlabelL1[T],label2,DBlabelL2[T],blank},Center],AxesLabel->{"Wave Length (A) ", "Density (#/cm^3/A)"}];
 SaveTheDBNl[fileName_,T_,numberOfPoints_] := Do[Write[fileName,CForm[SequenceForm[N[lMin+(j-1)*(lMax-lMin)/(numberOfPoints-1)],blank,DBNl[T,N[lMin+(j-1)*(lMax-lMin)/(numberOfPoints-1)]]]]],{j,1,numberOfPoints}];  
 ListPlotDBLamda[T_,numberOfPoints_]:=
  	{ n = numberOfPoints;
      ListDBNl = Table[{N[lMin+(j-1)*(lMax-lMin)/(n-1)],DBNl[T,N[lMin+(j-1)*(lMax-lMin)/(n-1)]]},{j,1,n}];
      ListDBNlFit = Interpolation[ListDBNl];
      ListDBLamda = NIntegrate[ListDBNlFit[x],{x,lMin,lMax}];
      ListPlot[ListDBNl,PlotJoined->True,PlotLabel->ColumnForm[{DBlabelL1[T],label2,DBLlabelL2,blank},Center],AxesLabel->{"Wave Length (A) ", "Density (#/cm^3/A)"}];
     };
  
 DBPhiE[T_,e1_] := N[ F  1/Tn^2 Sqrt[e1] IntDownE[T,e1]/(SigmaAe +Exp[e1/T]/Sqrt[e1]*IntDownEDB[T,e1])];
 DBFluxE[T_] := NIntegrate[ DBPhiE[T,e],{e,0,Eucn}, AccuracyGoal->5]; 
 PlotTheDBFluxDistribution[T_] := Plot[DBPhiE[T,x], {x,Eucn/100,Eucn},PlotLabel->ColumnForm[{DBlabelE2[T],label2,DBlabelE4[T],blank},Center],AxesLabel->{"Energy (K)", "Flux (#/cm^2/s/K)"}];
 SaveTheDBPhiE[fileName_,T_,numberOfPoints_] := Do[Write[fileName,CForm[SequenceForm[j*Eucn/numberOfPoints,blank,DBPhiE[T,j*Eucn/numberOfPoints]]]],{j,1,numberOfPoints}];  
 ListPlotDBFlux[T_,numberOfPoints_]:=
  	{ n = numberOfPoints;
      ListDBPhiE = Table[{j*Eucn/numberOfPoints,DBPhiE[T,j*Eucn/numberOfPoints]},{j,1,n}];
      ListDBPhiEFit = Interpolation[ListDBPhiE];
      ListDBFluxE = NIntegrate[ListDBPhiEFit[x],{x,0,Eucn}];
      ListPlot[ListDBPhiE,PlotJoined->True,PlotLabel->ColumnForm[{DBlabelE2[T],label2,DBLlabelE4,blank},Center],AxesLabel->{"Energy (K)", "Flux (#/cm^2/s/K)"}];
     };
 
(* Steady State, Non Equilibrium*) 
  Nv[T_,v_] := N[2 F 10^-2 1/alpha[Tn]^4 DeToDv v IntDownVE[T,v]/(SigmaAv + v*DeToDv*IntUpVE[T,v])];
  DensityV[T_] := NIntegrate[Nv[T,v] ,{v,0,Vucn}, AccuracyGoal->5];
  PlotTheDensityVDistribution[T_] := Plot[Nv[T,v], {v,Vucn/100,Vucn},PlotLabel->ColumnForm[{labelV1[T],label2,labelV3[T],blank},Center],AxesLabel->{"Velocity (m/s)", "Density (#/cm^3/(m/s))"}];
  SaveTheNv[fileName_,T_,numberOfPoints_] := Do[Write[fileName,CForm[SequenceForm[j*Vucn/numberOfPoints,blank,Nv[T,j*Vucn/numberOfPoints]]]],{j,1,numberOfPoints}];  
  ListPlotDensityV[T_,numberOfPoints_]:=
  	{ n = numberOfPoints;
      ListNv = Table[{j*Vucn/numberOfPoints,Nv[T,j*Vucn/numberOfPoints]},{j,1,n}];
      ListNvFit = Interpolation[ListNv];
      ListDensityV = NIntegrate[ListNvFit[x],{x,0,Vucn}];
      ListPlot[ListNv,PlotJoined->True,PlotLabel->ColumnForm[{labelV1[T],label2,LlabelV3,blank},Center],AxesLabel->{"Velocity (m/s)", "Density (#/cm^3/(m/s))"}];
     };
      
  Ne[T_,e_] := N[Sqrt[mass/2/kbzm] F/Tn^2 10^-2 IntDownE[T,e]/(SigmaAe + Sqrt[e]*IntUpE[T,e])];
  DensityE[T_] := NIntegrate[Ne[T,ee] ,{ee,Eucn/100,Eucn}, AccuracyGoal->5]; 
  ProdE[T_,e_] := N[ F/Tn^2  NucDen IntDownE[T,e]];
  ProdEF[T_,e_] :=  N[ NucDen*IntDownEF[T,e]];
  ProdRate[T_] := NIntegrate[ProdE[T,ee] ,{ee,Eucn/100,Eucn}, AccuracyGoal->5];
  ProdRateF[T_] := NIntegrate[ProdEF[T,ee] ,{ee,Eucn/100,Eucn}, AccuracyGoal->5];
  PlotTheDensityEDistribution[T_] := Plot[Ne[T,e], {e,Eucn/100,Eucn},PlotLabel->ColumnForm[{labelE1[T],label2,labelE3[T],blank},Center],AxesLabel->{"Energy (K)", "Density (#/cm^3/K)"}];
  SaveTheNe[fileName_,T_,numberOfPoints_] := Do[Write[fileName,CForm[SequenceForm[j*Eucn/numberOfPoints,blank,Ne[T,j*Eucn/numberOfPoints]]]],{j,1,numberOfPoints}];  
  ListPlotDensityE[T_,numberOfPoints_]:=
  	{ n = numberOfPoints;
      ListNe = Table[{j*Eucn/numberOfPoints,Ne[T,j*Eucn/numberOfPoints]},{j,1,n}];
      ListNeFit = Interpolation[ListNe];
      ListDensityE = NIntegrate[ListNeFit[x],{x,0,Eucn}];
      ListPlot[ListNe,PlotJoined->True,PlotLabel->ColumnForm[{labelE1[T],label2,LlabelE3,blank},Center],AxesLabel->{"Energy (K)", "Density (#/cm^3/K)"}];
     };
 Nl[T_,lu_] := N[4 F 1/alpha[Tn]^4 10^-2 lVConstant^3 lEConstant^2/lu^3  IntDownLE[T,lu]/(SigmaAl + IntUpLE[T,lu]*2*lEConstant^2/lu)];
  lamda[T_] := NIntegrate[Nl[T,lz], {lz,lMin,lMax},AccuracyGoal->5];
  PlotTheLamdaDistribution[T_] := Plot[Nl[T,x], {x,lMin,lMax},PlotLabel->ColumnForm[{labelL1[T],label2,labelL2[T],blank},Center],AxesLabel->{"Wave Length (A) ", "Density (#/cm^3/A)"}];
  SaveTheNl[fileName_,T_,numberOfPoints_] := Do[Write[fileName,CForm[SequenceForm[N[lMin+(j-1)*(lMax-lMin)/(numberOfPoints-1)],blank,Nl[T,N[lMin+(j-1)*(lMax-lMin)/(numberOfPoints-1)]]]]],{j,1,numberOfPoints}];  
  ListPlotLamda[T_,numberOfPoints_]:=
  	{ n = numberOfPoints;
      ListNl = Table[{N[lMin+(j-1)*(lMax-lMin)/(n-1)],Nl[T,N[lMin+(j-1)*(lMax-lMin)/(n-1)]]},{j,1,n}];
      ListNlFit = Interpolation[ListNl];
      ListLamda = NIntegrate[ListNlFit[x],{x,lMin,lMax}];
      ListPlot[ListNl,PlotJoined->True,PlotLabel->ColumnForm[{labelL1[T],label2,LlabelL2,blank},Center],AxesLabel->{"Wave Length (A) ", "Density (#/cm^3/A)"}];
     };
     
  NlCold[T_,lu_] := N[4 F 1/alpha[Tn]^4 10^-2 lVConstant^3 lEConstant^2/lu^3  IntProdLE[T,lu]/(SigmaAl + IntScatLE[T,lu]*2*lEConstant^2/lu + 2/Ws/lu/NucDen)];
   (*NlCold[T_,lu_] := N[4 F 1/alpha[Tn]^4 10^-2 lVConstant^3 lEConstant^2/lu^3  IntProdLE[T,lu]/(SigmaAl + IntScatLE[T,lu]*2*lEConstant^2/lu )];*)
  BrightL[T_,l_] := NlCold[T,l] lVConstant/l 100/4/N[Pi]; (* Brightness in n/cm^2/s/Å/sr at l Å*)
  ProdL[T_,lu_] := N[4 NucDen F 1/alpha[Tn]^4  lVConstant^4 lEConstant^2/lu^3  IntProdLE[T,lu]];
MinBrightL[T_,l_] := ProdL[T,l]*Ws/24/N[Pi];
                        
  						
  PhiE[T_,e1_] := N[ F  1/Tn^2 Sqrt[e1] IntDownE[T,e1]/(SigmaAe +Sqrt[e1]*IntUpE[T,e1])];
PhiECold[T_,e1_] := N[ F  1/Tn^2 Sqrt[e1] IntDownE[T,e1]/(SigmaAe +Sqrt[e1]*IntUpE[T,e1]+2*Sqrt[e1]/Ws/NucDen)];(*INCLUDE SIZE EFFECTS FOR UCN*)
 (*FluxE[T_] := NIntegrate[ PhiE[T,e],{e,0,Eucn}, AccuracyGoal->5];*)
  FluxE[T_] := PhiE[T,Eucn]*Eucn/2;
  PlotTheFluxDistribution[T_] := Plot[PhiE[T,x], {x,Eucn/100,Eucn},PlotLabel->ColumnForm[{labelE2[T],label2,labelE4[T],blank},Center],AxesLabel->{"Energy (K)", "Flux (#/cm^2/s/K)"}];
  SaveThePhiE[fileName_,T_,numberOfPoints_] := Do[Write[fileName,CForm[SequenceForm[j*Eucn/numberOfPoints,blank,PhiE[T,j*Eucn/numberOfPoints]]]],{j,1,numberOfPoints}];  
  ListPlotFlux[T_,numberOfPoints_]:=
  	{ n = numberOfPoints;
      ListPhiE = Table[{j*Eucn/numberOfPoints,PhiE[T,j*Eucn/numberOfPoints]},{j,1,n}];
      ListPhiEFit = Interpolation[ListPhiE];
      ListFluxE = NIntegrate[ListPhiEFit[x],{x,0,Eucn}];
      ListPlot[ListPhiE,PlotJoined->True,PlotLabel->ColumnForm[{labelE2[T],label2,LlabelE4,blank},Center],AxesLabel->{"Energy (K)", "Flux (#/cm^2/s/K)"}];
     };
  
  SaveTheDensityEGainCurve[fileName_,Tmax_,Tmin_,numberOfPoints_] := Do[Write[fileName,CForm[SequenceForm[N[Tmin+(j-1)*(Tmax-Tmin)/(numberOfPoints-1)],blank,DensityE[N[Tmin+(j-1)*(Tmax-Tmin)/(numberOfPoints-1)]]]]],{j,1,numberOfPoints}];  
  SaveTheFluxGainCurve[fileName_,Tmax_,Tmin_,numberOfPoints_] := Do[Write[fileName,CForm[SequenceForm[N[Tmin+(j-1)*(Tmax-Tmin)/(numberOfPoints-1)],blank,FluxE[N[Tmin+(j-1)*(Tmax-Tmin)/(numberOfPoints-1)]]/ReferenceFlux]]],{j,1,numberOfPoints}]; 

ReadFluxDbn[filename_]:={Clear[FluxData,instream2,FluxEDbn];
                         instream2 := OpenRead[filename];
                          FluxData  := ReadList[instream2,{Number,Number}];
                          FluxEDbn = Interpolation[FluxData,InterpolationOrder -> 1];
                          instream2:=Close[filename];
                         }; 
(*******************************************************************************************************************)
(* enter the commands here*)

(* the following command will calculate the production rate in SD2 at 5 K for 
 a Maxwellian neutron spectrum with a density of rho0 and a temperature of Tn,
 both specified in the input list above.*)
 
 (*ProdRate[5]*)

(* the following 2 commands will read in a cold neutron spectrum
 and calculate the UCN production rate in SD2 at 5 K*)

ReadFluxDbn["ColdFlux.txt"](*this is a text file containing the neutron energy spectrum.
                            Where column 1 = energy in K at the center of 
                            the bin, space, column2=n/cm^2/microC/K.
                            These are derived from an MCNP(X) tally (histogram)
                            of E(Mev)at the bottom of the bin and n/cm^2/p in the bin*)
ProdRateF[5](* for a flux spectrum as defined in the previous comment, the units of ProdRateF are
               UCN/cm^3/microC*)