Added tokamak field test and cleaned up old code

fusion_toolbox/Bfield.py

@@ -165,7 +165,6 @@ def make_PFset(self,R,Z,I):
         for i in range(N):
             self.coils.append(PFCoil(R[i],Z[i],I[i]))
 
-
     def make_TFset(self,phi,I):
         """
         Adds in a set of TF coils given a list of their toroidal angle phi and currents.
@@ -182,26 +181,6 @@ def make_TFset(self,phi,I):
         for i in range(N):
             self.coils.append(TFCoil(self.R,self.a,phi[i],I[i]))
 
-'''
-thetas = np.linspace(0,2*np.pi,100)
-X = 1*np.cos(thetas)
-Y = 1*np.sin(thetas)
-Z = 0*X
-XYZ = np.vstack((X,Y,Z)).T
-
-a = Coil(XYZ,1)
-field = a.B(np.array([[0,0,0],[0,1,2]]))
-print(field)
-'''
-
-a = PFCoil(1,0,1)
-field = a.B(np.array([[0,0,0],[0,1,2]]))
-print(field)
-
-#CMod = Tokamak(4,1)
-
-#coils = CMod.make_PFset(4,1,1)
-#coils = CMod.make_PFset(4,-1,1)
 def gen_plane_pts(axis,r,w1,w2,n1,n2):
     """
     Generates a 2D grid of points on a plane in 3D which is normal to one of the three axes.

tests/test_magnetics.py

@@ -23,6 +23,8 @@ def generate_pf_coil_xyz(R, N = 100, z = 0):
 
 def analytic_B_center_of_pf_coil(I, R): 
     return np.array([0,0,MU_0 * I / (2 * R)])
+def analytic_B_axis_of_pf_coil(I,R,z):
+    return np.array([0,0,MU_0*2*pi*R**2*I/(4*pi*(z**2+R**2)**(3/2))])
 
 @pytest.mark.parametrize(
     "I, R",
@@ -43,3 +45,25 @@ def test_B_center_of_pf_circular_coil(I, R):
     # Compare to analytic calculation
     B_analytic = [analytic_B_center_of_pf_coil(I,R)]
     npt.assert_almost_equal(B_center, B_analytic, decimal = 4)
+
+@pytest.mark.parametrize(
+    "R1, Z1, I1, R2, Z2, I2",
+    [[1,1,1,1,-1,1],[1,1,1,1,-1,-1]
+    ])
+def test_tokamak_B_total_center(R1,Z1,I1,R2,Z2,I2):
+    '''
+    R1: Radius of first loop [m]
+    Z1: Z offset of first loop [m]
+    I1: Current of first loop [A]
+    R2: Radius of second loop [m]
+    Z2: Z offset of second loop [m]
+    I2: Current of second loop [A]
+    Checks the field exactly centered between two parallel planar coil loops against analytic sol.
+    '''
+    #Generate tokamak with test coil pair (major and minor radius are arbitrary)
+    tok = Tokamak(1,1)
+    tok.make_PFset([R1,R2],[Z1,Z2],[I1,I2])
+    B = tok.get_B_from_coils(np.array([[0.,0,0]]))
+    B_center_analytic = analytic_B_axis_of_pf_coil(I1,R1,Z1-(Z1+Z2)/2)\
+                       +analytic_B_axis_of_pf_coil(I2,R2,-Z1+(Z1+Z2)/2)
+    npt.assert_almost_equal(B[0],B_center_analytic, decimal = 9)