Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Moffat: maxK #1208

Closed
jecampagne opened this issue Jan 21, 2023 · 2 comments · Fixed by #1210
Closed

Moffat: maxK #1208

jecampagne opened this issue Jan 21, 2023 · 2 comments · Fixed by #1210
Labels
bug report Bug report numerics Involving details of the numerical algorithms for performing some calculation(s) optimization/performance Related to the speed and/or memory consumption of some aspect of the code
Milestone

Comments

@jecampagne
Copy link

jecampagne commented Jan 21, 2023

Hi,
In the SBMoffat.cpp : SBMoffat::SBMoffatImpl::maxK() , I have found this comment (in the case of trunc=0, aka no truncation)

                // f(k) = 4 K(beta-1,k) (k/2)^beta / Gamma(beta-1)
                //
                // The asymptotic formula for K(beta-1,k) is
                //     K(beta-1,k) ~= sqrt(pi/(2k)) exp(-k)
                //
                // So f(k) becomes
                //
                // f(k) ~= 2 sqrt(pi) (k/2)^(beta-1/2) exp(-k) / Gamma(beta-1)
                //
                // Solve for f(k) = maxk_threshold

Well, the point is that the code reflect this comment, so it rises me a question.

  1. First takes the Moffat profile as
    f(r) = C (1+(r/rd)^2)^(-beta)
  2. Let takes the non-truncated Moffat, then in F is the flux that the user provides, it yields that the constant C is
    C = F (beta-1)/( pi rd^2)
  3. Taking the 2D Fourier Tansform needs first to fix the same definition that leads to the Gaussian power spectrum. In Mathematica this leads to fix {a,b}={1,-1} for instance, then it yields
    f(k) = 2pi rd^2 C \int_0^\infty J_0(krd x) (1+x^2)^(-beta) xdx = (4 F)/(2^beta Gamma(beta-1)) * (k rd)^(beta-1) K[beta-1, k rd]

So as you notice besides some constant what is puzzling to me is the k power scaling in the comment ie k^beta which is not what I find ie k^(beta-1).

  1. Moreover when computing the maxK value, the asymptotic behavior of f(k) is

f(k) \appros F \sqrt{pi}/Gamma(\beta-1) e^{k rd} (k rd/2)^(beta-3/2)

and then ones needs to solve w/o the F factor
f(maxK) = maxk_threshold

BUT what is used in the code is
" 2 sqrt(pi) (k/2)^(beta-1/2) exp(-k) / Gamma(beta-1) = maxk_threshold"

missing the correct k-scaling (concerning the Flux factor F, the Gaussian case does not takes it onto account so I guess this is the definition of maxk_threshold)

Anyway I puzzled concerning the k-scaling, do I make a mistake or a misunderstanding of the code?

@rmjarvis
Copy link
Member

rmjarvis commented Feb 2, 2023

Thanks Jean-Eric. I finally had some time to look into this.

First, you're quite correct that that bit of code is missing a power of k^-1. However, even fixing that, this approximation, using the high k asymptotic behavior, seems not to be a very good one for finding maxk. Especially for larger beta values, where that asymptote doesn't kick in until much higher k than what we need for maxk.

So I ended up scrapping that and just doing a proper fit for maxk using the real formula with the Bessel K function. It's plenty fast enough for this purpose, so there's not really a big need to do an approximate calculation.

I'll have a PR for this shortly.

@rmjarvis rmjarvis added numerics Involving details of the numerical algorithms for performing some calculation(s) bug report Bug report optimization/performance Related to the speed and/or memory consumption of some aspect of the code labels Feb 2, 2023
@rmjarvis rmjarvis added this to the v2.5 milestone Feb 2, 2023
@rmjarvis rmjarvis linked a pull request Feb 2, 2023 that will close this issue
@jecampagne
Copy link
Author

Hi, I'm glad to make GalSim better :)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug report Bug report numerics Involving details of the numerical algorithms for performing some calculation(s) optimization/performance Related to the speed and/or memory consumption of some aspect of the code
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants