human
I am working on the file "ElectroOptics/Main.ipynb".
The first task is:
2. Curve Fitting
Langevin-Debye Model
$$ U(\theta, \phi, E, \rho_0) = -\rho_0E\:sin\:\theta\:cos\:\phi\:(1+\alpha E\:cos\:\phi) $$
$$ tan(2\psi) =  \frac{\langle sin\:2\theta\:cos\:{\phi}\rangle}{\langle{cos}^2{\theta} - {sin}^2{\theta}\:{cos}^2{\phi}\rangle} $$
$$ \langle X \rangle = \int_{\theta_{min}}^{\theta_{max}} \int_0^{2\pi} X(\theta,\phi) f(\theta, \phi) sin(\theta) d\theta d\phi   $$
$$ f(\theta,\phi) = \frac{e^{\frac{-U}{k_b T}}}{\int_{\theta_{min}}^{\theta_{max}} \int_0^{2\pi} e^{\frac{-U}{k_b T}} \:sin\:{\theta}\: d\theta d\phi}  $$
$$ tan(2\psi) = \frac{\int_{\theta_{min}}^{\theta_{max}} \int_0^{2\pi} {sin\:{2\theta}\:cos\:{\phi}}\:e^{\frac{-U}{k_bT}}\:sin\:{\theta}\: d\theta d\phi}{\int_{\theta_{min}}^{\theta_{max}} \int_0^{2\pi} ({{cos}^2{\theta} - {sin}^2{\theta}\:{cos}^2{\phi}})\:e^{\frac{-U}{k_bT}}\:sin\:{\theta}\: d\theta d\phi}  $$
Computing
Can you write Python program for it?

Can you provide a the Python program you requested?
assistant
#Here is the Python program you requested.

#Test values at 40C
E_array = np.array(Fields[-1])
T = TempsK[-1]
PSI = np.array(Angles[-1])

#Phenomenological points from Shen paper at 40C
p0 = 619 #Debye
alpha = .338 # Micron/Volt

#Convert to SI
p0 = p0 * DebyeToSI
alpha = .338 * 1e-6
E = E_array[1]*1e6

#Combine p0 and kb to a single more reasonable number
p0k = p0/kb

#Scaling try: Puts E, alpha, and p0k all on the same order of magnitude ~200.
E = E*1e-4
alpha = alpha*1e9
p0k = p0k*1e6
E_array = E_array*1e2

def compute_psi(E,p0k,alpha):
    """Computes the tilt angle(psi) via the tan(2psi) equation at one temperature.
    Returns
    -------
    Float:
    The statistical tilt angle with conditions T,p0k,alpha,E"""
    
    def Boltz(theta,phi,T,p0k,alpha,E):
        """Compute the integrand for the Boltzmann factor.
        Returns
        -------
        A function of theta,phi,T,p0k,alpha,E to be used within dblquad
        """
        return np.exp((1/T)*p0k*E*np.sin(theta)*np.cos(phi)*((1*1e-2)+(alpha*E*np.cos(phi)*1e-7)))*np.sin(theta)
        #Factors: 1e-2 and 1-e7 and to cancel out the rescaling, of E,alpha, and p0k, done in the above cell
    
    
    def numerator(theta,phi,T,p0k,alpha,E):
        boltz = Boltz(theta,phi,T,p0k,alpha,E)
        return np.sin(2*theta)*np.cos(phi)*boltz

    def denominator(theta,phi,T,p0k,alpha,E):
        boltz = Boltz(theta,phi,T,p0k,alpha,E)
        return (np.cos(theta)**2 - ((np.sin(theta)*np.cos(phi))**2))*boltz
    
    avg_numerator, avg_numerator_error = dblquad(numerator, 0, 2*np.pi, lambda theta1: thetamin, lambda theta2: thetamax,args=(T,p0k,alpha,E))
    
    avg_denominator, avg_denominator_error = dblquad(denominator, 0, 2*np.pi, lambda theta3: thetamin, lambda theta4: thetamax,args=(T,p0k,alpha,E))
    
    psi = np.arctan(avg_numerator / avg_denominator) * (180 /(2*np.pi)) #Convert to degrees from radians and divide by two
    return psi

compute_psi(E,p0k,alpha)

plt.plot(E_array*1e-2,[compute_psi(entry,p0k,alpha) for entry in E_array],label = "Compute...