```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. ```