astro-group-bristol / Gradus.jl

Extensible spacetime agnostic general relativistic ray-tracing (GRRT).
https://astro-group-bristol.github.io/Gradus.jl/dev/
GNU General Public License v3.0
18 stars 2 forks source link

Can we do transfer functions better? #106

Closed fjebaker closed 1 year ago

fjebaker commented 1 year ago

I'm really not happy with how the transfer functions are currently defined and used. Since the integration domain is over $g$ when calculating the line profiles anyway, the conversion to $g^\ast$ is really inconvenient, as it introduces a number of numerical issues when calculating the transfer functions themselves. Consider the definition

$$ f := \frac{g}{\pi r} \sqrt{g^\ast (1 - g^\ast)} \left\lvert \frac{\partial(\alpha \beta}{\partial(r, g^\ast)} \right\rvert, $$

which is then used to find the flux through

$$ \text{d}F = g^2 I \frac{\pi r}{ \sqrt{g^\ast (1 - g^\ast)} } f \text{d}r \text{d} g^\ast. $$

This is itself used in the integral

$$ F = \int_{r=r0}^{r=r\text{max}}\int_{g^\ast=0}^{g^\ast=1} \text{d}F. $$

It's pretty obvious that most of the terms cancel. Care when integrating comes from the regions where $g^\ast \rightarrow { 0,1 }$, but that is because by the transfer functions go to zero at those points, but $\text{d}F \rightarrow \infty$ numerically, even though there is an analytic value. The prescription is then to have a buffer $h > 0$, where special edge integrations are performed instead.

The transfer functions could far more conveniently be defined

$$ f := g\left\lvert \frac{\partial(\alpha \beta}{\partial(r, g)} \right\rvert, $$

such that

$$ \text{d}F = g^2 I f \text{d}r \text{d} g, $$

and the integration is then performed directly over $g$ instead of $g^\ast$, which also then elides the need for a $g \rightarrow g^\ast$ Jacobian.

Possible gripes can easily be fixed in this formalism too:

It seems the transfer functions could be so much easier to use when defined this way.

fjebaker commented 1 year ago

Something I failed to account for here is that the Jacobian has the behaviour that

$$ \lim{g\rightarrow g\text{min}, g_\text{max}}\left\lvert \frac{\partial(\alpha, \beta)}{\partial(r, g)} \right\rvert \rightarrow \infty, $$

which is supressed by including the $\sqrt{g^\ast ( 1 - g^\ast)}$ term. Reparameterising $\text{d} r$ and $\text{d} \theta$ on the disc would solve this problem, since when $\text{d} \alpha / \text{d} r$ (respectively $\beta$) diverges, the orthogonal derivative w.r.t. $\theta$ on the disc would go to zero, keeping the Jacobian finite.

Of course the issue with this is then having to store more information (angles in addition to redshift values), and sampling representative $\theta$ such that integrating over $\text{d} \theta$ still converges correctly.

fjebaker commented 1 year ago

Closing, since this is not really that pertinent anymore.