Download this page as a Jupyter notebook

Equivalent width fitting

This tutorial shows how to use Korg to fit the abundances of a star using equivalent widths. We'll use the data from Melendez 2014, which are available in the Table 1.dat file, and recreate some of the analysis from that paper.

First, we'll import the necessary packages. These examples use the PythonPlot package, which provides a nice Julia interface to Matplotlib, but of course you can plot however you like. We'll also use CSV and DataFrames to read in the data from the paper.

    using Korg, PythonPlot, CSV, DataFrames

Creating a custom linelist

Korg can read in linelists in a variety of formats, including the VALD, MOOG, Kurucz, and ExoMol formats, but sometimes you'll want to use one in another format, e.g. from a table in a paper. We'll use the data from Table 1 of Melendez 2014, which are available in the Table 1.dat file.

First, we read the data into a DataFrame called lines. This is a convenient way to work with tabular data.

lines = CSV.read("../../assets/Table 1.dat", DataFrame; skipto=25, delim=' ', ignorerepeated=true,
                 header=["wl", "species", "ExPot", "log_gf", "C6", "EW_18Sco", "EW_Sun"]);

Next, we convert the numbers in the species column to Korg.Species objects. The Korg.Species constructor takes a string and returns a Korg.Species object, and supports nearly all of the formats for specifying species found in the wild.

The . applies the function to each element of the series. This is called "broadcasting".

lines.species = Korg.Species.(lines.species)

323-element Vector{Korg.Species}:
 Fe I
 Fe I
 Fe I
 Fe I
 Fe I
 Fe I
 Fe I
 Fe I
 Fe I
 Fe I
 ⋮
 Gd II
 Gd II
 Dy II
 Dy II
 Dy II
 Dy II
 Dy II
 Dy II
 Yb II

Let's look at Fe lines only, and sort them by wavelength.

filter!(lines) do row
    #careful, get_atom throws an error when applied to a molecular species
    Korg.get_atom(row.species) == 26 # atomic number of Fe
end

sort!(lines, :wl)
lines[1:4, :] # look at the first few rows

To pass the lines to Korg, we need to turn each row of the lines DataFrame into a Korg.Line object. Korg will use reasonable defaults for the broadening parameters, but see the Korg.Line documentation for details on how to specify them if you need to do that.

linelist = Korg.Line.(lines.wl, # can be in either cm or Å (like these), but NOT nm
                      lines.log_gf,
                      lines.species, # needs to be a Korg.Species, which we handled in the cell above
                      lines.ExPot) # excitation potential, i.e. lower level energy (must be in eV)

98-element Vector{Korg.Line{Float64, Float64, Float64, Float64, Float64, Float64}}:
 Fe I 5044.211 Å (log gf = -2.06, χ = 2.85 eV)
 Fe I 5054.642 Å (log gf = -1.92, χ = 3.64 eV)
 Fe I 5127.359 Å (log gf = -3.31, χ = 0.92 eV)
 Fe I 5127.679 Å (log gf = -6.12, χ = 0.05 eV)
 Fe II 5197.577 Å (log gf = -2.22, χ = 3.23 eV)
 Fe I 5198.711 Å (log gf = -2.13, χ = 2.22 eV)
 Fe I 5225.525 Å (log gf = -4.79, χ = 0.11 eV)
 Fe II 5234.625 Å (log gf = -2.18, χ = 3.22 eV)
 Fe I 5242.491 Å (log gf = -0.97, χ = 3.63 eV)
 Fe I 5247.05 Å (log gf = -4.95, χ = 0.09 eV)
 ⋮
 Fe I 6750.152 Å (log gf = -2.62, χ = 2.42 eV)
 Fe I 6752.707 Å (log gf = -1.2, χ = 4.64 eV)
 Fe I 6793.259 Å (log gf = -2.33, χ = 4.08 eV)
 Fe I 6806.845 Å (log gf = -3.11, χ = 2.73 eV)
 Fe I 6810.263 Å (log gf = -0.99, χ = 4.61 eV)
 Fe I 6837.006 Å (log gf = -1.69, χ = 4.59 eV)
 Fe I 6839.83 Å (log gf = -3.35, χ = 2.56 eV)
 Fe I 6843.656 Å (log gf = -0.86, χ = 4.55 eV)
 Fe I 6858.15 Å (log gf = -0.93, χ = 4.61 eV)

We can use lines directly with Korg at this point, but if we want to save it for later use, we can save it to an hdf5 file.

#read this back in with Korg.read_linelist("Fe_lines.h5"; format="korg")
Korg.save_linelist("Fe_lines.h5", linelist);

Equivalent width fitting

Now that we have our linelist, let's perform equivalent width fitting to determine abundances. We'll compare the Sun and 18 Sco using the equivalent widths from the paper.

First, we'll define the parameters for the Sun and 18 Sco.

#solar params
sun_Teff, sun_logg, sun_Fe_H, sun_vmic = 5777, 4.44, 0.0, 1.0
#vector of abundances for the sun
sun_A_X = Korg.format_A_X(sun_Fe_H)
#model atmosphere for the sun
sun_atm = Korg.interpolate_marcs(sun_Teff, sun_logg, sun_A_X)

Korg.PlanarAtmosphere{Float64, Float64, Float64, Float64, Float64} with 56 layers

and likewise for 18 Sco

sco_teff, sco_logg, sco_fe_h, sco_vmic = (5823, 4.45, 0.054, sun_vmic + 0.02)
sco_A_X = Korg.format_A_X(sco_fe_h)
sco_atm = Korg.interpolate_marcs(sco_teff, sco_logg, sco_A_X)

Korg.PlanarAtmosphere{Float64, Float64, Float64, Float64, Float64} with 56 layers

Now we can calculate abundances from the EWs for each star.

A_sun = Korg.Fit.ews_to_abundances(sun_atm, linelist, sun_A_X, lines.EW_Sun; vmic=sun_vmic)
A_18Sco = Korg.Fit.ews_to_abundances(sco_atm, linelist, sco_A_X, lines.EW_18Sco; vmic=sco_vmic);

Let's plot the abundance differences as a function of excitation potential to check for non-LTE effects or other systematic issues.

neutrals = [spec.charge == 0 for spec in lines.species] #bitmask for the lines of Fe I vs Fe II
χ = [l.E_lower for l in linelist] # the excitation potential for each line

scatter(χ[neutrals], (A_18Sco-A_sun)[neutrals]; label="Fe I")
scatter(χ[.!neutrals], (A_18Sco-A_sun)[.!neutrals]; label="Fe II")
ylabel(L"A(Fe)_\mathrm{18 Sco} - A(Fe)_\mathrm{sun}")
xlabel("χ [eV]")
legend()