import datetime
import numpy as np
from numba import njit, prange
from gtsimulation.common import Units, Regions
from gtsimulation.magnetic_field import AbsBfield
from gtsimulation.magnetic_field.heliosphere.Functions import transformations
[docs]
class Parker(AbsBfield):
"""
Parker spiral model of the interplanetary magnetic field with optional
heliospheric current sheet (HCS), corotating interaction regions (CIR),
and turbulent fluctuations (slab + 2D components).
The regular part is the classical Parker field with a possible HCS
whose tilt angle varies with the solar cycle. Turbulence is simulated
as a superposition of slab and 2D modes following the approach described
in [1]_.
Parameters
----------
date : datetime.date or int, optional
Initial epoch. If int, interpreted as seconds from an arbitrary
reference; if datetime.date, converted to Unix time (seconds since
1970-01-01). Default 0.
magnitude : float, default=2.09
Magnetic field magnitude at 1 AU (nT).
polarity : int, default=-1
Polarity sign of the field (+1 or -1).
tilt_angle : float or None, optional
Constant tilt angle of the heliospheric current sheet (radians).
If None (default), the tilt angle varies with the 11‑year solar cycle
according to an empirical fit.
use_reg : bool, default=True
Whether to include the regular Parker field (including HCS if enabled).
use_hcs : bool, default=True
Whether to include the heliospheric current sheet modulation.
use_cir : bool, default=False
Whether to include CIR effects (not implemented yet).
use_noise : bool, default=False
Whether to add turbulent fluctuations.
use_slab : bool, default=True
Whether to include the slab component of turbulence.
use_2d : bool, default=True
Whether to include the 2D component of turbulence.
noise_num : int, default=256
Number of modes used in the turbulence spectrum.
log_kmin : float, default=1
Decimal logarithm of the minimum wave number (AU⁻¹).
log_kmax : float, default=6
Decimal logarithm of the maximum wave number (AU⁻¹).
coeff_noise : float, default=0.47
Scaling factor for the turbulent component.
coeff_2d : float, default=2.9
Scaling factor for the 2D component of turbulence.
**kwargs
Additional arguments passed to the base class `AbsBfield`.
Notes
-----
The regular Parker field is given by:
B_r = A0 / r² * HCS
B_φ = -A0 / r² * ((r - rs) ω / v_sw) sinθ * HCS
where HCS is the current sheet modulation (tanh profile).
Turbulence is simulated as a sum of slab and 2D modes following the
synthetic turbulence model described in [1]_. The implementation
closely follows the formulas provided there.
References
----------
.. [1] Engelbrecht, N. E., et al. (2022). ApJ, 941(2), 168.
doi:10.3847/1538-4357/aca892
Examples
--------
>>> from gtsimulation.magnetic_field.heliosphere import Parker
Create a model with default parameters (regular field only, no noise).
Compute the field at Earth's orbit (1 AU along the X-axis).
>>> model = Parker()
>>> x, y, z = 1.0, 0.0, 0.0 # coordinates in AU
>>> Bx, By, Bz = model.CalcBfield(x, y, z)
>>> print(f"B = ({Bx:.2f}, {By:.2f}, {Bz:.2f}) nT")
Override the HCS tilt angle with a constant value (e.g., 0.5 rad).
>>> model_tilt = Parker(tilt_angle=0.5)
>>> Bx_t, By_t, Bz_t = model_tilt.CalcBfield(-x, y, z)
>>> print(f"B = ({Bx:.2f}, {By:.2f}, {Bz:.2f}) nT")
"""
ToMeters = Units.AU2m
rs = 0.0232523
omega = 2 * np.pi / 2160000
years11 = 347133600
km2AU = 1 / Units.AU2km
def __init__(self, date: datetime.date | int = 0, magnitude=2.09, polarity=-1, tilt_angle=None,
use_reg=True, use_hcs=True, use_cir=False, use_noise=False, use_slab=True, use_2d=True,
noise_num=256, log_kmin=1, log_kmax=6, coeff_noise=0.47, coeff_2d=2.9, **kwargs):
super().__init__(**kwargs)
self.Region = Regions.Heliosphere
self.ModelName = "Parker"
self.Units = "AU"
self.magnitude = magnitude
self.polarity = polarity
self.tilt_angle = tilt_angle
self.use_reg = use_reg
self.use_hcs = use_hcs
self.use_cir = use_cir
self.use_slab = use_slab
self.use_2d = use_2d
self.coeff_noise = coeff_noise
self.coeff_2d = coeff_2d
self.__set_time(date)
self.__set_noise(use_noise, noise_num, log_kmin, log_kmax)
def __set_time(self, date: int | datetime.datetime):
self.Date = date
self.t = 2488320 # To have a correct phase with ace data
if isinstance(date, int):
self.t += date
return
year = date.year
doy = date.timetuple().tm_yday
hour = date.hour
minute = date.minute
second = date.second
self.t += (((year * 365.25 + doy) * 24 + hour) * 60 + minute) * 60 + second
[docs]
def CalcBfield(self, x, y, z, **kwargs):
if kwargs.get("t") is not None:
self.t = kwargs.get("t")
A0 = self.magnitude * self.polarity
t = self.t
r, R, theta, phi = transformations.Cart2Sphere(x, y, z)
v_wind = self.v_wind(theta, self.km2AU)
omega = self.omega
rs = self.rs
use_hcs = self.use_hcs
Bx = np.zeros_like(r)
By = np.zeros_like(r)
Bz = np.zeros_like(r)
if self.use_reg:
years11 = self.years11
if self.tilt_angle is None:
alpha = self.CalcTiltAngle(t)
dalpha = np.sign(self.CalcTiltAngle(t + 1) - self.CalcTiltAngle(t - 1))
else:
alpha = self.tilt_angle
dalpha = 0.0
Br, Bphi = self._calc_regular(A0, t, r, theta, phi, v_wind, omega, rs, years11, alpha, dalpha, use_hcs)
# alpha -= np.pi / years11 * (r - rs) / v_wind * dalpha
#
# theta0 = np.pi / 2 - np.arctan(-np.tan(alpha) * np.sin(phi + omega * (r - rs) / v_wind -
# omega * t))
#
# HCS = self.HCS(theta, theta0, r) * (r >= rs)
# Br = A0 / r ** 2 * HCS
# Bphi = -A0 / r ** 2 * (((r - rs) * omega) / v_wind) * np.sin(theta) * HCS
Bx, By, Bz = transformations.Sphere2Cart(Br, 0, Bphi, theta, phi)
if self.use_cir:
# TODO: add CIR
pass
if not self.use_noise:
return Bx, By, Bz
# coeff2d = 1.4
# coeffslab = coeff2d / 2
a = v_wind / omega
A_rad = self.A_rad
alpha_rad = self.alpha_rad
delta_rad = self.delta_rad
A_azimuth = self.A_azimuth
alpha_azimuth = self.alpha_azimuth
delta_azimuth = self.delta_azimuth
A_2D = self.A_2D
alpha_2D = self.alpha_2D
delta_2D = self.delta_2D
k = self.k
dk = self.dk
Bx_n, By_n, Bz_n = self._calc_noise(r, theta, phi, a,
A_rad, alpha_rad, delta_rad,
A_azimuth, alpha_azimuth, delta_azimuth,
A_2D, alpha_2D, delta_2D,
rs, k, dk, self.use_slab, self.use_2d, self.coeff_2d)
Bx += self.magnitude * self.coeff_noise * Bx_n
By += self.magnitude * self.coeff_noise * By_n
Bz += self.magnitude * self.coeff_noise * Bz_n
return Bx, By, Bz
[docs]
def UpdateState(self, new_date):
self.__set_time(new_date)
@staticmethod
@njit(fastmath=True)
def _calc_regular(A0, t, r, theta, phi, v_wind, omega, rs, years11, alpha, dalpha, use_hcs):
HCS = 1.
if use_hcs:
alpha_n = alpha - np.pi / years11 * (r - rs) / v_wind * dalpha
theta0 = np.pi / 2 - np.arctan(-np.tan(alpha_n) * np.sin(phi + omega * (r - rs) / v_wind -
omega * t))
L = 0.0002
dt = r * (theta - theta0) / L
HCS = -np.tanh(dt)
Br = A0 / r ** 2 * HCS
Bphi = -A0 / r ** 2 * (((r - rs) * omega) / v_wind) * np.sin(theta) * HCS
return Br, Bphi
[docs]
@staticmethod
@njit(fastmath=True)
def CalcTiltAngle(t):
a0 = 0.7502
a1 = 0.02332
b1 = -0.01626
a2 = -0.3268
b2 = 0.2016
a3 = -0.02814
b3 = 0.0005215
a4 = -0.08341
b4 = -0.04852
w = 9.318e-09
alpha = (a0 +
a1 * np.cos(t * w) + b1 * np.sin(t * w) +
a2 * np.cos(2 * t * w) + b2 * np.sin(2 * t * w) +
a3 * np.cos(3 * t * w) + b3 * np.sin(3 * t * w) +
a4 * np.cos(4 * t * w) + b4 * np.sin(4 * t * w))
return alpha
[docs]
@staticmethod
@njit(fastmath=True)
def v_wind(theta, km2AU):
return (300 + 475 * (1 - np.sin(theta) ** 8)) * km2AU
[docs]
@classmethod
def a(cls, theta):
return cls.v_wind(theta) / cls.omega
[docs]
@staticmethod
@njit(fastmath=True)
def HCS(theta, theta0, r):
L = 0.0002
dt = r * (theta - theta0) / L
return -np.tanh(dt)
def __set_noise(self, use_noise, noise_num, log_kmin, log_kmax):
self.use_noise = use_noise
if self.use_noise is False:
return
self.noise_num = noise_num
self.log_kmin = log_kmin
self.log_kmax = log_kmax
self.k = np.logspace(log_kmin, log_kmax, self.noise_num)[:, np.newaxis]
self.dk = self.k * (10 ** ((log_kmax - log_kmin) / (self.noise_num - 1)) - 1)
self.A_2D = np.random.randn(self.noise_num, 1) / 130
self.alpha_2D = np.random.rand(self.noise_num, 1) * 2 * np.pi
n = np.trunc(np.sin(self.alpha_2D) * self.k)
self.alpha_2D = np.real(np.arcsin(n / self.k) * (np.cos(self.alpha_2D) > 0) +
(np.pi - np.arcsin(n / self.k)) * (np.cos(self.alpha_2D) < 0))
self.delta_2D = np.random.rand(self.noise_num, 1) * 2 * np.pi
self.A_rad = np.random.randn(self.noise_num, 1) / 1.5
self.alpha_rad = np.random.rand(self.noise_num, 1) * 2 * np.pi
n = np.trunc(np.sin(self.alpha_rad) * self.k)
self.alpha_rad = np.real(np.arcsin(n / self.k) * (np.cos(self.alpha_rad) > 0) +
(np.pi - np.arcsin(n / self.k)) * (np.cos(self.alpha_rad) < 0))
self.delta_rad = np.random.rand(self.noise_num, 1) * 2 * np.pi
self.A_azimuth = np.random.randn(self.noise_num, 1) / 4.5
self.alpha_azimuth = np.random.rand(self.noise_num, 1) * 2 * np.pi
n = np.trunc(np.sin(self.alpha_azimuth) * self.k)
self.alpha_azimuth = np.real(np.arcsin(n / self.k) * (np.cos(self.alpha_azimuth) > 0) +
(np.pi - np.arcsin(n / self.k)) * (np.cos(self.alpha_azimuth) < 0))
self.delta_azimuth = np.random.rand(self.noise_num, 1) * 2 * np.pi
@staticmethod
@njit(fastmath=True)
def _calc_noise(r, theta, phi, a,
A_rad, alpha_rad, delta_rad,
A_azimuth, alpha_azimuth, delta_azimuth,
A_2d, alpha_2d, delta_2d,
rs, k, dk, use_slab, use_2d, component_2d):
"""
doi.org/10.3847/1538-4357/aca892/meta
"""
q_slab = 5 / 3
q_2d = 8 / 3
p = 0
gamma = 3
cospsi = 1. / np.sqrt(1 + ((r - rs) * np.sin(theta) / a) ** 2)
sinpsi = ((r - rs) * np.sin(theta) / a) / np.sqrt(1 + ((r - rs) * np.sin(theta) / a) ** 2)
cospsi_ = 1. / np.sqrt(1 + ((r - rs) / a) ** 2)
sinpsi_ = ((r - rs) / a) / np.sqrt(1 + ((r - rs) / a) ** 2)
lam_2d = 0.04 * (r / (rs / 5)) ** 0.8 * (rs / 5)
dlamd_2d = 0.032 * (rs / (5 * r)) ** 0.2
lam_slab = 2 * lam_2d
dlam_slab = 2 * dlamd_2d
Br, Btheta, Bphi = 0., 0., 0.
# TODO: calculation is point wise
for mod in prange(len(k)):
numer_slab = dk[mod, 0] * k[mod, 0] ** p
# Radial spectrum
B_rad = A_rad[mod, 0] * r ** (-gamma / 2)
brk_rad = lam_slab * k[mod, 0] / np.sqrt(a * r)
denom_rad = (1 + brk_rad ** (p + q_slab))
spectrum_rad = np.sqrt(numer_slab / denom_rad)
deltaB_rad = 2 * B_rad * spectrum_rad * cospsi_ * r * np.sqrt(r * a)
# Azimuthal spectrum
B_azimuth = A_azimuth[mod, 0] * r ** (-gamma / 2)
brk_azimuth = lam_slab * k[mod, 0] / r
denom_azimuth = (1 + brk_azimuth ** (p + q_slab))
spectrum_azimuth = np.sqrt(numer_slab / denom_azimuth)
deltaB_azimuth = B_azimuth * spectrum_azimuth
dspectrum_azimuth = -spectrum_azimuth * (p + q_2d) * (denom_azimuth - 1) * (r * dlam_slab - lam_slab) / (
denom_azimuth * 2 * r * lam_slab)
ddeltaB_azimtuth = B_azimuth * dspectrum_azimuth + spectrum_azimuth * B_azimuth * (-gamma / (2 * r))
# 2d spectrum
B_2d = A_2d[mod, 0] * r ** (-gamma / 2)
brk_2d = lam_2d * k[mod, 0] / r
denom_2d = (1 + brk_2d ** (p + q_2d))
numer_2d = dk[mod, 0] * k[mod, 0] ** (p + 1)
spectrum_2d = np.sqrt(2 * np.pi * numer_2d / denom_2d)
deltaB_2d = B_2d * spectrum_2d
dspectrum_2d = -spectrum_2d * (p + q_2d) * (denom_2d - 1) * (r * dlamd_2d - lam_2d) / (
denom_2d * 2 * r * lam_2d)
ddeltaB_2d = B_2d * dspectrum_2d + spectrum_2d * B_2d * (-gamma / (2 * r))
# Radial polarization and phase
phase_rad = k[mod, 0] * np.sqrt(r / a) + delta_rad[mod, 0]
# Azimuthal polarization and phase
phase_azimuth = k[mod, 0] * phi + delta_azimuth[mod, 0]
# 2d polarization and phase
phase_2d = k[mod, 0] * ((r / a + phi) * np.sin(alpha_2d[mod, 0]) + theta * np.cos(alpha_2d[mod, 0])) + \
delta_2d[mod, 0]
# Radial field
Br_rad = 0
Btheta_rad = -deltaB_rad * a * np.sin(alpha_rad[mod, 0]) * np.cos(phase_rad) / (
2 * r * np.sin(theta) * np.sqrt(a * r))
Bphi_rad = deltaB_rad * a * np.cos(alpha_rad[mod, 0]) * np.cos(phase_rad) / (
2 * r * np.sin(theta) * np.sqrt(a * r))
# Azimuthal field
Br_az = -deltaB_azimuth * sinpsi_ * np.cos(alpha_azimuth[mod, 0]) * np.cos(phase_azimuth)
Btheta_az = deltaB_azimuth * sinpsi_ * np.sin(alpha_azimuth[mod, 0]) * np.cos(phase_azimuth)
Bphi_az = 1/k[mod, 0] * (np.sin(theta) * np.sin(phase_azimuth) * np.cos(alpha_azimuth[mod, 0]) *
(2*deltaB_azimuth*sinpsi_ + r/a * deltaB_azimuth * cospsi_ + r * sinpsi_ * ddeltaB_azimtuth) -
np.cos(theta)*deltaB_azimuth*np.sin(phase_azimuth)*sinpsi_*np.sin(alpha_azimuth[mod, 0]))
# 2d field
Br_2d = -deltaB_2d / (r * k[mod, 0] ) * (np.sin(phase_2d)*sinpsi*np.tan(theta)**(-1) +
k[mod, 0]*np.cos(alpha_2d[mod, 0])*np.cos(phase_2d)*sinpsi +
np.sin(phase_2d)*sinpsi*cospsi**2*np.tan(theta)**(-1))
Btheta_2d = deltaB_2d / (r * np.sin(theta)) * cospsi * np.sin(alpha_2d[mod, 0] * np.cos(phase_2d)) \
- np.sin(theta) * cospsi / (a * r * k[mod, 0]) * (ddeltaB_2d * r * (r - rs) * np.sin(phase_2d) +
deltaB_2d * np.sin(phase_2d) * (
2 * r - rs - r * sinpsi ** 2) +
k[mod, 0] * r * (r - rs) / a *
np.sin(alpha_2d[mod, 0]) * np.cos(phase_2d) *
deltaB_2d)
Bphi_2d = -deltaB_2d / (r * k[mod, 0]) * (cospsi * k[mod, 0] * np.cos(alpha_2d[mod, 0]) * np.cos(phase_2d) -
(np.tan(theta))**(-1) * np.sin(phase_2d) * cospsi * sinpsi**2)
# Total field
coeff_slab = 0
coeff_2d = 0
if use_slab:
coeff_slab = 1
if use_2d:
coeff_2d = component_2d
Br += coeff_2d*Br_2d + coeff_slab * (Br_az + Br_rad)
Btheta += coeff_2d*Btheta_2d + coeff_slab * (Btheta_az + Btheta_rad)
Bphi += coeff_2d*Bphi_2d + coeff_slab * (Bphi_az + Bphi_rad)
# B = np.zeros((3, *Br_2d.shape))
# B[0] = Br_2d + coeff_slab * (Br_az + Br_rad)
# B[1] = Btheta_2d + coeff_slab * (Btheta_az + Btheta_rad)
# B[2] = Bphi_2d + coeff_slab * (Bphi_az + Bphi_rad)
# B_s = np.sum(B, axis=1)
# Br, Btheta, Bphi = B_s[0], B_s[1], B_s[2]
# Br = np.sum(Br_2d + coeff_slab * (Br_az + Br_rad), axis=0)
# Btheta = np.sum(Btheta_2d + coeff_slab * (Btheta_az + Btheta_rad), axis=0)
# Bphi = np.sum(Bphi_2d + coeff_slab * (Bphi_az + Bphi_rad), axis=0)
Bx, By, Bz = transformations.Sphere2Cart(Br, Btheta, Bphi, theta, phi)
return Bx * (r > rs), By * (r > rs), Bz * (r > rs)
[docs]
def to_string(self):
s = f"""Parker
Regular: {self.use_reg}
Magnitude: {self.magnitude}
HCS: {self.use_hcs}
CIR: {self.use_cir}
Polarity: {self.polarity}
Noise: {self.use_noise}
"""
if self.use_noise:
s += f"""
Min wave length: {self.log_kmin}
Max wave length: {self.log_kmax}
Number of waves: {self.noise_num}
Coeff_Noise: {self.coeff_noise}
Coeff_2d: {self.coeff_2d}
Using Slab: {self.use_slab}
Using 2D: {self.use_2d}"""
return s
