\documentclass[a4paper,dvips,12pt]{article}
\usepackage[latin1]{inputenc}
\usepackage[english]{babel}
\usepackage[T1]{fontenc}
\usepackage{amsmath,amssymb}
\begin{document}

\section*{Fjeder Kanon Simulation}
Last week's øvelse resulted in a program (at least I hope it did!) which throws a ball with a given (by the user) initial velocity and with a constant acceleration. Many of you will realize that this is very similar to the 'fjeder kanon' we experimented with a couple of weeks ago. Today, we will try to make our simulation more realistic. Start by making a new directory and copy your program from last week into this directory. First, instead of pretending that our particle (question: why do I call this a particle? what role do the mass and the size as defined in your program play? And, yes, I do expect to see a discussion of these questions in your e-log, so open it now and start writing!) is massless and subject to some mysterious acceleration, let us calculate the net force on the object, and then calculate the acceleration from that. (nb: this is probably overkill for this problem as we are approaching it, but it is useful anyway as it will more easily allow you to modify the program later to include, say, a drag term).

In principle, the force could be time dependent, so let us calculate the force inside our loop over time: 
\begin{equation*}
\text{\texttt{\small{NetForce=mass1*g}}}
\end{equation*}
and note that if there were other forces acting on the mass than gravity, we would include them here. Note that since {\texttt{\small{g}}} has been defined as a vector, the variable {\texttt{\small{NetForce}}} is also now defined as a vector. We now calculate the net acceleration:
\begin{equation*}
\text{\texttt{\small{netA=NetForce/mass1 }}}
\end{equation*}
and we admit to ourselves that this is a very inefficient way to do this for a constant force, but that if there were a time dependent force we would need to do this (questions: when might we need to use time dependent forces? If we had more than one object, would we need to approach this differently? How?). We must now also slightly modify the equations of motion to take account of our new definition of acceleration.
\\ \\
Having done this, we run the program and see that it works exactly as before. Now, modify the program to reflect the measurements you made two weeks ago. That is to say, put your measured values in your program where appropriate. Once you have done this, compare the result of your program to your experimental results. What do you find?
\\ \\
By now, you have realized that one obtains a better estimate of, for example, the mean value of a measured quantity by taking more and more measurements (Is this also true of the results of your program? Why or why not?). One way of answering this question is to try to run the program many times and see how things change (how would you expect the answer you get to change if you run your program 10 times? 100 times? 1000000 times?).

Actually, we don't have to run the program many times, we can make another loop which will do this for us. Now, let's say we want to simulate 1000 shots of our kanon, how should we approach this? One way would be to define a new ``counter'' (like {\texttt{\small{t}}} is the counter in our time loop) we can call, for example, {\texttt{\small{nShots}}}, which will go from {\texttt{\small{0}}} to {\texttt{\small{maxShots}}}: 
\begin{align*}
& \text{\texttt{\small{maxShots=1000}}} \\
& \text{\texttt{\small{nShots=0}}} \\
& \text{\texttt{\small{while nShots < maxShots:}}} \\
& \text{\texttt{\small{  $\#$here you should input everything that changes from shot to shot}}} \\
& \text{\texttt{\small{  while mass1.pos.y>0 $\#$nb: this could be different in your program!}}} \\
& \text{\texttt{\small{    mass1.pos=...}}} \\
& \text{\texttt{\small{    mass1.vel=...}}} \\
& \text{\texttt{\small{    t=t+dt}}} \\
& \text{\texttt{\small{  nShots=nShots+1}}} 
\end{align*}
Now you can quickly loop over many, many shots. This gives us another problem, though. How do we keep track of the results of each shot? And which results, exactly, are we interested in?

Well, we are interested in comparing what we calculate to what we measured. We measured (or should have!) the distance from the canon to where the cannon ball landed. We can do that in our simulation, as well. The distance to the point of impact is just the position of {\texttt{\small{mass1}}} when it "hits" the ``ground''. (How do we define the ``ground''? How do we determine the position at this time?)

One way we can do this (and it is probably neither the only way nor the best way) is to simply get the value of the position when the time loop stops. We can, for example, use an {\texttt{\small{while...else}}} structure to do this (are there better ways?) as follows: 
\begin{align*}
& \text{\texttt{\small{while mass1.pos.y>0 $\#$nb: this could be different in your program!}}} \\
& \text{\texttt{\small{  mass1.pos=...}}} \\
& \text{\texttt{\small{  mass1.vel=...}}} \\
& \text{\texttt{\small{  t=t+dt}}} \\
& \text{\texttt{\small{else: $\#$that is, what if mass1.pos.y is not > 0?}}} \\
& \text{\texttt{\small{print nShots, mass1.pos }}} 
\end{align*}
This will, for each loop over shots, print the value of {\texttt{\small{nShots}}} and the position vector. Clearly, this is of limited utility, and we'd like to somehow store this information in a more useful way than to simply write it to the screen. One way of doing this (again, not necessarily the best way) is to create an array in which to store the values we want to keep. We'll need to write the final value of the position for each of our trials, so we'll need to define the array before the loop over ``shots'' and then add the value for the {\texttt{\small{n}}'th shot to the array after shot {\texttt{\small{n}}}. This will look something like this (although there is probably another, better way to do this, this way should work):
\begin{align*}
& \text{\texttt{\small{xdist=[] $\#$create a 1-d array for x-distance}}} \\
& \text{\texttt{\small{zdist=[] $\#$create a 1-d array for z-distance}}} 
\end{align*}
To write values into the list, one should use: 
\begin{align*}
& \text{\texttt{\small{xdist.append(mass1.pos.x) $\#$write value of x component}}} \\
& \text{\texttt{\small{zdist.append(mass1.pos.z) $\#$write value of z component}}} 
\end{align*}
It is up to you to figure out where these lines of code should be written to get the value in which you are interested (does it matter where in the program these go? what happens if you update the array inside the time loop? what happens if you create the array inside the loop over shots?)
\\ \\
Now that you have stored the data you generate using your simulation, you can use it to calculate, for example, the mean (middelværdi displacement in the plane in which the ball lands, the RMS (spredning), and the undertainty of the mean. Doing so, you will find that all of your statistics are wasted. Why? What can we do about this?
\\ \\
As you no doubt now realize, we need to somehow make our program mimic the uncertianty in our experiment. In our experiment, the ball leaves the canon with a different initial velocity vector each time we fire the canon (not too bad with a little spring canon, but something to be avoided in real canons!). One way we could try to model this is to try to understand in detail all the interactions between the spring, the canon barrel, and the projectile. In principle, we could do this, but it might be somewhat beyond the scope of a 3 hour lab, even for students as fantastic as you. So, we will do something else. We will assume a functional form for the uncertianty in the velocity and use that to mimic all the effects (what might these be?) which go into the uncertainty. That is to say, we sum our ignorance of the details of the system by adding an uncertianty to the result. You will also see that this uncertianty on the velocity is (or should be!) what you have calculated from your experimental results. To include this in your simulation, you need to do a couple of things. The first is that you need to include the {\texttt{\small{random}}} module in your code in the same way you included the {\texttt{\small{visual}}} module earlier. Having done this, you can add the following two lines of code to initialize the velocity and its uncertainty:
\begin{align*}
& \text{\texttt{\small{velocity=vector(vx,vy,vz)      $\#$nb: v is what you have measured! }}} \\
& \text{\texttt{\small{dvx=gauss(meanX,sigmaX)        $\#$what should meanX and sigmaX be? }}} \\
& \text{\texttt{\small{dvy=gauss(meanY,sigmaY)        $\#$what should meanY and sigmaY be? }}} \\
& \text{\texttt{\small{dvz=gauss(meanZ,sigmaZ)        $\#$what should meanX and sigmaX be? }}} \\
& \text{\texttt{\small{deltaV=vector(dvx,dvy,dvz)     $\#$uncertianty on the velocity }}} \\
& \text{\texttt{\small{mass1.velocity=velocity+deltaV $\#$intialize the velocity of a shot }}} 
\end{align*}
The function {\texttt{\small{gauss}}} comes from the module random and draws a random number from a Gaussian distribution:
\begin{equation*}
\text{\texttt{\small{http://en.wikipedia.org/wiki/Gaussian$\_$function}}}
\end{equation*}
It takes two arguments, the mean of the distribution and the variance. It is up to you to figure out where in your program this code should be placed. You will also note that this code can be made more compact; I have chosen this way to make explicity what it is we are trying to do.

Now we have a list of ranges in the two dimensions measured (nb: the height should be the same, or should it?), so we can perform the same mathematical operations on them as we did for the experimental data. Here is one possible way of doing this, but look around and see if you can't figure out a better way. 
\begin{align*}
& \text{\texttt{\small{middelX= distXSum/countMax}}} \\
& \text{\texttt{\small{middelZ= distZSum/countMax}}} \\
& \text{\texttt{\small{icount=1}}} \\
& \text{\texttt{\small{devX=0}}} \\
& \text{\texttt{\small{devX2=0}}} \\
& \text{\texttt{\small{while icount < maxShots:}}} \\
& \text{\texttt{\small{  devX=devX+(dxlist[icount]-middelX)}}} \\
& \text{\texttt{\small{  devX2=devX2+(dxlist[icount]-middelX)**2}}} \\
& \text{\texttt{\small{  icount=icount+1}}} \\
& \text{\texttt{\small{devX=devX/maxShots}}} \\
& \text{\texttt{\small{devX2=devX2/maxShots $\#$standard deviation squared (variance)}}} \\
& \text{\texttt{\small{sigma=pow(devX2,0.5) $\#$standard deviation}}} \\
& \text{\texttt{\small{errMeanX=sigma/(pow(maxShots,0.5)) $\#$error on mean}}} 
\end{align*}
You should also do this for the other dimension. How do the calculated values of the mean and the error on the mean compare to your measured values? Are they the same? Why or why not?

Finally, we will play a little bit with histograms and scatter plots. These are ways of visualizing the data in a more appealing and immediately understandable way than by mearly giving a list of numbers or even a mean and the error on the mean. In order to use VPython's histogramming capability, we need to add yet another module, this time we add the module {\texttt{\small{visual.graph}}}. This will allow us to make histograms (among other things) of our data. We will make two different plots (and you can look into the documentation to see what the following commands do, after all, you are all now experienced programmers and don't need to be shown every single thing anymore. right? RIGHT!?). Here is the code to initialize the graphs: 
\begin{align*}
& \text{\texttt{\scriptsize{graphwindow=gdisplay(xtitle='range x',ytitle='N',ymax=maxShots/4)}}} \\
& \text{\texttt{\scriptsize{distanceX=ghistogram(bins=arange(10,40,1),color=color.red,accumulate=1,average=1)}}} \\
& \text{\texttt{\scriptsize{graphwindow2=gdisplay(xtitle='range X',ytitle='range Z',xmax=50,ymax=50)}}} \\
& \text{\texttt{\scriptsize{hits=gdots(gdisplay=graphwindow2, color=color.green)}}} 
\end{align*}
To update these plots ``event-by-event'' you can use the following code (and you should be able to figure out where it goes, too!):
\begin{equation*}
 \text{\texttt{\small{distanceX.plot(data=xdist) hits.plot(pos=(mass1.pos.x,mass1.pos.z))}}}
\end{equation*}
\\ \\
Next, I'd like you to make a corresponding histogram for the z-direction (unless you measured r and theta, then use those two variables) and compare the histograms resulting from your simulation with that you produced from your data.

As a final, kind of fun thing to do, we can add air resistance (drag) to our simulation and see how the results compare with the ideal case that we have just considered. I will not discuss the model I will use for drag, but you should be able to find this. And since I know that only about 5$\%$ of you attend Lektie Cafe, I know you have some extra time to do this. So, I will assume that I have a resistive force acting on the canon ball, which I call the drag. This force acts to retard the motion, and so must be in the oposite direction of the velocity. Further, I will assume that the magnitude of the force is proportional to the velocity (when is this assumption valid? Is it a reasonable assumption in this case? Would it be a reasonable assumption if we were modelling a professional cyclist? An automobile? A bullet? Our ball falling through a very viscous oil?) It seems reasonable to assume (and I shall) that this drag term should also be proportional to the cross sectional area of the projectile (is this really reasonable, or am I just being sloppy?) and the density of the air. So, I define the constant part of this force as follows: 
\begin{align*}
& \text{\texttt{\small{A=3.14*mass1.radius**2 $\#$cross sectional area of canonball}}} \\
& \text{\texttt{\small{rho=1.29 $\#$density of air}}} \\
& \text{\texttt{\small{C=0.5 $\#$constant you can find the meaning of yourself}}} \\
& \text{\texttt{\small{D=(rho*C*A)/2}}} 
\end{align*}
All that remains now is to incorporate (properly!) this term into the equations of motion for the cannon ball. But now, you are all so good at this that you can figure it out more quickly than I can write it down.

Next week: HotWheels!
\\ \\
Ian Bearden 


\end{document}