from numpy import *
%pylab inline
yocto = 1e-24
mole = 6.022e23
yocto*mole
There are less than one atom in a yoctomole, thus zero.
nano = 1e-9
secs_per_year = 60*60*24*365
years_per_cent = 100
secs_per_year * years_per_cent * nano
Looks like $\pi$!
th_cd = .001 #thickness of cd, meters
bytes_per_cd = 737e6 #number of bytes on a cd
peta = 1e15
h = th_cd * peta / bytes_per_cd
meters_per_khalifa = 829.8
print "%.2f meters, %.2f times the Burj Khalifa"%(h,h/meters_per_khalifa)
There are $10^{80}$ atoms in the universe. The decimal value of the base-2 number which is a series of $n$ ones is $2^n$, so the number written by the atoms is $2^{10^{80}}$.
We can estimate the number of digits in this number using the relation $2^{10k} \approx 10^{3k}$. In this case, $k=10^{79}$, so there are about $3\times 10^{79} +1$ digits.
If we care more precisely, we can use the real formula for the number of decimal digits in $2^n$ of $1+\left \lfloor{n \frac{\log 2}{\log 10} }\right \rfloor $
#we can compare the three ways: explicit, formula, estimate
for i in range(7):
print len("%d"%2**(10**i)), 1+int( (10**i)*log(2)/log(10) ), int(3*10**(i-1)+1)
G = 6.67e-11 #m^3/kg/s/s
M_earth = 5.98e24 #kg
r_earth = 6.378e6 #m
def grav_accel(m,r):
return G*m/(r**2)
def db(s1,s2):
return 20*log(s1/s2)/log(10)
g = grav_accel(M_earth, r_earth)
g_kg = grav_accel(1.,1.)
print 'Earth accel at surface: %f m/s^2 \n One kilogram at one meter: %e m/s^2 \n Decibels of the ratio %f'%(g, g_kg, db(g,g_kg))
Approximately estimate the chemical energy in a ton of TNT. You can assume that nitrogen is the primary component; think about what kind of energy is released in a chemical reaction, where it is stored, and how much there is.
kg_per_ton = 907. #kg
energy_per_bond = 1 #eV
j_per_ev = 1.6e-19
atomic_weight_nitrogen = 14.
mole = 6e23
atoms_nitrogen_per_g = mole/atomic_weight_nitrogen
bonds_per_atom = 1
def J_per_kg_TNT():
return 1e3*atoms_nitrogen_per_g*bonds_per_atom*energy_per_bond*j_per_ev
def J_per_ton_TNT():
return J_per_kg_TNT()*kg_per_ton
print "%.1f MJ per kg of TNT"%(J_per_kg_TNT()*1e-6)
print "%.1f GJ per ton of TNT"%(J_per_ton_TNT()*1e-9)
This is roughly the right order of magnitude for this number. http://physics.nist.gov/Pubs/SP811/appenB8.html
atomic_weight_uranium = 238.
atoms_uranium_per_g = mole/atomic_weight_nitrogen
energy_per_nuclear_bond = 1e6 #eV
bonds_per_atom = 1
def J_per_kg_uranium():
return 1e3*atoms_uranium_per_g*bonds_per_atom*energy_per_nuclear_bond*j_per_ev
def J_per_ton_uranium():
return J_per_kg_uranium()*kg_per_ton
T = 1e4 #tons of TNT
print "%d tons of TNT has the same energy as %.1f kg of uranium"%(T, T*J_per_ton_TNT()/J_per_kg_uranium() )
c=3e8 #m/s
def rest_mass_energy(m):
return m*c*c
m=9.1
print "The rest mass energy of %.1f kg of uranium is %.2f EJ, compared to our estimate of %.2f EJ for the nuclear energy released."%(m, rest_mass_energy(m)*1e-12, m*J_per_kg_uranium()*1e-12)
h = 6.6e-34 #Js, plank constant
m_baseball = .15 #kg, mass of baseball
v_baseball = 45 #m/s, velocity of baseball
def debroglie_wavelength(m,v):
return h/m/v
print "Debroglie wavelength of a baseball is %e meters."%(debroglie_wavelength(m_baseball,v_baseball))
This is on the order of the Planck length!
T_room = 300. #Kelvin, room temp
P_room = 101000 #Pa, atmostpheric pressure
B = 1.38e-34 #J/K, Boltzmann constant
mass_N2 = 2*14*1e3/mole
def v_N2(T):
#velocity of N2 particle
return sqrt(2*T*B/mass_N2)
print "Debroglie wavelength of a N2 atom at room temp and pressure is %.2e meters."%(debroglie_wavelength(mass_N2,v_N2(T_room)))
rho_N2 = 1.2 #kg/m^3, density at room temperature and pressure
N2_per_m3 = rho_N2/mass_N2
side_length_per_N2 = (1./N2_per_m3)**.3333
#assume a cubic volume for each molecule
print "Rough distance between molecules is %.2e"%(side_length_per_N2)
T_d = linspace(.1,100,100)
plot(T_d,debroglie_wavelength(mass_N2,v_N2(T_d)))
xlabel("Temperature (K)")
ylabel("DeBroglie Wavelength (m)")
T_d = 3e-1
print "%.2f degrees K gives a DeBroglie wavelength of %.2e meters."%(T_d,debroglie_wavelength(mass_N2,v_N2(T_d)))
def U_G(m,M,r):
#gravitational potential energy
return -G*M*m/r
def escape_velocity(M,r):
return sqrt(2*G*M/r)
The radius at which nothing (not even light) can escape is $\frac{2GM}{c^2}$.
If the rest energy of a mass $M$ is converted into a photon, the wavelength is $\frac{h}{Mc}$, where $h$ is Planck's constant.
The mass at which the two are equal is $\sqrt{\frac{hc}{2G}}\approx 3.85\times10^{-8} kg$.
The corresponding size is $\frac{1}{c} \sqrt{\frac{2Gh}{c}} \approx 5.71\times10^{-35} m$. Roughly Planck distance
The corresponding energy is $c^2\sqrt{\frac{hc}{2G}} \approx 3.46 \times10^9 J \approx 2\times10^{28} eV$. Roughly the Planck energy.
The corresponding period is $\frac{\lambda}{c} = \frac{1}{c^2} \sqrt{\frac{2Gh}{c}} \approx 1.90\times10^{-43}$. Roughly Planck time.
Pyramid problem¶
Pyramid of height $H$ and square base of side length $2L$. A sphere is centered at the center of base, tangent to all edges of the pyramid.
What is $H$?
We can consider the triangle made of of the diagonal of the base, with the corresponding edges leading to the apex. The base length of this triangle is $2L\sqrt{2}$. The circle is tangent to the edges leading to the apex (with radius $L$), so $H=L\sqrt{2}$.
What is the common volume to the pyramid and sphere?
To tackle this problem, we consider cross sections normal to the height of the pyramid. In this cross section, we will have a circular section of the sphere and a square section of the pyramid. If we integrate the common area of these shapes over the range of $z$ values, we will get the desired volume. We define functions $r(z)$ and $s(z)$ giving the radius of the circular cross section of the sphere and the side length of the square cross section of the pyramid respectively. We calculate $r(z) = \sqrt{L^2-z^2}$ and $s(z) = 2L-\sqrt{2}z$. To get a better idea of the geometry, we plot these shape trajectories below:
L=1.
def r(z):
#radius of circular section of sphere
return sqrt(L**2-z**2)
def s(z):
#side length of pyramidal section
#return 2*L*(1-z/(L*sqrt(2)))
return 2*L-sqrt(2)*z
def plot_pyramid(z,artists=None):
#ax=subplot(aspect='equal')
if artists:
circ = artists[0]
circ.radius = r(z)
rect = artists[1]
rect.set_xy((-.5*s(z),-.5*s(z)))
rect.set_height(s(z))
rect.set_width(s(z))
else:
circ=Circle((0,0),r(z),color='r',alpha=0.1)
rect=Rectangle((-.5*s(z),-.5*s(z)),s(z),s(z),color='b',alpha=0.1)
return circ,rect
#need to evaluate video embedding cell at bottom first!
from matplotlib import animation
N_frames = 200
fig = figure(figsize=(8,8))
circ,rect = plot_pyramid(0.)
fig.gca().add_artist(circ)
fig.gca().add_artist(rect)
fig.gca().set_aspect('equal')
ylim([-L,L])
xlim([-L,L]);
def init():
return plot_pyramid(0.)
def animate(i):
return plot_pyramid(L*i/N_frames,artists=(circ,rect))
anim = animation.FuncAnimation(fig, animate,
init_func=init,
frames=N_frames,
interval=5,
blit=True)
display_animation(anim)
We see there are two distinct cases of cross section geometry: before the circle is completely contained in the square and after it is. To calculate the z value for the crossing point between these two phases, we solve $s(z) = 2r(z)$ to get $z = \frac{2L\sqrt{2}}{3}$. We can then split the desired integral into two regions separated by this value of $z$.
We need now to calculate the shared area. After the critical $z$ value, this is simple the area of the circular cross section: $\pi r(z)^2$. Before the critical value, the expression for shared area is more complicated.
$\int \frac{1}{2}(2 L-\sqrt{2} z) \sqrt{z (\sqrt{2} L-z/2)} dz = \frac{1}{6} (-z (z-2 \sqrt{2} L))^{3/2}$
import sympy
z = sympy.Symbol('z', real=True)
L = sympy.Symbol('L', real=True, positive=True)
def r(z): return sympy.sqrt(L*L-z*z)
def s(z): return 2*L-sympy.sqrt(2)*z
def alpha(z): return sympy.acos(s(z)/(2*r(z)))
def beta(z): return sympy.pi/2-2*alpha(z)
def area_tri(z): return (2*L-sympy.sqrt(2)*z)/2*sympy.sqrt(z*(sympy.sqrt(2)*L-z/2))
def area_sector(z): return r(z)*r(z)*beta(z)/2
sympy.integrate(area_tri(z),(z,0,2*L*sympy.sqrt(2)/3))
So the area of the triangular region integrated over the region of interest is $\frac{32 L^3}{81}$.
sympy.integrate(area_sector(z),(z,0,2*L*sympy.sqrt(2)/3))
sympy.integrate( sympy.pi*r(z)**2, (z, 2*L*sympy.sqrt(2)/3, L))
So the area of the top cap of the sphere is $\frac{2\pi}{3}L^3 - \frac{38}{81}\pi L^3\sqrt{2}$.
from tempfile import NamedTemporaryFile
#import base64
VIDEO_TAG = """<video controls>
<source src="data:video/x-m4v;base64,{0}" type="video/mp4">
Your browser does not support the video tag.
</video>"""
def anim_to_html(anim):
if not hasattr(anim, '_encoded_video'):
with NamedTemporaryFile(suffix='.mp4') as f:
anim.save(f.name, writer='mencoder', fps=20) #extra_args=['-vcodec', 'libx264'])
video = open(f.name, "rb").read()
anim._encoded_video = video.encode("base64")
return VIDEO_TAG.format(anim._encoded_video)
from IPython.display import HTML
def display_animation(anim):
plt.close(anim._fig)
return HTML(anim_to_html(anim))
circ,rect = plot_pyramid(0.15)
fig = figure()
fig.gca().add_artist(circ)
fig.gca().add_artist(rect)
fig.gca().set_aspect('equal')
ylim([-L,L])
xlim([-L,L]);