Zeeman splitting (broadening?)

.gitignore

@@ -17,6 +17,7 @@ misc/d.pdf
 misc/*.png
 docs/build/
 docs/src/index.md
+docs/src/generated/
 *.DS_Store
 results_*.json
 results_*.json.tmp

docs/.gitignore

@@ -0,0 +1 @@
+.CondaPkg

src/constants.jl

@@ -8,6 +8,8 @@ const electron_charge_cgs = 4.80320425e-10 # statcoulomb or cm^3/2 * g^1/2 / s
 const amu_cgs = 1.6605402e-24 #g
 const bohr_radius_cgs = 5.29177210903e-9 # cm 2018 CODATA recommended value
 
+const bohr_magneton_cgs = 9.2740100783e-21 # erg G
+
 const solar_mass_cgs = 1.9884e33  #g
 const G_cgs = 6.67430e-8
 

src/line_absorption.jl

@@ -27,8 +27,9 @@ other arguments:
     is multithreaded over the lines in `linelist`.
   - `verbose` (deprecated): no longer used.
 """
-function line_absorption!(α, linelist, λs::Wavelengths, temps, nₑ, n_densities, partition_fns, ξ,
-                          α_cntm; cutoff_threshold=3e-4, verbose=false, tasks_per_thread=1)
+function line_absorption!(α, linelist, lande_g_factors, magnetic_field, λs::Wavelengths, temps,
+                          nₑ, n_densities, partition_fns, ξ, α_cntm; cutoff_threshold=3e-4,
+                          verbose=false, tasks_per_thread=1)
     if length(linelist) == 0
         return zeros(length(λs))
     end
@@ -46,6 +47,7 @@ function line_absorption!(α, linelist, λs::Wavelengths, temps, nₑ, n_densiti
     n_chunks = tasks_per_thread * Threads.nthreads()
     chunk_size = max(1, length(linelist) ÷ n_chunks + (length(linelist) % n_chunks > 0))
     linelist_chunks = partition(linelist, chunk_size)
+
     tasks = map(linelist_chunks) do linelist_chunk
         # Each chunk of your data gets its own spawned task that does its own local, sequential work
         # and then returns the result
@@ -62,6 +64,8 @@ function line_absorption!(α, linelist, λs::Wavelengths, temps, nₑ, n_densiti
             levels_factor = Vector{eltype(α)}(undef, size(temps))
             ρ_crit = Vector{eltype(α)}(undef, size(temps))
             inverse_densities = Vector{eltype(α)}(undef, size(temps))
+            ΔE_zeeman = Vector{eltype(α)}(undef, size(temps))
+            Δλ_zeeman = Vector{eltype(α)}(undef, size(temps))
 
             for line in linelist_chunk
                 m = get_mass(line.species)
@@ -94,15 +98,43 @@ function line_absorption!(α, linelist, λs::Wavelengths, temps, nₑ, n_densiti
                 doppler_line_window = maximum(inverse_densities)
                 inverse_densities .= inverse_lorentz_density.(ρ_crit, γ)
                 lorentz_line_window = maximum(inverse_densities)
-                window_size = sqrt(lorentz_line_window^2 + doppler_line_window^2)
-                lb = searchsortedfirst(λs, line.wl - window_size)
-                ub = searchsortedlast(λs, line.wl + window_size)
-                # not necessary, but is faster as of 8f979cc2c28f45cd7230d9ee31fbfb5a5164eb1d
-                if lb > ub
-                    continue
+
+                # calculate using LS approximation
+                if !ismissing(line.J_even) && !ismissing(line.J_odd)
+                    d = line.J_even * (line.J_even + 1) - line.J_odd * (line.J_odd + 1)
+                    LS_g_eff = (line.lande_g_even + line.lande_g_odd) / 2 +
+                               (line.lande_g_even - line.lande_g_odd) / 4 * d
+                    if LS_g_eff == 0.0
+                        Δλ_zeeman .= 0.0
+                    else
+                        ΔE_zeeman .= bohr_magneton_cgs * LS_g_eff * magnetic_field
+                        # TODO upper or lower E
+                        Δλ_zeeman .= @. line.wl^2 / (c_cgs * hplanck_cgs) * ΔE_zeeman
+                        ∂λ_∂E = -line.wl^2 / (c_cgs * hplanck_cgs)
+                        Δλ_zeeman = ∂λ_∂E .* ΔE_zeeman
+                    end
+                else
+                    Δλ_zeeman .= 0.0
                 end
 
-                α_task[:, lb:ub] .+= line_profile.(line.wl, σ, γ, amplitude, view(λs, lb:ub)')
+                window_size = sqrt(lorentz_line_window^2 + doppler_line_window^2) +
+                              maximum(Δλ_zeeman) - minimum(Δλ_zeeman) # make sure it's big enough to include all Zeeman shifts
+                for shift in (-Δλ_zeeman, 0, Δλ_zeeman) # TODO selection rules
+                    λ = @. line.wl + shift
+
+                    # calculate the window center from the nominal line center, not the shifted line
+                    # centers
+                    # at present, this line is allocating. Would be good to fix that.
+                    lb = searchsortedfirst(λs, line.wl - window_size)
+                    ub = searchsortedlast(λs, line.wl + window_size)
+                    # not necessary, but is faster as of 8f979cc2c28f45cd7230d9ee31fbfb5a5164eb1d
+                    if lb > ub
+                        continue
+                    end
+
+                    # TODO selection rules scaling
+                    α_task[:, lb:ub] .+= line_profile.(λ, σ, γ, amplitude, view(λs, lb:ub)')
+                end
             end
             return α_task
         end

src/linelist.jl

@@ -14,10 +14,18 @@ struct Line{F1,F2,F3,F4,F5,F6}
     # ABO theory
     vdW::Tuple{F6,F6}
 
+    # TODO
+    lande_g_even # Landé g factor for even level
+    lande_g_odd  # Landé g factor for odd level
+    J_even
+    J_odd
+
     @doc """
         Line(wl::F, log_gf::F, species::Species, E_lower::F,
              gamma_rad::Union{F, Missing}=missing, gamma_stark::Union{F, Missing}=missing,
-             vdw::Union{F, Tuple{F, F}, Missing}, missing) where F <: Real
+             vdw::Union{F, Tuple{F, F}, Missing}=missing,
+             lande_g_even::Union{F, Missing}=missing, lande_g_odd::Union{F, Missing}=missing
+        ) where F <: Real
 
     Arguments:
      - `wl`: wavelength (Assumed to be in cm if < 1, otherwise in Å)
@@ -34,6 +42,8 @@ struct Line{F1,F2,F3,F4,F5,F6}
          - A fudge factor for the Unsoeld approximation, assumed if between 0 and 20
          - The [ABO parameters](https://github.com/barklem/public-data/tree/master/broadening-howto)
            as packed float (assumed if >= 20) or a `Tuple`, `(σ, α)`.
+     - `lande_g_even`: Landé g factor for even level (optional)
+     - `lande_g_odd`: Landé g factor for odd level (optional)
 
         This behavior is intended to mirror that of Turbospectrum as closely as possible.
 
@@ -51,8 +61,10 @@ struct Line{F1,F2,F3,F4,F5,F6}
     """
     function Line(wl::F1, log_gf::F2, species::Species, E_lower::F3,
                   gamma_rad::Union{F4,Missing}=missing, gamma_stark::Union{F5,Missing}=missing,
-                  vdW::Union{F6,Tuple{F6,F6},Missing}=missing) where {F1<:Real,F2<:Real,F3<:Real,
-                                                                      F4<:Real,F5<:Real,F6<:Real}
+                  vdW::Union{F6,Tuple{F6,F6},Missing}=missing, lande_g_even=missing,
+                  lande_g_odd=missing, J_even=missing,
+                  J_odd=missing) where {F1<:Real,F2<:Real,F3<:Real,F4<:Real,F5<:Real,
+                                        F6<:Real}
         if wl >= 1
             wl *= 1e-8 #convert Å to cm
         end
@@ -92,14 +104,19 @@ struct Line{F1,F2,F3,F4,F5,F6}
 
         new{F1,F2,F3,typeof(gamma_rad),typeof(gamma_stark),eltype(vdW)}(wl, log_gf, species,
                                                                         E_lower, gamma_rad,
-                                                                        gamma_stark, vdW)
+                                                                        gamma_stark, vdW,
+                                                                        lande_g_even,
+                                                                        lande_g_odd,
+                                                                        J_even, J_odd)
     end
 end
 # constructor to allow for copying a line and modifying some values (see docstring)
 function Line(line::Line; wl=line.wl, log_gf=line.log_gf, species=line.species,
               E_lower=line.E_lower, gamma_rad=line.gamma_rad, gamma_stark=line.gamma_stark,
-              vdW=line.vdW)
-    Line(wl, log_gf, species, E_lower, gamma_rad, gamma_stark, vdW)
+              vdW=line.vdW, lande_g_odd=line.lande_g_odd, lande_g_even=line.lande_g_even,
+              J_even=line.J_even, J_odd=line.J_odd)
+    Line(wl, log_gf, species, E_lower, gamma_rad, gamma_stark, vdW, lande_g_even, lande_g_odd,
+         J_even, J_odd)
 end
 
 function Base.show(io::IO, ::MIME"text/plain", line::Line)
@@ -383,6 +400,7 @@ function parse_kurucz_linelist(f; isotopic_abundances=nothing, vacuum=false, ver
         row == "" && continue #skip empty lines
 
         #some linelists have a missing column in the wavelenth region
+        # TODO make this more robust!  I think this was the problem!
         if length(row) == 159
             row = " " * row
         end
@@ -394,6 +412,15 @@ function parse_kurucz_linelist(f; isotopic_abundances=nothing, vacuum=false, ver
             abs(parse(Float64, s)) * c_cgs * hplanck_eV
         end
 
+        # TODO I HAD TO EDIT gf0640.10 !!!
+        # I inserted a space at the begining of every line to make it match the spec
+
+        # Parse lande g factors (columns 28 and 29, chars 141:145 and 146:150)
+        lande_g_even = parse(Float64, strip(row[145:149])) / 1000
+        lande_g_odd = parse(Int64, strip(row[150:154])) / 1000
+
+        quantum_g_even = parse(Float64, strip(row[38:41]))
+        quantum_g_odd = parse(Float64, strip(row[66:69]))
         species = Species(row[19:24])
 
         # log(gf) + adjustment for hyperfine structure
@@ -428,7 +455,11 @@ function parse_kurucz_linelist(f; isotopic_abundances=nothing, vacuum=false, ver
                    min(E_levels...),
                    tentotheOrMissing(parse(Float64, row[81:86])),
                    tentotheOrMissing(parse(Float64, row[87:92])),
-                   idOrMissing(parse(Float64, row[93:98]))))
+                   idOrMissing(parse(Float64, row[93:98])),
+                   lande_g_even,
+                   lande_g_odd,
+                   quantum_g_even,
+                   quantum_g_odd))
     end
     lines
 end
@@ -760,6 +791,45 @@ function read_korg_linelist(path)
     end
 end
 
+"""
+    get_lande_g_factor(filename)
+
+Get the Landé g-factors for lines from a VALD long-format file. Returns a vector of tuples
+(`mean`, `lower`, `upper`)  containing the Landé g-factors for each line.
+
+# Arguments
+
+  - `filename`: Path to the VALD long-format file
+
+# Returns
+
+  - Vector of named tuples with fields `lower`, `upper`, and `mean` containing the Landé g-factors
+"""
+function get_lande_g_factor(filename, nlines)
+    lines = filter!(collect(eachline(filename))) do line
+        length(line) > 0 && line[1] != '#' # remove comments and empty lines
+    end
+
+    # Find the header line that contains "Lande factors"
+    header_idx = findfirst(line -> occursin("Lande factors", line), lines)
+    isnothing(header_idx) && error("Could not find Landé factors header in file")
+
+    body = lines[header_idx+2:4:header_idx+2+4*(nlines-1)]
+
+    # Parse each data line to extract Landé g-factors
+    map(body) do line
+        # Split by commas and remove quotes
+        parts = split(strip(line, ['\'', ' ']), ',')
+        # Extract the three Landé g-factors (they are in positions 8, 9, and 10)
+        lower = parse(Float64, strip(parts[8]))
+        upper = parse(Float64, strip(parts[9]))
+        mean = parse(Float64, strip(parts[10]))
+        (mask_99(mean), mask_99(lower), mask_99(upper))
+    end
+end
+
+mask_99(x) = x == 99.0 ? 0.0 : x
+
 """
     get_VALD_solar_linelist()
 

src/synthesize.jl

@@ -132,6 +132,8 @@ result = synthesize(atm, linelist, A_X, 5000, 5100)
 """
 function synthesize(atm::ModelAtmosphere, linelist, A_X::AbstractVector{<:Real},
                     wavelength_params...;
+                    lande_g_factors=nothing,
+                    magnetic_field=nothing,
                     vmic=1.0,
                     line_buffer::Real=10.0,
                     cntm_step::Real=1.0,
@@ -162,6 +164,10 @@ function synthesize(atm::ModelAtmosphere, linelist, A_X::AbstractVector{<:Real},
         @warn "if you are synthesizing at wavelengths longer than 15000 Å (e.g. for APOGEE), setting use_MHD_for_hydrogen_lines=false is recommended for the most accurate synthetic spectra. This behavior may become the default in Korg 1.0."
     end
 
+    if I_scheme == "linear_flux_only" && mu_values != 20
+        @warn "mu_values was set to a non-default value ($mu_values), which will have no effect when I_scheme is set to \"linear_flux_only\"."
+    end
+
     # Add wavelength bounds check (Rayleigh scattering limitation)
     # we should really have an upper bound as well
     min_allowed_wavelength = 1300.0 * 1e-8  # cm
@@ -248,9 +254,8 @@ function synthesize(atm::ModelAtmosphere, linelist, A_X::AbstractVector{<:Real},
     # line contributions to α5
     if tau_scheme == "anchored"
         α_cntm_5 = [_ -> a for a in copy(α_ref)] # lambda per layer
-        line_absorption!(view(α_ref, :, 1), linelist5, Korg.Wavelengths([5000]), get_temps(atm),
-                         nₑs,
-                         number_densities,
+        line_absorption!(view(α_ref, :, 1), linelist5, lande_g_factors, magnetic_field,
+                         Wavelengths([5000]), get_temps(atm), nₑs, number_densities,
                          partition_funcs, vmic * 1e5, α_cntm_5;
                          cutoff_threshold=line_cutoff_threshold)
         interpolate_molecular_cross_sections!(view(α_ref, :, 1), molecular_cross_sections,
@@ -278,8 +283,9 @@ function synthesize(atm::ModelAtmosphere, linelist, A_X::AbstractVector{<:Real},
         end
     end
 
-    line_absorption!(α, linelist, wls, get_temps(atm), nₑs, number_densities, partition_funcs,
-                     vmic * 1e5, α_cntm; cutoff_threshold=line_cutoff_threshold, verbose=verbose)
+    line_absorption!(α, linelist, lande_g_factors, magnetic_field, wls, get_temps(atm), nₑs,
+                     number_densities, partition_funcs, vmic * 1e5, α_cntm;
+                     cutoff_threshold=line_cutoff_threshold, verbose=verbose)
     interpolate_molecular_cross_sections!(α, molecular_cross_sections, wls, get_temps(atm), vmic,
                                           number_densities)