# Physics 420/580 Final Exam 
## December 19, 2019 2pm-5pm  L1-150

Do three of the four following problems.  If you attempt all four, your best three will be used for your mark.  Use the Jupyter notebook, inserting your code and any textual answers/explanations in cells between the questions. (Feel free to add additional cells!)  Marks will be given based on how clearly you demonstrate your understanding.    

There are no restrictions on downloading from the internet, eclass,  or the use of books, notes, or any other widely available computing resources.  However, **you are not allowed** to communicate with each other or collaborate in any way and uploading to the internet or sending or receiving direct communications is not appropriate.  You are required to sit the exam in L1-150, and to upload it before you leave.     

When you are finished, upload the jupyter notebook to eclass. 

Also be careful to save the notebook periodically and double check **that you upload the file that has your answers in it!**



```python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import scipy
from scipy.integrate import odeint
from scipy.integrate import quad
from scipy.optimize import minimize
from scipy.optimize import fsolve
from scipy.optimize import least_squares
from scipy.interpolate import interp1d
```


```python
mpl.rc('figure',dpi=250)
mpl.rc('text',usetex=True)
```


```python
def add_labels(xlabel, ylabel, title):
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    plt.title(title)
    plt.legend()
```

# The Normal Distribution

 Write a simple program that:
 
 a.  Calculates and plots the function  
 $$G(x;\mu,\sigma)=1000e^{-\frac{(x-\mu)^2}{2\sigma^2}}$$ with an average value $\mu=10$, width $\sigma=2$.  This is proportional to the normal distribution.


```python
"""
Plan:
-Here we just graph G, and plot it
"""
```


```python
mu = 10
sigma = 2
def G(x, mu, sigma, A):
    return A*np.exp(-(x - mu)**2/(2*sigma**2))
```


```python
x = np.linspace(-25, 25)
plt.plot(x, G(x, mu, sigma, 1000), label=f"mu = {mu}, sigma = {sigma}, A = {1000}")
add_labels('x', 'G(x)', 'G(x)')

```

The three dimensional normal distribution is the product of normal distributions in the three coordinates:
\begin{equation} G(\vec{x}; \vec{x_0},\sigma, A)=Ae^{-\frac{(\vec{x}-\vec{x_0})^2}{2\sigma^2}}\end{equation}

Set $\vec{x_0}=(6,8,0)$ and throw 100000 random points by throwing the three coordinates independently with $\sigma=2$ for each coordinate.

Histogram and plot the distributions for each coordinate.  Comment on how you pick the right number of bins for your histogram.





```python
"""
Plan:
We could use our normal function G, and sample from it: just plug it into generate samples. However, but we realize that in numpy this is already done for us: 
Numpy.random.normal. It's the same thing, just that np.random.normal just much more efficient and I trust it more than myself. It runs faster too.


"""
```




    "\nPlan:\nWe could use our normal function G, and sample from it: just plug it into generate samples. However, but we realize that in numpy this is already done for us: \nNumpy.random.normal. It's the same thing, just that np.random.normal just much more efficient and I trust it more than myself. It runs faster too.\n\n\n"



```
def generate_samples(func, bounds, N): 
    """
    Assignment 5: function for sampling from a custom 1d distribution.
    
    Input:  func - any function
            bounds - tuple (a, b)
    
    Output: 1d samples of size N
    """
    a, b = bounds
    x = np.linspace(a, b)
    
    f_max = max(func(x))
    
    x_rand = np.random.uniform(a, b, size=N)
    y_rand = np.random.uniform(0, f_max, size=N)
        
    samples = x_rand[y_rand < func(x_rand)]

    return samples
```


```python
x0 = np.array([6, 8, 0])
N = 100000
sigma = 2
```


```python
throws = np.zeros((N, 3))

for i in range(N):
    throws[i, 0] = np.random.normal(x0[0], sigma)
    throws[i, 1] = np.random.normal(x0[1], sigma)
    throws[i, 2] = np.random.normal(x0[2], sigma)
    
throws.shape

# I picked the bins based on the domain and through just visual comparison. The bins are 100 linearly spaced along the domain. 
plt.figure()
plt.hist(throws[:, 0], bins=100, label=f'mu = {x0[0]}, sigma = 2')
add_labels('x', 'counts', 'distribution of x coordinate')
plt.figure()
plt.hist(throws[:, 1], bins=100, label=f'mu = {x0[1]}, sigma = 2')
add_labels('y', 'counts', 'distribution of y coordinate')
plt.figure()
plt.hist(throws[:, 2], bins=100, label=f'mu = {x0[2]}, sigma = 2')
add_labels('z', 'counts', 'distribution of z coordinate')
```

Now histogram and plot the distribution of $r=|\vec{x}-\vec{x_0}|$.





```python
radial = np.linalg.norm(throws - x0, axis = 1)
plt.hist(radial, bins=100)
add_labels('r', 'counts', 'distribution of radial coordinate with respect to mean')
```

This should look like the analytic 3d normal distribution in spherical coordinates:
\begin{equation}G(r;\sigma,A)=4\pi A r^2 e^{-\frac{r^2}{2\sigma^2}}.\end{equation}  We derived this equation by converting $G(\vec{x}; \vec{x_0},\sigma, A)$ into spherical coordinates and integrating over $\theta$ and $\phi.$

Determine A for 100000 points by numerically integrating:
\begin{equation}100000=\int_0^{20} G(r;\sigma,A)dr =A\int_0^{20} G(r;\sigma,1)dr\end{equation}

With a known A, you can now replot your histogram and overlay the calculation of $G(r;\sigma,A)$.  To have $G(r;\sigma,A)$ overlay properly you will need to include the bin widths, which you can do either by the exact integral or the approximation: $\int_{r_i}^{r_{i+1}}G(r;\sigma,A)dr \approx G(r;\sigma,A)\Delta r$, where $r_i$ is the start of the i-th bin and $r_{i+1}$ is the start of the next bin, and $\Delta r=r_{i+1}-r_i$.  




```python
def radial_Gaussian(r, sigma):
    return 4*np.pi*r**2*np.exp(-r**2/(2*sigma**2))
```


```python
sigma = 2
integral = quad(radial_Gaussian, 0, 20, args=sigma)

A = 100000/integral[0]
print(f'A = {A} is the normalization factor')
```

    A = 793.6704491780121 is the normalization factor



```python
r = np.linspace(0, 10)
plt.plot(r, A*radial_Gaussian(r, sigma))
```


```python
plt.hist(radial, bins=100)
add_labels('r', 'counts', 'distribution of radial coordinate with respect to mean')
plt.plot(r, A*radial_Gaussian(r, sigma))
```


```python
"""I have to normalize it now, and I have to fudge with some bins but if I have time after the exam I will do that. Right now, I need to solve the other 3 questions"""
```




    'I have to normalize it now, and I have to fudge with some bins but if I have time after the exam I will do that. Right now, I need to solve the other 3 questions'



# Monte Carlo

A radioactive source is placed into a large water Cerenkov detector (which is able to detect high energy electrons/gammas in pure water by looking at visible photons from the Cerenkov process.)  Before the source is deployed, the detector measures background radiation at a rate of 200 events/second.  The source produces two types of events:  
1.  single gamma rays with a rate of 100 events/second 
2. events that produce both a gamma ray and a neutron at a rate of 100 events/second.  The gamma ray is detected immediately after creation with an efficiency of 50% (for both gamma rays from the source- the singles and the gamma/neutron coincidences).  The neutron is detected after a capture time, which is quantified with a mean of 250 microseconds. In math, this says that the time distribution of neutron captures is $$P(t)=250\times 10^{-6} e^{-t/(250\times 10^{-6})},$$ where t is the time difference (in seconds) between gamma/neutron creation and neutron capture. The neutron detection efficiency is. 40%.    The time distribution of neutron captures is $$P(t)=250\times 10^{-6} e^{-t/(250\times 10^{-6})},$$ where t is the time difference (in seconds) between gamma/neutron creation and neutron capture.

Write a simple MC to model 1000 seconds of running time.   Generate the times for background events, the single gamma rays from the source and the gamma/neutron coincidences.  Then generate the neutron capture times and finally apply the detection efficiencies for each type of event.

Histogram and plot:  

1.  The times of each type of event (on an x-axis with a scale of 0-1000 seconds)
2.  The intervals between a detected event and the next detected event (from 0-10 ms).  To do this you might need to make an array with the times of all detections, sort it and then calculate the time from one event to the next.
3.  Repeat plots 1 and 2 for the case where the background rate is 1000 events/second.


```python
"""
Plan:
We need the times of each type of event. 
Let's say the source produces events perfectly. 
Sample from the uniform distribution then find the interavls
"""
```




    "\nPlan:\nWe need the times of each type of event. \nLet's say the source produces events perfectly. \nSample from the uniform distribution then find the interavls\n"




```python
N = 100000

efficiency_test = np.random.rand(N)
times = np.random.uniform(0, 1000, size=N)


plt.hist(times, bins=50, label='source/detector events')
add_labels('time of event', 'counts', 'times for event 1 (gamma)')
```


```python
intervals = []
times = np.sort(times)
for i in range(len(times) - 1):
    intervals.append(times[i+1] - times[i])
    
plt.hist(intervals, bins=50)
add_labels('interval between events', 'counts', 'time interval for event 1 (gamma)')
```


```python
N = 100000

efficiency_test = np.random.rand(N)
efficiency_test_2 = np.random.rand(N)
times = np.random.uniform(0, 1000, size=N) #times for event 2.
mask = efficiency_test > 0.5
mask2 = efficiency_test_2 > 0.4

plt.hist(times, bins=50, label='source events',alpha=0.3)
plt.hist(times[mask], bins=50, label='gamma detector events', alpha=0.7)
plt.hist(times[mask2], bins=50, label='neutron detector events', alpha=0.7)
add_labels('time of event', 'counts', 'times for event 2 (gamma+neutron)')
```


```python
meanTime = 250*10**-6

times = np.random.uniform(0, 1000, size=N) #times for event 2.
mask = efficiency_test > 0.5
mask2 = efficiency_test_2 > 0.4

gamma_events = np.sort(times[mask])
neutron_events = np.sort(times[mask2])


detector_events = neutron_events + np.random.exponential(meanTime, size = len(neutron_events)) 

intervals_gamma = []
intervals_neutron = []
detector_events = np.sort(detector_events)
for i in range(len(detector_events) - 1):
    intervals_neutron.append(detector_events[i+1] - detector_events[i])
    
for i in range(len(gamma_events) - 1):
    intervals_gamma.append(gamma_events[i+1] - gamma_events[i])
    
plt.hist(intervals_neutron, bins=50, label='neutron detector events', alpha = 0.7)
plt.hist(intervals_gamma, bins=50, label='gamma detector events', alpha=0.7)
add_labels('interval between events', 'counts', 'intervals between times for event 2 (gamma+neutron)')
```


```python
"""Will repeat for the background rate later. Gotta do the last two questions"""
```




    'Will repeat for the background rate later. Gotta do the last two questions'




```python

```

# Solving and Fitting

An object is viewed by 6 theodolites (telescopes that allow you to measure angles with respect to the vertical($\theta$) and with respect to a horizontal($\phi$).  The locations of the telescopes and the measured angles are given in the code cell below (all with respect to a single coordinate system,) which you can run in order to avoid cutting and pasting. 

Since each theodolite gives you two quantities we have 12 measurements.  There are three unknowns- the x, y, and z coordinates of a point on an object.  Consequently the system of equations is overdetermined, and we would typically be unable to solve.  

In a case like this, we "solve" the system of equations by fitting.  

Write a function that calculates the theta and phi angles (spherical coordinates) given an arbitrary point position and an theodolite position.

Using this function, write a second function that calculates the residuals of the calculation with respect to the data. Then find the position of the target by fitting the calculations to the data.  

Extract the uncertainties on the position.  In order to do this, you will need to estimate the size of the uncertainty in the angles to incorporate into the residual calculation.  One way to do this is by histogramming the residuals ((data-calculation)/uncertainty) and getting a sense of how big they are- if the uncertainties are correct the RMS/width of the residuals distribution should be about 1.    



```python
import numpy as np
positions =  np.array([[-0.0114519 ,  9.99999344,  6.69130799],
       [ 6.33116335, -7.74056655,  3.59009715],
       [-9.85887257, -1.67410623,  4.56138703],
       [ 4.09466118,  9.12325325,  7.32437587],
       [-3.37815507, -9.41212348,  3.67927878],
       [-5.4274145 , -8.39899827,  6.8865273 ]])
thetaMeasures = np.array( [2.3164460126412245, 1.6313464103926718, 1.6970719353033004, 2.533624050255211, 1.6289243873579955, 1.8159707640506253] )
phiMeasures = np.array( [-0.9844385514247201, 1.8432675248986576, 0.5831709514637226, -2.132723211863188, 1.2215622389804748, 1.0853416750730511] )

print('Theodolite x  \t Theodolite y \t Theodolite z \t Meas.theta \t Meas. phi') 
for i in range(0,6):
    print('%f  \t %f \t %f \t %f \t %f '%(positions[i,0],positions[i,1],positions[i,2],thetaMeasures[i],phiMeasures[i]))

```

    Theodolite x  	 Theodolite y 	 Theodolite z 	 Meas.theta 	 Meas. phi
    -0.011452  	 9.999993 	 6.691308 	 2.316446 	 -0.984439 
    6.331163  	 -7.740567 	 3.590097 	 1.631346 	 1.843268 
    -9.858873  	 -1.674106 	 4.561387 	 1.697072 	 0.583171 
    4.094661  	 9.123253 	 7.324376 	 2.533624 	 -2.132723 
    -3.378155  	 -9.412123 	 3.679279 	 1.628924 	 1.221562 
    -5.427415  	 -8.398998 	 6.886527 	 1.815971 	 1.085342 



```python
def xyz_to_spherical(vector):
    r = np.linalg.norm(vector)
    theta = np.arccos(vector[2]/r)
    phi = np.arctan2(vector[1], vector[0])
    return np.array([r, theta , phi])
```


```python
def find_theta_phi(point, telescope):
    difference = point - telescope
    return xyz_to_spherical(difference)


def find_residuals(point):
    theta_res = np.zeros((6))
    phi_res = np.zeros((6))
    for i in range(6):
        angles = find_theta_phi(point, positions[i, :3])
        theta_res[i] = thetaMeasures[i] - angles[1]
        phi_res[i] = phiMeasures[i] - angles[2]
    
    
    return np.concatenate([theta_res, phi_res], axis=0)
```


```python
from scipy.optimize import least_squares

params0 = np.array([1, 1, 1])
popt_minimization = least_squares(find_residuals, params0)
popt_minimization

```




     active_mask: array([0., 0., 0.])
            cost: 3.653545714200531e-19
             fun: array([-3.29253513e-10, -9.01634323e-12,  5.41442446e-11,  5.54386315e-10,
           -2.39801290e-10,  1.09330989e-10, -1.75378601e-10,  9.24476051e-11,
            7.91083865e-11,  3.49196672e-10,  6.83249013e-11,  2.65174105e-10])
            grad: array([-9.43717308e-11,  3.27579658e-11,  1.00605068e-11])
             jac: array([[ 0.06381229, -0.09606042,  0.12487658],
           [-0.00110712,  0.0039622 ,  0.06786023],
           [ 0.00711269,  0.00469231,  0.06712066],
           [-0.07762667, -0.12328795,  0.10138615],
           [ 0.00117931,  0.00323844,  0.05922442],
           [ 0.00656647,  0.01244677,  0.05624395],
           [-0.1927295 , -0.12802894,  0.        ],
           [ 0.06559699,  0.01832914,  0.        ],
           [ 0.03755727, -0.05692997,  0.        ],
           [-0.26295946,  0.16556903,  0.        ],
           [ 0.05583779, -0.02033393,  0.        ],
           [ 0.05285997, -0.027887  ,  0.        ]])
         message: '`gtol` termination condition is satisfied.'
            nfev: 6
            njev: 6
      optimality: 9.437173084897083e-11
          status: 1
         success: True
               x: array([2.38, 6.4 , 2.7 ])




```python
uncertainty = 5e-10
plt.hist(find_residuals(popt_minimization.x)/uncertainty)
print(f'the uncertainty is about {uncertainty}')
```

# Ordinary Differential Equations

Consider the mechanics problem shown in the diagram below.  We have a rigid bar $l=4$ meters long, with mass $m=20$ kg, at an angle $\theta$ to the floor and leaning against a wall.  The center of mass of the bar is at coordinates (x,y).  The normal force against the wall we designate as $N_1,$ while the normal force to the floor is $N_2.$  The coefficient of (sliding) friction with the floor is $\mu_2=0.2$, while the coefficient of friction with the wall is $\mu_1=0.4.$ 




Calculate how long it takes for the bar to fall to the floor.  The acceleration due to gravity is $g.$  Start with an initial condition of $\theta_0= 1.15$ radians, and $\omega_0=\dot{\theta}(0)=0.$

From the mechanics perspective, this involves setting up Newton's second law for the acceleration of the center of mass, and using the torque equation for the rotation of the bar around the center of mass:
\begin{align}\\
-mg+N_2+\mu_1 N_1&=m\ddot{y}\\
N_1-\mu_2 N_2=m\ddot{x}\\
(N_1+\mu_2 N_2)\frac{l}{2} \sin\theta - (N_2-\mu_1 N_1) \frac{l}{2}\cos\theta &=I\ddot{\theta}\\
\end{align}

For a bar, the moment of inertia $I=\frac{1}{12}ml^2.$

The geometry of the situation means that we can't simply plug the equation above into an ODE integrator.  Instead, we need to relate $x$ and $y$ to the angle $\theta:$
\begin{align}\\
x&=\frac{l}{2}\cos\theta\\
\dot{x}&=-\frac{l}{2}\sin\theta \dot{\theta}\\
\ddot{x}&=-\frac{l}{2}\sin\theta \ddot{\theta}-\frac{l}{2}\cos\theta \dot{\theta}^2\\
y&=\frac{l}{2}\sin\theta\\
\dot{y}&=\frac{l}{2}\cos\theta \dot{\theta}\\
\ddot{y}&=\frac{l}{2}\cos\theta \ddot{\theta}-\frac{l}{2}\sin\theta \dot{\theta}^2\\
\end{align}

Now substitute these equations for $\ddot{x}$ and $\ddot{y}$ into the original equation for Newton's laws.  
Solve the ODE numerically.  For each step, you will need to solve the set of simultaneous equations for $\ddot{\theta}$, $N_1,$ and $N_2$.  

Plot $\theta$ as a function of time $t$.  Graphically find the time the bar hits the floor ($\theta=0.$)
Plot $N_1$ and $N_2$ as a function of time $t$.  Discuss any mathematical assumptions in our model that might not reflect reality and how the model would need to change to accommodate them.  




```python
m = 20
l = 4
mu2 = 0.2
mu1 = 0.4
I = (1/12)*m*l**2
g=9.81
```


```python
global N1
global N2
global ddtheta
```


```python
#x = l/2*np.cos(theta)
#y = l/2*np.sin(theta)


```


```python
"""This is a coupled system. 
Hence we need to update both the linear newtons law
and the angular newtons law. we know the linear newtons law is 0, -mg.

I can't seem to get the height of the ground

"""
```


```python
N1 = []
N2 = []

def main(y, t):
    "Main thing to integrate"
    
    """
    y[0] - theta
    y[1] - omega
    y[2:4] - x,y 
    y[4:6] - dx, dy
    """
    
    theta = y[0:2]

    positions = y[2:4]
    velocities = y[4:6]
    
    ddx = 0
    ddy = -m*g
    
    #Solving for N. Need ddx, ddy
    A = np.array([[mu1, -m*g],[1, -mu2]])
    b = np.array([m*ddy, m*ddx])
    N = np.linalg.solve(A, b)
    
    N1.append(N[0])
    N2.append(N[1])
    
    dydt = np.empty(6)
    
    dx = -l/2*np.sin(y[0])*y[1]
    dy = l/2*np.sin(y[0])*y[1]
    
    dydt[0] = y[1] #omega
    dydt[1] = (1/I)*((N[0] + mu2*N[1])*(l/2)*np.sin(y[0]) - (N[1] - mu1*N[0])*(l/2)*np.cos(y[0]))
    dydt[2:4] = np.array([dx, dy])
    dydt[4:6] = np.array([ddx, ddy])
 
    return dydt 
```


```python
t = np.linspace(0, 10, 10000)
y0 = np.array([1.15, 0, l/2*np.cos(1.15), l/2*np.sin(1.15), 0, 0], dtype='float64')

y = odeint(main, y0, t)

"""Grab all physical numbers: the bar cannot fall into the floor"""
mask = (y[:, 3] > 0) & (y[:, 0] > 0)
t = t[mask]
y = y[mask]
```


```python

```

When you are finished, upload your completed exam onto eClass.  Double check that the uploaded exam contains your latest work!


```python

```


```python
plt.plot(t, y[:, 0]/np.pi*180)
add_labels('t','$\\theta$ (degrees)','angle of bar')
```


```python
t = np.linspace(0, 10, len(N1))
plt.plot(t, N1, label="N1")
plt.plot(t, N2, label="N2")
add_labels('time', 'N1 N2', 'N1, N2 as a function of time')
```

## Problems with the Model

N1, N2 are roughly constant which I think is not really correct, to get N1, N2 I had to let ddy = -mg which is gravity on the CoM. Otherwise there is no way to solve the problem due to recursion. We don't have enough variables.

I think that an good simulation should have friction be proportional to the velocity

If we plot the height of the bar, this is also not really correct as well. I think what's causing the problem is probably this assumption that friction is constant, when in reality, it's far from constant. So we know this model is wrong but I think I did the integration correctly.



```python
t = np.linspace(0, 10, 10000)
y0 = np.array([1.15, 0, l/2*np.cos(1.15), l/2*np.sin(1.15), 0, 0], dtype='float64')

y = odeint(main, y0, t)

"""Grab all physical numbers: the bar cannot fall into the floor"""
mask = (y[:, 3] > 0) & (y[:, 0] > 0)
t = t[mask]
y = y[mask]

plt.plot(t, y[:, 3]/np.pi*180)
add_labels('t','y','height of bar')
```


```python

```

5ed77d8fe639247917e65e8bee223452add87238
ipynb
Jupyter Notebook
Final2019.ipynb
hanzhihua72/phys-420
748d29b55d57680212b15bb70879a24b79cb16a9
MIT
Final2019.ipynb
hanzhihua72/phys-420
748d29b55d57680212b15bb70879a24b79cb16a9
MIT
Final2019.ipynb
hanzhihua72/phys-420
748d29b55d57680212b15bb70879a24b79cb16a9
MIT
Qwen/Qwen-72B
1. YES
2. YES

