Homec4science
Diffusion JuNoTE_git_repository (master) Browse

R9442/Interactive Cells/Induction_source.py
e4482c463143master

Last Change
Raw File

Induction_source.py
View Options

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
from matplotlib.widgets import Slider
from mpl_toolkits.mplot3d import Axes3D
def update_position_record(record, update, N_max=100):
''' Update record of positions.
record (list of lists): list of the records for every position to be updated.
update (list): list of the updates for every dimension.
N_max (int): Maximum number of elements to hold in records.'''
for i in range(len(update)):
if len(record[i]) >= N_max:
record[i].pop(0)
record[i].append(update[i])
return record
def get_overlap_area(coords, widths):
x1, y1, x2, y2 = coords
wx1, wy1, wx2, wy2 = widths
if x1+wx1 > x2 :
x_width = min([wx2, abs((x1 + wx1) - x2)])
y_width = min([wy2, abs((y1 + wy1) - y2), abs((y2+wy2) - y1)])
return x_width, y_width, x_width * y_width
else:
return 0, 0, 0
def plot_simple_coils(big_x, big_y, small_x, small_y, record, fig=None):
constant_B_magnitude = 1
big_width = 6
big_height = 8
small_width = 3
small_height = 4
coords = (big_x, big_y, small_x, small_y)
fig = fig or plt.figure(figsize=(15, 3))
ax0 = fig.add_subplot(131)
ax1 = fig.add_subplot(132)
ax2 = fig.add_subplot(133)
overlap_x, overlap_y, overlap_area = get_overlap_area(
coords, (big_width, big_height, small_width, small_height))
record = update_position_record(
record, (*coords, overlap_area * constant_B_magnitude))
ax0.set_xlim(0, 10)
ax0.set_ylim(0, 10)
# Draw rectangles
big_rect = patches.Rectangle((big_x, big_y),
big_width,
big_height,
linewidth=1,
edgecolor='r',
facecolor='none')
small_rect = patches.Rectangle((small_x, small_y),
small_width,
small_height,
linewidth=1,
edgecolor='b',
facecolor='none')
shaded_rect = patches.Rectangle((small_x, max(small_y, big_y)),
overlap_x,
overlap_y,
facecolor='k',
alpha=0.2)
ax0.add_artist(big_rect)
ax0.add_artist(small_rect)
ax0.add_artist(shaded_rect)
# Draw Phi
#ax1.clear()
line, = ax1.plot(record[-1], color='r')
#Draw FEM
#ax2.clear()
if len(record[-1]) >=2:
ax2.plot(-np.gradient(np.array(record[-1]))[:-1], color='b')
else:
ax2.plot(np.zeros((2)), color='b')
[a.set_xlabel('t (au)') for a in [ax1, ax2]]
ax0.set_xlabel('x'); ax0.set_ylabel('y')
ax1.set_ylabel('Magnetic flux')
ax2.set_ylabel('FEM')
plt.show()
plt.close(fig)
return big_rect, small_rect, shaded_rect
class Ring():
def __init__(self, center, radius, orientation):
self.center_x = center[0]
self.center_y = center[1]
self.center_z = center[2]
self.center = center
self.radius = radius
self.axis_th = orientation[0]
self.axis_ph = orientation[1]
self.norm_vector = np.array([
np.sin(self.axis_th) * np.cos(self.axis_ph) - self.center_x,
np.sin(self.axis_th) * np.sin(self.axis_ph) - self.center_y,
np.cos(self.axis_th) - self.center_z
])
self.area = np.pi * self.radius**2
self.plane_vector = np.array([
np.sin(self.axis_th - np.pi / 2) * np.cos(self.axis_ph) -
self.center_x,
np.sin(self.axis_th - np.pi / 2) * np.sin(self.axis_ph) -
self.center_y,
np.cos(self.axis_th - np.pi / 2) - self.center_z
])
self.m = np.array([0, 0, 1]) * 1e5
self.m_r0 = np.array([2, 0, 0])
self.flux_record = [self.magnetic_flux()]
def dipole_field(self,r, r0=[0,0,0],m=np.array([0, 1, 0]), c=1):
"""Magnetic field of a dipole with moment m along the y axis (default).
r (1d array): [x, y, z]"""
# Using Gaussian units & c=1
# mu_0/4 pi = 1/c = 1
r = np.array(r)
r0 = np.array(r0)
m = np.array(m)
r=np.subtract(r.T, r0).T
r_norm = np.sqrt(np.einsum('i...,i...', r, r))
r_unit = r / r_norm
B = (1 / r_norm**3) * (3 * r_unit *
(np.einsum('i...,i...', r_unit, m))) - np.einsum(
'...,i->i...',1 / r_norm**3, m)
return B
def update_m_r0(self, new_x, new_y, new_z, m_mag, m_th, m_ph):
self.m_r0 = np.array([new_x, new_y, new_z])
self.m = 1e5 * m_mag * np.array([
np.sin(m_th) * np.cos(m_ph),
np.sin(m_th) * np.sin(m_ph),
np.cos(m_th)
])
def plot(self, ax):
t = np.linspace(0, 2 * np.pi, 100)
cross = np.cross(self.norm_vector, self.plane_vector)
coords = self.radius * np.einsum(
'i,j', np.cos(t), self.plane_vector) + self.radius * np.einsum(
'i,j', np.sin(t), cross) + self.center
ax.plot(coords[:, 0], coords[:, 1], coords[:, 2])
ax.plot([self.center_x, self.norm_vector[0]],
[self.center_y, self.norm_vector[1]],
[self.center_z, self.norm_vector[2]])
ax.plot([self.center_x, self.plane_vector[0]],
[self.center_y, self.plane_vector[1]],
[self.center_z, self.plane_vector[2]])
# Plot monopole
ax.quiver(
self.m_r0[0],
self.m_r0[1],
self.m_r0[2], # <-- starting point of vector
1e-5*self.m[0],
1e-5*self.m[1],
1e-5*self.m[2], # <-- directions of vector
color='red',
alpha=.8,
lw=3,
)
ax.plot([self.m_r0[0]], [self.m_r0[1]], [self.m_r0[2]],
'ob',
label='Dipole')
ax.legend(loc=0)
def plot_flux(self, axes, flux_record):
self.line_flux, = axes[0].plot(flux_record)
if len(flux_record) >= 3:
self.line_fem, = axes[1].plot(
-np.gradient(np.array(flux_record))[:-1])
else:
self.line_fem, = axes[1].plot(np.zeros((3)))
def update_plots(self, new_x, new_y, new_z, m_mag, m_th, m_ph, fig=None):
Bfunc = self.dipole_field
fig = fig or plt.figure(figsize=(20, 5))
ax0 = fig.add_subplot(131, projection='3d')
ax1 = fig.add_subplot(132)
ax2 = fig.add_subplot(133)
ax0.set_xlim(-3, 3)
ax0.set_ylim(-3, 3)
ax0.set_zlim(-3, 3)
ax0.set_xlabel('x')
ax0.set_ylabel('y')
ax0.set_zlabel('z')
ax0.set_title('System')
ax1.set_title('Magnetic Flux')
ax2.set_title('EMF')
ax = [ax0, ax1, ax2]
#ax0, ax1, ax2 = ax
self.update_m_r0(new_x, new_y, new_z, m_mag, m_th, m_ph)
self.flux_record.append(self.magnetic_flux())
self.plot_flux(ax[1:], self.flux_record)
self.plot(ax[0])
self.line_flux.set_data(range(len(self.flux_record)), self.flux_record)
self.line_fem.set_data(range(len(self.flux_record)),
-np.gradient(self.flux_record))
fig.canvas.draw_idle()
def diff_magnetic_flux(self, z, y):
Bfunc = self.dipole_field
B = Bfunc([self.center_x, y, z], r0=self.m_r0, m=self.m)
dPhi = B.dot(self.norm_vector)
return dPhi
def magnetic_flux(self):
Bfunc = self.dipole_field
ybounds = (-self.radius, self.radius)
zbounds = (lambda y: -np.sqrt(self.radius**2 - y**2),
lambda y: np.sqrt(self.radius**2 - y**2))
from scipy.integrate import dblquad
integral, err = dblquad(
lambda z, y: self.diff_magnetic_flux(z, y), ybounds[0],
ybounds[1], zbounds[0], zbounds[1])
return integral
c4science · Help