6.3. Integrations (integr_los.h)

In the end, all calculations amount to multi-dimensional integrations repeated a great many number of times. Because the quantities to integrate are mostly power-laws of monotonic functions, we can optimise the integration for speed. Particularly sophisticated optimisations have been developed for the line-of-sight integration over extended (“resolved”) haloes. This is the core of the integr_los.h library and explained in this section.

6.3.1. Formal equations

For many quantities calculated along a given line of sight (l.o.s.) over a solid angle \(\Delta\Omega\) (J-factor for the smooth and clumpy component, number of clumps in a FOV, etc.), we have to perform the following integration \(I_{l.o.s.}\) on the generic function \(f\) (valid for spherical or triaxial haloes):

\[I_{l.o.s.}(\psi,\theta, \Delta\Omega) = \int_{\Delta\Omega} {\rm d}\Omega\int_{l.o.s.} f(l,\psi,\theta) \,{\rm d}l\,,\]

where the l.o.s. is defined by the angles \((\psi,\theta)\) w.r.t. the Galactic centre direction, as shown in Fig. 6.2. The integration \(\int_{\Delta\Omega}\) over the solid angle is described by the two angles \((\alpha,\beta)\). In practice, \((l,\alpha,\beta)\) are spherical coordinates over which the integration is performed, and for which the origin of the framework is Earth location (also see Fig. 6.2). Hence, using

\[\int_{\Delta\Omega} {\rm d}\Omega= \int_{0}^{2\pi} d\beta \int_{0}^{\alpha_{\rm int}} \sin(\alpha) {\rm d}\alpha\, ,\]

where \(\Delta\Omega\) and \(\alpha_{\rm int}\) are related by

\[\Delta\Omega = 2\pi\times[1-\cos(\alpha_{\rm int})] \, ,\]

we can rewrite \(I_{l.o.s.}\) to be

(6.2)\[I_{l.o.s.}(\psi,\theta, \Delta\Omega) = \int_{0}^{2\pi} F_\beta \,{\rm d}\beta\,,\]


(6.3)\[F_\beta = \int_{0}^{\alpha_{\rm int}} F_{(\beta,\alpha)} \,{\rm d}\alpha\]


(6.4)\[F_{(\beta,\alpha)} = \sin(\alpha) \int_{l_{\rm min}}^{l_{\rm max}} f(l,\beta,\alpha; \psi,\theta) \,{\rm d}l\,.\]

In practice, this amounts to performing three integrations: over \(\beta\) in (6.2), over \(\alpha\) in (6.3), and over \(l\) in (6.4). CLUMPY relies on a Simpson integration scheme in linear or logarithmic adaptive steps, and to optimise the computation time, the integration along \(l\) is broken down into several parts, as described in the following:.

6.3.2. Strategy and tricks

The DM halo density profile can be very cuspy at its centre, and thus it is hard to integrate over the l.o.s. with a regular integration step. We can take advantage of the quasi-sphericity of the clumps to tackle this problem. In any case, integrating a clump or the smooth DM halo requires, most of the time, a similar strategy.


For the smooth component, when we integrate over the solid angle, the distance \(l_{\rm max}\) (distance from Earth to the Galaxy DM halo boundary) slightly changes. However, for instance, for \(\Delta\Omega<10^{-3}\), the change is smaller than a few \(\rm pc\): as in practice there is no clear boundary of the DM halo, this is unimportant, and for simplicity, we keep the same \(l_{\rm max}\) for the integration over \((\alpha,\beta)\).

General integration over \(\alpha\) and \(\beta\)

The integration range is usually small on these parameters and a standard adaptive Simpson routine using linear steps is used. Indeed, the integral is not varying much over this range, unless, e.g., we wish to integrate slightly offset from the density peak centre, in which case the linear integration scheme used on these variable is less efficient. But still, the calculation converges, and we did not make any attempt to improve the situation in this very small l.o.s. region.

Strategy for the l.o.s. integration

DM profiles are decreasing power-law functions of their radius \(r\) (see Halo profiles (profiles.h)). To discuss our strategy, we shall consider the following symbolic integration:

\[[F_{(\beta,\alpha)}=] \,F \, \equiv \int_{0}^{l_{\rm max}} f(l,\phi) \,{\rm d}l\]

Note that tackling the integration over the range \([l_{min}-l_{max}]\) instead (for instance for a clump), implies more if/else statement in the code, but are irrelevant for the discussion below. For the sake of brevity, we define the quantities (see also figure below):

  • l.o.s.: line of sight;
  • l.o.i.: line of integration (for a given \(\alpha,\beta\) direction around the l.o.s.);
  • F.O.I.: field of integration (usually corresponds to the solid angle defined by the instrument resolution);
  • \(r_{\rm trick}\): largest radius (measured from the clump centre) for which \(\rho(r<r_{\rm trick})\) still falls in the F.O.I.
  • Special case 1: If \(\phi_{l.o.i.}>\pi/2\) (where \(\phi_{l.o.i.}\) is the angle between the GC and the l.o.i.) and if this is not a clump, \(\rho^2\) monotonically decreases with \(l\). A mere log-step integration suffices:

    \[F = \int_{\epsilon}^{l_{\rm max}} l\;\times\; f(l,\phi) \,{\rm d}\ln l\]
  • Special case 2: If \(\phi_{l.o.i.}<\pi/2\) (or if we integrate a clump anywhere in the Galaxy), \(\rho^2\) is no longer monotonic. Based on a simple geometric argument (halos are quasi-sphericals), the density peaks, for the Galactic halo, at \(l_{\rho_{\rm max}} = \cos(\phi)/R_\odot\). Hence, the integral can be broken down in two parts, both of them being now monotonic. An offset applied to \(l\) allows to take advantage of the Simpson log-step integral (as the integrand is a almost power-law decreasing function from both sides of \(l_{\rho_{\rm max}}\)):

    \[F = I_1 + I_2\,,\]


    \[I_1 = \int_{\epsilon}^{l_{\rm max}-l_{\rho_{\rm max}}} l'\;\times\; f_{\rm of\!fset}(l',\phi) \,{\rm d}\ln l' \quad {\rm having} \quad l' = l_{\rho_{\rm max}}-l\]


    \[I_2 = \int_{\epsilon}^{l_{\rho_{\rm max}}} l'\;\times\; f_{\rm of\!fset}(l',\phi) \,{\rm d}\ln l' \quad {\rm having} \quad l' = l-l_{\rho_{\rm max}}\,.\]

    This simple and efficient strategy is explained in the Fig. 6.4 below (it is for a clump not coincident with the GC, but the spirit is the same would the DM smooth component only be considered). The two-part split at \(l_{\rho_{\rm max}}\) is seen on the l.o.i. in dotted blue on the left:


    Fig. 6.4 Integration strategy along the l.o.s. Click to enlarge the figure.

Remark 1: If the integration region encompasses the position of the DM density peak, we treat the central part as a point-like contribution. We define \(r_{\rm trick}\) to be the maximum radius for which we are still in the solid angle of the integration and for which the point-like criterion is met (i.e. \(r_{\rm trick}\leq f \;\times; l_{\rm cl}\), where \(f\sim 10^{-3}\), in which case \(l\) can be considered constant and factors out of the integration). The rest is a question of geometry to get correctly the various quantities to test and skip the region within \(r_{\rm trick}\) to avoid double-counting.

Remark 2: If the field of integration crosses the clump centre (this is rare, but not impossible), the Simpson routine again is stuck. There is no simple and clean fix to this problem. For the time being we simply integrate on a slightly smaller field of integration to get rid of this stiff point.

Back to \(\alpha\) and \(\beta\) for a clump out of the l.o.s.

When the l.o.s. falls completely out of the clump (but chunks of it are still encompassed in the integration angle), most of the integration angles \(\alpha\) and \(\beta\) contribute to zero and the simpson routine gets stuck. Fig. 6.5 shows this situation (projected on the sky).


Fig. 6.5 Integration strategy for a clump out of the l.o.s.: face-on view

  • Strategy for \(\alpha\):

    The centre of the F.O.I. corresponds to \(\alpha=0\). We replace this lower integration boundary by \(\alpha_{\rm min} = \psi_{l.o.s.,\,\rm cl_{centre}} - \alpha_{\rm cl}\) where \(\alpha_{\rm cl} = \sin^{-1}(R_{\rm cl}/l_{\rm cl})\) is the angular size of the clump in the sky (this obviously equals 0 is the centre of the F.O.I. touches exactly the clump boundary).

  • Strategy for \(\beta\):

    We calculate the minimal integration range required \(\Delta\beta\), which from Fig. 6.5, is given by

    \[\Delta\beta = \min(\beta_{\rm view},\beta_{\rm intersect})\,.\]

    On the right-hand side, the first angle is the half angle which encompasses the clump as seen from a point at a distance \(l_{\rm cl}\) on the l.o.s.. This angle is indeed maximal when the distance on the l.o.s. equals the distance of the clump, i.e. when the distance between the l.o.s. and the clump-centre, denoted \(d_{(l.o.s.)-(\rm cl)}\) is minimal:

    \[\beta_{\rm view} = \sin^{-1}(R_{\rm cl}/d_{(l.o.s.)-(\rm cl)})\,.\]

    The second angle is where the F.O.I. intersects the clump boundary. It is given by the Al-Kashi theorem

    \[\beta_{\rm intersect} = \cos^{-1}\left( \frac{R_{\rm cl}^2 + R_{\rm F.O.I.}^2 - d_{(l.o.s.)-(\rm cl)}^2}{2 \; R_{\rm cl} \times R_{\rm F.O.I.}} \right)\]

    where \(R_{\rm F.O.I.} = \sin(\alpha_{\rm int}) \times l_{\rm cl}\).

    We finally have (symbolic notation)

    \[\int_{0}^{2\pi} F_\beta\, {\rm d}\beta = 2 \int_{\beta_{\rm ref}}^{\beta_{\rm ref}+\beta_{\rm integr}} F_\beta\, {\rm d}\beta = 2 \int_{\beta_{\rm ref}-\beta_{\rm integr}}^{\beta_{\rm ref}} F_\beta\, {\rm d}\beta\,,\]

    where \(\beta_{\rm ref}\) is the angle between \(\beta=0\) (i.e. the \(z'\) axis, see Fig. 6.2) and the centre of the clump.

Fig. 6.6 finally shows a gprof2dot caller flowchart for the -h5 module example from the Quick start tutorial (for every called function, it is shown the time spent in this function and all its children, in parentheses the time spent in the function alone, and the total number of function calls). Unsurprisingly, most of the run time is spent in integrating (\(\tt simpson\_log\_adapt\) and \(\tt trapeze\_log\_refine\), see integr.h for more details about CLUMPY’s integrators). Note how these integrals, being mostly line-of-sight integrals, are called by the integr_los.h core function \(\tt los\_integral\_mix()\) optimising the calculation.

Fig. 6.6 Caller flowchart