zblz / naima

Derivation of non-thermal particle distributions through MCMC spectral fitting
http://naima.readthedocs.io
BSD 3-Clause "New" or "Revised" License
45 stars 54 forks source link

Add synchrotron emission in turbulent magnetic field #135

Open tibaldo opened 7 years ago

tibaldo commented 7 years ago

See example Kelner, Aharonian, Khangulyan 2013 ApJ 774:61 Sec 4. Options could be power-law turbulence spectrum or tabulated values for the magnetic field PDF. I will branch and give it a go.

zblz commented 7 years ago

@tibaldo - that would be quite useful! I already have working code that implements Kelner et al. 2013, but its is quite slow at the moment. If you want to give it a go you could use it as a starting point: https://gist.github.com/zblz/35c9a635f710f20428eff0ca90a24a2c

This could also be a good opportunity to see how we can refactor the Synchrotron class so both this modification and proton synchrotron #134 can be included without much code duplication.

tibaldo commented 7 years ago

@zblz the code was very useful! I managed to reduce the time for the calculation of the spectrum by ~40% by invoking cbrt only once in Gtilda

            cb = cbrt(x)
            gt1 = 1.808 * cb / np.sqrt(1 + 3.4 * cb ** 2.)
            gt2 = 1 + 2.210 * cb ** 2. + 0.347 * cb ** 4.
            gt3 = 1 + 1.353 * cb ** 2. + 0.217 * cb ** 4.
            return gt1 * (gt2 / gt3) * np.exp(-x)

Original code (based on your snippet)

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
      700   98.579    0.141   98.579    0.141 magnetic_dist_3.py:249(Gtilde)

modified code with the tweak described above

     700   56.143    0.080   56.143    0.080 magnetic_dist_4.py:123(Gtilde)

Do you feel like this makes it fast enough? I did not find any other places where the calculation can be made faster.

zblz commented 7 years ago

That's brilliant, should probably add that modification to the regular Synchrotron class.

It will probably never be a super fast calculation, as it has to do quite a lot of work anyway, but I was thinking being a bit more clever about the integration over the magnetic field spectrum could reduce the number of single-B synchrotron computations needed.

tibaldo commented 7 years ago

@zblz Okay, I'll look more into the integration scheme, so far I have found that a large-number of single B computations (> 100) is needed for a power-law B distribution in order to get accurate results. This, and the overall refactoring of the class, may take some time to be done. Do you want me to send a pull request for the speeding of Gtilde now?

zblz commented 7 years ago

Yes, doing a first PR with the Gtilde speedup would be great.

With the code from the snippet my experience was similar, at least a 100 points in the B spectrum need to be evaluated to avoid the final spectrum being too different from an analytical calculation. However, each of these 100 spectra is very similar to the adjacent ones, so it feels like some interpolation should be possible before the integration if a lower number of single-B spectra are computed. I have not actually tested any of this so it might be completely wrong.