diff --git a/Coil_geometry/00_Pre_FINAL_coil_geometry.py b/Coil_geometry/00_Pre_FINAL_coil_geometry.py index fcbed74..4c0e09b 100644 --- a/Coil_geometry/00_Pre_FINAL_coil_geometry.py +++ b/Coil_geometry/00_Pre_FINAL_coil_geometry.py @@ -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()) diff --git a/Coil_geometry/AHH/08_FINAL_AHH.py b/Coil_geometry/AHH/08_FINAL_AHH.py index ce17cf6..1bc305b 100644 --- a/Coil_geometry/AHH/08_FINAL_AHH.py +++ b/Coil_geometry/AHH/08_FINAL_AHH.py @@ -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)) diff --git a/Coil_geometry/Additional coils/00_Test_coil_around viewport.py b/Coil_geometry/Additional coils/00_Test_coil_around viewport.py index e8df969..dd52074 100644 --- a/Coil_geometry/Additional coils/00_Test_coil_around viewport.py +++ b/Coil_geometry/Additional coils/00_Test_coil_around viewport.py @@ -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() diff --git a/Coil_geometry/HH/04_Iterative_Testing.py b/Coil_geometry/HH/04_Iterative_Testing.py index e56aa65..731a45a 100644 --- a/Coil_geometry/HH/04_Iterative_Testing.py +++ b/Coil_geometry/HH/04_Iterative_Testing.py @@ -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]}") diff --git a/Coil_geometry/HH/05_try_diff_geometry_HH1.py b/Coil_geometry/HH/05_try_diff_geometry_HH1.py index 06e65ab..fe698ba 100644 --- a/Coil_geometry/HH/05_try_diff_geometry_HH1.py +++ b/Coil_geometry/HH/05_try_diff_geometry_HH1.py @@ -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") diff --git a/Coil_geometry/HH/07_02_testing B_tot.py b/Coil_geometry/HH/07_02_testing B_tot.py index 4b3c63f..6dd6c43 100644 --- a/Coil_geometry/HH/07_02_testing B_tot.py +++ b/Coil_geometry/HH/07_02_testing B_tot.py @@ -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) diff --git a/Coil_geometry/HH/11_Final_HH.py b/Coil_geometry/HH/11_Final_HH.py index d12707e..3c7ae7f 100644 --- a/Coil_geometry/HH/11_Final_HH.py +++ b/Coil_geometry/HH/11_Final_HH.py @@ -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") diff --git a/Coil_geometry/HH/12_Final_Plotting.py b/Coil_geometry/HH/12_Final_Plotting.py index 5671c75..ea10e54 100644 --- a/Coil_geometry/HH/12_Final_Plotting.py +++ b/Coil_geometry/HH/12_Final_Plotting.py @@ -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) diff --git a/Coil_geometry/HH/15_new_HH_coil_iterative_winding_decision.py b/Coil_geometry/HH/15_new_HH_coil_iterative_winding_decision.py index d93d1ed..04ac4bc 100644 --- a/Coil_geometry/HH/15_new_HH_coil_iterative_winding_decision.py +++ b/Coil_geometry/HH/15_new_HH_coil_iterative_winding_decision.py @@ -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() diff --git a/Coil_geometry/New intermediate pair of coils/01_Coil_calculations.py b/Coil_geometry/New intermediate pair of coils/01_Coil_calculations.py index 412a9b2..b92cf95 100644 --- a/Coil_geometry/New intermediate pair of coils/01_Coil_calculations.py +++ b/Coil_geometry/New intermediate pair of coils/01_Coil_calculations.py @@ -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() diff --git a/Coil_geometry_AHH/Optimum_distance.py b/Coil_geometry_AHH/Optimum_distance.py index c931d71..124f69e 100644 --- a/Coil_geometry_AHH/Optimum_distance.py +++ b/Coil_geometry_AHH/Optimum_distance.py @@ -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) diff --git a/Coil_geometry_AHH/Theo_Opt_Dist.py b/Coil_geometry_AHH/Theo_Opt_Dist.py new file mode 100644 index 0000000..2e852e6 --- /dev/null +++ b/Coil_geometry_AHH/Theo_Opt_Dist.py @@ -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() diff --git a/Cooling/05_POWER_estimation_water_cooling.py b/Cooling/05_POWER_estimation_water_cooling.py index e5ce273..d239ee0 100644 --- a/Cooling/05_POWER_estimation_water_cooling.py +++ b/Cooling/05_POWER_estimation_water_cooling.py @@ -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 diff --git a/FINAL_Coil/FINAL_COIL_configuration.py b/FINAL_Coil/FINAL_COIL_configuration.py index 5944f3e..5d72ff0 100644 --- a/FINAL_Coil/FINAL_COIL_configuration.py +++ b/FINAL_Coil/FINAL_COIL_configuration.py @@ -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)) + diff --git a/FINAL_Coil/Final_Field_quality.py b/FINAL_Coil/Final_Field_quality.py index feb39e7..d5e1dbc 100644 --- a/FINAL_Coil/Final_Field_quality.py +++ b/FINAL_Coil/Final_Field_quality.py @@ -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]}") diff --git a/FINAL_Coil/Final_Gradient_quality.py b/FINAL_Coil/Final_Gradient_quality.py new file mode 100644 index 0000000..60cd229 --- /dev/null +++ b/FINAL_Coil/Final_Gradient_quality.py @@ -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 ############################################################################ +############################################################################### +############################################################################### +""" diff --git a/FINAL_Coil/Test_different_diameters.py b/FINAL_Coil/Test_different_diameters.py new file mode 100644 index 0000000..3368b9b --- /dev/null +++ b/FINAL_Coil/Test_different_diameters.py @@ -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 ############################################################################ +############################################################################### +############################################################################### +""" diff --git a/Noise/00_Simple_testing.py b/Noise/00_Simple_testing.py index 01731e0..b8b08b7 100644 --- a/Noise/00_Simple_testing.py +++ b/Noise/00_Simple_testing.py @@ -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(" ") """ diff --git a/Plots_Erlangen/Plots.py b/Plots_Erlangen/Plots.py index ce9fc1d..8315a55 100644 --- a/Plots_Erlangen/Plots.py +++ b/Plots_Erlangen/Plots.py @@ -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]") diff --git a/Plots_Erlangen/Plots_midterm_presentation.py b/Plots_Erlangen/Plots_midterm_presentation.py index ed316fa..e31de90 100644 --- a/Plots_Erlangen/Plots_midterm_presentation.py +++ b/Plots_Erlangen/Plots_midterm_presentation.py @@ -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() \ No newline at end of file +Coil1.plot_raster() diff --git a/RF_coil/frequency_response.py b/RF_coil/frequency_response.py index aa4e2c3..53d75a4 100644 --- a/RF_coil/frequency_response.py +++ b/RF_coil/frequency_response.py @@ -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) diff --git a/Stern_gerlach_separation/01_Calculate_trap_frequency.py b/Stern_gerlach_separation/01_Calculate_trap_frequency.py index fb3d60c..dc7cbc4 100644 --- a/Stern_gerlach_separation/01_Calculate_trap_frequency.py +++ b/Stern_gerlach_separation/01_Calculate_trap_frequency.py @@ -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)) diff --git a/Stern_gerlach_separation/03_Displacement_2DMOT.py b/Stern_gerlach_separation/03_Displacement_2DMOT.py new file mode 100644 index 0000000..a6fa6ae --- /dev/null +++ b/Stern_gerlach_separation/03_Displacement_2DMOT.py @@ -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 )) + + + + + diff --git a/Thesis_Plots/Coil_Design/Fig_Coil_configuration.py b/Thesis_Plots/Coil_Design/Fig_Coil_configuration.py new file mode 100644 index 0000000..181049b --- /dev/null +++ b/Thesis_Plots/Coil_Design/Fig_Coil_configuration.py @@ -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() diff --git a/Thesis_Plots/Field plots.py b/Thesis_Plots/Field plots.py new file mode 100644 index 0000000..a14549d --- /dev/null +++ b/Thesis_Plots/Field plots.py @@ -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() \ No newline at end of file diff --git a/Thesis_Plots/Final_design_figs.py b/Thesis_Plots/Final_design_figs.py new file mode 100644 index 0000000..2f012a2 --- /dev/null +++ b/Thesis_Plots/Final_design_figs.py @@ -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() diff --git a/Thesis_Plots/New_iterative_final_coil.py b/Thesis_Plots/New_iterative_final_coil.py new file mode 100644 index 0000000..639c327 --- /dev/null +++ b/Thesis_Plots/New_iterative_final_coil.py @@ -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(" ") diff --git a/Thesis_Plots/feshbach_res.py b/Thesis_Plots/feshbach_res.py new file mode 100644 index 0000000..6fd4ed8 --- /dev/null +++ b/Thesis_Plots/feshbach_res.py @@ -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() \ No newline at end of file diff --git a/src/coil_class.py b/src/coil_class.py index 3af2e23..715cbe8 100644 --- a/src/coil_class.py +++ b/src/coil_class.py @@ -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() diff --git a/src/physical_constants.py b/src/physical_constants.py index cd26612..e05dd99 100644 --- a/src/physical_constants.py +++ b/src/physical_constants.py @@ -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 diff --git a/time_response/03_FINAL_time_response.py b/time_response/03_FINAL_time_response.py index 5d3160a..7686a2b 100644 --- a/time_response/03_FINAL_time_response.py +++ b/time_response/03_FINAL_time_response.py @@ -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) diff --git a/time_response/04_Intermediate_time_response.py b/time_response/04_Intermediate_time_response.py new file mode 100644 index 0000000..a1b0037 --- /dev/null +++ b/time_response/04_Intermediate_time_response.py @@ -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) diff --git a/time_response/Lenny_Data.py b/time_response/Lenny_Data.py new file mode 100644 index 0000000..09f7422 --- /dev/null +++ b/time_response/Lenny_Data.py @@ -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() diff --git a/time_response_estimation.png b/time_response_estimation.png index f318e6a..49625c0 100644 Binary files a/time_response_estimation.png and b/time_response_estimation.png differ