Solutions to Exercises

# Solution 6.1
# Data
V = np.array([25.8, 26.2, 26.0, 26.5, 25.8, 26.1, 25.8, 26.3])

# Edit
Vmean = np.mean(V)
Vstd = np.std(V,ddof=1)
Verr = Vstd/np.sqrt(len(V))

#Print
print('(a) The average is: {} ml' .format(round(Vmean,2)))
print('(b) The standard deviation is: {} ml'.format(round(Vstd,2)))
print('(c) The uncertainty in the mean of the volume: {} ml'.format(round(Verr,2)))
print('(d) The best estimate of the true volume is: {:.2f} \u00B1 {:.2f} ml'.format(Vmean,Verr))

# Solution 6.1*
# add your calculations here

#calculate mean and std 'manually'
Tree_av = np.sum(L)/len(L)
Tree_std = np.sqrt(np.sum((L-Tree_av)**2)/(len(L)-1))

print('The best estimate for the height %.0f +/- %.0f m' %(Tree_av, Tree_std/np.sqrt(len(L))))

#Compare values with numpy functions, if true, nothing happens
assert Tree_av == np.mean(L), 'The results for the mean are different'
assert Tree_std == np.std(L,ddof=1), 'The results for the standard deviation are different'

# Solution 6.2
opg2_1 = 'iii'
opg2_2 = 'vii'
print('The measurements adequately displayed are: {} en {}'.format(opg2_1,opg2_2))

#Solution 6.3
#Data
numbers = np.array([602.20,0.00135,0.0225,1.60219,91.095])

#Print
print('The correct notation of (a) is: {:.1e}'.format(numbers[0]))
print('The correct notation of (b) is: {:.1e}'.format(numbers[1]))
print('The correct notation of (c) is: {:.1e}'.format(numbers[2]))
print('The correct notation of (d) is: {:.1e}'.format(numbers[3]))
print('The correct notation of (e) is: {:.1e}'.format(numbers[4]))

#Solution 6.4
data = np.loadtxt('experimental_results.csv')

N = len(data)
mean_val = np.mean(data)

plt.hist(data,bins=np.linspace(0,2*mean_val,50))
plt.xlabel('Value (-)')
plt.ylabel('Counts (-)')
plt.show()

# Calculation of mean and standard deviation
mean_array = np.zeros(N)
sigma_array = np.zeros(N)

for i in range(N):
    mean_array[i] = np.mean(data[:i+1])
    sigma_array[i] = np.std(data[:i+1], ddof=1)/(i+1)
    
# Plotting
plt.figure()
plt.plot(mean_array)
plt.axhline(mean_val,linestyle='dashed',color='grey')
plt.grid()
plt.xlabel('Data length (-)')
plt.ylabel('Average (-)')
plt.show()

plt.figure()
plt.loglog(sigma_array)
plt.grid()
plt.xlabel('Data length (-)')
plt.ylabel('Uncertainty (-)')
plt.show()

#Solution 6.5

# Input
mu = 502                   # mean
sigma = 14                 # standard deviation

# Values to be measured
x1 = 500
x2 = 530

# Calculation
Pa = norm.cdf(x1,mu,sigma)      # calc: Pa = 0.44320150
Pb = 1 - norm.cdf(x2,mu,sigma)  # calc: Pb = 1 - 0.97724987
aantal = Pb * 1000

#Print
print('(a) P = {:.4f}'.format(Pa))
print('(b) Expected amount = {:.2f}'.format(aantal))

#Solution 6.6

# a. Each series of measurements consists of five repeated measurements.
Measurements = np.array([[4.560,4.470,4.620,4.570,4.570],\
                     [4.570,4.540,4.700,4.510,4.150],\
                     [4.540,4.570,4.550,4.430,4.460],\
                     [4.470,5.130,4.490,5.150,4.530],\
                     [4.400,5.050,4.490,5.080,4.530],\
                     [4.530,5.130,4.560,4.910,4.560]]) #s

# b
# Random error: reaction time, teacher walks in front of the setup
# Systematic error: inertia of lamp in turning off and on (the lamps keeps glowing for a certain amount of time)

# c. The average of the measurement series
Measurements_av = np.array([np.mean(Measurements[0,:]),np.mean(Measurements[1,:]),\
                       np.mean(Measurements[2,:]),np.mean(Measurements[3,:]),\
                       np.mean(Measurements[4,:]),np.mean(Measurements[5,:])])


print('The average times, in s, are: ', np.round(Measurements_av,1))


# Standerd deviation of the measurements series
Measurements_std = np.array([np.std(Measurements[0,:],ddof=1),np.std(Measurements[1,:],ddof=1),\
                         np.std(Measurements[2,:],ddof=1),np.std(Measurements[3,:],ddof=1),\
                         np.std(Measurements[4,:],ddof=1),np.std(Measurements[0,:],ddof=1)])

print('The standard deviation in the times, in s, are: ', np.round(Measurements_std,1))
print('')


# d. Standard deviation in the mean
std_van_gem = np.std(Measurements_av,ddof=1)
print('The standard deviation in the mean is: ', np.round(std_van_gem,1), 's')


# Uncertainty in the measurement series
Measurements_unc = Measurements_std/np.sqrt(5)
print('The uncertainty in each measurement series: ', np.round(Measurements_unc,2), 's')
print('The mean uncertainty: ', np.round(np.mean(Measurements_unc),2), 's')
print()


# e. Final answer
print('The time during which the lamp is on: ', round(np.mean(Measurements),2), '+/-', round(np.std(Measurements,ddof=1)/np.sqrt(30),2),'s')
#Note: We act like the the 6 measurement series are just one big data set

# Solution 6.7

E = 14.5
u_E = 0.2
L = 14.9
u_L = 0.1

if (E-L) > 2*np.sqrt(u_E**2+u_L**2):
    print('No agreement')
else:
    print('The values are in agreement')

#Solution 6.8

#Data
Q = np.array([45.7, 53.2, 48.4, 45.1, 51.4, 62.1, 49.3])
N = len(Q)        # Total number of measurements

#Number to be determined
n = 6             # Number of measurement

#Editing
Qmean = np.mean(Q)
Qstd = np.std(Q,ddof=1)
Qerr = Qstd/np.sqrt(N)
Pout = 1-norm.cdf(Q[n-1],Qmean,Qstd)

#Print
print('The average is: {:.4f} uC'.format(Qmean))
print('The standard deviation is: {:.4f} uC'.format(Qstd))
print('The standard error is: {:.4f} uC'.format(Qerr))
print('Pout = {:.4f}'.format(Pout))
if 2*Pout*N < 0.5:
    print('Measurement {} has to be excluded'.format(n))
else:
    print('Measurement {} cannot be excluded'.format(n))

#Solution 6.8*
print('The 4.5 m may not be a faulty measurement, and can therefore not be discarded.')
print('The 4.5 m deviates from the other tree heights and may be another type or may be planted later.')

#Solution 6.9

#Data
N = np.array([0,1,2,3,4,5,6,7,8,9,10,11,12,13])
f = np.array([0,1,0,2,3,6,9,11,8,8,6,2,1,1])

#Calculations (a)
N_tot = sum(N*f)              # (a) Total number of counts
print('(a) The total number of counts is {:.0f}'.format(Ntot))
print()

#Calculations (b)
N_mean = sum(N*f)/sum(f)        # (b) Average number of counts
print('(b) The average number of counts is {:.2f}'.format(N_mean))
print()

#Calculations (c)
Nstd = np.sqrt(N_mean)              # (c) Standard deviation
Nerr = Nstd / np.sqrt(N_tot)

print('(c) The standard deiation is {:.2f}'.format(Nstd))
print()

#Calculations (d)
P=0
for i in range(6):
    P += N_mean**i*np.exp(-N_mean)/math.factorial(i)
print('(d) The expected number is {:.0f}'.format(P*58))

#Solution 6.9*

persons = np.genfromtxt("person_count.dat") #! aantal_peronen = persons
persons_mean = np.mean(persons)
persons_std = np.std(persons, ddof=1)

print('Average: %.3f \n Standard deviation: %.3f' %(persons_mean,persons_std))

num_bins = int(max(persons)-min(persons))

x = np.linspace(min(persons),max(persons),num_bins+1)
prob = poisson(persons_mean)

f,g,_ = plt.hist(persons,bins=num_bins,density=1)
plt.plot(x+.5,prob.pmf(x),'k.')
plt.xlabel('Number of visitors per minute')
plt.ylabel('Probability')
plt.show()

print('It looks like a Poisson distribution. It concerns a counting experiment, so that makes sense.')

# Solution 6.10

def func2a(a):
    z = np.sqrt(a)
    return z

def calc2a(a,ua):
    uz = np.abs( (1/2) * a**(-1/2) ) * ua
    return uz

def func2b(a):
    z = np.exp(a**2)
    return z

def calc2b(a,ua):
    uz = np.abs( 2 * a * np.exp(a**2) ) * ua
    return uz

a = 1.274
ua = 0.005

za = func2a(a)
uza_func = ( func2a(a+ua)-func2a(a-ua) ) / 2
uza_calc = calc2a(a,ua)
zb = func2b(a)
uzb_func = ( func2b(a+ua)-func2b(a-ua) ) / 2
uzb_calc = calc2b(a,ua)

# The number of significant figures is only given for demonstrative purposes
print('(a) functional: Z = {:.4f} \u00B1 {:.6f}'.format(za,uza_func))
print('      calculus: Z = {:.4f} \u00B1 {:.6f}'.format(za,uza_calc))
print('(b) functional: Z = {:.4f} \u00B1 {:.5f}'.format(zb,uzb_func))
print('      calculus: Z = {:.4f} \u00B1 {:.5f}'.format(zb,uzb_calc))

# Solution 6.11
# add your calculations here

# data
a = 12.3
ua = 0.4

b = 5.6
ub = 0.8

c = 89.0
uc = 0.2

# (a)
za = a - b
uza = np.sqrt( ua**2 + ub**2 )

# (b)
zb = a*b/c
uzb = zb * np.sqrt( (ua/a)**2 + (ub/b)**2 + (uc/c)**2 )

# (c)
zc = np.exp(a*b/c)
dzc_da = (b/c) * np.exp(a*b/c)
dzc_db = (a/c) * np.exp(a*b/c)
dzc_dc = (-a*b/c**2) * np.exp(a*b/c)
uzc = np.sqrt( dzc_da**2 * ua**2 + dzc_db**2 * ub**2 + dzc_dc**2 * uc**2 )

print('(a) Z = {:.1f} \u00B1 {:.1f}  (calculus)'.format(za,uza))
print('(b) Z = {:.1f} \u00B1 {:.1f}  (simplified calculus)'.format(zb,uzb))
print('(c) Z = {:.1f} \u00B1 {:.1f}  (calculus)'.format(zc,uzc))

# Solution 6.12
# add your calculations here
def func1(x):                      #first function
    return (x-1)/(x+1)

def func2(r):                      #second function
    return np.exp(r**2)

def functional(x,sig,f):           #function for the functional method
    return (f(x+sig) - f(x-sig))/2

#first function 
x = 3.2
sig = 0.2
error_functional = functional(x,sig,func1)
error_calculus = sig*((x+1)-(x-1))/(x+1)**2 

print('First function \n\
Value of Z: %.2f \n\
Error with functional approach: %.2f \n\
Error with calculus approach: %.2f \n\
' %(func1(x),error_functional,error_calculus))

#second function
r = 8.745
sig = 0.005
error_functional = functional(r,sig,func2)
error_calculus = sig*2*r*func2(r)  


print('Second function \n\
Value of Z: %.1e \n\
Error with functional approach: %.0e \n\
Error with calculus approach: %.0e \
' %(func2(x),error_functional,error_calculus))

#Solution 6.13
#import data
H_W_data = np.genfromtxt('weight-height.csv',delimiter=',',dtype=float,skip_header=1)
height = []
weight = []

for i in H_W_data:
    height.append(i[0]*2.54)   #convert to cm
    weight.append(i[1]*0.4536) #convert to kg

#Calculate mean and std
height_mean = np.mean(height)
height_std = np.std(height, ddof=1)
weight_mean = np.mean(weight)
weight_std = np.std(weight)

#Linear function to fit
def lin(x,a,b):
    return a*x + b

#Curvefit
val, cov = curve_fit(lin,height,weight)

x = np.linspace(145,205)
y = lin(x,val[0],val[1])

#Scatterplot
plt.figure(figsize=(12,4))
plt.plot(height,weight,'.',label='data')
plt.plot(x,y,label='fit')
plt.xlabel('Height (cm)')
plt.ylabel('Mass (kg)')
plt.xlim(min(x),max(x))
plt.show()

#Histogram plot
N=200
x2 = np.linspace(np.min(height),np.max(height),N)

plt.figure()
x1, y1, _ = plt.hist(height,bins=100)
dist = norm.pdf(x2,height_mean,height_std)
plt.plot(x2,dist*x1.max()/np.max(dist),'r')
plt.xlabel('Height (cm)')
plt.ylabel('Occurence')
plt.show()

#Histogram plot
N=50
x2 = np.linspace(np.min(weight),np.max(weight),N)

plt.figure()
x1, y1, _ = plt.hist(weight,bins=100)
dist = norm.pdf(x2,weight_mean,weight_std)
plt.plot(x2,dist*x1.max()/np.max(dist),'r')
plt.xlabel('Mass (kg)')
plt.ylabel('Occurence')
plt.show()

#Solution 6.14

m = np.array([50,100,150,200,250,300])
T = np.array([2.47,3.43,4.17,4.80,5.35,5.86])


#Define Function to fit
def Per(m,dm,C):
    return 2*np.pi*np.sqrt((m+dm)/C)

#Make figure
fig = plt.figure(figsize=(10,5))
ax = fig.add_subplot(121)
ax.errorbar(m,T,linestyle='none',marker='o')
ax.set_xlabel('m [kg]')
ax.set_ylabel('T [s]')

#Make Fit
vals, uncmat = curve_fit(Per,m,T)

#Make linspace
linm = np.linspace(0.7*np.min(m),1.1*np.max(m),200)

#Plot the fitted line
ax.plot(linm,Per(linm,*vals))

#Set the x-limit
ax.set_xlim(np.min(linm),np.max(linm))

#Find the uncertainties
u_dm, u_C = np.sqrt(np.diag(uncmat))


#Print the values gotten and its uncertainties
print('dm = {:.1f} +/- {:.1f} m'.format(vals[0],u_dm))
print('C = {:.1f} +/- {:.1f} kg m/s^2'.format(vals[1],u_C))

#Make a residual analysis
r = T - Per(m,*popt)

#Plot the risiduals
ax2 = fig.add_subplot(122)
ax2.plot(m,r,linestyle='none',marker='o')
ax2.set_xlabel('m [kg]')
ax2.set_ylabel('R [s]')

#Make it a bit pretier
plt.tight_layout()

#Solution 6.15
#function for gravitational force
def FG(G,m1,m2,r):
    return G * m1 * m2 /(r*r)

#values
G = 6.6759e-11
u_G = 0.00030e-11

m1 = 4.739e8
u_m1 = 0.154e8

m2 = 5.9722e24
u_m2 = 0.0006e24

r = 2.983e6
u_r = 0.037e6

#value of gravitatonal force
F_m = FG(G,m1,m2,r)

#calculus appraoch
r2 = r**2
u_F2 = (m1*m2/r2*u_G)**2 + (G*m2/r2*u_m1)**2 + (G*m1/r2*u_m2)**2 + (-2*G*m1*m2/(r2*r)*u_r)**2
u_F_calc = np.sqrt(u_F2)

#Func
U_F2 = ((FG(G+u_G,m1,m2,r)-FG(G-u_G,m1,m2,r))/2)**2 + ((FG(G,m1+u_m1,m2,r)-FG(G,m1-u_m1,m2,r))/2)**2 + ((FG(G,m1,m2+u_m2,r)-FG(G,m1,m2-u_m2,r))/2)**2 + ((FG(G,m1,m2,r+u_r)-FG(G,m1,m2,r-u_r))/2)**2
u_F_func = np.sqrt(U_F2)

print("F = %.2f +/- %.2f 10^10 N \n" %(F_m/1e10,u_F_func/1e10))

diff = np.abs(u_F_calc-u_F_func)

#onzekerheid in veel decimale cijfers om verschil te tonen
print("Uncertainties \n\
Caluculus approach: %f 10^8 N \n\
Functional approach: %f 10^8 N \n\
Difference: %f 10^8 N" %(u_F_calc/1e8,u_F_func/1e8,diff/1e8))

#Solution 6.16
data = np.genfromtxt('meting_massaveer.dat')
t = np.linspace(0,5,400)

def sine(x):                        #sine function to fit the data
      return 4.5*np.sin(x*5*np.pi)

#plot of data and sine function
plt.figure(figsize=(12,4))
plt.title('Data')
plt.plot(t,data, '.', label='data')
plt.plot(t,sine(t), label='$4.5\sin(5\pi x$)')
plt.xlim(min(t),max(t))
plt.ylabel('Position (cm)')
plt.xlabel('Time (s)')
plt.legend()
plt.show()

R = data - sine(t)  #Residuals

#plot of Residuals
plt.figure(figsize=(12,4))
plt.title('Residuals')
plt.plot(t,R)
plt.ylabel('Position (cm)')
plt.xlabel('Time (s)')
plt.show()

plt.hist(R[:100])
plt.show()

# Solution (PRE/POST EXERCISE)
I = np.array([3034,1494,756.1,384.9,199.5,100.6,39.93,20.11,10.23,5.00,2.556,1.269,0.601,0.295,0.137,0.067])*1e-6
a_I = np.array([4,2,0.8,0.4,0.3,0.2,0.05,0.03,0.02,0.01,0.008,0.007,0.007,0.006,0.006,0.006])*1e-6
V = np.array([670,636,604,572,541,505,465,432,400,365,331,296,257,221,181,145])*1e-3
a_V = np.array([4,4,4,4,4,4,3,3,3,3,3,2,2,2,2,2])*1e-3

def current(V,a,b):
    return a*(np.exp(b*V)-1)

popt, pcov = curve_fit(current,V,I)
x = np.linspace(0.1,0.7)
y = current(x,popt[0],popt[1])

plt.figure(figsize=(8,4))
plt.errorbar(V,I,xerr=a_V,yerr=a_I,linestyle='none',marker= '.',label='data')
plt.plot(x,y,label='fit')
plt.yscale('log')
plt.xlabel('Voltage (V)')
plt.ylabel('Current (A)')
#plt.legend()
plt.xlim(0.1,0.7)
plt.show()
print('The model is in good agreement with the measurements, it is therefore a valid model.')


# agreement analysis
print('a = %.3g +/- %.1g' %(popt[0], np.sqrt(pcov[0,0])))

a_2 = 2.0e-9
u_a_2 = 5e-10

if 2*np.sqrt(u_a_2**2 + pcov[0,0]) > np.abs(a_2 - popt[0]):
    print('The value of $a$ found in the experiment is in agreement with the value from theory,')
else:
    print('The value of $a$ found in the experiment does not agree with the value from theory,')