817 lines
30 KiB
Python
817 lines
30 KiB
Python
# -*- coding: utf-8 -*-
|
|
"""
|
|
Created on Tue Aug 17 15:17:44 2021
|
|
|
|
@author: Joschka
|
|
"""
|
|
import numpy as np
|
|
import matplotlib.pyplot as plt
|
|
import logging as log
|
|
from scipy import special as sp
|
|
|
|
import matplotlib
|
|
# matplotlib.use('Qt5Agg')
|
|
# matplotlib.use('Agg')
|
|
# get_ipython().run_line_magic('matplotlib', 'qt')
|
|
# get_ipython().run_line_magic('matplotlib', 'inline')
|
|
# import time
|
|
|
|
from src import physical_constants as cs
|
|
|
|
# logger = log.getLogger('example')
|
|
# log.setLevel(log.info)
|
|
log.basicConfig(level=log.WARNING, format='%(message)s')
|
|
|
|
|
|
# TODO: Docstrings for every function
|
|
|
|
# TODO: Implement conventions
|
|
|
|
class BCoil:
|
|
def __init__(self, HH, distance, radius, layers, windings, wire_width, wire_height, insulation_thickness=0.
|
|
, is_round=False, winding_scheme=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 (radial increase of thickness not diameter!)
|
|
:param is_round: True --> Round wire, False --> rectangular wire
|
|
:param winding_scheme: 0: standard, layer on layer, 1: with offset (for round wires primarily), 2: like 1, but alternatingly 8 --> 7 windings per layer
|
|
"""
|
|
if not is_round:
|
|
if winding_scheme == 1 or winding_scheme == 2:
|
|
log.warning('Is there a reason you want to wind a not round wire like that?')
|
|
if is_round:
|
|
if wire_width != wire_height:
|
|
log.error('Wire is round but width != height')
|
|
raise ValueError("Wire is round but width != height")
|
|
|
|
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
|
|
self.insulation_thickness = insulation_thickness * 1e-3
|
|
self.is_round = is_round
|
|
self.winding_scheme = winding_scheme
|
|
|
|
# Standard get/set functions
|
|
@property
|
|
def get_HH(self):
|
|
return self.HH
|
|
|
|
@property
|
|
def get_distance(self):
|
|
return self.distance
|
|
|
|
@property
|
|
def get_radius(self):
|
|
return self.radius
|
|
|
|
@property
|
|
def get_layers(self):
|
|
return self.layers
|
|
|
|
@property
|
|
def get_windings(self):
|
|
return self.windings
|
|
|
|
@property
|
|
def get_wire_width(self):
|
|
return wire_width
|
|
|
|
@property
|
|
def get_wire_height(self):
|
|
return self.wire_height
|
|
|
|
@property
|
|
def get_insulation_thickness(self):
|
|
return self.insulation_thickness
|
|
|
|
@property
|
|
def get_is_round(self):
|
|
return self.is_round
|
|
|
|
@property
|
|
def get_winding_scheme(self):
|
|
return self.winding_scheme
|
|
|
|
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_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()
|
|
|
|
def get_wire_area(self):
|
|
"""
|
|
calculates wire area in m^2
|
|
:return: float
|
|
"""
|
|
if self.is_round:
|
|
return np.pi * (self.wire_width / 2) ** 2
|
|
return self.wire_width * self.wire_height
|
|
|
|
def get_N(self):
|
|
"""
|
|
Calculates number of windings
|
|
"""
|
|
N = self.layers * self.windings
|
|
log.debug(f"N = {N}")
|
|
if self.winding_scheme == 2:
|
|
N -= self.layers // 2
|
|
log.debug(f"N = {N}")
|
|
return N
|
|
|
|
def get_wire_length(self):
|
|
"""
|
|
:return: Approximate length of wire per coil (m)
|
|
"""
|
|
return self.get_N() * 2 * self.radius * np.pi
|
|
|
|
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
|
|
|
|
def get_coil_height(self):
|
|
if self.winding_scheme == 1:
|
|
return self.get_tot_wire_height() * (self.windings + 0.5)
|
|
return self.get_tot_wire_height() * self.windings
|
|
|
|
def get_coil_width(self):
|
|
if self.winding_scheme == 1 or self.winding_scheme == 2:
|
|
if self.is_round:
|
|
log.info("Coil width: Be aware of the fact that this is an idealized situation of coil winding (slope "
|
|
"of windings changes each layer)")
|
|
return self.layers * self.get_tot_wire_width() - (self.layers - 1) * (
|
|
2 - np.sqrt(3)) * self.get_tot_wire_width() / 2 # width is reduced due to winding offset
|
|
return self.get_tot_wire_width() * self.layers
|
|
|
|
def winding_raster(self):
|
|
"""
|
|
generates raster of flowing currents
|
|
:param raster_value: wire height/raster_value is distance between rastered points in one wire
|
|
:return: 2 dim array [[z1,R1], [z2,R2], ...]
|
|
"""
|
|
outer_raster = np.zeros((self.get_N(), 2))
|
|
it = 0
|
|
z_start = self.get_zmin() + self.get_tot_wire_height() / 2 # (distance_coils/2 - windings * wire_height/2 + wire_height/2)*1e-3
|
|
R_start = self.get_R_inner() + self.get_tot_wire_width() / 2 # (R_inner + wire_width/2 )
|
|
log.debug(f"N = {self.get_N()}")
|
|
for xx in range(0, self.layers):
|
|
for zz in range(0, self.windings):
|
|
if self.winding_scheme == 2 and xx % 2 == 1 and zz == self.windings - 1:
|
|
continue # leave out every last winding in every second layer
|
|
z_pos = z_start + zz * self.get_tot_wire_height()
|
|
R_pos = R_start + xx * self.get_tot_wire_width()
|
|
|
|
# correct for different winding techniques and round wire
|
|
if self.winding_scheme == 1 or self.winding_scheme == 2:
|
|
if xx % 2 == 1:
|
|
z_pos += self.get_tot_wire_height() / 2
|
|
if self.is_round:
|
|
R_pos -= xx * ((2 - np.sqrt(3)) * self.get_tot_wire_width() / 2)
|
|
log.debug(f"lay = {xx}, wind = {zz}")
|
|
outer_raster[it] = [z_pos, R_pos]
|
|
# print(outer_raster[it])
|
|
it += 1
|
|
|
|
return outer_raster
|
|
|
|
def inner_raster(self, raster_value):
|
|
"""
|
|
Gives back inner raster for one wire around 0,0
|
|
Args:
|
|
raster_value (int): if N produces a N x N raster for rectangular and cut out of this for round
|
|
Returns: array containing raster around 0,0 [[z_pos_in_1,r_pos_in_1],...]
|
|
"""
|
|
# TODO: Less important, but rastering for round wires could be improved
|
|
if raster_value == 0:
|
|
raster_value = 1
|
|
log.info("raster value is 0 increased to 1")
|
|
|
|
if raster_value == 1:
|
|
return [0, 0]
|
|
|
|
if self.is_round & (raster_value % 2 == 0):
|
|
raster_value += 1
|
|
log.info(
|
|
f"for round wire rastering works better with uneven rastering value --> rastering value set to: {raster_value}")
|
|
|
|
inner_raster = np.zeros((raster_value ** 2, 2))
|
|
it = 0
|
|
for xx_in in range(0, raster_value):
|
|
for zz_in in range(0, raster_value):
|
|
z_pos_in = - self.wire_height / 2 + zz_in * self.wire_height / (raster_value - 1)
|
|
r_pos_in = - self.wire_width / 2 + xx_in * self.wire_width / (raster_value - 1)
|
|
|
|
inner_raster[it] = [z_pos_in, r_pos_in]
|
|
it += 1
|
|
|
|
# delete points out of round geometry
|
|
if self.is_round:
|
|
delete_list = []
|
|
for i in range(0, len(inner_raster)):
|
|
abs = np.sqrt(inner_raster[i, 0] ** 2 + inner_raster[i, 1] ** 2)
|
|
if abs > self.wire_width / 2:
|
|
delete_list.append(i)
|
|
inner_raster = np.delete(inner_raster, delete_list, axis=0)
|
|
|
|
return inner_raster
|
|
|
|
def full_raster(self, raster_value):
|
|
"""
|
|
|
|
Args:
|
|
raster_value: rastering value
|
|
|
|
Returns: [ wire 1:[[z_in1,r_in2], [z_in2,r_in,2],...],wire 2:[ [,] , ]...]...]
|
|
|
|
"""
|
|
outer_ras = self.winding_raster()
|
|
|
|
inner_ras = self.inner_raster(raster_value)
|
|
full_ras = np.zeros((len(outer_ras), len(inner_ras), 2))
|
|
|
|
for i in range(0, len(full_ras)):
|
|
full_ras[i] = outer_ras[i] + inner_ras
|
|
|
|
return full_ras
|
|
|
|
def plot_raster(self, raster_value=100):
|
|
full_structure = self.full_raster(raster_value) * 1e3
|
|
if self.get_coil_width() > self.get_coil_height():
|
|
extension = self.get_coil_width()
|
|
else:
|
|
extension = self.get_coil_height()
|
|
extension *= 1e3
|
|
|
|
plt.figure(77, figsize=(5, 5))
|
|
plt.scatter(full_structure[:, :, 1], full_structure[:, :, 0], linewidths=0.001)
|
|
plt.xlabel("radius [mm]")
|
|
plt.ylabel("z position [mm]")
|
|
plt.axvline(x=self.get_R_inner() * 1e3 - 0.1, lw=5, color="red")
|
|
plt.xlim(1e3 * self.get_R_inner() - 0.5, 1e3 * self.get_R_inner() + extension + 0.5)
|
|
plt.ylim(1e3 * self.get_zmin() - 0.5, 1e3 * self.get_zmin() + extension + 0.5)
|
|
|
|
plt.show()
|
|
|
|
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:.3f} mm, z_max = {self.get_zmax() * 1e3:.3f} mm")
|
|
print(
|
|
f"Radius = {self.radius * 1e3:.3f} mm, Radius_inner = {self.get_R_inner() * 1e3:.3f} mm, Radius_outer = {self.get_R_outer() * 1e3:.3f} mm")
|
|
print(
|
|
f"Coil width = {self.get_coil_width() * 1e3:.3f} mm, Coil height = {self.get_coil_height() * 1e3:.3f} mm"
|
|
)
|
|
print(
|
|
f"layers = {self.layers}, windings = {self.windings}, wire_width = {self.wire_width * 1e3:.3f}, wire_height = {self.wire_height * 1e3:.3f} mm, round wire = {self.is_round} ")
|
|
print(" ")
|
|
|
|
def print_basic_info(self):
|
|
if self.HH == 1:
|
|
print("Offset coil:")
|
|
else:
|
|
print("Gradient coil:")
|
|
print(
|
|
f" layers = {self.layers}, windings = {self.windings} \n "
|
|
f" wire: core = {self.wire_height * 1e3:.2f} mm, tot = {self.get_tot_wire_height() * 1e3:.2f} mm, total length for pair = {2 * self.get_wire_length():.0f} m"
|
|
)
|
|
print(
|
|
f" Distance = {self.distance * 1e3:.2f} mm, Radius = {self.radius * 1e3:.2f} mm"
|
|
)
|
|
print(
|
|
f" Coil width = {self.get_coil_width() * 1e3:.2f} mm, Coil height = {self.get_coil_height() * 1e3:.2f} mm"
|
|
)
|
|
print("")
|
|
|
|
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 * (
|
|
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))
|
|
B *= 1e4 # conversion Gauss
|
|
return B
|
|
|
|
@staticmethod
|
|
def make_axis(lim,precision,two_axis=True):
|
|
"""
|
|
Gives back two arrays x and z in mm, symmetric around zero
|
|
Args:
|
|
lim: positive (and negative if lim2 = None) limit of axis
|
|
precision: nr_points per mm
|
|
two_axis: If True --> gives back two arrays with given parameter
|
|
|
|
Returns: array
|
|
|
|
"""
|
|
nr_points = (2 * lim) * precision + 1
|
|
nr_points = int(nr_points)
|
|
z = np.linspace(-lim, lim, nr_points)
|
|
|
|
if two_axis:
|
|
x = np.linspace(-lim, lim, nr_points)
|
|
return x, z
|
|
return z
|
|
|
|
@staticmethod
|
|
def mm_it(x, x_pos):
|
|
"""
|
|
Takes argument in mm returns position in x-array
|
|
Args:
|
|
x: array of interest
|
|
x_pos: position to find position of
|
|
|
|
Returns: int of position of x_pos in array x
|
|
|
|
"""
|
|
precision = len(x)//(2*np.abs(x[0]))
|
|
it = int(len(x)//2 + precision * x_pos)
|
|
return it
|
|
|
|
@staticmethod
|
|
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)
|
|
|
|
@staticmethod
|
|
def __B_z_loop(I_current, r_loop, z_loop, r_pos, z_pos):
|
|
"""
|
|
calculate z-component of B-field at position r and z for each individual loop
|
|
Args:
|
|
I_current: Current in [A]
|
|
r_loop: Radial position of current loop
|
|
z_loop: Axial position of current loop along symmetry axis
|
|
r_pos: radial position of calculation in [m]
|
|
z_pos: axial position of calculation in [m]
|
|
|
|
Returns: z component of B-field at position r, z in G
|
|
|
|
"""
|
|
B_z = 2e-7 * I_current * 1 / ((r_loop + r_pos) ** 2 + (z_pos - z_loop) ** 2) ** (1 / 2) * (
|
|
sp.ellipk(BCoil.k_sq(r_loop, z_loop, r_pos, z_pos)) + sp.ellipe(BCoil.k_sq(r_loop, z_loop, r_pos, z_pos)) * (
|
|
r_loop ** 2 - r_pos ** 2 - (z_pos - z_loop) ** 2) / ((r_loop - r_pos) ** 2 + (z_pos - z_loop) ** 2))
|
|
B_z *= 1e4 # conversion to gauss
|
|
return B_z
|
|
|
|
@staticmethod
|
|
def __B_r_loop(I_current, r_loop, z_loop, r_pos, z_pos):
|
|
"""
|
|
calculate z-component of B-field at position r and z for each individual loop
|
|
Args:
|
|
I_current: Current in [A]
|
|
r_loop: Radial position of current loop
|
|
z_loop: Axial position of current loop along symmetry axis
|
|
r_pos: radial position of calculation in [m]
|
|
z_pos: axial position of calculation in [m]
|
|
|
|
Returns: z component of B-field at position r, z in G
|
|
|
|
"""
|
|
B_r = 2e-7 * I_current / r_pos * (z_pos - z_loop) / ((r_loop + r_pos) ** 2 + (z_pos - z_loop) ** 2) ** (1 / 2) * (
|
|
-sp.ellipk(BCoil.k_sq(r_loop, z_loop, r_pos, z_pos)) + sp.ellipe(BCoil.k_sq(r_loop, z_loop, r_pos, z_pos)) * (
|
|
r_loop ** 2 + r_pos ** 2 + (z_pos - z_loop) ** 2) / ((r_loop - r_pos) ** 2 + (z_pos - z_loop) ** 2))
|
|
B_r *= 1e4 # conversion to gauss
|
|
return B_r
|
|
|
|
def B_field_z(self,I_current):
|
|
pass
|
|
|
|
def B_field(self, I_current, x, z, raster = 7):
|
|
"""
|
|
Returns Bz along z-axis and B_x along x-axis,
|
|
HH = +1 --> Helmholtz configuration, HH = -1 --> Anti Helmholtz configuration
|
|
"""
|
|
|
|
x_SI = x * 1e-3
|
|
z_SI = z * 1e-3
|
|
|
|
if x[0] <= 0: # divide array into positive and negative values
|
|
x_negative = np.abs([el for el in x_SI if el < 0])
|
|
x_positive = 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_positive))
|
|
B_x_neg = np.zeros(len(x_negative))
|
|
B_x_zero = x_zero # 0 in x-array --> field is zero at this position
|
|
|
|
calc_raster = self.full_raster(raster)
|
|
|
|
if self.get_N() != len(calc_raster):
|
|
log.error("N is not equal length of raster")
|
|
|
|
rastering_value = len(calc_raster[0])
|
|
|
|
I_current /= rastering_value # divide current into smaller currents for mapping the whole wire
|
|
# start = time.time()
|
|
for wire in range(0, self.get_N()):
|
|
for ii in range(0, rastering_value):
|
|
# extract position information out of raster
|
|
z_pos = calc_raster[wire, ii, 0]
|
|
r_pos = calc_raster[wire, ii, 1]
|
|
|
|
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_positive, 0) + BCoil.__B_r_loop(
|
|
self.HH * I_current, r_pos, -z_pos, x_positive, 0)
|
|
B_x_neg += BCoil.__B_r_loop(I_current, r_pos, z_pos, x_negative, 0) + BCoil.__B_r_loop(
|
|
self.HH * I_current, r_pos, -z_pos, x_negative, 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))
|
|
# end = time.time()
|
|
# print(f"time = {end - start} s")
|
|
return B_z, B_x
|
|
|
|
def max_field(self, I_current, raster = 7):
|
|
B_z_0 = 0
|
|
calc_raster = self.full_raster(raster)
|
|
|
|
if self.get_N() != len(calc_raster):
|
|
log.error("N is not equal length of raster")
|
|
|
|
rastering_value = len(calc_raster[0])
|
|
|
|
I_current /= rastering_value # divide current into smaller currents for mapping the whole wire
|
|
|
|
for wire in range(0, self.get_N()):
|
|
for ii in range(0, rastering_value):
|
|
# extract position information out of raster
|
|
z_pos = calc_raster[wire, ii, 0]
|
|
r_pos = calc_raster[wire, ii, 1]
|
|
|
|
B_z_0 += BCoil.__B_z_loop(I_current, r_pos, z_pos, 0, 0) + BCoil.__B_z_loop(self.HH * I_current, r_pos,
|
|
-z_pos, 0, 0)
|
|
print(f" Max Field = {B_z_0:.2f} G")
|
|
return B_z_0
|
|
|
|
def max_gradient(self, i_current, raster = 7):
|
|
x, z = BCoil.make_axis(0.5, 1000)
|
|
Bz, Bx = self.B_field(i_current, x, z, raster = raster)
|
|
|
|
bz_grad = BCoil.grad(Bz, z)
|
|
bx_grad = BCoil.grad(Bx, x)
|
|
print(f" Max z Gradient = {bz_grad[len(z)//2]:.2f} G/cm")
|
|
print(f" Max x Gradient = {bx_grad[len(x)//2]:.2f} G/cm")
|
|
# return bx_grad[len(x)//2]
|
|
|
|
|
|
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
|
|
"""
|
|
|
|
# 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
|
|
calc_raster = self.full_raster(raster)
|
|
|
|
if self.get_N() != len(calc_raster):
|
|
log.error("N is not equal length of raster")
|
|
|
|
rastering_value = len(calc_raster[0])
|
|
|
|
I_current /= rastering_value # divide current into smaller currents for mapping the whole wire
|
|
# start = time.time()
|
|
for wire in range(0, self.get_N()):
|
|
for ii in range(0, rastering_value):
|
|
# extract position information out of raster
|
|
z_pos = calc_raster[wire, ii, 0]
|
|
r_pos = calc_raster[wire, ii, 1]
|
|
|
|
# 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)
|
|
# Bz 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)
|
|
B_tot_z = np.abs(B_tot_z)
|
|
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])
|
|
|
|
calc_raster = self.full_raster(raster)
|
|
if self.get_N() != len(calc_raster):
|
|
log.error("N is not equal length of raster")
|
|
rastering_value = len(calc_raster[0])
|
|
|
|
# TODO: why division by zero?
|
|
|
|
I_current /= rastering_value # divide current into smaller currents for mapping the whole wire
|
|
# start = time.time()
|
|
for el in range(0, len(x_SI)):
|
|
for wire in range(0, self.get_N()):
|
|
for ii in range(0, rastering_value):
|
|
# extract position information out of raster
|
|
z_pos = calc_raster[wire, ii, 0]
|
|
r_pos = calc_raster[wire, ii, 1]
|
|
|
|
# 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
|
|
|
|
@staticmethod
|
|
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)
|
|
print(x)
|
|
z = np.arange(-z_lim, z_lim, 5)
|
|
print(z)
|
|
x_m, z_m = np.meshgrid(x, z)
|
|
if self.is_round:
|
|
B = self.B_multiple_3d(I_current, x, z, raster=3)
|
|
else:
|
|
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[zz, xx] > 15:
|
|
B[zz, xx, :] /= B_tot[zz, xx] / 15
|
|
|
|
plt.figure()
|
|
plt.quiver(x_m, z_m, B[:, :, 0], B[:, :, 1])
|
|
plt.xlabel("x-axis [mm]")
|
|
plt.ylabel("z-axis [mm]")
|
|
plt.show()
|
|
|
|
@staticmethod
|
|
def curv(B_f, z):
|
|
return np.gradient(np.gradient(B_f, z), z) * 1e2
|
|
|
|
@staticmethod
|
|
def grad(B_f, z):
|
|
return np.gradient(B_f, z) * 1e1
|
|
|
|
def B_quick_plot(self, I_current, abs=True, x_lim=50, z_lim=50, nr_points=200):
|
|
x = np.linspace(-x_lim, x_lim, nr_points)
|
|
z = np.linspace(-z_lim, z_lim, nr_points)
|
|
if abs:
|
|
B_z, B_x = self.B_tot_along_axis(I_current, x, z)
|
|
else:
|
|
B_z, B_x = self.B_field(I_current, x, z)
|
|
|
|
plt.figure(11)
|
|
plt.plot(z, B_z, linestyle="solid", label=r"$B_{tot}$ along z-axis")
|
|
plt.plot(x, B_x, label=r"$B_{tot}$ along x-axis")
|
|
plt.title("B-field")
|
|
plt.ylabel(r"B-field [G]")
|
|
plt.xlabel("x-axis / z-axis [mm]")
|
|
plt.legend()
|
|
plt.show()
|
|
|
|
def B_grad_quick_plot(self, I_current, x_lim=50, z_lim=50, nr_points=200):
|
|
x = np.linspace(-x_lim, x_lim, nr_points)
|
|
z = np.linspace(-z_lim, z_lim, nr_points)
|
|
|
|
B_z, B_x = self.B_field(I_current, x, z)
|
|
# TODO: Think about Gradient calculation of tot B_field
|
|
B_z = BCoil.grad(B_z, z)
|
|
B_x = BCoil.grad(B_x, x)
|
|
|
|
plt.figure(12)
|
|
plt.plot(z, B_z, linestyle="solid", label=r"z grad of Bz along z-axis")
|
|
plt.plot(x, B_x, label=r"x Grad of B_x along x-axis")
|
|
plt.title("Gradient of B-field")
|
|
plt.ylabel(r"B-field [G/cm]")
|
|
plt.xlabel("x-axis / z-axis [mm]")
|
|
plt.legend()
|
|
plt.show()
|
|
|
|
def B_curv_quick_plot(self, I_current, x_lim=50, z_lim=50, nr_points=200):
|
|
x = np.linspace(-x_lim, x_lim, nr_points)
|
|
z = np.linspace(-z_lim, z_lim, nr_points)
|
|
|
|
B_z, B_x = self.B_field(I_current, x, z)
|
|
B_z = BCoil.curv(B_z, z)
|
|
B_x = BCoil.curv(B_x, x)
|
|
|
|
plt.figure(13)
|
|
plt.plot(z, B_z, linestyle="solid", label=r"$B_{tot}$ along z-axis")
|
|
plt.plot(x, B_x, label=r"$B_{tot}$ along x-axis")
|
|
plt.title("Curvature of B-field")
|
|
plt.ylabel(r"Curv. B-field [G/cm$^2$]")
|
|
plt.xlabel("x-axis / z-axis [mm]")
|
|
plt.legend()
|
|
plt.show()
|
|
|
|
def Bz_plot_HH_comp(self, Coil2, I_current, x, z):
|
|
B_z, B_x = self.B_field(I_current, x, z)
|
|
B_z_2, B_x_2 = Coil2.B_field(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 Bz along z-axis")
|
|
|
|
# Field plot
|
|
##########################
|
|
plt.subplot(2, 1, 1)
|
|
plt.plot(z, B_z, linestyle="solid", label=r"$Bz$")
|
|
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"$Bz$ [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 Bz$")
|
|
plt.plot(z, B_z_curvature_2, linestyle="solid", label=r"$\nabla_z^2 B_{z2}$")
|
|
|
|
plt.ylabel(r"$\nabla_z^2 Bz [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 for one coil
|
|
|
|
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.power(I_current, T)
|
|
|
|
Voltage = self.resistance(T) * I_current
|
|
|
|
#print("")
|
|
print("Cooling:")
|
|
print(f" Current I = {I_current} A")
|
|
print(f" current density = {j_dens:.2f} A/mm^2")
|
|
print(f" Power = {Power:.2f} W")
|
|
print(f" U = {Voltage:.2f} V")
|
|
|
|
def power(self, I_current, T):
|
|
"""
|
|
Power for one coil
|
|
"""
|
|
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 one rectangular coil via empirical formular by perry
|
|
|
|
"""
|
|
L = 4 * np.pi * (self.windings * self.layers) ** 2 * (1e2 * self.radius) ** 2 / (
|
|
0.2317 * 100 * self.radius + 0.44 * 100 * self.get_coil_height() + 0.39 * 100 * self.get_coil_width())
|
|
return L * 1e-9
|
|
|
|
@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
|
|
|
|
def resistance(self, tempr):
|
|
"""
|
|
Calculates resistance of one coil of configuration
|
|
|
|
Parameters
|
|
----------
|
|
tempr : double
|
|
temperature in degree Celsius
|
|
|
|
Returns
|
|
-------
|
|
double
|
|
resistance in Ohm
|
|
"""
|
|
|
|
return BCoil.resistivity_copper(tempr) * self.get_wire_length() / self.get_wire_area()
|
|
|
|
def tau(self, tempr = 22):
|
|
return self.induct_perry()/self.resistance(tempr)
|
|
|
|
|
|
def main():
|
|
Coil_1 = BCoil(HH=1, distance=54, radius=48, layers=8, windings=8, wire_height=0.5, wire_width=0.5,
|
|
insulation_thickness=0.05, is_round=True, winding_scheme=2)
|
|
# Coil_1.plot_raster()
|
|
res_m = BCoil.resistivity_copper(20) * 1 / Coil_1.get_wire_area()
|
|
print(res_m)
|
|
|
|
|
|
# main()
|
|
if __name__ == "__main__":
|
|
log.info("Run main Coil class")
|
|
main()
|