Exercises are from "Computational Physics" by Mark Newman.
Programmed by Sasha Bakker.
from numpy import *
A ball is again dropped from a tower of height h with initial velocity zero.
h = 100.0
g = 9.81 # m/s**2
vi = 0 # m/s
vf = sqrt(vi**2+2*g*h)
t = (vf-vi)/g
print(f"Tower height (m): {h}")
print(f"Ball drop time (s): {round(t,3)}")
A satellite is to be launched into a circular orbit around the Earth so that it orbits the planet once every T seconds.
The motion can be expressed by Kepler's Law of Periods:
Where G is Newton's gravitational constant, M is the mass of Earth, and R is its radius.
G = 6.67*10**-11 # m**3 kg**-1 s**-2
M = 5.97*10**24 # kg
R = 6371*1000 # m
periods = [86400, 86148, 5400, 2700] # s
for T in periods:
h = ((G*M*T**2)/(4*pi**2))**(1/3)-R
print(f"For a period of {T} seconds, the satellite altitude is {round(h/1000,3)} km.")
A satellite cannot orbit with a period of 45 minutes (2700 s). Therefore, a satellite's minimum altitude exists between a period of 90 minutes (5400 s) and 45 minutes. Technicaly, a geosynchoronous satellite is one that orbits the Earth once per sideral day, which is 23.93 hours (86148 s), not 24 hour (86400 s).
hours = array([24, 23.9344696])
length = hours*365
print(f"{round(length[0]-length[1],3)} hours difference in a year with solar days versus siderial days.")
"Because Earth orbits the Sun once a year, the sidereal time at any given place and time will gain about four minutes against local civil time, every 24 hours, until, after a year has passed, one additional sidereal "day" has elapsed compared to the number of solar days that have gone by." Source.
Write a program to perform the inverse operation to that of Example 2.2 for a 2-dimensional space.
x = 4.112
y = 6.05
r = sqrt(x**2 + y**2)
theta = arctan2(y,x)*(180/pi)
print(f"x = {x}, y = {y}")
print(f"r = {round(r,3)} units, theta = {round(theta,3)} degrees")
A spaceship travels from Earth in a straight line at relativistics speed v to another planet x light years away.
print("Choose x and v as a fraction of the speed of light c.")
x = 9.5
v = 0.95
c = 3e8 # m/s
t = x/v # dialated time (s)
t0 = t*sqrt(1-v**2) # proper time (s)
print(f"x = {x} light years, v = {v}c")
print(f"While {round(t0,2)} years pass on the space ship... {t} years pass on Earth!")
Read about time dialation here.
A well-known quantum mechanics problem involves a particle of mass m that encounters a one-dimensional potential step:
The particle with initial kinetic energy E and wavefactor $k_{1} = \sqrt{2mE}/\hbar$ enters from the left and encounters a sudden jump in potential energy of height V at position $x=0$. By solving the Schrödinger equation, one can show that when $E>V$ the particle may either (a) pass the step, in which case it has a lower kinetic energy of $E-V$ on the other side and a correspondingly smaller wavevector of $k_{2}=\sqrt{2m(E-V)}/\hbar$, or (b) it may be reflected, keeping all of its kinetic energy and unchanged wavevector but moving in the opposite direction. The probabilities T and R for transmission and reflection are given by:
Suppose we have a particle with mass equal to the electron mass m and energy 10 eV encountering a potential step of height 9 eV.
V = 9 # eV
E = 10 # eV
m = 9.11e-31 #kg
hbar = 6.5821e-16 # eV s
k1 = sqrt(2*m*E)/hbar
k2 = sqrt(2*m*(E-V))/hbar
T = (4*k1*k2)/(k1+k2)**2
R = ((k1-k2)/(k1+k2))**2
print(f"T = {round(T,3)}, R= {round(R,3)}")
print(f"T + R = {round(T+R,1)} because it represents the probability of all cases.")
The orbit in space of one body around another, such as a planet around the Sun, need not be circular. In general it t akes the form of an ellipse, with the body sometimes closer in and sometimes further out. If yo uare given the distance $l_{1}$ of closest approach that a planet makes to the Sun, also called its perihelion, and its linear velocity $v_{1}$ at perihelion, then any other property of the orbit can be calculated from these two as follows.
l1 = array([1.4710e11, 8.7830e10]) # m
v1 = array([3.0287e4, 5.4529e4]) # m/s
names = array(["Earth", "Halley's comet"])
G = 6.6738e-11 # m**3 k**-1 s**-2
M = 1.9891e30 # kg
c1 = (-2*G*M)/(v1*l1)
c2 = -(v1**2-2*G*M/l1)
roots = array([(-c1+sqrt(c1**2-4*c2))/2,(-c1-sqrt(c1**2-4*c2))/2])
v2 = array([min(roots[0]),min(roots[1])])
l2 = (l1*v1)/v2
a = 0.5*(l1+l2) # semi-major axis
b = sqrt(l1*l2) # semi-minor axis
T = (2*pi*a*b)/(l1*v1)
e = (l2-l1)/(l2+l1)
T = T/31536000
for i in range(len(names)):
print(f"{names[i]}: orbital period = {round(T[i],2)} years, eccentricity = {e[i]}")
The Catalan numbers $C_{n}$ are a sequence of integers 1, 1, 2, 5, 14, 42, 132... that play an important role in quantum mechanics and the theory of disordered systems. They are given by
C = 1
n = 0
print(C)
while C <= 1e9:
C = C*(4*n+2)/(n+2)
print(int(C))
n += 1
Suppose arrays a and b are defined. Compute:
print(b/a+1)print(b/(a+1))print(1/a)a = array([1,2,3,4],int)
b = array([2,4,6,8],int)
print(f"1. The result is {b/a+1}")
print(f"2. The result is {b/(a+1)}")
print(f"3. The result is {1/a}")
In condensed matter physics the Madelung constant gives the total electric potential felt by an atom in a solid. It depends on the chrages on the other atoms nearby and their locations. Consider for instance solid sodium chloride--table salt. The sodium chloride crystal has atoms arranged on a cubic lattice, but with alternating sodium and chlorine atoms, the sodium ones having a single positive charge +e and chlorine ones a single negative charge -e, wehter e is the charge on the electron. If we label each position on the lattice by three integer coordinates ( i , j ,k ), then the sodium atoms fall at positions where i + j + k is even, and the chlorine atoms at positions where i + j + k is odd.
Consider a sodium atom at the origin, i = j = k = 0, and let us calculate the Madelung constant. If the spacing of atoms on the lattice is a, then the distance from the origin to the atom at position ( i , j , k ) is
where M is the Madelung constant, at least approximately (assuming large L).
L = 100 # number of atoms in all directions
a = 1 # spacing between atoms (m)
c = 1.112e-10 # C**2 J**-1 m**-1 (represents 4*pi*epsilon-naught)
e = 1.6022e-19 # C
pos = arange(a,L+a,a)
neg = arange(-L,0,a)
M = 0
for i in range(-L,L+1):
for j in range(-L,L+1):
for k in range(-L,L+1):
if (i==j==k==0):
continue
M += (-1)**(i+j+k)/(sqrt(i**2+j**2+k**2))
print(f"The Madelung constant of NaCl is {abs(round(M,4))}")
In nuclear physics, the semi-empirical mass formula is a formula for calculating the approximate nuclear binding energy B of an atomic nucleu with atomic number Z and mass number A:
a1 = 15.8
a2 = 18.3
a3 = 0.714
a4 = 23.2
def a5calc(A,Z):
if (A%2)!=0:
return 0
elif (Z%2)==0:
return 12.0
else:
return -12.0
def Bcalc(A,Z,a5):
return a1*A-a2*A**(2/3)-a3*(Z**2/A**(1/3))-a4*(A-2*Z)**2/A+a5/A**(1/2)
A = 58
Z = 28
a5 = a5calc(A,Z)
B = Bcalc(A,Z,a5)
N = B/A
print(f"The binding energy of an atom with A = {A} and Z = {Z} has total binding energy "
f"B = {round(B,1)} MeV and binding energy per nucleon B/A = {round(N,2)} MeV")
def stable(Z):
A = arange(Z,3*Z+1,1)
a5 = list()
for a in A:
a5.append(a5calc(a,Z))
B = Bcalc(A,Z,a5)
N = list(B/A)
m = max(N)
i = N.index(max(N))
return (A[i],m)
Z = 28
print(f"The value of A that has the most stable nucleus for Z = {Z} is A = {stable(Z)[0]}")
Z = arange(1,101,1)
n = list()
a = list()
for z in Z:
x = stable(z)
a.append(x[0])
n.append(x[1])
print(f"The maximum binding energy per nucleon is {round(max(n),3)} MeV and occurs at Z = {Z[n.index(max(n))]}"
f" with A = {a[n.index(max(n))]}")
The binomial coefficient is an integer equal to:
def factorial(n):
f = 1.0
for k in range(1,n+1):
f *= k
return f
def binomial(n,k):
a = factorial(n)
b = factorial(k)
c = factorial(n-k)
return a/(b*c)
def pascal(N):
for n in range(0,N+1):
c = n+1
row = list()
for k in range(c):
row.append(int(binomial(n,k)))
print(row)
row.clear()
N = 12 # rows
print(f"The first {N} rows of Pascal's Triangle: \n")
pascal(N)
The probability that an unbiased coin, tossed n times, will come up heads k times is $\binom{n}{k}/2^{n}$.
def toss(n,k):
return binomial(n,k)/2**n
n = 100
k = 60
print(f"The probability that a coin tossed {n} times comes up heads exactly {k} times is p = {round(toss(n,k),5)}")
q = 0
for m in range(0,k):
q += toss(n,m)
p = 1-q
print(f"The probability that a coin tossed {n} times comes up heads {k} or more times is p = {round(p,5)}")
We can develop a fast program for prime numbers by making use of the following observations:
primes = [2]
N = 1000
for n in range(primes[0]+1,N+1):
i=0
boo = True
while primes[i] <= sqrt(n):
if (n%primes[i]) == 0:
boo = False
break
i += 1
if boo == True: primes.append(n)
print(primes)
A useful feature of user-defined functions is recursion, the ability of a function to call itself. Recall the definiton of Catalan numbers in Exercise 2.7.
def catalan(n):
if n == 0:
return 1
else:
return ((4*n+2)/(n+2))*catalan(n-1)
n = 100
print(f"the {n}th Catalan number is {catalan(n)}")
Euclid showed that the greatest common divisor $g(m,n)$ of two nonnegative integers m and n satisfies
def gcd(m,n):
if n == 0:
return m
else:
return gcd(n, m%n)
m = 108
n = 192
print(f"The greatest common divisor of {m} and {n} is {gcd(m,n)}")