Skip to content

Commit 33786b2

Browse files
authored
Minor optimizations on ellipsoids code (#610)
Allocate some square roots that are used multiple times as part of the analytic solutions of ellipsoid's gravity and magnetic forwards.
1 parent 38b6eb1 commit 33786b2

File tree

2 files changed

+10
-10
lines changed

2 files changed

+10
-10
lines changed

harmonica/_forward/ellipsoids/magnetic.py

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -394,16 +394,17 @@ def _demag_tensor_triaxial_internal(a: float, b: float, c: float):
394394
k = (a**2 - b**2) / (a**2 - c**2)
395395
ellipk = ellipkinc(phi, k)
396396
ellipe = ellipeinc(phi, k)
397+
a_c_sqrt = np.sqrt(a**2 - c**2)
397398

398-
nxx = (a * b * c) / (np.sqrt(a**2 - c**2) * (a**2 - b**2)) * (ellipk - ellipe)
399+
nxx = (a * b * c) / (a_c_sqrt * (a**2 - b**2)) * (ellipk - ellipe)
399400
nyy = (
400401
-1 * nxx
401-
+ ((a * b * c) / (np.sqrt(a**2 - c**2) * (b**2 - c**2))) * ellipe
402+
+ ((a * b * c) / (a_c_sqrt * (b**2 - c**2))) * ellipe
402403
- c**2 / (b**2 - c**2)
403404
)
404-
nzz = -1 * (
405-
(a * b * c) / (np.sqrt(a**2 - c**2) * (b**2 - c**2))
406-
) * ellipe + b**2 / (b**2 - c**2)
405+
nzz = -1 * ((a * b * c) / (a_c_sqrt * (b**2 - c**2))) * ellipe + b**2 / (
406+
b**2 - c**2
407+
)
407408

408409
return nxx, nyy, nzz
409410

@@ -426,9 +427,8 @@ def _demag_tensor_prolate_internal(a: float, b: float):
426427
if not m > 1:
427428
msg = f"Invalid aspect ratio for prolate ellipsoid: a={a}, b={b}, a/b={m}"
428429
raise ValueError(msg)
429-
nxx = (1 / (m**2 - 1)) * (
430-
((m / np.sqrt(m**2 - 1)) * np.log(m + np.sqrt(m**2 - 1))) - 1
431-
)
430+
sqrt = np.sqrt(m**2 - 1)
431+
nxx = (1 / (m**2 - 1)) * (((m / sqrt) * np.log(m + sqrt)) - 1)
432432
nyy = nzz = 0.5 * (1 - nxx)
433433
return nxx, nyy, nzz
434434

harmonica/_forward/ellipsoids/utils.py

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -347,8 +347,8 @@ def _get_elliptical_integrals_prolate(a, b, lmbda):
347347
sqrt_l2 = np.sqrt(b**2 + lmbda)
348348
log = np.log((sqrt_e + sqrt_l1) / sqrt_l2)
349349

350-
g1 = (2 / (e2 ** (3 / 2))) * (log - sqrt_e / sqrt_l1)
351-
g2 = (1 / (e2 ** (3 / 2))) * ((sqrt_e * sqrt_l1) / (b**2 + lmbda) - log)
350+
g1 = (2 / (sqrt_e**3)) * (log - sqrt_e / sqrt_l1)
351+
g2 = (1 / (sqrt_e**3)) * ((sqrt_e * sqrt_l1) / (b**2 + lmbda) - log)
352352
return g1, g2, g2
353353

354354

0 commit comments

Comments
 (0)