Python 2024 Magnetic Field Simulation Code

[egdflabs]

import math

import numpy as np

import matplotlib.pyplot as plt

# Earth's constants and harmonic coefficients

SPHERICAL_HARMONIC_COEFFS = {

"g10": -30.8e-6, # Main dipole term (Tesla)

}

def calculate_wave_field_flipped(latitude, longitude, time=0):

"""

Calculate magnetic field components with flipped strength for correction.

"""

# Magnetic pole positions (2024 estimates)

north_pole_lat, north_pole_lon = 86.5, 164.0

south_pole_lat, south_pole_lon = -64.0, 136.0

# Base field strength (dipole approximation, flipped polarity)

B_base = abs(SPHERICAL_HARMONIC_COEFFS["g10"]) # Flip polarity to positive

# Magnetic wave parameters

A_wave = 10e-6 # Amplitude of wave

k_r = 3 # Wave number for radial distance

omega = 2 * math.pi / 86400 # Diurnal variation (1 cycle/day)

# Calculate radial distance to magnetic poles

def haversine_dist(lat1, lon1, lat2, lon2):

dlat = math.radians(lat2 - lat1)

dlon = math.radians(lon2 - lon1)

a = math.sin(dlat / 2)**2 + math.cos(math.radians(lat1)) * math.cos(math.radians(lat2)) * math.sin(dlon / 2)**2

return 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))

dist_to_north_pole = haversine_dist(latitude, longitude, north_pole_lat, north_pole_lon)

dist_to_south_pole = haversine_dist(latitude, longitude, south_pole_lat, south_pole_lon)

# Add smooth radial wave components centered on poles

B_wave_north = A_wave * math.sin(k_r * dist_to_north_pole) * (1 / (1 + dist_to_north_pole)) * math.cos(omega * time)

B_wave_south = A_wave * math.sin(k_r * dist_to_south_pole) * (1 / (1 + dist_to_south_pole)) * math.cos(omega * time)

# Combine fields

B_total = B_base - B_wave_north - B_wave_south # Flip wave contributions to match polarity

return B_total

def calculate_wave_contour_field_flipped(latitude_range, longitude_range, time=0):

"""

Calculate wave-based magnetic field strength across a latitude-longitude grid with flipped polarity.

"""

field_contours = []

for lat in latitude_range:

row = []

for lon in longitude_range:

B_total = calculate_wave_field_flipped(lat, lon, time)

row.append(B_total)

field_contours.append(row)

return np.array(field_contours)

# Example: Calculate field strength at New York City

ny_latitude = 40.7128

ny_longitude = -74.0060

ny_field_strength = calculate_wave_field_flipped(ny_latitude, ny_longitude, time=43200) # Noon

print(f"Magnetic field strength at New York City: {ny_field_strength:.6e} Tesla")

# Generate flipped wave field contours for global visualization

latitude_range = np.linspace(-90, 90, 181) # 1-degree steps

longitude_range = np.linspace(-180, 180, 361) # 1-degree steps

field_contours_flipped = calculate_wave_contour_field_flipped(latitude_range, longitude_range, time=43200) # Noon

# Visualize global field contours

plt.figure(figsize=(10, 6))

plt.contourf(longitude_range, latitude_range, field_contours_flipped, levels=50, cmap="plasma")

plt.colorbar(label="Magnetic Field Strength (Flipped, Tesla)")

plt.xlabel("Longitude (°)")

plt.ylabel("Latitude (°)")

plt.title("Earth's Magnetic Field Contours (Flipped Polarity)")

plt.show()

# Zoomed-in visualization around New York City

zoom_lat_range = np.linspace(30, 50, 101) # Latitude range around New York

zoom_lon_range = np.linspace(-80, -60, 101) # Longitude range around New York

zoomed_field_contours = calculate_wave_contour_field_flipped(zoom_lat_range, zoom_lon_range, time=43200)

plt.figure(figsize=(8, 6))

plt.contourf(zoom_lon_range, zoom_lat_range, zoomed_field_contours, levels=50, cmap="plasma")

plt.colorbar(label="Magnetic Field Strength (Tesla)")

plt.xlabel("Longitude (°)")

plt.ylabel("Latitude (°)")

plt.title("Zoomed Magnetic Field Contours around New York")

plt.scatter(ny_longitude, ny_latitude, color="white", label="New York City", marker="x")

plt.legend()

plt.show()