1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45

from scipy.integrate import quad
from scipy import integrate
from scipy.optimize import fsolve
import pylab as pl
import numpy as np

# Variables.
boltzmann_const = 1.38e-23
planck_const = 6.62e-34
hbar = planck_const / ( 2 * np.pi )
transition_temp = 9.2
gap_energy_at_zero_kelvin = 3.528 / ( 2 * transition_temp * boltzmann_const )
debye_freq = ( 296 * boltzmann_const ) / hbar

# For subtracting from root_of_integral
a_const = np.log( ( 1.13 * hbar * debye_freq ) / ( boltzmann_const * transition_temp) )
# For simplifying function f.
b_const = ( hbar * debye_freq ) / ( 2 * boltzmann_const)


def f( coherence_length, temp ):
    # Defines the equation whose integral will have its roots found. Epsilon = coherence length. Delta = Gap energy.
    
    squareRoot = np.sqrt( coherence_length*coherence_length + gap_energy*gap_energy )
    return np.tanh( ( ( b_const / temp ) * squareRoot ) / squareRoot )


def integrate( coherence_length, temp ):
    # Integrates equation f with respect to E, between 0 and 1. 

    return integrate.quad( f, 0, 1, args = ( temp, ) )[0]


def root_of_integral( temp ):
    # Finds the roots of the integral with a guess of 0.01.
 
   return fsolve( integrate, 0.01, args = ( temp, ) )


def gap_energy_values( temp ):
    # Subtracts a_const from each root found, to obtain the gap_energy_values.

    return root_of_integral( temp ) - a_const

#Plot temperature values against the gap_energy_values obtained.