# Physics 420/580 Final Exam 
## December 19, 2018 9am-12am

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.  

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

```

# Monte Carlo

Write a Monte Carlo that calculates that calculates the area of a circle of radius 1 cm by randomly picking points from a square 2cm x 2 cm and checking to see if the points are within 1 cm of the center.




```python
import numpy as np
xy=2*np.random.rand(100000,2)-1
n=0
for aa in xy:
    if np.linalg.norm(aa)<=1:
        n=n+1
area= n/len(xy)*4
print(area)
import matplotlib.pyplot as plt
plt.figure()
plt.hist(xy[:,1])
plt.show()
```

    3.1486



    <Figure size 640x480 with 1 Axes>


Write a Monte Carlo that calculates the efficiency, $\epsilon$, for detecting cosmic ray muons with a set of three detectors as shown in the figure.  Here \begin{equation} \epsilon=\frac{N_3}{N_c},\end{equation} with $N_3$=number of muons that pass through all three detectors and $N_c$ is the number that pass through the central detector.  Assume the cosmic ray flux $I$ in units of $\frac{\rm{muons}}{\rm{m}^2\rm{s}}$ is uniform in position, and that the angular distribution is \begin{equation}I(\theta)=I_0\cos^2\theta \end{equation} where $\theta$ is the angle to the vertical ($\theta=0$)points straight down.  The actual value of $I_0$ doesn't matter here since we take a ratio- so you just need to throw a large number of muons with this angular distribution.  All three detectors are 0.25 m x 1 m.  




Generate and plot 100000 muon positions uniformly through the central counter. 


```python

```

Generate and plot thetas according to the distribution for the muons.


```python

```

Calculate the efficiency.


```python
#generate Nc muons passing through the central detector
Nc=100000
xyPhi=np.random.rand(Nc,4)*np.array([1,.25,0,2*np.pi]) #x,y,z Phi

#Now we want to generate the cos**2 (theta distribution)  Easiest just to throw and pick- so we will throw a theta, 
thetaTest=np.random.rand(4*Nc,2)*np.array([np.pi/2,1]) #Generate 2 numbers- theta, random number for picking theta,
ct2=np.cos(thetaTest[:,0])**2  

thetaArray=thetaTest[(np.where(ct2>thetaTest[:,1])),0][0]
print (thetaArray.shape, thetaArray)

x=np.linspace(0,np.pi/2,25)
plt.figure()
plt.hist(thetaArray,25)
plt.plot(x,np.cos(x)**2*1500)
plt.show()

xypt=zip(xyPhi,thetaArray[:Nc])
N3=0
for a in xypt:
    x0=np.array(a[0][0:3])
    theta=np.array(a[1])
    phi=a[0][3]
    n=np.array([np.sin(theta)*np.cos(phi),np.sin(theta)*np.sin(phi),np.cos(theta)])
    #find intersection with upper plane
    d=0.125/n[2]
    xt=x0+d*n
    if 0<xt[0]<1 and 0<xt[1]<0.125:  #passed through the upper plane
        x2=x0-d*n
        if 0<xt[0]<1 and 0<xt[1]<0.125:  #passed through the upper plane
            N3=N3+1

print('The efficiency is',N3/Nc)


```

# Fitting

A pulse of sound is emitted from a position (0,0), and travels with velocity $v=2000 \frac{\rm{m}}{\rm{s}}$.  The sound is reflected from a layer, which described by \begin{equation}(x,z)=(x,f(x))=(x,ax^3+bx^2+cx+d).\end{equation}  Here x is the horizontal position and z is the depth into the earth. 

The time it takes the sound to move from (0,0) on the surface to a point on the layer $\vec{p}$ is:\begin{equation}
t_1=\frac{|\vec{p}|}{v}\end{equation}

To do the reflection, we set up coordinates with x along the boundary at point $\vec{p}=(x,f(x))$ 

\begin{align}
\hat{i}=\frac{(1,df/dx)}{|(1,df/dx)|}\\
\hat{k}=-\frac{\vec{p}-({\vec{p}\cdot\hat{i}})\hat{i}}{|\vec{p}-({\vec{p}\cdot\hat{i}})\hat{i}|}
\end{align}

In this coordinate system the reflected wave just continues in the $\hat{i}$ direction, but reflects in the $\hat{k}$ direction.   So the final direction  \begin{equation}\hat{n}=\frac{(\vec{p}\cdot\hat{i})\hat{i}-(\vec{p}\cdot\hat{k})\hat{k}}{|(\vec{p}\cdot\hat{i})\hat{i}-(\vec{p}\cdot\hat{k})\hat{k}|}.\end{equation}  

To find the position $\vec{x_s}$ and time $t$ where the pulse returns to the surface we calculate
\begin{align} t_2=-\frac{f(x)}{v\cdot n_z}\\
\vec{x_s}=(x_s, z_s)=p+t_2v\hat{n}\\
t=t_1+t_2
\end{align}

Assume that the boundary is described by parameters $a-=0.0000002,b=0.00002, c=0.0001,d=-1000$.  

1.  Calculate the time of arrival of the pulse as a function of position x (the position on the layer that is reflecting the wave).

2. Plot t as a function of x.

3. Plot the $x_s$ (the x component of the vector $\vec{x_s}$) as a function of x. 

4.  Use the arrays you calculated in (2) and (3) to make a function that interpolates t as a function of $x_s$- and  then plot t as a function of $x_s$



```python

import numpy as np
import matplotlib.pyplot as plt
#d=np.poly1d([-0.0000002,.00002,0.0001,-1000])
d=np.poly1d([-0.00000010,-.00001,0.000,-1500])
deriv=np.polyder(d)
plt.figure()
x=np.linspace(-1000,1000,200)
plt.title('Position of boundary layer')
plt.plot(x,d(x))
v=2000 #m/s


def time(x):
#    import pdb
#    pdb.set_trace()
    p=np.array([x,d(x)])
    t1=np.linalg.norm(p)/v
    i=np.array([1,deriv(x)])
    
    i=i/np.linalg.norm(i)
    k=p-np.dot(p,i)*i
    knorm=np.linalg.norm(k)
    if knorm==0:
        k=np.array([0,1])
    else:
        k=k/knorm
    n=np.dot(p,i)*i-np.dot(p,k)*k
    n=n/np.linalg.norm(n)
    t2= -d(x)/(v*n[1])
    pf=p+v*t2*n
    return np.array((t1+t2,pf[0],pf[1]))


tvec=[]
xvec=[]
zvec=[]
for xx in x:
    t,xpos,zpos=time(xx)
    tvec.append(t)
    xvec.append(xpos)
    zvec.append(zpos)

from scipy import interpolate
ff=interpolate.interp1d(xvec,tvec)

    
plt.figure()
plt.title('Arrival time as a function of position on boundary layer')
plt.plot(x,tvec)
plt.show()

    
plt.figure()
plt.title('Time of arrival as a function of position on surface')
print(time(10))
plt.plot(xvec,tvec,'o')
plt.plot(x,ff(x))
plt.xlim(-1000,1000)
plt.ylim(1.4,1.8)

tsensor=[]
xsensor=[]
for a in range(-1000,1000,200):
    xsensor.append(a)
    tsensor.append(ff(a)+np.random.normal(0,0.0005))
    print(('|%f|%6.4f|'%(xsensor[-1],tsensor[-1])))

print(' '.join(['%6.0f,'%a for a in xsensor]))
print(' '.join(['%6.4f,'%a for a in tsensor]))
    #for aa in zip(xvec,tvec):
#    print(aa[0],aa[1])


```

Now we consider a layer with different, unknown polynomial coefficients.  
Sensors are placed every 200 meters along the surface of the ground, which measure the time of arrival of the pulse with an accuracy of 0.0005 seconds.



The data returned are:

|sensor x<br>position[m]|Time[s]|
|----|----|
|-1000.000000|1.5646|
|-800.000000|1.5456|
|-600.000000|1.5274|
|-400.000000|1.5124|
|-200.000000|1.5026|
|0.000000|1.5004|
|200.000000|1.5027|
|400.000000|1.5146|
|600.000000|1.5320|
|800.000000|1.5587|

Fit the 3rd order polynomial describing the boundary layer to these data, by defining a chisquare function and using scipy.optimize.minimize.

I have made arrays xsensor, tsensor which you may use to save cutting and pasting.


```python
xsensor=np.array([ -1000,   -800,   -600,   -400,   -200,      0,    200,    400,    600,    800])
tsensor=np.array([1.5646, 1.5456, 1.5274, 1.5124, 1.5026, 1.5004, 1.5027, 1.5146, 1.5320, 1.5587])


```


```python
v=2000 #m/s

from scipy import interpolate



def time(x):
#    import pdb
#    pdb.set_trace()
    p=np.array([x,d(x)])
    t1=np.linalg.norm(p)/v
    i=np.array([1,deriv(x)])
    i=i/np.linalg.norm(i)
    k=p-np.dot(p,i)*i
    knorm=np.linalg.norm(k)
    if knorm==0:
        k=np.array([0,1])
    else:
        k=k/knorm
    n=np.dot(p,i)*i-np.dot(p,k)*k
    n=n/np.linalg.norm(n)
    t2= -d(x)/(v*n[1])
    pf=p+v*t2*n
    return np.array((t1+t2,pf[0],pf[1]))


def chi2(pars):
    global d
    global deriv
#    import pdb
#    pdb.set_trace()
    d=np.poly1d(pars)
    deriv=np.polyder(d)

    tvec=[]
    xvec=[]
    zvec=[]
    for xx in x:
        t,xpos,zpos=time(xx)
        tvec.append(t)
        xvec.append(xpos)
        zvec.append(zpos)
    ff=interpolate.interp1d(xvec,tvec,fill_value='extrapolate')
    resid=(tsensor-ff(xsensor))**2
    return np.sum(resid)

params=np.array([-0.0000002,.00002,0.0001,-1450])
import scipy.optimize
scipy.optimize.minimize(chi2,params,method='Powell').x
#[-0.00000010,-.00001,0.000,-1500]


```




    array([-2.53300215e-07, -1.96962300e-04, -2.67924833e-02, -1.49343298e+03])




```python
from scipy import poly1d

x = np.linspace(-1000, 1000)

plt.plot(x, poly1d(np.array([-1.01825371e-07, -1.39292796e-05, -9.25237786e-04, -1.49995821e+03]))(x))
```

## Baseball Trajectory



A baseball, with diameter 7.4 cm and mass 146 g, has been hit a distance of about 150 meters.  After the ball leaves the bat, the only forces are those due to gravity and air resistance (we ignore the possibility of lift forces due to the spin of the ball).  Assuming that there is no wind, and the density of the air is $\rho=1.225 \frac{\rm{kg}}/{\rm{m}^3}$, and the magnitude of the turbulent drag, $F_t$, is
\begin{equation}
F_t=\frac{1}{2}\rho v^2AC_D
\end{equation}
where $A=\pi r^2$ is the cross sectional area of the sphere, $C_D$ is the drag coefficient.

The drag force is in the opposite direction to the velocity.


Use the following table to get the drag force as a function of Reynolds Number and interpolate.

The viscosity of air, $\eta=1.827\times 10^{-5}\rm{Pa-s}$ and Reynolds Number, $R=\frac{2\rho v r}{\eta}.$

|Reynolds Number|Drag Coefficient|
|---------------|----------------|
|1|24|
|4|8|
|10|4|
|40|1.7|
|100|1.1|
|400|0.6|
|1000|0.45|
|2000|0.4|
|4000|0.35|
|10,000|0.4|
|40,000|0.45|
|100,000|0.42|
|280,000|0.39|
|400,000|0.08|
|500,000|0.08|

In a coordinate system where x is along the direction of the ball, and z is up, we can integrate Newton's second law:
\begin{equation}
\vec{F}=mg(0,0,-1)-F_t \frac{\vec{v}}{|\vec{v}|}
\end{equation}

using odeint to determine the trajectory of the ball.  

Plot where the ball hits the ground (ie. has $z=0$) as a function of the initial speed $v$ and angle $\theta$,$\vec{v_0}=v(\cos\theta,0,\sin\theta)$; vary the range of these parameters to estimate the lowest necessary velocity for a range of 150 meters.

I have put the constants into the code block below to start you off.


```python
airDensity=1.225 # kg/m**3
airPressure=101.325 #kPa
radius=0.037
airViscosity=1.827e-5   #Pa-s  

RN=np.array([1,4,10,40,100,400,1000,2000,4000,10000,40000,100000,280000,400000,500000])
DC=np.array([24,8,4,1.7,1.1,0.6,0.45,0.4,0.35,0.4,0.45,0.42,0.39,0.08,0.08])


```


```python
import scipy.interpolate
from scipy.integrate import odeint

mass=0.146 # kg

logRN=np.log(RN)  #interpolation in log ReynoldsNumber will be more sensible
plotX=np.logspace(np.log10(RN[0]),np.log10(RN[-1]),100)
linearInterpolation=scipy.interpolate.interp1d(logRN,DC) #out x-axis is log Reynold's Number
splineInterpolation=scipy.interpolate.interp1d(logRN,DC,kind='cubic')
def DragCoefficientLinear(reynoldsNumber): 
    return linearInterpolation(np.log(reynoldsNumber))
def DragCoefficientSpline(reynoldsNumber):
    return splineInterpolation(np.log(reynoldsNumber))

def ReynoldsNumber(velocity):
    '''Calculates Reynolds number for a sphere.  Because it is used 
    in an logarithmic interpolation, we impose 1 as a lower cutoff'''
    R=airDensity*velocity*radius*2/airViscosity
    return np.maximum(R,1)


plt.figure()
plt.title("Drag Coefficient versus Reynolds Number")
plt.ylabel(r'$C_D$, drag coefficient')
plt.xlabel("Reynolds Number")
plt.semilogx(RN,DC,'o',color='blue',label='Data Points')
plt.plot(plotX,DragCoefficientLinear(plotX),color='red',label="Linear Interpolation")
plt.plot(plotX,DragCoefficientSpline(plotX),color='violet',label="Cubic Spline Interpolation")
plt.legend();

area=np.pi*radius**2
velocityAtR1=airViscosity/(2*airDensity*radius)
ceta= 6*np.pi*airViscosity*radius
ct=.5*airDensity*area
def fdrag(velocity):
    v=abs(velocity)
    if v<velocityAtR1:
        force=ceta*v
    else:
        force=ct*v**2*DragCoefficientSpline(ReynoldsNumber(v))
    return force

drag=np.vectorize(fdrag)  #so that drag can accept either a scalar or an array
print(drag([0,90,100]))

def accel(y,t):
    x,z,vx,vz=y
    pos=np.array([x,0,z])
    v=np.array([vx,0,vz])
    vmag=np.sqrt(np.dot(v,v))
    F=np.array([0,0,-mass*9.8])-fdrag(vmag)*v/vmag
    dydt=np.array([vx,vz,F[0]/mass,F[2]/mass])
#    theta, omega,energy = y
#    dydt = np.array([omega, - m*g*p.lcm/p.I*np.sin(theta)-drag(p.l*omega)/(p.m*p.l)*np.sign(omega),drag(p.l*omega)*p.l*np.abs(omega)])
    return dydt


t=np.linspace(0,15,100)  #  long enough so that it always works                     

def range(vinit,theta0):
    global ii
    y0 = np.array([0, 0.0,vinit*np.cos(theta0),vinit*np.sin(theta0)])
    sol = odeint(accel, y0, t)
    ii=np.where(sol[:,1]<0)[0][0]
#  Now we interpolate from x1,y1 to x2, y2  y=y1+(y2-y1)/(x2-x1)*(x-x1)=0
# x=-y1*(x2-x1)/(y2-y1)+x1
    ground=-sol[ii-1,1]*(sol[ii,0]-sol[ii-1,0])/(sol[ii,1]-sol[ii-1,1])+sol[ii-1,0]
    return ground

vinit=100
theta0=1/np.sqrt(2)
gr=range(vinit,theta0)
ix=ii
print(ix,sol[ix-1],sol[ix],gr)

plt.figure()
plt.plot(t,sol[:,1])
plt.figure()
plt.plot(t,sol[:,0])


v0=np.linspace(10,100,10)
r=[]
for vinit in v0:
    r.append(range(vinit,theta0))
plt.figure()
plt.plot(v0,r)

t0=np.linspace(0.55,0.7,10)
r=[]
for tinit in t0:
    r.append(range(67.7,tinit))
plt.figure()
plt.plot(t0,r)

print(67.7*3.6,ReynoldsNumber(67.7))
```

## Statistics and Monte Carlo

A beam of particles is incident on a target and produces radioactive atoms.  In typical operation it produces R atoms/second.  The way in which the atoms are produced is *independent*- that is, one atom is produced at a time and there is no connection between the rate at which atoms are produced doesn't depend on how many atoms are there already.  

This problem involves writing various Monte Carlo calculations that simulate the statistics of this situation.  

A simple way to approximate this is to break the time interval into very small time steps, so that in each step only $\lambda<<1$ atoms are produced.   We throw a random number for each time step, and generate an atom if that number is less than $\lambda$.

Consider the case when R=10 atoms/second.  Pick $\lambda=0.05,$ so the time interval will be $\delta=\lambda/R=0.005 \text{s}.$  Calculate (by MC) the number of atoms created in a second.  Run the 10000 experiments, and plot the distribution of the number of atoms created in a second interval.




```python
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import matplotlib as mpl
mpl.rc('figure',dpi=250)

R=10
lam=0.05
delta=lam/R
NperSec=int(1/delta)
Nexperiments=10000
nexp=np.random.rand(Nexperiments*NperSec) # We run the experiments together- one second is one experiment
tbirth=np.where(nexp<lam)[0]*delta #array containing times at which atoms are created
nexp=np.histogram(tbirth,Nexperiments,(0,Nexperiments))
plt.figure()
aa=plt.hist(nexp[0],25,(0,25))
plt.show()
print(aa)
```

Repeat with $\lambda=0.005$. 

For this second distribution, make a log plot which shows the number of events in a 0.005 second interval.  (The first simulation, by definition, always has 0 or 1 event in such an interval)


```python
lam=0.005
delta=lam/R
NperSec=int(1/delta)
Nexperiments=10000
nexp=np.random.rand(Nexperiments*NperSec) # We run the experiments together- one second is one experiment
tbirth=np.where(nexp<lam)[0]*delta #array containing times at which atoms are created
nexp=np.histogram(tbirth,Nexperiments,(0,Nexperiments))
plt.figure()
aa=plt.hist(nexp[0],25,(0,25))
plt.show()
print(aa)
```


```python
# to plot the number of atoms in an 0.005 second interval we histogram tbirth in 0.005 second bins
nexp2=np.histogram(tbirth,Nexperiments*200,(0,Nexperiments))
plt.figure()
plt.yscale('log')
aa=plt.hist(nexp2[0],10,(0,10))
plt.show()
print(aa)
```

To do the problem exactly, we generate an expontial time separation between adjacent events.  That is\begin{equation}P(t+dt)=P(t)(1-dt*R)\rightarrow\frac{dP}{dt}=-RP\rightarrow P=e^{-Rt}\end{equation}

Thus, the probability distribution for time intervals between events is $P=e^{-Rt}.$  Generate 10000 1 second experiments with this distribution, and plot the number of events in a second.  Also plot the number of events in a 100,000 experiments with a 0.005 second interval.




```python
#We do the same as before, and will generate intervals, then sum them to get the time, #
# and then histogram them in one second bins.
scale=1.0/R
t0=np.random.exponential(scale, int(1.1*Nexperiments*NperSec)) # We run the experiments together- one second is one experiment
tbirth=np.cumsum(t0) #array containing times at which atoms are created
nexp=np.histogram(tbirth,Nexperiments,(0,Nexperiments))
plt.figure()
aa=plt.hist(nexp[0],25,(0,25))
plt.show()
print(aa)



```


```python
# to plot the number of atoms in an 0.005 second interval we histogram tbirth in 0.005 second bins
Nexperiments=100000
nexp2=np.histogram(tbirth,Nexperiments,(0,int(Nexperiments/200)))
plt.figure()
plt.yscale('log')
aa=plt.hist(nexp2[0],10,(0,10))
plt.show()
print(aa)
```


```python

```


```python

```

## Upload the exam

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

c4e69846962d5c936cb6a2d34817765e41c093ee
ipynb
Jupyter Notebook
Assignments/Final2018Master.ipynb
hanzhihua72/phys-420
748d29b55d57680212b15bb70879a24b79cb16a9
MIT
Assignments/Final2018Master.ipynb
hanzhihua72/phys-420
748d29b55d57680212b15bb70879a24b79cb16a9
MIT
Assignments/Final2018Master.ipynb
hanzhihua72/phys-420
748d29b55d57680212b15bb70879a24b79cb16a9
MIT
Qwen/Qwen-72B
1. YES
2. YES