initial

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 B_z 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 B_z 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 B_z 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
+
+

src/coil_class.py

@@ -0,0 +1,442 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Aug 17 15:17:44 2021
+
+@author: Joschka
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy import special as sp
+
+from src import physical_constants as cs
+from src import coil_class_add_functions as cf
+
+
+class BCoil:
+    def __init__(self, HH, distance, radius, layers, windings, wire_width, wire_height, layers_spacing=0,
+                 windings_spacing=0):
+        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.layers_spacing = layers_spacing * 1e-3
+        self.windings_spacing = windings_spacing * 1e-3
+
+    def get_wire_area(self):
+        return self.wire_width * self.wire_height
+
+    def get_N(self):
+        return self.layers * self.windings
+
+    def get_wire_length(self):
+        return self.get_N() * 2 * self.radius * np.pi
+
+    def get_coil_height(self):
+        return self.windings * self.wire_height + (self.windings - 1) * self.windings_spacing
+
+    def get_coil_width(self):
+        return self.layers * self.wire_width + (self.layers - 1) * self.layers_spacing
+
+    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_radius(self):
+        return self.radius
+
+    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_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 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} mm, z_max = {self.get_zmax() * 1e3} mm")
+        print(
+            f"Radius = {self.radius * 1e3} mm, Radius_inner = {self.get_R_inner() * 1e3} mm, Radius_outer = {self.get_R_outer() * 1e3:.2f} mm")
+        print(
+            f"layers = {self.layers}, windings = {self.windings}, wire_width = {self.wire_width * 1e3}, wire_height = {self.wire_height * 1e3} mm ")
+        print(" ")
+
+    def print_basic_info(self):
+        print(
+            f"HH = {self.HH}, Distance = {self.distance * 1e3} mm, Radius = {self.radius * 1e3} mm,layers = {self.layers}, windings = {self.windings}")
+
+    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
+
+    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(BCoil.k_sq(R_loop, z_loc, r, z)) + sp.ellipe(BCoil.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(BCoil.k_sq(R_loop, z_loc, r, z)) + sp.ellipe(BCoil.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(self, I_current, x, z, raster=10):
+        """
+        Returns B_z along z-axis and B_x along x-axis,
+        HH = +1 --> Helmholtz configuration, HH = -1 --> Anti Helmholtz configuration
+        """
+
+        z_start = self.get_zmin() + self.wire_height / 2  # (distance_coils/2 - windings * wire_height/2 + wire_height/2)*1e-3
+        R_start = self.get_R_inner() + self.wire_width / 2  # (R_inner + wire_width/2 )*1e-3
+
+        # 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_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
+
+        I_current /= raster ** 2
+
+        for xx in range(0, self.layers):
+            for zz in range(0, self.windings):
+                for xx_in in range(0, raster):
+                    for zz_in in range(0, raster):
+                        z_pos = z_start + zz * (
+                                    self.wire_height + self.windings_spacing) - self.wire_height / 2 + zz_in * self.wire_height / (
+                                            raster - 1)
+                        R_pos = R_start + xx * (
+                                    self.wire_width + self.layers_spacing) - self.wire_width / 2 + xx_in * self.wire_width / (
+                                            raster - 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_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)
+
+        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_tot_along_axis(self, I_current, x, z, raster=10):
+        """
+        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
+        """
+
+        z_start = self.get_zmin() + self.wire_height / 2  # (distance_coils/2 - windings * wire_height/2 + wire_height/2)*1e-3
+        R_start = self.get_R_inner() + self.wire_width / 2  # (R_inner + wire_width/2 )*1e-3
+
+        # 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
+
+        I_current /= raster ** 2
+
+        for xx in range(0, self.layers):
+            for zz in range(0, self.windings):
+                for xx_in in range(0, raster):
+                    for zz_in in range(0, raster):
+                        z_pos = z_start + zz * (
+                                    self.wire_height + self.windings_spacing) - self.wire_height / 2 + zz_in * self.wire_height / (
+                                            raster - 1)
+                        R_pos = R_start + xx * (
+                                    self.wire_width + self.layers_spacing) - self.wire_width / 2 + xx_in * self.wire_width / (
+                                            raster - 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)
+                        # B_z 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)
+
+        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])
+
+        I_current /= raster ** 2
+
+        for el in range(0, len(x_SI)):
+            for xx in range(0, self.layers):
+                for zz in range(0, self.windings):
+                    for xx_in in range(0, raster):
+                        for zz_in in range(0, raster):
+                            z_pos = z_start + zz * (
+                                        self.wire_height + self.windings_spacing) - self.wire_height / 2 + zz_in * self.wire_height / (
+                                                raster - 1)
+                            R_pos = R_start + xx * (
+                                        self.wire_width + self.layers_spacing) - self.wire_width / 2 + xx_in * self.wire_width / (
+                                                raster - 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
+
+    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)
+        z = np.arange(-z_lim, z_lim, 5)
+        x_m, z_m = np.meshgrid(x, z)
+
+        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[xx, zz] > 15:
+                    B[xx, zz, :] /= B_tot[xx, zz] / 15
+
+        plt.figure(33)
+        plt.quiver(x_m, z_m, B[:, :, 0], B[:, :, 1])
+        plt.xlabel("x-axis [mm]")
+        plt.ylabel("z-axis [mm]")
+        plt.show()
+
+    def Bcurv(B_field, z):
+        return np.gradient(np.gradient(B_field, z), z) * 1e2
+
+    def Bgrad(B_field, z):
+        return np.gradient(B_field, z) * 1e1
+
+    def Bz_plot_HH(self, I_current, x, z):
+        B_z, B_x = self.B_multiple(I_current, x, z)
+
+        B_z_curvature = np.gradient(np.gradient(B_z, z), z)
+
+        plt.figure(100, figsize=(13, 10))
+
+        # plt.rcParams.update({'font.size': 15})
+        plt.suptitle("Helmholtz coil field B_z along z-axis")
+
+        # Field plot
+        ##########################
+        plt.subplot(2, 1, 1)
+        plt.plot(z, B_z, linestyle="solid", label=r"$B_z$: Result via elliptic integrals")
+
+        # plt.xlim(-0.01,0.01)
+        plt.title("B-field")
+
+        plt.ylabel(r"$B_z$ [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 B_z$: Result via elliptic integrals")
+
+        plt.ylabel(r"$\nabla_z^2 B_z [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 Bz_plot_HH_comp(self, Coil2, I_current, x, z):
+        B_z, B_x = self.B_multiple(I_current, x, z)
+        B_z_2, B_x_2 = Coil2.B_multiple(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 B_z along z-axis")
+
+        # Field plot
+        ##########################
+        plt.subplot(2, 1, 1)
+        plt.plot(z, B_z, linestyle="solid", label=r"$B_z$")
+        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"$B_z$ [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 B_z$")
+        plt.plot(z, B_z_curvature_2, linestyle="solid", label=r"$\nabla_z^2 B_{z2}$")
+
+        plt.ylabel(r"$\nabla_z^2 B_z [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
+
+        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.resistance(T) * I_current ** 2
+
+        print(f"current density = {j_dens} A/mm^2")
+        print(f"Power = {Power} W")
+
+    def power(self, I_current, T):
+        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 rectangular coil via empirical formular by perry
+
+        """
+        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
+
+    def resistance(self, T):
+        """
+        Calculates resistance of coil configuration in series
+
+        Parameters
+        ----------
+        T : double
+            temperature in degree Celsius
+
+        Returns
+        -------
+        double
+            resistance in Ohm
+
+        """
+        return 2 * cf.resistivity_copper(T) * self.get_wire_length() / self.get_wire_area()
+
+
+HH_Coil = BCoil(HH=1, distance=54, radius=48, layers=4, windings=4, wire_height=1, wire_width=1, windings_spacing=0.25,
+                layers_spacing=0.25)
+HH_Coil.set_R_outer(49.3)
+HH_Coil.set_d_min(49.8)
+
+L = HH_Coil.induct_perry()
+
+HH_Coil.print_info()
+print(L)

src/coil_class_add_functions.py

@@ -0,0 +1,18 @@
+# -*- 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
+
+def resistivity_copper(T):
+    R_ref = cs.rho_copper_20
+    T_ref = 20 #degree celsius
+    alpha = 0.00393
+    R = R_ref *(1 + alpha * (T-T_ref))
+    return R
+
+print(resistivity_copper(30))

src/physical_constants.py

@@ -0,0 +1,31 @@
+# -*- 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.68e-8
+density_copper = 8.940e3
+