Spergel shoot photon for nu<0

[jecampagne]

Hello, Well I'm looking at the SBSpergel.cpp::SpergelInfo::shoot method for nu<0.

Let me write some formula to a little more precise in the following. The Spergel profile noted I_{nu}(r) = N_{Nu} * f(r/r0,nu) is such that (N_{nu} is a normalisation factor, see later)

$$\Large f(z,\nu)=z^\nu K_\nu(z)$$ 

as the following property for z->0 (or r->0 for I_{nu}(r))

$$
\begin{cases}\Large
 f(z,\nu>0) = 2^{\nu -1} \Gamma(\nu)\\
\Large f(z,\nu=0) \propto \log(1/z) \\
\Large f(z,\nu<0) \propto 1/z^{-2\nu}
\end{cases}
$$

So, I agree that I_{nu}(r) is diverging for nu<=0. But the cumulated integrated flux is proportional to

$$\Large \int_0^z  f(u,\nu) u du$$

so has no divergence at z=0 if nu>-1 (which is the original paper bound) and Galsim uses nu is in [-0.85, 4.0] so in principle there is no theoretical problem to generate a distribution with the cumulative flux function.

But, let say the current implementation of OneDimensionalDeviatesampling that is a general purpose tool needs to have something that f(r) is non diverging. Let us have a look at the method used in the SpergelInfo::shootfunction. The main idea is to replace the I_{nu}(r) expression by a linear one for r in [0, rmin] imposing a flux conservation. Now, rmin/r0 = zmin (with r0 the scale radius in the Galsim language) is defined such

\Large \int_0^{zmin }  N_{\nu} f(z,\nu)\ 2\pi z dz = flux\_target

This is coded as

double flux_target = _gsparams->shoot_accuracy;
double shoot_rmin = calculateFluxRadius(flux_target);

fine.

Now considering the linear approximation mentioned above, it is written (below I have replaced "r" by "z" as this is the reduced radius r/r0)

                // int(2 pi z (a + b z) dz, 0..zmin) = shoot_accuracy (= flux_target)
                // a + b zmin = K_nu(zmin) * zmin^nu

and I verify that it is what is really used to compute the (a,b)parameters as solutions of a 2eq-2var system. These paramters are used to defined a modified profile (ie. SpergelNuNegativeRadialFunction) such that

\Large \tilde{I}(z) = \begin{cases}
a+b z & z \leq zmin \\
f(z,\nu) & z> zmin 
\end{cases}

But, there is for me a problem: if the (a,b) are used to define \tilde{I}(z) then I would expect

\Large  \int_0^{zmin} (a + b z)\ 2\pi z dz = flux\_target / N_\nu

And then it leads to a different 2eq-2var system, leading to different numerical values.

So there is a missmatch of definition and as (a,b)are defined & computed according to the first one, then I suspect a problem in the definition of SpergelNuNegativeRadialFunction and in the total flux normalization.

Am I wrong in the math somewhere?

[rmjarvis]

Sorry, I didn't follow this. What do you think is being computed wrong? The point of the whole rmin, a, b thing is just so that the right number of photons get generated near r=0. Do you think that doesn't come out to the right number? Or that they are distributed wrongly within r<rmin? (Or something else?) If the latter, then it's definitionally ok, since rmin is determined by shoot_accuracy -- the level at which we are ok with a little bit of inaccuracy.

[jecampagne]

I have coded the linear approx with the normalization of (a,b) mentioned in the last post, and the pytest is passed. So one can have a linear approx with a correct expression that is ok.l