@@ -0,0 +1,1520 @@
+import numpy as np%0Aimport scipy.integrate%0A%0Adef trapz2(func, a, b, eps, npoints, args):%0A    eps = (b-a)*eps%0A    midpoint = 0.5*(a+b)%0A    xpoints = a + np.logspace(np.log10(eps),np.log10(midpoint-a),npoints/2)%0A    ypoints = func(xpoints, *args)%0A    part1 = scipy.integrate.trapz(ypoints,xpoints)%0A%0A    xpoints = b - np.logspace(np.log10(eps), np.log10(b-midpoint), npoints/2)%0A    ypoints = func(xpoints,*args)%0A    part2 = -1*scipy.integrate.trapz(ypoints,xpoints)%0A    return part1 + part2%0A%0Adef mb_integrand_1a(eps,delta,kbt,hf):%0A    f2 = 1/(np.exp((eps+hf)/kbt)+1)%0A    return ((1-2*f2) * np.abs(eps * (eps+hf) + delta**2)/%0A            (np.sqrt(eps**2-delta**2)*np.sqrt((eps+hf)**2-delta**2)))%0A%0Adef mb_integrand_1b(eps, delta, kbt, hf):%0A    f1 = 1/(np.exp(eps/kbt)+1)%0A    f2 = 1/(np.exp((eps+hf)/kbt)+1)%0A    return ((f1-f2) * (eps * (eps+hf) + delta**2)/%0A            (np.sqrt(eps**2-delta**2)*np.sqrt((eps+hf)**2-delta**2)))%0A%0Adef _mb1(delta, kbt, hf):%0A    s1b = 2*trapz2(mb_integrand_1b, delta, 20*delta, 1e-13, 1000, args=(delta,kbt,hf))/hf%0A    if hf %3E 2*delta:%0A        s1a = trapz2(mb_integrand_1a, delta-hf, -delta, 1e-13, 1000, args = (delta,kbt,hf))/hf%0A    else:%0A        s1a = 0%0A    return s1a + s1b%0A%0Adef mb_integrand_2(eps,delta,kbt,hf):%0A    f2 = 1/(np.exp((eps+hf)/kbt)+1)%0A    return ((1-2*f2) * np.abs(eps * (eps+hf) + delta**2)/%0A            (np.sqrt(delta**2-eps**2)*np.sqrt((eps+hf)**2-delta**2)))%0A%0Adef _mb2(delta, kbt, hf):%0A    return trapz2(mb_integrand_2, delta - hf, delta, 1e-10, 6000, args=(delta,kbt,hf))/hf

add basic M-B integral code, based on Peter Day's matlab code