Wavelength oddities

[andycasey]

I have encountered some strange behaviour with Korg v0.4.1. Maybe it is a bug in how Julia is handling floating points, or it might not be a bug: it might just be how things should be. This is low priority as I have found a workaround. A minimum reproducible example is below.

Let me describe the example. I want to synthesize between 15,000 and 15,500 Angstroms (air wavelengths), at 0.01 Angstrom spacings.

Case A Since Korg works in vacuum wavelengths I first convert 15,000 and 15,500 Angstroms to vacuum wavelengths and then supply those bounds to Korg, with 0.01 Angstrom spacings: 15004.099625020779:0.01:15504.245608304296.

Korg returns a spectrum with length 50,014 but I expected 50,001 (+/- 1, depending on how the array is constructed). For now I will just trim the edge of the spectrum to match the number of points.

Case B Now for the comparison. I will now do the synthesis where I supply fewer significant digits to Korg: 15004.09:0.01:15504.24. This will give me a few more flux calculations than what I need, but I can use it to interpolate or trim it down, etc.

Below you can see the plotted spectra for Case A (red, bad) and Case B (blue, good). I say that Case B is "good" because when I compare to other radiative transfer codes, Case B looks correct but Case A looks incorrect. You can see that both spectra agree in the blue end of the spectra, but by the right end of the spectra there is significant shifting going on.

Is this something to do with how Julia constructs the wavelength array differently to how Python would? Is there some rounding bug going on?

To compare with other radiative transfer codes I have switched to using Case B and interpolating the Korg spectrum to the comparison spectra.

Minimum reproducible example

Run example.jl then example.py for plotting. Files for the minimum reproducible example:

• atmosphere.txt
• transitions.txt

example.jl code:

using Pkg
using Korg

Pkg.status("Korg")

atm = Korg.read_model_atmosphere("atmosphere.txt")
linelist = Korg.read_linelist("transitions.txt", format="vald")

bad_spectrum = Korg.synthesize(atm, linelist, 15004.099625020779:0.01:15504.245608304296; metallicity=0)
bad_continuum = Korg.synthesize(atm, [], 15004.099625020779:0.01:15504.245608304296; metallicity=0)

open("bad_spectrum.txt", "w") do fp
    for (flux, cont) in zip(bad_spectrum.flux, bad_continuum.flux)
        println(fp, flux / cont)
    end
end

spectrum = Korg.synthesize(atm, linelist, 15004.09:0.01:15504.24; metallicity=0)
continuum = Korg.synthesize(atm, [], 15004.09:0.01:15504.24; metallicity=0)

open("spectrum.txt", "w") do fp
    for (flux, cont) in zip(spectrum.flux, continuum.flux)
        println(fp, flux / cont)
    end
end

example.py code:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator   

def air_to_vacuum(lambdas):
    s = 10**4 / lambdas
    n = 1 + 0.00008336624212083 + 0.02408926869968 / (130.1065924522 - s**2) + 0.0001599740894897 / (38.92568793293 - s**2)
    return lambdas * n

def vacuum_to_air(lambdas):
    s = 10**4 / lambdas
    n = 1 + 0.0000834254 + 0.02406147 / (130 - s**2) + 0.00015998 / (38.9 - s**2)
    return (lambdas / n) 

air_wavelength = np.arange(15000, 15500 + 0.01, 0.01)

bad_wavelength = air_to_vacuum(air_wavelength)
assert bad_wavelength[0] == 15004.099625020779
assert bad_wavelength[-1] == 15504.245608304296

bad_spectrum = np.loadtxt("bad_spectrum.txt")

print(f"bad_wavelength has shape {bad_wavelength.shape} and bad_spectrum has shape {bad_spectrum.shape}")

bad_wavelength = bad_wavelength[:len(bad_spectrum)]

wavelength = np.arange(15004.09, 15504.24 + 0.01, 0.01)
spectrum = np.loadtxt("spectrum.txt")

print(f"wavelength has shape {wavelength.shape} and spectrum has shape {spectrum.shape}")

fig, axes = plt.subplots(1, 2, figsize=(12, 6))

# Since bad_wavelength is shape 50002 and bad_spectrum is shape 50014, we need to trim the arrays to the same length.
for ax in axes:
    ax.plot(vacuum_to_air(bad_wavelength), bad_spectrum[:len(bad_wavelength)], label="Bad spectrum", c="tab:red")
    ax.plot(vacuum_to_air(wavelength), spectrum, label="Good spectrum", c="tab:blue", ls="-.")

    ax.xaxis.set_major_locator(MaxNLocator(6, integer=True))
    ax.ticklabel_format(useOffset=False)

    ax.set_xlabel("Wavelength (air) [Angstroms]")

axes[0].legend(loc="lower left", frameon=False)
axes[-1].set_yticklabels([])
axes[0].set_xlim(15000 - 1, 15010)
axes[-1].set_xlim(15490, 15500 + 1)
fig.tight_layout()

fig.savefig("example.png", dpi=300)

[ajwheeler]

I'm fairly sure this is a direct consequence of the nonlinear relationship between air and vacuum wavelengths, i.e. working as intended. The range object 15004.099625020779 : 0.01 : 15504.245608304296 itself has 50015 elements just by the nature of its bounds and step size.

For very small wavelength ranges, converting starting and stopping wavelengths and keeping the step the same is a reasonable approximation. For large ones it will fail, so for the best possible comparisons between synthesized air and vacuum spectra, interpolation is required at some level. However, for the 15000-15500 \AA range, we are actually in a middle regime where the relationship between air and vacuum wavelengths is linear-ish, as long you adapt both the step size and the range endpoints.

I've opened a PR, which I will try to merge soon, that will do this conversion automatically for you, raising a warning if the approximation becomes inaccurate. This overlaps with some other things I've been thinking about, and will make it easier to change the internals later, e.g. to use frequencies instead of wavelengths.

[andycasey]

That's great! I thought that could be the case, but I could not find a .wavelength property (or similar) in the outputted spectrum object to check for sure. Sorry for the noise! I'll leave this open until it gets closed by your pull request.