Refactor axial_deritvatives_radial_ring

traceon/backend/__init__.py

@@ -154,7 +154,7 @@ def __init__(self, eff, *args, **kwargs):
     'potential_radial_ring': (dbl, dbl, dbl, dbl, dbl, vp), 
     'dr1_potential_radial_ring': (dbl, dbl, dbl, dbl, dbl, vp), 
     'dz1_potential_radial_ring': (dbl, dbl, dbl, dbl, dbl, vp), 
-    'axial_derivatives_radial_ring': (None, arr(ndim=2), charges_2d, jac_buffer_2d, pos_buffer_2d, sz, z_values, sz),
+    'axial_derivatives_radial': (None, arr(ndim=2), charges_2d, jac_buffer_2d, pos_buffer_2d, sz, z_values, sz),
     'potential_radial': (dbl, v2, charges_2d, jac_buffer_2d, pos_buffer_2d, sz),
     'potential_radial_derivs': (dbl, v2, z_values, arr(ndim=3), sz),
     'flux_density_to_charge_factor': (dbl, dbl),
@@ -403,14 +403,14 @@ def trace_particle_3d_derivs(position, velocity, bounds, atol, z, electrostatic_
 dr1_potential_radial_ring = lambda *args: backend_lib.dr1_potential_radial_ring(*args, None)
 dz1_potential_radial_ring = lambda *args: backend_lib.dz1_potential_radial_ring(*args, None)
 
-def axial_derivatives_radial_ring(z, charges, jac_buffer, pos_buffer):
+def axial_derivatives_radial(z, charges, jac_buffer, pos_buffer):
     derivs = np.zeros( (z.size, DERIV_2D_MAX) )
 
     assert jac_buffer.shape == (len(charges), N_QUAD_2D)
     assert pos_buffer.shape == (len(charges), N_QUAD_2D, 2)
     assert charges.shape == (len(charges),)
 
-    backend_lib.axial_derivatives_radial_ring(derivs,charges, jac_buffer, pos_buffer, len(charges), z, len(z))
+    backend_lib.axial_derivatives_radial(derivs,charges, jac_buffer, pos_buffer, len(charges), z, len(z))
     return derivs
 
 def potential_radial(point, charges, jac_buffer, pos_buffer):

traceon/backend/radial.c

@@ -1,5 +1,4 @@
 
-#define DERIV_2D_MAX 9
 
 EXPORT const int DERIV_2D_MAX_SYM = DERIV_2D_MAX;
 
@@ -17,26 +16,22 @@ struct effective_point_charges_2d {
 
 
 EXPORT void
-axial_derivatives_radial_ring(double *derivs_p, double *charges, jacobian_buffer_2d jac_buffer, position_buffer_2d pos_buffer, size_t N_lines, double *z, size_t N_z) {
+axial_derivatives_radial(double *derivs_p, double *charges, jacobian_buffer_2d jac_buffer, position_buffer_2d pos_buffer, size_t N_lines, double *z, size_t N_z) {
 
 	double (*derivs)[DERIV_2D_MAX] = (double (*)[DERIV_2D_MAX]) derivs_p;	
 
 	for(int i = 0; i < N_z; i++) 
 	for(int j = 0; j < N_lines; j++)
 	for(int k = 0; k < N_QUAD_2D; k++) {
 		double z0 = z[i];
-		double r = pos_buffer[j][k][0], z = pos_buffer[j][k][1];
-
-		double R = norm_2d(z0-z, r);
+		double r = pos_buffer[j][k][0];
+		double z = pos_buffer[j][k][1];
+
+		double D[DERIV_2D_MAX];
+
+		axial_derivatives_radial_ring(z0, r, z, D);
 
-		double D[DERIV_2D_MAX] = {0.}; // Derivatives of the currently considered line element.
-		D[0] = 1/R;
-		D[1] = -(z0-z)/pow(R, 3);
-
-		for(int n = 1; n+1 < DERIV_2D_MAX; n++)
-			D[n+1] = -1./pow(R,2) *( (2*n + 1)*(z0-z)*D[n] + pow(n,2)*D[n-1]);
-
-		for(int l = 0; l < DERIV_2D_MAX; l++) derivs[i][l] += jac_buffer[j][k] * charges[j] * r/2 * D[l];
+		for(int l = 0; l < DERIV_2D_MAX; l++) derivs[i][l] += jac_buffer[j][k] * charges[j] * D[l];
 	}
 }
 

traceon/backend/radial_ring.c

@@ -1,5 +1,5 @@
 
-
+#define DERIV_2D_MAX 9
 
 INLINE double flux_density_to_charge_factor(double K) {
 // There is quite some derivation to this factor.
@@ -103,4 +103,18 @@ current_field_radial_ring(double x0, double y0, double x, double y, double resul
 	result[1] = 1./M_PI * 1/(2*sqrt(A)) * ( (pow(a,2) - pow(z,2) - pow(r,2))/B * ellipe(k) + ellipk(k) );
 }
 
+EXPORT void
+axial_derivatives_radial_ring(double z0, double r, double z, double derivs[DERIV_2D_MAX]) {
+
+	double R = norm_2d(z0-z, r);
+
+	derivs[0] = 1/R;
+	derivs[1] = -(z0-z)/pow(R, 3);
+
+	for(int n = 1; n+1 < DERIV_2D_MAX; n++)
+		derivs[n+1] = -1./pow(R,2) *( (2*n + 1)*(z0-z)*derivs[n] + pow(n,2)*derivs[n-1]);
+
+	for(int n = 0; n < DERIV_2D_MAX; n++)
+		derivs[n] *= r/2;
+}
 

traceon/solver.py

@@ -816,7 +816,7 @@ def get_electrostatic_axial_potential_derivatives(self, z):
         charges = self.electrostatic_point_charges.charges
         jacobians = self.electrostatic_point_charges.jacobians
         positions = self.electrostatic_point_charges.positions
-        return backend.axial_derivatives_radial_ring(z, charges, jacobians, positions)
+        return backend.axial_derivatives_radial(z, charges, jacobians, positions)
 
     def get_magnetostatic_axial_potential_derivatives(self, z):
         """
@@ -836,7 +836,7 @@ def get_magnetostatic_axial_potential_derivatives(self, z):
         jacobians = self.magnetostatic_point_charges.jacobians
         positions = self.magnetostatic_point_charges.positions
 
-        derivs_magnetic = backend.axial_derivatives_radial_ring(z, charges, jacobians, positions)
+        derivs_magnetic = backend.axial_derivatives_radial(z, charges, jacobians, positions)
         derivs_current = self.get_current_axial_potential_derivatives(z)
         return derivs_magnetic + derivs_current