commit d0b900120e03e6982ffefb0c269fe085f94a4815 Author: bukehu Date: Sun Mar 10 02:23:34 2024 +0000 initial diff --git a/coil_sim.py b/coil_sim.py new file mode 100644 index 0000000..a737d2c --- /dev/null +++ b/coil_sim.py @@ -0,0 +1,431 @@ +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()) \ No newline at end of file