Add normal state function

src/wl/owenScalapinoFreeEnergy.wl

@@ -7,7 +7,8 @@ ATerm::usage = "Calculates A term from notes";
 BTerm::usage = "B term";
 CTerm::usage = "C term";
 DTerm::usage = "D term";
-Fs::usage "Free energy for sc case";
+Fs::usage = "Free energy for sc case";
+Fn::usage = "Free energy for n case";
 
 Begin["`Private`"];
 
@@ -17,28 +18,37 @@ Begin["`Private`"];
  *)
 fs[energy_, delta_, temp_, mustar_] := 1 / (Exp[(Sqrt[energy^2 + delta^2] - mustar) / temp] + 1);
 
+ffree[energy_, temp_] := 1 / (Exp[energy/temp] + 1);
+
 (* Entropy calculation *)
 entropyPart[delta_?NumericQ, temp_?NumericQ, mustar_?NumericQ,
 	energy_?NumericQ] :=
 		With[{f =
 				fs[energy, delta, temp, mustar]}, (f *
 				Log[f]) + (1 - f) * (Log[1 - f])];
+
 entropy[delta_?NumericQ, temp_?NumericQ, mustar_?NumericQ,
 	N0_?NumericQ] := - 2 * N0 *
 		NIntegrate[entropyPart[delta, temp, mustar, e], {e, 0, Infinity}];
 
+entropyPartNormal[temp_?NumericQ, energy_?NumericQ] :=
+		With[{f = ffree[energy, temp]}, (f * Log[f]) + (1 - f) *(Log[1 - f])];
+
+entropyNormal[temp_?NumericQ, N0_?NumericQ] := - 2 * N0 *
+		NIntegrate[entropyPartNormal[temp, e], {e, 0, Infinity}];
+
 (* A from notes *)
 ATermIntegrand[energy_?NumericQ, delta_?NumericQ, temp_?NumericQ,
 	mustar_?NumericQ] :=
-		Abs[energy] * ((
+		Abs[energy] * (
 			Abs@energy /
 					Sqrt[energy^2 + delta^2]) (2 * fs[energy, delta, temp, mustar] -
-				1) + 1);
+				1);
 
 ATerm[delta_?NumericQ, temp_?NumericQ, mustar_?NumericQ,
 	debyeFreq_?NumericQ] :=
-		2 * NIntegrate[
-			ATermIntegrand[en, delta, temp, mustar], {en, 0, debyeFreq}];
+		NIntegrate[
+			ATermIntegrand[en, delta, temp, mustar], {en, -debyeFreq, debyeFreq}];
 
 
 (* B term *)
@@ -67,6 +77,9 @@ Fs[delta_?NumericQ, temp_?NumericQ, mustar_?NumericQ,
 				CTerm[delta, temp, mustar, N0] +
 				DTerm[delta, temp, mustar, debyeFreq, N0];
 
+Fn[temp_?NumericQ, debyeFreq_?NumericQ, N0_?NumericQ] :=
+		(2 * N0) *NIntegrate[en * ffree[en, temp], {en, -debyeFreq, debyeFreq}] - temp * entropyNormal[temp, 2 * N0];
+
 End[]; (* `Private` *)
 
 EndPackage[];

src/wltest/BasicTest.wlt

@@ -17,7 +17,7 @@ VerificationTest[
 VerificationTest[
 	ATerm[1.3, .5, .2, 100]
 	,
-	7.910233149747325`
+	-9992.089994992271`
 	, TestID -> "A term test"
 ];
 
@@ -45,8 +45,15 @@ VerificationTest[
 VerificationTest[
 	Fs[.5, .5, .247777, 50, 1, .2]
 	,
-	1.7128117685372855
+	-2498.287187445874
 	, TestID -> "Fs Test"
 ];
 
+VerificationTest[
+	Fs[0, .5, 0, 100, 1, 1]
+	, Fn[.5, 100, 1]
+	, TestID -> "Fs == Fn test"
+	, SameTest -> approxEqual
+];
+
 EndTestSection[];