The RT modules

mag_ray

This module does operations on ray_t objects.

Note that the casting of a ray itself is a part of the ray_vol_t objects.

A ray_t is and extension of a func_t object.

Importantly, it has two important quantities:

  • S_a and S_b: source function at the start and end of a ray segment

  • chi_a and chi_b: \(\kappa\rho\) at the start and end of a eay segment.

A ray_t object has the following methods:

  • tau: return piecewise constant optical depth \(\tau\) at the start of each segment, with \(\tau=0\) at the origin of the ray (ray_t%r_0). The resulting vector has ray_t%n +1 elements.

  • I: This wrapper method will act differently if it is given one or two parameters.

  • I_1(I_star): Returns \(I(tau=0)\) for the ray. I_star is the intensity at the back of the ray (so \(I_\star\) or 0)

  • I_n(I_star, profile): Returns \(I_\lambda(tau=0)\). ‘profile’ is the photospheric profile as a function of \(\lambda\). Note that the source function is grey in this calculation.

  • I_X: This is something added by Vero (needs a better implementation) to handle the Xray transfer. BEWARE: It assumes that the quantity named S is in fact \(\eta \rho_h\).

  • I_n: for halpha line profiles – more on this later on.

More background information about the I_X implementation in mag_ray

In principle, we can create a source function for X-ray radiative transfer. After all:

\[ \frac{dI}{ds} = -I \kappa\rho + j\rho \]

which, using \(d\tau=-\kappa \rho ds\) can be recast as:

\[ \frac{dI}{\kappa\rho ds} = -I \frac{\kappa\rho}{\kappa\rho} + \frac{j\rho}{\kappa\rho} \]

In the classical picture where each grid cell in the model contains only a single density (at a single temperature and thus has a single \(\kappa\) and \(\rho\)), we can simply cancel out the terms in the right-hand side, and define \(S=j/\kappa\).

However, in the case of e.g. the adm model, this was problematic. The hot component is conceptually separated from the cold component (upflow and downflow). This could be fixed by setting:

\[ \frac{dI}{\kappa\rho_c ds} = -I \frac{\kappa\rho_c}{\kappa_c\rho} + \frac{j\rho_h}{\kappa_c\rho} \]

And the the source function could be described as \(\frac{j\rho_h}{\kappa\rho_c}\).

The problem here is that this is not numerically stable. In the limit where \(\kappa_c\) goes to zero, the source function calculated this way blows up. A work around is to use a taylor expension.

The formal solution for piece-wise constant source function is:

\[I_a = I_b \mathrm{e}^{\tau_a - \tau_b} + S[1 - \mathrm{e}^{\tau_a-\tau_b}]\]

Noting that \(\tau_a - \tau_b = -\kappa\rho_c ds = \chi ds\) and defining \(\eta = j\rho_h\), we have:

\[I_a = I_b \mathrm{e}^{\tau_a - \tau_b} + \eta ds \left[ \frac{1 - \mathrm{e}^{\tau_a-\tau_b}}{\chi ds}\right]\]

The left term in bracket is analytically 1.0 when \(\chi ds\) goes to zero, but will not be numerically stable.

However, noting that

\[\frac{1-\mathrm{e}^x}{x} \approx 1 - \frac{x}{2} + \frac{x^2}{6} - \frac{x^3}{24} + ...\]

and given the fact that using the 4 first terms of this expension is rather good for \(x< 1\), Vero made a separate routine do calculate the intensity by doing:

\[I_a = I_b \mathrm{e}^{\tau_a - \tau_b} + \eta ds \mathcal{f}(\chi ds)\]

where \(\mathcal{f}\) is the exact expression with source function if \((\chi ds) > 1\) and the taylor approximation if \((\chi ds)< 1\).

Practically, we are only interested in \(I(\tau=0)\). Thus we can simplify piece-wise calculation analytically first. Starting at the back of the ray, we can analitically write down the intensity at the end of the first segment. Then we can use this expression as the intensity entering in the next segment, etc. Doing this n-times, simplifies to:

\[I(\tau=0) = I_\star \mathrm{e}^{-tau_{n+1}} + \sum_i \eta_i ds_i \mathcal{f}(\chi_i ds_i) \mathrm{e}^{-\tau_i}\]

This implies that the values S and \(chi\) in the ray_t definition do not have their original values to use this function. The way it is coded right now, it assumes that the quantity named S is in fact \(\eta \rho_h\). This is not great coding. We should create a new object that is a modified ray for which the name of the quantity reflects what they actually are.

mag_ray_vol

This module does operations on ray_vol_t.

A ray_vol_t is an expension of a volume_t object. In addition to all the volume_t information and operation, a ray_vol_t object contains a star_t object, and two 3D arrays:

  • S: The value of the source function at all grid points (BEWARE – different for the Xrays)

  • Chi: kappa*rho at all grid points.

The module also defined the following operations:

  • cast_ray: is the procedure that will trace a ray through the ray_vol_t, resulting in a ray_t object, whose values of S and chi have been interpolated through the volume. This is used in the programs that e.g. builds optical depth maps – the rays are cast through the ray_vol_t at a certain viewing angle in a loop, and for each rays the ray_t tau procedure is used to get the optical depth of the rays.

  • eval_S: returns the value of S at a given location.

  • eval_chi: returns the value of chi at a given location.

Creating ray_volumes objects

There are a few programs that will create ray_volumes for various situations.

  • build_ray_volume: generic program that will create a ray_volume for different inputs:

  • ARRM: from a arrm_volume file (TODO: list the options)

  • ADM: from a adm_volume file – although there is a separate program for this (see below), which I think was created afterward with more options (and this is the one in which I added the Xrays)

  • ETA: I cannot remember what an eta_volume is on the top of my head – need to check.

  • build_ray_volume_adm_vero: program that has more options for building a ray_volume from an adm file. Vero added a specific option to build a ray_vol where the quantity named “S” has the x-ray emmissivity in it. This is not great coding – we should do this properly later on (by adding a ray_t that has an appropiately named value in it.)

  • make_ray_vol: this is a program that builds a ray_volume from scratch, with chi and S being a pair of spheres at x=+/- 3 – this is useful to test the orientation of the coordinate systems in the following program that performs the RT using ray_volumes.

Creating maps

There are some programs that will create a maps in the plane of the sky. For all of these program, the inclination of the rotational axis with respect to the line of sight can be specified, along with the number of desired rotational phases.

The resulting file is a set of sky maps.

  • calc_tau_map: makes a map of optical depth

  • calc_emission_map: makes a map of “intensity” (result from I_1 in the ray_t procedure)

  • calc_Xray_map: makes a map of “intensity” using the result from I_X in the ray_t procedure.

workflow

This is a good place to explain how all of the classes are used together in these map programs.

  1. The ray_volume object is read from a file.

  2. A xy map of the sky is made, with defines the origin of the rays in the LOS coordinate system

  3. In a loop over sky xy:

    1. The direction of the ray in the ray_volume object is determined (with a coordinate transformation between LOS and the rotational coordinate system, using the inclination of the rotation axis and the rotational phase)

    2. We cast a ray_t object through the ray_volume (this is a procedure of the ray_volume class, which determines the values of S and chi along the ray by interpolating through the ray_volume)

    3. We get the value of the X-ray intensity through the I_X procedure of the ray_t class.

  4. Repeat for other rotational phases, and save the maps to file