import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import os
os.chdir('..')
import proppy as pp

Simulation setup

We use the same simulation from Tutorial 1 with a different propagation module.

sim = pp.Simulation()

start simulation

Source

First, we have to define a source of particles that we want to propagate. The simples source type is the point source that emmits particles isotropically. The only user-specified parameters are:

• energy: Energy of the particles in eV.
• source posistion: The position of the point source.
• number of particles: Number of particles that should be emitted from this source.

The source can be easily added to the simulation. Afterwards, calling the description on the source of the simulation prints all relevant information and the values of the source parameters.

nr_particles = 10**3
source_pos = np.array([0.0, 0.0, 0.0], dtype=np.float32)
energy = 3*10**15 # eV

source = pp.PointSourceIsotropicPhi(energy, source_pos, nr_particles)
sim.add_source(source)
sim.source.get_description()

Description Source:
                The source defines the start conditions of the particles 
                and covers the position, direction, energy, etc

position:  [0. 0. 0.]
number particles:  1000
energy:  3000000000000000  eV
source tpye: PointSourceIsotropic

Propagator

Propagates particles via a correlated rrandom walk. The overall behaviour is governed by a generalized telegraph equation: $\frac{\partial f}{\partial t} = \sum_i \tau_i \frac{\partial^2 f}{\partial t^2} = \sum_i \kappa_i \frac{\partial^2 f}{\partial x_i^2}$.

Here, we use isotropic diffusion for simplicity that propagates particles using a correlated random walk in Cartesian coordinates. Isotropic diffusion is given when the turbulence is isotropic and there is no background field. In the following tutorials also the anisotropic diffusion is used. The only user-specified parameters for the simple case of isotropic diffusion are:

• mean-free paths: The mean-free paths $\lambda_i$ can be derived from the diffusion coefficients $\kappa_i$ via $\lambda_i = 3 \kappa_i/c$. Therfore, the diagonal elements of the diagonalized diffusion tensor are needed to determine the characteristics of the transport. For isotropic diffusion all diffusion coefficients are equal $\kappa = [\kappa_{xx}, \kappa_{yy}, \kappa_{zz}]$ with $\kappa_{xx}=\kappa_{yy}=\kappa_{zz}$. Typical free-mean paths of charged particles in plasmoids in AGN jets are $10^{12}$m (see Reichherzer et al. (2021)).
• number of steps: The number of simulation steps for each individual particle.
• step size: Size of an individual step. Together with the parameter number of steps, the step size determines the trajectory length of the particles.

The propagator can be easily added to the simulation. Afterwards, calling the description on the propagator of the simulation prints all relevant information and the values of the propagation parameters.

nr_steps = 1*10**5
step_size = 0.1*10**10 # [m]
speed_of_light = 3*10**8 # [m/s]
diffusion_coefficient_perp = 1.3*10**18 # [m^2/s]
diffusion_coefficient_para = 1.4*10**20 # [m^2/s]
mfp_perp = 3*diffusion_coefficient_perp/speed_of_light*2
mfp_para = 3*diffusion_coefficient_para/speed_of_light
mfp = np.array([mfp_perp, mfp_perp, mfp_para], dtype=np.float32)
rms = 1 # Gaus
magnetic_field = pp.OrderedBackgroundField(rms, [0,0,1]).magnetic_field

propagator = pp.AnisotropicPropagator(magnetic_field, mfp, nr_steps, step_size)
sim.add_propagator(propagator)
sim.propagator.get_description()

Propagator initialized
Description Propagator:
                The propagator is responsible for the movement of the particles. 
                It performs the change of direction and the movement in the respective direction.
                There are two phases:
                 - change direction with probability (see below)
                 - move in all directions
                The movement takes place according to the correlated random walk (CRW).

propagation tpye: AnisotropicPropagator
Description Propagator:
                The propagator is responsible for the movement of the particles. 
                It performs the change of direction and the movement in the respective direction.
                There are two phases:
                 - change direction with probability (see below)
                 - move in all directions
                The movement takes place according to the correlated random walk (CRW).

coordinates: Cylindrical coordinates
dimensions:  3
pitch angle: constant
m/s
speed:  299800000.0 m/s
number steps:  100000
step size:  1000000000.0 m
step duration:  3.335556983947754  s
total distance:  100000000000000.0 m
total duration:  333555.703802535  s
probability to change directions in step:  [1.923077   1.923077   0.03571429] %

Observer

The Observer determines during the simulation when and what data to write out (observe).

In each simulation step, the current particle state is evaluated by the Observer to check if one of the observing contions is satisfyed. The conditions to observe can be based on the time (-> step) or the coordinates of the particle.

The conditions to observe can be based on the time (or step) or the coordinates of the particle.

• step number [unique_steps] -> time (TimeEvolutionObservers)
• radius of observer sphere [shperes] -> sphere around source (SphericalObservers)
• cartesian coordinates [box_dimensions] -> box around source (BoxObserver) (not yet implemented)

All special observer will create an Observer object and specify the relevant parameters for the observation conditions (unique_steps, shperes, box_dimensions)

substeps = [False, False, True] # observe only steps (no substeps)
min_step = 1
max_step = nr_steps
nr_obs_steps = 600

observer = pp.TimeEvolutionObserverLog(min_step, max_step, nr_obs_steps, substeps)

sim.add_observer(observer)
sim.observer.get_description()

Observer initialized
Description Observer:
                The observer defines the conditions for when to write data to the output.

observer tpye: TimeEvolutionObserverLog
spheres:  [-1.]
steps [0:10]:  [ 1  2  3  4  5  6  7  8  9 10]
steps [-11:]:  [ 82513  84115  85747  87411  89107  90837  92600  94397  96228  98096
 100000]
nr steps:  445
substeps:  [False False  True]
all_steps:  False

%time sim.run_simulation()
sim.save_data('data/data_tut_2')

CPU times: user 3min 33s, sys: 644 ms, total: 3min 34s
Wall time: 3min 34s
saved data file data/data_tut_2.pkl

Analyze statistics

Diffusion coefficients

As the particles propagate via a random walk, statistical properties of many particles are interesting, such as the diffusion coefficients and particle distributions. These quantities can be compared to analytical predictions.

Running diffusion coefficients are computed with the mean-squared displacement method:

$\kappa_{ii}(t) = \frac{<\Delta x_i^2>}{2t}$.

df = pd.read_pickle("data/data_tut_2.pkl")
df_time_evolution_observer = df.loc[df['radius'] == -1.0]

sta = pp.Statistics(df_time_evolution_observer)
isotropic = False # diffusion is anisotropic
errors = False # don't show error bars
df_kappas = sta.plot_diffusion_coefficients(isotropic, errors)
print('input kappa_perp:', f"{float(diffusion_coefficient_perp):.3}", 'm²/s')
print('input kappa_para:', f"{float(diffusion_coefficient_para):.3}", 'm²/s')
n = 100
print('kappa_{xx}:', f"{np.mean(df_kappas['kappa_xx'][-n:]):.3}", 'm²/s', '+-', f"{np.std(df_kappas['kappa_xx'][-n:]):.3}", 'm²/s')
print('kappa_{yy}:', f"{np.mean(df_kappas['kappa_yy'][-n:]):.3}", 'm²/s', '+-', f"{np.std(df_kappas['kappa_yy'][-n:]):.3}", 'm²/s')
print('kappa_{zz}:', f"{np.mean(df_kappas['kappa_zz'][-n:]):.3}", 'm²/s', '+-', f"{np.std(df_kappas['kappa_zz'][-n:]):.3}", 'm²/s')

init statistics plotting class

diffusion coefficients computed between 1.49e+13m and 1.00e+14m with 100 data points
kappa_{xx}: 1.26e+18 m²/s +- 4.08e+16 m²/s
kappa_{yy}: 1.29e+18 m²/s +- 5.23e+16 m²/s
kappa_{zz}: 1.28e+20 m²/s +- 5.09e+18 m²/s
input kappa_perp: 1.3e+18 m²/s
input kappa_para: 1.4e+20 m²/s
kappa_{xx}: 1.26e+18 m²/s +- 4.08e+16 m²/s
kappa_{yy}: 1.29e+18 m²/s +- 5.23e+16 m²/s
kappa_{zz}: 1.28e+20 m²/s +- 5.09e+18 m²/s