Ph265: Take-home final exam
(deadline: at the end of the document)

Problem 1. Modify the program birthdays.py we used in class so that it:
  1. Calculates the approximate probablity that three individuals in a randomly selected group of N people have their birthday on the very same day;
  2. Calculates the approximate probability that in a group of N randomly selected people there are two pairs who celebrate their birthdays on the same day (example: Mrs. Smith and Mr. Brown both on June 4, and Ms. White and Mr. Jackson both on Nov. 28).
Problem 2. Write a program simulating an individual act of rolling a die. You may use the random.random() function that generates random numbers from the (0,1) range. If you multiply it by an appropriate factor and convert the result to an integer using the int(decarg) function (where "decarg" is a decimal argument), then you may get a sequence of random integers from 1 to 6.
  1. First, check whether your program indeed generates these numbers with a uniform probability. For instance, simulate 600 000 "dice rolling acts", add all "ones", all "twos", etc., each to a separate variable (you will have to use "if" many times!), and print the results to see whether all the six sums agree within a reasonable "statistical error margin" (mathematical statistics theory says that in such counting the standard deviation of a sum is given simply by its square root; so, for 100 000 events the fluctuations should be plus minus 316, so that two numbers should not differ much more than twice this value; if the differences between the numbers you obtain are much larger than, say, 1000, it would mean that something is wrong.

  2. Next, when you are already sure that your simulation of dice rolling works OK, write the next part of the program. The problem is the following: you toss N dice and count all "dots" -- what is the probability that the sum equals n? It does not matter whether you have N dice an toss them all ("parallel rolling"), or have only one die and toss is N times ("sequential rolling") -- the probability of obtaining a given n value is the same in both these methods. But sequential rolling, in my opinion, is easier to think of.

    The value of n cannot be lower than N, or larger than 6N. So for all n values outside this range the probability is zero. For n = N or n = 6N the answer is not too difficult, either, because the total number of all sequences one can obtain is 6N, and there is only one way of obtaining n = N (all ones), and n = 6N (all sixes). So, the probability for theses two n-s is simply p(N) = p(6N) = 1/6N. For n = N+1 the result can still be readily obtained, because there are only N such sequences possible (211111..., 1211111..., 1121111..., etc.), so the probability is now p(N+1) = N/6N, and the same is true for p(6N-1). But if you want to continue in this way, then the problem becomes a real nightmare!

    Monte Carlo method is a viable tool for obtaining pretty accurate values of the probability for n values that are not too far from the most probable value, i.e., n = 7N/2. You simply simulate rolling a die N times, add all results, and check whether the sum is equal to a given n. You repeat such a proces M times, where M is a reasonably large number, and count all instances in which you get a positive result. Then, simply divide that number by M, and you get your probability.

    So, I want you to write a program that prompts for N, and for n, and calculates the probability p(n). Use a "reasonable" M value -- meaning large one, but not such large that you would have to wait forever for the result.

Problem 3. Write a program that calculates the motion of a simple pendulum consisting of a point mass suspended on a weightless string of lenght L. When the pendulum swings, the angle of its maximum displacement from its equilibrium position is thetamax.

The program should prompt for the L and thetamax values (the latter in degrees), and then it should calculate the pendulum's motion (using either the Euler's or the Verlet's algorithm, whichever you like better). The results should be displayed in a graphic form, as a theta(t) function, where theta is the displacement angle, and t is time (the program should also prompt for the total time period to be graphed).

It is well known that a simple pendulum can be thought of as a harmonic oscillator of angular frequency omega = sqrt(g/L). However, thisis a good approximation only for small swinging amplitudes (i.e., for small thetamax valuse). What happens for larger amplitudes? Does the oscillation period T become longer, or shorter than the "ideal" value of 2pi/omega? I want you to find the answer using a simple trick: namely, on the same field you graph the calculated theta(t) function, plot an "idealized" function thetamax*cos(omega*t). For small amplitudes the two functions should be identical -- and you will immediately see "what's going on" for larger thetamax values!

In this website you may find the results of the exact analythical theory of a simple pendulum. You may find in it how much the period T differs from the "ideal" period for a few "larger" amplitudes -- and you may check, if you wish, whether your program results are consistent with those data.

Hints: In the Euler's or Verlat's algorithm, my advice is to use not the angular coordinate theta (and definitely not the xy Carthesian coordinates!), but a linear coordinate l (small L -- sorry, it does not show up very well on the screen)), defined as the displacement from the equilibrium position measured along the circular path of the pendulum's end. The acceleration one should use in such calculations is the g vector component tangential to this circular path -- which can be readily calculated using just one trigonometric function if you know the value of the displacement angle theta, right? But note that the relation between your l (again, "small L") coordinate and theta is very simple: theta = l/L. Of course, if theta is expressed in radians, not in degrees!

I asked you that the input value of thetamx should be in degrees -- so you need to convert it to radians right away, and use radians throughout the program. But the theta(t) values in the graph should be again in degrees. Which is definitely not a big problem, as converting degrees to radians and radians back to degrees is quite an easy task.

Problem (or Task) 4. In class, we briefly discussed the Monte Carlo program for calculating the critical mass of uranium. The program prompts for the radius of the 235U sphere, and the "compression coefficient". When you enter 1, the program makes calculations for uranium of normal density. If 2, for "compressed" uranium of density twice as high as that of normal density uranium.

After you enter the above two values, the calculates the mass of the sphere. Then, it simulates the path of a neutron created by a fission, and displays the results in a graphical form. The green dot shows where the neutron is "born", and the red line is the calculated path, "zig-zacking" due to scattering processes. A neutron colliding with with a uranium nucleus may either get scattered, or may trigger a fission (releasing 2 or 3 "next generation" neutrons). Such a process is indicated by an orange dot, and the two or three orange lines branching out from it symbolize the beginning of the paths of the "next generation" neutrons.

The probability that a neutron gets scattered is about 2.5 times larger than the probability that it triggers a fission. Therefore, most neutrons, after some "zig-zacking", reach the surfaca and escape form the sphere. The "escape point" is shown by a blue dot, and the blue line is the neutron's flight path after the escape.

However, even though most neutrons escape, those who trigger fission events release such a number of neutrons that their overall population inside may increase from generation to generation. This is the infamous "chain reaction".

The program not only simulates and displays the single-neutron paths, but it also may simulate the course of events when a larger number of neutrons are simultaneously "born" inside the sphere. If you are fed up with wathing the "single-neutron paths", enter "n", and then the other part will start. The initial "portion" of neutrons in the sphere is 1000. For all of them the program performs an analogous simulations as before, but now it "memorizes" the location of every fission event that has occurred, and how many "next generation" neutrons were released in each. Then, it starts pocessing those "next generation" neutrons in the same way. And so on, for as many generations as you wish (the program prompts for that). The number of neutrons created in every generation is plotted in the graph that appears on the screen.

For uncompressed uranium, as you may check, the "critical" radius of the sphere is about 8.73 cm. For smaller values, the number of neutrons in consecutive generations shows a clearly decreasing tendency -- but above this value a real chain reaction starts. If the value is much higher, the increase is so fast that the program may crash -- it can handle only 10 000 neutrons.

The simulation is quite realistic -- well, Dr. Tom has spent about 40 years of his professional cereer studying neutron scattering processes, so he knows the underlying physics pretty well. The "critical radius" value of R = 8.73 cm that the program yields corresponds to a mass of about 53 kilograms, which appears to be pretty close to the critical mass value of Uranium 235 that one can find in literature (many such sources can be readily found in the Web).

But it turns out that if the sphere is "squeezed", the critical mass changes. You can easily check that if you enter a "compression factor" larger than one. Start, e.g., with 1.2, and see what happens. You will find that the "critical radius" is markedly smaller. And, of course, so is the corresponding mass.

So, what I want you to is to find out the dependence of the critical mass on the sphere density. Enter a compression factor 1.4142 (square root of two) and by the trial-and error method determine the approximate value of the critical mass. The do the same with 1.732 (square root of 3) and then with 2. And it will probably be enough for figuring out what the dependence in question is!

I hope that you will get the answer in less than 30 minutes. If not, don't spend more time on it!

Problem (or Task) 5. I want you to play a while with the "docking" program. I prepared a version of the program -- click on this to see it -- that records all relevant inputs and outputs in a file named "docking.log".

Please run the program. Let the initial distance between the spacecraft be 1 second. Then try to get the spacecrafts as close to each other as you can, by using short "bursts of thrust", and then letting the spacecrafts to mo mowe freely for some time. The thrust is such that it produces an acceleration of 10 m/s2, and its direction may be chosen by inputing the angle:
  1. zero degrees is a burst sent "backward", adding speed;
  2. 90 degrees is a burst sent "downward", toward Earth, pushing the spacecraft up;
  3. 180 degrees is a "reverse thrust", decelerating the spacecraft;
  4. and 270 degrees is a burst sent "up", pushing the spacecraft down.
Usually people who start using the program miserably fail at the beginning, but then they quickly learn how to maneuver the crafts efficiently. If A-B distance becomes large -- larger than the initial distance -- then my advice is not to try desperately to "resque" the situation, but to terminate the run, and start another one.

My best result was some 70 meters -- can you get a better result? Please send me the "docking.log" file recording your most successful trial, as the "solution" of the present Problem.

Please keep in mind that each time you start a new run, the previous "docking.log" file is erased. So, if you get a result that is really good, but you want to try to get a result that is even better -- first rename the "docking.log" file, and only then re-run the program, because you may never be sure that you'll get the "even better" result!
Because the exam was posted later than intended, the dealine is extended to Wednesday (by the end of the day).