From 8b5d2dfcd4f46bec44d000f247d9b6e9875bb061 Mon Sep 17 00:00:00 2001 From: Karl Laundal Date: Wed, 22 Apr 2026 14:23:52 +0200 Subject: [PATCH] added function to calculate equivalent current --- lompe/model/model.py | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/lompe/model/model.py b/lompe/model/model.py index e4654c6..051a1d3 100644 --- a/lompe/model/model.py +++ b/lompe/model/model.py @@ -1097,6 +1097,41 @@ def FAC(self, lon = None, lat = None): return ju.reshape(shape) + @check_input + def equivalent_current_function(self, lon = None, lat = None): + """ + Calculate the equivalent current function psi, where j_equivalent_current = u_hat cross grad(psi) + + Requires the model vector to be defined. + + Parameters + ---------- + lon : array, optional + Longitudes [degrees] of the evaluation points, default is center of exterior grid points (grid_E). + Must have same shape as lat + lat : array, optional + Latitudes [degrees] of the evaluation points, default is center of exterior grid points (grid_E). + Must have same shape as lon + + Returns + ------- + psi : array + Equivalent current function, unit A/m + + Note + ---- + psi is calculated by evaluating the SECS G matrix for a potential, with a sign reversal + + """ + + shape = np.broadcast(lon, lat).shape + + S_df = self._B_df_matrix(return_poles = True) + G_psi = -get_SECS_J_G_matrices(lat, lon, self.lat_J, self.lon_J, current_type = 'potential', RI = self.R, singularity_limit = self.secs_singularity_limit) + + return G_psi.dot(S_df) + + @check_input def get_SECS_currents(self, lon = None, lat = None): """