Example 3: Generate mock observed spectrum

About

This example will show you how to manipulate the raw data using functions in breidablik.analysis.read and breidablik.analysis.tools.

Example

We start with our simulated spectra for the sun with 1.1 dex ⁷Li, sun.txt, this can be downloaded from https://zenodo.org/records/10428804

First, we’ll read in the simulated spectra and wavelengths.

# read spectra
import numpy as np
flux, _ = np.loadtxt('sun.txt', unpack = True)

# read wavelengths
from breidablik.analysis import read
wl = read.get_wavelengths()
# all of the simulated spectra have the same wavelength points

Next, we’ll narrow down the spectra to the region we care about, the original spectra ranges from about 600.0 nm to 820.0 nm.

# cut spectra
from breidablik.analysis import tools
wl_cut, flux_cut = tools.cut(wl, flux, center = 670.9659,
                             upper = 0.4, lower = 0.4)

Lastly, to mimic real observations, we’ll need spectra observed at roughly equidistant points in wavelength, and some noise.

# re-interpolate to equidistant wavelength points
from scipy.interpolate import CubicSpline
num_points = 500
wl_equidist = np.linspace(wl_cut[0], wl_cut[-1], num_points)
flux_equidist = CubicSpline(wl_cut, flux_cut)(wl_equidist)

# add Gaussian noise
from numpy.random import normal, seed
seed(19) # set seed for replicability
noise_scale = 0.003 # control how noisy the spectrum is
noise = normal(scale = noise_scale, size = num_points)
flux_noise = flux_equidist + noise
# generate array for error in each pixel
flux_err = np.full(len(flux_noise), noise_scale)

Now we can write this spectra to a file

# write to file
name = 'example_spec.txt'
data = np.array([wl_equidist, flux_noise, flux_err]).T
header = 'wavelength (nm) \t normalised flux \t flux error'
np.savetxt(name, data, fmt = '%.5e', header = header)