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

Several analysis

Browse Source
This commit is contained in:
Joschka Schöner 2022-09-02 13:30:37 +02:00
parent dd2a0f93b3
commit 7d1180ed95
34 changed files with 1772 additions and 129 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Coil_geometry/00_Pre_FINAL_coil_geometry.py
View File

@ -26,7 +26,7 @@ HH_Coil = BC.BCoil(HH = 1, distance = 54, radius = 48, layers = 8, windings = 8,
print(f"HH N = {HH_Coil.get_N()}")
I = 1.33 # 64 / HH_Coil.get_N() * 1.25
I = 1.33 # 64 / AHH_Coil.get_N() * 1.25
print(HH_Coil.resistance(20))
print(HH_Coil.induct_perry())

8
Coil_geometry/AHH/08_FINAL_AHH.py
View File

@ -29,14 +29,14 @@ print(f"HH N = {HH_Coil.get_N()}")
I = 64 / HH_Coil.get_N() * 1.25
# set radius plus distance
#HH_Coil.set_R_outer(50.5 - HH_Coil.get_tot_wire_width()*1e3)
#AHH_Coil.set_R_outer(50.5 - AHH_Coil.get_tot_wire_width()*1e3)
HH_Coil.set_R_inner(45.6)
HH_Coil.set_d_min(45.8+2*0.8)
# HH_Coil.B_quick_plot(I)
# HH_Coil.B_curv_quick_plot(I)
# HH_Coil.plot_raster()
# AHH_Coil.B_quick_plot(I)
# AHH_Coil.B_curv_quick_plot(I)
# AHH_Coil.plot_raster()
HH_Coil.print_info()
D_max = 2 * (HH_Coil.get_R_inner()*1e3 - 1) * np.tan(np.radians(41.11))

4
Coil_geometry/Additional coils/00_Test_coil_around viewport.py
View File

@ -67,6 +67,6 @@ print(f"Diff B 15 mm: {Bz[30000] - Bz[15000]}, relative: {(Bz[30000] - Bz[15000]
HH_Coil.cooling(I_current,25)
HH_Coil.print_info()
#HH_Coil.B_curv_quick_plot(I_current)
#HH_Coil.plot_raster()
#AHH_Coil.B_curv_quick_plot(I_current)
#AHH_Coil.plot_raster()

67
Coil_geometry/HH/04_Iterative_Testing.py
View File

@ -11,8 +11,8 @@ from src import B_field_calculation as bf
from src import coil_class as BC
from IPython import get_ipython
get_ipython().run_line_magic('matplotlib', 'qt')
#from IPython import get_ipython
#get_ipython().run_line_magic('matplotlib', 'qt')
#get_ipython().run_line_magic('matplotlib', 'inline')
x = np.linspace(-1, 1, 11)
@ -22,7 +22,7 @@ I_current = 5*16
HH_Coil = HH_Coil_comp = BC.BCoil(HH = 1, distance = 54 ,radius = 37,layers = 1, windings = 1,wire_width = 8, wire_height = 8)
HH_Coil.set_R_outer(49.3)
HH_Coil.set_d_min(49.8)
HH_Coil.set_d_min(47.4)
HH_Coil.print_info()
Bz, Bx = HH_Coil.B_field(I_current, x, z, raster = 50)
@ -31,16 +31,17 @@ HH_Coil.cooling(I_current,30)
print(f"B_z(0) = {Bz[1]:.2f} G")
print(f"B_z_curvature(0) = {Bz_curv[1]:.4f} G/cm^2")
# %%
B = []
Curv = []
array_width = np.arange(0.2,11,0.1)
array_width = np.arange(2.5,6,0.1)
#array_width = [5.7]
for width in array_width:
height = 20/width
height = 16/width
HH_Coil = HH_Coil_comp = BC.BCoil(HH = 1, distance = 54 ,radius = 37,layers = 1, windings = 1,wire_width = width, wire_height = height)
HH_Coil.set_R_outer(49.3)
HH_Coil.set_d_min(49.8)
#HH_Coil.print_info()
HH_Coil.set_R_outer(50)
HH_Coil.set_d_min(47.4)
#AHH_Coil.print_info()
Bz, Bx = HH_Coil.B_field(I_current, x, z, raster = 30)
Bz_curv = BC.BCoil.curv(Bz, z)
@ -57,3 +58,53 @@ plt.ylabel("curvature")
plt.xlabel("total width [mm]")
plt.show()
# %%
HH_Coil = HH_Coil_comp = BC.BCoil(HH = 1, distance = 54 ,radius = 37,layers = 1, windings = 1,wire_width = 8, wire_height = 8)
HH_Coil.set_R_outer(49.3)
HH_Coil.set_d_min(47.4)
#set up axis
x = np.linspace(-15, 15, 30001)
z = np.linspace(-15, 15, 30001)
# New coil
Wire_1 = [0.5, 0.568]
#Wire_1 = [0.45, 0.514]
#I_current = 0.94
HH_Coil = HH_Coil_comp = BC.BCoil(HH = 1, distance = 54 ,radius = 37,layers = 1, windings = 1,wire_width = 4, wire_height = 4)
HH_Coil.set_R_outer(50)
HH_Coil.set_d_min(47.4)
HH_Coil.print_info()
R = HH_Coil.resistance(22.5)
print(f"U = {1 * R}")
I_current = 55
# 0.4 to get from +-30300
HH_Coil.print_info()
#Bz, Bx = AHH_Coil.B_field(I_current, x, z, raster = 7)
Bz, Bx_tot = HH_Coil.B_tot_along_axis(I_current, x, z, raster = 7)
Bz_curv = BC.BCoil.curv(Bz, z)
#AHH_Coil.cooling(I_current,28)
print(f"Bz(0) = {Bz[15000]} G")
print(f"B_z_curvature(0) = {Bz_curv[15000]:.10f} G/cm^2")
#print(f"Bz(1 μm) = {Bz[15001]}")
#print(f"Bz(1 mm) = {Bz[16000]}")
print(f"Diff B +/- 1 μm: {Bz[15001] - Bz[15000]}, relative: {(Bz[15001] - Bz[15000])/Bz[15000]}")
print(f"Diff B +/- 0.5 mm: {Bz[15500] - Bz[15000]}, relative: {(Bz[15500] - Bz[15000])/Bz[15000]}")
print(f"Diff B +/- 1 mm: {Bz[16000] - Bz[15000]}, relative: {(Bz[16000] - Bz[15000])/Bz[15000]}")

12
Coil_geometry/HH/05_try_diff_geometry_HH1.py
View File

@ -38,14 +38,14 @@ print(f"B_z_curvature(0) = {Bz_curv[150]:.4f} G/cm^2")
print(x[500])
# I_current = 5*16
# HH_Coil = BC.BCoil(HH = 1, distance = 54 ,radius = 48 , layers = 1, windings = 1, wire_height = 10, wire_width = 6)
# HH_Coil.set_R_outer(49.3)
# HH_Coil.set_d_min(49.8)
# AHH_Coil = BC.BCoil(HH = 1, distance = 54 ,radius = 48 , layers = 1, windings = 1, wire_height = 10, wire_width = 6)
# AHH_Coil.set_R_outer(49.3)
# AHH_Coil.set_d_min(49.8)
# HH_Coil.print_info()
# Bz, Bx = HH_Coil.B_field(I_current,x,z,raster = 50)
# AHH_Coil.print_info()
# Bz, Bx = AHH_Coil.B_field(I_current,x,z,raster = 50)
# Bz_curv = BC.BCoil.curv(Bz, z)
# HH_Coil.cooling(I_current)
# AHH_Coil.cooling(I_current)
# print(f"Bz(0) = {Bz[150]:.2f} G")
# print(f"B_z_curvature(0) = {Bz_curv[150]:.4f} G/cm^2")

2
Coil_geometry/HH/07_02_testing B_tot.py
View File

@ -30,7 +30,7 @@ HH_Coil.set_R_outer(49.3)
HH_Coil.set_d_min(49.8)
HH_Coil.print_info()
#Bz, Bx = HH_Coil.B_field(I_current,x,z,raster = 10)
#Bz, Bx = AHH_Coil.B_field(I_current,x,z,raster = 10)
B_tot_z, B_tot_x = HH_Coil.B_tot_along_axis(I_current, x, z,raster = 8)
Bz_curv = BC.BCoil.curv(B_tot_z, z)

8
Coil_geometry/HH/11_Final_HH.py
View File

@ -32,25 +32,25 @@ print(f"U = {1 * R}")
I_current = 64 / HH_Coil.get_N() * 1.25
#set radius plus distance
#HH_Coil.set_R_outer(50.5 - HH_Coil.get_tot_wire_width()*1e3 - 0.5)
#AHH_Coil.set_R_outer(50.5 - AHH_Coil.get_tot_wire_width()*1e3 - 0.5)
HH_Coil.set_R_outer(49.93)
additional_space = 0.3
HH_Coil.set_R_outer(49.93-additional_space)
HH_Coil.set_R_inner(45.6)
#HH_Coil.set_R_inner(45.9-0.1)
#AHH_Coil.set_R_inner(45.9-0.1)
# 0.4 to get from +-30300
HH_Coil.set_d_min(47.15+0.4+ 2*additional_space)
print(f"height = {HH_Coil.get_coil_height()}")
HH_Coil.print_info()
#Bz, Bx = HH_Coil.B_field(I_current, x, z, raster = 7)
#Bz, Bx = AHH_Coil.B_field(I_current, x, z, raster = 7)
Bz, B_tot_x = HH_Coil.B_tot_along_axis(I_current, x, z, raster = 7)
Bz_curv = BC.BCoil.curv(Bz, z)
#HH_Coil.cooling(I_current,28)
#AHH_Coil.cooling(I_current,28)
print(f"Bz(0) = {Bz[15000]} G")
print(f"B_z_curvature(0) = {Bz_curv[15000]:.10f} G/cm^2")

2
Coil_geometry/HH/12_Final_Plotting.py
View File

@ -37,7 +37,7 @@ print(HH_Coil.resistance(22))
print(HH_Coil.induct_perry())
HH_Coil.print_info()
#Bz, Bx = HH_Coil.B_field(I_current,x,z,raster = 10)
#Bz, Bx = AHH_Coil.B_field(I_current,x,z,raster = 10)
B_tot_z, B_tot_x = HH_Coil.B_tot_along_axis(I_current, x, z,raster = 8)
Bz_curv = BC.BCoil.curv(B_tot_z, z)

10
Coil_geometry/HH/15_new_HH_coil_iterative_winding_decision.py
View File

@ -30,19 +30,19 @@ def mu_it(x_pos):
Wires = [[0.45, 0.514],[0.475, 0.543],[0.5, 0.568]]
for i in range(0,2):
for i in [2]:
Wire_1 = Wires[i]
print(f"Wire core = {Wire_1[0]} mm:")
print(" ")
for ll in [8,10]:
for ww in [8,9]:
for ll in [6,10]:
for ww in [7,9]:
print(f"layers = {ll}, windings = {ww}")
Coil_1 = BC.BCoil(HH = 1, distance = 54, radius = 48, layers = ll, windings = ww, wire_height = Wire_1[0], wire_width = Wire_1[0], insulation_thickness= (Wire_1[1] - Wire_1[0])/2, is_round = True, winding_scheme= 2)
#set radius plus distance
Coil_1.set_R_outer(50.5-Coil_1.get_tot_wire_width()*1e3)
Coil_1.set_R_outer(49.5)
Coil_1.set_d_min(47.15)
Coil_1.set_d_min(47.4)
#Coil_1.print_info()

10
Coil_geometry/New intermediate pair of coils/01_Coil_calculations.py
View File

@ -28,7 +28,7 @@ AHH_Coil = BC.BCoil(HH=-1, distance=100.336, radius=85.016, layers=10, windings=
winding_scheme=2)
#%%
AHH_Coil.print_basic_info()
I = 1
I = 0.4
AHH_Coil.max_gradient(I)
HH_Coil.max_field(I)
@ -42,11 +42,11 @@ AHH_Coil.print_basic_info()
AHH_Coil.plot_raster()
I = 0.75
#HH_Coil.B_quick_plot(I)
#HH_Coil.B_curv_quick_plot(I)
#AHH_Coil.B_quick_plot(I)
#AHH_Coil.B_curv_quick_plot(I)
# Power
#HH_Coil.cooling(I, 23)
#AHH_Coil.cooling(I, 23)
#AHH_Coil.cooling(I, 23)
print(f"resistance = {HH_Coil.resistance(23)} Ohm")
@ -95,5 +95,5 @@ I = 1
# AHH_Coil.B_curv_quick_plot(I, nr_points = scale)
# HH_Coil.print_info()
# AHH_Coil.print_info()
# AHH_Coil.print_info()

4
Coil_geometry_AHH/Optimum_distance.py
View File

@ -3,6 +3,8 @@ import numpy as np
import matplotlib
#matplotlib.use('Qt5Agg')
from src import coil_class as BC
# %%
HH_Coil = BC.BCoil(HH = 1, distance = 54, radius = 48, layers = 8, windings = 8, wire_height = 0.5,
wire_width = 0.5, insulation_thickness = (0.546-0.5)/2, is_round = True,
winding_scheme= 2)
@ -22,6 +24,8 @@ d = AHH_Coil.get_radius()
print(d)
AHH_Coil.set_d(d + 5)
# %%
z = np.linspace(-50,50,500)
I = 1
Bz, Bx = AHH_Coil.B_field(I,z,z)

33
Coil_geometry_AHH/Theo_Opt_Dist.py Normal file
View File

@ -0,0 +1,33 @@
import matplotlib.pyplot as plt
import numpy as np
import matplotlib
#matplotlib.use('Qt5Agg')
from src import coil_class as BC
# %%
z = np.linspace(-10,10,100)
x = np.linspace(-10,10,100)
I = 1
for i in range(0, 10):
AHH_Coil = BC.BCoil(HH = -1, distance=10+i, radius = 10, layers = 1, windings=1,
wire_height = 0.1, wire_width=0.1, insulation_thickness=0,
is_round = True)
Bz, Bx = AHH_Coil.B_field(I,z,z)
Bz_grad = BC.BCoil.grad(Bz,z)
Bx_grad = BC.BCoil.grad(Bx,x)
#plt.plot(z,Bz)
plt.plot(z, Bx_grad, label=f"d={10+i}mm")
plt.legend()
#plt.plot(z, Bz_curv)
#plt.plot(z, Bz_4)
plt.show()

6
Cooling/05_POWER_estimation_water_cooling.py
View File

@ -19,12 +19,12 @@ def Q_heat(flow,d_T):
def main():
d_T = 1.5
d_T = 2
flow = 0.3/5#/5 #m/s
flow = 1#/5 #m/s
#flow *=2
print(f"flow = {flow}m/s")
V_t = 4.8 * 3.2 * 1e-6 * flow
V_t = 4.7 * 3.2 * 1e-6 * flow
print(f"Volume rate = {V_t * 1e6} mL/s")
m = cs.water_dens * V_t
Q = m * cs.water_c_p * d_T

97
FINAL_Coil/FINAL_COIL_configuration.py
View File

@ -10,25 +10,33 @@ HH_Coil = BC.BCoil(HH = 1, distance = 54, radius = 48, layers = 8, windings = 8,
winding_scheme= 2)
HH_Coil.set_R_inner(45.6)
HH_Coil.set_d_min(2*24.075)
HH_Coil.set_d_min(47.4)
HH_Coil.print_info()
print(f"inductivity = {HH_Coil.induct_perry()*2} H")
AHH_Coil = BC.BCoil(HH = -1, distance = 54, radius = 48, layers = HH_Coil.get_layers, windings=2 * HH_Coil.get_windings,
wire_height = 0.5, wire_width=0.5, insulation_thickness=(0.546-0.5)/2,
is_round = True, winding_scheme= 2)
is_round = True, winding_scheme= 2, red_N=1)
AHH_Coil.set_R_inner(45.6)
AHH_Coil.set_d_min(HH_Coil.get_zmax()*2 * 1e3 + 4)
AHH_Coil.set_d_max(77.6)
AHH_Coil.print_info()
R = HH_Coil.resistance(25)
I = 5
AHH_Coil.cooling(I, 30)
# %%
AHH_Coil.plot_raster()
print(AHH_Coil.get_N())
# %%
I = 1
AHH_Coil.cooling(I, 22.5)
AHH_Coil.max_gradient(I)
I_2 = 0.92
HH_Coil.cooling(I_2,22.5)
HH_Coil.max_field(I_2)
# %%
AHH_Coil.B_quick_plot(I)
@ -37,19 +45,84 @@ AHH_Coil.print_info()
print("All values for pair in series")
print(f"inductivity AHH: {AHH_Coil.induct_perry()*2*1e3} mH" )
print(f"inductivity HH: {HH_Coil.induct_perry() * 2*1e3} mH" )
print(f"resistance AHH: {2*AHH_Coil.resistance(22)} Ω")
print(f"resistance HH: {2*HH_Coil.resistance(22)} Ω")
print(f"resistance AHH: {AHH_Coil.resistance(22.5)*2} Ω")
print(f"resistance HH: {HH_Coil.resistance(22.5)*2} Ω")
# %%
I = 1.3
I = 1
AHH_Coil.cooling(I,22)
AHH_Coil.max_gradient(I)
I_HH = 1.4
HH_Coil.cooling(I,1.2)
I_HH = 1
HH_Coil.cooling(I,22)
HH_Coil.max_field(I)
#%%
# field quality
HH_Coil.B_quality(z_pos=1)
HH_Coil.B_quality(z_pos=19)
# %%
AHH_Coil.Grad_quality(z_pos=1)
#%% gradient quality
x = np.linspace(-20, 20, 40001)
z = np.linspace(-20, 20, 40001)
B_z, B_x_tot = AHH_Coil.B_field(1, x, z, raster = 7)
Bz = BC.BCoil.grad(B_z, z)
#AHH_Coil.cooling(I_current,28)
print(f"Bz_grad(0) = {Bz[15000]} G")
#print(f"B_z_curvature(0) = {Bz_curv[15000]:.10f} G/cm^2")
#print(f"Bz(1 μm) = {Bz[15001]}")
#print(f"Bz(1 mm) = {Bz[16000]}")
print(f"Diff B +/- 1 μm: {Bz[15001] - Bz[15000]}, relative: {(Bz[15001] - Bz[15000])/Bz[15000]}")
print(f"Diff B +/- 0.5 mm: {Bz[15500] - Bz[15000]}, relative: {(Bz[15500] - Bz[15000])/Bz[15000]}")
print(f"0 = {z[20000]}, 1 = {z[21000]}, 19 = {z[39000]}")
print(f"Diff Grad B +/- 1 mm: relative: {(Bz[20000] - Bz[21000])/Bz[20000]}")
print(f"Diff Grad B +/- 19 mm: relative: {(Bz[20000] - Bz[39000])/Bz[20000]}")
# %%
plt.plot(z,Bz)
plt.show()
# %%
rho = BC.BCoil.resistivity_copper(50)
l = 48e-3*2*np.pi * 64
A = 0.25e-3**2 *np.pi
R_22 = BC.BCoil.resistivity_copper(16) *l/A
R_50 = BC.BCoil.resistivity_copper(50)*l/A
print(R_22)
print(R_50)
print((R_50-R_22)/R_50)
print(BC.BCoil.resistivity_copper(28))
# SF coil
# %%
SF_Coil = BC.BCoil(HH = 1, distance = 54, radius = 48, layers = 1, windings = 1, wire_height = 0.5,
wire_width = 0.5, insulation_thickness = (0.546-0.5)/2, is_round = True,
winding_scheme= 2)
SF_Coil.set_R_inner(38.4)
SF_Coil.set_d_min(47.4)
SF_Coil.print_info()
print(f"inductivity = {SF_Coil.induct_perry()*2} H")
print(f"resistance SF: {SF_Coil.resistance(22.5)*2} Ω")
SF_Coil.max_field(1)
SF_Coil.B_quality()
SF_Coil.B_quality( z_pos=19)
SF_Coil.cooling(1, 22.5)
print(SF_Coil.power(1, 22.5))

7
FINAL_Coil/Final_Field_quality.py
View File

@ -41,13 +41,13 @@ I_current = 64 / HH_Coil.get_N() * 1.25
# 0.4 to get from +-30300
HH_Coil.print_info()
#Bz, Bx = HH_Coil.B_field(I_current, x, z, raster = 7)
#Bz, Bx = AHH_Coil.B_field(I_current, x, z, raster = 7)
Bz, B_x_tot = HH_Coil.B_tot_along_axis(I_current, x, z, raster = 7)
Bz_curv = BC.BCoil.curv(Bz, z)
#HH_Coil.cooling(I_current,28)
#AHH_Coil.cooling(I_current,28)
print(f"Bz(0) = {Bz[15000]} G")
print(f"B_z_curvature(0) = {Bz_curv[15000]:.10f} G/cm^2")
@ -62,6 +62,9 @@ print(f"Diff B +/- 0.5 mm: {Bz[15500] - Bz[15000]}, relative: {(Bz[15500] - Bz[1
print(f"Diff B +/- 1 mm: {Bz[16000] - Bz[15000]}, relative: {(Bz[16000] - Bz[15000])/Bz[15000]}")
# %%
HH_Coil.B_quality()
#print(f"Diff B +/- 15 mm: {Bz[30000] - Bz[15000]}, relative: {(Bz[30000] - Bz[15000])/Bz[15000]}")

113
FINAL_Coil/Final_Gradient_quality.py Normal file
View File

@ -0,0 +1,113 @@
# -*- coding: utf-8 -*-
"""
Created on Mon Aug 23 17:40:37 2021
@author: Joschka
"""
import matplotlib.pyplot as plt
import numpy as np
#from src import B_field_calculation as bf
from src import coil_class as BC
#from IPython import get_ipython
#get_ipython().run_line_magic('matplotlib', 'qt')
#get_ipython().run_line_magic('matplotlib', 'inline')
#set up axis
x = np.linspace(-15, 15, 30001)
z = np.linspace(-15, 15, 30001)
# New coil
Wire_1 = [0.5, 0.568]
#Wire_1 = [0.45, 0.514]
#I_current = 0.94
AHH_Coil = BC.BCoil(HH = -1, distance = 69.622, radius = 47.528, layers = 8, windings=16,
wire_height = 0.5, wire_width=0.5, insulation_thickness=(0.546-0.5)/2,
is_round = True, winding_scheme= 2)
AHH_Coil.print_info()
R = AHH_Coil.resistance(22.5)
print(f"U = {1 * R}")
I_current = 1
# 0.4 to get from +-30300
AHH_Coil.print_info()
#Bz, Bx = AHH_Coil.B_field(I_current, x, z, raster = 7)
B_z, B_x_tot = AHH_Coil.B_field(I_current, x, z, raster = 7)
Bz = BC.BCoil.grad(B_z, z)
Bz_curv = BC.BCoil.curv(Bz, z)
#AHH_Coil.cooling(I_current,28)
print(f"Bz_grad(0) = {Bz[15000]} G")
#print(f"B_z_curvature(0) = {Bz_curv[15000]:.10f} G/cm^2")
#print(f"Bz(1 μm) = {Bz[15001]}")
#print(f"Bz(1 mm) = {Bz[16000]}")
print(f"Diff B +/- 1 μm: {Bz[15001] - Bz[15000]}, relative: {(Bz[15001] - Bz[15000])/Bz[15000]}")
print(f"Diff B +/- 0.5 mm: {Bz[15500] - Bz[15000]}, relative: {(Bz[15500] - Bz[15000])/Bz[15000]}")
print(f"Diff B +/- 1 mm: {Bz[16000] - Bz[15000]}, relative: {(Bz[16000] - Bz[15000])/Bz[15000]}")
#print(f"Diff B +/- 15 mm: {Bz[30000] - Bz[15000]}, relative: {(Bz[30000] - Bz[15000])/Bz[15000]}")
"""
plt.figure(300)
#Field plot
##########################
plt.subplot(2,1,1)
plt.plot(z,Bz,linestyle = "solid", label = r"$Bz along z-axis")
plt.plot(z,B_tot_z, linestyle = "dashed", label = "New B_tot along z-axis")
#plt.plot(x,B_tot_x, label = "B_tot along x-axis")
#plt.plot(z,B_z_comp,linestyle = "solid", label = r"$B_{z,1}$, d = 54 mm, R = 48.8 mm, I = 5 A, 4 x 4")
#plt.plot(z,B_tot,linestyle = "solid", label = r"$B_{z,1} + B_{z,2}$")
#plt.xlim(-0.01,0.01)
plt.title("B-field" )
plt.ylabel(r"$Bz$ [G]")
plt.xlabel("z-axis [mm]")
plt.legend()#bbox_to_anchor=(1.05, 1), loc='upper left')
plt.subplot(2,1,2)
plt.plot(z,Bz_curv,linestyle = "solid", label = r"$\nabla_z^2 B_{ref}$, d = 44 mm, R = 44 mm")
#plt.plot(z,B_z_comp_curv,linestyle = "solid", label = r"$\nabla_z^2 B_{z,1}$, d = 54 mm, R = 48.8 mm, I = 5 A")
#plt.plot(z,B_tot_curv,linestyle = "solid", label = r"$\nabla_z^2 B_{z,1} + B_{z,2}$")
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()#bbox_to_anchor=(1.05, 1), loc='upper left')
#plt.savefig("output/first_compensation_idea.png")
plt.show()
"""
"""
AHH ############################################################################
###############################################################################
###############################################################################
"""

168
FINAL_Coil/Test_different_diameters.py Normal file
View File

@ -0,0 +1,168 @@
# -*- coding: utf-8 -*-
"""
Created on Mon Aug 23 17:40:37 2021
@author: Joschka
"""
import matplotlib.pyplot as plt
import numpy as np
#from src import B_field_calculation as bf
from src import coil_class as BC
#from IPython import get_ipython
#get_ipython().run_line_magic('matplotlib', 'qt')
#get_ipython().run_line_magic('matplotlib', 'inline')
# %%
#set up axis
x = np.linspace(-15, 15, 30001)
z = np.linspace(-15, 15, 30001)
# New coil
Wire_1 = [0.5, 0.568]
#Wire_1 = [0.45, 0.514]
#I_current = 0.94
HH_Coil = BC.BCoil(HH = 1, distance = 54, radius = 48, layers = 8, windings = 8, wire_height = 0.5,
wire_width = 0.5, insulation_thickness = (0.546-0.5)/2, is_round = True,
winding_scheme= 2)
HH_Coil.set_R_inner(45.6)
HH_Coil.set_d_min(2*24.075)
HH_Coil.print_info()
R = HH_Coil.resistance(22.5)
print(f"U = {1 * R}")
# %%
I_current = 1
# 0.4 to get from +-30300
#Bz, Bx = AHH_Coil.B_field(I_current, x, z, raster = 7)
Bz, B_x_tot = HH_Coil.B_field(I_current, x, z, raster = 7)
Bz_as, B_x_as = HH_Coil.B_field_asym(I_current, x, z, raster = 7)
Bz_curv = BC.BCoil.curv(Bz, z)
Bz_curv_as = BC.BCoil.curv(Bz_as, z)
#AHH_Coil.cooling(I_current,28)
# %%
print("Symmetric coil configuration")
print(f"Bz(0) = {Bz[15000]} G")
print(f"B_z_curvature(0) = {Bz_curv[15000]:.10f} G/cm^2")
#print(f"Bz(1 μm) = {Bz[15001]}")
#print(f"Bz(1 mm) = {Bz[16000]}")
print(f"Diff B +/- 1 μm: {Bz[15001] - Bz[15000]}, relative: {(Bz[15001] - Bz[15000])/Bz[15000]}")
print(f"Diff B +/- 0.5 mm: {Bz[15500] - Bz[15000]}, relative: {(Bz[15500] - Bz[15000])/Bz[15000]}")
print(f"Diff B +/- 1 mm: {Bz[16000] - Bz[15000]}, relative: {(Bz[16000] - Bz[15000])/Bz[15000]}")
print("")
print("Asymmetric coil configuration (one with 200 μm smaller diameter)")
print(f"Bz(0) = {Bz_as[15000]} G")
print(f"B_z_curvature(0) = {Bz_curv_as[15000]:.10f} G/cm^2")
#print(f"Bz(1 μm) = {Bz[15001]}")
#print(f"Bz(1 mm) = {Bz[16000]}")
print(f"Diff B + 1 μm: {Bz_as[15001] - Bz_as[15000]}, relative: {(Bz_as[15001] - Bz_as[15000])/Bz_as[15000]}")
print(f"Diff B - 1 μm: {Bz_as[14999] - Bz_as[15000]}, relative: {(Bz_as[14999] - Bz_as[15000])/Bz_as[15000]}")
print("")
print(f"Diff B + 0.5 mm: {Bz_as[15500] - Bz_as[15000]}, relative: {(Bz_as[15500] - Bz_as[15000])/Bz_as[15000]}")
print(f"Diff B - 0.5 mm: {Bz_as[14500] - Bz_as[15000]}, relative: {(Bz_as[14500] - Bz_as[15000])/Bz_as[15000]}")
print("")
print(f"Diff B + 1 mm: {Bz_as[16000] - Bz_as[15000]}, relative: {(Bz_as[16000] - Bz_as[15000])/Bz_as[15000]}")
print(f"Diff B - 1 mm: {Bz_as[14000] - Bz_as[15000]}, relative: {(Bz_as[14000] - Bz_as[15000])/Bz_as[15000]}")
print("")
print("Comparison Bx")
print(f"Symmetric Bx = {B_x_tot[15001]}")
print(f"Bx(0) = {B_x_as[15001]}")
print(f"Diff B - 1 mm: {B_x_as[14000] - B_x_as[15000]}, relative: {(B_x_as[14000] - B_x_as[15000])/B_x_as[15000]}")
#print(f"Diff B +/- 15 mm: {Bz[30000] - Bz[15000]}, relative: {(Bz[30000] - Bz[15000])/Bz[15000]}")
_
# %%
HH_Coil = BC.BCoil(HH = 1, distance = 54, radius = 48, layers = 8, windings = 8, wire_height = 0.5,
wire_width = 0.5, insulation_thickness = (0.546-0.5)/2, is_round = True,
winding_scheme= 2)
HH_Coil.set_R_inner(45.6-0.1)
HH_Coil.set_d_min(2*24.075)
HH_Coil.print_info()
# %%
I_current = 1
# 0.4 to get from +-30300
#Bz, Bx = AHH_Coil.B_field(I_current, x, z, raster = 7)
Bz_2, B_x_2 = HH_Coil.B_field(I_current, x, z, raster = 7)
Bz_curv_2 = BC.BCoil.curv(Bz_2, z)
print("")
print(f"Diff B +/- 1 mm: {Bz_2[16000] - Bz_2[15000]}, relative: {(Bz_2[16000] - Bz_2[15000])/Bz_2[15000]}")
"""
plt.figure(300)
#Field plot
##########################
plt.subplot(2,1,1)
plt.plot(z,Bz,linestyle = "solid", label = r"$Bz along z-axis")
plt.plot(z,B_tot_z, linestyle = "dashed", label = "New B_tot along z-axis")
#plt.plot(x,B_tot_x, label = "B_tot along x-axis")
#plt.plot(z,B_z_comp,linestyle = "solid", label = r"$B_{z,1}$, d = 54 mm, R = 48.8 mm, I = 5 A, 4 x 4")
#plt.plot(z,B_tot,linestyle = "solid", label = r"$B_{z,1} + B_{z,2}$")
#plt.xlim(-0.01,0.01)
plt.title("B-field" )
plt.ylabel(r"$Bz$ [G]")
plt.xlabel("z-axis [mm]")
plt.legend()#bbox_to_anchor=(1.05, 1), loc='upper left')
plt.subplot(2,1,2)
plt.plot(z,Bz_curv,linestyle = "solid", label = r"$\nabla_z^2 B_{ref}$, d = 44 mm, R = 44 mm")
#plt.plot(z,B_z_comp_curv,linestyle = "solid", label = r"$\nabla_z^2 B_{z,1}$, d = 54 mm, R = 48.8 mm, I = 5 A")
#plt.plot(z,B_tot_curv,linestyle = "solid", label = r"$\nabla_z^2 B_{z,1} + B_{z,2}$")
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()#bbox_to_anchor=(1.05, 1), loc='upper left')
#plt.savefig("output/first_compensation_idea.png")
plt.show()
"""
"""
AHH ############################################################################
###############################################################################
###############################################################################
"""

2
Noise/00_Simple_testing.py
View File

@ -29,7 +29,7 @@ print(diff)
Bz2, Bx = HH_Coil.B_field(I+ diff, x, z)
Bz2, Bx = AHH_Coil.B_field(I+ diff, x, z)
print(Bz2[1500]-Bz1[1500])
print(" ")
"""

24
Plots_Erlangen/Plots.py
View File

@ -94,7 +94,7 @@ I = 1
Max_field = HH_Coil.max_field(I)
def I_t_cut(time, coil, I_end, U_0):
I = U_0 / coil.resistance(22) * (1 - np.exp(- time / coil.tau()))
I = U_0 / coil.resistance(22)/2 * (1 - np.exp(- time / coil.tau()))
if I >= I_end:
I = I_end
return I
@ -102,7 +102,7 @@ def I_t_cut(time, coil, I_end, U_0):
def I_current(Coil, I_0, t):
L = Coil.induct_perry()
R = Coil.resistance(22.5)
R = Coil.resistance(22.5)*2
print(f"L={L}")
print(f" R= {R}")
tau = L / R
@ -119,25 +119,27 @@ t = np.linspace(0, 0.002, 10000)
i_to_B = 10.64
fig, ax = plt.subplots(figsize = (11,7))
#fig.suptitle(f"Time response HH-coil: I_max = {I} A --> Max Field = {Max_field:.2f} G \n \n I(t) = U(t) / R * (1 - exp(- R/L * t))")
ax.set_title("Time Response Offset Coil", y = 1.05)
ax.set_title("Expected Time Response Final Offset Coil", y = 1.05)
ax.text(0.6, 5, r'$I(t) = \frac{U(t)}{R} - \frac{L}{R} \cdot \frac{dI(t)}{dt} $', fontsize=34)
# ax.text(0.6, 5, r'$I(t) = \frac{U(t)}{R} - \frac{L}{R} \cdot \frac{dI(t)}{dt} $', fontsize=34)
ax.plot(t * 1e3, i_to_B * I_current(HH_Coil, I, t), label=f"U(t) = 1.5 V", zorder=1, color=my_colors['pastel_blue'])
# ax.plot(t * 1e3, i_to_B * I_current(HH_Coil, I, t), label=f"U(t) = 1.5 V", zorder=1, color=my_colors['pastel_blue'])
U_0 = 28
ax.plot(t * 1e3, i_to_B * I_t_cut_vec(t, HH_Coil, I, U_0), label=f"U(t) regulated via PI feedback loop", zorder=1, color=my_colors['light_green'])
plt.vlines(3.1e-2, 0, 10.64, zorder=2, linestyles=(0, (1.5, 3.06)),color=my_colors['orange'], label='t = 30 μs')
ax.plot(t * 1e3, i_to_B * I_t_cut_vec(t, HH_Coil, I, U_0), zorder=1, color=my_colors['pastel_blue'])
plt.vlines(6.2e-2, 0, 10.64, zorder=2, linestyles=(0, (1.5, 3.06)),color=my_colors['orange'], label='t = 60 μs')
# for scaling in np.arange(2,5,0.5):
# ax.plot(t * 1e3, I_t_exp(t, HH_Coil, I, 15, scaling), label=f"Exponential decay U")
# ax.plot(t * 1e3, I_t_exp(t, AHH_Coil, I, 15, scaling), label=f"Exponential decay U")
ax.set_xlabel("time [ms]")
ax.set_ylabel("Magnetic field [G]")
ax.set_ylim(ylim)
ax.set_xlim(-0.09,2)
ax.set_xlim(-0.009,0.25)
ax.legend()
plt.show()
# %%
# Comparison different number of windings
Coil1 = BC.BCoil(HH=1, distance=54, radius=48, layers=8, windings=8, wire_height=0.5,
@ -195,10 +197,10 @@ ax.text(0.6, 5, r'$I(t) = \frac{U(t)}{R} - \frac{L}{R} \cdot \frac{dI(t)}{dt} $'
ax.plot(t * 1e3, i_to_B * I_current(Coil1, I, t), label=f"U(t) = 1.5 V", zorder=1, color=my_colors['pastel_blue'])
ax.plot(t * 1e3, i_to_B * I_current(Coil2, I , t), label=f"U(t) = 1.5 V", zorder=1, linestyle='- -', color=my_colors['light_green'])
U_0 = 28
#ax.plot(t * 1e3, i_to_B * I_t_cut_vec(t, HH_Coil, I, U_0), label=f"U(t) regulated via PI feedback loop", zorder=1, color=my_colors['light_green'])
#ax.plot(t * 1e3, i_to_B * I_t_cut_vec(t, AHH_Coil, I, U_0), label=f"U(t) regulated via PI feedback loop", zorder=1, color=my_colors['light_green'])
#plt.vlines(3.1e-2, 0, 10.64, zorder=2, linestyles=(0, (1.5, 3.06)),color=my_colors['orange'], label='t = 30 μs')
# for scaling in np.arange(2,5,0.5):
# ax.plot(t * 1e3, I_t_exp(t, HH_Coil, I, 15, scaling), label=f"Exponential decay U")
# ax.plot(t * 1e3, I_t_exp(t, AHH_Coil, I, 15, scaling), label=f"Exponential decay U")
ax.set_xlabel("time [ms]")
ax.set_ylabel("Magnetic field [G]")

111
Plots_Erlangen/Plots_midterm_presentation.py
View File

@ -1,6 +1,7 @@
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import scipy.io.wavfile as wavfile
from src import coil_class as BC
# %%
@ -21,7 +22,6 @@ mpl.rcParams['ytick.major.width'] = 3
mpl.rcParams['ytick.minor.size'] = 10
mpl.rcParams['ytick.minor.width'] = 3
mpl.rcParams.update({'font.size': 22, 'axes.linewidth': 3, 'lines.linewidth': 3})
# %%
@ -65,7 +65,7 @@ my_colors = {'light_green': '#97e144',
'light_blue': '#71C8F4',
'purple': '#7c588c'}
c_field = my_colors['light_green']
c_field = my_colors['orange']
c_grad = my_colors['purple']
fig, ax1 = plt.subplots(figsize=(11, 6))
@ -88,19 +88,30 @@ plt.plot(z, np.abs(AHH_B_grad_z), color=c_grad)
ax2.set_ylim(2.1, 5.5)
plt.show()
# %% import data
samplingratePI, dataPI = wavfile.read(
"C:/Users/Joschka/Desktop/Midterm_presentation/Plot_time_response/SquareWave300HzWithPI.wav")
dataPI = dataPI[:2200000] / (2.5 * 2 ** 16)
samplingrateNoPI, dataNoPI = wavfile.read(
"C:/Users/Joschka/Desktop/Midterm_presentation/Plot_time_response/SquareWave300HzNoPI.wav")
dataNoPI = dataNoPI[:2200000] / (2.5 * 2 ** 16)
# %%
I = 1
Max_field = HH_Coil.max_field(I)
def I_t_cut(time, coil, I_end, U_0):
I = U_0 / coil.resistance(22) * (1 - np.exp(- time / coil.tau()))
if I >= I_end:
I = I_end
return I
def I_current(Coil, I_0, t):
L = Coil.induct_perry()
# L = 300e-6
R = Coil.resistance(22.5)
print(f"L={L}")
@ -110,54 +121,64 @@ def I_current(Coil, I_0, t):
I = I_0 * (1 - np.exp(-R / L * t))
return I
I_t_cut_vec = np.vectorize(I_t_cut)
I_t_cut_vec = np.vectorize(I_t_cut)
fig, ax1 = plt.subplots(figsize=(11, 8))
ylim = (0, 11.5)
t = np.linspace(0, 0.002, 10000)
i_to_B = 10.64
fig, ax = plt.subplots(figsize = (11,7))
#fig.suptitle(f"Time response HH-coil: I_max = {I} A --> Max Field = {Max_field:.2f} G \n \n I(t) = U(t) / R * (1 - exp(- R/L * t))")
ax.set_title("Time Response Offset Coil", y = 1.05)
fig, ax = plt.subplots(figsize=(11, 7))
# fig.suptitle(f"Time response HH-coil: I_max = {I} A --> Max Field = {Max_field:.2f} G \n \n I(t) = U(t) / R * (1 - exp(- R/L * t))")
ax.set_title("Time Response Offset Coil", y=1.05)
ax.text(0.6, 5, r'$I(t) = \frac{U(t)}{R} - \frac{L}{R} \cdot \frac{dI(t)}{dt} $', fontsize=34)
ax.plot(t * 1e3, i_to_B * I_current(HH_Coil, I, t), label=f"U(t) = 1.5 V", zorder=1, color=my_colors['pastel_blue'])
U_0 = 28
ax.plot(t * 1e3, i_to_B * I_t_cut_vec(t, HH_Coil, I, U_0), label=f"U(t) regulated via PI feedback loop", zorder=1, color=my_colors['light_green'])
plt.vlines(3.1e-2, 0, 10.64, zorder=2, linestyles=(0, (1.5, 3.06)),color=my_colors['orange'], label='t = 30 μs')
ax.plot(t * 1e3, i_to_B * I_t_cut_vec(t, HH_Coil, I, U_0), label=f"U(t) regulated via PI feedback loop", zorder=1,
color=my_colors['light_green'])
plt.vlines(3.1e-2, 0, 10.64, zorder=2, linestyles=(0, (1.5, 3.06)), color=my_colors['orange'], label='t = 30 μs')
# for scaling in np.arange(2,5,0.5):
# ax.plot(t * 1e3, I_t_exp(t, HH_Coil, I, 15, scaling), label=f"Exponential decay U")
# ax.plot(t * 1e3, I_t_exp(t, AHH_Coil, I, 15, scaling), label=f"Exponential decay U")
ax.set_xlabel("time [ms]")
ax.set_ylabel("Magnetic field [G]")
ax.set_ylim(ylim)
ax.set_xlim(-0.09,2)
ax.legend()
ax.set_xlim(-0.001, 2)
#ax.legend()
time = np.linspace(0,5030/samplingratePI,5030)*1e3
plt.plot(time,i_to_B *dataPI[5470:10500]/0.1,color="limegreen",label="U(t) regulated via PI feedback loop")#/dataPI[10500])
plt.plot(time,i_to_B *dataNoPI[6004:11034]/0.136,color="royalblue", label = "U(t) = 1.5 V")#/dataNoPI[11034])
plt.xlabel("Time [ms]")
plt.ylabel("Current [A]")
plt.vlines(0.052,0,1.01, color="orange", linestyle = "--", label="t = 52 s$")
#plt.ylim(0,1.05)
#plt.legend()
plt.show()
# %%
# Comparison different number of windings
Coil1 = BC.BCoil(HH=1, distance=54, radius=48, layers=8, windings=8, wire_height=0.5,
wire_width=0.5, insulation_thickness=(0.546 - 0.5)/2, is_round=True,
winding_scheme=2)
wire_width=0.5, insulation_thickness=(0.546 - 0.5) / 2, is_round=True,
winding_scheme=1)
Coil1.set_R_inner(45.6)
Coil1.set_d_min(2 * 24.075)
Coil1.print_info()
factor = 2
Coil2 = BC.BCoil(HH=1, distance=54, radius=48, layers=8//factor, windings=8//factor, wire_height=0.5*factor,
wire_width=0.5*factor, insulation_thickness=(0.546 - 0.5)*factor/2, is_round=True,
winding_scheme=1)
Coil2 = BC.BCoil(HH=1, distance=54, radius=48, layers=8 // factor, windings=8 // factor, wire_height=0.5 * factor,
wire_width=0.5 * factor, insulation_thickness=(0.546 - 0.5) * factor / 2, is_round=True,
winding_scheme=1)
Coil2.set_R_inner(45.6)
Coil2.set_d_min(2 * 24.075)
Coil2.print_info()
full_structure = Coil1.full_raster(100) * 1e3
if Coil1.get_coil_width() > Coil1.get_coil_height():
extension = Coil1.get_coil_width()
@ -171,27 +192,41 @@ plt.scatter(full_structure[:, :, 1], full_structure[:, :, 0], linewidths=0.001)
plt.xlabel("radius [mm]")
plt.ylabel("z position [mm]")
plt.axvline(x=Coil1.get_R_inner() * 1e3 - 0.1, lw=5, color="red")
plt.xlim(45, 50.4)
plt.ylim(23.5, 28.9)
plt.xlim(45, 52)
plt.ylim(23.5, 30.5)
plt.show()
# plt.figure(1,figsize=[8,5])
time = np.linspace(0, 5030 / samplingratePI, 5030) * 1e3
plt.plot(time, dataPI[5470:10500] / 0.1, color="limegreen",
label="U(t) regulated via PI feedback loop") # /dataPI[10500])
plt.plot(time, dataNoPI[6004:11034] / 0.136, color="royalblue", label="U(t) = 1.5 V") # /dataNoPI[11034])
plt.xlabel("Time [ms]")
plt.ylabel("Current [A]")
plt.vlines(0.052, 0, 1.01, color="orange", linestyle="--", label="t = 52 s$")
plt.ylim(0, 1.05)
plt.legend()
# plt.show()
plt.show()
# %%
Coil2.plot_raster()
factor = Coil1.get_N()//Coil2.get_N()
factor = Coil1.get_N() // Coil2.get_N()
I = 1
I2 = I * factor
Max_field = HH_Coil.max_field(I)
def I_t_cut(time, coil, I_end, U_0):
I = U_0 / coil.resistance(22) * (1 - np.exp(- time / coil.tau()))
if I >= I_end:
I = I_end
return I
def I_current(Coil, I_0, t):
L = Coil.induct_perry()
@ -203,45 +238,47 @@ def I_current(Coil, I_0, t):
I = I_0 * (1 - np.exp(-R / L * t))
return I
I_t_cut_vec = np.vectorize(I_t_cut)
I_t_cut_vec = np.vectorize(I_t_cut)
fig, ax1 = plt.subplots(figsize=(11, 8))
ylim = (0, 11.5)
t = np.linspace(0, 0.002, 10000)
i_to_B = Max_field/I
fig, ax = plt.subplots(figsize = (11,7))
#fig.suptitle(f"Time response HH-coil: I_max = {I} A --> Max Field = {Max_field:.2f} G \n \n I(t) = U(t) / R * (1 - exp(- R/L * t))")
i_to_B = Max_field / I
fig, ax = plt.subplots(figsize=(11, 7))
# fig.suptitle(f"Time response HH-coil: I_max = {I} A --> Max Field = {Max_field:.2f} G \n \n I(t) = U(t) / R * (1 - exp(- R/L * t))")
# ax.set_title("Time Response", y = 1.05)
# ax.text(0.4, 4, r"$ I(t) = \frac{U_{0}}{R} (1 - e^{-\frac{R}{L} t})$", fontsize=34)
ax.plot(t * 1e3, i_to_B * I_current(Coil1, I, t), label=f"N = {Coil1.get_N()}, I$_{{end}}$ = {I} A", zorder=1, color=my_colors['pastel_blue'])
ax.plot(t * 1e3, i_to_B/factor * I_current(Coil2, I2 , t), label=f"N = {Coil2.get_N()}, I$_{{end}}$ = {I2} A", zorder=1, linestyle=(0, (4, 4)), color=my_colors['light_green'])
ax.plot(t * 1e3, i_to_B * I_current(Coil1, I, t), label=f"N = {Coil1.get_N()}, I$_{{end}}$ = {I} A", zorder=1,
color=my_colors['pastel_blue'])
ax.plot(t * 1e3, i_to_B / factor * I_current(Coil2, I2, t), label=f"N = {Coil2.get_N()}, I$_{{end}}$ = {I2} A",
zorder=1, linestyle=(0, (4, 4)), color=my_colors['light_green'])
U_0 = 28
#ax.plot(t * 1e3, i_to_B * I_t_cut_vec(t, HH_Coil, I, U_0), label=f"U(t) regulated via PI feedback loop", zorder=1, color=my_colors['light_green'])
#plt.vlines(3.1e-2, 0, 10.64, zorder=2, linestyles=(0, (1.5, 3.06)),color=my_colors['orange'], label='t = 30 μs')
# ax.plot(t * 1e3, i_to_B * I_t_cut_vec(t, AHH_Coil, I, U_0), label=f"U(t) regulated via PI feedback loop", zorder=1, color=my_colors['light_green'])
# plt.vlines(3.1e-2, 0, 10.64, zorder=2, linestyles=(0, (1.5, 3.06)),color=my_colors['orange'], label='t = 30 μs')
# for scaling in np.arange(2,5,0.5):
# ax.plot(t * 1e3, I_t_exp(t, HH_Coil, I, 15, scaling), label=f"Exponential decay U")
# ax.plot(t * 1e3, I_t_exp(t, AHH_Coil, I, 15, scaling), label=f"Exponential decay U")
ax.set_xlabel("time [ms]")
ax.set_ylabel("Magnetic field [G]")
ax.set_ylim(ylim)
ax.set_xlim(-0.09,2)
ax.set_xlim(-0.09, 2)
ax.legend()
plt.show()
#%%
# %%
Coil1 = BC.BCoil(HH=1, distance=54, radius=48, layers=8, windings=8, wire_height=0.5,
wire_width=0.5, insulation_thickness=(0.546 - 0.5)/2, is_round=True,
winding_scheme=2)
wire_width=0.5, insulation_thickness=(0.546 - 0.5) / 2, is_round=True,
winding_scheme=2)
Coil1.set_R_inner(45.6)
Coil1.set_d_min(2 * 24.075)
Coil1.print_info()
mpl.rcParams.update(mpl.rcParamsDefault)
print(f"Cross_section = {Coil1.get_cross_section()}")
fill_factor = Coil1.get_wire_area() * Coil1.get_N()/Coil1.get_cross_section()
fill_factor = Coil1.get_wire_area() * Coil1.get_N() / Coil1.get_cross_section()
print(f"fill_factor = {fill_factor}")
Coil1.plot_raster()
Coil1.plot_raster()

2
RF_coil/frequency_response.py
View File

@ -14,7 +14,7 @@ R = RF_Coil.resistance(22.5)
L = RF_Coil.induct_perry()
L = 0.26e-6
#L = 0.26e-6
f = R/(2*np.pi* L)

47
Stern_gerlach_separation/01_Calculate_trap_frequency.py
View File

@ -57,3 +57,50 @@ print(f"Fall without magnetic field = {-0.5 * 9.81 * dt**2 *1e3} mm")
print(f"s7 = {s7 * 1e3}mm")
print(f"diff = {(s8 - s7) * 1e3} mm")
# %%
# For thermal cloud size from paper
dt = 15e-3
sigma = 0.5e-3
T = 3.2e-6
sigma = np.sqrt(cs.k_B*T/cs.m_Dy_164) * dt
dz = 2 * sigma
print(f"cloud size: {sigma*1e3:.3f} mm")
#dz = 250e-6
Grad_Bz = 2 * dz * cs.m_Dy_164/(dt**2 * 1.241 * cs.mu_B)
print(" ")
print("For Stern-Gerlach separation:")
print(f"dBz/dz = {Grad_Bz*1e4*1e-2:.4f} G/cm")
print(" ")
# %%
dt1 = 10e-3
dt2 = 10e-3
dt = dt1 + dt2
sigma = np.sqrt(0.3e-3**2 + (np.sqrt(cs.k_B*T/cs.m_Dy_164) * dt)**2)
T = 3.2e-6
#sigma = np.sqrt(cs.k_B*T/cs.m_Dy_164) * dt
def cloud_sep(Grad_Bz, dt1 = 10e-3, dt2 = 0):
Grad_Bz*=1e-2
a = cs.mu_B * 1.241*Grad_Bz/cs.m_Dy_164
s1 = 0.5 * a *dt1**2
v1 = a * dt1
s2 = s1 + v1 * dt2
s3 = 0.5 * 9.81 *(dt1+dt2)**2
#print(f"fall:{s3*1e3}mm")
return s2*1e3
print(f"2 x sigma = {2*sigma*1e3}")
print(cloud_sep(8, dt1 = dt1, dt2 = dt2))
print(cloud_sep(30, dt1 = dt1, dt2 = dt2))

29
Stern_gerlach_separation/03_Displacement_2DMOT.py Normal file
View File

@ -0,0 +1,29 @@
# -*- coding: utf-8 -*-
"""
Created on Fri Aug 27 15:14:48 2021
@author: Joschka
"""
from src import physical_constants as cs
import numpy as np
mu = 9.9* cs.mu_B/8
def cloud_pos(mF,Grad_Bz, dt = 10e-3):
a = mF*mu*Grad_Bz/cs.m_Dy_164 - 9.81
s = 0.5 * a *dt**2
return s
Grad = 40*1e-2
t = 0.2/100
print(t)
print(cloud_pos(-8,Grad,dt=t ))

241
Thesis_Plots/Coil_Design/Fig_Coil_configuration.py Normal file
View File

@ -0,0 +1,241 @@
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from matplotlib.patches import Ellipse
from matplotlib import patches
from src import coil_class as BC
# %%
my_colors = {'light_green': '#97e144',
'orange': '#FF914D',
'light_grey': '#545454',
'pastel_blue': '#1b64d1',
'light_blue': '#71C8F4',
'purple': '#7c588c'}
# %%
def myarrow(ax,posX,posY,direction=1):
h=0.3
w=0.35 * direction
ax.arrow(posX,posY,w,h,linewidth=3,head_width=0,head_length=0)
ax.arrow(posX,posY,w,-h,linewidth=3,head_width=0,head_length=0)
# %%
AHH = BC.BCoil(HH=-1, distance=3*np.sqrt(3), radius=3, layers=1,windings=1,wire_width=0.1
,wire_height=0.1,is_round=True)
ylim = 4
xlim = 3
x = np.arange(-xlim, xlim+0.1, 0.5)
z = np.arange(-ylim, ylim+0.1, 0.5)
x_m_AHH, z_m_AHH = np.meshgrid(x, z)
I_current = 1
B = AHH.B_multiple_3d(I_current, x, z, raster=1)
B_tot = BC.BCoil.B_tot_3d(B)
for xx in range(0, len(x)):
for zz in range(0, len(z)):
if B_tot[zz, xx] > 1.5:
B[zz, xx, :] /= B_tot[zz, xx] / 1.5
# %%
HH = BC.BCoil(HH=1, distance=3, radius=3, layers=1,windings=1,wire_width=0.1
,wire_height=0.1,is_round=True)
ylim = 3
xlim = 3
x = np.arange(-xlim, xlim+0.1, 0.5)
z = np.arange(-ylim, ylim+0.1, 0.5)
x_m, z_m = np.meshgrid(x, z)
I_current = 1
B_HH = HH.B_multiple_3d(I_current, x, z, raster=1)
B_tot_HH = BC.BCoil.B_tot_3d(B_HH)
for xx in range(0, len(x)):
for zz in range(0, len(z)):
if B_tot_HH[zz, xx] > 1.5:
B_HH[zz, xx, :] /= B_tot_HH[zz, xx] / 1.5
# %%
mpl.rcParams.update({'font.size':11})
fig = plt.figure(figsize=(5,2.5),dpi=600)
ax = fig.add_subplot(121,aspect=1)
ypos = 1.5
width = 6
height = 1.3
lwidth = 4
#axes
ylim = 4
xlim = 3.5
e1 = Ellipse((0,ypos), width, height, linewidth = lwidth, fill = False, capstyle='projecting', zorder=2)
e2 = Ellipse((0,-ypos), width, height, linewidth = lwidth, fill = False, capstyle='projecting', zorder=2)
ax.add_patch(e1)
ax.add_patch(e2)
#add arrows to ellipse
myarrow(ax,0.5,0.86)
myarrow(ax,0.5,-2.14, direction=1)
# make axis
x_min = 3.3
ax.arrow(-x_min,0,2*x_min,0,head_width=0.2, length_includes_head=True, facecolor='black',zorder=0)
ax.text(x_min + 0.16 ,-0.15,'y',alpha=1)
y_min = 3.2
ax.arrow(0,-y_min,0, 2*y_min,head_width=0.2, length_includes_head=True, facecolor='black', zorder=0)
ax.text(-0.2, y_min +0.25, 'z')
ax.set_ylim(-ylim,ylim)
ax.set_xlim(-xlim,xlim)
# plot field
x_m, z_m = np.meshgrid(x, z)
arr_color = my_colors['orange']
zord=3
lim1=0
lim2=3
ax.quiver(x_m[lim1:lim2, :], z_m[lim1:lim2, :], B_HH[lim1:lim2, :, 0], B_HH[lim1:lim2, :, 1], color=arr_color,zorder=1)
lim1=3
lim2=6
ax.quiver(x_m[lim1:lim2, :], z_m[lim1:lim2, :], B_HH[lim1:lim2, :, 0], B_HH[lim1:lim2, :, 1], color=arr_color,zorder=3)
lim1=6
lim2=9
ax.quiver(x_m[lim1:lim2, :], z_m[lim1:lim2, :], B_HH[lim1:lim2, :, 0], B_HH[lim1:lim2, :, 1], color=arr_color,zorder=1)
lim1=9
lim2=12
ax.quiver(x_m[lim1:lim2, :], z_m[lim1:lim2, :], B_HH[lim1:lim2, :, 0], B_HH[lim1:lim2, :, 1], color=arr_color,zorder=3)
# add distance + radius annotation
d = patches.FancyArrowPatch((-width/2 - 0.5, -ypos-0.2), (-width/2 - 0.5, ypos+0.2), arrowstyle='<->',mutation_scale=10,
linewidth=1, color=my_colors['pastel_blue'])
ax.text(-width/2-1, -0.18, 'd', color=my_colors['pastel_blue'])
y_off = -3.4
R = patches.FancyArrowPatch((-width/2 , y_off), (0.13, y_off), arrowstyle='<->',mutation_scale=10,
linewidth=1, color=my_colors['pastel_blue'])
ax.text(-width/4-0.1,y_off-0.5, 'r', color=my_colors['pastel_blue'])
ax.add_patch(d)
ax.add_patch(R)
# general plot settings
ax.set_ylim(-ylim-0.5,ylim+0.5)
ax.set_xlim(-xlim-0.5,xlim+0.5)
ax.axis('off')
#ax.title.set_text('(a)',pos='left')
ax.set_title('(a)',loc='left')
# plt.savefig("Out/magnetic_configuration.png")
##########################################################################################################
ax = fig.add_subplot(122, aspect=1)
width = 6
scale = width/4 *np.sqrt(3) - ypos
ypos = ypos + scale
height = 1.3
print(ypos)
#axes
e1 = Ellipse((0,ypos), width, height, linewidth = lwidth, fill = False, capstyle='projecting', zorder=2)
e2 = Ellipse((0,-ypos), width, height, linewidth = lwidth, fill = False, capstyle='projecting', zorder=2)
ax.add_patch(e1)
ax.add_patch(e2)
#add arrows to ellipse
myarrow(ax,0.5,0.86+scale)
myarrow(ax,-0.5,-2.14-scale, direction=-1)
# make axis
x_min = 3.2
ax.arrow(-x_min,0,2*x_min,0,head_width=0.2, length_includes_head=True, facecolor='black',zorder=0)
ax.text(x_min + 0.16 ,-0.15,'y',alpha=1)
y_min = 4.1
ax.arrow(0,-y_min,0, 2*y_min,head_width=0.2, length_includes_head=True, facecolor='black', zorder=0)
ax.text(-0.2, y_min +0.25, 'z')
ax.set_ylim(-ylim,ylim)
ax.set_xlim(-xlim,xlim)
ax.axis('off')
#ax = fig.add_subplot(222,aspect=1)
x_m = x_m_AHH
z_m = z_m_AHH
print(np.shape(B))
arr_color = my_colors['orange']
zord=3
lim1=1
lim2=3
ax.quiver(x_m[lim1:lim2, :], z_m[lim1:lim2, :], B[lim1:lim2, :, 0], B[lim1:lim2, :, 1], color=arr_color,zorder=1)
ax.quiver(x_m[3, 1:-1], z_m[3, 1:-1], B[3, 1:-1, 0], B[3, 1:-1, 1], color=arr_color,zorder=1)
ax.quiver(x_m[3, 0], z_m[3, 0], B[3, 0, 0], B[3, 0, 1], color=arr_color,zorder=3)
ax.quiver(x_m[3, -1], z_m[3, -1], B[3, -1, 0], B[3, -1, 1], color=arr_color,zorder=3)
lim1=4
lim2=6
ax.quiver(x_m[lim1:lim2, :], z_m[lim1:lim2, :], B[lim1:lim2, :, 0], B[lim1:lim2, :, 1], color=arr_color,zorder=3)
lim1=6
lim2=13
ax.quiver(x_m[lim1:lim2, :], z_m[lim1:lim2, :], B[lim1:lim2, :, 0], B[lim1:lim2, :, 1], color=arr_color,zorder=1)
lim1=14
lim2=len(x_m)-1
ax.quiver(x_m[lim1:lim2, :], z_m[lim1:lim2, :], B[lim1:lim2, :, 0], B[lim1:lim2, :, 1], color=arr_color,zorder=3)
ax.quiver(x_m[13, 1:-1], z_m[13, 1:-1], B[13, 1:-1, 0], B[13, 1:-1, 1], color=arr_color,zorder=3)
ax.quiver(x_m[13, 0], z_m[13, 0], B[13, 0, 0], B[13, 0, 1], color=arr_color,zorder=1)
ax.quiver(x_m[13, -1], z_m[13, -1], B[13, -1, 0], B[13, -1, 1], color=arr_color,zorder=1)
# Add d + R annotations
d = patches.FancyArrowPatch((-width/2 - 0.5, -ypos-0.2), (-width/2 - 0.5, ypos+0.2), arrowstyle='<->',mutation_scale=10,
linewidth=1, color=my_colors['pastel_blue'])
ax.text(-width/2-1, -0.18, 'd', color=my_colors['pastel_blue'])
y_off = -4.2
R = patches.FancyArrowPatch((-width/2 , y_off), (0.13, y_off), arrowstyle='<->',mutation_scale=10,
linewidth=1, color=my_colors['pastel_blue'])
ax.text(-width/4-0.1,y_off-0.5, 'r', color=my_colors['pastel_blue'])
ax.add_patch(d)
ax.add_patch(R)
ax.set_ylim(-ylim-0.5,ylim+0.5)
ax.set_xlim(-xlim-0.5,xlim+0.5)
ax.axis('off')
ax.set_title('(b)',loc='left')
plt.savefig("C:/Users/Joschka/Desktop/Magnetic_field_coils_project/Thesis_Plots/Coil_Design/Out/magnetic_configuration.png")
plt.show()

200
Thesis_Plots/Field plots.py Normal file
View File

@ -0,0 +1,200 @@
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from matplotlib.patches import Ellipse
from matplotlib import patches
from src import coil_class as BC
from matplotlib.ticker import AutoMinorLocator
# %%
my_colors = {'light_green': '#97e144',
'orange': '#FF914D',
'light_grey': '#545454',
'pastel_blue': '#1b64d1',
'light_blue': '#71C8F4',
'purple': '#7c588c'}
mpl.rcParams.update({'font.size': 10, 'axes.linewidth': 1, 'lines.linewidth': 2, 'text.usetex': False, 'font.family': 'arial'})
mpl.rcParams['xtick.direction'] = 'in'
mpl.rcParams['ytick.direction'] = 'in'
mpl.rcParams['xtick.top'] = True
mpl.rcParams['ytick.right'] = True
def a_scat(B, a_bg = 1, B_0 = 0, DeltaB = 1):
return a_bg * (1 - (DeltaB/(B - B_0)))
mm = 1/25.4
# %%
HH_Coil = BC.BCoil(HH = 1, distance = 54, radius = 48, layers = 8, windings = 8, wire_height = 0.5,
wire_width = 0.5, insulation_thickness = (0.546-0.5)/2, is_round = True,
winding_scheme= 2)
HH_Coil.set_R_inner(45.6)
HH_Coil.set_d_min(47.4)
HH_Coil.print_info()
AHH_Coil = BC.BCoil(HH = -1, distance = 54, radius = 48, layers = HH_Coil.get_layers, windings=2 * HH_Coil.get_windings,
wire_height = 0.5, wire_width=0.5, insulation_thickness=(0.546-0.5)/2,
is_round = True, winding_scheme= 2, red_N=1)
AHH_Coil.set_R_inner(45.6)
AHH_Coil.set_d_max(77.6)
AHH_Coil.print_info()
# %%
x = np.linspace(-80, 80, 1000)
z = np.linspace(-80, 80, 1000)
I = 1
Bz, Bx = HH_Coil.B_tot_along_axis(I, x, z, raster = 7)
Bz_grad = np.abs(BC.BCoil.grad(Bz, z))
Bx_grad = np.abs(BC.BCoil.grad(Bx, x))
# %%
#HH plots
fig, axes = plt.subplots(2, 2, figsize=(5.8, 4), gridspec_kw={'height_ratios':[12,5]},dpi=600)
#fig.tight_layout()
lim1 = 55
lim1 =(-lim1, lim1)
lim2= 5
lim2 = (-lim2, lim2)
axs = axes[0, 0]
axs.plot(x, Bx, color=my_colors['pastel_blue'], label=r'$x/y-$axis')
axs.plot(z, Bz, color=my_colors['orange'], label=r'$z-$axis')
axs.set_ylabel(r'$B_{tot}$ (G/cm)')
axs.set_xlim(lim1)
axs.set_xticks(np.arange(-40,42,20))
axs.legend(loc=8)
axs = axes[1, 0]
axs.plot(x, Bx, color=my_colors['pastel_blue'])
axs.plot(z, Bz, color=my_colors['orange'])
axs.set_ylabel(r'$B_{tot}$ (G)')
axs.set_xlim(lim2)
axs.set_ylim(10.725,10.775)
axs.set_xlabel('distance to center (mm)')
axs.set_yticks([10.73,10.75,10.77])
axs.set_xticks(np.arange(-4,5,2))
axs = axes[0, 1]
axs.plot(x, Bx_grad, color=my_colors['pastel_blue'])
axs.plot(z, Bz_grad, color=my_colors['orange'])
axs.set_xlim(lim1)
axs.set_ylabel(r'$| ∇ B_{tot}|$ (G/cm)')
axs.set_xticks(np.arange(-40,42,20))
axs = axes[1, 1]
axs.plot(x, Bx_grad, color=my_colors['pastel_blue'])
axs.plot(z, Bz_grad, color=my_colors['orange'])
axs.set_ylabel(r'$| ∇ B_{tot}|$ (G/cm)')
axs.set_xlim(lim2)
axs.set_ylim(-0.01,0.09)
axs.set_yticks([0.0,0.04, 0.08])
axs.set_xticks(np.arange(-4,5,2))
axs.set_xlabel('distance to center (mm)')
for i in [0,1]:
for j in [0,1]:
axes[1,j].xaxis.set_minor_locator(AutoMinorLocator(2))
axes[0, j].xaxis.set_minor_locator(AutoMinorLocator(2))
axes[i,j].yaxis.set_minor_locator(AutoMinorLocator(2))
axes[i, 0].yaxis.set_label_coords(-0.19, 0.5)
axes[i, 1].yaxis.set_label_coords(-0.17, 0.5)
fig.tight_layout()
fig.savefig('Thesis_Plots/Coil_Design/Out/HH_field.png')
plt.show()
# %%
x = np.linspace(-80, 80, 1000)
z = np.linspace(-80, 80, 1000)
I = 1
Bz, Bx = AHH_Coil.B_field(I, x, z, raster = 7)
Bz_grad = np.abs(BC.BCoil.grad(Bz, z))
Bx_grad = np.abs(BC.BCoil.grad(Bx, x))
Bz = np.abs(Bz)
Bx = np.abs(Bx)
# %%
#AHH plots
fig, axes = plt.subplots(2, 2, figsize=(5.8, 4), gridspec_kw={'height_ratios':[12,5]},dpi=600)
#fig.tight_layout()
lim1 = 55
lim1 =(-lim1, lim1)
lim2= 5
lim2 = (-lim2, lim2)
axs = axes[0, 0]
axs.plot(x, Bx, color=my_colors['pastel_blue'], label=r'$x/y-$axis')
axs.plot(z, Bz, color=my_colors['orange'], label=r'$z-$axis')
axs.set_ylabel(r'$B_{tot}$ (G/cm)')
axs.set_xlim(lim1)
axs.set_xticks(np.arange(-40,42,20))
axs = axes[1, 0]
axs.plot(x, Bx, color=my_colors['pastel_blue'])
axs.plot(z, Bz, color=my_colors['orange'])
axs.set_ylabel(r'$B_{tot}$ (G)')
axs.set_xlim(lim2)
axs.set_ylim(-0.2, 3.5)
axs.set_xlabel('distance to center (mm)')
axs.set_yticks([0,1,2,3])
axs.set_xticks(np.arange(-4,5,2))
axs = axes[0, 1]
axs.plot(x, Bx_grad, color=my_colors['pastel_blue'], label=r'$x/y-$axis')
axs.plot(z, Bz_grad, color=my_colors['orange'], label=r'$z-$axis')
axs.set_xlim(lim1)
axs.set_ylabel(r'$| ∇ B_{tot}|$ (G/cm)')
axs.set_xticks(np.arange(-40,42,20))
axs.legend()
axs = axes[1, 1]
axs.plot(x, Bx_grad, color=my_colors['pastel_blue'])
axs.plot(z, Bz_grad, color=my_colors['orange'])
axs.set_ylabel(r'$| ∇ B_{tot}|$ (G/cm)')
axs.set_xlim(lim2)
axs.set_ylim(1.5,6.4)
#axs.set_yticks([0.0,0.04, 0.08])
axs.set_xticks(np.arange(-4,5,2))
axs.set_xlabel('distance to center (mm)')
for i in [0,1]:
for j in [0,1]:
axes[1,j].xaxis.set_minor_locator(AutoMinorLocator(2))
axes[0, j].xaxis.set_minor_locator(AutoMinorLocator(2))
axes[i,j].yaxis.set_minor_locator(AutoMinorLocator(2))
axes[i,0].yaxis.set_label_coords(-0.15,0.5)
axes[i, 1].yaxis.set_label_coords(-0.08, 0.5)
fig.tight_layout()
fig.savefig('Thesis_Plots/Coil_Design/Out/AHH_field.png')
plt.show()

65
Thesis_Plots/Final_design_figs.py Normal file
View File

@ -0,0 +1,65 @@
import matplotlib.pyplot as plt
import numpy as np
import matplotlib as mpl
# matplotlib.use('Qt5Agg')
from src import coil_class as BC
# %%
HH_Coil = BC.BCoil(HH=1, distance=54, radius=48, layers=8, windings=8, wire_height=0.5,
wire_width=0.5, insulation_thickness=(0.546 - 0.5) / 2, is_round=True,
winding_scheme=2)
HH_Coil.set_R_inner(45.6)
HH_Coil.set_d_min(47.4)
HH_Coil.print_info()
AHH_Coil = BC.BCoil(HH=-1, distance=54, radius=48, layers=HH_Coil.get_layers, windings=2 * HH_Coil.get_windings,
wire_height=0.5, wire_width=0.5, insulation_thickness=(0.546 - 0.5) / 2,
is_round=True, winding_scheme=2)
AHH_Coil.set_R_inner(45.6)
AHH_Coil.set_d_max(77.6)
AHH_Coil.print_info()
# %%
mpl.rcParams['xtick.direction'] = 'in'
mpl.rcParams['ytick.direction'] = 'in'
mpl.rcParams['xtick.top'] = True
mpl.rcParams['ytick.right'] = True
# mpl.rcParams['xtick.major.size'] = 10
# mpl.rcParams['xtick.major.width'] = 3
# mpl.rcParams['xtick.minor.size'] = 10
# mpl.rcParams['xtick.minor.width'] = 3
# mpl.rcParams['ytick.major.size'] = 10
# mpl.rcParams['ytick.major.width'] = 3
# mpl.rcParams['ytick.minor.size'] = 10
# mpl.rcParams['ytick.minor.width'] = 3
mpl.rcParams['font.size'] = 11
# mpl.rcParams.update({'font.size': 22, 'axes.linewidth': 3, 'lines.linewidth': 3})
# %%
#HH_Coil.plot_raster()
raster_value = 400
full_structure = HH_Coil.full_raster(raster_value, plot_radius_offset=0.00007) * 1e3
fig, ax = plt.subplots(figsize=(3.8, 3.8),dpi=600)
ax.scatter(full_structure[:, :, 1], full_structure[:, :, 0], linewidths=0.001)
ax.set_xlabel("radius ϕ [mm]")
ax.set_ylabel("z-axis [mm]")
ax.axvline(x=HH_Coil.get_R_inner() * 1e3 - 0.1, lw=5, color="red")
ext = 0.2
ax.set_xlim(1e3 * HH_Coil.get_R_inner() - ext- 0.2, 1e3 * (HH_Coil.get_R_inner()+HH_Coil.get_coil_height() )+ ext)
ax.set_ylim(1e3 * HH_Coil.get_zmin() - ext, 1e3 * HH_Coil.get_zmax() + ext)
plt.savefig('C:/Users/Joschka/Desktop/Magnetic_field_coils_project/Thesis_Plots/Coil_Design/Out/winding_scheme_2.png')
plt.show()

61
Thesis_Plots/New_iterative_final_coil.py Normal file
View File

@ -0,0 +1,61 @@
# -*- coding: utf-8 -*-
"""
Created on Tue Sep 7 13:18:18 2021
@author: Joschka
"""
import matplotlib.pyplot as plt
import numpy as np
import matplotlib
# matplotlib.use('Qt5Agg')
from src import coil_class as BC
# %%
scale = 1000
lim = 1
nr_points = (2 * lim) * scale + 1
x = np.linspace(-lim, lim, nr_points)
z = np.linspace(-lim, lim, nr_points)
# Wire_1 = [0.45, 0.514]
# Wire_1 = [0.475, 0.543]
# Wire_1 = [0.5, 0.568]
Wires = [[0.45, 0.514], [0.475, 0.543], [0.5, 0.568]]
for i in [0]:
Wire_1 = Wires[i]
print(f"Wire core = {Wire_1[0]} mm:")
print(" ")
for ll in [6, 8, 10]:
for ww in [7, 8, 9]:
print(f"layers = {ll}, windings = {ww}")
Coil_1 = BC.BCoil(HH=1, distance=54, radius=48, layers=ll, windings=ww, wire_height=Wire_1[0],
wire_width=Wire_1[0], insulation_thickness=(Wire_1[1] - Wire_1[0]) / 2, is_round=True,
winding_scheme=2)
# set radius plus distance
Coil_1.set_R_outer(49.5)
Coil_1.set_d_min(47.4)
# Coil_1.print_info()
print(f" Coil crossection area = {Coil_1.get_coil_width() * 1e3 * Coil_1.get_coil_height() * 1e3} mm2")
if (Coil_1.get_coil_width() * 1e3 * Coil_1.get_coil_height() * 1e3) < 16:
print("")
continue
#print(f" H = {Coil_1.get_coil_height() * 1e3}, W = {Coil_1.get_coil_width() * 1e3}")
I = 55 / Coil_1.get_N()
# I = 1
#print(f" Current needed for 13.8 G: I = {I} A")
# Coil_1.cooling(I, 22.5)
print(f" power = {Coil_1.power(I,22.5)}")
Coil_1.B_quality()
#print(f"B(0) = {Coil_1.max_field(I):.4f}")
print(" ")

72
Thesis_Plots/feshbach_res.py Normal file
View File

@ -0,0 +1,72 @@
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from matplotlib.patches import Ellipse
from matplotlib import patches
from src import coil_class as BC
# %%
my_colors = {'light_green': '#97e144',
'orange': '#FF914D',
'light_grey': '#545454',
'pastel_blue': '#1b64d1',
'light_blue': '#71C8F4',
'purple': '#7c588c'}
mpl.rcParams.update({'font.size': 10, 'axes.linewidth': 1, 'lines.linewidth': 2, 'text.usetex': False, 'font.family': 'arial'})
mpl.rcParams['xtick.direction'] = 'in'
mpl.rcParams['ytick.direction'] = 'in'
mpl.rcParams['xtick.top'] = True
mpl.rcParams['ytick.right'] = True
def a_scat(B, a_bg = 1, B_0 = 0, DeltaB = 1):
return a_bg * (1 - (DeltaB/(B - B_0)))
mm = 1/25.4
#%%
xlim = 3.9
ylim = 3.9
x_neg = np.linspace(-7, -0.0001, 500)
x_pos = np.linspace(0.0001, 5, 500)
color_scat = '#FF914D'
color_add = '#71C8F4'
color_grey = '#545454'
fig, ax = plt.subplots(1, 1, figsize=(65 * mm, 55 * mm),dpi=600)
plt.subplots_adjust(bottom=0.17, left=0.17)
ax.grid(alpha= 0.4)
ax.hlines(132/80, -xlim, xlim, color=color_add, linestyles=(2, (2, 2)), zorder=2)
ax.plot(x_neg, a_scat(x_neg), color=color_scat, zorder=1)
ax.plot(x_pos, a_scat(x_pos), color=color_scat, zorder=1)
# plt.ylim(-ylim, ylim)
ax.vlines(0, -ylim, ylim, color='grey', linestyles='dotted', zorder=0)
ax.hlines(0, -xlim, xlim, color='grey', linestyles='dotted', zorder=0)
ax.text(2.6, 132/80 + 0.3, r'a$_{\rm dd}$', color=color_add)
ax.arrow(0, 0, 1, 0, length_includes_head=True, width=0.15, head_length = 0.3, color=color_grey, zorder=2)
ax.text(0.14, 0.36, r'$\Delta_{\rm FR}$', color=color_grey, fontweight='bold')
mpl.patches.Rectangle((-xlim,0), 5, 1.2)
# plt.xticks()
ax.set_xlim(-xlim, xlim)
ax.set_ylim(-ylim, ylim)
ax.set_xticks(np.arange(-3,4,2))
ax.set_yticks(np.arange(-3, 4, 2))
ax.set_ylabel(r'a$_{\rm S}$/a$_{\rm bg}$')
ax.set_xlabel(r'(B-B$_0$)/$\Delta_{\rm FR}$')
ax.yaxis.set_label_coords(-0.12,0.5)
ax.xaxis.set_label_coords(0.5,-0.12)
plt.savefig('scattering.svg')
#plt.tight_layout()
plt.show()

114
src/coil_class.py
View File

@ -29,7 +29,7 @@ log.basicConfig(level=log.WARNING, format='%(message)s')
class BCoil:
def __init__(self, HH, distance, radius, layers, windings, wire_width, wire_height, insulation_thickness=0.
, is_round=False, winding_scheme=False):
, is_round=False, winding_scheme=False, red_N=0):
"""
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
@ -42,6 +42,7 @@ class BCoil:
: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
:param red_N: account for lower number of windings, reduces N by red_N
"""
if not is_round:
if winding_scheme == 1 or winding_scheme == 2:
@ -61,6 +62,7 @@ class BCoil:
self.insulation_thickness = insulation_thickness * 1e-3
self.is_round = is_round
self.winding_scheme = winding_scheme
self.red_N = red_N
# Standard get/set functions
@property
@ -153,7 +155,7 @@ class BCoil:
if self.winding_scheme == 2:
N -= self.layers // 2
log.debug(f"N = {N}")
return N
return N-self.red_N
def get_wire_length(self):
"""
@ -201,6 +203,10 @@ class BCoil:
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
if self.red_N !=0:
if xx == self.layers-1 and zz >= self.windings - self.red_N - 1:
continue
z_pos = z_start + zz * self.get_tot_wire_height()
R_pos = R_start + xx * self.get_tot_wire_width()
@ -217,11 +223,12 @@ class BCoil:
return outer_raster
def inner_raster(self, raster_value):
def inner_raster(self, raster_value, plot_radius_offset=0.):
"""
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
plot_radius_offset: reduce radius of wire for plotting
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
@ -252,24 +259,25 @@ class BCoil:
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:
if abs > (self.wire_width / 2 - plot_radius_offset):
delete_list.append(i)
inner_raster = np.delete(inner_raster, delete_list, axis=0)
return inner_raster
def full_raster(self, raster_value):
def full_raster(self, raster_value, plot_radius_offset=0.):
"""
Args:
raster_value: rastering value
plot_radius_offset: reduce radius of wire for plotting
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)
inner_ras = self.inner_raster(raster_value, plot_radius_offset=plot_radius_offset)
full_ras = np.zeros((len(outer_ras), len(inner_ras), 2))
for i in range(0, len(full_ras)):
@ -277,8 +285,8 @@ class BCoil:
return full_ras
def plot_raster(self, raster_value=100):
full_structure = self.full_raster(raster_value) * 1e3
def plot_raster(self, raster_value=100, plot_radius_offset=0.):
full_structure = self.full_raster(raster_value, plot_radius_offset=plot_radius_offset) * 1e3
if self.get_coil_width() > self.get_coil_height():
extension = self.get_coil_width()
else:
@ -474,6 +482,74 @@ class BCoil:
# print(f"time = {end - start} s")
return B_z, B_x
def B_field_asym(self, I_current, x, z, raster = 7, rad_diff=0.1e-3):
"""
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-rad_diff,
-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-rad_diff, -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-rad_diff, -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 B_quality(self,z_pos=1, raster = 7):
z = np.array([0, z_pos])
Bz, Bx = self.B_tot_along_axis(1, z, z, raster=raster)
print(f" Diff B +/- {z_pos} mm, relative: {(Bz[1] - Bz[0]) / Bz[0]}")
def Grad_quality(self, z_pos=1, I=1, raster = 7):
"""Does not work!!!"""
z = np.arange(-1, z_pos + 1, 0.001)
x = np.array([0])
z_0 = 1000
z_pos_ind = z_0 + 1000 * z_pos
print(z[z_0])
print(z[z_pos_ind])
B_z_tot, B_x_tot = self.B_field(I, x, z, raster=raster)
print(B_z_tot[z_0])
Bz_grad = BCoil.grad(B_z_tot, z)
print(Bz_grad[z_0])
print(f"rel Diff Grad B +/- {z_pos} mm relative: {(Bz_grad[z_pos_ind] - Bz_grad[z_0]) / Bz_grad[z_0]}")
def max_field(self, I_current, raster = 7):
B_z_0 = 0
@ -540,7 +616,7 @@ class BCoil:
# return bx_grad[len(x)//2]
def B_tot_along_axis(self, I_current, x, z, raster=10):
def B_tot_along_axis(self, I_current, x, z, raster=7):
"""
return B_tot_z, B_tot_x
Returns B_tot_z along z-axis and B_tot_x along x-axis,
@ -806,7 +882,7 @@ class BCoil:
def induct_perry(self):
"""
Calculates inductance of one rectangular coil via empirical formular by perry
Calculates inductance of one rectangular coil via empirical formular by perry in μH
"""
L = 4 * np.pi * (self.windings * self.layers) ** 2 * (1e2 * self.radius) ** 2 / (
@ -815,6 +891,14 @@ class BCoil:
@staticmethod
def resistivity_copper(T):
"""
For annealed copper
Args:
T: at temperature T
Returns: resistivity in Ohm * meter
"""
R_ref = cs.rho_copper_20
T_ref = 20 # degree celsius
alpha = 0.00393
@ -843,11 +927,11 @@ class BCoil:
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)
AHH_Coil = BCoil(HH=-1, distance=54, radius=48, layers=8, windings=16,
wire_height=0.5, wire_width=0.5, insulation_thickness=(0.546 - 0.5) / 2,
is_round=True, winding_scheme=2, red_N=1)
AHH_Coil.B_grad_quick_plot(1)
AHH_Coil.Grad_quality()
# main()

2
src/physical_constants.py
View File

@ -26,7 +26,7 @@ f_626 = 478.839e12
lambda_626 = 626.082e-9
#specific constants
rho_copper_20 = 1.68e-8
rho_copper_20 = 1.72e-8
density_copper = 8.940e3
# water

20
time_response/03_FINAL_time_response.py
View File

@ -15,7 +15,7 @@ HH_Coil = BC.BCoil(HH = 1, distance = 54, radius = 48, layers = 8, windings = 8,
winding_scheme= 2)
HH_Coil.set_R_inner(45.6)
HH_Coil.set_d_min(2*24.075)
# HH_Coil.print_info()
# AHH_Coil.print_info()
I = 1
@ -26,9 +26,9 @@ HH_Coil.cooling(I, 22)
# AHH_Coil = BC.BCoil(-1, 82 , 47.3 , 4, 6, wire_width= 1, wire_height= 1.5 ,layers_spacing = 0.25, windings_spacing= 0.25)
def I_current(Coil, I_0, t):
L = Coil.induct_perry()
L = Coil.induct_perry()*2
R = Coil.resistance(22.5)
R = Coil.resistance(22.5)*2
print(f"L={L}")
print(f" R= {R}")
tau = L / R
@ -50,10 +50,10 @@ def U_t(t, coil, U_0, I_end, scaling):
return U
# t = np.linspace(0, 0.005, 1000)
# plt.plot(t* 1e3,U_t(t, HH_Coil, 15, I))
# plt.plot(t* 1e3,U_t(t, AHH_Coil, 15, I))
# plt.xlabel("t [ms]")
# plt.show()
# print(U_t(0, HH_Coil, 15, I))
# print(U_t(0, AHH_Coil, 15, I))
@ -66,7 +66,9 @@ def I_t_exp(time, coil, I_end, U_0, scaling):
return I
def I_t_cut(time, coil, I_end, U_0):
I = U_0 / coil.resistance(22) * (1 - np.exp(- time / coil.tau()))
R = 2 * coil.resistance(22)
print(f"resistance I_cut = {R} Ω")
I = U_0 / R * (1 - np.exp(- time / coil.tau()))
if I >= I_end:
I = I_end
return I
@ -88,7 +90,7 @@ for U_0 in [10, 15, 28, 58]:
ax.plot(t * 1e3, I_t_cut_vec(t, HH_Coil, I, U_0), label=f"U_0 = {U_0} V, U_end = 1.5 V")
# for scaling in np.arange(2,5,0.5):
# ax.plot(t * 1e3, I_t_exp(t, HH_Coil, I, 15, scaling), label=f"Exponential decay U")
# ax.plot(t * 1e3, I_t_exp(t, AHH_Coil, I, 15, scaling), label=f"Exponential decay U")
ax.set_xlabel("time [ms]")
ax.set_ylabel("current I [A]")
@ -101,11 +103,11 @@ ax.plot(t * 1e6, I_current(HH_Coil, I, t), label=f"U_0 = U_end = 1.5 V")
for U_0 in [10, 15, 28, 58]:
ax.plot(t * 1e6, I_t_cut_vec(t, HH_Coil, I, U_0), label=f"U_0 = {U_0} V, U_end = 1.5 V")
# for scaling in np.arange(2,5,0.5):
# ax.plot(t * 1e3, I_t_exp(t, HH_Coil, I, 15, scaling), label=f"Exponential decay U")
# ax.plot(t * 1e3, I_t_exp(t, AHH_Coil, I, 15, scaling), label=f"Exponential decay U")
ax.set_xlabel("time [μs]")
ax.set_ylabel("current I [A]")
ax.set_xlim(0,120)
ax.set_xlim(0,300)
ax.legend()
plt.savefig('time_response_estimation.png', dpi=400)

126
time_response/04_Intermediate_time_response.py Normal file
View File

@ -0,0 +1,126 @@
# -*- coding: utf-8 -*-
"""
Created on Tue Sep 7 13:18:18 2021
@author: Joschka
"""
import matplotlib.pyplot as plt
import numpy as np
from src import coil_class as BC
HH_Coil = BC.BCoil(HH=1, distance=79.968, radius=80.228, layers=8, windings=8, wire_height=0.5,
wire_width=0.5, insulation_thickness=0.046 / 2, is_round=True,
winding_scheme=2)
AHH_Coil = BC.BCoil(HH=-1, distance=100.336, radius=85.016, layers=10, windings=10, wire_height=1.18,
wire_width=1.18, insulation_thickness=0.06, is_round=True,
winding_scheme=2)
I = 1
Max_field = HH_Coil.max_field(I)
HH_Coil.cooling(I, 22)
# AHH_Coil = BC.BCoil(-1, 82 , 47.3 , 4, 6, wire_width= 1, wire_height= 1.5 ,layers_spacing = 0.25, windings_spacing= 0.25)
def I_current(Coil, I_0, t):
L = Coil.induct_perry()*2
R = Coil.resistance(22.5)*2
print(f"L={L}")
print(f" R= {R}")
tau = L / R
print(f" τ = {tau}")
I = I_0 * (1 - np.exp(-R / L * t))
return I
def I_current_exp(I_0, R, L, t):
print("")
print(L / R)
I = I_0 * (1 - np.exp(-R / L * t))
return I
#%%
def U_t(t, coil, U_0, I_end, scaling):
U_end = I_end * coil.resistance(22)
U = (U_0-U_end) * np.exp(-t/coil.tau() * scaling) + U_end
return U
# t = np.linspace(0, 0.005, 1000)
# plt.plot(t* 1e3,U_t(t, AHH_Coil, 15, I))
# plt.xlabel("t [ms]")
# plt.show()
# print(U_t(0, AHH_Coil, 15, I))
U_vec = np.vectorize(U_t)
def I_t_exp(time, coil, I_end, U_0, scaling):
I = U_vec(time, coil, U_0, I_end, scaling) / coil.resistance(22) * (1 - np.exp(- time/coil.tau()))
return I
def I_t_cut(time, coil, I_end, U_0):
R = 2 * coil.resistance(22)
print(f"resistance I_cut = {R} Ω")
I = U_0 / R * (1 - np.exp(- time / coil.tau()))
if I >= I_end:
I = I_end
return I
I_t_cut_vec = np.vectorize(I_t_cut)
# %%
t = np.linspace(0, 0.003, 10000)
fig, axs = plt.subplots(2, 1, figsize = (7,7.5))
fig.suptitle(f"Time response HH-coil: I_max = {I} A --> Max Field = {Max_field:.2f} G \n \n I(t) = U(t) / R * (1 - exp(- R/L * t))")
ax = axs[0]
# ax.title("test")
ax.plot(t * 1e3, I_current(HH_Coil, I, t), label=f"U_0 = U_end = 1.5 V")
for U_0 in [10, 15, 28, 58]:
ax.plot(t * 1e3, I_t_cut_vec(t, HH_Coil, I, U_0), label=f"U_0 = {U_0} V, U_end = 1.5 V")
# for scaling in np.arange(2,5,0.5):
# ax.plot(t * 1e3, I_t_exp(t, AHH_Coil, I, 15, scaling), label=f"Exponential decay U")
ax.set_xlabel("time [ms]")
ax.set_ylabel("current I [A]")
ax.legend()
ax = axs[1]
ax.plot(t * 1e6, I_current(HH_Coil, I, t), label=f"U_0 = U_end = 1.5 V")
for U_0 in [10, 15, 28, 58]:
ax.plot(t * 1e6, I_t_cut_vec(t, HH_Coil, I, U_0), label=f"U_0 = {U_0} V, U_end = 1.5 V")
# for scaling in np.arange(2,5,0.5):
# ax.plot(t * 1e3, I_t_exp(t, AHH_Coil, I, 15, scaling), label=f"Exponential decay U")
ax.set_xlabel("time [μs]")
ax.set_ylabel("current I [A]")
ax.set_xlim(0,300)
ax.legend()
plt.savefig('time_response_estimation.png', dpi=400)
plt.show()
# %%
print(HH_Coil.max_field(0.4))
print(f"L = {HH_Coil.induct_perry() * 1e6} μH")
print(f"R = {HH_Coil.resistance(22)}")
V_max = 15
I_set = 5e-3
L = HH_Coil.induct_perry()
f_drop = V_max/(2*np.pi*L*I_set)
print(f_drop)

232
time_response/Lenny_Data.py Normal file
View File

@ -0,0 +1,232 @@
import numpy as np
import matplotlib.pyplot as plt
import scipy.io.wavfile as wavfile
from scipy.signal import savgol_filter
import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np
from scipy.optimize import curve_fit
from src import coil_class as BC
plt.rcParams["font.size"] = 12
samplingrateinput, GradientTenHzinput = wavfile.read("Z:/Users/Lenny/Power Supply/TimeSeriesOP549/FinalMeasurements/SquareWave10HzGradientCoilsWithPIinput.wav")
GradientTenHzinput = GradientTenHzinput*0.4/(2**14) /1.2093 +0.00247
samplingratesensingNoPI, GradientTenHzsensingNoPI = wavfile.read("Z:/Users/Lenny/Power Supply/TimeSeriesOP549/FinalMeasurements/SquareWave10HzGradientCoilsNoPIsensing.wav")
GradientTenHzsensingNoPI = GradientTenHzsensingNoPI*0.4/(2**14)/1.2093 +0.00247
samplingratesensingWithPI, GradientTenHzsensingWithPI = wavfile.read("Z:/Users/Lenny/Power Supply/TimeSeriesOP549/FinalMeasurements/SquareWave10HzGradientCoilsWithPIsensing.wav")
GradientTenHzsensingWithPI = GradientTenHzsensingWithPI*0.4/(2**14)/1.2093 +0.00247
t1 = np.arange(-199875/samplingrateinput, (len(GradientTenHzinput)-199876)/samplingrateinput, 1/samplingrateinput)*1000
plt.plot(t1[115491::],GradientTenHzsensingWithPI[115491::],label="With PI-controller")
#plt.plot(t1[115491::],GradientTenHzsensingNoPI[100000:-15491:],label="Without PI-controller")
plt.plot(t1[115491::],GradientTenHzinput[115491::],label="Input signal")
#plt.plot(t1[115491::],savgol_filter((GradientTenHzinput[115491::]-GradientTenHzsensingWithPI[115491::])/np.mean(GradientTenHzinput[0:1000]),50,3))
#plt.plot(t1[115491::],savgol_filter((GradientTenHzinput[115491::]-GradientTenHzsensingNoPI[100000:-15491:])/np.mean(GradientTenHzinput[0:1000]),50,3))
#plt.plot(t1[115491::],GradientTenHzinput[115491::]/np.mean(GradientTenHzinput[0:1000]))
#plt.xlim(t1[115491],t1[-1])
#plt.ylim(0.3, 0.65)
#plt.xlim(-0.05,0.3)
#plt.xlim(-0.05, 0.3)
plt.grid("both")
plt.xlabel("Time [ms]")
plt.ylabel("Current [A]")
handles, labels = plt.gca().get_legend_handles_labels()
order = [1,0]
plt.legend([handles[idx] for idx in order],[labels[idx] for idx in order], loc="upper right")
plt.show()
sampligrateinput, OffsetTenHzinput = wavfile.read("Z:/Users/Lenny/Power Supply/TimeSeriesOP549/FinalMeasurements/SquareWave10HzOffsetCoilsWithPIinput.wav")
OffsetTenHzinput = OffsetTenHzinput*0.4/(2**14)/1.2093 +0.00247
sampligratesensingNoPI, OffsetTenHzsensingNoPI = wavfile.read("Z:/Users/Lenny/Power Supply/TimeSeriesOP549/FinalMeasurements/SquareWave10HzOffsetCoilsNoPIsensing.wav")
OffsetTenHzsensingNoPI = OffsetTenHzsensingNoPI*0.4/(2**14)/1.2093 +0.00247
sampligratesensingWithPI, OffsetTenHzsensingWithPI = wavfile.read("Z:/Users/Lenny/Power Supply/TimeSeriesOP549/FinalMeasurements/SquareWave10HzOffsetCoilsWithPIsensing.wav")
OffsetTenHzsensingWithPI = OffsetTenHzsensingWithPI*0.4/(2**14)/1.2093 +0.00247
t2 = np.arange(-199069/samplingrateinput, (len(OffsetTenHzinput)-199070-833)/samplingrateinput, 1/samplingrateinput)*1000
plt.plot(t2[100000:],OffsetTenHzsensingWithPI[100833:], label="With PI-controller")
plt.plot(t2[100000:],OffsetTenHzsensingNoPI[100000:-833], label = "Without PI-controller")
plt.plot(t2[100000:],OffsetTenHzinput[100833:], label="Input signal")
#plt.plot(t2[100000:-150000],OffsetTenHzinput[100833:-150000]/np.mean(OffsetTenHzinput[0:1000]))
#plt.plot(t2[100000:-150000],savgol_filter((OffsetTenHzinput[100833:-150000]-OffsetTenHzsensingWithPI[100833:-150000])/np.mean(OffsetTenHzinput[0:1000]),50,3))
#plt.plot(t2[100000:-150000],savgol_filter((OffsetTenHzinput[100833:-150000]-OffsetTenHzsensingNoPI[100000:-150833])/np.mean(OffsetTenHzinput[0:1000]),50,3))
# My data
t = np. linspace (0,0.003,1000)
U_0 = 30
plt.plot(t * 1e3, I_t_cut_vec(t, int_HH_Coil, I, U_0) - 0.4, label=f"U_0 = {U_0} V, U_end = 1.5 V")
plt.plot(t * 1e3, 0.8 * I_current(int_HH_Coil, I, t) - 0.4, label=f"U_0 = U_end = 1.5 V")
#plt.ylim(0.3, 0.65)
# plt.xlim(-0.05,0.3)
plt.xlim(-0.05, 4 )
plt.grid("both")
plt.xlabel("Time [ms]")
plt.ylabel("Current [A]")
handles, labels = plt.gca().get_legend_handles_labels()
order = [1,0]
plt.legend([handles[idx] for idx in order],[labels[idx] for idx in order], loc="upper right")
plt.show()
# %%
def I_current(Coil, I_0, t):
L = Coil.induct_perry()*2
R = Coil.resistance(22.5)*2
print(f"L={L}")
print(f" R= {R}")
tau = L / R
print(f" τ = {tau}")
I = I_0 * (1 - np.exp(-R / L * t))
return I
def I_t_cut(time, coil, I_end, U_0):
R = 2 * coil.resistance(22)
#print(f"resistance I_cut = {R} Ω")
I = U_0 / R * (1 - np.exp(- time / coil.tau()))
if I >= I_end:
I = I_end
return I
I_t_cut_vec = np.vectorize(I_t_cut)
# %%
my_colors = {'light_green': '#97e144',
'orange': '#FF914D',
'light_grey': '#545454',
'pastel_blue': '#1b64d1',
'light_blue': '#71C8F4',
'purple': '#7c588c',
'pastel_red': '#D42024'}
#%%
mpl.rcParams['xtick.direction'] = 'in'
mpl.rcParams['ytick.direction'] = 'in'
mpl.rcParams['xtick.top'] = True
mpl.rcParams['ytick.right'] = True
mpl.rcParams['xtick.major.size'] = 10
mpl.rcParams['xtick.major.width'] = 3
mpl.rcParams['xtick.minor.size'] = 10
mpl.rcParams['xtick.minor.width'] = 3
mpl.rcParams['ytick.major.size'] = 10
mpl.rcParams['ytick.major.width'] = 3
mpl.rcParams['ytick.minor.size'] = 10
mpl.rcParams['ytick.minor.width'] = 3
mpl.rcParams.update({'font.size': 22, 'axes.linewidth': 3, 'lines.linewidth': 3})
# %%
HH_Coil = BC.BCoil(HH = 1, distance = 54, radius = 48, layers = 8, windings = 8, wire_height = 0.5,
wire_width = 0.5, insulation_thickness = (0.546-0.5)/2, is_round = True,
winding_scheme= 2)
HH_Coil.set_R_inner(45.6)
HH_Coil.set_d_min(2*24.075)
int_HH_Coil = BC.BCoil(HH=1, distance=79.968, radius=80.228, layers=8, windings=8, wire_height=0.5,
wire_width=0.5, insulation_thickness=0.046 / 2, is_round=True,
winding_scheme=2)
fig, ax1 = plt.subplots(figsize=(11, 8))
ylim = (0, 11.5)
t = np.linspace(0, 0.002, 10000)
i_to_B_inter = int_HH_Coil.max_field(1)
fig, ax = plt.subplots(figsize = (11,7))
#fig.suptitle(f"Time response HH-coil: I_max = {I} A --> Max Field = {Max_field:.2f} G \n \n I(t) = U(t) / R * (1 - exp(- R/L * t))")
ax.set_title("Time Response Intermediate Coils", y = 1.05)
# ax.text(0.6, 5, r'$I(t) = \frac{U(t)}{R} - \frac{L}{R} \cdot \frac{dI(t)}{dt} $', fontsize=34)
I = 0.8
ax.plot(t * 1e3, i_to_B_inter * (I_current(int_HH_Coil, I, t) - I / 2), label=f"Calculation without PI", zorder=1, color=my_colors['pastel_blue'])
U_0 = 28
ax.plot(t * 1e3, i_to_B_inter * (I_t_cut_vec(t, int_HH_Coil, I, U_0) - I / 2), label=f"Calculation with PI", zorder=1, color=my_colors['light_blue'])
i_to_B = HH_Coil.max_field(1)
I2 = 1
#ax.plot(t * 1e3, i_to_B * (I_t_cut_vec(t, HH_Coil, I2, U_0) - I2 / 2), label=f"Intermediate coils: calculation with PI", zorder=1, color=my_colors['light_blue'])
#plt.vlines(3.1e-2, 0, 10.64, zorder=2, linestyles=(0, (1.5, 3.06)),color=my_colors['orange'], label='t = 30 μs')
# for scaling in np.arange(2,5,0.5):
# ax.plot(t * 1e3, I_t_exp(t, AHH_Coil, I, 15, scaling), label=f"Exponential decay U")
#ax.plot(t*1e3, i_to_B_inter * I_current_fit(t*1e3, *popt))
#data
plt.plot(t2[100000:],i_to_B_inter * OffsetTenHzsensingWithPI[100833:], label="Measurement with PI", color=my_colors['pastel_red'])
plt.plot(t2[100000:],i_to_B_inter * OffsetTenHzsensingNoPI[100000:-833], label = "Measurement without PI", color = my_colors['orange'])
ax.set_xlabel("time [ms]")
ax.set_ylabel("Magnetic field [G]")
ax.set_ylim(-5,5)
ax.set_xlim(-0.09,1.8)
ax.legend()
plt.show()
# %% Fit
def I_current_fit(t, L1, L2, I_0, Coil=int_HH_Coil):
# L = Coil.induct_perry()*2
t_s = t * 1e-3
R = Coil.resistance(22.5)*2
#print(f"L={L}")
#print(f" R= {R}")
#tau = L / R
#print(f" τ = {tau}")
I = I_0 * ((1 - np.exp(-R / L1 * t_s))+ (1 - np.exp(-R / L2 * t_s)))
return I-I_0/2
# %%
def I_current_fit(t, L1, I_0, Coil=int_HH_Coil):
# L = Coil.induct_perry()*2
t_s = t * 1e-3
R = Coil.resistance(22.5)*2
#print(f"L={L}")
#print(f" R= {R}")
#tau = L / R
#print(f" τ = {tau}")
I = I_0 * (1 - np.exp(-R / L1 * t_s))
return I-I_0/2
# %%
p0 = [20e-3 , 0.8]
#bounds = ([-50, 0], [0, 5])
lim1 = 1.994e5
lim2 = 2.5e5
lim1 = int(lim1)
lim2 = int(lim2)
t_plot = np.copy(t2)
t_plot = t_plot[lim1:lim2]
I_plot = np.copy(OffsetTenHzsensingNoPI)
I_plot = I_plot[lim1:lim2]
print(t_plot)
print(I_plot)
popt, pcov = curve_fit(I_current_fit, t_plot, I_plot, p0=p0)
print(popt)
# %%
t = np.linspace(0,20, 10000)
#popt = [0.02,0.8]
fig, ax = plt.subplots(figsize = (11,7))
ax.plot(t, I_current_fit(t, *popt))
ax.plot(t2[100000:], OffsetTenHzsensingNoPI[100000:-833], label = "Intermediate coil: without PI", color = my_colors['orange'])
#ax.plot(t_plot, I_plot, label = "Intermediate coil: without PI", color = my_colors['orange'])
# ax.vlines(t2[lim1], -0.5, 0.5)
ax.vlines(t2[lim2], - 0.5, 0.5)
ax.set_xlabel("time [ms]")
ax.set_ylabel("Current [A]")
#ax.set_ylim(-0.5, 0.5)
ax.set_xlim(-0.09,5)
ax.legend()
plt.show()

BIN
time_response_estimation.png
View File

Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 401 KiB

After

Width:  |  Height:  |  Size: 407 KiB