In [4]:
import numpy as np
import matplotlib.pyplot as plt
09.01
In [5]:
np.logspace(-20,-10,11)
Out[5]:
array([1.e-20, 1.e-19, 1.e-18, 1.e-17, 1.e-16, 1.e-15, 1.e-14, 1.e-13,
1.e-12, 1.e-11, 1.e-10])
In [6]:
np.array(10*[2.])
#
2*np.ones(10)
Out[6]:
array([2., 2., 2., 2., 2., 2., 2., 2., 2., 2.])
In [7]:
np.zeros(5, dtype=np.float32)
Out[7]:
array([0., 0., 0., 0., 0.], dtype=float32)
09.02
In [8]:
x=np.random.poisson(100,10)
print (x)
print (x[1:] - x[:-1])
[ 82 109 83 112 83 79 94 88 84 93] [ 27 -26 29 -29 -4 15 -6 -4 9]
09.03
In [9]:
year, temp = np.loadtxt('data/munich_temperatures_average_with_bad_data.txt',
unpack=True)
temp_red=temp[(temp>-90) & (temp<90)]
print (len(temp)-len(temp_red),'elements were removed')
60 elements were removed
10.01
In [10]:
year, temp = np.loadtxt('data/munich_temperatures_average_with_bad_data.txt',
unpack=True)
keep=(temp>-90) & (temp<90)
In [11]:
plt.plot(year[keep],temp[keep], marker='.',linewidth=0)
Out[11]:
[<matplotlib.lines.Line2D at 0x7d769ef5b890>]
In [12]:
plt.plot(year[keep]%1,temp[keep], marker='.',linewidth=0)
Out[12]:
[<matplotlib.lines.Line2D at 0x7d769e263b10>]
10.02
In [13]:
# define normalized Gaussian with mean 'mu' and sigma 'sig'
def gaussian(x,mu,sig):
return (1/(np.sqrt(2*np.pi)*sig))*np.exp(-(x-mu)**2/(2*sig**2))
In [14]:
# Gaussian with mean 10 and sigma sqrt(10)
x=np.linspace(-5,30,100)
d=np.random.normal(7,np.sqrt(7),10000)
_=plt.hist(d,bins=np.linspace(-5,30,36))
plt.plot(x,10000*gaussian(x,7,np.sqrt(7)),color='red',linewidth=3)
Out[14]:
[<matplotlib.lines.Line2D at 0x7d769e2e7d90>]
In [15]:
# I plot here the sample from the Poisson distribution together with the Gaussian
x=np.linspace(-5,30,100)
d=np.random.poisson(7,10000)
_=plt.hist(d,bins=np.linspace(-5,30,36))
plt.plot(x,10000*gaussian(x,7,np.sqrt(7)),color='red',linewidth=3)
Out[15]:
[<matplotlib.lines.Line2D at 0x7d769ec4e0d0>]
10.03
In [16]:
import numpy as np
import matplotlib.pyplot as plt
In [17]:
x=np.arange(-10,10,0.1)
y=np.arange(-10,10,0.1)
# Note we calculate the meshgrid first here
# this makes the computations below
# more readable
xx,yy = np.meshgrid(x,y)
r1sq=xx**2+(yy-5)**2 # work out radii
r2sq=xx**2+(yy+5)**2
b1x,b1y = - (yy-5)/r1sq, xx/r1sq # calculate magnetic field vectors
b2x,b2y = + (yy+5)/r2sq, -xx/r2sq
plt.streamplot(x,y,b1x+b2x,b1y+b2y,color='k')
Out[17]:
<matplotlib.streamplot.StreamplotSet at 0x7d769eea9940>
Solution for the coil
In [18]:
a=5 #radius of the coil
n=90 #elements on the coil
x=np.arange(-10,10,0.05)
y=np.arange(-10,10,0.05)
dphi=2*np.pi/n
xx,yy=np.meshgrid(x,y) #generate coordinate grid
bx=0*xx
by=0*yy
for phi in np.arange(0,2*np.pi,dphi):
cos=np.cos(phi)
#separation coil element-observer
rp=np.sqrt((xx-a*cos)**2+a**2*(1.0-cos**2)+yy**2)
#magnetic field components (cross product)
bx+=-a*cos*yy*dphi/rp**3
by+=-a*(a-cos*xx)*dphi/rp**3
plt.streamplot(x,y,bx,by,color='k')
Out[18]:
<matplotlib.streamplot.StreamplotSet at 0x7d769c34ce10>
In [19]:
# Example "Helmholtz coils"
a=5 #radius of the coil
n=90 #elements on the coil
x=np.arange(-10,10,0.05)
y=np.arange(-10,10,0.05)
dphi=2*np.pi/n
xx,yy=np.meshgrid(x,y) #generate coordinate grid
bx=0*xx
by=0*yy
for phi in np.arange(0,2*np.pi,dphi):
cos=np.cos(phi)
#separation coil element-observer
rp=np.sqrt((xx-a*cos)**2+a**2*(1.0-cos**2)+(yy-a/2)**2)
rp2=np.sqrt((xx-a*cos)**2+a**2*(1.0-cos**2)+(yy+a/2)**2)
#magnetic field components (cross product)
bx+=-a*cos*(yy-a/2)*dphi/rp**3
by+=-a*(a-cos*xx)*dphi/rp**3
bx+=-a*cos*(yy+a/2)*dphi/rp2**3
by+=-a*(a-cos*xx)*dphi/rp2**3
plt.streamplot(x,y,bx,by,color='k')
Out[19]:
<matplotlib.streamplot.StreamplotSet at 0x7d769c726350>
11.1
In [20]:
import glob
import os
for file in glob.glob('data/*'):
print (str(int(os.path.getsize(file)/1024+0.5))+'kB', file)
1kB data/data.txt 2149kB data/ice_data_2019_area.npy 122kB data/munich_temperatures_average.txt 42kB data/germany.txt 9kB data/autofahrt.txt 0kB data/problem_1_codons.txt 0kB data/data_new.txt 2kB data/problem_1_question_6.fasta 2kB data/problem_1_question_4_new.fasta 129kB data/autofahrt2019.txt 5kB data/new_decay_data.txt 3kB data/UID_0031246_RVC_001.tbl 1kB data/problem_1_question_5.fasta 0kB data/riddle.txt 4kB data/off 9kB data/autofahrt_nohd.txt 3kB data/decay_data.txt 2149kB data/ice_data_area.npy 0kB data/columns.txt 123kB data/munich_temperatures_average_with_bad_data.txt 9kB data/autofahrt_data.txt 73kB data/autofahrt_2018.txt
12.1, 12.2
In [21]:
def date2str(tup):
year,month,day,hour,minute,second = tup
date="{0:04d}-{1:02d}-{2:02d}".format(year,month,day)
time="{0:02d}:{1:02d}:{2:02d}".format(hour,minute,second)
s=f'{date} {time}'
return s
def str2date(s):
date, time =s.split()
year,month,day = date.split('-')
hour,minute,second = time.split(':')
return (int(year),int(month),int(day),int(hour),int(minute),int(second))
# foward and backward
print (str2date(date2str((2017,4,5,17,43,12))))
print (date2str(str2date('2017-04-05 17:43:12')))
(2017, 4, 5, 17, 43, 12) 2017-04-05 17:43:12
... and we can now also more easily format the output from getsize()
In [23]:
for file in glob.glob('data/*'):
x=os.path.getsize(file)/1024
print (f'{x:8.1f}kB {file:s}')
0.5kB data/data.txt
2148.6kB data/ice_data_2019_area.npy
121.7kB data/munich_temperatures_average.txt
42.0kB data/germany.txt
9.2kB data/autofahrt.txt
0.4kB data/problem_1_codons.txt
0.0kB data/data_new.txt
2.1kB data/problem_1_question_6.fasta
2.0kB data/problem_1_question_4_new.fasta
129.0kB data/autofahrt2019.txt
5.3kB data/new_decay_data.txt
3.4kB data/UID_0031246_RVC_001.tbl
1.0kB data/problem_1_question_5.fasta
0.3kB data/riddle.txt
4.0kB data/off
9.2kB data/autofahrt_nohd.txt
3.2kB data/decay_data.txt
2148.5kB data/ice_data_area.npy
0.2kB data/columns.txt
122.5kB data/munich_temperatures_average_with_bad_data.txt
9.4kB data/autofahrt_data.txt
72.6kB data/autofahrt_2018.txt
Monte Carlo Error propagation¶
In [24]:
# Defined constant(s)
G = 6.67384e-11 # SI units
# Define values to sample from
mean_m_1 = 40.e4
sigma_m_1 = 0.05e4
mean_m_2 = 30.e4
sigma_m_2 = 0.1e4
mean_r = 3.2
sigma_r = 0.01
Standard error propagation
In [25]:
# Compute mean and error of force
mean_f = G * mean_m_1 * mean_m_2 / mean_r ** 2
sigma_f = mean_f * np.sqrt((sigma_m_1 / mean_m_1) ** 2
+ (sigma_m_2 / mean_m_2) ** 2
+ 4. * (sigma_r / mean_r) ** 2)
print(mean_f, sigma_f)
0.7820906249999999 0.005625407272639265
Generate random samples for m_1, m_2 and r and calculate the force
In [26]:
N = 1000000
m_1 = np.random.normal(mean_m_1, sigma_m_1, N)
m_2 = np.random.normal(mean_m_2, sigma_m_2, N)
r = np.random.normal(mean_r, sigma_r, N)
# calculate force
F = G * m_1 * m_2 / r ** 2
print(np.mean(F), np.std(F))
0.7821065448652971 0.0056252380186608995
Plot the Monte-Carlo histogram using the option density=True and the Gaussian from the error propagation
In [27]:
# Define range of output values for plotting
xmin = 0.75
xmax = 0.82
# use Gaussian definition from above!
x = np.linspace(xmin, xmax, 1000)
y = gaussian(x, mean_f, sigma_f)
plt.hist(F, bins=50, range=[xmin, xmax], density=True)
plt.plot(x, y, color='red', lw=3)
plt.xlabel("Force [N]")
plt.ylabel("Relative Probability")
plt.xlim(xmin, xmax)
Out[27]:
(0.75, 0.82)
Repeat this with larger error bars:
In [28]:
# Define values to sample from
mean_m_1 = 40.e4
sigma_m_1 = 2.e4
mean_m_2 = 30.e4
sigma_m_2 = 10.e4
mean_r = 3.2
sigma_r = 1.0
# Define range of output values for plotting
xmin = -1.
xmax = 5.
# Define number of samples
N = 1000000
## STANDARD ERROR PROPAGATION
# Compute mean and error of force
mean_f = G * mean_m_1 * mean_m_2 / mean_r ** 2
sigma_f = mean_f * np.sqrt((sigma_m_1 / mean_m_2) ** 2
+ (sigma_m_2 / mean_m_2) ** 2
+ 2. * (sigma_r / mean_r) ** 2)
# Define Gaussian function
x = np.linspace(xmin, xmax, 1000)
y = gaussian(x, mean_f, sigma_f)
## MONTE-CARLO ERROR PROPAGATION
# Sample from initial values
m_1 = np.random.normal(mean_m_1, sigma_m_1, N)
m_2 = np.random.normal(mean_m_2, sigma_m_2, N)
r = np.random.normal(mean_r, sigma_r, N)
# Compute final values
F = G * m_1 * m_2 / r ** 2
## PLOTTING
plt.hist(F, bins=50, range=[xmin, xmax], density=True)
plt.plot(x, y, color='red', lw=3)
plt.xlabel("Force [N]")
plt.ylabel("Relative Probability")
plt.xlim(xmin, xmax)
Out[28]:
(-1.0, 5.0)
