cherab / core

The core source repository for the Cherab project.
https://www.cherab.info
Other
44 stars 24 forks source link

Bremsstrahlung formula reference #409

Closed skuba31 closed 1 year ago

skuba31 commented 1 year ago

Currently, the reference for the bremsstrahlung formula is the Beurskens ITER LIDAR article, which uses an unexplained constant and is not in SI units. Looking into references I was able to find an article by Kadota, where the formula is also presented. However, without derivation or explanation and it is written in some form of CGS unit system as the vacuum permittivity is not present. This is the oldest English written reference. It is pointing to Kontinuierliche Spektren (in German), where I hit the limits of my language skills. My colleague recommended me Hutchinson's Principles of Plasma Diagnostics where the formula is derived and the constants are in SI units. Some information about the gaunt factor is also present. The downside is that it is written with respect to frequency. But this can be solved by substitution factor.

$$ P{\mathrm{ff}} \approx \int{f\mathrm{min}}^{f\mathrm{max}} j\mathrm{ff} (f) \mathrm{d}f = \int{\lambda\mathrm{max}}^{\lambda\mathrm{min}} - \frac{c}{\lambda^2} j\mathrm{ff}\left(\frac{c}{\lambda}\right) \mathrm{d}\lambda = \int{\lambda\mathrm{min}}^{\lambda\mathrm{max}} \frac{c}{\lambda^2} j\mathrm{ff}\left(\frac{c}{\lambda}\right) \mathrm{d}\lambda $$

I did simple test and the bremsstrahlung based on Hutchinson was 1.008937... times greater than that computed using current formula. This seems to be caused by the limited accuracy of the constant in Beurskens article $0.95\cdot10^{−19}$.

As the difference is less than a percent, I would propose to at least include Hutchinson as a reference to the documentation for better clarity and optionally implement the equation with pre_factor based on fundamental constants or change it to be based on Hutchinson's book.

vsnever commented 1 year ago

Hi, @skuba31. I totally agree. This factor, 0.95 • 10-19, should be replaced with an actual expression for accuracy and clarity, and included as a constant just like the EXP_FACTOR. The reference to Beurskens is no longer appropriate because we have already replaced their fit for the Gaunt factor with the much more general precalculated tables from Carson expanded by de Avillez & Breitschwerdt. The reference to Hutchinson is more appropriate. Can you make a PR?

skuba31 commented 1 year ago

Hi @vsnever, I can make the PR. Should I wait for #402 to be merged as it significantly changes the Bremsstrahlung class?

vsnever commented 1 year ago

Sorry for not replying earlier. Since I still need to add a test in #402, it will be easier for me to replace the factor with an expression in the same PR, so as not to rewrite the test later. I'll let you know when I make these changes so you can review them.

skuba31 commented 1 year ago

Hello Vlad, I pushed the class I made for testing should you want to have a look. The formulas and docstrings are modified. It is based on an older version of Cherab though. https://github.com/skuba31/cherab-core/blob/hutchinson/cherab/core/model/plasma/hutchinson.pyx

vsnever commented 1 year ago

Hello Vlad, I pushed the class I made for testing should you want to have a look. The formulas and docstrings are modified. It is based on an older version of Cherab though. https://github.com/skuba31/cherab-core/blob/hutchinson/cherab/core/model/plasma/hutchinson.pyx

Thanks! I replaced the coarse numerical constant with an expression and updated the reference in #402 following your implementation. I think we can give the equation for λ even if it's given for ν in the reference without describing the conversion. Also, in the future, users will be able to fill the local atomic data repository with custom data for Gaunt factors, so I don't think we should refer to de Avillies and Breitschwerdt's paper in the Bremsstrahlung class itself.

Could you take a look at #402 and check if it's ok?

skuba31 commented 1 year ago

I agree that the equation can be written with respect to wavelength, but I think that keeping at least a note about conversion would make it easier for a newcomer to understand. You are right about the de Avilles reference. It should not be in the class itself.

I checked the changes you propose. I have only one comment regarding the return statement.

return radiance / wvl

It feels confusing to have one more division on that line. I think readability would benefit if it would be put directly to prefactor calculation as that is the place where I would expect it based on the docstrings

pre_factor = BREMS_CONST * ni_gff_z2 * self.ne / (sqrt(self.te) * wvl * wvl)

or even reorder to better match the order in docstrings

pre_factor = BREMS_CONST / (sqrt(self.te) * wvl * wvl) * self.ne * ni_gff_z2
vsnever commented 1 year ago

Thank you, @skuba31, I agree. I changed the expression for pre_factor as you suggested and added a note about the ν→λ conversion in the docstring in the last commit.

skuba31 commented 1 year ago

Excellent, I think this issue can be closed now.