Surface radiation field
Extensions to compute specific intensities in a local comoving frame.
Overview
We require the surface radiation field model to be implemented in C extensions to Python if signal computation is to be sufficiently fast. Signal computation includes:
integrating over the image of a star on a distant observer’s sky, yielding a phase-energy resolved signal for likelihood function evaluation;
resolving the image of star on a distant observer’s sky as a function of phase and energy for visualisation purposes, and for more general simulation purposes.
The first type of signal computation is enabled by surface discretization,
whilst the latter is enabled via image-plane discretization. For surface
discretization there are several 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 (they can be selected using appropriate flags when creating
instances of the radiating regions); 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 (by default they are the same as hot_BB.pyx
).
These extension modules contain functions whose bodies enable evaluation of
specific intensities with respect to local comoving frames at the stellar surface.
The hot_*.pyx
modules define the radiation field inside surface hot regions
represented by HotRegion
instances; however, the same
modules are used by an instance of the Everywhere
class.
The elsewhere_user.pyx
extension module defines the customized radiation
field exterior to the HotRegion
instances on the surface.
In case no customized module for the exterior surface is needed, the surface
radiation field for exterior surface can be selected from one of the built-in
options for the hot regions.
For image-plane discretization, the hot_*.pyx
extensions define the local
surface radiation field properties with respect to a comoving frame. However,
to compute the variables required by this module to evaluate specific intensity,
another extension module is required with transforms a set of global
variables into local variables given spacetime coordinates of an event
belonging to a null geodesic connecting the surface to a radiation-detecting
distant instrument. Users and developers must customize the function bodies
in local_variables.pyx
in order to implement more advanced physics models.
For all of these modules, numerical data can be read from disk to define the
functions for specific intensity evaluation given local variables, and for
transformation of global variables to local variables. A common usage pattern
is to statistically model data assuming some numerical physics model for the
radiation emergent from an atmosphere. Statistical modeling requires many
likelihood function evaluations, and therefore to avoid unnecessarily reading
numerical data from disk every time the likelihood function is called, users
can preload the numerical atmosphere data in an encapsulation Python process
and point to the data structures in memory when computing signals; developers
can develop extension modules to work with preloaded numerical atmosphere data.
When preloading an atmosphere, the number of grid points in each dimension
is dynamically inferred, but the existing source code for multi-dimensional
interpolation (hot_Num4D.pyx
) is hard-coded for four dimensions. With
development, this can be made more sophisticated.
Developers can contribute extension modules to the
surface_radiation_field/archive
. Extensions from this archive must replace
the contents of the hot_user.pyx
, elsewhere_user.pyx
, and local_variables.pyx
modules, adhering to the function prototypes defined in the associated
header files and using existing examples in the archive for guidance.
Use the Python-exposed functions below to test the C implementations of surface radiation field models.
Wrappers
- xpsi.surface_radiation_field.intensity(double[::1] energies, double[::1] mu, double[:, ::1] local_variables, atmosphere=None, int stokesQ=0, region_extension=u'hot', atmos_extension=u'BB', beam_opt=0, size_t numTHREADS=1)
Evaluate the photon specific intensity using an extension module.
The intensities are computed directly from local variable values.
- Parameters:
energies (double[::1]) – A 1D
numpy.ndarray
of energies in keV to evaluate photon specific intensity at, in the local comoving frame.mu (double[::1]) – A 1D
numpy.ndarray
of angles at which to evaluate the photon specific intensity at. Specifically, the angle is the cosine of the angle of the ray direction to the local surface normal, in the local comoving frame.local_variables (double[:,::1]) – A 2D
numpy.ndarray
(with beam_opt=0) or 7Dnumpy.ndarray
(with beam_opt > 0) of local variables such as temperature, effective gravity, and beaming parameters. Rows correspond to the sequence of points in the space of energy and angle specified above, and columns contain the required variable values, one variable per column, in the expected order.atmosphere (tuple) –
A numerical atmosphere preloaded into buffers. The buffers must be supplied in an \(n\)-tuple whose \(n^{th}\) element is an \((n-1)\)-dimensional array flattened into a one-dimensional
numpy.ndarray
. The first \(n-1\) elements of the \(n\)-tuple must each be an ordered one-dimensionalnumpy.ndarray
of parameter values for the purpose of multi-dimensional interpolation in the \(n^{th}\) buffer. The first \(n-1\) elements must be ordered to match the index arithmetic applied to the \(n^{th}\) buffer. An example would be(logT, logg, mu, logE, buf)
, where:logT
is a logarithm of local comoving effective temperature;logg
is a logarithm of effective surface gravity;mu
is the cosine of the angle from the local surface normal;logE
is a logarithm of the photon energy; andbuf
is a one-dimensional buffer of intensities of size given by the product of sizes of the first \(n-1\) tuple elements.It is highly recommended that buffer preloading is used, instead of loading from disk in the customisable radiation field extension module, to avoid reading from disk for every signal (likelihood) evaluation. This can be a non-negligible waste of compute resources. By preloading in Python, the memory is allocated and references to that memory are not in general deleted until a sampling script exits and the kernel stops. The likelihood callback accesses the same memory upon each call without I/O.
StokesQ (int) – If StokesQ=1, Stokes Q intensities will be returned. Otherwise Stokes I intensities. Default is StokesQ=0.
region_extension (str) – Specify the radiating region extension module to invoke. Options are
'hot'
and'elsewhere'
.atmos_extension (str) – Used to determine which atmospheric extension to use. Options at the moment: “BB”: Analytical blackbody, “Num4D”: Numerical atmosphere using 4D-interpolation from the provided atmosphere data, “Pol_BB_burst”: Polarized analytical blackbody+burst approximation, “Pol_Num2D”: Polarized numerical atmosphere using 2D-interpolation from the provided atmosphere data, “user”: A user-provided extension which can be set up by replacing the contents of the file hot_user.pyx (and elsewhere_user.pyx if needed) and re-installing X-PSI (if not changed, “user” is the same as “BB”).
beam_opt (int) – Used to determine which atmospheric beaming modification model to use. Options at the moment: 0: No modification (default) 1: Original*beaming_correction without re-normalization 2: Original*beaming_correction with analytical re-normalization estimate 3: Original*beaming_correction with numerical re-normalization
numTHREADS (int) – Number of OpenMP threads to launch.
- Returns:
A 1D
numpy.ndarray
of the photon specific intensities in units of photons/s/keV/cm^2/sr.
- xpsi.surface_radiation_field.intensity_from_globals(double[::1] energies, double[::1] mu, double[::1] colatitude, double[::1] azimuth, double[::1] phase, double[::1] global_variables, double R_eq, double zeta, double epsilon, atmosphere=None, atmos_extension=u'BB', size_t numTHREADS=1)
Evaluate the photon specific intensity using extension modules.
The local variables are computed from global variables and spacetime coordinates of events at the surface.
- Parameters:
energies (double[::1]) – A 1D
numpy.ndarray
of energies in keV to evaluate photon specific intensity at, in the local comoving frame.mu (double[::1]) – A 1D
numpy.ndarray
of angles at which to evaluate the photon specific intensity at. Specifically, the angle is the cosine of the angle of the ray direction to the local surface normal, in the local comoving frame.colatitude (double[::1]) – A 1D
numpy.ndarray
of colatitudes in radians at which to evaluate photon specific intensity, in the local comoving frame. This is the coordinate of a fixed point relative to a distant static observer.azimuth (double[::1]) – A 1D
numpy.ndarray
of azimuths in radians at which to evaluate photon specific intensity, in the local comoving frame. This is the coordinate of a fixed point relative to a distant static observer.phase (double[::1]) – A 1D
numpy.ndarray
of rotational phases in radians at which to evaluate photon specific intensity, in the local comoving frame. The phase describes the rotation of the surface radiation field through the fixed points alluded to above. A localised hot region for instance rotates so that it’s azimuth changes relative to a point with fixed azimuth described above. Specifying zero phase means that the intensity is evaluated at the fixed points for the initial configuration of the hot region (or more general radiation field).global_variables (double[::1]) – A 1D
numpy.ndarray
of global variables controlling the surface radiation field.R_eq (double) – The equatorial radius in metres. Needed for effective gravity universal relation.
zeta (double) – See the
Spacetime
class property. Needed for effective gravity universal relation. It is recommended to supply this dimensionless variable by using an instance ofSpacetime
.epsilon (double) – See the
Spacetime
class property. Needed for effective gravity universal relation. It is recommended to supply this dimensionless variable by using an instance ofSpacetime
.atmosphere (tuple) –
A numerical atmosphere preloaded into buffers. The buffers must be supplied in an \(n\)-tuple whose \(n^{th}\) element is an \((n-1)\)-dimensional array flattened into a one-dimensional
numpy.ndarray
. The first \(n-1\) elements of the \(n\)-tuple must each be an ordered one-dimensionalnumpy.ndarray
of parameter values for the purpose of multi-dimensional interpolation in the \(n^{th}\) buffer. The first \(n-1\) elements must be ordered to match the index arithmetic applied to the \(n^{th}\) buffer. An example would be(logT, logg, mu, logE, buf)
, where:logT
is a logarithm of local comoving effective temperature;logg
is a logarithm of effective surface gravity;mu
is the cosine of the angle from the local surface normal;logE
is a logarithm of the photon energy; andbuf
is a one-dimensional buffer of intensities of size given by the product of sizes of the first \(n-1\) tuple elements.It is highly recommended that buffer preloading is used, instead of loading from disk in the customisable radiation field extension module, to avoid reading from disk for every signal (likelihood) evaluation. This can be a non-negligible waste of compute resources. By preloading in Python, the memory is allocated and references to that memory are not in general deleted until a sampling script exits and the kernel stops. The likelihood callback accesses the same memory upon each call without I/O.
atmos_extension (str) – Used to determine which atmospheric extension to use. Options at the moment: “BB”: Analytical blackbody, “Num4D”: Numerical atmosphere using 4D-interpolation from the provided atmosphere data, “user”: A user-provided extension which can be set up by replacing the contents of the file hot_user.pyx (and elsewhere_user.pyx if needed) and re-installing X-PSI (if not changed, “user” is the same as “BB”).
numTHREADS (int) – Number of OpenMP threads to launch.
- Returns:
A 1D
numpy.ndarray
of the photon specific intensities in units of photons/s/keV/cm^2/sr.
- xpsi.surface_radiation_field.effective_gravity(double[::1] cos_colatitude, double[::1] R_eq, double[::1] zeta, double[::1] epsilon)
Approximate local effective gravity using a universal relation.
See Morsink et al. (2007), AlGendy & Morsink (2014), and Bogdanov et al. (2019, ApJL, 887, L26).
- Parameters:
cos_colatitude (double[::1]) – A 1D
numpy.ndarray
of colatitudes, with respect to a rotational pole, at which to evaluate the local effective gravity. on the surface. Specifically, the angle is the cosine of the colatitude of the ray direction to the local surface normal, in the local comoving frame.R_eq (double[::1]) – A 1D
numpy.ndarray
of surface coordinate equatorial radii in metres.zeta (double[::1]) – A 1D
numpy.ndarray
of a dimensionless function of stellar properties. Seezeta
.epsilon (double[::1]) – A 1D
numpy.ndarray
of a dimensionless function of stellar properties. Seeepsilon
.
- Returns:
A 1D
numpy.ndarray
of the logarithms (base 10) of the effective gravities in units of cm/s^2.
A numerical atmosphere
The following is the extension module hot_Num4D.pyx
for surface radiation field evaluation
that may be found on path surface_radiation_field/hot_Num4D.pyx
.
This extension works with numerical atmosphere preloading to execute cubic polynomial interpolation
in four dimensions.
#cython: cdivision=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
from libc.math cimport M_PI, sqrt, sin, cos, acos, log10, pow, exp, fabs
from libc.stdio cimport printf, fopen, fclose, fread, FILE
from GSL cimport gsl_isnan, gsl_isinf
from libc.stdlib cimport malloc, free
from xpsi.global_imports import _keV, _k_B, _h_keV
cdef int SUCCESS = 0
cdef int ERROR = 1
cdef double erg = 1.0e-7
cdef double k_B = _k_B
cdef double keV = _keV
cdef double h_keV = _h_keV
cdef double k_B_over_keV = k_B / keV
cdef int VERBOSE = 0
ctypedef struct ACCELERATE:
size_t **BN # base node for interpolation
double **node_vals
double **SPACE
double **DIFF
double **INTENSITY_CACHE
double **VEC_CACHE
# Modify this struct if useful for the user-defined source radiation field.
# Note that the members of DATA will be shared by all threads and are
# statically allocated, whereas the members of ACCELERATE will point to
# dynamically allocated memory, not shared by threads.
ctypedef struct DATA:
const _preloaded *p
ACCELERATE acc
#----------------------------------------------------------------------->>>
# >>> 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_Num4D(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.
# Return NULL if dynamic memory is not required for the model
if preloaded == NULL :
printf("ERROR: The numerical atmosphere data were not preloaded, which are required by this extension.\n")
cdef DATA *D = <DATA*> malloc(sizeof(DATA))
D.p = preloaded
D.p.BLOCKS[0] = 64
D.p.BLOCKS[1] = 16
D.p.BLOCKS[2] = 4
cdef size_t T, i, j, k, l
D.acc.BN = <size_t**> malloc(numThreads * sizeof(size_t*))
D.acc.node_vals = <double**> malloc(numThreads * sizeof(double*))
D.acc.SPACE = <double**> malloc(numThreads * sizeof(double*))
D.acc.DIFF = <double**> malloc(numThreads * sizeof(double*))
D.acc.INTENSITY_CACHE = <double**> malloc(numThreads * sizeof(double*))
D.acc.VEC_CACHE = <double**> malloc(numThreads * sizeof(double*))
for T in range(numThreads):
D.acc.BN[T] = <size_t*> malloc(D.p.ndims * sizeof(size_t))
D.acc.node_vals[T] = <double*> malloc(2 * D.p.ndims * sizeof(double))
D.acc.SPACE[T] = <double*> malloc(4 * D.p.ndims * sizeof(double))
D.acc.DIFF[T] = <double*> malloc(4 * D.p.ndims * sizeof(double))
D.acc.INTENSITY_CACHE[T] = <double*> malloc(256 * sizeof(double))
D.acc.VEC_CACHE[T] = <double*> malloc(D.p.ndims * sizeof(double))
for i in range(D.p.ndims):
D.acc.BN[T][i] = 0
D.acc.VEC_CACHE[T][i] = D.p.params[i][1]
D.acc.node_vals[T][2*i] = D.p.params[i][1]
D.acc.node_vals[T][2*i + 1] = D.p.params[i][2]
j = 4*i
D.acc.SPACE[T][j] = 1.0 / (D.p.params[i][0] - D.p.params[i][1])
D.acc.SPACE[T][j] /= D.p.params[i][0] - D.p.params[i][2]
D.acc.SPACE[T][j] /= D.p.params[i][0] - D.p.params[i][3]
D.acc.SPACE[T][j + 1] = 1.0 / (D.p.params[i][1] - D.p.params[i][0])
D.acc.SPACE[T][j + 1] /= D.p.params[i][1] - D.p.params[i][2]
D.acc.SPACE[T][j + 1] /= D.p.params[i][1] - D.p.params[i][3]
D.acc.SPACE[T][j + 2] = 1.0 / (D.p.params[i][2] - D.p.params[i][0])
D.acc.SPACE[T][j + 2] /= D.p.params[i][2] - D.p.params[i][1]
D.acc.SPACE[T][j + 2] /= D.p.params[i][2] - D.p.params[i][3]
D.acc.SPACE[T][j + 3] = 1.0 / (D.p.params[i][3] - D.p.params[i][0])
D.acc.SPACE[T][j + 3] /= D.p.params[i][3] - D.p.params[i][1]
D.acc.SPACE[T][j + 3] /= D.p.params[i][3] - D.p.params[i][2]
D.acc.DIFF[T][j] = D.acc.VEC_CACHE[T][i] - D.p.params[i][1]
D.acc.DIFF[T][j] *= D.acc.VEC_CACHE[T][i] - D.p.params[i][2]
D.acc.DIFF[T][j] *= D.acc.VEC_CACHE[T][i] - D.p.params[i][3]
D.acc.DIFF[T][j + 1] = D.acc.VEC_CACHE[T][i] - D.p.params[i][0]
D.acc.DIFF[T][j + 1] *= D.acc.VEC_CACHE[T][i] - D.p.params[i][2]
D.acc.DIFF[T][j + 1] *= D.acc.VEC_CACHE[T][i] - D.p.params[i][3]
D.acc.DIFF[T][j + 2] = D.acc.VEC_CACHE[T][i] - D.p.params[i][0]
D.acc.DIFF[T][j + 2] *= D.acc.VEC_CACHE[T][i] - D.p.params[i][1]
D.acc.DIFF[T][j + 2] *= D.acc.VEC_CACHE[T][i] - D.p.params[i][3]
D.acc.DIFF[T][j + 3] = D.acc.VEC_CACHE[T][i] - D.p.params[i][0]
D.acc.DIFF[T][j + 3] *= D.acc.VEC_CACHE[T][i] - D.p.params[i][1]
D.acc.DIFF[T][j + 3] *= D.acc.VEC_CACHE[T][i] - D.p.params[i][2]
cdef double *address = NULL
# Cache intensity
for T in range(numThreads):
for i in range(4):
for j in range(4):
for k in range(4):
for l in range(4):
address = D.p.I + (D.acc.BN[T][0] + i) * D.p.S[0]
address += (D.acc.BN[T][1] + j) * D.p.S[1]
address += (D.acc.BN[T][2] + k) * D.p.S[2]
address += D.acc.BN[T][3] + l
D.acc.INTENSITY_CACHE[T][i * D.p.BLOCKS[0] + j * D.p.BLOCKS[1] + k * D.p.BLOCKS[2] + l] = address[0]
# Cast for generalised usage in integration routines
return <void*> D
cdef int free_hot_Num4D(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_hot()
# because data is expected to be NULL in this case
cdef DATA *D = <DATA*> data
cdef size_t T
for T in range(numThreads):
free(D.acc.BN[T])
free(D.acc.node_vals[T])
free(D.acc.SPACE[T])
free(D.acc.DIFF[T])
free(D.acc.INTENSITY_CACHE[T])
free(D.acc.VEC_CACHE[T])
free(D.acc.BN)
free(D.acc.node_vals)
free(D.acc.SPACE)
free(D.acc.DIFF)
free(D.acc.INTENSITY_CACHE)
free(D.acc.VEC_CACHE)
free(D)
return SUCCESS
#----------------------------------------------------------------------->>>
# >>> Cubic polynomial interpolation.
# >>> Improve acceleration properties... i.e. do not recompute numerical
# ... weights or re-read intensities if not necessary.
#----------------------------------------------------------------------->>>
cdef double eval_hot_Num4D(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
# This function must cast the void pointer appropriately for use.
cdef DATA *D = <DATA*> data
cdef:
size_t i = 0, ii
double I = 0.0, temp
double *node_vals = D.acc.node_vals[THREAD]
size_t *BN = D.acc.BN[THREAD]
double *SPACE = D.acc.SPACE[THREAD]
double *DIFF = D.acc.DIFF[THREAD]
double *I_CACHE = D.acc.INTENSITY_CACHE[THREAD]
double *V_CACHE = D.acc.VEC_CACHE[THREAD]
double vec[4]
double E_eff = k_B_over_keV * pow(10.0, VEC[0])
int update_baseNode[4]
int CACHE = 0
vec[0] = VEC[0]
vec[1] = VEC[1]
vec[2] = mu
vec[3] = log10(E / E_eff)
while i < D.p.ndims:
# if parallel == 31:
# printf("\nDimension: %d", <int>i)
update_baseNode[i] = 0
if vec[i] < node_vals[2*i] and BN[i] != 0:
# if parallel == 31:
# printf("\nExecute block 1: %d", <int>i)
update_baseNode[i] = 1
while vec[i] < D.p.params[i][BN[i] + 1]:
# if parallel == 31:
# printf("\n!")
# printf("\nvec i: %.8e", vec[i])
# printf("\nBase node: %d", <int>BN[i])
if BN[i] > 0:
BN[i] -= 1
elif vec[i] <= D.p.params[i][0]:
vec[i] = D.p.params[i][0]
break
elif BN[i] == 0:
break
node_vals[2*i] = D.p.params[i][BN[i] + 1]
node_vals[2*i + 1] = D.p.params[i][BN[i] + 2]
# if parallel == 31:
# printf("\nEnd Block 1: %d", <int>i)
elif vec[i] > node_vals[2*i + 1] and BN[i] != D.p.N[i] - 4:
# if parallel == 31:
# printf("\nExecute block 2: %d", <int>i)
update_baseNode[i] = 1
while vec[i] > D.p.params[i][BN[i] + 2]:
if BN[i] < D.p.N[i] - 4:
BN[i] += 1
elif vec[i] >= D.p.params[i][D.p.N[i] - 1]:
vec[i] = D.p.params[i][D.p.N[i] - 1]
break
elif BN[i] == D.p.N[i] - 4:
break
node_vals[2*i] = D.p.params[i][BN[i] + 1]
node_vals[2*i + 1] = D.p.params[i][BN[i] + 2]
# if parallel == 31:
# printf("\nEnd Block 2: %d", <int>i)
# if parallel == 31:
# printf("\nTry block 3: %d", <int>i)
if V_CACHE[i] != vec[i] or update_baseNode[i] == 1:
# if parallel == 31:
# printf("\nExecute block 3: %d", <int>i)
ii = 4*i
DIFF[ii] = vec[i] - D.p.params[i][BN[i] + 1]
DIFF[ii] *= vec[i] - D.p.params[i][BN[i] + 2]
DIFF[ii] *= vec[i] - D.p.params[i][BN[i] + 3]
DIFF[ii + 1] = vec[i] - D.p.params[i][BN[i]]
DIFF[ii + 1] *= vec[i] - D.p.params[i][BN[i] + 2]
DIFF[ii + 1] *= vec[i] - D.p.params[i][BN[i] + 3]
DIFF[ii + 2] = vec[i] - D.p.params[i][BN[i]]
DIFF[ii + 2] *= vec[i] - D.p.params[i][BN[i] + 1]
DIFF[ii + 2] *= vec[i] - D.p.params[i][BN[i] + 3]
DIFF[ii + 3] = vec[i] - D.p.params[i][BN[i]]
DIFF[ii + 3] *= vec[i] - D.p.params[i][BN[i] + 1]
DIFF[ii + 3] *= vec[i] - D.p.params[i][BN[i] + 2]
V_CACHE[i] = vec[i]
# if parallel == 31:
# printf("\nEnd block 3: %d", <int>i)
# if parallel == 31:
# printf("\nTry block 4: %d", <int>i)
if update_baseNode[i] == 1:
# if parallel == 31:
# printf("\nExecute block 4: %d", <int>i)
CACHE = 1
SPACE[ii] = 1.0 / (D.p.params[i][BN[i]] - D.p.params[i][BN[i] + 1])
SPACE[ii] /= D.p.params[i][BN[i]] - D.p.params[i][BN[i] + 2]
SPACE[ii] /= D.p.params[i][BN[i]] - D.p.params[i][BN[i] + 3]
SPACE[ii + 1] = 1.0 / (D.p.params[i][BN[i] + 1] - D.p.params[i][BN[i]])
SPACE[ii + 1] /= D.p.params[i][BN[i] + 1] - D.p.params[i][BN[i] + 2]
SPACE[ii + 1] /= D.p.params[i][BN[i] + 1] - D.p.params[i][BN[i] + 3]
SPACE[ii + 2] = 1.0 / (D.p.params[i][BN[i] + 2] - D.p.params[i][BN[i]])
SPACE[ii + 2] /= D.p.params[i][BN[i] + 2] - D.p.params[i][BN[i] + 1]
SPACE[ii + 2] /= D.p.params[i][BN[i] + 2] - D.p.params[i][BN[i] + 3]
SPACE[ii + 3] = 1.0 / (D.p.params[i][BN[i] + 3] - D.p.params[i][BN[i]])
SPACE[ii + 3] /= D.p.params[i][BN[i] + 3] - D.p.params[i][BN[i] + 1]
SPACE[ii + 3] /= D.p.params[i][BN[i] + 3] - D.p.params[i][BN[i] + 2]
# if parallel == 31:
# printf("\nEnd block 4: %d", <int>i)
i += 1
cdef size_t j, k, l, INDEX, II, JJ, KK
cdef double *address = NULL
# Combinatorics over nodes of hypercube; weight cgs intensities
for i in range(4):
II = i * D.p.BLOCKS[0]
for j in range(4):
JJ = j * D.p.BLOCKS[1]
for k in range(4):
KK = k * D.p.BLOCKS[2]
for l in range(4):
address = D.p.I + (BN[0] + i) * D.p.S[0]
address += (BN[1] + j) * D.p.S[1]
address += (BN[2] + k) * D.p.S[2]
address += BN[3] + l
temp = DIFF[i] * DIFF[4 + j] * DIFF[8 + k] * DIFF[12 + l]
temp *= SPACE[i] * SPACE[4 + j] * SPACE[8 + k] * SPACE[12 + l]
INDEX = II + JJ + KK + l
if CACHE == 1:
I_CACHE[INDEX] = address[0]
I += temp * I_CACHE[INDEX]
#if gsl_isnan(I) == 1:
#printf("\nIntensity: NaN; Index [%d,%d,%d,%d] ",
#<int>BN[0], <int>BN[1], <int>BN[2], <int>BN[3])
#printf("\nBase-nodes [%d,%d,%d,%d] ",
#<int>BN[0], <int>BN[1], <int>BN[2], <int>BN[3])
if I < 0.0:
return 0.0
return I * pow(10.0, 3.0 * vec[0])
cdef double eval_hot_norm_Num4D() 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 / h_keV