Run the Fleischman & Kuznetsov (2010) synchrotron code (pwkit.fk10)

This module helps you run the synchrotron code described in Fleischman & Kuznetsov (2010; hereafter FK10) [DOI:10.1088/0004-637X/721/2/1127]. The code is provided as a precompiled binary module. It’s meant to be called from IDL, but we won’t let that stop us!

The main interface to the code is the Calculator class. But before you can use it, you must install the code, as described below.

Installing the code

To do anything useful with this module, you must first obtain the precompiled module. This isn’t the sort of module you’d want to install into your system shared libraries, so for applications you’ll probably be storing it in some random directory. Therefore, all actions in this module start by specifying the path to the library.

The module can be downloaded from as part of a Supplementary Data archive attached to the journal paper. At the moment, the direct link is here, but that might change over time. The journal’s website for the paper should always have a link.

The archive contains compiled versions of the code for Windows, 32-bit Linux, and 64-bit Linux. It is quite worrisome that maybe one day these files will stop working, but that’s what we’ve got right now.

Anyway, you should download and unpack the archive and copy the desired file to wherever makes the most sense for your software environment and application. On 64-bit Linux, the file name is libGS_Std_HomSrc_CEH.so.64. Any variable named shlib_path that comes up in the API should be a path to this file. Note that relative paths should include a directory component (e.g. ./libGS_Std_HomSrc_CEH.so.64); the ctypes module treats bare filenames specially.

class pwkit.fk10.Calculator(shlib_path)

An interface to the FK10 synchrotron routines.

This class maintains state about the input parameters that can be passed to the routines, and can invoke them for you.

Constructor arguments

shlib_path

The path to the compiled FK10 code, as described in the module-level documentation.

Newly-constructed objects are initialized with:

self.set_hybrid_parameters(12, 12)
self.set_ignore_q_terms(False)
self.set_trapezoidal_integration(15)

Setting parameters

set_bfield(B_G)	Set the strength of the local magnetic field.
set_bfield_for_s0(s0)	Set B to probe a certain harmonic number.
set_edist_powerlaw(emin_mev, emax_mev, ...)	Set the energy distribution function to a power law.
set_edist_powerlaw_gamma(gmin, gmax, delta, ...)	Set the energy distribution function to a power law in the Lorentz factor
set_freqs(n, f_lo_ghz, f_hi_ghz)	Set the frequency grid on which to perform the calculations.
set_hybrid_parameters(s_C, s_WH[, do_renorm])	Set the hybrid/renormalization control parameters.
set_ignore_q_terms(ignore_q_terms)	Set whether "Q" terms are ignored.
set_obs_angle(theta_rad)	Set the observer angle relative to the field.
set_one_freq(f_ghz)	Set the code to calculate results at just one frequency.
set_padist_gaussian_loss_cone(boundary_rad, ...)	Set the pitch-angle distribution to a Gaussian loss cone.
set_padist_isotropic()	Set the pitch-angle distribution to be isotropic.
set_thermal_background(T_K, nth_cc)	Set the properties of the background thermal plasma.
set_trapezoidal_integration(n)	Set the code to use trapezoidal integration.

Running calculations

find_rt_coefficients([depth0])	Figure out emission and absorption coefficients for the current parameters.
find_rt_coefficients_tot_intens([depth0])	Figure out total-intensity emission and absorption coefficients for the current parameters.

set_bfield(B_G)

Set the strength of the local magnetic field.

Call signature

B_G

The magnetic field strength, in Gauss

Returns

self for convenience in chaining.

set_bfield_for_s0(s0)

Set B to probe a certain harmonic number.

Call signature

s0

The harmonic number to probe at the lowest frequency

Returns

self for convenience in chaining.

This just proceeds from the relation nu = s nu_c = s e B / 2 pi m_e c. Since s and nu scale with each other, if multiple frequencies are being probed, the harmonic numbers being probed will scale in the same way.

set_edist_powerlaw(emin_mev, emax_mev, delta, ne_cc)

Set the energy distribution function to a power law.

Call signature

emin_mev

The minimum energy of the distribution, in MeV

emax_mev

The maximum energy of the distribution, in MeV

delta

The power-law index of the distribution

ne_cc

The number density of energetic electrons, in cm^-3.

Returns

self for convenience in chaining.

set_edist_powerlaw_gamma(gmin, gmax, delta, ne_cc)

Set the energy distribution function to a power law in the Lorentz factor

Call signature

gmin

The minimum Lorentz factor of the distribution

gmax

The maximum Lorentz factor of the distribution

delta

The power-law index of the distribution

ne_cc

The number density of energetic electrons, in cm^-3.

Returns

self for convenience in chaining.

set_freqs(n, f_lo_ghz, f_hi_ghz)

Set the frequency grid on which to perform the calculations.

Call signature

n

The number of frequency points to sample.

f_lo_ghz

The lowest frequency to sample, in GHz.

f_hi_ghz

The highest frequency to sample, in GHz.

Returns

self for convenience in chaining.

set_hybrid_parameters(s_C, s_WH, do_renorm=True)

Set the hybrid/renormalization control parameters.

Call signature

s_C

The harmonic number above which the continuous approximation is used (with special behavior; see below).

s_WH

The harmonic number above which the Wild-Hill BEssel function approximations are used.

do_renorm (default True)

Whether to do any renormalization at all.

Returns

self for convenience in chaining.

FK10 uses frequency parameters f^C_cr and f^WH_cr to control some of its optimizations. This function sets these parameters as multiples of the electron cyclotron frequency (f_Be in FK10 notation): e.g., f^C_cr = s_C * f_Be.

At frequencies above f^C_cr, the “continuum” approximation is introduced, replacing the “exact” sum with an integral. At frequencies above f^WH_cr, the Wild-Hild approximations to the Bessel functions are used. In both cases, the activation of the optimizations can result in normalization shifts in the calculations. “Renormalization” computes these shifts (by doing both kinds of calculations at the transition frequencies) and attempts to correct them. (Some of the FK10 documentation seems to refer to renormalization as “R-optimization”.)

If f^C_cr is below the lowest frequency integrated, all calculations will be done in continuum mode. In this case, the sign of s_C sets whether Wild-Hill renormalization is applied. If s_C is negative and f^WH_cr is above the lowest frequency integration, renormalization is done. Otherwise, it is not.

The documentation regarding f^WH_cr is confusing. It states that f^WH_cr only matters if (1) s_WH < s_C or (2) s_C < 0 and f^WH_cr > f_0. It is not obvious to me why s_WH > s_C should only matter if s_C < 0, but that’s what’s implied.

In most examples in FK10, both of these parameters are set to 12.

set_ignore_q_terms(ignore_q_terms)

Set whether “Q” terms are ignored.

Call signature

ignore_q_terms

If true, ignore “Q” terms in the integrations.

Returns

self for convenience in chaining.

See Section 4.3 of FK10 and OnlineII.pdf in the Supplementary Data for a precise explanation. The default is to not ignore the terms.

set_obs_angle(theta_rad)

Set the observer angle relative to the field.

Call signature

theta_rad

The angle between the ray path and the local magnetic field, in radians.

Returns

self for convenience in chaining.

set_one_freq(f_ghz)

Set the code to calculate results at just one frequency.

Call signature

f_ghz

The frequency to sample, in GHz.

Returns

self for convenience in chaining.

set_padist_gaussian_loss_cone(boundary_rad, expwidth)

Set the pitch-angle distribution to a Gaussian loss cone.

Call signature

boundary_rad

The angle inside which there are no losses, in radians.

expwidth

The characteristic width of the Gaussian loss profile in direction-cosine units.

Returns

self for convenience in chaining.

See OnlineI.pdf in the Supplementary Data for a precise definition. (And note the distinction between α_c and μ_c since not everything is direction cosines.)

set_padist_isotropic()

Set the pitch-angle distribution to be isotropic.

Returns

self for convenience in chaining.

set_thermal_background(T_K, nth_cc)

Set the properties of the background thermal plasma.

Call signature

T_K

The temperature of the background plasma, in Kelvin.

nth_cc

The number density of thermal electrons, in cm^-3.

Returns

self for convenience in chaining.

Note that the parameters set here are the same as the ones that describe the thermal electron distribution, if you choose one of the electron energy distributions that explicitly models a thermal component (“thm”, “tnt”, “tnp”, “tng”, “kappa” in the code’s terminology). For the power-law-y electron distributions, these parameters are used to calculate dispersion parameters (e.g. refractive indices) and a free-free contribution, but their synchrotron contribution is ignored.

set_trapezoidal_integration(n)

Set the code to use trapezoidal integration.

Call signature

n

Use this many nodes

Returns

self for convenience in chaining.

find_rt_coefficients(depth0=None)

Figure out emission and absorption coefficients for the current parameters.

Argument

depth0 (default None)

A first guess to use for a good integration depth, in cm. If None, the most recent value is used.

Return value

A tuple (j_O, alpha_O, j_X, alpha_X), where:

j_O

The O-mode emission coefficient, in erg/s/cm^3/Hz/sr.

alpha_O

The O-mode absorption coefficient, in cm^-1.

j_X

The X-mode emission coefficient, in erg/s/cm^3/Hz/sr.

alpha_X

The X-mode absorption coefficient, in cm^-1.

The main outputs of the FK10 code are intensities and “damping factors” describing a radiative transfer integration of the emission from a homogeneous source. But there are times when we’d rather just know what the actual emission and absorption coefficients are. These can be backed out from the FK10 outputs, but only if the “damping factor” takes on an intermediate value not extremely close to either 0 or 1. Unfortunately, there’s no way for us to know a priori what choice of the “depth” parameter will yield a nice value for the damping factor. This routine automatically figures one out, by repeatedly running the calculation.

To keep things simple, this routine requires that you only be solving for coefficients for one frequency at a time (e.g., set_one_freq()).

find_rt_coefficients_tot_intens(depth0=None)

Figure out total-intensity emission and absorption coefficients for the current parameters.

Argument

depth0 (default None)

A first guess to use for a good integration depth, in cm. If None, the most recent value is used.

Return value

A tuple (j_I, alpha_I), where:

j_I

The total intensity emission coefficient, in erg/s/cm^3/Hz/sr.

alpha_I

The total intensity absorption coefficient, in cm^-1.

See find_rt_coefficients() for an explanation how this routine works. This version merely postprocesses the results from that method to convert the coefficients to refer to total intensity.