# Physics 420/580 Midterm Exam 
## October 19, 2017 1pm-2pm

Do the following problems.  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.  

When you are finished, upload the jupyter notebook to eclass. Eclass times out after 2:05 so make sure that you upload things before then.  Also be careful to save the notebook periodically and that you upload your final exam file.



```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
from scipy.interpolate import CubicSpline

```


```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()
```

## Graphics

Plot the two curves:
$$\begin{align} y&=4x^3-3x-2\\
x&=\sin(\frac{y^4}{4}-2y^2+2)
\end{align}
$$

How many intersections are there?  Read the x,y value corresponding to the intersections from the plots.


```python
fsolve(lambda x: x**3 - 1, 9)
```




    array([1.])




```python
"""
Plan: 

-Initialize an array of xdata, ydata. 
-Create two functions and then plot the two functions.
-Find the intersections by inspection, or if time allows, using fsolve/brute force.
- Read off the x, y values.
"""
```


```python
y1 = lambda x: 4*x**3 - 3*x - 2

xdata = np.linspace(-2, 2, 1000)
plt.plot(xdata, y1(xdata))

x1 = lambda y: np.sin(y**4/4 - 2*y**2 + 2)

ydata = np.linspace(-10, 10,1000)
plt.plot(x1(ydata), ydata)

curve1_x = xdata
curve1_y = y1(xdata)

curve2_x = x1(ydata)
curve2_y = ydata

tol = 1e-4
plt.ylim(-5, 2.5)
plt.xlim(-1.5, 1.5)
```


```python
from ipywidgets import interact, interact_manual


def plot(a):
    x = np.linspace(0, 5, 50)
    plt.plot(x, np.exp(a*x))
    plt.ylim(0, 5)
    
interact_manual(plot, a=5)
```

# Reconstruction

A instantaneous flash of light occurs in a large tank of water at time $t_0$ and at position $\vec{x_0}$  

The group velocity of light in water is about $2.2 \times 10^8$ m/s.  Four sensors detect the light flash, and report a measurement of the time at which light strikes the sensor.  The locations of the sensors and the time which each sensor was hit is recorded in the table below:

|Sensor #|	x|	y|	z|		Time|
|--------|------------|------------|------------|------------|
|	|[m]|	[m]|	[m]|		[s]|
|1	|0	|0|	10	|	6.3859E-08|
|2	|8.66025404|	0|	-5|		1.1032E-07|
|3	|-4.33012702|	7.5|	-5|		7.9394E-08|
|4	|-4.33012702|	-7.5|	-5|		1.0759E-07|

Calculate the initial time $t_0$ and the location of the flash.  

Here is a sketch of the general geometry.  Note that the gray region is not different from the rest of the water- it is just to show the tetrahedral arrangement of the light sensors.  The flash, shown as the yellow star is in an arbitrary location which you want to find.  




```python
"""
Plan: 

- Pick a random point called position. Calculate the distance to the each sensors, call it d
- The optimal point minimizes d - v*t for each sensor.
- Calculate the residuals 
- Run least_squares on the residuals
- Find the optimal point and print it.
"""
```


```python
"cost function solution"

sen1 = np.array([0, 0, 10, 6.3859E-08])
sen2 = np.array([8.66025404,	0,	-5, 1.1032E-07])
sen3 = np.array([-4.33012702,	7.5,	-5,	7.9394E-08])
sen4 = np.array([-4.33012702,	-7.5,	-5,	1.0759E-07])

v = 2.2*10**8
def cost(pos):
    
    error = 0
    for sen in [sen1, sen2, sen3, sen4]:
        error += (np.linalg.norm(pos - sen[0:3]) - v*sen[3])**2
    
    return error

minimize(cost, np.random.randn(3))
```




          fun: 2.889558609330956
     hess_inv: array([[ 2.7050527 ,  1.63590179,  0.27960445],
           [ 1.63590179,  1.60674196, -0.18068716],
           [ 0.27960445, -0.18068716,  0.51680401]])
          jac: array([ 1.04308128e-06,  8.94069672e-08, -1.07288361e-06])
      message: 'Optimization terminated successfully.'
         nfev: 85
          nit: 12
         njev: 17
       status: 0
      success: True
            x: array([-7.8833527 , 10.57964442, 10.80804805])




```python
from scipy.optimize import least_squares
"residuals solution"

sen1 = np.array([0, 0, 10, 6.3859E-08])
sen2 = np.array([8.66025404,	0,	-5, 1.1032E-07])
sen3 = np.array([-4.33012702,	7.5,	-5,	7.9394E-08])
sen4 = np.array([-4.33012702,	-7.5,	-5,	1.0759E-07])

def cost(pos):
    error = np.zeros(4)
    for i, sen in enumerate([sen1, sen2, sen3, sen4]):
        error[i] = (np.linalg.norm(pos - sen[0:3]) - v*sen[3])
    
    return error

least_squares(cost, np.random.randn(3))
```




     active_mask: array([0., 0., 0.])
            cost: 1.4447793059407599
             fun: array([-0.83046438,  0.93904664, -0.97413042,  0.6075762 ])
            grad: array([-1.42800140e-05, -1.02340696e-05, -8.34587431e-07])
             jac: array([[-0.59639614,  0.80035917,  0.06112968],
           [-0.65625106,  0.41966651,  0.62706827],
           [-0.21545151,  0.1867243 ,  0.95849604],
           [-0.14636443,  0.74470816,  0.651143  ]])
         message: '`ftol` termination condition is satisfied.'
            nfev: 15
            njev: 15
      optimality: 1.4280013981682327e-05
          status: 2
         success: True
               x: array([-7.88347184, 10.57956012, 10.80804357])




```python
'''
Alternative solution

'''

class Sensor():
    def __init__(self, position, time):
        self.position = position
        self.time = time
        self.velocity = 2.2*10**8
        
    def draw_radius(self):
        r = self.velocity*self.time
        possible_vals = []
        for phi in np.linspace(0, 2*np.pi):
            for theta in np.linspace(0, np.pi):
                x = np.array([np.sin(theta)*np.cos(phi), np.sin(theta)*np.sin(phi), np.cos(theta)])
                possible_vals.append(self.position + r*x)
        return np.array(possible_vals)
        
```


```python
sen1 = Sensor(np.array([0, 0, 10]), 6.3859E-08)
sen2 = Sensor(np.array([8.66025404,	0,	-5]), 1.1032E-07)
sen3 = Sensor(np.array([4.33012702,	7.5,	-5]),	7.9394E-08)
sen4 = Sensor(np.array([-4.33012702,	-7.5,	-5]),	1.0759E-07)

circ = []
for i in [sen1, sen2, sen3, sen4]:
    circ.append(i.draw_radius())
```


```python
tol = 18
mask = np.abs(circ[0]-circ[1]-circ[2]-circ[3]) < tol

for i, j in enumerate(mask):
    if j.all() == True:
        print(i, j)
        print(circ[0][i])
```

    1010 [ True  True  True]
    [-7.04232267  4.58404406 21.25904395]
    1460 [ True  True  True]
    [-7.04232267 -4.58404406 21.25904395]


## Trajectories

An alpha-particle is a helium nucleus, with mass 6.644×10−27 kg or 4.002 u, and charge equal to twice the charge of an electron (but positive).  Calculate the trajectory of a 5MeV=(5*1.609 e-13 J) alpha particle as it moves by a gold nucleus (mass 196.966 u), with charge 79e, as a function of  the impact parameter b, below.  Assume both the alpha and the gold are point particles, and ignore special relativity.  Plot the scattering angle $\theta$ and energy loss of the alpha as a function of $b$, for values of b between 1e-16 and 1e-9m.  

The force between two charge particles is given by the Coulomb potential:
$$\vec{F}=\frac{1}{4\pi\epsilon_0}\frac{q_1q_2 (\vec{r_2}-\vec{r_1})}{|\vec{r_2}-\vec{r_1}|^3}$$ with $\epsilon_0=8.85\times10^{-12}\frac{\rm{C}^2}{\rm{N\cdot m^2}}$







```python
"""
Plan:

- Write NEwtons law as system of coupled ODE and integrate. 
- F is couloumbs law, so the acceleration is F/m
- 

"""
```


```python
eps = 8.85*10**-12
elementary_charge = 1.60*10**-19

q1 = 2*elementary_charge
q2 = 79*elementary_charge

u = 1.66*10**-27
m1 = 4*u
m2 = 196.966*u


def CouloubsLaw(r1, r2, q1, q2): 
    distance = np.linalg.norm(r2-r1)
    direction = (r2-r1)/distance    
    F = (1/(4*np.pi*eps))*(q2/distance**2)*direction
    return F
```


```python
kinetic = 5*1.609e-13
velocity = np.sqrt(2*kinetic/m1)
velocity
```




    15566607.758546297




```python
import numpy as np
from scipy.integrate import odeint

def main(y, t, m1 = m1, m2 = m2): 
    '''
    Main function to be integrated with odeint.
    Input: y (array of elements 12), t (scalar)
    Output: dydt (array of elements 12)   
    '''
    
    particle1_x = y[0:2]
    particle1_dx = y[2:4]
    
    particle2_x = np.array([0, velocity])
    
    F = CouloubsLaw(particle1_x, particle2_x, q1, q2)

    dxdt = np.zeros(4)
    
    dxdt[0:2] = particle1_dx
    dxdt[2:4] = F/m1
    
    #print(y[0:2])
    
    return dxdt

b = 1e-16
y0 = np.array([0, b, velocity, 0])
t = np.linspace(0, 3, num=500)
y = odeint(main, y0, t)

```


```python
y

plt.plot(y[:,0], y[:,1])
```


```python
energy = 1/2*m1*y[:, 2]**2 + 1/2*m1*y[:, 3]**2
plt.plot(t, energy)
```


```python
for b in np.logspace(-16, -9, num=10):
    y0 = np.array([0, b, velocity, 0])
    t = np.linspace(0, 3, num=500)
    y = odeint(main, y0, t)
    plt.plot(y[:,0], y[:,1], label=b)
    
add_labels('x', 'y', 'Scattering for various b')
```


```python
for b in np.logspace(-16, -9, num=10):
    y0 = np.array([0, b, velocity, 0])
    t = np.linspace(0, 3, num=500)
    y = odeint(main, y0, t)

    energy = 1/2*m1*y[:, 2]**2 + 1/2*m1*y[:, 3]**2
    plt.plot(t, energy, label = b)
    add_labels('x', 'y', 'Scattering for various b')
```


```python

```


```python
def f(x):
    return np.array([x**2, np.sin(x)])

xdata = np.linspace(0, 10, num=10)
ydata = f(xdata)


g = scipy.interpolate.interp1d(xdata, ydata, kind='cubic')
```


```python
x = np.linspace(0, 10)

plt.plot(x, g(x)[0])
plt.plot(x, g(x)[1])
```


```python

```

# Don't forget to upload your work to eclass!


```python

```

56722069cd1d6e945bac1f368991e53a8823f819
ipynb
Jupyter Notebook
Assignments/Finished/.ipynb_checkpoints/Midterm 2018 (1)-checkpoint.ipynb
hanzhihua72/phys-420
748d29b55d57680212b15bb70879a24b79cb16a9
MIT
Assignments/Finished/.ipynb_checkpoints/Midterm 2018 (1)-checkpoint.ipynb
hanzhihua72/phys-420
748d29b55d57680212b15bb70879a24b79cb16a9
MIT
Assignments/Finished/.ipynb_checkpoints/Midterm 2018 (1)-checkpoint.ipynb
hanzhihua72/phys-420
748d29b55d57680212b15bb70879a24b79cb16a9
MIT
Qwen/Qwen-72B
1. YES
2. YES