From 75c1dda67cccc2574b7072a58ae6e7227a270079 Mon Sep 17 00:00:00 2001
From: Deepak <dmallubhotla+github@gmail.com>
Date: Sun, 22 Nov 2020 15:09:37 -0600
Subject: [PATCH] adds equilibrium gap solver, and test for zero temp warnings

---
 pysuperconductor/os_gap_calc.py | 13 +++++++++++++
 tests/test_os_gap_calc.py       | 21 +++++++++++++++++++++
 2 files changed, 34 insertions(+)

diff --git a/pysuperconductor/os_gap_calc.py b/pysuperconductor/os_gap_calc.py
index dd25fa1..db47a94 100644
--- a/pysuperconductor/os_gap_calc.py
+++ b/pysuperconductor/os_gap_calc.py
@@ -1,5 +1,6 @@
 import numpy
 import scipy.integrate as integrate
+import scipy.optimize
 
 
 def energy(freq, delta):
@@ -8,6 +9,8 @@ def energy(freq, delta):
 
 def gap_integrand_function(xi, temp, delta, mu_star):
 	big_e = energy(xi, delta)
+	if temp == 0:
+		return (1 / (2 * big_e)) * (2 * numpy.heaviside(big_e - mu_star, .5) - 1)
 	return numpy.tanh((big_e - mu_star) / (2 * temp)) / (2 * big_e)
 
 
@@ -17,3 +20,13 @@ def gap_integral(temp, delta, mu_star, debye_frequency):
 		# to cut the integral to zero to omega_debyeh
 		return 2 * gap_integrand_function(xi, temp, delta, mu_star)
 	return integrate.quad(integrand, 0, debye_frequency)[0]
+
+
+def equilibrium_gap(temp, debye_frequency, nv):
+	nv_inv = 1 / nv
+	sol = scipy.optimize.root(
+		lambda d: gap_integral(temp, d, 0, debye_frequency) - nv_inv,
+		x0=debye_frequency / (numpy.sinh(nv_inv))
+	)
+	print(sol)
+	return sol.x
diff --git a/tests/test_os_gap_calc.py b/tests/test_os_gap_calc.py
index 1844be9..f032dc0 100644
--- a/tests/test_os_gap_calc.py
+++ b/tests/test_os_gap_calc.py
@@ -1,5 +1,6 @@
 import pysuperconductor.os_gap_calc
 import numpy
+import pytest
 
 
 def test_gap_integrand_function():
@@ -10,9 +11,29 @@ def test_gap_integrand_function():
 	)
 
 
+def test_gap_integrand_function_zero_temp():
+	actual = None
+	with pytest.warns(None) as record:
+		actual = pysuperconductor.os_gap_calc.gap_integrand_function(10, 0, 5, 2)
+	numpy.testing.assert_almost_equal(
+		actual, 0.0447213595499958,
+		decimal=7, err_msg='gap integrand function fails at zero temp', verbose=True
+	)
+	assert not record
+
+
 def test_gap_integral():
 	actual = pysuperconductor.os_gap_calc.gap_integral(1, 5, 2, 10)
 	numpy.testing.assert_almost_equal(
 		actual, 1.390972295468881,
 		decimal=7, err_msg='gap integral is off', verbose=True
 	)
+
+
+def test_equilibrium_gap_zero_temp():
+	actual = pysuperconductor.os_gap_calc.equilibrium_gap(0, 100, .2)
+	numpy.testing.assert_almost_equal(
+		actual, 1.3476505830587,
+		decimal=7, err_msg='gap finding at zero temp equilibrium is wrong',
+		verbose=True
+	)