Introduction to Python

Taken from Computational Physics by Newman Chapter 2

1. Another ball dropped from a tower

A ball is dropped from a tower of height h with initial velocity zero. Write a program that asks the user to enter the height in meters of the tower and then calculates and prints the time the ball takes until it hits the ground, ignoring air resistance. Use your program to calculate the time for a ball dropped from a 100m high tower.

2. Altitude of a satellite

A satellite is to be launched into a circular orbit around the Earth so that it orbits the planet once every T seconds.

Show that the altitude h above the Earth’s surface that the satellite must have is

h = \biggl( {GMT^2\over4\pi^2} \biggr)^{1/3} - R

where G=6.67\times10^{-11}\,\textrm{m}^3\,\textrm{kg}^{-1}\,\textrm{s}^{-2} is Newton’s gravitational constant, M=5.97\times10^{24}\,kg is the mass of the Earth, and R=6371\,km is its radius.

Write a program that asks the user to enter the desired value of T and then calculates and prints out the correct altitude in meters.

Use your program to calculate the altitudes of satellites that orbit the Earth once a day (so-called “geosynchronous” orbit), once every 90 minutes, and once every 45 minutes. What do you conclude from the last of these calculations?

Technically a geosynchronous satellite is one that orbits the Earth once per sidereal day, which is 23.93 hours, not 24 hours. Why is this? And how much difference will it make to the altitude of the satellite?

3. Coordinate Tranformation

Write a program to perform the inverse operation to that of Example 2.2. That is, ask the user for the Cartesian coordinates x,y of a point in two-dimensional space, and calculate and print the corresponding polar coordinates, with the angle \theta given in degrees.

4. Space Travel

A spaceship travels from Earth in a straight line at relativistic speed v to another planet x light years away. Write a program to ask the user for the value of x and the speed v as a fraction of the speed of light c, then print out the time in years that the spaceship takes to reach its destination (a) in the rest frame of an observer on Earth and (b) as perceived by a passenger on board the ship. Use your program to calculate the answers for a planet 10 light years away with v=0.99c.

5. Quantum potential step

A well-known quantum mechanics problem involves a particle of mass m that encounters a one-dimensional potential step, like this:qstep
The particle with initial kinetic energy E and wavevector 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\”odinger 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 wavevectorof k_2=\sqrt{2m(E-V)}/\hbar, or (b)~it may be reflected, keeping all of its kinetic energy and an unchanged wavevector but moving in the opposite direction. The probabilities T and R for transmission and reflection are given by

T = {4k_1k_2\over(k_1+k_2)^2}
R = \biggl( {k_1-k_2\over k_1+k_2} \biggr)^2

Suppose we have a particle with mass equal to the electron mass m=9.11\times10^{-31}\,kg and energy 10eV encountering a potential step of height 9eV. Write a Python program to compute and print out the transmission and reflection probabilities using the formulas above.

6. Planetary orbits

The orbit in space of one body around another, such as a planet around the Sun, need not be circular. In general it takes the form of an ellipse, with the body sometimes closer in and sometimes further out. If you are given the distance \ell_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.
Kepler’s second law tells us that the distance \ell_2 and velocity v_2 of the planet at its most distant point, or aphelion, satisfy \ell_2 v_2 = \ell_1 v_1. At the same time the total energy, kinetic plus gravitational, of a planet with velocity v and distance r from the Sun is given by E = \frac{1}{2} m v^2 - G {mM\over r} where m is the planet’s mass, M=1.9891\times10^{30}kg is the mass of the Sun, and G=6.6738\times10^{-11} {m^3\,kg^{-1}\,s^{-2}} is Newton’s gravitational constant. Given that energy must be conserved, show that v_2 is the smaller root of the quadratic equation

v_2^2 - {2GM\over{v_1 \ell_1}} v_2 -[v_1^2 -{2GM\over{\ell_1}}] = 0

Once we have v_2 we can calculate \ell_2 using the relation \ell_2 = \ell_1 v_1/v_2.  Given the values of v_1, \ell_1, and \ell_2, other parameters of the orbit are given by simple formulas can that be derived from Kepler’s laws and the fact that the orbit is an ellipse:

Semi-major axis: a = \frac{1}{2}(\ell_1+\ell_2)
Semi-minor axis: b = \sqrt{\ell_1\ell_2}
Orbital period: T= {2\pi ab\over\ell_1 v_1}
Orbital eccentricity: e ={\ell_2-\ell_1\over\ell_2+\ell_1}

Write a program that asks the user to enter the distance to the Sun and velocity at perihelion, then calculates and prints the quantities \ell_2, v_2, T, and e.
Test your program by having it calculate the properties of the orbits of the Earth (for which \ell_1=1.4710\times10^{11}m and v_1=3.0287\times10^4\,\mathrm{m\,s^{-1}}) and Halley’s comet
(\ell_1=8.7830\times10^{10}\,m and v_1=5.4529\times10^4\,\mathrm{m\,s^{-1}}). Among other things, you should find that the orbital period of the Earth is one year and that of Halley’s comet is about 76 years.

10. The semi-empirical mass formula

In nuclear physics, the semi-empirical mass formula is a formula for calculating the approximate nuclear binding energy B of an atomic nucleus with atomic number Z and mass number A:
B = a_1 A - a_2 A^{2/3} - a_3 {Z^2\over A^{1/3}} - a_4 {(A - 2Z)^2\over A} + {a_5\over A^{1/2}}
where, in units of millions of electron volts, the constants are a_1=15.67, a_2=17.23, a_3=0.75, a_4=93.2, and a_5=0 if A is odd, a_5 = 12.0 if A and Z are even or a_5 = -12.0 if A is even and Z is odd.

Write a program that takes as its input the values of A and Z, and prints out the binding energy for the corresponding atom. Use your program to find the binding energy of an atom with A=58 and Z=28. (Hint: The correct answer is around 490MeV.) Modify your program to print out not the total binding energy B, but the binding energy per nucleon, which is B/A.Now modify your program so that it takes as input just a single value of the atomic number Z and then goes through all values of A from A=Z to A=3Z, to find the one that has the largest binding energy per nucleon. This is the most stable nucleus with the given atomic number. Have your program print out the value of A for this most stable nucleus
and the value of the binding energy per nucleon. Modify your program again so that, instead of taking Z as input, it runs through all values of Z from 1 to 100 and prints out the most
stable value of A for each one. At what value of Z does the maximum binding energy per nucleon occur? (The true answer, in real life, is Z=28, which is nickel. You should find that the semi-empirical mass formula gets the answer roughly right, but not exactly.)

11 Binomial coefficients

The binomial coefficient {n\choose k} is an integer equal to
{n\choose k} = {n!\over k!(n-k)!} = {n\times(n-1)\times(n-2)\times\ldots\times(n-k+1)\over 1\times2\times\ldots\times k}
when k\ge1, or {n\choose0}=1 when k=0.
Using this form for the binomial coefficient, write a user-defined function binomial(n,k) that calculates the binomial coefficient for given n and k. Make sure your function returns the answer in the form of an integer (not a float) and gives the correct value of 1 for the
case where k=0.  Using your function write a program to print out the first 20 lines
of “Pascal’s triangle.” The nth line of Pascal’s triangle contains n+1 numbers, which are the coefficients {n\choose 0}, {n\choose1}, and so on up to {n\choose n}. Thus the first few lines are
1 1
1 2 1
1 3 3 1
1 4 6 4 1
The probability that an unbiased coin, tossed n times, will come up heads k times is {n\choose k}/2^n. Write a program to calculate (a) the total probability that a coin tossed 100 times comes up heads exactly 60 times, and (b) the probability that it comes up heads 60 or more times.

12 Prime numbers

The program in Example 2.8 is not a very efficient way of calculating prime numbers: it checks each number to see if it is divisible by any number less than it. We can develop a much faster program for prime numbers by making use of the following observations:
A number n is prime if it has no prime factors less than n. Hence we only need to check if it is divisible by other primes.  If a number n  is non-prime, having a factor r, then n=rs, where s is also a factor. If r\ge\sqrt{n} then n = rs \ge \sqrt{n}s, which implies that s\le\sqrt{n}. In other words, any non-prime must have factors, and hence also prime factors, less than or equal to \sqrt{n}. Thus to determine if a number is prime we have to check its prime factors only up to and including \sqrt{n}—if there are none then the number is prime.  If we find even a single prime factor less than \sqrt{n} then we know that the number is non-prime, and hence there is no need to check any further—we can abandon this number and move on to something else.
Write a Python program that finds all the primes up to ten thousand. Create a list to store the primes, which starts out with just the one prime number 2 in it. Then for each number n from 3 to 10,000 check whether the number is divisible by any of the primes in the list up to and
including \sqrt{n}. As soon as you find a single prime factor you can stop checking the rest of them—you know n  is not a prime. If you find no prime factors \sqrt{n} or less then n is prime and you should add it to the list. You can print out the list all in one go at the end of the program, or you can print out the individual numbers as you find them.

13 Recursion

A useful feature of user-defined functions is recursion, the ability of a function to call itself. For example, consider the following definition of the factorial n! of a positive integer n:
$latex
n! = \biggl\lbrace\begin{array}{ll}
1 & \qquad\mbox{if $n=1$,} \\
n\times(n-1)! & \qquad\mbox{if $n>1$.}
\end{array}
$
This constitutes a complete definition of the factorial which allows us to
calculate the value of n! for any positive integer. We can employ this
definition directly to create a Python function for factorials, like this:

def factorial(n):
    if n==1:
       return 1
    else:
       return n*factorial(n-1)

Note how, if $n$ is not equal to~1, the function calls itself to calculate
the factorial of $n-1$. This is recursion. If we now say
“\verb|print(factorial(5))|” the computer will correctly print the
answer~120.
\begin{enumerate}\setlength{\itemsep}{0pt}
\item We encountered the Catalan numbers~$C_n$ previously in Exercise~2.7
on page~46. With just a little rearrangement, the definition given there
can be rewritten in the form
\begin{displaymath}
C_n = \left\lbrace\begin{array}{ll}
\rule[-9pt]{0pt}{10pt}1 & \qquad\mbox{if $n=0$,} \\
\dfrac{4n-2}{n+1}\,C_{n-1} & \qquad\mbox{if $n>0$.}
\end{array}\right.
\end{displaymath}
Write a Python function, using recursion, that calculates~$C_n$. Use your
function to calculate and print~$C_{100}$.
\item Euclid showed that the greatest common divisor~$g(m,n)$ of two
nonnegative integers $m$ and $n$ satisfies
\begin{displaymath}
g(m,n) = \biggl\lbrace\begin{array}{ll}
m & \qquad\mbox{if $n=0$,} \\
g(n,m\>\textrm{mod}\>n) & \qquad\mbox{if $n>0$.}
\end{array}
\end{displaymath}
Write a Python function \verb|g(m,n)| that employs recursion to calculate
the greatest common divisor of $m$ and $n$ using this formula. Use your
function to calculate and print the greatest common divisor of 108 and~192.
\end{enumerate}
Comparing the calculation of the Catalan numbers in part~(a) above with
that of Exercise~2.7, we see that it’s possible to do the calculation two
ways, either directly or using recursion. In most cases, if a quantity can
be calculated \emph{without} recursion, then it will be faster to do so,
and we normally recommend taking this route if possible. There are some
calculations, however, that are essentially impossible (or at least much
more difficult) without recursion. We will see some examples later in this
book.