from math import factorial
#Variable Declaration
aH = 40 #Number of heads
N = 100 #Total events
#Calculations
aT = 100 - aH
We = factorial(N)/(factorial(aT)*factorial(aH))
Wexpected = factorial(N)/(factorial(N/2)*factorial(N/2))
#Results
print 'The observed weight %5.2e compared to %5.2e'%(We,Wexpected)
from sympy import symbols, diff, log
#Varialbe declaration
n = 10000 #Total number of particles
#Calcualtions
def ster(i):
return i*log(i)-i
n1, n2, n3, W = symbols('n1 n2 n3 W',positive=True)
n2 = 5000 - 2*n3
n1 = 10000 - n2 -n3
logW = ster(n) - ster(n1) - ster(n2) - ster(n3)
fun = diff(logW, n3)
dfun = diff(fun, n3)
x0 = 10.0
err = 1.0
while err>0.001:
f = fun.subs(n3,x0)
df = dfun.subs(n3,x0)
xnew = x0 - f/df
err = abs(x0-xnew)/x0
x0 = xnew
x0 = int(x0)
N2 = n2.subs(n3,x0)
N3 = x0
n1 = n1.subs(n3,x0)
N1 = n1.subs(n2,N2)
lnW = logW.subs(n3,N3)
#Results
print 'At maximum value of ln(W)'
print 'Values of N1 : %4d, N2: %4d and N3: %4d '%(N1, N2,N3)
print 'Maximum value of ln(W)= %6d'%lnW
#Variable Declaration
p0 = 0.633 #Probabilities of Energy level 1,2,3
p1 = 0.233
p2 = 0.086
#Calculation
p4 = 1. -(p0+p1+p2)
#Results
print 'Probability of finding an oscillator at energy level of n>3 is %4.3f i.e.%4.1f percent'%(p4,p4*100)
#Variable Declaration
p0 = 0.394 #Probabilities of Energy level 1,2,3
p1by2 = 0.239
p2 = 0.145
#Calculation
p4 = 1. -(p0+p1by2+p2)
#Results
print 'Probability of finding an oscillator at energy level of n>3 is %4.3f'%(p4)
from math import exp
#Variable Declaration
I2 = 208 #Vibrational frequency, cm-1
T = 298 #Molecular Temperature, K
c = 3.00e10 #speed of light, cm/s
h = 6.626e-34 #Planks constant, J/K
k = 1.38e-23 #Boltzman constant, J/K
#Calculation
q = 1./(1.-exp(-h*c*I2/(k*T)))
p2 = exp(-2*h*c*I2/(k*T))/q
#Results
print 'Partition function is %4.3f'%(q)
print 'Probability of occupying the second vibrational state n=2 is %4.3f'%(p2)
#Variable Declaration
B = 1.45 #Magnetic field streangth, Teslas
T = 298 #Molecular Temperature, K
c = 3.00e10 #speed of light, cm/s
h = 6.626e-34 #Planks constant, J/K
k = 1.38e-23 #Boltzman constant, J/K
gnbn = 2.82e-26 #J/T
#Calculation
ahpbyahm = exp(-gnbn*B/(k*T))
#Results
print 'Occupation Number is %7.6f'%(ahpbyahm)