Merge pull request 'first commit' (#1) from master into main

Reviewed-on: #1

coil_design/src/B_field_calculation.py

@@ -0,0 +1,151 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Aug 16 11:51:36 2021
+
+@author: Joschka
+"""
+
+import numpy as np
+from scipy import special as sp
+from src import physical_constants as cs
+
+
+#HH = +1 --> Helmholtz configuration, HH = -1 --> Anti Helmholtz configuration
+def B_field_ideal(N_windings,I_current,R_radius,d_distance,HH,z_arg): 
+    """ Calculate B-field for ideal point like wires with R_radius and d_distance
+    HH = +1 --> Helmholtz configuration, HH = -1 --> Anti Helmholtz configuration
+    """
+    B = cs.mu_0*N_windings*I_current/2 * R_radius**2 * (1/(R_radius**2+(z_arg-d_distance/2)**2)**(3/2) +HH* 1/(R_radius**2+(z_arg+d_distance/2)**2)**(3/2))
+    B *= 1e4 #conversion Gauss
+    return  B
+
+#Argument for elliptic integrals
+def k_sq(R_loop, z_loc, r, z):
+    """ Argument for elliptic integrals"""
+    return (4*R_loop*r)/((R_loop+r)**2 + (z-z_loc)**2)
+
+
+
+def B_z_loop(I_current, R_loop, z_loc, r, z):
+    """calculate z-component of B-field at position r and z for each individual loop"""
+    B_z = 2e-7*I_current * 1/( (R_loop+r)**2 + (z-z_loc)**2 )**(1/2) *(sp.ellipk(k_sq(R_loop,z_loc,r,z))+ sp.ellipe(k_sq(R_loop,z_loc,r,z))*(R_loop**2 - r**2 - (z-z_loc)**2)/((R_loop-r)**2 + (z-z_loc)**2))
+    B_z *= 1e4 #conversion to gauss
+    return B_z
+
+
+def B_r_loop(I_current, R_loop, z_loc, r, z):
+    """calculate r-component of B-field at position r and z for each individual loop""" 
+    B_r = 2e-7*I_current/r * (z-z_loc)/( (R_loop+r)**2 + (z-z_loc)**2 )**(1/2) *(-sp.ellipk(k_sq(R_loop,z_loc,r,z))+ sp.ellipe(k_sq(R_loop,z_loc,r,z))*(R_loop**2 + r**2 + (z-z_loc)**2)/((R_loop-r)**2 + (z-z_loc)**2))
+
+    B_r *= 1e4 #conversion to gauss
+    return B_r
+
+def B_multiple(I_current, HH, R_inner, distance_coils, layers, windings, wire_width, wire_height, x, z):
+    """Returns Bz along z-axis and B_r along r-axis
+    HH = +1 --> Helmholtz configuration, HH = -1 --> Anti Helmholtz configuration"""
+    z_start = (distance_coils/2 - windings * wire_height/2 + wire_height/2)*1e-3
+    R_start = (R_inner + wire_width/2 )*1e-3
+
+    
+    if x[0] <=0:
+        x_neg = np.abs([el for el in x if el<0])
+        x_pos = np.array([el for el in x if el>0])
+        x_zero = np.array([el for el in x if el == 0])
+        
+        
+    B_z = np.zeros(len(z))
+    B_x_pos = np.zeros(len(x_pos))
+    B_x_neg = np.zeros(len(x_neg))
+    B_x_zero = x_zero #If there is a zero in the x-array it is assured that the x component of the field at this point is zero
+    
+    
+    
+    for xx in range(0,layers):
+        for zz in range(0,windings):
+            z_pos = z_start + zz * wire_height*1e-3
+            R_pos = R_start + xx * wire_width*1e-3
+
+            B_z += B_z_loop(I_current,R_pos,z_pos,0,z) + B_z_loop(HH*I_current,R_pos,-z_pos,0,z)            
+            B_x_pos += B_r_loop(I_current,R_pos,z_pos,x_pos,0) + B_r_loop(HH*I_current,R_pos,-z_pos,x_pos,0)
+            B_x_neg += B_r_loop(I_current,R_pos,z_pos,x_pos,0) + B_r_loop(HH*I_current,R_pos,-z_pos,x_pos,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))
+    return B_z,B_x
+
+
+def B_multiple_raster(I_current, HH, R_inner, distance_coils, layers, windings, wire_width, wire_height, x, z):
+    """Returns Bz along z-axis and B_r along r-axis
+    HH = +1 --> Helmholtz configuration, HH = -1 --> Anti Helmholtz configuration"""
+    z_start = (distance_coils/2 - windings * wire_height/2 + wire_height/2)*1e-3
+    R_start = (R_inner + wire_width/2 )*1e-3
+    
+    if x[0] <=0:
+        x_neg = np.abs([el for el in x if el<0])
+        x_pos = np.array([el for el in x if el>0])
+        x_zero = np.array([el for el in x if el == 0])
+        
+        
+    B_z = np.zeros(len(z))
+    B_x_pos = np.zeros(len(x_pos))
+    B_x_neg = np.zeros(len(x_neg))
+    B_x_zero = x_zero #If there is a zero in the x-array it is assured that the x component of the field at this point is zero
+    
+    
+    
+    #dividing each wire into 10 x 10 raster
+    raster = 10
+    I_current /=raster**2
+    for xx in range(0,layers):
+        for zz in range(0,windings):
+            for xx_in in range(0,raster):
+                for zz_in in range(0,raster):
+                    z_pos = z_start + (zz * wire_height - wire_height/2 + zz_in * wire_height/(raster-1) )*1e-3
+                    R_pos = R_start + (xx * wire_width - wire_width/2 + xx_in * wire_width/(raster-1) )*1e-3
+
+                    B_z += B_z_loop(I_current,R_pos,z_pos,0,z) + B_z_loop(HH*I_current,R_pos,-z_pos,0,z)            
+                    B_x_pos += B_r_loop(I_current,R_pos,z_pos,x_pos,0) + B_r_loop(HH*I_current,R_pos,-z_pos,x_pos,0)
+                    B_x_neg += B_r_loop(I_current,R_pos,z_pos,x_pos,0) + B_r_loop(HH*I_current,R_pos,-z_pos,x_pos,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))
+    return B_z,B_x
+
+def B_multiple_raster_test(I_current, HH, R_inner, distance_coils, layers, windings, wire_width, wire_height, x, z):
+    """Returns Bz along z-axis and B_r along r-axis
+    HH = +1 --> Helmholtz configuration, HH = -1 --> Anti Helmholtz configuration"""
+    z_start = (distance_coils/2 - windings * wire_height/2 + wire_height/2)*1e-3
+    R_start = (R_inner + wire_width/2 )*1e-3
+    
+    if x[0] <=0:
+        x_neg = np.abs([el for el in x if el<0])
+        x_pos = np.array([el for el in x if el>0])
+        x_zero = np.array([el for el in x if el == 0])
+        
+        
+    B_z = np.zeros(len(z))
+    B_x_pos = np.zeros(len(x_pos))
+    B_x_neg = np.zeros(len(x_neg))
+    B_x_zero = x_zero #If there is a zero in the x-array it is assured that the x component of the field at this point is zero
+    
+    
+    
+    #dividing each wire into 10 x 10 raster
+    raster = 10
+    I_current /=raster
+    for xx in range(0,layers):
+        for zz in range(0,windings):
+            for xx_in in range(0,1):
+                for zz_in in range(0,raster):
+                    z_pos = z_start + (zz * wire_height - wire_height/2 + zz_in * wire_height/(raster-1) )*1e-3
+                    R_pos = R_start + (xx * wire_width)*1e-3
+
+                    B_z += B_z_loop(I_current,R_pos,z_pos,0,z) + B_z_loop(HH*I_current,R_pos,-z_pos,0,z)            
+                    B_x_pos += B_r_loop(I_current,R_pos,z_pos,x_pos,0) + B_r_loop(HH*I_current,R_pos,-z_pos,x_pos,0)
+                    B_x_neg += B_r_loop(I_current,R_pos,z_pos,x_pos,0) + B_r_loop(HH*I_current,R_pos,-z_pos,x_pos,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))
+    return B_z,B_x
+
+

coil_design/src/coil_class.py

@@ -0,0 +1,940 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Aug 17 15:17:44 2021
+
+@author: Joschka
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+import logging as log
+from scipy import special as sp
+
+import matplotlib
+# matplotlib.use('Qt5Agg')
+# matplotlib.use('Agg')
+# get_ipython().run_line_magic('matplotlib', 'qt')
+# get_ipython().run_line_magic('matplotlib', 'inline')
+# import time
+
+from src import physical_constants as cs
+
+# logger = log.getLogger('example')
+# log.setLevel(log.info)
+log.basicConfig(level=log.WARNING, format='%(message)s')
+
+
+# TODO: Docstrings for every function
+
+# TODO: Implement conventions
+
+class BCoil:
+    def __init__(self, HH, distance, radius, layers, windings, wire_width, wire_height, insulation_thickness=0.
+                 , 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
+        :param distance: distance between center of coils
+        :param radius: radius to center of coil
+        :param layers: number of radial (horizontal) layers
+        :param windings: number of axial (vertical) layers
+        :param wire_width: width of conductive wire core
+        :param wire_height: height of conductive wire core
+        :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:
+                log.warning('Is there a reason you want to wind a not round wire like that?')
+        if is_round:
+            if wire_width != wire_height:
+                log.error('Wire is round but width != height')
+                raise ValueError("Wire is round but width != height")
+
+        self.HH = HH
+        self.distance = distance * 1e-3
+        self.radius = radius * 1e-3
+        self.layers = layers
+        self.windings = windings
+        self.wire_width = wire_width * 1e-3
+        self.wire_height = wire_height * 1e-3
+        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
+    def get_HH(self):
+        return self.HH
+
+    @property
+    def get_distance(self):
+        return self.distance
+
+    @property
+    def get_radius(self):
+        return self.radius
+
+    @property
+    def get_layers(self):
+        return self.layers
+
+    @property
+    def get_windings(self):
+        return self.windings
+
+    @property
+    def get_wire_width(self):
+        return self.wire_width
+
+    @property
+    def get_wire_height(self):
+        return self.wire_height
+
+    @property
+    def get_insulation_thickness(self):
+        return self.insulation_thickness
+
+    @property
+    def get_is_round(self):
+        return self.is_round
+
+    @property
+    def get_winding_scheme(self):
+        return self.winding_scheme
+
+    def get_zmin(self):
+        return self.distance / 2 - self.get_coil_height() / 2
+
+    def get_zmax(self):
+        return self.distance / 2 + self.get_coil_height() / 2
+
+    def get_R_inner(self):
+        return self.radius - self.get_coil_width() / 2
+
+    def get_R_outer(self):
+        return self.radius + self.get_coil_width() / 2
+
+    def set_R_outer(self, R_outer):
+        R_outer *= 1e-3
+        self.radius = R_outer - self.get_coil_width() / 2
+
+    def set_R_inner(self, R_inner):
+        R_inner *= 1e-3
+        self.radius = R_inner + self.get_coil_width() / 2
+
+    def set_d(self, d):
+        d *= 1e-3
+        self.distance = d
+
+    def set_d_min(self, d_min):
+        d_min *= 1e-3
+        self.distance = d_min + self.get_coil_height()
+
+    def set_d_max(self, d_max):
+        d_max *= 1e-3
+        self.distance = d_max - self.get_coil_height()
+
+    def get_wire_area(self):
+        """
+        calculates wire area in m^2
+        :return: float
+        """
+        if self.is_round:
+            return np.pi * (self.wire_width / 2) ** 2
+        return self.wire_width * self.wire_height
+
+    def get_N(self):
+        """
+        Calculates number of windings
+        """
+        N = self.layers * self.windings
+        log.debug(f"N = {N}")
+        if self.winding_scheme == 2:
+            N -= self.layers // 2
+            log.debug(f"N = {N}")
+        return N-self.red_N
+
+    def get_wire_length(self):
+        """
+        :return: Approximate length of wire per coil (m)
+        """
+        return self.get_N() * 2 * self.radius * np.pi
+
+    def get_tot_wire_height(self):
+        """ returns wire height incl. insulation"""
+        return self.wire_height + 2 * self.insulation_thickness
+
+    def get_tot_wire_width(self):
+        """ returns wire width incl. insulation"""
+        return self.wire_width + 2 * self.insulation_thickness
+
+    def get_coil_height(self):
+        if self.winding_scheme == 1:
+            return self.get_tot_wire_height() * (self.windings + 0.5)
+        return self.get_tot_wire_height() * self.windings
+
+    def get_coil_width(self):
+        if self.winding_scheme == 1 or self.winding_scheme == 2:
+            if self.is_round:
+                log.info("Coil width: Be aware of the fact that this is an idealized situation of coil winding (slope "
+                         "of windings changes each layer)")
+                return self.layers * self.get_tot_wire_width() - (self.layers - 1) * (
+                        2 - np.sqrt(3)) * self.get_tot_wire_width() / 2  # width is reduced due to winding offset
+        return self.get_tot_wire_width() * self.layers
+
+    def get_cross_section(self):
+        return self.get_coil_height() * self.get_coil_width()
+
+    def winding_raster(self):
+        """
+        generates raster of flowing currents
+        :param raster_value: wire height/raster_value is distance between rastered points in one wire
+        :return: 2 dim array [[z1,R1], [z2,R2], ...]
+        """
+        outer_raster = np.zeros((self.get_N(), 2))
+        it = 0
+        z_start = self.get_zmin() + self.get_tot_wire_height() / 2  # (distance_coils/2 - windings * wire_height/2 + wire_height/2)*1e-3
+        R_start = self.get_R_inner() + self.get_tot_wire_width() / 2  # (R_inner + wire_width/2 )
+        log.debug(f"N = {self.get_N()}")
+        for xx in range(0, self.layers):
+            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()
+
+                # correct for different winding techniques and round wire
+                if self.winding_scheme == 1 or self.winding_scheme == 2:
+                    if xx % 2 == 1:
+                        z_pos += self.get_tot_wire_height() / 2
+                    if self.is_round:
+                        R_pos -= xx * ((2 - np.sqrt(3)) * self.get_tot_wire_width() / 2)
+                log.debug(f"lay = {xx}, wind = {zz}")
+                outer_raster[it] = [z_pos, R_pos]
+                # print(outer_raster[it])
+                it += 1
+
+        return outer_raster
+
+    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
+        if raster_value == 0:
+            raster_value = 1
+            log.info("raster value is 0 increased to 1")
+
+        if raster_value == 1:
+            return [0, 0]
+
+        if self.is_round & (raster_value % 2 == 0):
+            raster_value += 1
+            log.info(
+                f"for round wire rastering works better with uneven rastering value --> rastering value set to: {raster_value}")
+
+        inner_raster = np.zeros((raster_value ** 2, 2))
+        it = 0
+        for xx_in in range(0, raster_value):
+            for zz_in in range(0, raster_value):
+                z_pos_in = - self.wire_height / 2 + zz_in * self.wire_height / (raster_value - 1)
+                r_pos_in = - self.wire_width / 2 + xx_in * self.wire_width / (raster_value - 1)
+
+                inner_raster[it] = [z_pos_in, r_pos_in]
+                it += 1
+
+        # delete points out of round geometry
+        if self.is_round:
+            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 - 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, 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, 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)):
+            full_ras[i] = outer_ras[i] + inner_ras
+
+        return full_ras
+
+    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:
+            extension = self.get_coil_height()
+        extension *= 1e3
+
+        plt.figure(77, figsize=(5, 5))
+        plt.scatter(full_structure[:, :, 1], full_structure[:, :, 0], linewidths=0.001, s=.1)
+        plt.xlabel("radius [mm]")
+        plt.ylabel("z position [mm]")
+        plt.axvline(x=self.get_R_inner() * 1e3 - 0.1, lw=5, color="red")
+        plt.xlim(1e3 * self.get_R_inner() - 0.5, 1e3 * self.get_R_inner() + extension + 0.5)
+        plt.ylim(1e3 * self.get_zmin() - 0.5, 1e3 * self.get_zmin() + extension + 0.5)
+
+        plt.show()
+
+    def print_info(self):
+        print(" ")
+        # print(f"{self.get_zmin()}")
+        # print(self.__name__)
+        print(
+            f"HH = {self.HH}, Distance = {self.distance * 1e3} mm, z_min = {self.get_zmin() * 1e3:.3f} mm, z_max = {self.get_zmax() * 1e3:.3f} mm")
+        print(
+            f"Radius = {self.radius * 1e3:.3f} mm, Radius_inner = {self.get_R_inner() * 1e3:.3f} mm, Radius_outer = {self.get_R_outer() * 1e3:.3f} mm")
+        print(
+            f"Coil width = {self.get_coil_width() * 1e3:.3f} mm, Coil height = {self.get_coil_height() * 1e3:.3f} mm"
+        )
+        print(
+            f"layers = {self.layers}, windings = {self.windings}, wire_width = {self.wire_width * 1e3:.3f}, wire_height = {self.wire_height * 1e3:.3f} mm, round wire = {self.is_round} ")
+        print(" ")
+
+    def print_basic_info(self):
+        if self.HH == 1:
+            print("Offset coil:")
+        else:
+            print("Gradient coil:")
+        print(
+            f"    layers = {self.layers}, windings = {self.windings} \n "
+            f"   wire: core = {self.wire_height * 1e3:.2f} mm, tot = {self.get_tot_wire_height() * 1e3:.2f} mm, total length for pair = {2 * self.get_wire_length():.0f} m"
+        )
+        print(
+            f"    Distance = {self.distance * 1e3:.2f} mm, Radius = {self.radius * 1e3:.2f} mm"
+        )
+        print(
+              f"    Coil width = {self.get_coil_width() * 1e3:.2f} mm, Coil height = {self.get_coil_height() * 1e3:.2f} mm"
+        )
+        print("")
+
+    def B_field_ideal_z(self, I_current, z_arg):
+        """ 
+        Calculate B-field for ideal point like wires with R_radius and d_distance
+        HH = +1 --> Helmholtz configuration, HH = -1 --> Anti Helmholtz configuration
+        z_arg in mm
+        """
+        z_SI = z_arg * 1e-3
+        N_windings = self.layers * self.windings
+
+        B = cs.mu_0 * N_windings * I_current / 2 * self.radius ** 2 * (
+                1 / (self.radius ** 2 + (z_SI - self.distance / 2) ** 2) ** (3 / 2) + self.HH * 1 / (
+                self.radius ** 2 + (z_arg + self.distance / 2) ** 2) ** (3 / 2))
+        B *= 1e4  # conversion Gauss
+        return B
+
+    @staticmethod
+    def make_axis(lim,precision,two_axis=True):
+        """
+        Gives back two arrays x and z in mm, symmetric around zero
+        Args:
+            lim: positive (and negative if lim2 = None) limit of axis
+            precision: nr_points per mm
+            two_axis: If True --> gives back two arrays with given parameter
+
+        Returns: array
+
+        """
+        nr_points = (2 * lim) * precision + 1
+        nr_points = int(nr_points)
+        z = np.linspace(-lim, lim, nr_points)
+
+        if two_axis:
+            x = np.linspace(-lim, lim, nr_points)
+            return x, z
+        return z
+
+    @staticmethod
+    def mm_it(x, x_pos):
+        """
+        Takes argument in mm returns position in x-array
+        Args:
+            x: array of interest
+            x_pos: position to find position of
+
+        Returns: int of position of x_pos in array x
+
+        """
+        precision = len(x)//(2*np.abs(x[0]))
+        it = int(len(x)//2 + precision * x_pos)
+        return it
+
+    @staticmethod
+    def k_sq(R_loop, z_loc, r, z):
+        """ Argument for elliptic integrals"""
+        return (4 * R_loop * r) / ((R_loop + r) ** 2 + (z - z_loc) ** 2)
+
+    @staticmethod
+    def __B_z_loop(I_current, r_loop, z_loop, r_pos, z_pos):
+        """
+        calculate z-component of B-field at position r and z for each individual loop
+        Args:
+            I_current: Current in [A]
+            r_loop: Radial position of current loop
+            z_loop: Axial position of current loop along symmetry axis
+            r_pos: radial position of calculation in [m]
+            z_pos: axial position of calculation in [m]
+
+        Returns: z component of B-field at position r, z in G
+
+        """
+        B_z = 2e-7 * I_current * 1 / ((r_loop + r_pos) ** 2 + (z_pos - z_loop) ** 2) ** (1 / 2) * (
+                sp.ellipk(BCoil.k_sq(r_loop, z_loop, r_pos, z_pos)) + sp.ellipe(BCoil.k_sq(r_loop, z_loop, r_pos, z_pos)) * (
+                r_loop ** 2 - r_pos ** 2 - (z_pos - z_loop) ** 2) / ((r_loop - r_pos) ** 2 + (z_pos - z_loop) ** 2))
+        B_z *= 1e4  # conversion to gauss
+        return B_z
+
+    @staticmethod
+    def __B_r_loop(I_current, r_loop, z_loop, r_pos, z_pos):
+        """
+         calculate z-component of B-field at position r and z for each individual loop
+         Args:
+             I_current: Current in [A]
+             r_loop: Radial position of current loop
+             z_loop: Axial position of current loop along symmetry axis
+             r_pos: radial position of calculation in [m]
+             z_pos: axial position of calculation in [m]
+
+         Returns: z component of B-field at position r, z in G
+
+         """
+        B_r = 2e-7 * I_current / r_pos * (z_pos - z_loop) / ((r_loop + r_pos) ** 2 + (z_pos - z_loop) ** 2) ** (1 / 2) * (
+                -sp.ellipk(BCoil.k_sq(r_loop, z_loop, r_pos, z_pos)) + sp.ellipe(BCoil.k_sq(r_loop, z_loop, r_pos, z_pos)) * (
+                r_loop ** 2 + r_pos ** 2 + (z_pos - z_loop) ** 2) / ((r_loop - r_pos) ** 2 + (z_pos - z_loop) ** 2))
+        B_r *= 1e4  # conversion to gauss
+        return B_r
+
+    def B_field_z(self,I_current):
+        pass
+
+    def B_field(self, I_current, x, z, raster = 7):
+        """
+        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,
+                                                                                             -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, -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, -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_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
+        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
+
+        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_0 += BCoil.__B_z_loop(I_current, r_pos, z_pos, 0, 0) + BCoil.__B_z_loop(self.HH * I_current, r_pos,
+                                                                                            -z_pos, 0, 0)
+        print(f"    Max Field = {B_z_0:.2f} G")
+        return B_z_0
+
+    def max_field_one_coil(self, I_current, raster = 7):
+        """
+        Returns and prints max (z-) field in Gauss at center of one coil of configuration pair if second one is not powered
+        Args:
+            I_current: Current in [A]
+            raster: rastering value for each wire (increasing increases accuracy), 7 is fine for most cases
+
+        Returns:
+            Max B-field at center position of coil
+
+        """
+        B_z_0 = 0
+        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
+
+        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_0 += BCoil.__B_z_loop(I_current, r_pos, z_pos, 0, z_pos)
+
+        print(f"    Max Field = {B_z_0:.2f} G")
+        return B_z_0
+
+    def max_gradient(self, i_current, raster = 7):
+        x, z = BCoil.make_axis(0.5, 1000)
+        Bz, Bx = self.B_field(i_current, x, z, raster = raster)
+
+        bz_grad = BCoil.grad(Bz, z)
+        bx_grad = BCoil.grad(Bx, x)
+        print(f"    Max z Gradient = {bz_grad[len(z)//2]:.2f} G/cm")
+        print(f"    Max x Gradient = {bx_grad[len(x)//2]:.2f} G/cm")
+        # return bx_grad[len(x)//2]
+
+
+    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,
+        HH = +1 --> Helmholtz configuration, HH = -1 --> Anti Helmholtz configuration
+        """
+
+        # Convert axis to SI units
+        x_SI = x * 1e-3
+        z_SI = z * 1e-3
+
+        if x[0] <= 0:
+            x_neg = np.abs([el for el in x_SI if el < 0])
+            x_pos = 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_tot_z = np.zeros(len(z))
+        B_x_pos = np.zeros(len(x_pos))
+        B_x_neg = np.zeros(len(x_neg))
+        B_x_zero = x_zero  # If there is a zero in the x-array it is assured that the x component of the field at this point is zero
+
+        B_z_x = np.zeros(len(x_SI))
+        # dividing each wire into 10 x 10 raster
+        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]
+
+                # z-field along z-axis (x-Field always zero)
+                B_tot_z += BCoil.__B_z_loop(I_current, r_pos, z_pos, 0, z_SI) + BCoil.__B_z_loop(self.HH * I_current,
+                                                                                                 r_pos, -z_pos, 0, z_SI)
+
+                # x-field along x-axis
+                B_x_pos += BCoil.__B_r_loop(I_current, r_pos, z_pos, x_pos, 0) + BCoil.__B_r_loop(self.HH * I_current,
+                                                                                                  r_pos, -z_pos, x_pos,
+                                                                                                  0)
+                B_x_neg += BCoil.__B_r_loop(I_current, r_pos, z_pos, x_neg, 0) + BCoil.__B_r_loop(self.HH * I_current,
+                                                                                                  r_pos, -z_pos, x_neg,
+                                                                                                  0)
+                # Bz along x-axis:
+                B_z_x += BCoil.__B_z_loop(I_current, r_pos, z_pos, x_SI, 0) + BCoil.__B_z_loop(self.HH * I_current,
+                                                                                               r_pos, -z_pos, x_SI, 0)
+
+        B_x_neg *= -1  # B_x_neg is pointing in opposite x_direction
+
+        # B_x_x = np.zeros(len(x_SI))
+        B_x_x = np.concatenate((B_x_neg, B_x_zero, B_x_pos))
+        B_tot_x = np.sqrt(B_x_x ** 2 + B_z_x ** 2)
+        B_tot_z = np.abs(B_tot_z)
+        return B_tot_z, B_tot_x
+
+    def B_multiple_3d(self, I_current, x, z, raster=4):
+        z_start = self.get_zmin() + self.wire_height / 2
+        R_start = self.get_R_inner() + self.wire_width / 2
+
+        x_SI = x * 1e-3
+        z_SI = z * 1e-3
+
+        B = np.zeros([len(z_SI), len(x_SI), 2])
+
+        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])
+
+        # TODO: why division by zero?
+
+        I_current /= rastering_value  # divide current into smaller currents for mapping the whole wire
+        # start = time.time()
+        for el in range(0, len(x_SI)):
+            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]
+
+                    # compute z-value of field
+                    B[:, el, 1] += BCoil.__B_z_loop(I_current, r_pos, z_pos, np.abs(x_SI[el]), z_SI) + BCoil.__B_z_loop(
+                        self.HH * I_current, r_pos, -z_pos, np.abs(x_SI[el]), z_SI)
+
+                    # compute x-value
+                    if x_SI[el] < 0:
+                        B[:, el, 0] -= BCoil.__B_r_loop(I_current, r_pos, z_pos, -x_SI[el], z_SI) + BCoil.__B_r_loop(
+                            self.HH * I_current, r_pos, -z_pos, -x_SI[el], z_SI)
+
+                    elif x_SI[el] > 0:
+                        B[:, el, 0] += BCoil.__B_r_loop(I_current, r_pos, z_pos, x_SI[el], z_SI) + BCoil.__B_r_loop(
+                            self.HH * I_current, r_pos, -z_pos, x_SI[el], z_SI)
+
+                    elif x_SI[el] == 0:
+                        B[:, el, 0] = 0
+
+        return B
+
+    @staticmethod
+    def B_tot_3d(B):
+        return np.sqrt(B[:, :, 0] ** 2 + B[:, :, 1] ** 2)
+
+    def plot_3d(self, I_current, x_lim, z_lim):
+        x = np.arange(-x_lim, x_lim, 5)
+        print(x)
+        z = np.arange(-z_lim, z_lim, 5)
+        print(z)
+        x_m, z_m = np.meshgrid(x, z)
+        if self.is_round:
+            B = self.B_multiple_3d(I_current, x, z, raster=3)
+        else:
+            B = self.B_multiple_3d(I_current, x, z, raster=2)
+        B_tot = BCoil.B_tot_3d(B)
+        for xx in range(0, len(x)):
+            for zz in range(0, len(z)):
+                if B_tot[zz, xx] > 15:
+                    B[zz, xx, :] /= B_tot[zz, xx] / 15
+
+        plt.figure()
+        plt.quiver(x_m, z_m, B[:, :, 0], B[:, :, 1])
+        plt.xlabel("x-axis [mm]")
+        plt.ylabel("z-axis [mm]")
+        plt.show()
+
+    @staticmethod
+    def curv(B_f, z):
+        return np.gradient(np.gradient(B_f, z), z) * 1e2
+
+    @staticmethod
+    def grad(B_f, z):
+        return np.gradient(B_f, z) * 1e1
+
+    def B_quick_plot(self, I_current, abs=True, x_lim=50, z_lim=50, nr_points=200):
+        x = np.linspace(-x_lim, x_lim, nr_points)
+        z = np.linspace(-z_lim, z_lim, nr_points)
+        if abs:
+            B_z, B_x = self.B_tot_along_axis(I_current, x, z)
+        else:
+            B_z, B_x = self.B_field(I_current, x, z)
+
+        plt.figure(11)
+        plt.plot(z, B_z, linestyle="solid", label=r"$B_{tot}$ along z-axis")
+        plt.plot(x, B_x, label=r"$B_{tot}$ along x-axis")
+        plt.title("B-field")
+        plt.ylabel(r"B-field [G]")
+        plt.xlabel("x-axis / z-axis [mm]")
+        plt.legend()
+        plt.show()
+
+    def B_grad_quick_plot(self, I_current, x_lim=50, z_lim=50, nr_points=200):
+        x = np.linspace(-x_lim, x_lim, nr_points)
+        z = np.linspace(-z_lim, z_lim, nr_points)
+
+        B_z, B_x = self.B_field(I_current, x, z)
+        # TODO: Think about Gradient calculation of tot B_field
+        B_z = BCoil.grad(B_z, z)
+        B_x = BCoil.grad(B_x, x)
+
+        plt.figure(12)
+        plt.plot(z, B_z, linestyle="solid", label=r"z grad of Bz along z-axis")
+        plt.plot(x, B_x, label=r"x Grad of B_x along x-axis")
+        plt.title("Gradient of B-field")
+        plt.ylabel(r"B-field [G/cm]")
+        plt.xlabel("x-axis / z-axis [mm]")
+        plt.legend()
+        plt.show()
+
+    def B_curv_quick_plot(self, I_current, x_lim=50, z_lim=50, nr_points=200):
+        x = np.linspace(-x_lim, x_lim, nr_points)
+        z = np.linspace(-z_lim, z_lim, nr_points)
+
+        B_z, B_x = self.B_field(I_current, x, z)
+        B_z = BCoil.curv(B_z, z)
+        B_x = BCoil.curv(B_x, x)
+
+        plt.figure(13)
+        plt.plot(z, B_z, linestyle="solid", label=r"$B_{tot}$ along z-axis")
+        plt.plot(x, B_x, label=r"$B_{tot}$ along x-axis")
+        plt.title("Curvature of B-field")
+        plt.ylabel(r"Curv. B-field [G/cm$^2$]")
+        plt.xlabel("x-axis / z-axis [mm]")
+        plt.legend()
+        plt.show()
+
+    def Bz_plot_HH_comp(self, Coil2, I_current, x, z):
+        B_z, B_x = self.B_field(I_current, x, z)
+        B_z_2, B_x_2 = Coil2.B_field(I_current, x, z)
+
+        B_z_curvature = np.gradient(np.gradient(B_z, z), z) * 1e2
+        B_z_curvature_2 = np.gradient(np.gradient(B_z_2, z), z) * 1e2
+
+        plt.figure(100, figsize=(13, 10))
+
+        # plt.rcParams.update({'font.size': 15})
+        plt.suptitle("Helmholtz coil field Bz along z-axis")
+
+        # Field plot
+        ##########################
+        plt.subplot(2, 1, 1)
+        plt.plot(z, B_z, linestyle="solid", label=r"$Bz$")
+        plt.plot(z, B_z_2, linestyle="solid", label=r"$B_{z2}$")
+        # plt.xlim(-0.01,0.01)
+        plt.title("B-field")
+
+        plt.ylabel(r"$Bz$ [G]")
+        plt.xlabel("z-axis [mm]")
+        plt.legend()
+
+        plt.subplot(2, 1, 2)
+        plt.plot(z, B_z_curvature, linestyle="solid", label=r"$\nabla_z^2 Bz$")
+        plt.plot(z, B_z_curvature_2, linestyle="solid", label=r"$\nabla_z^2 B_{z2}$")
+
+        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(loc='lower right')
+
+        plt.show()
+
+    def cooling(self, I_current, T):
+        """
+        Print current density and power for one coil
+
+        Parameters
+        ----------
+        I_current : current in A
+        T : temperature in degree Celsius
+
+        Returns
+        -------
+        None.
+
+        """
+
+        j_dens = I_current / self.get_wire_area() * 1e-6
+
+        Power = self.power(I_current, T)
+
+        Voltage = self.resistance(T) * I_current
+
+        #print("")
+        print("Cooling:")
+        print(f"    Current I = {I_current} A")
+        print(f"    current density = {j_dens:.2f} A/mm^2")
+        print(f"    Power = {Power:.2f} W")
+        print(f"    U = {Voltage:.2f} V")
+
+    def power(self, I_current, T):
+        """
+        Power for one coil
+        """
+        P = self.resistance(T) * I_current ** 2
+        return P
+
+    def inductivity(self):
+        return cs.mu_0 * (self.layers * self.windings) ** 2 * self.radius ** 2 * np.pi / self.get_coil_height()
+
+    def induct_perry(self):
+        """
+        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 / (
+                0.2317 * 100 * self.radius + 0.44 * 100 * self.get_coil_height() + 0.39 * 100 * self.get_coil_width())
+        return L * 1e-9
+
+    @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
+        R = R_ref * (1 + alpha * (T - T_ref))
+        return R
+
+    def resistance(self, tempr):
+        """
+        Calculates resistance of one coil of configuration
+
+        Parameters
+        ----------
+        tempr : double
+            temperature in degree Celsius
+
+        Returns
+        -------
+        double
+            resistance in Ohm
+        """
+
+        return BCoil.resistivity_copper(tempr) * self.get_wire_length() / self.get_wire_area()
+
+    def tau(self, tempr = 22):
+        return self.induct_perry()/self.resistance(tempr)
+
+
+def main():
+    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()
+if __name__ == "__main__":
+    log.info("Run main Coil class")
+    main()

coil_design/src/coil_class_add_functions.py

@@ -0,0 +1,13 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Sep 15 12:01:05 2021
+
+@author: Joschka
+"""
+import numpy as np
+
+from src import physical_constants as cs
+
+
+
+print(resistivity_copper(30))

coil_design/src/coil_helpers.py

@@ -0,0 +1,224 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from src import coil_class as BC
+
+import magpylib as magpy
+import matplotlib.ticker as ticker
+
+
+def circular_HH(HH, radius, distance, layers, windings,
+                I = 1, #A
+                wire_height = .5 *1e-3,
+                wire_width = .5 *1e-3,
+                insulation_thickness = (0.568-0.5)/2 *1e-3,
+                print_info = True):
+    """returns magpylib system for a circular HH-coil (parameters in m) and prints info on it if requested"""
+
+    #### give out coil parameter s####
+    HH_Coil = BC.BCoil(HH = HH, distance = distance*1e3, radius = radius*1e3, layers = layers, windings = windings,
+                    wire_height = wire_height*1e3, wire_width = wire_width*1e3, insulation_thickness = insulation_thickness*1e3,
+                    is_round = True, winding_scheme= 2)
+
+    if print_info:
+        print("Single coil parameters:")
+
+        HH_Coil.print_info()
+        print(f"N = {HH_Coil.get_N()}")
+        print(f"wire length = {HH_Coil.get_wire_length()} m")
+
+        HH_Coil.plot_raster()
+
+        #### total inductance from Lenny's code ####
+        def Self_Inductance_Abbot(alpha, beta, gamma, N):
+            """
+            Formula for self inductance of circular coil by J.J.Abbott
+            """
+            L = (7.9e-6 * alpha**2 * N**2) / (3*alpha + 9*beta + 10*gamma)
+            return L
+
+        alpha = 2*radius
+        beta = HH_Coil.get_coil_width()
+        gamma = HH_Coil.get_coil_height()
+        N = HH_Coil.get_N()
+        dist = distance
+
+        L_pair = Self_Inductance_Abbot(alpha, dist + beta, gamma, N) +Self_Inductance_Abbot(alpha, dist - beta, gamma, N) - 2 * Self_Inductance_Abbot(alpha, dist, gamma, N) + 2 * Self_Inductance_Abbot(alpha, beta, gamma, N)
+        R_pair = 2*HH_Coil.get_wire_length() * 0.08708
+
+        print("Pair in series:")
+        print(f"wire length = {2*HH_Coil.get_wire_length()} m")
+        print(f"R = {2*HH_Coil.resistance(25)} Ohm")
+        print(f"L = {L_pair*1e3} mH")
+        print(f"tau = {L_pair/R_pair*1e3} ms")
+
+        print(f"\nat {I} A:")
+        print(f"V = {R_pair*I} V")
+        print(f"P = {R_pair*I**2} W")
+
+    #### calculate field properties ####
+
+    #define coils
+    winding_offsets = HH_Coil.full_raster(raster_value=1)[:,0]
+
+    coil1 = magpy.Collection()
+    for coord in winding_offsets:
+        winding = magpy.current.Circle(
+            current=1,
+            diameter=2*coord[1],
+            position=(0,0,coord[0]),
+        )
+        coil1.add(winding)
+
+    coil2 = magpy.Collection()
+    for coord in winding_offsets:
+        winding = magpy.current.Circle(
+            current=HH*1,
+            diameter=2*coord[1],
+            position=(0,0,-coord[0]),
+        )
+        coil2.add(winding)
+
+    #move to test if slight changes in positioning have a large effect
+    #coil1.move([0.5e-3,0,0])
+
+    # Combine them into a collection (like multiple windings)
+    system = magpy.Collection(coil1, coil2)
+
+    return system, 2*HH_Coil.get_wire_length()
+
+
+def rect_current(current, a,b,position, corner_radius):
+    """returns rectangular mapy current with given dimensions (in m)"""
+    #first quadrant
+    angles = np.linspace(0,np.pi/2,50)
+    radii = np.zeros((len(angles),3))
+    radii[:,0] = np.cos(angles)
+    radii[:,1] = np.sin(angles)
+    radii *= corner_radius
+    vertices = np.array([a/2-corner_radius, b/2-corner_radius, 0]) + radii
+
+    #second quadrant
+    angles = np.linspace(np.pi/2,np.pi,50)
+    radii = np.zeros((len(angles),3))
+    radii[:,0] = np.cos(angles)
+    radii[:,1] = np.sin(angles)
+    radii *= corner_radius
+    vertices = np.append(vertices,np.array([-(a/2-corner_radius), b/2-corner_radius, 0]) + radii, axis=0)
+
+    #third quadrant
+    angles = np.linspace(np.pi,1.5*np.pi,50)
+    radii = np.zeros((len(angles),3))
+    radii[:,0] = np.cos(angles)
+    radii[:,1] = np.sin(angles)
+    radii *= corner_radius
+    vertices = np.append(vertices,np.array([-(a/2-corner_radius), -(b/2-corner_radius), 0]) + radii, axis=0)
+
+    #fourth quadrant
+    angles = np.linspace(1.5*np.pi, 2*np.pi,50)
+    radii = np.zeros((len(angles),3))
+    radii[:,0] = np.cos(angles)
+    radii[:,1] = np.sin(angles)
+    radii *= corner_radius
+    vertices = np.append(vertices,np.array([a/2-corner_radius, -(b/2-corner_radius), 0]) + radii, axis=0)
+
+    #close loop
+    vertices = np.append(vertices, [vertices[0]], axis=0)
+
+    # Create a rectangular current coil
+    coil1 = magpy.current.Polyline(position=position, vertices=vertices, current=current)
+    return coil1
+
+
+def rectangular_HH(HH, a, b, distance, corner_radius, layers, windings,
+                    I = 1, #A
+                    wire_height = .5 *1e-3,
+                    wire_width = .5 *1e-3,
+                    insulation_thickness = (0.568-0.5)/2 *1e-3,
+                    print_info = True):
+    """returns magpylib system for a rectangular HH-coil (parameters in m) and prints info on it if requested"""
+
+    #### give out coil parameter s####
+    HH_Coil = BC.BCoil(HH = HH, distance = 0, radius = 0, layers = layers, windings = windings,
+                    wire_height = wire_height*1e3, wire_width = wire_width*1e3, insulation_thickness = insulation_thickness*1e3,
+                    is_round = True, winding_scheme= 2)
+
+    if print_info:
+        N = HH_Coil.get_N()
+        print("Single coil parameters: \n")
+        #HH_Coil.print_info()
+        print(f"coil width = {HH_Coil.get_coil_width()*1e3} mm")
+        print(f"coil height = {HH_Coil.get_coil_height()*1e3} mm")
+        print(f"N = {HH_Coil.get_N()}")
+
+        wire_length = N*(2*a - 4*corner_radius + 2*b - 4*corner_radius + 2*np.pi*corner_radius)
+        print(f"wire length = {wire_length} m")
+
+        HH_Coil.plot_raster()
+
+
+        #### total inductance from Lenny's code ####
+        #convert to centimeters
+        a *= 100
+        b *= 100
+        distance *= 100
+        p1 = (a + np.sqrt(a ** 2 + distance ** 2)) * np.sqrt(b ** 2 + distance ** 2) / ((a + np.sqrt(a ** 2 + b ** 2 + distance ** 2)) * distance)
+        p2 = (b + np.sqrt(b ** 2 + distance ** 2)) * np.sqrt(a ** 2 + distance ** 2) / ((b + np.sqrt(a ** 2 + b ** 2 + distance ** 2)) * distance)
+        p3 = np.sqrt(a ** 2 + b ** 2 + distance ** 2) - np.sqrt(a ** 2 + distance ** 2) - np.sqrt(b ** 2 + distance ** 2) + distance
+        M = 4 * (a * np.log(p1) + b * np.log(p2)) + 8 * p3
+        M *= N**2 * 1e-9 #Conversion to SI and accounting for N wires
+
+        r = (wire_width/2 + insulation_thickness) * 100
+        diag = np.sqrt(a**2 + b**2)
+        L_self = 4 * ( (a+b) * np.log(2*a*b / r) - a * np.log(a+diag) - b * np.log(b+diag) - 7/4 * (a+b) + 2 * (diag + r) )
+        L_self *= N**2 * 1e-9 #Conversion to SI and accounting for N wires
+
+        L_pair = 2 * (L_self+M)
+
+        R_pair = 2*wire_length * 0.08708
+
+        print("Pair in series:")
+        print(f"wire length = {2*wire_length} m")
+        print(f"R = {R_pair} Ohm")
+        print(f"L = {L_pair*1e3} mH")
+        print(f"tau = {L_pair/R_pair *1e3} ms")
+
+        print(f"\nat {I} A:")
+        print(f"V = {R_pair*I} V")
+        print(f"P = {R_pair*I**2} W")
+
+    #### calculate field properties ####
+        #convert back to meters
+        a /= 100
+        b /= 100
+        distance /= 100
+
+    #define coils
+    winding_offsets = HH_Coil.full_raster(raster_value=1)[:,0]
+
+    coil1 = magpy.Collection()
+    for coord in winding_offsets:
+        winding = rect_current(
+            current=1,
+            a=a+coord[1], b=b+coord[1],
+            position=[0,0,distance/2+coord[0]],
+            corner_radius=corner_radius+coord[1]/2
+        )
+        coil1.add(winding)
+
+    coil2 = magpy.Collection()
+    for coord in winding_offsets:
+        winding = rect_current(
+            current=HH*1,
+            a=a+coord[1], b=b+coord[1],
+            position=[0,0,-(distance/2+coord[0])],
+            corner_radius=corner_radius+coord[1]/2
+        )
+        coil2.add(winding)
+
+    #move to test if slight changes in positioning have a large effect
+    #coil1.move([0.5e-3,0,0])
+
+    # Combine them into a collection (like multiple windings)
+    system = magpy.Collection(coil1, coil2)
+
+    return system, 2*wire_length

coil_design/src/physical_constants.py

@@ -0,0 +1,34 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Aug 16 11:54:40 2021
+
+@author: Joschka
+"""
+import numpy as np
+#once and forever
+
+mu_0 = 1.25663706212e-6
+
+mu_B = 9.2740100783e-24
+
+h = 6.62607015e-34
+h_bar = h/(2*np.pi)
+
+k_B = 1.38064852e-23
+
+sigma_B = 5.670374419e-8
+
+#Dysprosium
+m_Dy_164 = 2.72210785e-25
+
+Gamma_626 = 8.5e5
+f_626 = 478.839e12
+lambda_626 = 626.082e-9
+
+#specific constants
+rho_copper_20 = 1.72e-8
+density_copper = 8.940e3
+
+# water
+water_dens = 999.78  # kg/m^3
+water_c_p = 4190  # J/kg*K