hoenen/FieldCalculations
1
0
Fork 0
Code Issues Pull Requests repo.project_board Releases Wiki Activity

Initial upload of code to calculate the field of rectangular coils. See Lenny's BA Thesis for details.

Browse Source
This commit is contained in:
Hoenen 2022-10-14 12:44:13 +02:00
commit 1167ee64b8
1 changed files with 625 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

625
main.py Normal file
View File

@ -0,0 +1,625 @@
import numpy as np
import matplotlib.pyplot as plt
import logging as log
import scipy.constants as cs
from mpl_toolkits.axes_grid1 import host_subplot
import mpl_toolkits.axisartist as AA
# TODO:
# -Different wire winding schemes
# -Plot a cut through the Coil
# -Plotte den erlaubten unsicherheitsbereich in die B-Feld plots. (durchsichtige Box sollte reichen)
# -Möchte ich HHCoils als klasse oder funktion?
# --> Einmal B berechnen und dann funktionen druaf anwenden oder einmal objekt kreieren?
# --> Mit sinnvollen funktionen eigentlich kein unterschied
# -Coords in meter oder mm
"""
The following functions allow the calculation of the magnetic field of rectangular coils.
The most important important inputs the following:
-"distance" describes the distance between two coils in Helmholtz (HH) config in mm. In the case of a single coil,
it is the offset of its plane from zero in mm.
-"a" is the side length of the coil(s) in +/- x direction in mm, when the coil is lying in the z-plane.
-"b" is the side length of the coil(s) in +/- y direction in mm, when the coil is lying in the z-plane.
-"plane" describes the plane in which the coil lies. "plane = 'z'" describes a coil which extends in x and y
direction with its sides and has an offset of "distance" mm in z direction.
-The wire parameters layers, windings, wire_width and wire_height describe the winding of the coil. The side lengths
"a" and "b" are defined with regard to the center of the coil winding. That means, a coil with side length a has
wires extending from a-h to a+h, where h is the height of the accumulated winding layers.
-NO REAL WINDING SCHEME IS IMPLEMENTED.
Instead of crossing wires, the coil is made up of many single wire coils lying next to each other. In reality, a
layer of 5 wires is followed by a layer of 4 wires, as the next layer always sits in the groves of the layer below.
THIS IS NOT IMPLEMENTED. Here "layers =3, windings = 5" means that 15 wires in total sit perfectly on top of each
other. They are not filling any gaps and have a height of exactly 3*(wire_height+insulation thickness).
-The calculation is done once per wire, as if it where a coil without physical extend in the center of the wire.
-The parameter "center_offset" allows for coils that are not centered around an axis, but have some offset.
The center offset is to be given as a list of the three coords, that define the center of the coil (pair) in m(!).
Note that the "distance" is calculated from this new center. For non ideal/crooked HH coils, define 2 individual
coils with two individual centers.
-The calculation of the field is done a whichever coordinates are given as an input to the various field calculation
functions (mainly calc_B_field_coil() or HH_coil_calc_B_field). The coords have to be given in units of meters and
the output will be a set of three components of the field for each coordinate in Gauss. The field can also be given
in Tesla.
"""
class RectangularCoil:
def __init__(self, distance, a, b, plane, layers, windings, wire_width, wire_height, I, mu0=cs.mu_0,
center_offset=None, insulation_thickness=0., is_round=True, winding_offset=False, unit="Gauss"):
# if center is None:
# center = [0, 0, 0]
if not is_round:
log.warning("Non-round wires not implemented yet. Please choose round wires instead.")
if is_round:
if wire_width != wire_height:
log.error('Wire is round but width != height')
self.distance = distance * 1e-3
self.a = a * 1e-3
self.b = b * 1e-3
self.plane = plane
self.layers = layers
self.windings = windings
self.wire_width = wire_width * 1e-3
self.wire_height = wire_height * 1e-3
self.I = I
self.mu0 = mu0
self.insulation_thickness = insulation_thickness * 1e-3
self.is_round = is_round
self.winding_offset = winding_offset
self.unit = unit
if center_offset is None:
self.center_offset = [0,0,0]
else:
self.center_offset = center_offset
def get_N(self):
"""
Calculates number of windings
"""
return self.layers * self.windings
def wires(self):
"""
According to the wire geometry the winding is done and layers of wires are stacked. Each loop of wire can then be treated as a coil each.
"""
width = self.wire_width + 2 * self.insulation_thickness
coil_width = self.windings * width
coil_height = self.layers * width
Coilwires = np.zeros( (self.get_N(), 3) ) # Each wire of the coil will be stored in Coilwires with Params = (distance, a, b)
"""
This winding scheme just lays wires next to each other and on top of each other, but only their radius etc is used. The geometry is not yet accounted for.
"""
# Fill in different Wirewinding schemes later
"""
The Coil is centered around the parameters d,a,b
"""
for L in range(self.layers):
layer_offset = -coil_height / 2 + (2 * L + 1) * width / 2
a_layer = self.a - layer_offset
b_layer = self.b - layer_offset
for W in range(self.windings):
d_offset = -coil_width / 2 + (2 * W + 1) * width / 2
distance_wire = self.distance - d_offset
Coilwires[L * self.windings + W] = np.array([distance_wire, a_layer, b_layer])
return Coilwires
def calc_B_field_z_wire(self, coords, wire):
x = coords[0]
y = coords[1]
z = coords[2]
a_wire = wire[1] / 2 # The input is the side-length, but the calculation uses half-lengths of each side
b_wire = wire[2] / 2
c1 = a_wire + x
c2 = a_wire - x
c3 = -c2
c4 = -c1
d1 = y + b_wire
d2 = d1
d3 = y - b_wire
d4 = d3
r1 = np.sqrt((a_wire + x) ** 2 + (y + b_wire) ** 2 + z ** 2)
r2 = np.sqrt((a_wire - x) ** 2 + (y + b_wire) ** 2 + z ** 2)
r3 = np.sqrt((a_wire - x) ** 2 + (y - b_wire) ** 2 + z ** 2)
r4 = np.sqrt((a_wire + x) ** 2 + (y - b_wire) ** 2 + z ** 2)
const = self.mu0 * self.I / (4 * cs.pi)
Bx = const * (z / (r1 * (r1 + d1)) - z / (r2 * (r2 + d2)) + z / (r3 * (r3 + d3)) - z / (r4 * (r4 + d4)))
By = const * (z / (r1 * (r1 + c1)) - z / (r2 * (r2 - c2)) + z / (r3 * (r3 + c3)) - z / (r4 * (r4 - c4)))
Bz = const * (-1 * d1 / (r1 * (r1 + c1)) - c1 / (r1 * (r1 + d1)) + d2 / (r2 * (r2 - c2)) - c2 / (
r2 * (r2 + d2)) - d3 / (r3 * (r3 + c3)) - c3 / (r3 * (r3 + d3)) + d4 / (r4 * (r4 - c4)) - c4 / (
r4 * (r4 + d4)))
return Bx, By, Bz
def calc_B_field_any_wire(self, coords, wire):
x_coords = coords[0]
y_coords = coords[1]
z_coords = coords[2]
distance_wire = wire[0]
x_center = self.center_offset[0]
y_center = self.center_offset[1]
z_center = self.center_offset[2]
"""
The input coords are the actual physical coordinates, at which the BField is to be calculated.
The calculation is only done for Coils that lie in the z-plane. Therefore the input has to be rotated,
and shifted in z-direction (depending on the winding) and the resulting fields are then rotated back.
Therefore the parameter "plane" has to be given to determine in which plane the coil lies.
"""
if self.plane == "x":
x_center, y_center, z_center = -z_center, y_center, x_center
x_coords, y_coords, z_coords = -z_coords, y_coords, x_coords
elif self.plane == "y":
x_center, y_center, z_center = x_center, -z_center, y_center
x_coords, y_coords, z_coords = x_coords, -z_coords, y_coords
elif self.plane == "z":
pass
else:
log.error("Unexpected plane. Please choose 'x', 'y' or 'z'.")
"""
The rotated coil now lies in the z-plane.
By subtracting the distance_wire from the coords, the distnace of the wire to the z-plane is accounted for.
Then, the accordingly rotated coordinates of the center of the coils can be subtracted likewise from the coords, where B is calculated.
Therefore not the coil in the coordsys is shifted by +offset, but the coordsystem is shifted by -offset, which results in the same offset.
"""
offset_coords = np.array([x_coords - x_center, y_coords - y_center, (z_coords - distance_wire - z_center)])
x_field, y_field, z_field = self.calc_B_field_z_wire(offset_coords, wire)
if self.plane == "x":
x_field, y_field, z_field = z_field, y_field, -x_field
elif self.plane == "y":
x_field, y_field, z_field = x_field, z_field, -y_field
elif self.plane == "z":
pass
else:
log.error("Unexpected plane. Please choose 'x', 'y' or 'z'.")
field = np.array([x_field, y_field, z_field])
return field
def calc_B_field_coil(self, coords): # B Field of one Coil. No HH config
B_coil = np.zeros(np.shape(coords))
Coil_Wires = self.wires()
"""
To calculate B-Field of whole coil, calculate B-field once per wire
"""
for l in range(self.layers):
for w in range(self.windings):
B_coil += self.calc_B_field_any_wire(coords, Coil_Wires[l * self.windings + w])
if self.unit.lower() == "gauss":
B_coil = B_coil * 1e4
elif self.unit.lower() == "tesla":
pass
else:
log.warning(f"{self.unit} Not recognized. Please choose unit 'gauss' or 'tesla'.")
return B_coil
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.power(I_current, T)
Voltage = self.resistance(T) * I_current
print("")
print("Cooling:")
print(f" current density = {j_dens} A/mm^2")
print(f" Power = {Power} W")
print(f" U = {Voltage} V")
def power(self, I_current, T):
P = self.resistance(T) * I_current ** 2
return P
def resistance(self, T):
"""
Calculates resistance of one coil of configuration
Parameters
----------
T : double
temperature in degree Celsius
Returns
-------
double
resistance in Ohm
"""
return RectangularCoil.resistivity_copper(T) * self.get_wire_length() / self.get_wire_area()
def get_wire_length(self):
"""
:return: Approximate length of wire per coil (m)
"""
return self.get_N() * 2 * (self.a + self.b)
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
@staticmethod
def resistivity_copper(T):
R_ref = 1.68e-8 # Resistance at T=20
T_ref = 20 # degree celsius
alpha = 0.00393
R = R_ref * (1 + alpha * (T - T_ref))
return R
@staticmethod
def HHCoil_calc_B_field(Coords, d, a, b, plane, layers, windings, wire_width, wire_height, I, mu0=4e-7 * np.pi,
center_offset=None,
insulation_thickness=0., is_round=True, winding_offset=False, unit="Gauss"):
Coil1 = RectangularCoil(distance=d / 2, a=a, b=b, plane=plane, layers=layers, windings=windings,
wire_width=wire_width, wire_height=wire_height, I=I, mu0=mu0, center_offset=center_offset,
insulation_thickness=insulation_thickness, is_round=is_round,
winding_offset=winding_offset, unit=unit)
Coil2 = RectangularCoil(distance=-d / 2, a=a, b=b, plane=plane, layers=layers, windings=windings,
wire_width=wire_width, wire_height=wire_height, I=I, mu0=mu0, center_offset=center_offset,
insulation_thickness=insulation_thickness, is_round=is_round,
winding_offset=winding_offset, unit=unit)
B = Coil1.calc_B_field_coil(Coords)
B += Coil2.calc_B_field_coil(Coords)
return B
@staticmethod
def HHCoil_calc_total_B_field(Coords, d, a, b, plane, layers, windings, wire_width, wire_height, I,
mu0=4e-7 * np.pi,
center_offset=None,
insulation_thickness=0., is_round=True, winding_offset=False, unit="Gauss"):
Coil1 = RectangularCoil(distance=d / 2, a=a, b=b, plane=plane, layers=layers, windings=windings,
wire_width=wire_width, wire_height=wire_height, I=I, mu0=mu0, center_offset=center_offset,
insulation_thickness=insulation_thickness, is_round=is_round,
winding_offset=winding_offset, unit=unit)
Coil2 = RectangularCoil(distance=-d / 2, a=a, b=b, plane=plane, layers=layers, windings=windings,
wire_width=wire_width, wire_height=wire_height, I=I, mu0=mu0, center_offset=center_offset,
insulation_thickness=insulation_thickness, is_round=is_round,
winding_offset=winding_offset, unit=unit)
B = Coil1.calc_B_field_coil(Coords)
B += Coil2.calc_B_field_coil(Coords)
B = np.sqrt(B[0] ** 2 + B[1] ** 2 + B[2] ** 2)
return B
@staticmethod
def grad(B_f, z):
return np.gradient(B_f, z) # *1e3)
@staticmethod
def curv(B_f, z):
return np.gradient(np.gradient(B_f, z)) # *1e3), z*1e3)
@staticmethod
def calc_B_i_grad(Field_i, Coords_i):
"""
The '_i' is supposed to indicate, that this function only takes scalar inputs. Total field is possible as well.
"""
Grad_Field = RectangularCoil.grad(Field_i, Coords_i)
return Grad_Field
@staticmethod
def calc_B_i_curv(Field_i, Coords_i):
Curv_Field = RectangularCoil.curv(Field_i, Coords_i)
return Curv_Field
@staticmethod
def plot_B_i(Field_i, Coords_i, show=False):
plt.plot(Coords_i, Field_i, ls="-", label=r"$B_{tot}$")
plt.ylabel("Total Magnetic Field along axis")
plt.xlabel("axis [m]")
plt.ylabel("Magnetic field strength [G]")
plt.grid()
if show:
plt.show()
@staticmethod
def plot_B_grad_i(Field_i, Coords_i, show=False):
Grad_Field_i = RectangularCoil.calc_B_i_grad(Field_i, Coords_i)
plt.plot(Coords_i, Grad_Field_i, ls="-.", label=r"$grad(B_{tot})$")
plt.ylabel("Total Magnetic Field Gradient along axis")
plt.xlabel("axis [m]")
plt.ylabel("Magnetic field Gradient [G/m]")
plt.grid()
if show:
plt.show()
@staticmethod
def plot_B_curv_i(Field_i, Coords_i, show=False):
Curv_Field_i = RectangularCoil.calc_B_i_curv(Field_i, Coords_i)
plt.plot(Coords_i, Curv_Field_i, ls="-.", label=r"$curv(B_{tot})$")
plt.ylabel("Total Magnetic Field Curvature along axis")
plt.xlabel("axis [m]")
plt.ylabel("Magnetic field Curvature [G/m^2]")
plt.grid()
if show:
plt.show()
@staticmethod
def plot_B_everything_fancy(Field_i, Coords_i, f=0, show=True):
Grad_Field_i = RectangularCoil.grad(Field_i, Coords_i)
Curv_Field_i = RectangularCoil.curv(Field_i, Coords_i)
xmin, xmax = np.min(Coords_i) * (1 - f), np.max(Coords_i) * (1 + f)
y0min, y0max = np.min(Field_i) * (1 - f), np.max(Field_i) * (1 + f)
y1min, y1max = np.min(Grad_Field_i) * (1 - f), np.max(Grad_Field_i) * (1 + f)
y2min, y2max = np.min(Curv_Field_i[2:-2:]) * (1 - f), np.max(Curv_Field_i[2:-2:]) * (1 + f)
host = host_subplot(111, axes_class=AA.Axes)
plt.subplots_adjust(right=0.75)
par1 = host.twinx()
par2 = host.twinx()
new_fixed_axis = par2.get_grid_helper().new_fixed_axis
par1.axis["right"] = new_fixed_axis(loc="right", axes=par1, offset=(0, 0))
par1.axis["right"].toggle(all=True)
par2.axis["right"] = new_fixed_axis(loc="right", axes=par2, offset=(60, 0))
par2.axis["right"].toggle(all=True)
host.set_xlim(xmin, xmax)
host.set_ylim(y0min, y0max)
host.set_xlabel("axis [m]")
host.set_ylabel("Field strength [G]")
par1.set_ylabel("Field Gradient [G/m]")
par2.set_ylabel("Field Curvature [G/m^2]")
p1, = host.plot(Coords_i, Field_i, label="Srength")
p2, = par1.plot(Coords_i, Grad_Field_i, label="Gradient")
p3, = par2.plot(Coords_i, Curv_Field_i, label="Curvature")
par1.set_ylim(y1min, y1max)
par2.set_ylim(y2min, y2max)
host.legend()
host.axis["left"].label.set_color(p1.get_color())
par1.axis["right"].label.set_color(p2.get_color())
par2.axis["right"].label.set_color(p3.get_color())
host.grid()
par1.grid(alpha=0.3)
par2.grid(alpha=0.3)
plt.draw()
if show:
plt.show()
@staticmethod
def plot_B_everything_simple(Field_i, Coords_i, f=1.1):
RectangularCoil.plot_B_i(Field_i, Coords_i)
RectangularCoil.plot_B_grad_i(Field_i, Coords_i)
RectangularCoil.plot_B_curv_i(Field_i, Coords_i)
plt.legend()
plt.show()
@staticmethod
def homogenity_1D(Field_i, Coords_i, Range, perc=False):
"""
first check whether 0 is in the list of coords or not
if not, calc B(0).
"""
clipped_coords = Coords_i[(-Range <= Coords_i)]
clipped_coords = clipped_coords[(clipped_coords <= Range)]
near0 = np.min(np.abs(Coords_i))
index = np.where(np.abs(Coords_i) == near0)
if near0 != 0:
print(
f"Point (0,0,0) not part of given coordinate list. Closest point is {near0} and will be used for calculation.")
else:
pass
index1 = int(np.where(clipped_coords[0] == Coords_i)[0])
index2 = int(np.where(clipped_coords[-1] == Coords_i)[0])
homogenity = np.abs((Field_i - Field_i[index]) / Field_i[index])
if perc:
homogenity = homogenity[index1:index2 + 1:] * 100
print(f"The max. inhomogenity in the range of +/- {Range}m was {np.max(np.abs(homogenity)) * 100} %")
else:
homogenity = homogenity[index1:index2 + 1:]
print(f"The max. inhomogenity in the range of +/- {Range}m was {np.max(np.abs(homogenity)) }")
return clipped_coords, homogenity
def Mutual_Inductance_Rover(self):
"""
Mutual inductance of two equal parallel rectangles according to Rosa & Grover
p1,p2,p3 only introduced for readability
"""
a = self.a * 100 # Converion to [cm]
b = self.b * 100
d = self.distance * 100
p1 = (a + np.sqrt(a ** 2 + d ** 2)) * np.sqrt(b ** 2 + d ** 2) / ((a + np.sqrt(a ** 2 + b ** 2 + d ** 2)) * d)
p2 = (b + np.sqrt(b ** 2 + d ** 2)) * np.sqrt(a ** 2 + d ** 2) / ((b + np.sqrt(a ** 2 + b ** 2 + d ** 2)) * d)
p3 = np.sqrt(a ** 2 + b ** 2 + d ** 2) - np.sqrt(a ** 2 + d ** 2) - np.sqrt(b ** 2 + d ** 2) + d
M = 4 * (a * np.log(p1) + b * np.log(p2)) + 8 * p3
M *= self.get_N()**2 * 1e-9 #Conversion to SI and accounting for N wires
return M
def Self_Inductance_Rover(self):
a = self.a * 100 #Converion to [cm]
b = self.b * 100
r = self.wire_width/2 *100
diag = np.sqrt(a**2 + b**2)
L = 4 * ( (a+b) * np.log(2*a*b / r) - a * np.log(a+diag) - b * np.log(b+diag) - 7/4 * (a+b) + 2 * (diag + r) )
L *= self.get_N()**2 * 1e-9 #Conversion to SI and accounting for N wires
return L
@staticmethod
def Self_Inductance_Abbot(alpha, beta, gamma, N):
"""
Formula for self inductance of SQUARE coil by J.J.Abbott
"""
L = (1e-5 * alpha**2 * N**2) / (3*alpha + 9*beta + 10*gamma)
return L
def Total_Inductance_Abbot(self):
"""
Above Formula for self inductance used, but with modified alpha to roughly represent rect. coils.
Total inductance then according to Grover
"""
alpha = (self.a + self.b) / 2
width = self.wire_width + 2 * self.insulation_thickness
hight = self.wire_height + 2 * self.insulation_thickness
beta = self.windings * width
gamma = self.layers * hight
N = self.get_N()
dist = self.distance
L_tot = RectangularCoil.Self_Inductance_Abbot(alpha, dist + gamma, gamma, N) + RectangularCoil.Self_Inductance_Abbot(alpha, dist - gamma, gamma, N) - 2 * RectangularCoil.Self_Inductance_Abbot(alpha, dist, gamma, N) + 2 * RectangularCoil.Self_Inductance_Abbot(alpha, beta, gamma, N)
return L_tot
def Total_Inductance_Rover(self):
L = self.Self_Inductance_Rover()
M = self.Mutual_Inductance_Rover()
L_tot = 2 * (L+M)
return L_tot
@staticmethod
def I_time_response(U, R, L, t):
I = np.abs(U) / R * (1 - np.exp(-R * t*1e-3 / L))
return I
def main():
"""
Up to 1e5 steps for up to 100 wires take a reasonable amount of time
"""
x0 = np.linspace(-0.001, 0.001, 1000)
y0 = np.ones(np.shape(x0)) * 0 / 2
z0 = np.ones(np.shape(y0)) * 0 / 2
Coordsx = np.array([x0, y0, z0])
#Coordsx2 = np.array([x0-0.01, y0, z0])
y1 = np.linspace(-0.001, 0.001, 1000)
x1 = np.ones(np.shape(y1)) * 0 / 2
z1 = np.ones(np.shape(y1)) * 0 / 2
Coordsy = np.array([x1, y1, z1])
z2 = np.linspace(-0.001, 0.001, 1000)
y2 = np.ones(np.shape(z2)) * 0 / 2
x2 = np.ones(np.shape(y2)) * 0 / 2
Coordsz = np.array([x2, y2, z2])
Rec_off = RectangularCoil( 78.4 / 2, 129.4, 206.4, "z", 10, 10, 0.5, 0.5, 1, center_offset=[0, 0.000, 0])
Rec_off_2 = RectangularCoil(-78.4 / 2, 129.4, 206.4, "z", 10, 10, 0.5, 0.5, 1, center_offset=[0, -0.000, 0])
BLA = Rec_off.calc_B_field_coil(Coordsz)
BLA2 = Rec_off_2.calc_B_field_coil(Coordsz)
FIELD_off = BLA + BLA2
Rec = RectangularCoil( 86 / 2, 142.2, 206.4, "z", 10, 10, 0.5, 0.5, 1)
Rec2 = RectangularCoil(-86 / 2, 142.2, 206.4, "z", 10, 10, 0.5, 0.5, 1)
BLA3 = Rec.calc_B_field_coil(Coordsx)
BLA4 = Rec2.calc_B_field_coil(Coordsx)
FIELD = BLA3 + BLA4
homx = RectangularCoil.homogenity_1D(FIELD[0], Coordsy[1], 0.001)
homy = RectangularCoil.homogenity_1D(FIELD[1], Coordsy[1], 0.001)
homz = RectangularCoil.homogenity_1D(FIELD[2], Coordsy[1], 0.001)
#plt.plot(Coordsx[0], FIELD[0], label=r"$\Delta$B_x")
#plt.plot(Coordsx[0], FIELD[1], label=r"$\Delta$B_y")
#plt.plot(Coordsx[0], FIELD[2], label=r"$\Delta$B_z")
plt.plot(Coordsx[0], FIELD[0], label="x_off", alpha=0.6)
plt.plot(Coordsx[0], FIELD[1], label="y_off", alpha=0.6)
plt.plot(Coordsx[0], FIELD[2], label="z_off", alpha=0.6)
plt.legend()
plt.show()
"""
Field = RectangularCoil.HHCoil_calc_B_field(Coordsx, 78.4, 129.4, 206.4, "z", 10, 10, 0.5, 0.5, 1)
plt.plot(Coordsx[0], Field[0], label="x2")
plt.plot(Coordsx[0], Field[1], label="2y")
plt.plot(Coordsx[0], Field[2], label="z2")
plt.show()
Rec.cooling(1,20)
t = np.linspace(0, 12, 1000)
L = Rec.Total_Inductance_Abbot()
L2 = Rec.Total_Inductance_Rover()
R = RectangularCoil.resistance(Rec, T=20)
I = RectangularCoil.I_time_response(U=4.82, R=R, L=L, t=t)
I2 = RectangularCoil.I_time_response(U=4.82, R=R, L=L2, t=t)
print("Inductance Abbott: ",L,"\n inductance rosa: ",L2,R)
print("Resistance: ",R)
tau = L/R
tau2 = L2/R
plt.figure(1,figsize=[8,6])
#plt.title("Time Response for a=145, b=207, d=87.3 with 10x10 windings and 6 Volts")
plt.plot(t, I, label=r"Abbott")#: time constant $ \tau $"+f" = {tau:.2e} s")
plt.plot(t, I2, label=r"Rover")#: time constant $ \tau $"+f" = {tau:.2e} s")
plt.hlines(1/np.sqrt(2), t[0],t[-1], color="tab:green")
plt.legend(fontsize=12)
plt.ylabel("Current I [A]",fontsize=12)
plt.xlabel("Time t [ms]",fontsize=12)
plt.xlim(t[0],t[-1])
plt.ylim(0,1.1)
plt.grid()
plt.show()"""
"""
axis = np.linspace(-0.001, 0.001, 1000)
zeros = np.ones(np.shape(axis)) * 0
Coordsx = np.array([axis, zeros, zeros])
Coordsy = np.array([zeros, axis, zeros])
Coordsz = np.array([zeros, zeros, axis])
Bxaxis = RectangularCoil.HHCoil_calc_B_field(Coordsx, 86, 142.2, 206.4, "z", 10, 10, 0.546, 0.546, 1)
Byaxis = RectangularCoil.HHCoil_calc_B_field(Coordsy, 86, 142.2, 206.4, "z", 10, 10, 0.546, 0.546, 1)
Bzaxis = RectangularCoil.HHCoil_calc_B_field(Coordsz, 86, 142.2, 206.4, "z", 10, 10, 0.546, 0.546, 1)
#BDiagaxis = RectangularCoil.HHCoil_calc_B_field(CoordsDiag, 78.4, 129.4, 206.4, "z", 10, 10, 0.546, 0.546, 1)
Btot = RectangularCoil.HHCoil_calc_total_B_field(Coordsz, 78.4, 129.4, 206.4, "z", 10, 10, 0.546, 0.546, 1)
#hxaxis = RectangularCoil.homogenity_1D(Bxaxis[0], Coordsx[0], 0.001, perc=False)
#hyaxis = RectangularCoil.homogenity_1D(Byaxis[2], Coordsy[1], 0.001, perc=False)
#hzaxis = RectangularCoil.homogenity_1D(Bzaxis[2], Coordsz[2], 0.001, perc=False)
#H = np.array([hxaxis,hyaxis,hzaxis])
plt.plot(axis, Bxaxis[2], label="B_z along x-axis")
plt.plot(axis, Byaxis[2], label="B_z along y-axis")
plt.plot(axis, Bzaxis[2], label="B_z along z-axis")
plt.legend()
plt.grid()
plt.ylabel(r"$B(\vec{r})$ [G]")
plt.xlabel("Position [m]")
plt.show()"""
main()