emp.model

Copyright (C) 2023 by The RAND Corporation See LICENSE and README.md for information on usage and licensing

Contains the EmpModel class for simulating EMP effects, as well as the EmpLosResult dataclass for storing results.

Classes

EmpLosResult(time_points, E_theta_at_ground, ...)

Result dataclass for EMP line-of-sight calculations.

EmpModel(burst_point, target_point[, ...])

The EMP model class.

class emp.model.EmpLosResult(time_points, E_theta_at_ground, E_phi_at_ground, E_norm_at_ground, model_params, burst_point_dict, target_point_dict)[source]

Result dataclass for EMP line-of-sight calculations.

Contains time series data for electric field components and magnitude at ground level, along with model parameters used for the calculation.

time_points: ndarray[Any, dtype[float64]]
E_theta_at_ground: List[float]
E_phi_at_ground: List[float]
E_norm_at_ground: List[float]
model_params: Dict[str, Any]
burst_point_dict: Dict[str, float]
target_point_dict: Dict[str, float]
save(filepath)[source]

Save the EmpLosResult to a JSON file.

Parameters:

filepath (Union[str, Path]) – Path where to save the result file.

Return type:

None

classmethod load(filepath)[source]

Load an EmpLosResult from a JSON file.

Parameters:

filepath (Union[str, Path]) – Path to the result file to load.

Returns:

Loaded result object.

Return type:

EmpLosResult

get_max_field_magnitude()[source]

Get the maximum field magnitude from the result.

Return type:

float

get_max_field_time()[source]

Get the time at which maximum field occurs.

Return type:

float

__init__(time_points, E_theta_at_ground, E_phi_at_ground, E_norm_at_ground, model_params, burst_point_dict, target_point_dict)
class emp.model.EmpModel(burst_point, target_point, total_yield_kt=5.0, gamma_yield_fraction=0.05, Compton_KE=1.28, pulse_param_a=0.01, pulse_param_b=0.37, rtol=0.0001, numerical_integration_method='Radau', magnetic_field_model='dipole', magnetic_field_date=None, time_max=100.0, num_time_points=300)[source]

The EMP model class.

__init__(burst_point, target_point, total_yield_kt=5.0, gamma_yield_fraction=0.05, Compton_KE=1.28, pulse_param_a=0.01, pulse_param_b=0.37, rtol=0.0001, numerical_integration_method='Radau', magnetic_field_model='dipole', magnetic_field_date=None, time_max=100.0, num_time_points=300)[source]

Init method.

Parameters:
  • burst_point (geometry.Point) – Nuclear burst location.

  • target_point (geometry.Point) – Target point where EMP effects are calculated.

  • total_yield_kt (float, optional) – Total yield in kilotons. By default DEFAULT_total_yield_kt.

  • gamma_yield_fraction (float, optional) – Fraction of yield deposited in prompt Gamma radiation. By default DEFAULT_gamma_yield_fraction.

  • Compton_KE (float, optional) – Kinetic energy of Compton electron in MeV. By default DEFAULT_Compton_KE.

  • pulse_param_a (float, optional) – Pulse parameter in 1/ns. By default, DEFAULT_pulse_param_a.

  • pulse_param_b (float, optional) – Pulse parameter in 1/ns. By default, DEFAULT_pulse_param_b.

  • rtol (float, optional) – Relative tolerance for ODE integration. By default, DEFAULT_rtol.

  • numerical_integration_method (str, optional) – Integration method used for solve_ivp. By default, ‘Radau’, which is well-suited for stiff problems.

  • magnetic_field_model (Union[str, MagneticFieldModel], optional) – Magnetic field model to use (‘dipole’ or ‘igrf’). By default ‘dipole’.

  • magnetic_field_date (datetime, optional) – Date for the magnetic field model, if applicable.

  • time_max (float, optional) – Max time to integrate to in ns. By default 100.0 ns.

  • num_time_points (int, optional) – Number of time points to compute. By default 300.

Yields:

None – No returns.

RCompton(r)[source]

The Compton electron stopping distance.

Parameters:

r (float) – Radius from burst to point of interest, in km.

Returns:

Stopping distance.

Return type:

float

TCompton(r)[source]

The (scaled) Compton electron lifetime.

Parameters:

r (float) – Radius from burst to point of interest, in km.

Returns:

Compton electron lifetime, in ns.

Return type:

float

f_pulse(t)[source]

Normalized gamma pulse, for the difference of exponential form used by Seiler. This parameterization of the pulse profile has the benefit that it allows many integrals to be analytically solved. In general, other parameterizations are possible but they will change the equations derived by Seiler.

Parameters:

t (Union[float, ndarray]) – Time, in ns.

Returns:

Gamma flux pulse profile, in 1/ns.

Return type:

Union[float, ndarray]

rho_divided_by_rho0(r)[source]

Ratio of air density at radius r to air density at sea level.

Parameters:

r (float) – Radius from burst to point of interest, in km.

Returns:

Air density ratio, dimensionless.

Return type:

float

mean_free_path(r)[source]

Compton electron mean free path.

Parameters:

r (float) – Radius from burst to point of interest, in km.

Returns:

Mean free path, in km.

Return type:

float

gCompton(r)[source]

The g function for the creation of primary electrons, as introduced in Karzas, Latter Eq (4). The radius is measured from the burst.

Parameters:

r (float) – Radius from burst to point of interest, in km.

Returns:

Compton g function, in km^(-3)

Return type:

float

gCompton_numerical_integration(r)[source]

The g function for the creation of primary electrons The radius is measured from the burst. Here g is computed using a numerical integration. TO DO: - can we vectorize the numerical integral w/ scipy? - should we increase the precision?

Parameters:

r (float) – Radius from burst to point of interest, in km.

Returns:

Compton g function, in km^(-3)

Return type:

Union[float, NDArray[np.float64]]

electron_collision_freq_at_sea_level(E, t)[source]

Electron collision frequency at sea level. Defined in Seiler eq 55. I modified his expression so t is in units of ns.

Parameters:
  • E (float) – Norm of the electric field, in V/m.

  • t (float) – Evaluation retarded time, in ns.

Returns:

Electron collision frequency, in 1/ns.

Return type:

float

conductivity(r, t, nuC_0)[source]

The conductivity computed using Seiler’s approximation. Defined in Seiler eq. 67 and 70. The expression in eq. 70 has been simlified.

Parameters:
  • r (float) – Radius from burst to point of interest, in km.

  • t (float) – Evaluation retarded time, in ns.

  • nuC_0 (float) – Ground-level electron collision frequency, in units of 1/ns.

Returns:

Conductivity, in units of C^2 s m^{-3} kg^{-1}.

Return type:

float

JCompton_theta(r, t)[source]

The theta component of the Compton current, computed using Seiler’s approximation. Defined in Seiler eq. 65 and 68.

Parameters:
  • r (float) – Radius from burst to point of interest, in km.

  • t (float) – Evaluation retarded time, in ns.

Returns:

J_theta, in units of C s^{-1} m^{-2}.

Return type:

float

JCompton_phi(r, t)[source]

The azimuthal component of the Compton current, computed using Seiler’s approximation. Defined in Seiler eq. 66 and 69.

Parameters:
  • r (float) – Radius from burst to point of interest, in km.

  • t (float) – Evaluation retarded time, in ns.

Returns:

J_phi, in units of C s^{-1} m^{-2}.

Return type:

float

JCompton_theta_KL(r, t)[source]

The theta component of the Compton current, computed using the KL approximation. Defined in KL eq. 15.

Parameters:
  • r (float) – Radius from burst to point of interest, in km.

  • t (float) – Evaluation retarded time, in ns.

Returns:

J_theta, in units of C s^{-1} m^{-2}.

Return type:

float

JCompton_phi_KL(r, t)[source]

The phi component of the Compton current, computed using the KL approximation. Defined in KL eq. 16.

Parameters:
  • r (float) – Radius from burst to point of interest, in km.

  • t (float) – Evaluation retarded time, in ns.

Returns:

J_phi, in units of C s^{-1} m^{-2}.

Return type:

float

conductivity_KL(r, t, nuC_0)[source]

The conductivity computed using the KL approximation. Defined in KL eq. 53 and 13.

Parameters:
  • r (float) – Radius from burst to point of interest, in km.

  • t (float) – Evaluation retarded time, in ns.

  • nuC_0 (float) – Ground-level electron collision frequency, in units of 1/ns.

Returns:

Conductivity, in units of C^2 s m^{-3} kg^{-1}.

Return type:

float

F_theta_Seiler(E, r, t, nuC_0)[source]

The theta-component of the Maxwell equations, expressed as dE/dr = F(E).

Parameters:
  • E (float) – Electric field, in V/m.

  • r (float) – Radius from burst to point of interest, in km.

  • t (float) – Evaluation retarded time, in ns.

  • nuC_0 (Callable[[float], float]) – Ground-level electron collision frequency function, in units of 1/ns.

Returns:

Returns the RHS of the ODE. Units are in (V/m/km). The factor of 1e3 is to account for the fact that r is measured in km.

Return type:

float

F_phi_Seiler(E, r, t, nuC_0)[source]

The phi-component of the Maxwell equations, expressed as dE/dr = F(E).

Parameters:
  • E (float) – Electric field, in V/m.

  • r (float) – Radius from burst to point of interest, in km.

  • t (float) – Evaluation retarded time, in ns.

  • nuC_0 (Callable[[float], float]) – Ground-level electron collision frequency function, in units of 1/ns.

Returns:

Returns the RHS of the ODE. Units are in (V/m/km). The factor of 1e3 is to account for the fact that r is measured in km.

Return type:

float

solve_KL_ODEs(t, nuC_0)[source]

Solve the angular components of the KL ODEs.

TO DO: consider adding radial component

Parameters:
  • t (float) – Evaluation retarded time, in ns.

  • nuC_0 (Callable[[float], float]) – Ground-level electron collision frequency function, in units of 1/ns.

Returns:

A result tuple for each component (theta, phi).

Return type:

Tuple[OdeResult, OdeResult]

run()[source]

Solve the KL equations using the Seiler approximations for the source terms for a range of retarded times and return the angular components of the electric field for r = r_target (at the Earth’s surface).

Returns:

Result object containing time series of the components and norm of the E-field evaluated at a target point on the Earth’s surface.

Return type:

EmpLosResult

to_yaml(filepath)[source]

Export the EmpModel configuration to a YAML file.

Parameters:

filepath (Union[str, Path]) – Path where to save the YAML configuration file.

Return type:

None

classmethod from_yaml(filepath)[source]

Create an EmpModel instance from a YAML configuration file.

Parameters:

filepath (Union[str, Path]) – Path to the YAML configuration file.

Returns:

New EmpModel instance configured from the YAML file.

Return type:

EmpModel

Raises:
  • FileNotFoundError – If the configuration file doesn’t exist.

  • KeyError – If required configuration parameters are missing.

  • ValueError – If configuration parameters have invalid values.