schoener/DyLab_3D_MOT
1
0
Fork 0
Code Issues Pull Requests repo.project_board Releases Wiki Activity
b7b30c250f
DyLab_3D_MOT/src/coil_class.py

514 lines
20 KiB
Python
Raw Normal View History

initial
2021-10-01 14:25:09 +02:00
# -*- coding: utf-8 -*-
"""
Created on Tue Aug 17 15:17:44 2021
@author: Joschka
"""
import numpy as np
import matplotlib.pyplot as plt
Added round wire, checked class up to line 96
2021-10-06 21:58:19 +02:00
import logging as log
initial
2021-10-01 14:25:09 +02:00
from scipy import special as sp
from src import physical_constants as cs
Incomplete class modifications for round wire --> Must be continued
2021-10-07 17:49:17 +02:00
#TODO: Docstrings for every function
#TODO: Implement conventions
initial
2021-10-01 14:25:09 +02:00
class BCoil:
Incomplete class modifications for round wire --> Must be continued
2021-10-07 17:49:17 +02:00
def __init__(self, HH, distance, radius, layers, windings, wire_width, wire_height, insulation_thickness=0
, is_round = False, winding_offset = False):
"""
Creates object that represents a configuration of two coils, with symmetry axis = z
:param HH: HH = 1 --> current in same direction, homogeneous field, HH = -1 --> current in opposite direction, quadrupole field
:param distance: distance between center of coils
:param radius: radius to center of coil
:param layers: number of radial (horizontal) layers
:param windings: number of axial (vertical) layers
:param wire_width: width of conductive wire core
:param wire_height: height of conductive wire core
:param insulation_thickness: thickness of wire insulation
:param is_round: True --> Round wire, False --> rectangular wire
:param winding_offset: layer offset by half winding (especially for round wires)
"""
if not is_round:
if winding_offset:
log.warning('Is there a reason you want to wind a not round wire like that?')
Added round wire, checked class up to line 96
2021-10-06 21:58:19 +02:00
if is_round:
if wire_width != wire_height:
log.error('Wire is round but width != height')
Incomplete class modifications for round wire --> Must be continued
2021-10-07 17:49:17 +02:00
initial
2021-10-01 14:25:09 +02:00
self.HH = HH
self.distance = distance * 1e-3
self.radius = radius * 1e-3
self.layers = layers
self.windings = windings
self.wire_width = wire_width * 1e-3
self.wire_height = wire_height * 1e-3
Incomplete class modifications for round wire --> Must be continued
2021-10-07 17:49:17 +02:00
self.insulation_thickness = insulation_thickness * 1e-3
Added round wire, checked class up to line 96
2021-10-06 21:58:19 +02:00
self.is_round = is_round
Incomplete class modifications for round wire --> Must be continued
2021-10-07 17:49:17 +02:00
self.winding_offset = winding_offset
initial
2021-10-01 14:25:09 +02:00
def get_wire_area(self):
Incomplete class modifications for round wire --> Must be continued
2021-10-07 17:49:17 +02:00
"""
calculates wire area in m^2
:return: float
"""
Added round wire, checked class up to line 96
2021-10-06 21:58:19 +02:00
if self.is_round:
return np.pi * (self.wire_width/2)**2
initial
2021-10-01 14:25:09 +02:00
return self.wire_width * self.wire_height
def get_N(self):
Incomplete class modifications for round wire --> Must be continued
2021-10-07 17:49:17 +02:00
"""
Calculates number of windings
:return: int
"""
initial
2021-10-01 14:25:09 +02:00
return self.layers * self.windings
def get_wire_length(self):
Incomplete class modifications for round wire --> Must be continued
2021-10-07 17:49:17 +02:00
"""
:return: Approximate length of wire per coil (m)
"""
initial
2021-10-01 14:25:09 +02:00
return self.get_N() * 2 * self.radius * np.pi
Incomplete class modifications for round wire --> Must be continued
2021-10-07 17:49:17 +02:00
def get_tot_wire_height(self):
""" returns wire height incl. insulation"""
return self.wire_height + 2 * self.insulation_thickness
def get_tot_wire_width(self):
""" returns wire width incl. insulation"""
return self.wire_width + 2 * self.insulation_thickness
initial
2021-10-01 14:25:09 +02:00
def get_coil_height(self):
Incomplete class modifications for round wire --> Must be continued
2021-10-07 17:49:17 +02:00
if self.winding_offset:
return self.get_tot_wire_height() * (self.windings + 0.5)
return self.get_tot_wire_height() * self.windings
initial
2021-10-01 14:25:09 +02:00
def get_coil_width(self):
Incomplete class modifications for round wire --> Must be continued
2021-10-07 17:49:17 +02:00
if self.winding_offset:
if self.is_round:
return self.layers * self.get_tot_wire_width() -
#Todo : continue
return self.get_tot_wire_height() * self.windings
initial
2021-10-01 14:25:09 +02:00
def get_zmin(self):
return self.distance / 2 - self.get_coil_height() / 2
def get_zmax(self):
return self.distance / 2 + self.get_coil_height() / 2
def get_radius(self):
return self.radius
def get_R_inner(self):
return self.radius - self.get_coil_width() / 2
def get_R_outer(self):
return self.radius + self.get_coil_width() / 2
def set_R_outer(self, R_outer):
R_outer *= 1e-3
self.radius = R_outer - self.get_coil_width() / 2
def set_R_inner(self, R_inner):
R_inner *= 1e-3
self.radius = R_inner + self.get_coil_width() / 2
def set_d_min(self, d_min):
d_min *= 1e-3
self.distance = d_min + self.get_coil_height()
def set_d_max(self, d_max):
d_max *= 1e-3
self.distance = d_max - self.get_coil_height()
Incomplete class modifications for round wire --> Must be continued
2021-10-07 17:49:17 +02:00
def raster(self,raster_value):
"""
generates raster of flowing currents
:param raster_value: wire height/raster_value is distance between rastered points in one wire
:return: 2 dim array with z values in colums and R_values in rows
"""
initial
2021-10-01 14:25:09 +02:00
def print_info(self):
print(" ")
# print(f"{self.get_zmin()}")
# print(self.__name__)
print(
f"HH = {self.HH}, Distance = {self.distance * 1e3} mm, z_min = {self.get_zmin() * 1e3} mm, z_max = {self.get_zmax() * 1e3} mm")
print(
f"Radius = {self.radius * 1e3} mm, Radius_inner = {self.get_R_inner() * 1e3} mm, Radius_outer = {self.get_R_outer() * 1e3:.2f} mm")
print(
Incomplete class modifications for round wire --> Must be continued
2021-10-07 17:49:17 +02:00
f"layers = {self.layers}, windings = {self.windings}, wire_width = {self.wire_width * 1e3}, wire_height = {self.wire_height * 1e3} mm, round wire = {self.is_round} ")
initial
2021-10-01 14:25:09 +02:00
print(" ")
def print_basic_info(self):
print(
Incomplete class modifications for round wire --> Must be continued
2021-10-07 17:49:17 +02:00
f"HH = {self.HH}, Distance = {self.distance * 1e3} mm, Radius = {self.radius * 1e3} mm,layers = {self.layers}, windings = {self.windings}, round wire = {self.is_round}")
initial
2021-10-01 14:25:09 +02:00
def B_field_ideal_z(self, I_current, z_arg):
"""
Calculate B-field for ideal point like wires with R_radius and d_distance
HH = +1 --> Helmholtz configuration, HH = -1 --> Anti Helmholtz configuration
z_arg in mm
"""
z_SI = z_arg * 1e-3
N_windings = self.layers * self.windings
B = cs.mu_0 * N_windings * I_current / 2 * self.radius ** 2 * (
more stuff
2021-10-01 14:37:07 +02:00
1 / (self.radius ** 2 + (z_SI - self.distance / 2) ** 2) ** (3 / 2) + self.HH * 1 / (
self.radius ** 2 + (z_arg + self.distance / 2) ** 2) ** (3 / 2))
initial
2021-10-01 14:25:09 +02:00
B *= 1e4 # conversion Gauss
return B
Incomplete class modifications for round wire --> Must be continued
2021-10-07 17:49:17 +02:00
@staticmethod
initial
2021-10-01 14:25:09 +02:00
def k_sq(R_loop, z_loc, r, z):
""" Argument for elliptic integrals"""
return (4 * R_loop * r) / ((R_loop + r) ** 2 + (z - z_loc) ** 2)
Incomplete class modifications for round wire --> Must be continued
2021-10-07 17:49:17 +02:00
@staticmethod
initial
2021-10-01 14:25:09 +02:00
def B_z_loop(I_current, R_loop, z_loc, r, z):
"""calculate z-component of B-field at position r and z for each individual loop"""
B_z = 2e-7 * I_current * 1 / ((R_loop + r) ** 2 + (z - z_loc) ** 2) ** (1 / 2) * (
more stuff
2021-10-01 14:37:07 +02:00
sp.ellipk(BCoil.k_sq(R_loop, z_loc, r, z)) + sp.ellipe(BCoil.k_sq(R_loop, z_loc, r, z)) * (
R_loop ** 2 - r ** 2 - (z - z_loc) ** 2) / ((R_loop - r) ** 2 + (z - z_loc) ** 2))
initial
2021-10-01 14:25:09 +02:00
B_z *= 1e4 # conversion to gauss
return B_z
Incomplete class modifications for round wire --> Must be continued
2021-10-07 17:49:17 +02:00
@staticmethod
initial
2021-10-01 14:25:09 +02:00
def B_r_loop(I_current, R_loop, z_loc, r, z):
"""calculate r-component of B-field at position r and z for each individual loop"""
B_r = 2e-7 * I_current / r * (z - z_loc) / ((R_loop + r) ** 2 + (z - z_loc) ** 2) ** (1 / 2) * (
more stuff
2021-10-01 14:37:07 +02:00
-sp.ellipk(BCoil.k_sq(R_loop, z_loc, r, z)) + sp.ellipe(BCoil.k_sq(R_loop, z_loc, r, z)) * (
R_loop ** 2 + r ** 2 + (z - z_loc) ** 2) / ((R_loop - r) ** 2 + (z - z_loc) ** 2))
initial
2021-10-01 14:25:09 +02:00
B_r *= 1e4 # conversion to gauss
return B_r
def B_multiple(self, I_current, x, z, raster=10):
"""
Returns B_z along z-axis and B_x along x-axis,
HH = +1 --> Helmholtz configuration, HH = -1 --> Anti Helmholtz configuration
"""
z_start = self.get_zmin() + self.wire_height / 2 # (distance_coils/2 - windings * wire_height/2 + wire_height/2)*1e-3
R_start = self.get_R_inner() + self.wire_width / 2 # (R_inner + wire_width/2 )*1e-3
# Convert axis to SI units
x_SI = x * 1e-3
z_SI = z * 1e-3
if x[0] <= 0:
x_neg = np.abs([el for el in x_SI if el < 0])
x_pos = np.array([el for el in x_SI if el > 0])
x_zero = np.array([el for el in x_SI if el == 0])
B_z = np.zeros(len(z))
B_x_pos = np.zeros(len(x_pos))
B_x_neg = np.zeros(len(x_neg))
B_x_zero = x_zero # If there is a zero in the x-array it is assured that the x component of the field at this point is zero
# dividing each wire into 10 x 10 raster
I_current /= raster ** 2
for xx in range(0, self.layers):
for zz in range(0, self.windings):
for xx_in in range(0, raster):
for zz_in in range(0, raster):
z_pos = z_start + zz * (
more stuff
2021-10-01 14:37:07 +02:00
self.wire_height + self.windings_spacing) - self.wire_height / 2 + zz_in * self.wire_height / (
raster - 1)
initial
2021-10-01 14:25:09 +02:00
R_pos = R_start + xx * (
more stuff
2021-10-01 14:37:07 +02:00
self.wire_width + self.layers_spacing) - self.wire_width / 2 + xx_in * self.wire_width / (
raster - 1)
Incomplete class modifications for round wire --> Must be continued
2021-10-07 17:49:17 +02:00
#TODO: fix z_pos, R_pos, either by raster or by spacing
initial
2021-10-01 14:25:09 +02:00
B_z += BCoil.B_z_loop(I_current, R_pos, z_pos, 0, z_SI) + BCoil.B_z_loop(self.HH * I_current,
R_pos, -z_pos, 0, z_SI)
B_x_pos += BCoil.B_r_loop(I_current, R_pos, z_pos, x_pos, 0) + BCoil.B_r_loop(
self.HH * I_current, R_pos, -z_pos, x_pos, 0)
B_x_neg += BCoil.B_r_loop(I_current, R_pos, z_pos, x_neg, 0) + BCoil.B_r_loop(
self.HH * I_current, R_pos, -z_pos, x_neg, 0)
B_x_neg *= -1 # B_x_neg is pointing in opposite x_direction
B_x = np.concatenate((B_x_neg, B_x_zero, B_x_pos))
return B_z, B_x
def B_tot_along_axis(self, I_current, x, z, raster=10):
"""
return B_tot_z, B_tot_x
Returns B_tot_z along z-axis and B_tot_x along x-axis,
HH = +1 --> Helmholtz configuration, HH = -1 --> Anti Helmholtz configuration
"""
z_start = self.get_zmin() + self.wire_height / 2 # (distance_coils/2 - windings * wire_height/2 + wire_height/2)*1e-3
R_start = self.get_R_inner() + self.wire_width / 2 # (R_inner + wire_width/2 )*1e-3
# Convert axis to SI units
x_SI = x * 1e-3
z_SI = z * 1e-3
if x[0] <= 0:
x_neg = np.abs([el for el in x_SI if el < 0])
x_pos = np.array([el for el in x_SI if el > 0])
x_zero = np.array([el for el in x_SI if el == 0])
B_tot_z = np.zeros(len(z))
B_x_pos = np.zeros(len(x_pos))
B_x_neg = np.zeros(len(x_neg))
B_x_zero = x_zero # If there is a zero in the x-array it is assured that the x component of the field at this point is zero
B_z_x = np.zeros(len(x_SI))
# dividing each wire into 10 x 10 raster
I_current /= raster ** 2
for xx in range(0, self.layers):
for zz in range(0, self.windings):
for xx_in in range(0, raster):
for zz_in in range(0, raster):
z_pos = z_start + zz * (
more stuff
2021-10-01 14:37:07 +02:00
self.wire_height + self.windings_spacing) - self.wire_height / 2 + zz_in * self.wire_height / (
raster - 1)
initial
2021-10-01 14:25:09 +02:00
R_pos = R_start + xx * (
more stuff
2021-10-01 14:37:07 +02:00
self.wire_width + self.layers_spacing) - self.wire_width / 2 + xx_in * self.wire_width / (
raster - 1)
Incomplete class modifications for round wire --> Must be continued
2021-10-07 17:49:17 +02:00
# TODO: fix z_pos, R_pos, either by raster or by spacing
initial
2021-10-01 14:25:09 +02:00
# z-field along z-axis (x-Field always zero)
B_tot_z += BCoil.B_z_loop(I_current, R_pos, z_pos, 0, z_SI) + BCoil.B_z_loop(
self.HH * I_current, R_pos, -z_pos, 0, z_SI)
# x-field along x-axis
B_x_pos += BCoil.B_r_loop(I_current, R_pos, z_pos, x_pos, 0) + BCoil.B_r_loop(
self.HH * I_current, R_pos, -z_pos, x_pos, 0)
B_x_neg += BCoil.B_r_loop(I_current, R_pos, z_pos, x_neg, 0) + BCoil.B_r_loop(
self.HH * I_current, R_pos, -z_pos, x_neg, 0)
# B_z along x-axis:
B_z_x += BCoil.B_z_loop(I_current, R_pos, z_pos, x_SI, 0) + BCoil.B_z_loop(self.HH * I_current,
R_pos, -z_pos, x_SI,
0)
B_x_neg *= -1 # B_x_neg is pointing in opposite x_direction
B_x_x = np.zeros(len(x_SI))
B_x_x = np.concatenate((B_x_neg, B_x_zero, B_x_pos))
B_tot_x = np.sqrt(B_x_x ** 2 + B_z_x ** 2)
return B_tot_z, B_tot_x
def B_multiple_3d(self, I_current, x, z, raster=4):
z_start = self.get_zmin() + self.wire_height / 2
R_start = self.get_R_inner() + self.wire_width / 2
x_SI = x * 1e-3
z_SI = z * 1e-3
B = np.zeros([len(z_SI), len(x_SI), 2])
I_current /= raster ** 2
for el in range(0, len(x_SI)):
for xx in range(0, self.layers):
for zz in range(0, self.windings):
for xx_in in range(0, raster):
for zz_in in range(0, raster):
z_pos = z_start + zz * (
more stuff
2021-10-01 14:37:07 +02:00
self.wire_height + self.windings_spacing) - self.wire_height / 2 + zz_in * self.wire_height / (
raster - 1)
initial
2021-10-01 14:25:09 +02:00
R_pos = R_start + xx * (
more stuff
2021-10-01 14:37:07 +02:00
self.wire_width + self.layers_spacing) - self.wire_width / 2 + xx_in * self.wire_width / (
raster - 1)
Incomplete class modifications for round wire --> Must be continued
2021-10-07 17:49:17 +02:00
# TODO: fix z_pos, R_pos, either by raster or by spacing
initial
2021-10-01 14:25:09 +02:00
# compute z-value of field
B[:, el, 1] += BCoil.B_z_loop(I_current, R_pos, z_pos, np.abs(x_SI[el]),
z_SI) + BCoil.B_z_loop(self.HH * I_current, R_pos, -z_pos,
np.abs(x_SI[el]), z_SI)
# compute x-value
if x_SI[el] < 0:
B[:, el, 0] -= BCoil.B_r_loop(I_current, R_pos, z_pos, -x_SI[el],
z_SI) + BCoil.B_r_loop(self.HH * I_current, R_pos, -z_pos,
-x_SI[el], z_SI)
elif x_SI[el] > 0:
B[:, el, 0] += BCoil.B_r_loop(I_current, R_pos, z_pos, x_SI[el], z_SI) + BCoil.B_r_loop(
self.HH * I_current, R_pos, -z_pos, x_SI[el], z_SI)
elif x_SI[el] == 0:
B[:, el, 0] = 0
return B
more stuff
2021-10-01 14:37:07 +02:00
@staticmethod
initial
2021-10-01 14:25:09 +02:00
def B_tot_3d(B):
return np.sqrt(B[:, :, 0] ** 2 + B[:, :, 1] ** 2)
def plot_3d(self, I_current, x_lim, z_lim):
x = np.arange(-x_lim, x_lim, 5)
z = np.arange(-z_lim, z_lim, 5)
x_m, z_m = np.meshgrid(x, z)
B = self.B_multiple_3d(I_current, x, z, raster=2)
B_tot = BCoil.B_tot_3d(B)
for xx in range(0, len(x)):
for zz in range(0, len(z)):
if B_tot[xx, zz] > 15:
B[xx, zz, :] /= B_tot[xx, zz] / 15
plt.figure(33)
plt.quiver(x_m, z_m, B[:, :, 0], B[:, :, 1])
plt.xlabel("x-axis [mm]")
plt.ylabel("z-axis [mm]")
plt.show()
more stuff
2021-10-01 14:37:07 +02:00
def curv(B_field, z):
initial
2021-10-01 14:25:09 +02:00
return np.gradient(np.gradient(B_field, z), z) * 1e2
def Bgrad(B_field, z):
return np.gradient(B_field, z) * 1e1
def Bz_plot_HH(self, I_current, x, z):
B_z, B_x = self.B_multiple(I_current, x, z)
B_z_curvature = np.gradient(np.gradient(B_z, z), z)
plt.figure(100, figsize=(13, 10))
# plt.rcParams.update({'font.size': 15})
plt.suptitle("Helmholtz coil field B_z along z-axis")
# Field plot
##########################
plt.subplot(2, 1, 1)
plt.plot(z, B_z, linestyle="solid", label=r"$B_z$: Result via elliptic integrals")
# plt.xlim(-0.01,0.01)
plt.title("B-field")
plt.ylabel(r"$B_z$ [G]")
plt.xlabel("z-axis [mm]")
plt.legend()
plt.subplot(2, 1, 2)
plt.plot(z, B_z_curvature, linestyle="solid", label=r"$\nabla_z^2 B_z$: Result via elliptic integrals")
plt.ylabel(r"$\nabla_z^2 B_z [G/cm^2]$")
plt.xlabel("z-axis [mm]")
plt.xlim(-10, 10)
plt.title("Curvature of B-field")
plt.legend(loc='lower right')
plt.show()
def Bz_plot_HH_comp(self, Coil2, I_current, x, z):
B_z, B_x = self.B_multiple(I_current, x, z)
B_z_2, B_x_2 = Coil2.B_multiple(I_current, x, z)
B_z_curvature = np.gradient(np.gradient(B_z, z), z) * 1e2
B_z_curvature_2 = np.gradient(np.gradient(B_z_2, z), z) * 1e2
plt.figure(100, figsize=(13, 10))
# plt.rcParams.update({'font.size': 15})
plt.suptitle("Helmholtz coil field B_z along z-axis")
# Field plot
##########################
plt.subplot(2, 1, 1)
plt.plot(z, B_z, linestyle="solid", label=r"$B_z$")
plt.plot(z, B_z_2, linestyle="solid", label=r"$B_{z2}$")
# plt.xlim(-0.01,0.01)
plt.title("B-field")
plt.ylabel(r"$B_z$ [G]")
plt.xlabel("z-axis [mm]")
plt.legend()
plt.subplot(2, 1, 2)
plt.plot(z, B_z_curvature, linestyle="solid", label=r"$\nabla_z^2 B_z$")
plt.plot(z, B_z_curvature_2, linestyle="solid", label=r"$\nabla_z^2 B_{z2}$")
plt.ylabel(r"$\nabla_z^2 B_z [G/cm^2]$")
plt.xlabel("z-axis [mm]")
plt.xlim(-10, 10)
plt.title("Curvature of B-field")
plt.legend(loc='lower right')
plt.show()
def cooling(self, I_current, T):
"""
Print current density and power
Parameters
----------
I_current : current in A
T : temperature in degree Celsius
Returns
-------
None.
"""
j_dens = I_current / self.get_wire_area() * 1e-6
Power = self.resistance(T) * I_current ** 2
print(f"current density = {j_dens} A/mm^2")
print(f"Power = {Power} W")
def power(self, I_current, T):
P = self.resistance(T) * I_current ** 2
return P
def inductivity(self):
return cs.mu_0 * (self.layers * self.windings) ** 2 * self.radius ** 2 * np.pi / self.get_coil_height()
def induct_perry(self):
"""
Calculates inductance of rectangular coil via empirical formular by perry
"""
L = 4 * np.pi * (self.windings * self.layers) ** 2 * (1e2 * self.radius) ** 2 / (
more stuff
2021-10-01 14:37:07 +02:00
0.2317 * 100 * self.radius + 0.44 * 100 * self.get_coil_height() + 0.39 * 100 * self.get_coil_width())
initial
2021-10-01 14:25:09 +02:00
return L * 1e-9
Incomplete class modifications for round wire --> Must be continued
2021-10-07 17:49:17 +02:00
@staticmethod
def resistivity_copper(T):
R_ref = cs.rho_copper_20
T_ref = 20 # degree celsius
alpha = 0.00393
R = R_ref * (1 + alpha * (T - T_ref))
return R
initial
2021-10-01 14:25:09 +02:00
def resistance(self, T):
"""
Incomplete class modifications for round wire --> Must be continued
2021-10-07 17:49:17 +02:00
Calculates resistance of one coil of configuration
initial
2021-10-01 14:25:09 +02:00
Parameters
----------
T : double
temperature in degree Celsius
Returns
-------
double
resistance in Ohm
"""
Incomplete class modifications for round wire --> Must be continued
2021-10-07 17:49:17 +02:00
return cf.resistivity_copper(T) * self.get_wire_length() / self.get_wire_area()
initial
2021-10-01 14:25:09 +02:00
Incomplete class modifications for round wire --> Must be continued
2021-10-07 17:49:17 +02:00
HH_Coil = BCoil(HH=1, distance=54, radius=48, layers=8, windings=8, wire_height=0.4,wire_width=0.4, insulation_thickness = 0.1 ,is_round= False, winding_offset= True)
initial
2021-10-01 14:25:09 +02:00
HH_Coil.set_R_outer(49.3)
HH_Coil.set_d_min(49.8)
Incomplete class modifications for round wire --> Must be continued
2021-10-07 17:49:17 +02:00
HH_Coil.print_info()
Added round wire, checked class up to line 96
2021-10-06 21:58:19 +02:00
print(HH_Coil.get_wire_area()*1e6)
print(HH_Coil.get_wire_length())
print(HH_Coil.get_coil_height())
Reference in New Issue Copy Permalink