Merge pull request #147 from fjebaker/fergus/brems

BremsStrahlung Julia implementation

Project.toml

@@ -20,6 +20,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -31,3 +32,5 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Surrogates = "6fc51010-71bc-11e9-0e15-a3fcc6593c49"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
+[compat]
+Polynomials = "4.0.12"

docs/src/models/models.md

@@ -13,10 +13,11 @@ Order = [:type]
 ### Additive
 
 ```@docs
-PowerLaw
 BlackBody
-GaussianLine
+BremsStrahlung
 DeltaLine
+GaussianLine
+PowerLaw
 ```
 
 ### Multiplicative

src/SpectralFitting.jl

@@ -18,6 +18,7 @@ using FileIO
 using Interpolations
 import DataInterpolations
 using SpecialFunctions
+using Polynomials
 using PreallocationTools
 using EnumX
 

src/julia-models/additive.jl

@@ -110,6 +110,152 @@ end
     end
 end
 
+"""
+    BremsStrahlung
+
+$(FIELDS)
+
+# Example
+
+```julia
+energy = collect(range(0.1, 20.0, 100))
+invokemodel(energy, BremsStrahlung())
+```
+
+```
+                   BremsStrahlung               
+     ┌────────────────────────────────────────┐ 
+   3 │.                                       │ 
+     │:                                       │ 
+     │:                                       │ 
+     │:                                       │ 
+     │:                                       │ 
+     │:                                       │ 
+     │:                                       │ 
+     │:                                       │ 
+     │:                                       │ 
+     │:                                       │ 
+     │:                                       │ 
+     │:                                       │ 
+     │:                                       │ 
+     │ :                                      │ 
+   0 │ ':.....................................│ 
+     └────────────────────────────────────────┘ 
+      0                                     20  
+                       E (keV)                  
+```
+"""
+struct BremsStrahlung{T} <: AbstractSpectralModel{T,Additive}
+    "Normalisation."
+    K::T
+    "Temperature (keV)."
+    kT::T
+    "Helium to hydrogen ratio"
+    ab::T
+end
+function BremsStrahlung(; K = FitParam(1.0), kT = FitParam(7.0), ab = FitParam(0.085))
+    BremsStrahlung{typeof(K)}(K, kT, ab)
+end
+@inline function invoke!(flux, energy, model::BremsStrahlung)
+    let kT = model.kT, ab = model.ab
+        integration_kernel!(flux, energy) do E, δE
+            gaunt = (
+                _Internal_BremsStrahlung.gaunt(E, kT, 1) +
+                4 * ab * _Internal_BremsStrahlung.gaunt(E, kT, 2)
+            )
+            gaunt * exp(-E / kT) * δE / (E * sqrt(kT))
+        end
+    end
+end
+
+# scope these implementation specifics in a module so they don't collide in the
+# namespace
+module _Internal_BremsStrahlung
+
+using ..Polynomials
+
+const A1 = [
+    1.001; 1.004; 1.017; 1.036; 1.056; 1.121;; 1.001; 1.005; 1.017; 1.046; 1.073; 1.115;; 0.9991; 1.005; 1.03; 1.055; 1.102; 1.176;; 0.997; 1.005; 1.035; 1.069; 1.134; 1.186;; 0.9962; 1.004; 1.042; 1.1; 1.193; 1.306;; 0.9874; 0.9962; 1.047; 1.156; 1.327; 1.485;; 0.9681; 0.9755; 1.020; 1.208; 1.525; 1.965;;;
+    0.3029; 0.1616; 0.04757; 0.013; 0.0049; -0.0032;; 0.4905; 0.2155; 0.08357; 0.02041; 0.00739; 0.00029;; 0.654; 0.2833; 0.08057; 0.03257; 0.00759; -0.00151;; 1.029; 0.391; 0.1266; 0.05149; 0.01274; 0.00324;; 0.9569; 0.4891; 0.1764; 0.05914; 0.01407; -0.00024;; 1.236; 0.7579; 0.326; 0.1077; 0.028; 0.00548;; 1.327; 1.017; 0.6017; 0.205; 0.0605; 0.00187;;;
+    -1.323; -0.254; -0.01571; -0.001; -0.000184; 0.00008;; -4.762; -0.3386; -0.03571; -0.001786; -0.0003; 0.00001;; -6.349; -0.4206; -0.02571; -0.003429; -0.000234; 0.00005;; -13.231; -0.59; -0.04571; -0.005714; -0.000445; -0.00004;; -7.672; -0.6852; -0.0643; -0.005857; -0.00042; 0.00004;; -7.143; -0.9947; -0.12; -0.01007; -0.000851; -0.00004;; -3.175; -1.116; -0.2270; -0.01821; -0.001729; 0.00023
+]
+
+const γ2 = [0.7783, 1.2217, 2.6234, 4.3766, 20.0, 70.0]
+const γ3 = [1.0, 1.7783, 3.0, 5.6234, 10.0, 30.0]
+const Alo1 =
+    [-0.57721566, 0.4227842, 0.23069756, 0.0348859, 0.00262698, 0.0001075, 0.0000074]
+const Alo2 = [1.0, 3.5156229, 3.089942, 1.2067492, 0.2659732, 0.0360768, 0.0045813]
+const Ahi1 =
+    [1.25331414, 0.07832358, 0.02189568, 0.01062446, 0.00587872, 0.0025154, 0.00053208]
+
+function gaunt(E::T, kT::T, z::Integer) where {T<:Real}
+    γ = 0.01358 * z^2 / kT
+    if kT == 0 || E > 50 * kT || E == 0
+        #  Temperature or Energy is out of range
+        gaunt = zero(T)
+    elseif γ > 0.1
+        #  Low kT regime, use Kurucz's algorithm
+        gaunt = kurucz(E / kT, γ)
+    else
+        #  Calculate Born approximation
+        u = E / kT / 4
+        born =
+            0.5513 *
+            exp(2u) *
+            (
+                u <= 1 ? Polynomial(Alo1)(u^2) - log(u) * Polynomial(Alo2)((4u / 7.5)^2) :
+                Polynomial(Ahi1)(-1 / u) / exp(2u) / sqrt(2u)
+            )
+        if min(1000 * γ, 100) < 1
+            #  Use Born approximation if valid
+            gaunt = born
+        else
+            u, γ1 = max(u, 0.003), min(1000 * γ, 100.0)
+            n, m = N(u), M(γ1)
+            p = (γ1 - γ3[m]) / γ2[m]
+            G1 = [A1[n, m, 1], A1[n, m, 2], A1[n, m, 3]]
+            G2 = [A1[n, m+1, 1], A1[n, m+1, 2], A1[n, m+1, 3]]
+            gaunt = ((1.0 - p) * Polynomial(G1)(u) + p * Polynomial(G2)(u)) * born
+        end
+    end
+    return gaunt
+end
+
+N(x) = x <= 0.03 ? 1 : x <= 0.3 ? 2 : x <= 1.0 ? 3 : x <= 5.0 ? 4 : x <= 15.0 ? 5 : 6
+M(x) = x <= 1.773 ? 1 : x <= 3.0 ? 2 : x <= 5.6234 ? 3 : x <= 10.0 ? 4 : x <= 30.0 ? 5 : 6
+
+const A2 = [
+    5.40; 5.25; 5.00; 4.69; 4.48; 4.16; 3.85;;
+    4.77; 4.63; 4.40; 4.13; 3.87; 3.52; 3.27;;
+    4.15; 4.02; 3.80; 3.57; 3.27; 2.98; 2.70;;
+    3.54; 3.41; 3.22; 2.97; 2.70; 2.45; 2.20;;
+    2.94; 2.81; 2.65; 2.44; 2.21; 2.01; 1.81;;
+    2.41; 2.32; 2.19; 2.02; 1.84; 1.67; 1.50;;
+    1.95; 1.90; 1.80; 1.68; 1.52; 1.41; 1.30;;
+    1.55; 1.56; 1.51; 1.42; 1.33; 1.25; 1.17;;
+    1.17; 1.30; 1.32; 1.30; 1.20; 1.15; 1.11;;
+    0.86; 1.00; 1.15; 1.18; 1.15; 1.11; 1.08;;
+    0.59; 0.76; 0.97; 1.09; 1.13; 1.10; 1.08;;
+    0.38; 0.53; 0.76; 0.96; 1.08; 1.09; 1.09
+]
+
+tkur(γ::Real, j::Integer) = 2log10(γ) + 3 - j
+ukur(μ::Real, k::Integer) = 2log10(μ) + 9 - k
+
+function kurucz(μ::T, γ::T) where {T<:Real}
+    #  Algorithm for low kT
+    j = round(Integer, tkur(γ, 0))
+    k = max(1, round(Integer, ukur(μ, 0)))
+    t, u = tkur(γ, j), ukur(μ, k)
+    return j >= 7 || j < 1 || k >= 12 ? zero(T) :
+           (1 - t) * (1 - u) * A2[j, k] +
+           t * (1 - u) * A2[j+1, k] +
+           (1 - t) * u * A2[j, k+1] +
+           t * u * A2[j+1, k+1]
+end
+
+end # module _Internal_BremsStrahlung
+
 """
     GaussianLine
 
@@ -229,4 +375,4 @@ Reflection.get_parameter_symbols(model::Type{<:DeltaLine}) = fieldnames(model)[2
     invoke!(flux, energy, GaussianLine(promote(model.K, model.E, model._width)...))
 end
 
-export PowerLaw, BlackBody, GaussianLine, DeltaLine
+export PowerLaw, BlackBody, BremsStrahlung, GaussianLine, DeltaLine

test/models/test-julia-models.jl

@@ -3,7 +3,7 @@ using SpectralFitting
 
 include("../utils.jl")
 
-ALL_JULIA_MODELS = [PowerLaw, BlackBody]
+ALL_JULIA_MODELS = [PowerLaw, BlackBody, BremsStrahlung]
 
 # has data requirements, so skip on the CI
 @ciskip push!(ALL_JULIA_MODELS, PhotoelectricAbsorption)