Implement potential_radial_ring in terms of delta_r and delta_z

tests/test_radial_ring.py

@@ -67,19 +67,19 @@ class TestRadialRing(unittest.TestCase):
 
     def test_potential_radial_ring(self):
         for r0, z0, r, z in np.random.rand(10, 4):
-            assert np.isclose(B.potential_radial_ring(r0, z0, r, z)/epsilon_0, potential_of_ring_arbitrary(1., r0, z0, r, z))
+            assert np.isclose(B.potential_radial_ring(r0, z0, r-r0, z-z0)/epsilon_0, potential_of_ring_arbitrary(1., r0, z0, r, z))
 
     def test_potential_field_radial_single_ring(self):
         rs = 1 # Source point
         zs = 0
 
         r = np.linspace(1.1, 2, 10000)
-        pot = [B.potential_radial_ring(r_, zs, rs, zs) for r_ in r]
+        pot = [B.potential_radial_ring(r_, zs, rs-r_, 0.) for r_ in r]
         deriv = [B.dr1_potential_radial_ring(r_, zs, rs-r_, 0.) for r_ in r]
         assert np.allclose(CubicSpline(r, pot)(r, 1), deriv, atol=0., rtol=1e-9)
 
         z = np.linspace(0.1, 1, 10000)
-        pot = [B.potential_radial_ring(rs, z_, rs, zs) for z_ in z]
+        pot = [B.potential_radial_ring(rs, z_, 0., zs-z_) for z_ in z]
         deriv = [B.dz1_potential_radial_ring(rs, z_, 0., zs-z_) for z_ in z]
         assert np.allclose(CubicSpline(z, pot)(z, 1), deriv, atol=0., rtol=1e-9)
 
@@ -108,7 +108,7 @@ def test_axial(self):
         correct = k * Q / np.sqrt(z**2 + r**2)
 
         pot = [potential_of_ring_arbitrary(dq, 0., z0, r, 0.) for z0 in z]
-        traceon = [dq/epsilon_0*B.potential_radial_ring(0., z0, r, 0.) for z0 in z]
+        traceon = [dq/epsilon_0*B.potential_radial_ring(0., z0, r, 0.-z0) for z0 in z]
 
         assert np.allclose(pot, correct)
         assert np.allclose(traceon, correct)

traceon/backend/radial.c

@@ -112,7 +112,7 @@ potential_radial(double point[3], double* charges, jacobian_buffer_2d jacobian_b
 	for(int i = 0; i < N_vertices; i++) {  
 		for(int k = 0; k < N_QUAD_2D; k++) {
 			double *pos = &position_buffer[i][k][0];
-			double potential = potential_radial_ring(r0, z0, pos[0], pos[1], NULL);
+			double potential = potential_radial_ring(r0, z0, pos[0]-r0, pos[1]-z0, NULL);
 			sum_ += charges[i] * jacobian_buffer[i][k] * potential;
 		}
 	}  
@@ -186,13 +186,13 @@ EXPORT double self_potential_radial(double alpha, double line_points[4][3]) {
 	double *v3 = line_points[3];
 	double *v4 = line_points[1];
 
-	double pos[2], jac;
-	position_and_jacobian_radial(alpha, v1, v2, v3, v4, pos, &jac);
+	double delta_pos[2], jac;
+	delta_position_and_jacobian_radial(alpha, v1, v2, v3, v4, delta_pos, &jac);
 
 	double target[2], jac2;
 	position_and_jacobian_radial(0, v1, v2, v3, v4, target, &jac2);
 
-	return jac*potential_radial_ring(target[0], target[1], pos[0], pos[1], NULL);
+	return jac*potential_radial_ring(target[0], target[1], delta_pos[0], delta_pos[1], NULL);
 }
 
 struct self_field_dot_normal_radial_args {
@@ -280,7 +280,7 @@ EXPORT void fill_matrix_radial(double *matrix,
 
 					double *pos = pos_buffer[j][k];
 					double jac = jacobian_buffer[j][k];
-					matrix[i*N_matrix + j] += jac * potential_radial_ring(target[0], target[1], pos[0], pos[1], NULL);
+					matrix[i*N_matrix + j] += jac * potential_radial_ring(target[0], target[1], pos[0]-target[0], pos[1]-target[1], NULL);
 				}
             }
 		}

traceon/backend/radial_ring.c

@@ -34,9 +34,10 @@ EXPORT double dr1_potential_radial_ring(double r0, double z0, double delta_r, do
     return 1./M_PI * (ellipkm1_term + ellipem1_term) / denominator;
 }
 
-EXPORT double potential_radial_ring(double r0, double z0, double r, double z, void *_) {
-    double delta_z = z - z0;
-    double delta_r = r - r0;
+EXPORT double potential_radial_ring(double r0, double z0, double delta_r, double delta_z, void *_) {
+	double r = r0 + delta_r;
+	double z = z0 + delta_z;
+
     double t = (pow(delta_z, 2) + pow(delta_r, 2)) / (pow(delta_z, 2) + pow(delta_r, 2) + 4 * r0 * delta_r + 4 * pow(r0, 2));
     return 1./M_PI * ellipkm1(t) * (delta_r + r0) / sqrt(pow(delta_z, 2) + pow((delta_r + 2 * r0), 2));
 }