adds equilibrium gap solver, and test for zero temp warnings

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

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
+	)