Issue with prism_magnetic - unwanted asymmetry in solution

[lruepke]

Thanks so much for adding prism_magnetic!

Description of the problem: I am playing with the new functions and fail to reproduce the classic 2D solutions for a burried rod in Heitzler et al. 1962 using Harmonica's new prism_magnetic function.

Setup: Infinite rod of 10 km width and 15km vertical extent buried at 2 km depth that causes an induced field anomaly.

Expected behavior: For an inclination of 45degree (in the northern hemisphere) and an angle of 90degree between strike of the rod and geographic north, I should get a symmetric total field anomaly with a low to the North and a high to the South. For a magnetization of 0.001emu (unit in heitzler, I guess it means emu/cm^3 equal to 1 A/m), the peaks are at +145 and -145-ish gamma in Heitzler.

Harmonica result: I am using the development version of Harmonica from github. Version: v0.6.0.post33+gc677dfa I use a rod that is aligned with the easting direction. A setup of declination=0 should reproduce the strike=90 of Heitzler.

I get total field anomalies as functions of strike and inclination that are very similar to the Heitzler results. But, at inclination=45 and strike=90 (declination 0 in the harmonica setup) the predicted anomaly is not fully symmetric and the magnitude is also different at +303nT and -292 nT. For other symmetric combinations of strike and inclination, I also get asymmetries.

Possible solutions: Most likely I am doing something wrong - any help would be highly appreciated! Might also point to an issue in the kernel functions or wrappers - that's why it made it a bug report. Hope that's ok. I am unsure about the unit conversions but as far as I understand code, the 4s and pis are taken care of inside choclo and harmonica. I am also unsure about the projections of b_e, b_n, b_v to the total field but can't see anything wrong there.

Thanks! lars

Full code that generated the error

import numpy as np
import harmonica as hm # this is the development version build from source
import matplotlib.pyplot as plt

easting = np.linspace(-13e3, 13e3, 201)
northing = np.linspace(-13e3, 13e3, 201)
easting, northing = np.meshgrid(easting, northing)

# in radians
inclination = np.deg2rad(45)  # in radians
declination = np.deg2rad(0)  # in radians

# Compute induced magnetization
magnetization_strength = 1 #A/m
magnetization = [
    magnetization_strength * np.cos(inclination) * np.sin(declination),
    magnetization_strength * np.cos(inclination) * np.cos(declination),
    -magnetization_strength * np.sin(inclination) # minus because vertical axis is positive upwards in harmonica
]

prisms = np.zeros((1,6))
prisms[0, 0] = -50e3 # west - make large horizontal extent
prisms[0, 1] = 50e3  # east
prisms[0, 2] = -5e3  # south
prisms[0, 3] = 5e3  # north
prisms[0, 4] = -17e3  # bottom
prisms[0, 5] = -2e3  # top

# Compute magnetic anomaly
coordinates = [easting.ravel(), northing.ravel(), np.full_like(easting, 0).ravel()]
anomaly = hm.prism_magnetic(coordinates, prisms, magnetization)
shape = easting.shape

# get total field anomaly from b_e, b_n, b_v
anomaly_bu =  anomaly[0]*np.cos(inclination) * np.sin(declination) + \
   anomaly[1]*np.cos(inclination) * np.cos(declination) - anomaly[2]*np.sin(inclination)

# Reshape the anomaly to match the original grid shape
anomaly_reshaped = anomaly_bu.reshape(shape)

# Plot
fig, axs = plt.subplots(1, 2, figsize=(16, 6))

pc = axs[0].pcolormesh(easting, northing, anomaly_reshaped, shading='auto', cmap='RdBu_r')
fig.colorbar(pc, ax=axs[0], label='Magnetic Anomaly (nT)')
axs[0].set_xlabel('Easting (m)')
axs[0].set_ylabel('Northing (m)')
axs[0].set_title('Computed total Magnetic Anomaly')

cross_section_index = len(easting) // 2
vertical_slice = anomaly_reshaped[:, cross_section_index]  

axs[1].plot(vertical_slice, northing[:, cross_section_index])
axs[1].set_xlabel('Magnetic Anomaly (nT)')
axs[1].set_ylabel('Northing (m)')
axs[1].set_title('Max: {:.1f} nT Min: {:.1f} nT'.format(max(vertical_slice), min(vertical_slice)))

plt.show()

Full error message

System information

• Operating system: MacOS 13.6
• Python installation (Anaconda, system, ETS): anaconda,
• Version of Python: 3.11.5
• Version of this package: v0.6.0.post33+gc677dfa
• If using conda, paste the output of conda list below:

Output of conda list

PASTE OUTPUT OF CONDA LIST HERE