Modeling (without statistics)

In this tutorial we build a star object for signal computation, but streamline by omitting any demonstration of statistical modeling (for which there exist other tutorials).

%matplotlib inline

import os
import numpy as np
import math
import time

from matplotlib import pyplot as plt
from matplotlib import rcParams
from matplotlib.ticker import MultipleLocator, AutoLocator, AutoMinorLocator
from matplotlib import gridspec
from matplotlib import cm
rcParams['text.usetex'] = False
rcParams['font.size'] = 14.0

import xpsi

from xpsi.global_imports import _c, _G, _dpr, gravradius, _csq, _km, _2pi
from xpsi.utilities import PlottingLibrary as XpsiPlot

/=============================================\
| X-PSI: X-ray Pulse Simulation and Inference |
|---------------------------------------------|
|                Version: 3.0.0               |
|---------------------------------------------|
|      https://xpsi-group.github.io/xpsi      |
\=============================================/

Imported emcee version: 3.1.6
Imported PyMultiNest.
Imported UltraNest.
Imported GetDist version: 1.5.3
Imported nestcheck version: 0.2.1

Constructing a star

We need to build our star. The basic units for building a star are:

• the Spacetime class;

• the Photosphere class;

• the HotRegion class;

• the Elsewhere class (optionally used in conjuction with HotRegion instances);

• the Everywhere class (cannot be used with HotRegion instances);

• and four low-level user-modifiable routines for evaluation of a parametrised specific intensity model.

For this demonstration we will assume that the surface radiation field elsewhere (other than the hot regions) can be ignored in the soft X-ray regime our model instrument is sensitive to. We will not be utilizing the Elsewhere class, for which their exists a distinct tutorial (Global surface emission). For more advanced modelling, we can simply write custom derived classes, and instantiate those derived classes to construct objects for our model. In particular, a common pattern will be to subclass the HotRegion class. Let’s start with the Spacetime class.

The ambient spacetime

xpsi.Spacetime#? # uncomment to query

xpsi.Spacetime.Spacetime

We can instantiate the spacetime by defining all parameters with a given value. In this case, all parameters will be fixed because the bounds are not specified (empty dictionary). Note that all parameters must be defined at least once in ``bounds`` or ``values``.

values = dict(distance =0.150,          # (Earth) distance
                mass = 1.5,             # mass
                radius = 12.0,          # equatorial radius
                cos_inclination = 0.5,  # (Earth) inclination to rotation axis
                frequency = 300. )      # spin frequency

spacetime = xpsi.Spacetime(bounds={}, values=values)

Creating parameter:
    > Named "frequency" with fixed value 3.000e+02.
    > Spin frequency [Hz].
Creating parameter:
    > Named "mass" with fixed value 1.500e+00.
    > Gravitational mass [solar masses].
Creating parameter:
    > Named "radius" with fixed value 1.200e+01.
    > Coordinate equatorial radius [km].
Creating parameter:
    > Named "distance" with fixed value 1.500e-01.
    > Earth distance [kpc].
Creating parameter:
    > Named "cos_inclination" with fixed value 5.000e-01.
    > Cosine of Earth inclination to rotation axis.

Alternatively we can specify bounds manually for the free parameters. Using (None, None) will set the default bounds. The fixed parameters are defined in the values dictionary. If both the bounds and a value are defined for a parameters, the value will not be fixed but instead will be set as an initial value.

bounds = dict(distance = (None, None),                   # Default bounds for (Earth) distance
                mass = (None, None),                     # Default bounds for mass
                radius = (None, None),                   # Default bounds for equatorial radius
                cos_inclination = (None, None))          # Default bounds for (Earth) inclination to rotation axis

values = dict(frequency=300.0,                           # Fixed spin frequency
              mass=1.4)                                  # mass with initial value

spacetime = xpsi.Spacetime(bounds=bounds, values=values)

Creating parameter:
    > Named "frequency" with fixed value 3.000e+02.
    > Spin frequency [Hz].
Creating parameter:
    > Named "mass" with bounds [1.000e-03, 3.000e+00] and initial value 1.400e+00.
    > Gravitational mass [solar masses].
Creating parameter:
    > Named "radius" with bounds [1.000e+00, 2.000e+01].
    > Coordinate equatorial radius [km].
Creating parameter:
    > Named "distance" with bounds [1.000e-02, 3.000e+01].
    > Earth distance [kpc].
Creating parameter:
    > Named "cos_inclination" with bounds [-1.000e+00, 1.000e+00].
    > Cosine of Earth inclination to rotation axis.

To display the free parameters:

for p in spacetime:
    print(p)

Gravitational mass [solar masses].
Coordinate equatorial radius [km].
Earth distance [kpc].
Cosine of Earth inclination to rotation axis.

For the following example in this tutorial, we will specify the bounds of all parameters.

bounds = dict(distance = (0.1, 1.0),                     # (Earth) distance
                mass = (1.0, 3.0),                       # mass with default bounds
                radius = (3.0 * gravradius(1.0), 16.0),  # equatorial radius
                cos_inclination = (0.0, 1.0))            # (Earth) inclination to rotation axis

values = dict(frequency=300.0 )                          # Fixed spin frequency

spacetime = xpsi.Spacetime(bounds=bounds, values=values)

Creating parameter:
    > Named "frequency" with fixed value 3.000e+02.
    > Spin frequency [Hz].
Creating parameter:
    > Named "mass" with bounds [1.000e+00, 3.000e+00].
    > Gravitational mass [solar masses].
Creating parameter:
    > Named "radius" with bounds [4.430e+00, 1.600e+01].
    > Coordinate equatorial radius [km].
Creating parameter:
    > Named "distance" with bounds [1.000e-01, 1.000e+00].
    > Earth distance [kpc].
Creating parameter:
    > Named "cos_inclination" with bounds [0.000e+00, 1.000e+00].
    > Cosine of Earth inclination to rotation axis.

The photosphere and its constituent regions

It is not necessary for us to write a custom derived class for the photosphere object, so we will simply instantiate a Photosphere object. However, we first need an instance of HotRegion to instantiate the photosphere, and we need to implement a low-level parametrised model for the specific intensity emergent from the photosphere in a local comoving frame.

The neutron star atmosphere is assumed to be geometrically thin. In the applications thus far, the surface local-comoving-frame radiation field as being described by a single free parameter: the effective temperature. The radiation field is also dependent on the local effective gravitational acceleration, however this is a derived parameter in the model. The parametrised radiation field as a function of energy and angle subtended to the normal to the (plane-parallel) atmosphere in a local comoving frame is provided as numerical model data for multi-dimensional interpolation.

In X-PSI, integration over the surface radiation field is performed via calls to low-level C routines. To reduce likelihood evaluation times, the atmosphere interpolator is written in C, and calls to that interpolator are from C routine. In other words, in X-PSI, we do not use Python callback functions for evaluation of specific intensities, but C functions which are compiled when the X-PSI package is built. Unfortunately this means that to change the parametrised surface radiation field you need to get your hands a little dirty; on the bright side, the body of these functions can be implemented almost completely in the Cython language, so syntactically there is some similarity to Python because the language syntax is somewhat of a hybrid/superset. Beware, however, that the body of these functions must not contain calls to the Python API, and only to external C libraries if required: the code must evaluate to pure C, and not require the Python/C API. Note that the Python global interpreter lock is deactivated during integration to enable OpenMP multi-threading in some applications of the integrator; thus the code needs to be thread safe and nogil (not require the global interpreter lock, although a context manager could in principle be used to reacquire the lock within the integrator). Also note that if external C libraries are required, that you include a Cython .pxd (header) file in the package which externs the required library components; the library also needs to be included and linked in setup.py at package build-time.

If you have ideas for making this model specification more user-friendly, without, crucially, increasing signal integration time, please submit a pull request.

There are several atmosphere extension modules for users and developers to work with; hot_BB.pyx, and hot_Num4D.pyx are the built-in options for isotropic blackbody radiation field and numerical 4D-interpolation from a preloaded atmosphere table, and they can be selected with the atm_ext option when creating instances of the radiating regions (from HotRegion, Elsewhere, or Everywhere classes); hot_user.pyx and elsewhere_user.pyx are additional extensions which can be replaced by user-developed versions, and used after re-installing X-PSI and setting atm_ext="user" (by default they are the same as hot_BB.pyx).

Note that in X-PSI versions before 2.1.0 the selection of the atmosphere extension needed to be done when installing X-PSI (using installation flags), and not using the atm_ext option as now.

The following is the contents of the hot_user.pxd file which the X-PSI integrators use as the header file for including other C functions in the package (if using the non-modified “user” option).

from .preload cimport _preloaded

cdef double eval_hot_user(size_t THREAD,
                     double E,
                     double mu,
                     const double *const VEC,
                     void *const data) nogil

cdef double eval_hot_norm_user() nogil

cdef void* init_hot_user(size_t numThreads, const _preloaded *const preloaded) nogil

cdef int free_hot_user(size_t numThreads, void *const data) nogil

You are free to modify these functions in the associated hot_user.pyx implementation file, and you have almost complete control over the function bodies, but not the signatures. By default the package includes an isotropic blackbody model:

#cython: cdivision=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False

from libc.math cimport exp, pow
from libc.stdio cimport printf

from xpsi.global_imports import _keV, _k_B

cdef int SUCCESS = 0

cdef double erg = 1.0e-7
cdef double Planck_dist_const = 5.040366110812353e22

cdef double k_B = _k_B
cdef double keV = _keV
cdef double k_B_over_keV = k_B / keV

#----------------------------------------------------------------------->>>
# >>> User modifiable functions.
# >>> Note that the user is entirely free to wrap thread-safe and
# ... non-parallel external C routines from an external library.
# >>> Thus the bodies of the following need not be written explicitly in
# ... the Cython language.
#----------------------------------------------------------------------->>>
cdef void* init_hot_user(size_t numThreads, const _preloaded *const preloaded) nogil:
    # This function must match the free management routine free_hot()
    # in terms of freeing dynamically allocated memory. This is entirely
    # the user's responsibility to manage.

    if preloaded != NULL :
        printf("WARNING: Numerical atmosphere data were preloaded, even though those are not used by this atmosphere extension.\n")

    # Return NULL if dynamic memory is not required for the model.
    return NULL

cdef int free_hot_user(size_t numThreads, void *const data) nogil:
    # This function must match the initialisation routine init_hot()
    # in terms of freeing dynamically allocated memory. This is entirely
    # the user's responsibility to manage.
    # The void pointer must be appropriately cast before memory is freed --
    # only the user can know this at compile time.
    # Just use free(<void*> data) iff no memory was dynamically
    # allocated in the function:
    #   init_local_hot()
    # because data is expected to be NULL in this case

    #printf("\nNo data to be freed.")

    return SUCCESS

cdef double eval_hot_user(size_t THREAD,
                     double E,
                     double mu,
                     const double *const VEC,
                     void *const data) nogil:
    # Arguments:
    # E = photon energy in keV
    # mu = cosine of ray zenith angle (i.e., angle to surface normal)
    # VEC = variables such as temperature, effective gravity, ...
    # data = numerical model data required for intensity evaluation

    cdef double temp = k_B_over_keV * pow(10.0, VEC[0])

    return E * E * E / ( exp(E / temp) - 1.0 )

cdef double eval_hot_norm_user() nogil:
    # Source radiation field normalisation which is independent of the
    # parameters of the parametrised model -- i.e. cell properties, energy,
    # and angle.
    # Writing the normalisation here reduces the number of operations required
    # during integration.
    # The units of the specific intensity need to be J/cm^2/s/keV/steradian.

    return erg * Planck_dist_const

In most use-cases we apply the hot_Num4D.pyx extension to handle the numerical atmosphere data. This extension is presented in surface radiation field extension module, and is essentially the same as used by Riley et al. (2019) (and the follow-up papers) to implement a numerical atmosphere generated by the NSX atmosphere code (Ho & Lai (2001); Ho & Heinke (2009)), courtesy of W.C.G. Ho for NICER modeling efforts. A fully-ionized hydrogen atmosphere was used in Riley et al. (2019); also see Bogdanov et al. (2021) and Riley et al. (2021).

We now instantiate hot region objects. We can find the required and optional parameter names in the class docstring:

xpsi.HotRegion#? # uncomment to query

xpsi.HotRegion.HotRegion

The names can also be found as class attributes as follows (refer to the HotRegion class docstring and Riley et al. (2019) for the meaning of super):

xpsi.HotRegion.required_names

['super_colatitude',
 'super_radius',
 'phase_shift',
 'super_temperature (if no custom specification)']

The condition if no custom specification means that this name is required if we do not supply custom parameters for the radiation field in the superseding member of the hot region. If we supply custom parameters, we also need to subclass xpsi.HotRegion and overwrite the __compute_cellParamVecs method to handle our parameters.

xpsi.HotRegion.optional_names

['omit_colatitude',
 'omit_radius',
 'omit_azimuth',
 'cede_colatitude',
 'cede_radius',
 'cede_azimuth',
 'cede_temperature']

For the purpose of illustration, we tie the temperatures of the hot regions together. There is more than one way to achieve this, but we will use the most powerful approach.

bounds = dict(super_colatitude = (None, None),
              super_radius = (None, None),
              phase_shift = (0.0, 0.1),
              super_temperature = (None, None))

# a simple circular, simply-connected spot
primary = xpsi.HotRegion(bounds=bounds,
                            values={}, # no initial values and no derived/fixed
                            symmetry=True,
                            omit=False,
                            cede=False,
                            concentric=False,
                            sqrt_num_cells=32,
                            min_sqrt_num_cells=10,
                            max_sqrt_num_cells=64,
                            num_leaves=100,
                            num_rays=200,
                            do_fast=False,
                            atm_ext="BB",#default blackbody, other options: "Num4D" or "user"
                            prefix='p') # unique prefix needed because >1 instance

Creating parameter:
    > Named "super_colatitude" with bounds [0.000e+00, 3.142e+00].
    > The colatitude of the centre of the superseding region [radians].
Creating parameter:
    > Named "super_radius" with bounds [0.000e+00, 1.571e+00].
    > The angular radius of the (circular) superseding region [radians].
Creating parameter:
    > Named "phase_shift" with bounds [0.000e+00, 1.000e-01].
    > The phase of the hot region, a periodic parameter [cycles].
Creating parameter:
    > Named "super_temperature" with bounds [3.000e+00, 7.600e+00].
    > log10(superseding region effective temperature [K]).

Note that since the atmospheric local-comoving-frame effective temperature is uniform everywhere within the hot region boundaries, we can use the default value of the symmetry keyword, True. All other arguments determine the numerical resolution, and have defaults which have been (somewhat arbitrarily) chosen to be result in a likelihood evaluation time of \(\mathcal{O}(1)\) s.

Let’s take a look at the xpsi.Derive docstring for guidance:

xpsi.Derive#? # uncomment to query

xpsi.Parameter.Derive

class derive(xpsi.Derive):
    def __init__(self):
        """
        We can pass a reference to the primary here instead
        and store it as an attribute if there is risk of
        the global variable changing.

        This callable can for this simple case also be
        achieved merely with a function instead of a magic
        method associated with a class.
        """
        pass

    def __call__(self, boundto, caller = None):
        # one way to get the required reference
        global primary # unnecessary, but for clarity
        return primary['super_temperature'] - 0.2

bounds['super_temperature'] = None # declare fixed/derived variable

secondary = xpsi.HotRegion(bounds=bounds, # can otherwise use same bounds
                              values={'super_temperature': derive()}, # create a callable value
                              symmetry=True,
                              omit=False,
                              cede=False,
                              concentric=False,
                              sqrt_num_cells=32,
                              min_sqrt_num_cells=10,
                              max_sqrt_num_cells=100,
                              num_leaves=100,
                              num_rays=200,
                              do_fast=False,
                              is_antiphased=True,
                              atm_ext="BB",
                              prefix='s') # unique prefix needed because >1 instance

Creating parameter:
    > Named "super_colatitude" with bounds [0.000e+00, 3.142e+00].
    > The colatitude of the centre of the superseding region [radians].
Creating parameter:
    > Named "super_radius" with bounds [0.000e+00, 1.571e+00].
    > The angular radius of the (circular) superseding region [radians].
Creating parameter:
    > Named "phase_shift" with bounds [0.000e+00, 1.000e-01].
    > The phase of the hot region, a periodic parameter [cycles].
Creating parameter:
    > Named "super_temperature" that is derived from ulterior variables.
    > log10(superseding region effective temperature [K]).

The description derived from ulterior variables means that when we lookup the value, it is calculated dynamically from the values of other (ulterior) model parameters. We clearly expect the temperature of the secondary hot region to behave in this way. A few other varibles do to because of keyword arguments passed upon instantiation of the hot regions. For example, the colatitude of the zero-radii omission and ceding regions (omit=False and cede=False) are equivalent to the colatitude of the centre of the superseding region. The azimuths are relative to the superseding region, and are thus listed as being fixed at zero azimuthal separation. If one of omit or cede was True, and concentric=True, a similar setup is performed, but with the radius of omit or cede being free, fixed (at finite value, but zero achieves the same as False for both omit and cede keyword arguments), or derived.

We now need to encapsulate the hot region instances in a container with properties expected by the Photosphere class.

from xpsi import HotRegions

hot = HotRegions((primary, secondary))

Let’s check out the hot regions:

hot.objects[0] # 'p'

Free parameters
---------------
p__phase_shift: The phase of the hot region, a periodic parameter [cycles].
p__super_colatitude: The colatitude of the centre of the superseding region [radians].
p__super_radius: The angular radius of the (circular) superseding region [radians].
p__super_temperature: log10(superseding region effective temperature [K]).

Derived/fixed parameters
------------------------
p__cede_colatitude: The colatitude of the centre of the ceding region [radians].
p__cede_radius: The angular radius of the (circular) ceding region [radians].
p__cede_azimuth: The azimuth of the centre of the ceding region relative to the
centre of the superseding region [radians].
p__omit_colatitude: The colatitude of the centre of the omission region [radians].
p__omit_radius: The angular radius of the (circular) omission region [radians].
p__omit_azimuth: The azimuth of the centre of the omission region relative to the
centre of the superseding region [radians].

hot.objects[1] # 's'

Free parameters
---------------
s__phase_shift: The phase of the hot region, a periodic parameter [cycles].
s__super_colatitude: The colatitude of the centre of the superseding region [radians].
s__super_radius: The angular radius of the (circular) superseding region [radians].

Derived/fixed parameters
------------------------
s__super_temperature: log10(superseding region effective temperature [K]).
s__cede_colatitude: The colatitude of the centre of the ceding region [radians].
s__cede_radius: The angular radius of the (circular) ceding region [radians].
s__cede_azimuth: The azimuth of the centre of the ceding region relative to the
centre of the superseding region [radians].
s__omit_colatitude: The colatitude of the centre of the omission region [radians].
s__omit_radius: The angular radius of the (circular) omission region [radians].
s__omit_azimuth: The azimuth of the centre of the omission region relative to the
centre of the superseding region [radians].

A list of names, with the prefix, can also be accessed as follows:

h = hot.objects[0]
h.names

['p__phase_shift',
 'p__super_colatitude',
 'p__super_radius',
 'p__super_temperature',
 'p__cede_colatitude',
 'p__cede_radius',
 'p__cede_azimuth',
 'p__omit_colatitude',
 'p__omit_radius',
 'p__omit_azimuth']

h.prefix

'p'

h.get_param('phase_shift')

The phase of the hot region, a periodic parameter [cycles]

Let’s set a value for the temperature of the primary hot region:

hot['p__super_temperature'] = 6.0 # equivalent to ``primary['super_temperature'] = 6.0``

Now let’s lookup the temperature of the secondary:

secondary['super_temperature']

5.8

No value was set explicitly for this secondary temperature: it is looked up dynamically from that of the primary hot region.

Note that the following access will not work:

# hot['s__super_temperature']

The reason for this is because the temperature of the secondary hot region is a fixed/derived variable. Only free model parameters are merged into larger spaces. A fixed/derived variable needs to be accessed via the subspace that directly encapsulates a reference to it.

We can now instantitate the photosphere (assuming a blackbody atmosphere, see https://github.com/xpsi-group/xpsi/tree/main/examples/examples_modeling_tutorial/TestRun_Num.py for an example when using numerical atmosphere):

class CustomPhotosphere(xpsi.Photosphere):
    """ Implement method for imaging."""

    @property
    def global_variables(self):
        """ This method is needed if we also want to invoke the image-plane signal simulator. """

        return np.array([self['p__super_colatitude'],
                          self['p__phase_shift'] * 2.0 * math.pi,
                          self['p__super_radius'],
                          0.0, #self['p__cede_colatitude'],
                          0.0, #self['p__phase_shift'] * 2.0 * math.pi - self['p__cede_azimuth'],
                          0.0, #self['p__cede_radius'],
                          self['s__super_colatitude'],
                          (self['s__phase_shift'] + 0.5) * 2.0 * math.pi,
                          self['s__super_radius'],
                          0.0, #self['s__cede_colatitude'],
                          0.0, #(self['s__phase_shift'] + 0.5) * 2.0 * math.pi - self['s__cede_azimuth'],
                          0.0, #self['s__cede_radius'],
                          self['p__super_temperature'],
                          0.0, #self['p__cede_temperature'],
                          self.hot.objects[1]['s__super_temperature'],
                          0.0]) #self['s__cede_temperature']])

photosphere = CustomPhotosphere(hot = hot, elsewhere = None,
                                values=dict(mode_frequency = spacetime['frequency']))

Creating parameter:
    > Named "mode_frequency" with fixed value 3.000e+02.
    > Coordinate frequency of the mode of radiative asymmetry in the
photosphere that is assumed to generate the pulsed signal [Hz].

photosphere['mode_frequency'] == spacetime['frequency']

True

We do not model the surface radiation field elsewhere, and we thus leave the elsewhere keyword as None (the default). Elsewhere means on the surface, exterior to radiating hot regions that are typically expected to span a smaller angular extent; in the current version, the radiation from elsewhere, if explicitly computed is assumed to be time-invariant supposing the hot regions were not present. To account for radiation from elsewhere, a time-invariant signal is first computed, meaning an axisymmetric surface local-comoving-frame radiation field is assumed. The time-dependent signals from the hot regions are then computed, and modified by subtracting the specific intensity that would otherwise be generated by the surface local-comoving-frame radiation field from elsewhere (i.e., in place of the hot regions).

Star

We can now combine the ambient spacetime, spacetime, and the embedded photosphere, photosphere, into a model star represented by an instance of Star. We do not need to extend this class, so we can simply construct and instantiate a star as follows:

star = xpsi.Star(spacetime = spacetime, photospheres = photosphere)

Let’s check out the star object, which merged parameter subspaces associated with objects lower in the hierarchy:

star

Free parameters
---------------
mass: Gravitational mass [solar masses].
radius: Coordinate equatorial radius [km].
distance: Earth distance [kpc].
cos_inclination: Cosine of Earth inclination to rotation axis.
p__phase_shift: The phase of the hot region, a periodic parameter [cycles].
p__super_colatitude: The colatitude of the centre of the superseding region [radians].
p__super_radius: The angular radius of the (circular) superseding region [radians].
p__super_temperature: log10(superseding region effective temperature [K]).
s__phase_shift: The phase of the hot region, a periodic parameter [cycles].
s__super_colatitude: The colatitude of the centre of the superseding region [radians].
s__super_radius: The angular radius of the (circular) superseding region [radians].

Note that only the free parameters are merged into a subspace higher in the object hierarchy. The reason for this is that there is not a clear and common pattern (at present) for accessing fixed/derived variables outside of the primary subspace to which they belong. If you try hard enough, of course, you can still get at them.

len(star)

11

Let’s specify a parameter vector to the star subspace with the true model parameter values that we injected to generate the synthetic data rendered above:

p = [1.4,
     12.5,
     0.2,
     math.cos(1.25),
     0.0,
     1.0,
     0.075,
     6.2,
     0.025,
     math.pi - 1.0,
     0.2]

star(p)

The order of the parameters is as follows:

star.params

[Gravitational mass [solar masses] = 1.400e+00,
 Coordinate equatorial radius [km] = 1.250e+01,
 Earth distance [kpc] = 2.000e-01,
 Cosine of Earth inclination to rotation axis = 3.153e-01,
 The phase of the hot region, a periodic parameter [cycles] = 0.000e+00,
 The colatitude of the centre of the superseding region [radians] = 1.000e+00,
 The angular radius of the (circular) superseding region [radians] = 7.500e-02,
 log10(superseding region effective temperature [K]) = 6.200e+00,
 The phase of the hot region, a periodic parameter [cycles] = 2.500e-02,
 The colatitude of the centre of the superseding region [radians] = 2.142e+00,
 The angular radius of the (circular) superseding region [radians] = 2.000e-01]

If the secondary temperature was free, we would extend the vector p by one element, passing the injected value of 6.0:

secondary['super_temperature']

6.0

Inspecting functionality

When we simulate a pulsation signal, it is stored as the signal property of the photosphere object. The simulated 1D profile in photosphere are obtained by summing the count-rate over output instrument channels. This is done inside the plot_1d_pulse method.

We need some energies (of photons incident in our rest frame) to integrate at:

energies = np.logspace(-1.0, np.log10(3.0), 128, base=10.0)

We first need to update the star given the parameter vector, and we also show another way to set the parameter value:

star['cos_inclination'] = math.cos(2.0)

star.update()
photosphere.integrate(energies, threads=1) # the number of OpenMP threads to use

The photosphere object does not contain the value of the phase so we use the ones in the HotRegions by setting photosphere_phases=hot.phases_in_cycles.

XpsiPlot.plot_1d_pulse(signal=None,
                       photosphere=photosphere,
                       photosphere_phases=hot.phases_in_cycles)

These pulse profiles are the signals incident on the telescope, before operating on them with the response model. The markers, linearly spaced in phase, denote the phase resolution.

The star.update method which in-turn calls the photosphere.embed method. We then call the photosphere.integrate method, passing the energies. Here we sum over incident specific photon flux pulses as an approximation to integrating over energy.

star['cos_inclination'] = math.cos(1.0)

star.update()
photosphere.integrate(energies, threads=1)

XpsiPlot.plot_1d_pulse(signal=None,
                       photosphere=photosphere,
                       photosphere_phases=hot.phases_in_cycles)

Notice the solid pulses without markers are unchanged from the plot a few cells above, and can be used to guide the eye to the effect of a change in Earth inclination.

Below we print crude representations of the cell meshes spanning each hot region. The elements of a mesh cell-area array which are finite are not all identical: at the boundary of a hot region the proper area elements are smaller because of partial coverage by radiating material. The sum of all finite proper areas effectively equals the total proper area within a hot-region boundary.

XpsiPlot.plot_meshes(regions=hot.objects)

Plotting the meshes for 2 regions, with 1 temperatures

Note that the lowest colatitude row is at zero on the y-axis.

Now, let’s compute the incident specific flux signal normalised to the maximum and plot a pulse in two dimensions.

Note that if want to obtain the specific flux in units of photons/cm\(^{2}\)/s/keV instead, the output of photosphere.signal needs to be divided by distance squared, where distance is measured in meters.

Also note that we can interpolate the signal in phase as follows.

from xpsi.tools import phase_interpolator

# Interpolation in phase with 500 phase bins
phases = np.linspace(0.0, 1, 500)

# Interpolating and summing spot1 (in photosphere.signal[0]) and spot2 (in photosphere.signal[1])
# with their respective signal phase shifts.

# SPOT 1 + SPOT 2
signal_photosphere = phase_interpolator(phases, hot.phases_in_cycles[0], photosphere.signal[0][0], hot['p__phase_shift']) +\
                     phase_interpolator(phases, hot.phases_in_cycles[0], photosphere.signal[1][0], hot['s__phase_shift'])

XpsiPlot.plot_2d_pulse(signal_photosphere,
                       x=phases,
                       y=energies,
                       rotations=5,
                       ylabel=r'Energy (keV)',
                       cbar_label=r'Signal (normalized by maximum)',
                       yticks=[0.2, 0.5, 1.0, 2.0],
                       normalized=True,
                      )

Now we increase the phase resolution, and plot a single rotational pulse:

for obj in hot.objects:
    obj.set_phases(num_leaves = 1024)

star.update()
photosphere.integrate(energies, threads=1)

# Interpolation in phase with 500 phase bins
phases = np.linspace(0.0, 1, 500)

# Interpolating and summing spot1 (in photosphere.signal[0]) and spot2 (in photosphere.signal[1])
# with their respective signal phase shifts.

# SPOT 1 + SPOT 2
signal_photosphere = phase_interpolator(phases, hot.phases_in_cycles[0], photosphere.signal[0][0], hot['p__phase_shift']) +\
                     phase_interpolator(phases, hot.phases_in_cycles[0], photosphere.signal[1][0], hot['s__phase_shift'])

XpsiPlot.plot_2d_pulse(signal_photosphere,
                       x=phases,
                       y=energies,
                       rotations=1,
                       ylabel=r'Energy (keV)',
                       cbar_label=r'Signal (normalized by maximum)',
                       yticks=[0.2, 0.5, 1.0, 2.0],
                       normalized=True,
                      )

Let’s iterate over a monotonically increasing set of values of the hot-region angular radius.

For objects that derive from ParameterSubspace we can get the current parameter vector in several ways:

star()

[1.4,
 12.5,
 0.2,
 0.5403023058681398,
 0.0,
 1.0,
 0.075,
 6.2,
 0.025,
 2.141592653589793,
 0.2]

star() == star.vector # possible shortcut to save some characters

True

Finally, let’s play with some geometric parameters:

fig = plt.figure(figsize=(7,7))
ax = fig.add_subplot(111)
ax.set_ylabel('photons/cm$^2$/s/keV (normalised by maxima)')
ax.set_xlabel('Phase [cycles]')

for obj in hot.objects:
    obj.set_phases(num_leaves = 256)

# let's play with the angular radius of the primary hot region
angular_radii = np.linspace(0.01, 1.0, 10)

star['cos_inclination'] = math.cos(1.25)

for angular_radius in angular_radii:
    star['p__super_radius'] = angular_radius # play time
    star.update()
    photosphere.integrate(energies, threads=3)
    temp = np.sum(photosphere.signal[0][0] + photosphere.signal[1][0], axis=0)
    _ = ax.plot(hot.phases_in_cycles[0], temp/np.max(temp), 'k-', lw=0.5)

star['cos_inclination'] = math.cos(1.0)

for angular_radius in angular_radii:
    star['p__super_radius'] = angular_radius
    star.update()
    photosphere.integrate(energies, threads=3)
    temp = np.sum(photosphere.signal[0][0] + photosphere.signal[1][0], axis=0)
    _ = ax.plot(hot.phases_in_cycles[0], temp/np.max(temp), 'r-', lw=0.5)

star['cos_inclination'] = math.cos(0.5)

for angular_radius in angular_radii:
    star['p__super_radius'] = angular_radius
    star.update()
    photosphere.integrate(energies, threads=3)
    temp = np.sum(photosphere.signal[0][0] + photosphere.signal[1][0], axis=0)
    _ = ax.plot(hot.phases_in_cycles[0], temp/np.max(temp), 'b-', lw=0.5)

XpsiPlot.veneer((0.05,0.2), (0.05,0.2), ax)

Imaging

We can also image the photosphere to improve visualisation of the signal generation process. This is very brief: a more details are provided in the Global surface emission tutorial.

spacetime.a = 0.0 # spacetime spin parameter (~angular momentum)
spacetime.q = 0.0 # spacetime mass quadrupole moment

sky_map_kwargs = {'panel_indices': (0,1,2,3,4,5),
                  'num_levels': 100, # in intensity field rendering
                  'colormap': cm.Purples_r,   # a colormap like this will make the lowest finite
                  'phase_average': False,     # intensity distinct from the black zero-intensity
                  'annotate_energies': True,  # background from the surface and behind the star
                  'energy_annotation_format': '[%.2f keV]',
                  'annotate_location': (0.025,0.025)}

animate_kwargs = {'cycles': 8, 'fps': 32}

# later we will choose a higher phase resolution in
# order to phase-average, because it reduces artefacts
# that manifest due to the hard boundary of the hot regions
for obj in hot.objects:
    obj.set_phases(num_leaves = 96) # animation will have 3 second period

star['p__super_radius'] = 0.1

photosphere.image(reimage = True,
                  cache_intensities = 1.0, # cache size limit in GBs
                  energies = np.array([0.01,0.1,0.5,1.0,2.0,5.0]),
                  phases = hot.phases_in_cycles[0] * _2pi,
                  sqrt_num_rays = 400,
                  threads = 4,
                  max_steps = 100000,
                  epsrel_ray = 1.0e-12,
                  plot_sky_maps = True, # activate if you want to plot frames
                  sky_map_kwargs = sky_map_kwargs,
                  animate_sky_maps = True, # activate if you want to animate
                  free_memory = True, # try to free ray-map/intensity cache memory before animating
                  animate_kwargs = animate_kwargs)

Imaging the star...
Warning: a ray map has not been cached... tracing new ray set...
Commencing ray tracing and imaging...
Imaging the star and computing light-curves...
Calculating image plane boundary...
Image plane semi-major: 1.00
Image plane semi-minor: 1.00
Thread 0 is tracing annulus #0 of rays.
Thread 0 is tracing annulus #100 of rays.
Thread 0 is tracing annulus #200 of rays.
Thread 0 is tracing annulus #300 of rays.
Warning: cosine aberrated surface zenith angle = -1.55991567e-03...
Warning: forcing the aberrated ray to be tangential.
Warning: cosine aberrated surface zenith angle = -9.39012537e-04...
Warning: forcing the aberrated ray to be tangential.
Warning: cosine aberrated surface zenith angle = -8.31628268e-05...
Warning: forcing the aberrated ray to be tangential.
Warning: cosine aberrated surface zenith angle = -6.39549950e-04...
Warning: forcing the aberrated ray to be tangential.
Warning: cosine aberrated surface zenith angle = -5.11681429e-04...
Warning: forcing the aberrated ray to be tangential.
Warning: cosine aberrated surface zenith angle = -3.07033381e-03...
Warning: forcing the aberrated ray to be tangential.
Warning: cosine aberrated surface zenith angle = -3.07255372e-03...
Warning: forcing the aberrated ray to be tangential.
Warning: cosine aberrated surface zenith angle = -5.17969130e-04...
Warning: forcing the aberrated ray to be tangential.
Warning: cosine aberrated surface zenith angle = -6.48373520e-04...
Warning: forcing the aberrated ray to be tangential.
Warning: cosine aberrated surface zenith angle = -8.50120359e-05...
Warning: forcing the aberrated ray to be tangential.
Warning: cosine aberrated surface zenith angle = -1.02009401e-03...
Warning: forcing the aberrated ray to be tangential.
Warning: cosine aberrated surface zenith angle = -1.69951233e-03...
Warning: forcing the aberrated ray to be tangential.

Global ray map computed.
Coordinates transformed from Boyer-Lindquist to spherical.
Commencing imaging...Ray tracing complete.
Ray set cached.
Intensity caching complete.
Plotting intensity sky maps...
Normalising each sky map panel separately...
Normalised sky map panels separately.
Rendering image numbers [1, 10]...
Rendering image numbers (10, 20]...
Rendering image numbers (20, 30]...
Rendering image numbers (30, 40]...
Rendering image numbers (40, 50]...
Rendering image numbers (50, 60]...
Rendering image numbers (60, 70]...
Rendering image numbers (70, 80]...
Rendering image numbers (80, 90]...
Rendering image numbers (90, 100]...
Intensity sky maps plotted.
Animating intensity sky maps...
Writing to disk: ./images/skymap_animated.mp4...
Intensity sky maps animated.
Star imaged.

If getting ValueError: unknown file extension: .mp4 from the code block above, try installing ffmpeg e.g., with conda install conda-forge/label/gcc7::ffmpeg" and then re-start the kernel and re-run the notebook.

%%HTML
<div align="middle">
<video width="100%" controls loop>
    <source src="images/skymap_animated.mp4" type="video/mp4">
</video></div>

!mkdir images/frames_twospots
!mv images/*.png images/frames_twospots/.

mkdir: cannot create directory ‘images/frames_twospots’: File exists

from IPython.display import Image

Image("./images/frames_twospots/skymap_0.png")

We freed memory above before animating, and thus lost the cached ray map and intensity cache. The ray map will now be recomputed, and the intensities cached again:

sky_map_kwargs['phase_average'] = True
# only one set of panels so why not choose higher res.?
sky_map_kwargs['num_levels'] = 500

for obj in hot.objects:
    obj.set_phases(num_leaves = 260)

photosphere.image(reimage = True,
                  cache_intensities = 1.0, # cache size limit in GBs
                  energies = np.array([0.01,0.1,0.5,1.0,2.0,5.0]),
                  phases = hot.phases_in_cycles[0] * _2pi,
                  sqrt_num_rays = 400,
                  threads = 4,
                  max_steps = 100000,
                  epsrel_ray = 1.0e-12,
                  plot_sky_maps = True,
                  sky_map_kwargs = sky_map_kwargs)

Imaging the star...
Warning: a ray map has not been cached... tracing new ray set...
Commencing ray tracing and imaging...
Imaging the star and computing light-curves...
Calculating image plane boundary...
Image plane semi-major: 1.00
Image plane semi-minor: 1.00
Thread 0 is tracing annulus #0 of rays.
Thread 0 is tracing annulus #100 of rays.
Thread 0 is tracing annulus #200 of rays.
Thread 0 is tracing annulus #300 of rays.
Warning: cosine aberrated surface zenith angle = -1.55991567e-03...
Warning: forcing the aberrated ray to be tangential.
Warning: cosine aberrated surface zenith angle = -9.39012537e-04...
Warning: forcing the aberrated ray to be tangential.
Warning: cosine aberrated surface zenith angle = -8.31628268e-05...
Warning: forcing the aberrated ray to be tangential.
Warning: cosine aberrated surface zenith angle = -6.39549950e-04...
Warning: forcing the aberrated ray to be tangential.
Warning: cosine aberrated surface zenith angle = -5.11681429e-04...
Warning: forcing the aberrated ray to be tangential.
Warning: cosine aberrated surface zenith angle = -3.07033381e-03...
Warning: forcing the aberrated ray to be tangential.
Warning: cosine aberrated surface zenith angle = -3.07255372e-03...
Warning: forcing the aberrated ray to be tangential.
Warning: cosine aberrated surface zenith angle = -5.17969130e-04...
Warning: forcing the aberrated ray to be tangential.
Warning: cosine aberrated surface zenith angle = -6.48373520e-04...
Warning: forcing the aberrated ray to be tangential.
Warning: cosine aberrated surface zenith angle = -8.50120359e-05...
Warning: forcing the aberrated ray to be tangential.
Warning: cosine aberrated surface zenith angle = -1.02009401e-03...
Warning: forcing the aberrated ray to be tangential.
Warning: cosine aberrated surface zenith angle = -1.69951233e-03...
Warning: forcing the aberrated ray to be tangential.

Global ray map computed.
Coordinates transformed from Boyer-Lindquist to spherical.
Commencing imaging...Ray tracing complete.
Ray set cached.
Intensity caching complete.
Plotting intensity sky maps...
Averaging (specific) intensity over rotational phase...
Averaged (specific) intensity over rotational phase.
Normalising each sky map panel separately...
Normalised sky map panels separately.
Rendering phase-averaged images...
Intensity sky maps plotted.
Star imaged.

!mv images/skymap_0.png images/frames_twospots/skymap_phase_averaged.png

Image("images/frames_twospots/skymap_phase_averaged.png")

Artefacts are not very discernable, but can be further reduced with a higher phase resolution at higher memory and time cost.

Summary

In this notebook we constructed a model star and simulated some pulsed X-ray signals.