Kirori Mal College , Delhi University
Statistical Mechanics Programs
INDEX
[Link]. Program Name Date Sign.
Plot the probability of various macrostates in
1 coin-tossing experiment Vs no. of heads. 10/01/2023
Plot the Maxwell speed distribution function at
2 different temperatures in a 3-dimension system 17/01/2023
Plot Specific Heat of Solids w.r.t temperature
a) Dulong-Petit law,
3 b) Einstein distribution function
13/02/2023
c) Debye distribution function
Plot Planck’s law of Black Body radiation w.r.t
4 wavelength and Compare it with Rayleigh-Jean 14/02/2023
and Wein’s distribution law.
Plot the following functions with energy at
different temperatures
5 a) Maxwell-Boltzmann distribution 20/02/2023
b) Fermi-Dirac distribution
c) Bose-Einstein distribution
Plot the distribution of particles w.r.t. energy
(dN/de versus e) in 3 Dimensions for
a) Relativistic and non-relativistic bosons both at
6 high and low temperature.
13/03/2023
b) Relativistic and non-relativistic fermions both
at high and low temperature.
Name -SUNIT SINGH
College Roll No. -2030113 Submitted to Dr. Shyam
Exam Roll No. -20036567131
0
Page
Semester -6th (G-5) Sign…………………….
Experiment 1:
Aim: Plot the probability of various macrostates in coin-tossing experiment versus
number of heads with 4, 8, 16 coins etc.
1
Page
Program:
//Probability of macrostate vs heads in coin tossing experiment.
//sunit singh 2030113
clc;clear;clf;
N= input("Enter the number of seperate experiments ")
for i=1:N
n=input("Enter the No. of coins ")
nom=2^n//microstates
disp('n(h) P(h)')
for j=0:1:n
//No of microstates for a particular macrostate
ns=factorial(n)/(factorial(j)*factorial(n-j));
P(j+1)=ns/nom;//probablity
h(j+1)=j;
disp([h(j+1) P(j+1)])
end
subplot(1,N,i);plot(h,P,'o-','linewidth',3)
title("Number of coins: "+string(n),'fontsize',4);xlabel('Number of heads,
n(h)','fontsize',4);ylabel('Probablity , P(h)','fontsize',4)
end
Output
2
Page
Result:
3
Page
Experiment 2
Plot the Maxwell speed distribution function at different temperatures in a 3-
dimension system. Calculate the average speed, root mean square and most
probable speed
4
Page
Program
//Maxwell speed distributino function at different temperature
//sunit singh
clc;clear;
k=1.38e-23;N=6e23;pi=3.14
v=0:10:3000// speed in m/s
T=300:300:900//temperature in kelvin
n=input("Enter the no. of gases to be considered ")
for r=1:n
name=input("Enter the name of gase: ","string")
M=input("Enter the molar mas in g/mol of : ")
m=M/(N*1000); //molecular weigt (kg)
disp('Temp(K) Vmp Vav Vrms')
for j=1:length(T)
a=m/(2*k*T(j))
for i=1:length(v)
f(j,i)=(4*pi)*((a/pi)^1.5)*(v(i)^2)*exp(-a*v(i)^2);
end
[p,q]=max(f(j,:))
vmp=v(q);
vav=sqrt(4/pi)*vmp;
vrms=sqrt(3/2)*vmp
disp([T(j) vmp vav vrms])
end
subplot(n,1,r)
plot(v',f')
xlabel('v (m/s)','fontsize',4)
ylabel('f(v)','fontsize',4)
legend('T = '+string(T)+' K')
title('Maxwell speed distribution function for ' +string(name),'fontsize',4)
end
Output
5
Page
Result:
6
Page
Experiment 3
Plot Specific Heat of Solids w.r.t temperature
a) Dulong-Petit law, b) Einstein distribution function c) Debye distribution
function
7
Page
Program
//To plot specific heat with temperature using Doulong-Petite , Einstein, and Debye
theory.
//sunit singh 2030113
clc;clear;clf;
k=1.38e-23;N=6e23;
name =input('Enter the value of solid: ','string')
Te=input('Enter the Einstien Temp: ')
Td=input('Enter the Debye Temp: ')
T=0:2:Td//kelvin
for i=1:length(T)
dp(i)=3*N*k//Doulong Petit law
if T(i)==0 then
e(i)=0;d(i)=0
else
x=(Te/T(i))
e(i)=(3*N*k*(x^2)*exp(x))/((exp(x)-1)^2);//einstien model
if T(i)>Td/50
then
I=integrate('((y^4)*exp(y))/((exp(y)-1)^2)','y',0,Td/T(i))
d(i)=9*N*k*I*((T(i)/Td)^3);//Debye mode;
else
d(i)=((12*3.14^4)/5)*N*k*((T(i)/Td)^3);//Debye model (low temp)
end
end
end
plot(T',[dp e d],'linewidth',3)
legend('Doulong Petite law', 'Eisntein mode', 'Debye model',4)
title(string(name)+' ( Td='+string(Td)+'K and Te = '+string(Te)+'K )','fontsize',4)
ylabel('Specific Heat','fontsize',4)
xlabel('Temperature ','fontsize',4)
Output
8
Page
Graphic window
Result:
9
Page
Experiment 4
Plot Planck’s law of Black body radiation w.r.t. wavelength/frequency at different
temperatures. Compare it with
a. Rayleigh-Jeans Law and
b. Wien’s distribution law
for a given temperature.
10
Page
Program
//sunit singh 2030113
clc;clear;
clf;
h=6.626e-34 //planks constant
k=1.38e-23 //boltzmann constant
c=3e8 //speed of light
t=[1000,1200,1400] //temp (KELVIN)
n=1400
r=linspace(0.1e-6,15e-6,n) //wavelength
r0=linspace(5e-6,15e-6,n) //wavelength
for i=1:n
for j=1:3
u(j,i)=(8*%pi*h*c)/((r(i))^5*(exp(h*c/(r(i)*k*t(j)))-1)) //planck's law
11
m(j,i)=(8*%pi*k*t(j))/(r0(i))^4 //rayleigh'Jeans law
v(j,i)=(8*%pi*h*c)/((r(i))^5*(exp(h*c/(r(i)*k*t(j)))))//wein distribution law
Page
end
end
subplot(2,2,1)
plot(r,u(1,:),'r')
plot(r,u(2,:),'g')
plot(r,u(3,:),'b')
xlabel("Wavelength--->","fontsize",3)
ylabel("Energy--->","fontsize",3)
title("Planks law","fontsize",3)
legend("T1=1000","T2=1200","T3=1400")
xgrid
subplot(2,2,2)
plot(r0,m(1,:),'r')
plot(r0,m(2,:),'g')
plot(r0,m(3,:),'b')
xlabel("Wavelength--->","fontsize",3)
ylabel("Energy--->","fontsize",3)
title("Rayleigh law","fontsize",3)
legend("T1=1000","T2=1200","T3=1400")
xgrid
subplot(2,2,3)
plot(r,u(3,:),'b')
plot(r0,m(3,:),'g')
xlabel("Wavelength--->","fontsize",3)
ylabel("Energy --->","fontsize",3)
title("Compare b/w Planks and Rayleigh law","fontsize",3)
legend("(PLANKS)","(RAYLEIGH)")
xstring(1.2e-05,500,"T=1200 K")
xgrid
subplot(2,2,4)
plot(r,u(3,:),'r')
plot(r0,m(3,:),'g')
plot(r,v(3,:),'b--')
xlabel("Wavelength --->","fontsize",3)
ylabel("Energy --->","fontsize",3)
title("Compare b/w Planks,Rayleigh law and Weins law ","fontsize",3)
legend("Planks","Rayleigh Jeans","Weins")
xstring(1.2e-05,450,"T=1200 K")
xgrid
Output
12
Page
Result
13
Page
Experiment 5
Plot the following functions with energy at different temperatures a) Maxwell-
Boltzmann distribution b) Fermi-Dirac distribution c) Bose-Einstein distribution
14
Page
Program
//plot a) Maxwell-Boltzmann b) Fermi-Dirac and c) Bose-Einstein distribution
//Sunit singh 2030113
clc;clear;clf;
e=1.6e-19;k=1.38e-23
u=0;//chemical potential (eV)
E=-0.5:0.001:0.5//energy in eV
T=100:200:1100//Kelvin
To=500//Temp for comparison
dist=['Bose Einstein','Maxwell Boltzman','Fermi Dirac']
a=-1;
for n=1:3
for j=1:length(T)
if To==T(j)
q=j
end
for i=1:length(E)
f(j,i)=1/(exp(((E(i)-u)*e)/(k*T(j)))+a)//general distribution
end
end
a=a+1
subplot(2,2,n)
plot(E',f','linewidt',3)
ylabel('f(E)','fontsize',4)
xlabel('Energy (eV)','fontsize',4)
legend ('T = '+string(T)+' k')
title(string(dist(n))+' distribution for u = '+string(u),'fontsize',4)
C(n,:)=f(q,:)
end
subplot(2,2,4)
plot(E',C','linewidth',4)
ylabel('f(E)','fontsize',4)
xlabel('Energy (eV)','fontsize',4)
legend(string(dist)+' distribution',3)
15
title('Temperature = '+string(T(q))+' K','fonsize',4)
Page
Output
Result
16
Page
Program 6
Plot the distribution of particles w.r.t. energy (dN/de versus e) in 3 Dimensions for
a) Relativistic and non-relativistic bosons both at high and low temperature.
b) Relativistic and non-relativistic fermions both at high and low temperature.
17
Page
Non-Relativistic Bosons
//dN/dE for non-relavastic bosons
//Sunit singh 2030113
clf;clear;clc;
//define constants
e=1.6e-19;kb=1.38e-23;h=6.626e-34;s=1;u=-1;V=1;m=4*1.66e-27
E=0:0.001:0.5//Energy in eV
T=[100 1000]//Temperature in Kelvin
Cn=(2*s+1)*(2*3.14*V*(2*m)^1.5)/(h^3);
for j=1:length(T)
b=1/(kb*T(j));
for i=1:length(E)
g(i)=Cn*(E(i))^0.5;
n(j,i)=1/(exp((E(i)-u)*e*b)-1);
f(j,i)=g(i)*n(j,i)
end
subplot(2,3,j*j)
plot(E',g,'linewidth',3)
ylabel('g(E)','fontsize',4)
xlabel('Energy(eV)','fontsize',4)
subplot(2,3,j*j+1)
plot(E',n(j,:)','linewidth',3);legend(string(T(j))+' K');
ylabel('n(E)','fontsize',4)
xlabel('Energy(eV)','fontsize',4)
title('Non-relativistic Bosons
(s='+string(s)+',u='+string(u)+'eV)','fontsize',5);
subplot(2,3,j*j+2)
plot(E',f(j,:)','linewidth',3)
legend(string(T(j))+' K')
ylabel('dN/dE','fontsize',4)
xlabel('Energy(eV)','fontsize',4)
end
18
Page
Output
Relativistic Bosons
//dN/dE for relavastic bosons
//Sunit singh 2030113
clf;clear;clc;
//define constants
e=1.6e-19;kb=1.38e-23;h=6.626e-34;s=1;u=-1;V=1;c=3e8
E=0:0.001:6//Energy in MeV
T=[10^9 10^10]//Temperature in Kelvin
Cr=(2*s*4*3.14*V)/((h^3)*(c^3));
for j=1:length(T)
b=1/(kb*T(j));
for i=1:length(E)
g(i)=Cr*(E(i))^2;
19
n(j,i)=1/(exp((E(i)-u)*10^6*e*b)-1);
f(j,i)=g(i)*n(j,i)
Page
end
subplot(2,3,j*j)
plot(E',g,'linewidth',3)
ylabel('g(E)','fontsize',4)
xlabel('Energy(MeV)','fontsize',4)
subplot(2,3,j*j+1)
plot(E',n(j,:)','linewidth',3);legend(string(T(j))+' K');
ylabel('n(E)','fontsize',4)
xlabel('Energy(MeV)','fontsize',4)
title('Relativistic Bosons (s='+string(s)+',u='+string(u)+'eV)','fontsize',5);
subplot(2,3,j*j+2)
plot(E',f(j,:)','linewidth',3)
legend(string(T(j))+' K')
ylabel('dN/dE','fontsize',4)
xlabel('Energy(MeV)','fontsize',4)
end
Output
20
Page
Non-Relativistic Fermions
//dN/dE for non-relavastic fermions
//Sunit singh 2030113
clf;clear;clc;
//define constants
e=1.6e-19;kb=1.38e-23;h=6.626e-34;s=0.5;u=1;V=1;m=9.1e-31
E=0:0.001:2//Energy in eV
T=[100 1000]//Temperature in Kelvin
Cn=(2*s+1)*(2*3.14*V*(2*m)^1.5)/(h^3);
for j=1:length(T)
b=1/(kb*T(j));
for i=1:length(E)
g(i)=Cn*(E(i))^0.5;
n(j,i)=1/(exp((E(i)-u)*e*b)+1);
f(j,i)=g(i)*n(j,i)
end
end
subplot(1,3,1)
plot(E',g,'linewidth',3)
ylabel('g(E)','fontsize',4)
xlabel('Energy(eV)','fontsize',4)
subplot(1,3,2)
plot(E',n,'linewidth',3);legend(string(T)+' K');
ylabel('n(E)','fontsize',4)
xlabel('Energy(eV)','fontsize',4)
title('Non-relativistic fermions
(s='+string(s)+',u='+string(u)+'eV)','fontsize',5);
subplot(1,3,3)
plot(E',f','linewidth',3);legend(string(T)+' K')
ylabel('dN/dE','fontsize',4)
xlabel('Energy(eV)','fontsize',4)
Output
21
Page
Relativistic Fermions
//dN/dE for relavastic fermions
//Sunit singh 2030113
clf;clear;clc;
//define constants
e=1.6e-19;kb=1.38e-23;h=6.626e-34;s=0.5;u=1;V=1;c=3e8
E=0:0.001:2//Energy in MeV
T=[10^8 10^9]//Temperature in Kelvin
Cr=(2*s*4*3.14*V)/((h^3)*(c^3));
for j=1:length(T)
b=1/(kb*T(j));
for i=1:length(E)
g(i)=Cr*(E(i))^2;
22
n(j,i)=1/(exp((E(i)-u)*10^6*e*b)+1);
Page
f(j,i)=g(i)*n(j,i)
end
end
subplot(1,3,1)
plot(E',g,'linewidth',3)
ylabel('g(E)','fontsize',4)
xlabel('Energy(MeV)','fontsize',4)
subplot(1,3,2)
plot(E',n','linewidth',3);legend(string(T)+' K');
ylabel('n(E)','fontsize',4)
xlabel('Energy(MeV)','fontsize',4)
title('Relativistic fermions
(s='+string(s)+',u='+string(u)+'eV)','fontsize',5);
subplot(1,3,3)
plot(E',f','linewidth',3);legend(string(T)+' K')
ylabel('dN/dE','fontsize',4)
xlabel('Energy(MeV)','fontsize',4)
23
Page