Adds c term and tests

do.sh

@@ -5,14 +5,17 @@
 # shellcheck disable=SC2039
 set -Eeuo pipefail # -e "Automatic exit from bash shell script on error"  -u "Treat unset variables and parameters as errors"
 
+_install_location() {
+	# shellcheck disable=SC2016
+	wolframscript -c 'FileNameJoin[{StringReplace[$UserBaseDirectory, "\\" -> "/"], "Applications", "owenScalapinoFreeEnergy"}, OperatingSystem -> "Unix"]' | tr -d '\r\n'
+}
+
 install() {
 	# shellcheck disable=SC2039
 	local install_loc
-	# shellcheck disable=SC2016
-	install_loc=$(wolframscript -c 'FileNameJoin[{StringReplace[$UserBaseDirectory, "\\" -> "/"], "Applications", "owenScalapinoFreeEnergy"}, OperatingSystem -> "Unix"]' | tr -d '\r\n')
+	install_loc=$(_install_location)
 	echo "Installing to [$install_loc]"
 	mkdir -p "$install_loc"
-	echo $install_loc | hexdump
 	cp src/wl/owenScalapinoFreeEnergy.wl "$install_loc"
 }
 

src/wl/owenScalapinoFreeEnergy.wl

@@ -2,6 +2,10 @@ BeginPackage["owenScalapinoFreeEnergy`"];
 (* Exported symbols added here with SymbolName::usage *)
 
 fs::usage = "fs[energy, delta, temp, mustar] is the distribution function for quasiparticles.";
+entropy::usage = "entropy[delta, temp, mustar, n0] returns superconducting entropy";
+ATerm::usage = "Calculates A term from notes";
+BTerm::usage = "B term";
+CTerm::usage = "C term";
 
 Begin["`Private`"];
 
@@ -9,7 +13,45 @@ Begin["`Private`"];
  * fs is the quasiparticle distribution function for the superconducting state.
  * Here, temp, energy, delta and mustar are all given in the same units.
  *)
-fs[energy_, delta_, temp_, mustar_] := 1 / (Exp[(Sqrt[energy^2 + delta^2] - mustar)/temp] + 1);
+fs[energy_, delta_, temp_, mustar_] := 1 / (Exp[(Sqrt[energy^2 + delta^2] - mustar) / 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}];
+
+(* A from notes *)
+ATermIntegrand[energy_?NumericQ, delta_?NumericQ, temp_?NumericQ,
+	mustar_?NumericQ] :=
+		Abs[energy] * ((
+			Abs@energy /
+					Sqrt[energy^2 + delta^2]) (2 * fs[energy, delta, temp, mustar] -
+				1) + 1);
+
+ATerm[delta_?NumericQ, temp_?NumericQ, mustar_?NumericQ,
+	debyeFreq_?NumericQ] :=
+		2 * NIntegrate[
+			ATermIntegrand[en, delta, temp, mustar], {en, 0, debyeFreq}];
+
+
+(* B term *)
+BTermIntegrand[energy_?NumericQ, delta_?NumericQ, temp_?NumericQ,
+	mustar_?NumericQ] := (1 - 2 * fs[energy, delta, temp, mustar]) /
+		Sqrt[energy^2 + delta^2];
+
+BTerm[delta_?NumericQ, temp_?NumericQ, mustar_?NumericQ,
+	debyeFreq_?NumericQ, N0_?NumericQ,
+	V0_?NumericQ] := (V0) * (delta * N0 *
+		NIntegrate[
+			BTermIntegrand[en, delta, temp, mustar], {en, 0, debyeFreq}] )^2;
+
+CTerm[delta_?NumericQ, temp_?NumericQ, mustar_?NumericQ,
+  N0_?NumericQ] := -entropy[delta, temp, mustar, N0] * temp;
 
 End[]; (* `Private` *)
 

src/wltest/BasicTest.wlt

@@ -7,4 +7,32 @@ VerificationTest[(* 1 *)
 	, TestID -> "f_s test"
 ];
 
+VerificationTest[
+	entropy[2, 1, .5, 1]
+	,
+	2.6166163281810166`
+	, TestID -> "Entropy_s test"
+];
+
+VerificationTest[
+	ATerm[1.3, .5, .2, 100]
+	,
+	7.910233149747325`
+	, TestID -> "A term test"
+];
+
+VerificationTest[
+	BTerm[1.3, .5, .2, 100, 1, .2]
+	,
+	8.058538888403167`
+	, TestID -> "B term"
+];
+
+VerificationTest[
+	CTerm[1.3, .5, .2, 1]
+	,
+	-0.4560502473422899`
+	, TestID -> "C term"
+];
+
 EndTestSection[];