from dataclasses import dataclass import magpylib as magpy import matplotlib.pyplot as plt import numpy as np import random import json from concurrent.futures import ProcessPoolExecutor import concurrent.futures @dataclass class AqpBookerParameters: """ Contains the default parameters for QT3 atomic quantum processor magnet geometry Base units for magpylib are millimeters, milliteslas, degrees, and Amperes """ # assume x, y, z coils are all round, may consider racetrack coil later # assume x and y axis coils are symmetric with respect to the sample center, current has same magnitude, in the opposite direction # assume z axis coil are not symmetric # shape of the coils, can be 'round' or 'racetrack' # race track means the coil is a rectangle with rounded corners # x, y and z can be 0, which means the coil is round with r "radius" # # x/y/z # ___________ radius # ◜ ◝ # | | # | | # x/y/z | | # | | # | | # | | # ◟___________◞ # x_axis_coil_y_length: float = 0 x_axis_coil_z_length: float = 60 x_axis_coil_corner_radius: float = 18 x_axis_coil_winding_number: int = 200 x_axis_coil_current_a: float = -1.5 x_axis_coil_current_b: float = 1.5 x_axis_coil_center_a: tuple = (38, 0, 0) x_axis_coil_center_b: tuple = (-38, 0, 0) y_axis_coil_x_length: float = 0 y_axis_coil_z_length: float = 60 y_axis_coil_corner_radius: float = 18 y_axis_coil_winding_number: int = 200 y_axis_coil_current_a: float = -1.5 y_axis_coil_current_b: float = 1.5 y_axis_coil_center_a: tuple = (0, 38, 0) y_axis_coil_center_b: tuple = (0, -38, 0) # z axis coil not used for 2D MOT z_axis_coil_x_length: float = 27 z_axis_coil_y_length: float = 27 z_axis_coil_corner_radius: float = 10 z_axis_coil_winding_number: int = 1 z_axis_coil_current_a: float = -1*0 z_axis_coil_current_b: float = 0 z_axis_coil_center_a: tuple = (0, 0, 80) z_axis_coil_center_b: tuple = (0, 0, -80) sample_center: tuple = (0, 0, 0) x_wire_length: float = ((x_axis_coil_y_length + x_axis_coil_z_length)*2 + x_axis_coil_corner_radius*2*np.pi) * x_axis_coil_winding_number /1000 # meters y_wire_length: float = ((y_axis_coil_x_length + y_axis_coil_z_length)*2 + y_axis_coil_corner_radius*2*np.pi) * y_axis_coil_winding_number /1000 # meters z_wire_length: float = ((z_axis_coil_x_length + z_axis_coil_y_length)*2 + z_axis_coil_corner_radius*2*np.pi) * z_axis_coil_winding_number /1000 # meters wire_selection: int = 24 # AWG wire_diameter: float = 0.127 * (92 ** ((36 - wire_selection)/39)) # mm # R = rho * L / A x_wire_resistance: float = (1.7241*10**(-8)) * x_wire_length / (np.pi*((wire_diameter/2/1000) ** 2)) # ohm y_wire_resistance: float = (1.7241*10**(-8)) * y_wire_length / (np.pi*((wire_diameter/2/1000) ** 2)) # ohm z_wire_resistance: float = (1.7241*10**(-8)) * z_wire_length / (np.pi*((wire_diameter/2/1000) ** 2)) # ohm wire_weight: float = x_wire_length * np.pi * ((wire_diameter/2000) ** 2) * 8941 # winding_layer_number: int = 3 # x_axis_coil_thickness: float = x_axis_coil_winding_number * wire_diameter / winding_layer_number # y_axis_coil_thickness: float = y_axis_coil_winding_number * wire_diameter / winding_layer_number # z_axis_coil_thickness: float = z_axis_coil_winding_number * wire_diameter / winding_layer_number # using fixed thickness x_axis_coil_thickness: float = 20 y_axis_coil_thickness: float = 20 z_axis_coil_thickness: float = 20 position_tolerance: float = 1 # mm angle_tolerance: float = 1 # degree def get_default_aqp_parameters(): p = AqpBookerParameters() return p def build_magpy_collection(p: AqpBookerParameters = get_default_aqp_parameters()) -> magpy.Collection: # assume the coils are all rounds. x_axis_coll = magpy.Collection() y_axis_coll = magpy.Collection() z_axis_coll = magpy.Collection() # field_sensor = magpy.Sensor(position=p.sample_center, style_label='field_sensor') # coll = magpy.Collection(x_axis_coll, y_axis_coll, z_axis_coll, field_sensor) # coll = magpy.Collection(x_axis_coll, y_axis_coll, z_axis_coll) coll = magpy.Collection(x_axis_coll, y_axis_coll) # sample_number = 1000 # ts = np.linspace(-1 * p.x_axis_coil_winding_number/2, p.x_axis_coil_winding_number/2, sample_number) x_vertices = shape_racetrack(p.x_axis_coil_winding_number, p.x_axis_coil_corner_radius, p.x_axis_coil_y_length, p.x_axis_coil_z_length, p.x_axis_coil_thickness) x_axis_coil_a = magpy.current.Line(current=p.x_axis_coil_current_a, vertices=x_vertices, position=p.x_axis_coil_center_a, style_label='x_axis_coil_a', style_color = 'r') x_axis_coil_b = magpy.current.Line(current=p.x_axis_coil_current_b, vertices=x_vertices, position=p.x_axis_coil_center_b, style_label='x_axis_coil_b', style_color = 'g') x_axis_coil_a.rotate_from_angax(90, 'y') x_axis_coil_b.rotate_from_angax(90, 'y') x_axis_coil_a.rotate_from_angax(90, 'x') x_axis_coil_b.rotate_from_angax(90, 'x') x_axis_coll.add(x_axis_coil_a, x_axis_coil_b) # y_axis_coil_a = magpy.current.Loop(current = p.y_axis_coil_current_a, diameter=p.y_axis_coil_diameter_a,\ # position= p.y_axis_coil_center_a, style_label='y_axis_coil_a') # y_axis_coil_b = magpy.current.Loop(current = p.y_axis_coil_current_b, diameter=p.y_axis_coil_diameter_b,\ # position= p.y_axis_coil_center_b, style_label='y_axis_coil_b') y_vertices = shape_racetrack(p.y_axis_coil_winding_number, p.y_axis_coil_corner_radius, p.y_axis_coil_x_length, p.y_axis_coil_z_length, p.y_axis_coil_thickness) y_axis_coil_a = magpy.current.Line(current=p.y_axis_coil_current_a, vertices=y_vertices, position=p.y_axis_coil_center_a, style_label='y_axis_coil_a', style_color='b') y_axis_coil_b = magpy.current.Line(current=p.y_axis_coil_current_b, vertices=y_vertices, position=p.y_axis_coil_center_b, style_label='y_axis_coil_b', style_color='y') y_axis_coil_a.rotate_from_angax(90, 'x') y_axis_coil_b.rotate_from_angax(90, 'x') y_axis_coll.add(y_axis_coil_a, y_axis_coil_b) # z_axis_coil_a = magpy.current.Loop(current = p.z_axis_coil_current_a, diameter=p.z_axis_coil_diameter_a,\ # position= p.z_axis_coil_center_a, style_label='z_axis_coil_a') # z_axis_coil_b = magpy.current.Loop(current = p.z_axis_coil_current_b, diameter=p.z_axis_coil_diameter_b,\ # position= p.z_axis_coil_center_b, style_label='z_axis_coil_b') # z_vertices = shape_racetrack(p.z_axis_coil_winding_number, p.z_axis_coil_corner_radius, p.z_axis_coil_x_length, p.z_axis_coil_y_length, p.z_axis_coil_thickness) # z_axis_coil_a = magpy.current.Line(current=p.z_axis_coil_current_a, vertices=z_vertices, position=p.z_axis_coil_center_a, style_label='z_axis_coil_a', style_color='m') # z_axis_coil_b = magpy.current.Line(current=p.z_axis_coil_current_b, vertices=z_vertices, position=p.z_axis_coil_center_b, style_label='z_axis_coil_b', style_color='c') # z_axis_coil_a.rotate_from_angax(90, 'z') # z_axis_coil_b.rotate_from_angax(90, 'z') # z_axis_coll.add(z_axis_coil_a, z_axis_coil_b) return coll def simulate_mag_field(p: AqpBookerParameters = get_default_aqp_parameters()): coll = build_magpy_collection(p) # field_sensor = coll.sensors[0] print("x direction gradient is", get_gradient(coll, np.asarray(p.sample_center), 'x')) print("y direction gradient is", get_gradient(coll, np.asarray(p.sample_center), 'y')) print("z direction gradient is", get_gradient(coll, np.asarray(p.sample_center), 'z')) image = magpy.show(coll) def shape_racetrack(winding_number, corner_radius, edge_1_length, edge_2_length, coil_thickness) -> np.ndarray: sample_number = 10000 ts = np.linspace(-1 * winding_number/2, winding_number/2, sample_number) thick = np.linspace(-1 * coil_thickness/2, coil_thickness/2, sample_number) vertices = np.c_[corner_radius*np.cos(ts*2*np.pi), corner_radius*np.sin(ts*2*np.pi), thick] for i in range(len(vertices)): if (np.cos(ts[i] * 2 * np.pi) >= 0 and np.sin(ts[i] * 2 * np.pi) >= 0): vertices[i][0] += edge_1_length/2 vertices[i][1] += edge_2_length/2 elif (np.cos(ts[i] * 2 * np.pi) < 0 and np.sin(ts[i] * 2 * np.pi) >= 0): vertices[i][0] -= edge_1_length/2 vertices[i][1] += edge_2_length/2 elif (np.cos(ts[i] * 2 * np.pi) < 0 and np.sin(ts[i] * 2 * np.pi) < 0): vertices[i][0] -= edge_1_length/2 vertices[i][1] -= edge_2_length/2 elif (np.cos(ts[i] * 2 * np.pi) >= 0 and np.sin(ts[i] * 2 * np.pi) < 0): vertices[i][0] += edge_1_length/2 vertices[i][1] -= edge_2_length/2 vertices = np.append(vertices, [[vertices[0][0], vertices[0][1], vertices[-1][2]]], axis=0) # print(vertices) return vertices def get_gradient(coll: magpy.Collection, location: np.ndarray, direction="xyz"): diff = 0.001 # mm B0 = coll.getB(location) if direction == "xyz": B1 = coll.getB(location + np.array([diff, 0, 0])) x_gradient = (B1[0] - B0[0])/diff # mT/mm B1 = coll.getB(location + np.array([0, diff, 0])) y_gradient = (B1[1] - B0[1])/diff B1 = coll.getB(location + np.array([0, 0, diff])) z_gradient = (B1[2] - B0[2])/diff gradient = np.array([x_gradient/100, y_gradient/100, z_gradient/100]) else: if direction == 'x': B1 = coll.getB(location + np.array([diff, 0, 0])) gradient = (B1[0] - B0[0])/diff # mT/mm elif direction == 'y': B1 = coll.getB(location + np.array([0, diff, 0])) gradient = (B1[1] - B0[1])/diff elif direction == 'z': B1 = coll.getB(location + np.array([0, 0, diff])) gradient = (B1[2] - B0[2])/diff else: raise ValueError("direction must be one of 'x', 'y', 'z'") gradient = gradient * 100 # convert to G/cm if B0 is None or B1 is None: raise ValueError("B0 or B1 is None") return gradient def heat_generation(p: AqpBookerParameters = get_default_aqp_parameters()): # heat map for winding number and current winding_min = 50 winding_max = 300 current_min = 1 current_max = 3 sample_points = 50 winding_number = np.linspace(winding_min, winding_max, sample_points) current = np.linspace(current_min, current_max, sample_points) data = np.zeros((len(winding_number), len(current))) data_heat = np.zeros((len(winding_number), len(current))) for i in range(len(winding_number)): for j in range(len(current)): print("calculating: ", i, ", ", j) p.x_axis_coil_winding_number = winding_number[i] p.y_axis_coil_winding_number = winding_number[i] p.x_axis_coil_current_a = -1 * current[j] p.x_axis_coil_current_b = current[j] p.y_axis_coil_current_a = -1 * current[j] p.y_axis_coil_current_b = current[j] coll = build_magpy_collection(p) data[i][j] = get_gradient(coll, np.asarray(p.sample_center), 'x') if np.abs(data[i][j]) > 15: x_heat_production = p.x_wire_resistance * p.x_axis_coil_current_a**2 wire_weight = p.wire_weight delta_temp = x_heat_production * 1 / 0.385 / (wire_weight * 1000) data_heat[i][j] = delta_temp # export the data to csv np.savetxt("data.csv", data, delimiter=",") np.savetxt("data_heat.csv", data_heat, delimiter=",") plt.subplot(1, 2, 1) plt.imshow(data, cmap='hot', interpolation='nearest', extent=[current_min, current_max, winding_max, winding_min], aspect='auto') plt.colorbar() plt.subplot(1, 2, 2) plt.imshow(data_heat, cmap='hot', interpolation='nearest', extent=[current_min, current_max, winding_max, winding_min], aspect='auto') plt.colorbar() plt.show() def B_field_plot(p: AqpBookerParameters = get_default_aqp_parameters(), coll: magpy.Collection = None): if coll is None: coll = build_magpy_collection(p) fig, (ax1, ax2, ax3) = plt.subplots(3, 1) # plot along x axis x_axis_field = [] for i in np.linspace(-1 * p.z_axis_coil_x_length/2, 1*p.z_axis_coil_x_length/2, 100): B = coll.getB((i, p.x_axis_coil_center_a[1], p.x_axis_coil_center_a[2])) x_axis_field.append(np.abs(B[0])) ax1.plot(np.linspace(-1 * p.z_axis_coil_x_length/2, 1*p.z_axis_coil_x_length/2, 100), x_axis_field) ax1.set_xlabel('x (mm)') ax1.set_ylabel('B (mT)') ax1.set_title('B field along x axis') # plot along y axis y_axis_field = [] for i in np.linspace(-1 * p.z_axis_coil_y_length/2, 1* p.z_axis_coil_y_length/2, 100): B = coll.getB((p.y_axis_coil_center_a[0], i, p.y_axis_coil_center_a[2])) y_axis_field.append(np.abs(B[1])) ax2.plot(np.linspace(-1 * p.z_axis_coil_y_length/2, 1*p.z_axis_coil_y_length/2, 100), y_axis_field) ax2.set_xlabel('y (mm)') ax2.set_ylabel('B (mT)') ax2.set_title('B field along y axis') # plot along z axis z_axis_field = [] for i in np.linspace(-1 * p.x_axis_coil_z_length/2, 1*p.x_axis_coil_z_length/2, 100): B = coll.getB((p.z_axis_coil_center_a[0], p.z_axis_coil_center_a[1], i)) z_axis_field.append(np.abs(B[2])) ax3.plot(np.linspace(-1 * p.x_axis_coil_z_length/2, 1*p.x_axis_coil_z_length/2, 100), z_axis_field) ax3.set_xlabel('z (mm)') ax3.set_ylabel('B (mT)') ax3.set_title('B field along z axis') plt.tight_layout() plt.show() def resistance_sim(p:AqpBookerParameters = get_default_aqp_parameters()): print("coil diameter is", p.wire_diameter, "mm") print("coil thickness is", p.x_axis_coil_thickness, "mm") print("x coil resistance is", p.x_wire_resistance, "ohm") print("y coil resistance is", p.y_wire_resistance, "ohm") print("z coil resistance is", p.z_wire_resistance, "ohm") print("x coil wire length is", p.x_wire_length, "m") def heat_sim(p:AqpBookerParameters = get_default_aqp_parameters()): x_heat_production = p.x_wire_resistance * p.x_axis_coil_current_a**2 print("Each x coil heat production is", p.x_wire_resistance * p.x_axis_coil_current_a**2, "W") print("Each y coil heat production is", p.y_wire_resistance * p.y_axis_coil_current_a**2, "W") # print("z coil heat production is", p.z_wire_resistance * p.z_axis_coil_current_a**2, "W") wire_weight = p.x_wire_length * np.pi * ((p.wire_diameter/2000) ** 2) * 8941 # kg print("x coil wire weight is", wire_weight, "kg") delta_temp = x_heat_production * 1 / 0.385 / (wire_weight * 1000) # assume only 1 second on time print("x coil temperature rise is", delta_temp, "K") def find_lowest_heat_gen(p:AqpBookerParameters = get_default_aqp_parameters()): winding_min = 50 winding_max = 300 current_min = 0.5 current_max = 3.0 sample_points = 50 gradient_limit = 18 # G/cm winding_number = np.linspace(winding_min, winding_max, sample_points, dtype=int) current = np.linspace(current_min, current_max, sample_points) gradient_list = [] for i in range(len(winding_number)): for j in range(len(current)): print("calculating: ", winding_number[i], " windings, ", current[j], " A\n") p.x_axis_coil_winding_number = winding_number[i] p.y_axis_coil_winding_number = winding_number[i] p.x_axis_coil_current_a = -1 * current[j] p.x_axis_coil_current_b = current[j] p.y_axis_coil_current_a = -1 * current[j] p.y_axis_coil_current_b = current[j] coll = build_magpy_collection(p) gradient = get_gradient(coll, np.asarray(p.sample_center), 'x') if np.abs(gradient) > gradient_limit: x_heat_production = p.x_wire_resistance * p.x_axis_coil_current_a**2 gradient_list.append([winding_number[i], current[j], gradient, x_heat_production]) break # export the data to csv np.savetxt("gradient_list_24awg.csv", gradient_list, delimiter=",") print("done") def find_min_B_field(coll: magpy.Collection): search_range = [-1,1] sample_number = 51 x_min = 999 y_min = 999 z_min = 999 x_min_location = (0,0,0) y_min_location = (0,0,0) z_min_location = (0,0,0) for i in np.linspace(search_range[0], search_range[1], sample_number): for j in np.linspace(search_range[0], search_range[1], sample_number): for k in np.linspace(search_range[0], search_range[1], sample_number): B = coll.getB((i, j, k)) if np.abs(B[0]) < np.abs(x_min): x_min = B[0] x_min_location = (i, j, k) if np.abs(B[1]) < np.abs(y_min): y_min = B[1] y_min_location = (i, j, k) if np.abs(B[2]) < np.abs(z_min): z_min = B[2] z_min_location = (i, j, k) print("x min is", x_min, "at", x_min_location) print("y min is", y_min, "at", y_min_location) print("z min is", z_min, "at", z_min_location) return {"x_min": x_min, "y_min": y_min, "z_min": z_min, "x_min_location": x_min_location, "y_min_location": y_min_location, "z_min_location": z_min_location} def run_simulation(iteration, p:AqpBookerParameters = get_default_aqp_parameters()): direction_1 = random.choice(['x', 'y', 'z']) direction_2 = random.choice(['x', 'y', 'z']) direction_3 = random.choice(['x', 'y', 'z']) direction_4 = random.choice(['x', 'y', 'z']) side_1 = random.choice([-1, 1]) side_2 = random.choice([-1, 1]) coll = build_magpy_collection(p) coll.children[0][0].rotate_from_angax(p.angle_tolerance, direction_1) coll.children[0][1].rotate_from_angax(side_1*p.angle_tolerance, direction_2) coll.children[1][0].rotate_from_angax(p.angle_tolerance, direction_3) coll.children[1][1].rotate_from_angax(side_1*p.angle_tolerance, direction_4) result = find_min_B_field(coll) return { "iteration": iteration, "directions": [direction_1, direction_2, direction_3, direction_4], "sides": [side_1, side_2], "result": result } def error_sim_parallel(p:AqpBookerParameters = get_default_aqp_parameters()): iterations = 10 results = [] with ProcessPoolExecutor() as executor: futures = [executor.submit(run_simulation, i, p) for i in range(iterations)] for future in concurrent.futures.as_completed(futures): results.append(future.result()) # Write results to file with open("error_sim_result.txt", "w") as f: for result in results: f.write(json.dumps(result) + "\n") print("Simulation completed.") if __name__ == "__main__": # resistance_sim() # heat_sim() # simulate_mag_field() # heat_generation() # B_field_plot() # find_lowest_heat_gen() error_sim_parallel(get_default_aqp_parameters())