Computational Programs

10 interactive simulations covering solar wind dynamics, magnetospheric structure, ionospheric physics, and space weather modeling. Programs use Python and Fortran with real-time execution.

1. CME Drag-Based Propagation

Model CME propagation using the aerodynamic drag equation (Vršnak et al. 2013).

CME Propagation

Python

Parameters

Base CME Speed (km/s)

Solar Wind Speed (km/s)

Drag Parameter (km⁻¹)

script.py70 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

V_CME_0 = 1500
V_CME = V_CME_0
V_SW_0 = 400
V_SW = V_SW_0
GAMMA_0 = 5e-8
GAMMA_DRAG = GAMMA_0

AU_km = 1.496e8

def drag_model(t, y):
    v, r = y
    dv = -GAMMA_DRAG * (v - V_SW) * abs(v - V_SW)
    dr = v
    return [dv, dr]

# Solve for multiple initial speeds
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
fig.patch.set_facecolor('#0f172a')
for ax in axes:
    ax.set_facecolor('#1e293b')
    ax.tick_params(colors='#cbd5e1')
    for spine in ax.spines.values():
        spine.set_color('#334155')

colors = ['#f87171', '#60a5fa', '#34d399', '#fbbf24', '#a78bfa']
speeds = [V_CME*0.5, V_CME*0.75, V_CME, V_CME*1.25, V_CME*1.5]

for i, v0 in enumerate(speeds):
    t_span = (0, 5*86400)
    t_eval = np.linspace(*t_span, 3000)
    sol = solve_ivp(drag_model, t_span, [v0, 0], t_eval=t_eval, max_step=600)
    t_hr = sol.t / 3600
    v = sol.y[0]
    r_au = sol.y[1] / AU_km

axes[0].plot(t_hr, v, color=colors[i], linewidth=1.5, label=f'{v0:.0f} km/s')
    axes[1].plot(t_hr, r_au, color=colors[i], linewidth=1.5)

idx_1au = np.argmin(np.abs(r_au - 1.0))
    if r_au[idx_1au] > 0.9:
        axes[2].bar(i, t_hr[idx_1au], color=colors[i], alpha=0.8)
        axes[2].text(i, t_hr[idx_1au]+1, f'{t_hr[idx_1au]:.1f}h', ha='center', color=colors[i], fontsize=9)

axes[0].axhline(V_SW, color='white', linestyle='--', linewidth=0.8, alpha=0.5, label=f'v_sw={V_SW} km/s')
axes[0].set_xlabel('Time (hours)', color='#cbd5e1')
axes[0].set_ylabel('Speed (km/s)', color='#cbd5e1')
axes[0].set_title('CME Speed Evolution', color='#f87171', fontweight='bold')
axes[0].legend(fontsize=7, facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')

axes[1].axhline(1.0, color='white', linestyle='--', linewidth=0.8, alpha=0.5)
axes[1].set_xlabel('Time (hours)', color='#cbd5e1')
axes[1].set_ylabel('Distance (AU)', color='#cbd5e1')
axes[1].set_title('Heliocentric Distance', color='#60a5fa', fontweight='bold')

axes[2].set_xlabel('CME Index', color='#cbd5e1')
axes[2].set_ylabel('Transit Time (hours)', color='#cbd5e1')
axes[2].set_title('1 AU Arrival Time', color='#34d399', fontweight='bold')

fig.suptitle(f'CME Drag Model: γ = {GAMMA_DRAG:.1e} km⁻¹, v_sw = {V_SW} km/s', color='white', fontsize=13, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a')
plt.close()
print(f"CME drag propagation model computed for {len(speeds)} initial speeds")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

2. Parker Solar Wind Solutions

Solve Parker's isothermal wind equation at multiple coronal temperatures.

Parker Wind

Python

Parameters

Base Coronal Temperature (K)

script.py67 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy.optimize import brentq

T_CORONA = 1.5e6
T_VAL = T_CORONA

k_B = 1.381e-23; m_p = 1.673e-27; G = 6.674e-11
M_sun = 1.989e30; R_sun = 6.957e8

temperatures = [T_VAL * 0.5, T_VAL * 0.75, T_VAL, T_VAL * 1.25, T_VAL * 1.5]

fig, axes = plt.subplots(1, 2, figsize=(14, 6))
fig.patch.set_facecolor('#0f172a')
for ax in axes:
    ax.set_facecolor('#1e293b')
    ax.tick_params(colors='#cbd5e1')
    for spine in ax.spines.values():
        spine.set_color('#334155')

colors = ['#60a5fa', '#34d399', '#f87171', '#fbbf24', '#a78bfa']

for idx, T in enumerate(temperatures):
    c_s = np.sqrt(k_B * T / (0.5 * m_p))
    r_c = G * M_sun / (2 * c_s**2)

def parker_eq(u, r):
        xi = r / r_c
        return u**2 - np.log(u**2) - 4*np.log(xi) - 4/xi + 3

r = np.linspace(1.05 * R_sun, 250 * R_sun, 3000)
    v_wind = np.zeros_like(r)
    for i, ri in enumerate(r):
        try:
            if ri < r_c:
                v_wind[i] = brentq(parker_eq, 1e-4, 0.999, args=(ri,)) * c_s
            else:
                v_wind[i] = brentq(parker_eq, 1.001, 30.0, args=(ri,)) * c_s
        except:
            v_wind[i] = np.nan

r_Rs = r / R_sun
    axes[0].plot(r_Rs, v_wind/1e3, color=colors[idx], linewidth=1.5, label=f'T = {T/1e6:.2f} MK')
    axes[1].plot(r_Rs, v_wind/1e3, color=colors[idx], linewidth=1.5)
    axes[1].axvline(r_c/R_sun, color=colors[idx], linestyle=':', linewidth=0.8, alpha=0.5)

axes[0].set_xlabel('r / R_sun', color='#cbd5e1')
axes[0].set_ylabel('Speed (km/s)', color='#cbd5e1')
axes[0].set_title('Parker Wind: Near Sun', color='#f87171', fontweight='bold')
axes[0].set_xlim(1, 30)
axes[0].legend(fontsize=8, facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')

axes[1].set_xlabel('r / R_sun', color='#cbd5e1')
axes[1].set_ylabel('Speed (km/s)', color='#cbd5e1')
axes[1].set_title('Parker Wind: Full Range', color='#60a5fa', fontweight='bold')
axes[1].set_xlim(1, 250)
axes[1].axvline(215, color='white', linestyle='--', linewidth=0.8, alpha=0.5, label='~1 AU')
axes[1].legend(fontsize=8, facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')

fig.suptitle('Parker Solar Wind Solutions at Different Coronal Temperatures', color='white', fontsize=13, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a')
plt.close()
print(f"Parker wind computed for {len(temperatures)} temperatures")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

3. Magnetopause & Bow Shock Shape (Fortran)

Compute magnetopause and bow shock positions using the Shue et al. (1998) model, compiled and run in Fortran.

Magnetopause Shape

Python

Parameters

Dynamic Pressure (nPa)

IMF Bz (nT)

script.py117 lines

import subprocess, os, tempfile, sys
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

DP_VAL = 2.0
Dp = DP_VAL
BZ_VAL = -5.0
Bz = BZ_VAL

fortran_code = f"""
program magnetopause_shape
  implicit none
  integer, parameter :: N = 360
  real(8), parameter :: pi = 3.141592653589793d0
  real(8) :: Dp, Bz, r0, alpha_val, theta, r_mp, r_bs, delta_R
  integer :: i

Dp = {Dp}d0
  Bz = {Bz}d0

! Shue et al. (1998) model
  r0 = (10.22d0 + 1.29d0*tanh(0.184d0*(Bz + 8.14d0))) * Dp**(-1.0d0/6.6d0)
  alpha_val = (0.58d0 - 0.007d0*Bz) * (1.0d0 + 0.024d0*log(Dp))

! Bow shock standoff (Fairfield 1971 / Spreiter)
  delta_R = 1.1d0 * r0 * 0.28d0

write(*,'(A)') '# theta(deg) r_mp(Re) x_mp(Re) y_mp(Re) r_bs(Re) x_bs(Re) y_bs(Re)'
  write(*,'(A,F8.3)') '# r0 = ', r0
  write(*,'(A,F8.3)') '# alpha = ', alpha_val

do i = 0, N
    theta = -pi + 2.0d0*pi*dble(i)/dble(N)
    r_mp = r0 * (2.0d0 / (1.0d0 + cos(theta)))**alpha_val
    r_bs = r_mp + delta_R * (2.0d0 / (1.0d0 + cos(theta)))**0.6d0
    write(*,'(7F12.4)') theta*180.0d0/pi, r_mp, r_mp*cos(theta), r_mp*sin(theta), &
                         r_bs, r_bs*cos(theta), r_bs*sin(theta)
  end do
end program
"""

with tempfile.TemporaryDirectory() as tmpdir:
    src = os.path.join(tmpdir, 'mpause.f90')
    exe = os.path.join(tmpdir, 'mpause')
    with open(src, 'w') as f:
        f.write(fortran_code)

comp = subprocess.run(['gfortran', '-O2', '-o', exe, src], capture_output=True, text=True)
    if comp.returncode != 0:
        print("Compilation error:", comp.stderr)
        sys.exit(1)

result = subprocess.run([exe], capture_output=True, text=True)
    lines = result.stdout.strip().split('\n')

header_lines = [l for l in lines if l.startswith('#')]
data_lines = [l for l in lines if not l.startswith('#')]

for hl in header_lines:
    print(hl)

data = np.array([[float(x) for x in line.split()] for line in data_lines if line.strip()])

theta_deg = data[:,0]
x_mp = data[:,2]
y_mp = data[:,3]
x_bs = data[:,5]
y_bs = data[:,6]

# Left: Full 2D shape
axes[0].plot(x_mp, y_mp, color='#f87171', linewidth=2, label='Magnetopause (Shue)')
axes[0].plot(x_bs, y_bs, color='#60a5fa', linewidth=2, linestyle='--', label='Bow Shock')
axes[0].plot(0, 0, 'o', color='#34d399', markersize=12, label='Earth')
axes[0].set_xlabel('X_GSM (R_E)', color='#cbd5e1')
axes[0].set_ylabel('Y_GSM (R_E)', color='#cbd5e1')
axes[0].set_title(f'Magnetopause & Bow Shock (Dp={Dp} nPa, Bz={Bz} nT)', color='#f87171', fontweight='bold', fontsize=10)
axes[0].legend(fontsize=8, facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')
axes[0].set_aspect('equal')
axes[0].set_xlim(-40, 20)
axes[0].grid(True, alpha=0.15, color='#475569')
axes[0].axhline(0, color='#475569', linewidth=0.5)
axes[0].axvline(0, color='#475569', linewidth=0.5)

# Right: Parametric study
Bz_vals = [-10, -5, 0, 5, 10]
colors_bz = ['#ef4444', '#f87171', '#fbbf24', '#34d399', '#60a5fa']
for bz, col in zip(Bz_vals, colors_bz):
    r0_t = (10.22 + 1.29*np.tanh(0.184*(bz + 8.14))) * Dp**(-1/6.6)
    alpha_t = (0.58 - 0.007*bz) * (1 + 0.024*np.log(Dp))
    th = np.linspace(-np.pi, np.pi, 300)
    r_t = r0_t * (2 / (1 + np.cos(th)))**alpha_t
    axes[1].plot(r_t*np.cos(th), r_t*np.sin(th), color=col, linewidth=1.5, label=f'Bz={bz} nT')

axes[1].plot(0, 0, 'o', color='#34d399', markersize=10)
axes[1].set_xlabel('X_GSM (R_E)', color='#cbd5e1')
axes[1].set_ylabel('Y_GSM (R_E)', color='#cbd5e1')
axes[1].set_title(f'Magnetopause vs IMF Bz (Dp={Dp} nPa)', color='#a78bfa', fontweight='bold')
axes[1].legend(fontsize=8, facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')
axes[1].set_aspect('equal')
axes[1].set_xlim(-40, 20)
axes[1].grid(True, alpha=0.15, color='#475569')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a')
plt.close()
print(f"Magnetopause model computed for Dp={Dp} nPa, Bz={Bz} nT")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

4. Magnetospheric Convection Patterns

Volland-Stern convection electric field + corotation potential, showing equipotential drift paths.

Convection

Python

Parameters

Cross-Polar Cap Potential (kV)

script.py85 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

CPCP_VAL = 60.0
CPCP = CPCP_VAL

# Volland-Stern convection + corotation model
# Phi_conv = A * L^gamma * sin(phi)  +  Phi_corot
# A chosen so that CPCP matches desired value

R_E = 6.371e6
B0 = 3.12e-5
Omega_E = 7.27e-5
gamma_VS = 2  # Volland-Stern exponent

# Grid in L-shell and MLT
L = np.linspace(1.5, 10, 200)
phi = np.linspace(0, 2*np.pi, 200)
L2d, phi2d = np.meshgrid(L, phi)

# Corotation potential (kV)
Phi_cor = -Omega_E * B0 * R_E**2 / L2d * 1e-3  # in kV

# Convection potential
A = CPCP / (2 * 10**gamma_VS)  # kV / L^gamma
Phi_conv = A * L2d**gamma_VS * np.sin(phi2d)

Phi_total = Phi_conv + Phi_cor

# Convert to Cartesian for plotting
x = L2d * np.cos(phi2d)  # noon at top
y = L2d * np.sin(phi2d)

fig, axes = plt.subplots(1, 3, figsize=(18, 6))
fig.patch.set_facecolor('#0f172a')
for ax in axes:
    ax.set_facecolor('#1e293b')
    ax.tick_params(colors='#cbd5e1')
    for spine in ax.spines.values():
        spine.set_color('#334155')

# Panel 1: Corotation
c1 = axes[0].contour(y, x, Phi_cor, levels=15, colors='#60a5fa', linewidths=0.8)
axes[0].set_title('Corotation Potential (kV)', color='#60a5fa', fontweight='bold')
axes[0].set_aspect('equal')
circle = plt.Circle((0,0), 1, color='#34d399', fill=True, zorder=5)
axes[0].add_patch(circle)
axes[0].set_xlabel('Dawn-Dusk (R_E)', color='#cbd5e1')
axes[0].set_ylabel('Noon-Midnight (R_E)', color='#cbd5e1')

# Panel 2: Convection
c2 = axes[1].contour(y, x, Phi_conv, levels=15, colors='#f87171', linewidths=0.8)
axes[1].set_title(f'Convection Potential (CPCP={CPCP:.0f} kV)', color='#f87171', fontweight='bold')
axes[1].set_aspect('equal')
circle2 = plt.Circle((0,0), 1, color='#34d399', fill=True, zorder=5)
axes[1].add_patch(circle2)
axes[1].set_xlabel('Dawn-Dusk (R_E)', color='#cbd5e1')

# Panel 3: Total with plasmapause
c3 = axes[2].contour(y, x, Phi_total, levels=20, cmap='RdBu_r', linewidths=0.8)
axes[2].set_title('Total Potential (Convection + Corotation)', color='#fbbf24', fontweight='bold')
axes[2].set_aspect('equal')
circle3 = plt.Circle((0,0), 1, color='#34d399', fill=True, zorder=5)
axes[2].add_patch(circle3)
axes[2].set_xlabel('Dawn-Dusk (R_E)', color='#cbd5e1')

# Find approximate plasmapause (last closed equipotential)
Phi_pp = -Omega_E * B0 * R_E**2 * 1e-3 * A**(1.0/gamma_VS) if A > 0 else 0
axes[2].contour(y, x, Phi_total, levels=[Phi_pp*0.8], colors='#fbbf24', linewidths=2, linestyles='--')

for ax in axes:
    ax.set_xlim(-10, 10)
    ax.set_ylim(-10, 10)
    ax.text(0, 10.5, 'Noon', ha='center', color='#cbd5e1', fontsize=8)
    ax.text(0, -10.5, 'Midnight', ha='center', color='#cbd5e1', fontsize=8)
    ax.text(10.5, 0, 'Dusk', ha='center', color='#cbd5e1', fontsize=8, rotation=90)
    ax.text(-10.5, 0, 'Dawn', ha='center', color='#cbd5e1', fontsize=8, rotation=90)

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a')
plt.close()
print(f"Convection pattern computed for CPCP = {CPCP} kV")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

5. Harris Current Sheet Structure

Compute magnetic field, current density, and pressure profiles for a Harris equilibrium.

Harris Sheet

Python

Parameters

Lobe Field B₀ (nT)

Half-thickness δ (R_E)

script.py58 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

B0_VAL = 20.0
B0 = B0_VAL  # nT
DELTA_VAL = 2.0
delta = DELTA_VAL  # R_E, half-thickness

mu0 = 4 * np.pi * 1e-7
R_E = 6.371e6

z = np.linspace(-5*delta, 5*delta, 500)

Bx = B0 * np.tanh(z / delta)
Jy = (B0 * 1e-9) / (mu0 * delta * R_E) * 1.0 / np.cosh(z / delta)**2  # A/m^2
Jy_nA = Jy * 1e9  # nA/m^2
P = (B0*1e-9)**2 / (2*mu0) * 1.0 / np.cosh(z / delta)**2  # Pa
P_nPa = P * 1e9

fig, axes = plt.subplots(1, 3, figsize=(16, 5))
fig.patch.set_facecolor('#0f172a')
for ax in axes:
    ax.set_facecolor('#1e293b')
    ax.tick_params(colors='#cbd5e1')
    for spine in ax.spines.values():
        spine.set_color('#334155')

axes[0].plot(z, Bx, color='#f87171', linewidth=2)
axes[0].axhline(0, color='#475569', linewidth=0.5)
axes[0].axhline(B0, color='#f87171', linewidth=0.8, linestyle='--', alpha=0.5)
axes[0].axhline(-B0, color='#f87171', linewidth=0.8, linestyle='--', alpha=0.5)
axes[0].fill_between(z, Bx, alpha=0.1, color='#f87171')
axes[0].set_xlabel('z (R_E)', color='#cbd5e1')
axes[0].set_ylabel('Bx (nT)', color='#cbd5e1')
axes[0].set_title('Magnetic Field B_x(z)', color='#f87171', fontweight='bold')

axes[1].plot(z, Jy_nA, color='#60a5fa', linewidth=2)
axes[1].fill_between(z, Jy_nA, alpha=0.15, color='#60a5fa')
axes[1].set_xlabel('z (R_E)', color='#cbd5e1')
axes[1].set_ylabel('J_y (nA/m²)', color='#cbd5e1')
axes[1].set_title('Current Density j_y(z)', color='#60a5fa', fontweight='bold')

axes[2].plot(z, P_nPa, color='#34d399', linewidth=2)
axes[2].fill_between(z, P_nPa, alpha=0.15, color='#34d399')
axes[2].axhline((B0*1e-9)**2/(2*mu0)*1e9, color='#fbbf24', linestyle='--', linewidth=1, label='B²/2μ₀ (lobe)')
axes[2].set_xlabel('z (R_E)', color='#cbd5e1')
axes[2].set_ylabel('Pressure (nPa)', color='#cbd5e1')
axes[2].set_title('Plasma Pressure P(z)', color='#34d399', fontweight='bold')
axes[2].legend(fontsize=8, facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')

fig.suptitle(f'Harris Current Sheet: B₀ = {B0} nT, δ = {delta} R_E', color='white', fontsize=13, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a')
plt.close()
print(f"Harris sheet: B0={B0} nT, delta={delta} R_E, peak J={np.max(Jy_nA):.2f} nA/m²")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

6. Dst Storm Evolution (Fortran)

Simulate geomagnetic storm Dst index using the O'Brien-McPherron model in Fortran.

Dst Storm

Python

Parameters

Solar Wind Speed (km/s)

IMF Bz (nT)

Southward Duration (hours)

script.py112 lines

import subprocess, os, tempfile, sys
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

VSW_VAL = 500.0
v_sw = VSW_VAL
BZ_STORM = -15.0
Bz = BZ_STORM
DURATION_HR = 6.0
dur = DURATION_HR

fortran_code = f"""
program dst_evolution
  implicit none
  integer, parameter :: N = 7200
  real(8) :: dt, t, Dst, Ey, Q, tau_val
  real(8) :: v_sw, Bz, dur_hr
  integer :: i

v_sw = {v_sw}d0
  Bz = {Bz}d0
  dur_hr = {dur}d0

dt = 0.01d0  ! hours
  Dst = 0.0d0

write(*,'(A)') '# t(hr) Dst(nT) Ey(mV/m) Q(nT/hr) tau(hr)'

do i = 0, N
    t = dble(i) * dt

! IMF Bz profile: southward for dur hours, then northward
    if (t < dur_hr) then
      Ey = v_sw * abs(Bz) * 1.0d-3  ! mV/m
    else
      Ey = 0.0d0
    end if

! O'Brien-McPherron (2000) model
    if (Ey > 0.5d0) then
      Q = -4.4d0 * (Ey - 0.5d0)
    else
      Q = 0.0d0
    end if

tau_val = 2.40d0 * exp(9.74d0 / (4.69d0 + Ey))

! Euler step
    Dst = Dst + (Q - Dst/tau_val) * dt

if (mod(i, 10) == 0) then
      write(*,'(5F12.4)') t, Dst, Ey, Q, tau_val
    end if
  end do
end program
"""

with tempfile.TemporaryDirectory() as tmpdir:
    src = os.path.join(tmpdir, 'dst.f90')
    exe = os.path.join(tmpdir, 'dst')
    with open(src, 'w') as f:
        f.write(fortran_code)
    comp = subprocess.run(['gfortran', '-O2', '-o', exe, src], capture_output=True, text=True)
    if comp.returncode != 0:
        print("Compilation error:", comp.stderr)
        sys.exit(1)
    result = subprocess.run([exe], capture_output=True, text=True)
    lines = result.stdout.strip().split('\n')

data_lines = [l for l in lines if not l.startswith('#')]
data = np.array([[float(x) for x in line.split()] for line in data_lines if line.strip()])

t = data[:,0]; Dst_arr = data[:,1]; Ey = data[:,2]; Q = data[:,3]; tau = data[:,4]

fig, axes = plt.subplots(2, 2, figsize=(14, 10))
fig.patch.set_facecolor('#0f172a')
for row in axes:
    for ax in row:
        ax.set_facecolor('#1e293b')
        ax.tick_params(colors='#cbd5e1')
        for spine in ax.spines.values():
            spine.set_color('#334155')

axes[0,0].plot(t, Ey, color='#fbbf24', linewidth=2)
axes[0,0].set_ylabel('Ey (mV/m)', color='#cbd5e1')
axes[0,0].set_title('Dawn-Dusk E-field (vBs)', color='#fbbf24', fontweight='bold')

axes[0,1].plot(t, Q, color='#f87171', linewidth=2)
axes[0,1].set_ylabel('Q (nT/hr)', color='#cbd5e1')
axes[0,1].set_title('Injection Rate', color='#f87171', fontweight='bold')

axes[1,0].plot(t, Dst_arr, color='#60a5fa', linewidth=2)
axes[1,0].axhline(-50, color='#fbbf24', linestyle='--', linewidth=0.8, alpha=0.5, label='Moderate storm')
axes[1,0].axhline(-100, color='#f87171', linestyle='--', linewidth=0.8, alpha=0.5, label='Intense storm')
axes[1,0].set_xlabel('Time (hours)', color='#cbd5e1')
axes[1,0].set_ylabel('Dst (nT)', color='#cbd5e1')
axes[1,0].set_title(f'Dst Index (min = {np.min(Dst_arr):.0f} nT)', color='#60a5fa', fontweight='bold')
axes[1,0].legend(fontsize=8, facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')

axes[1,1].plot(t, tau, color='#34d399', linewidth=2)
axes[1,1].set_xlabel('Time (hours)', color='#cbd5e1')
axes[1,1].set_ylabel('τ (hours)', color='#cbd5e1')
axes[1,1].set_title('Ring Current Decay Time', color='#34d399', fontweight='bold')

fig.suptitle(f'Geomagnetic Storm: v_sw={v_sw:.0f} km/s, Bz={Bz:.0f} nT for {dur:.0f} hr', color='white', fontsize=13, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a')
plt.close()
print(f"Storm simulation: min Dst = {np.min(Dst_arr):.0f} nT")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

7. Trapped Particle Trajectories in Dipole

Compute mirror points, bounce/drift periods, and loss cone for protons trapped in a dipole field.

Trapped Particles

Python

Parameters

Proton Energy (keV)

L-shell

Equatorial Pitch Angle (°)

script.py130 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

ENERGY_KEV = 100.0
E_keV = ENERGY_KEV
L_SHELL = 5.0
L = L_SHELL
PITCH_DEG = 45.0
alpha_eq = PITCH_DEG

# Constants
q = 1.602e-19; m_p = 1.673e-27
R_E = 6.371e6; B0 = 3.12e-5  # dipole moment

E_J = E_keV * 1e3 * q
v = np.sqrt(2*E_J/m_p)

# Dipole field
def B_dipole(r, theta):
    Br = -2*B0*(R_E/r)**3 * np.cos(theta)
    Bt = -B0*(R_E/r)**3 * np.sin(theta)
    return np.sqrt(Br**2 + Bt**2), Br, Bt

B_eq = B0 / L**3
v_perp = v * np.sin(np.radians(alpha_eq))
v_par = v * np.cos(np.radians(alpha_eq))

# Compute mirror point
mu = 0.5 * m_p * v_perp**2 / B_eq
B_mirror = 0.5 * m_p * v**2 / mu

# Find mirror latitude
from scipy.optimize import brentq
def mirror_eq(lam):
    r = L * R_E * np.cos(lam)**2
    B_mag = B0 * np.sqrt(1 + 3*np.sin(lam)**2) / (L**3 * np.cos(lam)**6)
    return B_mag - B_mirror

try:
    lam_m = brentq(mirror_eq, 0.001, np.pi/2 - 0.01)
    lam_m_deg = np.degrees(lam_m)
except:
    lam_m_deg = 90
    lam_m = np.pi/2

# Loss cone
B_atm = B0 * np.sqrt(1 + 3) / 1.0  # at surface, L=1
B_loss = B0 * np.sqrt(4) / 1  # simplified
alpha_LC = np.degrees(np.arcsin(np.sqrt(B_eq / (B0 * L**3 * np.sqrt(4)))))

# Bounce period approximation
tau_b = 4 * L * R_E / v * (1.30 - 0.56 * np.sin(np.radians(alpha_eq)))

# Drift period (gradient-curvature)
tau_d = 2 * np.pi * q * B0 * R_E**2 / (3 * E_J * L)

# Plot field line and particle info
fig, axes = plt.subplots(1, 3, figsize=(18, 6))
fig.patch.set_facecolor('#0f172a')
for ax in axes:
    ax.set_facecolor('#1e293b')
    ax.tick_params(colors='#cbd5e1')
    for spine in ax.spines.values():
        spine.set_color('#334155')

# Panel 1: Dipole field lines and mirror points
theta_fl = np.linspace(-np.pi/2 + 0.05, np.pi/2 - 0.05, 500)
for Li in [2, 3, 4, 5, 6, 7, 8]:
    r_fl = Li * np.cos(theta_fl)**2
    x_fl = r_fl * np.cos(theta_fl)
    y_fl = r_fl * np.sin(theta_fl)
    col = '#f87171' if Li == L else '#475569'
    lw = 2 if Li == L else 0.5
    axes[0].plot(x_fl, y_fl, color=col, linewidth=lw)

# Mirror points
r_m = L * np.cos(lam_m)**2
axes[0].plot(r_m*np.cos(lam_m), r_m*np.sin(lam_m), 'o', color='#fbbf24', markersize=10, zorder=5, label=f'Mirror λ={lam_m_deg:.1f}°')
axes[0].plot(r_m*np.cos(lam_m), -r_m*np.sin(lam_m), 'o', color='#fbbf24', markersize=10, zorder=5)

earth = plt.Circle((0,0), 1, color='#34d399', fill=True, zorder=10)
axes[0].add_patch(earth)
axes[0].set_aspect('equal')
axes[0].set_xlim(-1, 10)
axes[0].set_ylim(-6, 6)
axes[0].set_xlabel('x (R_E)', color='#cbd5e1')
axes[0].set_ylabel('z (R_E)', color='#cbd5e1')
axes[0].set_title(f'Dipole Field Lines (L={L})', color='#f87171', fontweight='bold')
axes[0].legend(fontsize=8, facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')

# Panel 2: B vs latitude along field line
lam_arr = np.linspace(0, lam_m*1.2, 200)
B_arr = B0 * np.sqrt(1 + 3*np.sin(lam_arr)**2) / (L**3 * np.cos(lam_arr)**6)
axes[1].plot(np.degrees(lam_arr), B_arr*1e9, color='#a78bfa', linewidth=2)
axes[1].axhline(B_mirror*1e9, color='#fbbf24', linestyle='--', linewidth=1, label=f'B_mirror = {B_mirror*1e9:.1f} nT')
axes[1].axhline(B_eq*1e9, color='#34d399', linestyle='--', linewidth=1, label=f'B_eq = {B_eq*1e9:.1f} nT')
axes[1].set_xlabel('Latitude (°)', color='#cbd5e1')
axes[1].set_ylabel('|B| (nT)', color='#cbd5e1')
axes[1].set_title('Field Strength vs Latitude', color='#a78bfa', fontweight='bold')
axes[1].legend(fontsize=8, facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')

# Panel 3: Summary text
info = [
    f"Proton Energy: {E_keV} keV",
    f"L-shell: {L}",
    f"Eq. Pitch Angle: {alpha_eq}°",
    f"Speed: {v/1e6:.2f} × 10⁶ m/s",
    f"",
    f"Mirror Latitude: {lam_m_deg:.1f}°",
    f"Loss Cone: {alpha_LC:.1f}°",
    f"Trapped: {'Yes' if alpha_eq > alpha_LC else 'NO - in loss cone!'}",
    f"",
    f"Bounce Period: {tau_b:.2f} s",
    f"Drift Period: {tau_d:.0f} s ({tau_d/60:.1f} min)",
    f"Gyroradius: {m_p*v_perp/(q*B_eq):.0f} m",
]
for i, line in enumerate(info):
    col = '#f87171' if 'NO' in line else '#cbd5e1'
    axes[2].text(0.1, 0.95 - i*0.07, line, transform=axes[2].transAxes, fontsize=11, color=col, family='monospace')
axes[2].set_title('Particle Parameters', color='#34d399', fontweight='bold')
axes[2].axis('off')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a')
plt.close()
print(f"Trapped {E_keV} keV proton at L={L}: mirror at {lam_m_deg:.1f}°, bounce={tau_b:.2f} s, drift={tau_d:.0f} s")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

8. Ionospheric Conductivity Profiles

Compute Pedersen, Hall, and parallel conductivities as functions of altitude.

Conductivity

Python

Parameters

Solar Zenith Angle (°)

F10.7 Solar Flux (sfu)

script.py93 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

CHI_DEG = 30.0
chi = CHI_DEG
F107_VAL = 150.0
F107 = F107_VAL

# Altitude grid
z = np.linspace(60, 500, 500)  # km

# Neutral density (exponential atmosphere)
H_n = 7.0  # neutral scale height (km), simplified
n_neutral = 2e19 * np.exp(-(z - 100)/H_n)  # m^-3

# Electron density (Chapman layer approximation for E and F regions)
def chapman(z_arr, z0, H, NmF2):
    y = (z_arr - z0) / H
    return NmF2 * np.exp(0.5*(1 - y - 1/np.cos(np.radians(chi))*np.exp(-y)))

ne_E = chapman(z, 110, 10, 1e11)  # E region
ne_F1 = chapman(z, 170, 20, 5e10 * (F107/150))
ne_F2 = chapman(z, 300, 60, 3e11 * (F107/150))
ne = ne_E + ne_F1 + ne_F2

# Ion-neutral and electron-neutral collision frequencies
nu_in = 2.6e-15 * n_neutral * 1e-6  # simplified, Hz
nu_en = 5.4e-16 * n_neutral * 1e-6 * np.sqrt(300)

# Gyrofrequencies
B = 5e-5  # T (50000 nT)
q = 1.602e-19; m_e = 9.109e-31; m_i = 16 * 1.673e-27  # O+
Omega_e = q * B / m_e
Omega_i = q * B / m_i

# Conductivities
sigma_P = ne * q / B * (nu_in**2/(nu_in**2 + Omega_i**2) + nu_en**2/(nu_en**2 + Omega_e**2))
sigma_H = ne * q / B * (Omega_e*nu_en/(nu_en**2 + Omega_e**2) - Omega_i*nu_in/(nu_in**2 + Omega_i**2))
sigma_0 = ne * q**2 / (m_e * nu_en)  # parallel

fig, axes = plt.subplots(1, 3, figsize=(16, 7))
fig.patch.set_facecolor('#0f172a')
for ax in axes:
    ax.set_facecolor('#1e293b')
    ax.tick_params(colors='#cbd5e1')
    for spine in ax.spines.values():
        spine.set_color('#334155')

# Panel 1: Electron density
axes[0].semilogx(ne/1e6, z, color='#60a5fa', linewidth=2)
axes[0].semilogx(ne_E/1e6, z, color='#34d399', linewidth=1, linestyle='--', label='E region')
axes[0].semilogx(ne_F2/1e6, z, color='#f87171', linewidth=1, linestyle='--', label='F2 region')
axes[0].set_xlabel('Electron Density (cm⁻³)', color='#cbd5e1')
axes[0].set_ylabel('Altitude (km)', color='#cbd5e1')
axes[0].set_title(f'Electron Density (χ={chi}°, F10.7={F107})', color='#60a5fa', fontweight='bold')
axes[0].legend(fontsize=8, facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')
axes[0].set_ylim(60, 500)

# Panel 2: Collision vs gyro frequencies
axes[1].semilogx(nu_in, z, color='#f87171', linewidth=2, label='ν_in')
axes[1].semilogx(nu_en, z, color='#60a5fa', linewidth=2, label='ν_en')
axes[1].axvline(Omega_i, color='#fbbf24', linestyle='--', linewidth=1, label=f'Ω_i = {Omega_i:.0f} Hz')
axes[1].axvline(Omega_e, color='#a78bfa', linestyle='--', linewidth=1, label=f'Ω_e = {Omega_e:.0e} Hz')
axes[1].set_xlabel('Frequency (Hz)', color='#cbd5e1')
axes[1].set_ylabel('Altitude (km)', color='#cbd5e1')
axes[1].set_title('Collision vs Gyrofrequencies', color='#fbbf24', fontweight='bold')
axes[1].legend(fontsize=7, facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')
axes[1].set_ylim(60, 500)

# Panel 3: Conductivities
axes[2].semilogx(sigma_P, z, color='#f87171', linewidth=2, label='Pedersen σ_P')
axes[2].semilogx(np.abs(sigma_H), z, color='#60a5fa', linewidth=2, label='Hall |σ_H|')
axes[2].semilogx(sigma_0, z, color='#34d399', linewidth=2, linestyle='--', label='Parallel σ₀')
axes[2].set_xlabel('Conductivity (S/m)', color='#cbd5e1')
axes[2].set_ylabel('Altitude (km)', color='#cbd5e1')
axes[2].set_title('Ionospheric Conductivities', color='#f87171', fontweight='bold')
axes[2].legend(fontsize=8, facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')
axes[2].set_ylim(60, 300)

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a')
plt.close()

# Height-integrated conductances
dz_m = (z[1]-z[0]) * 1e3
Sigma_P = np.sum(sigma_P) * dz_m
Sigma_H = np.sum(np.abs(sigma_H)) * dz_m
print(f"Height-integrated Pedersen conductance: {Sigma_P:.2f} S")
print(f"Height-integrated Hall conductance: {Sigma_H:.2f} S")
print(f"Hall/Pedersen ratio: {Sigma_H/Sigma_P:.2f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

9. Auroral Emission Altitude Profiles

Model volume emission rates for principal auroral lines as functions of precipitating electron energy.

Aurora

Python

Parameters

Characteristic Energy (keV)

Number Flux (cm⁻²s⁻¹)

script.py92 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

E_CHAR = 5.0
E_char_keV = E_CHAR
FLUX_VAL = 1e9
flux = FLUX_VAL  # electrons/cm²/s

# Altitude grid
z = np.linspace(80, 400, 500)

# MSIS-like neutral density profiles (simplified)
n_O = 5e11 * np.exp(-(z-200)/50)  # cm^-3, atomic oxygen
n_N2 = 1e12 * np.exp(-(z-120)/15)  # cm^-3, N2
n_O2 = 3e11 * np.exp(-(z-110)/10)  # cm^-3, O2

# Energy deposition (Rees 1963 parameterization)
# Peak altitude: z_peak ≈ 150 - 20*ln(E_keV)
z_peak = 150 - 20 * np.log(E_char_keV)
H_dep = 15  # deposition scale height

# Normalized energy deposition profile
q_dep = np.exp(-0.5*((z - z_peak)/H_dep)**2)
q_dep = q_dep / np.max(q_dep)

# Volume emission rates (simplified)
# 557.7 nm: O(1S->1D), from O excitation
k_557 = 1.3e-17  # excitation cross section cm^2
VER_557 = flux * E_char_keV * k_557 * n_O * q_dep / 35  # photons/cm³/s

# 630.0 nm: O(1D->3P), from O excitation
k_630 = 2.2e-18
VER_630 = flux * E_char_keV * k_630 * n_O * q_dep / 35
# Quenching below ~200 km
tau_630 = 110  # s lifetime
n_quench = 2e11  # quenching density
quench_factor = 1.0 / (1.0 + n_N2 * 2e-17 * tau_630)
VER_630 = VER_630 * quench_factor

# 427.8 nm: N2+(1NG) from N2 ionization
k_428 = 3.5e-17
VER_428 = flux * E_char_keV * k_428 * n_N2 * q_dep / 35

# Panel 1: Energy deposition profile
axes[0].plot(q_dep, z, color='#f87171', linewidth=2)
axes[0].axhline(z_peak, color='#fbbf24', linestyle='--', linewidth=1, label=f'Peak: {z_peak:.0f} km')
axes[0].fill_betweenx(z, 0, q_dep, alpha=0.15, color='#f87171')
axes[0].set_xlabel('Normalized Deposition', color='#cbd5e1')
axes[0].set_ylabel('Altitude (km)', color='#cbd5e1')
axes[0].set_title(f'Energy Deposition (E = {E_char_keV} keV)', color='#f87171', fontweight='bold')
axes[0].legend(fontsize=9, facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')

# Panel 2: Volume emission rates
axes[1].semilogx(VER_557, z, color='#34d399', linewidth=2, label='557.7 nm (green)')
axes[1].semilogx(VER_630, z, color='#ef4444', linewidth=2, label='630.0 nm (red)')
axes[1].semilogx(VER_428, z, color='#818cf8', linewidth=2, label='427.8 nm (blue)')
axes[1].set_xlabel('VER (photons/cm³/s)', color='#cbd5e1')
axes[1].set_ylabel('Altitude (km)', color='#cbd5e1')
axes[1].set_title('Volume Emission Rates', color='#34d399', fontweight='bold')
axes[1].legend(fontsize=8, facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')

# Panel 3: Multiple energies comparison
energies = [0.5, 1, 5, 10, 30]
colors = ['#ef4444', '#f87171', '#fbbf24', '#34d399', '#60a5fa']
for E, col in zip(energies, colors):
    zp = 150 - 20 * np.log(E)
    q_e = np.exp(-0.5*((z - zp)/H_dep)**2)
    axes[2].plot(q_e, z, color=col, linewidth=1.5, label=f'{E} keV')

axes[2].set_xlabel('Normalized Deposition', color='#cbd5e1')
axes[2].set_ylabel('Altitude (km)', color='#cbd5e1')
axes[2].set_title('Deposition vs Electron Energy', color='#fbbf24', fontweight='bold')
axes[2].legend(fontsize=8, facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')

for ax in axes:
    ax.set_ylim(80, 350)

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a')
plt.close()
print(f"Auroral emission profiles for E = {E_char_keV} keV, peak at {z_peak:.0f} km")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

10. Tsyganenko-like Magnetic Field Model (Fortran)

Compute magnetospheric field in the noon-midnight plane using dipole + tail + ring current modules in Fortran.

Tsyganenko

Python

Parameters

Dst Index (nT)

Dynamic Pressure (nPa)

script.py157 lines

import subprocess, os, tempfile, sys
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

DST_VAL = -50.0
Dst = DST_VAL
DP_T = 2.0
Dp = DP_T

fortran_code = f"""
program tsyganenko_field
  implicit none
  integer, parameter :: NX = 200, NZ = 200
  real(8), parameter :: pi = 3.141592653589793d0
  real(8) :: B0, RE, x, z_val, r, theta, Bx_d, Bz_d, Bx_t, Bz_t
  real(8) :: Dst, Dp, B_tail, R_tail, delta_t
  real(8) :: Bx_mp, Bz_mp, R_mp, Bx_rc, Bz_rc
  real(8) :: Bx_tot, Bz_tot
  integer :: i, j

B0 = 3.12d-5
  RE = 6.371d6
  Dst = {Dst}d0
  Dp = {Dp}d0

! Magnetopause standoff
  R_mp = 10.22d0 * Dp**(-1.0d0/6.6d0)

! Tail parameters
  B_tail = 20.0d-9  ! nT -> T
  R_tail = 10.0d0  ! characteristic scale
  delta_t = 3.0d0  ! tail half-thickness in RE

write(*,'(A)') '# x(RE) z(RE) Bx(nT) Bz(nT) |B|(nT)'

do i = 1, NX
    do j = 1, NZ
      x = -30.0d0 + 40.0d0 * dble(i-1) / dble(NX-1)
      z_val = -15.0d0 + 30.0d0 * dble(j-1) / dble(NZ-1)

r = sqrt(x**2 + z_val**2)
      if (r < 1.0d0) cycle

theta = atan2(z_val, x)

! 1. Dipole field
      Bx_d = -3.0d0 * B0 * 1d9 / r**3 * x * z_val / r**2
      Bz_d = B0 * 1d9 / r**3 * (3.0d0 * z_val**2 / r**2 - 1.0d0)

! 2. Tail current contribution (Harris-like)
      Bx_t = B_tail * 1d9 * tanh(z_val / delta_t)
      if (x < -5.0d0) then
        Bx_t = Bx_t * (1.0d0 - exp((x + 5.0d0)/R_tail))
      else
        Bx_t = 0.0d0
      end if
      Bz_t = 0.0d0

! 3. Ring current (simplified as Bz depression)
      Bx_rc = 0.0d0
      Bz_rc = 0.0d0
      if (r > 2.0d0 .and. r < 8.0d0 .and. abs(z_val) < 2.0d0) then
        Bz_rc = Dst * 1d-9 * 1d9 * exp(-0.5d0*((r-4.5d0)/2.0d0)**2) * exp(-0.5d0*(z_val/1.0d0)**2)
      end if

! 4. Magnetopause currents (simplified compression)
      Bx_mp = 0.0d0
      Bz_mp = 0.0d0
      if (x > 0 .and. r > R_mp * 0.8d0) then
        Bz_mp = -Bz_d * 0.3d0 * exp(-(r - R_mp)**2 / 4.0d0)
      end if

Bx_tot = Bx_d + Bx_t + Bx_rc + Bx_mp
      Bz_tot = Bz_d + Bz_t + Bz_rc + Bz_mp

write(*,'(5F12.4)') x, z_val, Bx_tot, Bz_tot, sqrt(Bx_tot**2 + Bz_tot**2)
    end do
  end do
end program
"""

with tempfile.TemporaryDirectory() as tmpdir:
    src = os.path.join(tmpdir, 'tsyg.f90')
    exe = os.path.join(tmpdir, 'tsyg')
    with open(src, 'w') as f:
        f.write(fortran_code)
    comp = subprocess.run(['gfortran', '-O2', '-o', exe, src], capture_output=True, text=True)
    if comp.returncode != 0:
        print("Compilation error:", comp.stderr)
        sys.exit(1)
    result = subprocess.run([exe], capture_output=True, text=True, timeout=60)
    lines = result.stdout.strip().split('\n')

data_lines = [l for l in lines if not l.startswith('#') and l.strip()]
data = np.array([[float(v) for v in line.split()] for line in data_lines])

x = data[:,0]; z = data[:,1]; Bx = data[:,2]; Bz = data[:,3]; Bmag = data[:,4]

# Reshape
nx = len(np.unique(x)); nz = len(np.unique(z))
x2d = x.reshape(nx, nz); z2d = z.reshape(nx, nz)
Bx2d = Bx.reshape(nx, nz); Bz2d = Bz.reshape(nx, nz)
Bmag2d = Bmag.reshape(nx, nz)

fig, axes = plt.subplots(1, 2, figsize=(16, 7))
fig.patch.set_facecolor('#0f172a')
for ax in axes:
    ax.set_facecolor('#1e293b')
    ax.tick_params(colors='#cbd5e1')
    for spine in ax.spines.values():
        spine.set_color('#334155')

# Panel 1: Field magnitude
c1 = axes[0].pcolormesh(x2d, z2d, np.log10(np.clip(Bmag2d, 0.1, 1e4)), cmap='inferno', shading='auto')
cb1 = fig.colorbar(c1, ax=axes[0], shrink=0.8)
cb1.set_label('log₁₀|B| (nT)', color='#cbd5e1')
cb1.ax.yaxis.set_tick_params(color='#cbd5e1')
plt.setp(plt.getp(cb1.ax.axes, 'yticklabels'), color='#cbd5e1')

earth = plt.Circle((0,0), 1, color='#34d399', fill=True, zorder=10)
axes[0].add_patch(earth)
axes[0].set_xlabel('X_GSM (R_E)', color='#cbd5e1')
axes[0].set_ylabel('Z_GSM (R_E)', color='#cbd5e1')
axes[0].set_title(f'|B| Field (Dst={Dst:.0f} nT, Dp={Dp:.1f} nPa)', color='#f87171', fontweight='bold')
axes[0].set_aspect('equal')

# Panel 2: Field lines (streamplot)
# Need uniform grid for streamplot
from scipy.interpolate import griddata
xi = np.linspace(-30, 10, 150)
zi = np.linspace(-15, 15, 150)
Xi, Zi = np.meshgrid(xi, zi)
Bxi = griddata((x, z), Bx, (Xi, Zi), method='linear', fill_value=0)
Bzi = griddata((x, z), Bz, (Xi, Zi), method='linear', fill_value=0)

speed = np.sqrt(Bxi**2 + Bzi**2)
lw = 0.5 + 1.5 * speed / np.percentile(speed[speed>0], 95)
lw = np.clip(lw, 0.3, 2)

axes[1].streamplot(xi, zi, Bxi, Bzi, color=np.log10(np.clip(speed, 0.1, 1e4)),
                   cmap='plasma', linewidth=lw, density=2, arrowsize=0.8)
earth2 = plt.Circle((0,0), 1, color='#34d399', fill=True, zorder=10)
axes[1].add_patch(earth2)
axes[1].set_xlabel('X_GSM (R_E)', color='#cbd5e1')
axes[1].set_ylabel('Z_GSM (R_E)', color='#cbd5e1')
axes[1].set_title('Magnetic Field Lines (Noon-Midnight Plane)', color='#a78bfa', fontweight='bold')
axes[1].set_aspect('equal')
axes[1].set_xlim(-30, 10)
axes[1].set_ylim(-15, 15)

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a')
plt.close()
print(f"Tsyganenko-like field computed: Dst={Dst} nT, Dp={Dp} nPa")

Click Run to execute the Python code

Code will be executed with Python 3 on the server