import numpy as np import math import scipy.constants as sc from scipy.optimize import root import matplotlib.pyplot as plt from scipy.optimize import root from scipy.misc import derivative from scipy import optimize

def expected_particle_number(n, μ, T): nsq = n[0]*n[0] + n[1]*n[1] + n[2]*n[2] return 1/(np.exp((nsq-μ)/(T))-1)

delta = 1e-8 def N_overall(μ, T, tol): nmax = np.sqrt((np.log(1/tol+1)*T+μ)/3) if np.isnan(nmax): nmax = np.nan_to_num(nmax) nmax = int(nmax.real) +1 nval = [expected_particle_number(np.array([i+1,j+1,k+1]),μ,T) for i in range (nmax) for j in range(nmax) for k in range(nmax)] return np.sum(nval)

xval=np.linspace(-2,3.1,1000) temp=[0.1,0.5,1,2,5] yval=[ [N_overall(x,T,delta) for x in xval] for T in temp] plt.style.use('fivethirtyeight') fig = plt.figure(figsize=(15, 7)) ax = fig.add_subplot(1,1,1) plt.title("Expected Particle Number $ N(\mu) $", fontsize=18) plt.xlabel("Chemical Potential") plt.ylabel("N") plt.ylim(0,4) plt.xlim(-2,3.1) for i in range(len(temp)): plt.plot(xval,yval[i],label=temp[i]) plt.grid(which='major') plt.grid(which='minor', linestyle = ':') plt.grid() plt.legend() plt.show()

def μ(n,T): tol=1e-8 guess=np.linspace(1,2.99,10) guess=np.flipud(guess) for x0 in guess: rs =root(lambda μ: N_overall(μ,T,tol)-n, x0, method="lm") if rs.success: return rs.x[0] print("not found")

def tbe(n): return (4*n/(2.315*np.pi))**(2/3)

nvalues = np.arange(1,9,1) tvalues = [np.linspace(0.01*tbe(n), 2*tbe(n) ,30) for n in nvalues] μvalues = np.zeros((nvalues.size, tvalues[0].size))

for i in range (nvalues.size): print(f"T_BE(n={nvalues[i]})=",tbe(nvalues[i])) μvalues[i] = [μ(nvalues[i],T) for T in tvalues[i]] plt.plot(tvalues[i]/tbe(nvalues[i]), μvalues[i], label=f"n={nvalues[i]}",marker=".") plt.style.use('ggplot') plt.title("Chem. Potential", fontsize=18) plt.xlabel("$ T/T_{BE} $") plt.ylabel("$\mu(N,T/T_{BE})$") plt.grid(which='major') plt.grid(which='minor', linestyle = ':') plt.grid() plt.legend() plt.show()

lowest=[expected_particle_number([1,1,1],μvalues[i],tvalues[i]) for i in range(len(μvalues))] secondlowest=[expected_particle_number([2,1,1],μvalues[i],tvalues[i]) for i in range(len(μvalues))]

fig = plt.figure(figsize=(15, 7)) plt.style.use('fivethirtyeight') ax = fig.add_subplot(1,2,1) plt.title("Lowest State $ N_{1,1,1} $", fontsize=18) plt.xlabel("$T/T_{BE}$") plt.ylabel("$ N_{1,1,1} $") plt.ylim(0,max(nvalues)+0.1) for i in range(nvalues.size): plt.plot(tvalues[i]/tbe(nvalues[i]),lowest[i],label=f"n={nvalues[i]}") plt.grid(which='major') plt.grid(which='minor', linestyle = ':') plt.grid() plt.legend() ax = fig.add_subplot(1,2,2) plt.title("Lowest State $ N_{2,1,1} $", fontsize=18) plt.xlabel("$T/T_{BE}$") plt.ylabel("$ N_{1,1,1} $") for i in range(nvalues.size): plt.plot(tvalues[i]/tbe(nvalues[i]),secondlowest[i],label=f"n={nvalues[i]}") plt.grid(which='major') plt.grid(which='minor', linestyle = ':') plt.grid() plt.legend() plt.show()

xvalues = np.linspace(1,100,100) for T in [0.5, 1.5]: plt.plot( xvalues, [expected_particle_number(n=[1,1,1], μ=μ(n=N,T=T*tbe(N)), T=T*tbe(N))/N for N in xvalues], label=f"T = {T}" ) #plt.plot(xvalues, [expected_particle_number()/N for N in xvalues], label=f"N={nvalues[i]})") plt.vlines(x = 0.5, ymin = 0, ymax = 1) #plt.vlines(x = 1.5, ymin = 0, ymax = 1) plt.xlabel("$ T/T_{BE} $") plt.ylabel("$ N_0/N $") plt.title("$ N_{(0,0,0)}/N $") plt.grid(which='major') plt.grid(which='minor', linestyle = ':') plt.xlim(1,100) plt.grid() plt.legend() plt.show()

def energy(n): return np.dot(n,n)

def system_energy(μ, T): tolerance = 1e-5 N = 0 nmax = 2 N_new = 1 difference = 1 while (difference>tolerance): n_values = [energy([i,j,k])*expected_particle_number([i,j,k], μ, T) for i in range(1,nmax+1) for j in range(1,nmax+1) for k in range(1,nmax+1)] N_new = np.sum(n_values) difference = N_new-N nmax = nmax+1 N = N_new return N

system_energy(-3,1)#test

def my_heat_capacity(T, N): def f(T): return system_energy(μ(N,T), T) return derivative(func=f, x0=T, dx=1e-5)

n_values=np.arange(2,6,1) t_values=np.linspace(0.1,10,50) heat_capacities = [] for N in n_values: TBE = tbe(N) heat_capacities.append([my_heat_capacity(temp * TBE, N) for temp in t_values])

plt.style.use('fivethirtyeight') plt.figure(figsize=(30, 14)) plt.xlabel('$T/T_{BE}$') plt.ylabel('$C_N$') plt.title('Heat Capacity $C_N$') for i in range(n_values.size): plt.plot( t_values, heat_capacities[i], label = f"N ={n_values[i]}" ) plt.legend() plt.show()

specific_heat_capacities = [] for N in n_values: TBE = tbe(N) specific_heat_capacities.append([my_heat_capacity(temp*TBE, N)/N for temp in t_values]