diff --git a/.gitignore b/.gitignore index d3ff13e00..b050b05b6 100644 --- a/.gitignore +++ b/.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 diff --git a/docs/.gitignore b/docs/.gitignore new file mode 100644 index 000000000..06e0f9666 --- /dev/null +++ b/docs/.gitignore @@ -0,0 +1 @@ +.CondaPkg diff --git a/src/constants.jl b/src/constants.jl index 41ded3d28..bddb34237 100644 --- a/src/constants.jl +++ b/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 diff --git a/src/line_absorption.jl b/src/line_absorption.jl index 7aac10f43..98b109d02 100644 --- a/src/line_absorption.jl +++ b/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 diff --git a/src/linelist.jl b/src/linelist.jl index f9efeb292..3e854ece4 100644 --- a/src/linelist.jl +++ b/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() diff --git a/src/synthesize.jl b/src/synthesize.jl index 8e869f680..d27b344b4 100644 --- a/src/synthesize.jl +++ b/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)