import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
09.01
np.logspace(-20, -10, 11)
np.array(10*[2.])
#
2*np.ones(10)
np.zeros(5, dtype = np.float32)
09.02
#np.random.poisson?
x = np.random.poisson(100,10)
print(x, x[1:], x[:-1])
dx = x[1:] - x[:-1]
print(dx)
09.03
year, temp = np.loadtxt('../data/munich_temperatures_average_with_bad_data.txt',
unpack=True)
print(temp.min(), temp.max())
temp_good = temp[(temp!=-99) & (temp!=99)]
print(len(temp)-len(temp_good), 'elements were removed')
10.01
date, temperature = np.loadtxt('../data/munich_temperatures_average_with_bad_data.txt', unpack=True)
keep = np.abs(temperature) < 99
date = date[keep]
temperature = temperature[keep]
plt.plot(date, temperature, marker='.', linewidth=0)
frac = date%1
plt.plot(frac, temperature, marker='.', linewidth=0)
10.02
# define normalized Gaussian with mean 'mu' and sigma 'sig'
def gaussian(x, mu, sig):
return 1./(np.sqrt(2.*np.pi)*sig) * np.exp(-0.5*((x-mu)/sig)**2)
# Gaussian with mean 10 and sigma sqrt(10)
x = np.linspace(-5, 30, 100)
d = np.random.normal(10, np.sqrt(10), 10000)
n_bins = int(10000**(1./3.))
h = plt.hist(d, bins=n_bins, density=True)
plt.plot(x, gaussian(x, 10, np.sqrt(10)), color='red', linewidth=3)
# I plot here the sample from the Poisson distribution together with the Gaussian
N = 10000
mu = 12
x = np.linspace(-5, 50, 100)
d2 = np.random.poisson(mu, N)
h2 = plt.hist(d2, bins=np.linspace(-5, 50, 56), density=True)
plt.plot(x, gaussian(x, mu, np.sqrt(mu)), color='red', linewidth=3)
plt.xlim([-10,40])
What you observe here is the central limit theorem in action. For large mean values mu, the Poisson distribution approaches a Gaussian distribution (it becomes more and more symmetric). Also, if we draw $N$ times from a Poisson distribution with mean $\lambda$, the sum of the draws will be distributed approximately normal with mean $\mu=N\lambda$ and variance $\sigma^2=N\lambda$ (since $\lambda$ is the variance of the Poisson distribution).
10.03
First, it is useful to realise that ($\vec{r}=x\vec{e}_\mathrm{x} + y\vec{e}_\mathrm{y} + z\vec{e}_\mathrm{z}$) $$\vec{e}_\varphi = \vec{e}_\mathrm{z} \times \vec{e}_\mathrm{r} = \vec{e}_\mathrm{z} \times \frac{\vec{r}}{|\vec{r}|} = -y\vec{e}_\mathrm{x} + x\vec{e}_\mathrm{y}.$$ For the magnetic field of each wire, we then have $$\vec{B} \propto -\frac{y}{r^2}\vec{e}_\mathrm{x} + \frac{x}{r^2}\vec{e}_\mathrm{y}$$ with the proportionality constant given by $\mu_0 I/2\pi$ and $r^2=|\vec{r}|^2$.
We then apply a shift along the y-axis by an amount $a$ for each wire and add the contributions of both wires to obtain the total magnetic field strength.
x = np.arange(-10, 10, 0.1)
y = np.arange(-10, 10, 0.1)
xx, yy = np.meshgrid(x, y)
a = 5
r1_sq = xx**2 + (yy-a)**2
r2_sq = xx**2 + (yy+a)**2
b1x, b1y = -(yy-a)/r1_sq, +xx/r1_sq
b2x, b2y = -(yy+a)/r2_sq, +xx/r2_sq
plt.axes().set_aspect('equal')
plt.streamplot(x, y, b1x+b2x, b1y+b2y, color='k')
We have to integrate the B-fields $\mathrm{d}\vec{B}$ produced by wire elements $\mathrm{d}\vec{l}$. We therefore define the polar angle $\phi$ in the xz-plane such the coil of radius $a$ has a position vector $\vec{l}=a\cos\phi\vec{e}_\mathrm{x} + 0\vec{e}_\mathrm{y} + a\sin\phi\vec{e}_\mathrm{z}$. The distance vector $\vec{r}$ from any position in space to the coil is then $$\vec{r}=(x-a\cos\phi)\vec{e}_\mathrm{x} + y\vec{e}_\mathrm{y} + (z-a\sin\phi)\vec{e}_\mathrm{z}.$$ The absolute distance from any point in space to one wire element at a certain $\phi$ is thus $r^2=|\vec{r}|^2=(x-a\cos\phi)^2 + y^2 + (z-a\sin\phi)^2$. In the Biot-Savart law, the wire element $\mathrm{d}\vec{l}$ points along the path of the wire such that $$\mathrm{d}\vec{l} = \left[\vec{l}\times\vec{e}_\mathrm{y}\right]\mathrm{d}\phi=\left[-a\sin\phi\vec{e}_\mathrm{x} + a\cos\phi\vec{e}_\mathrm{z}\right]\mathrm{d}\phi.$$ For the cross product $\mathrm{d}\vec{l}\times\vec{r}$ in the magnetic field $\mathrm{d}\vec{B}$, we then have $$\mathrm{d}\vec{l}\times\vec{r} = \left[-ay\cos\phi\vec{e}_\mathrm{x} - a(a-x\cos\phi)\vec{e}_\mathrm{y} - ay\sin\phi\vec{e}_\mathrm{z}\right]\mathrm{d}\phi.$$ What is left now is to sum all the contributions from the wire elememts for all possible $\phi$-values from $0$ to $2\pi$ and plot the resulting fieldlines in the xy-plane.
a = 5 # radius of the coil
n = 90 # number of steps for integration along the wire
x = np.arange(-10, 10, 0.05)
y = np.arange(-10, 10, 0.05)
dphi = 2*np.pi/n
# init. the grid and zero values for the B-field
xx, yy = np.meshgrid(x, y)
bx = 0*xx
by = 0*yy
for phi in np.arange(0, 2*np.pi, dphi):
cos = np.cos(phi)
sin = np.sin(phi)
# separation coil-element to observer
rp = np.sqrt((xx-a*cos)**2 + yy**2 + (0-a*sin)**2)
# magnetic field components in x and y direction
bx += -a*cos*yy*dphi/rp**3
by += -a*(a-cos*xx)*dphi/rp**3
plt.axes().set_aspect('equal')
plt.streamplot(x, y, bx, by, color='k')
11.01
import glob
import os
for file in glob.glob('*'):
print(file, str(int(os.path.getsize(file)/1024+0.5))+'kB')
12.01 and 12.02
def date2str(li):
[year, month, day, hour, minute, second] = li
date = "{0:04d}-{1:02d}-{2:02d}".format(year, month, day)
time = "{0:02d}:{1:02d}:{2:02d}".format(hour, minute, second)
s = "{0} {1}".format(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([2019, 2, 27, 15, 37, 12])))
print(date2str(str2date('2019-02-27 15:37:12')))
# Define constant(s) in SI units
G = 6.67384e-11
# Define values to sample from
mean_m_1 = 40.0e4
sigma_m_1 = 0.05e4
mean_m_2 = 30.0e4
sigma_m_2 = 0.1e4
mean_r = 3.2
sigma_r = 0.01
Standard error propagation gives
# 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
+ (2. * sigma_r / mean_r)**2)
print(mean_f, sigma_f)
Generate random samples for $M_1$, $M_2$ and $r$ and calculate the force
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))
Plot the Monte-Carlo histogram using the option density=True and the Gaussian from the error propagation
xmin = 0.75
xmax = 0.82
# use Gaussian definition from above!
#def gaussian(x, mu, sig):
# return 1./(np.sqrt(2.*np.pi)*sig) * np.exp(-0.5*((x-mu)/sig)**2)
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("Probability density")
plt.xlim(xmin, xmax)
Repeat this with larger uncertainties:
# Define values to sample from
mean_m_1 = 40.0e4
sigma_m_1 = 2.0e4
mean_m_2 = 30.0e4
sigma_m_2 = 10.0e4
mean_r = 3.2
sigma_r = 1.0
# Define range of output values for plotting
xmin = -1.0
xmax = 5.0
# 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_1)**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("Probability density")
plt.xlim(xmin, xmax)
Standard error propagation assumes a linear expansion of a function. So for small uncertainties, the linear expansion is likely a good approximation. For larger uncertainties or highly non-linear functions, this must not be the case anymore. Hence, standard error propagation breaks down in such cases and the Monte-Carlo method is the method of choice.