gaia-benchmarks-hermes.py

# coding: utf-8

""" Analyse Gaia 'Benchmark' stars observed with the HERMES spectrograph. """

__author__ = "Andy Casey <andy@astrowizici.st>"

import os
from glob import glob

import smh

hermes_configuration_filename = os.path.join(
    os.path.dirname(__file__), "smh-hermes.json")

gaia_benchmark_stars = {
    "hd22879":      "HD22879_lowres/*ex3.ex.fits",
    "hd49933":      "HD49933_lowres/*ex4.ex.fits",
    "procyon":      "Procyon_lowres/*ex3.ex.fits",
    "alpha_tau":    "alpha_Tau_lowres/*ex3.ex.fits",
    "beta_hyi":     "beta_Hyi_lowres/*ex.fits",
    "delta_eri":    "delta_Eri/*.ex.fits", 
    "epsilon_eri":  "epsilon_Eri/*.ex.fits",
    "epsilon_for":  "epsilon_For/*.ex.fits",
    "sol": "../sol/uvessun*.txt"
}

# We will use the Sun as a cross-correlation template
sol_blue = smh.specutils.Spectrum1D.load(os.path.join(
    os.path.dirname(__file__), "sol", "uvessun1.txt"))

for i, (star, filename_wildmask) in enumerate(gaia_benchmark_stars.iteritems()):

    spectra = glob(os.path.join(os.path.dirname(__file__),
        "benchmark-stars", filename_wildmask))

    session = smh.Session()
    session.load_spectra(spectra)
    session.load_settings(hermes_configuration_filename)

    # Since our normalisation arguments includes regions to exclude
    # -- and these regions are at rest-frame -- we will measure the
    # radial velocity first, and place the spectrum at rest. In the
    # final GALAH pipeline the normalisation and radial velocity is
    # expected to be handled by the 'GAP' before spectral analysis,
    # so this is a temporary step.

    v_rad, v_err = smh.specutils.cross_correlate(
        observed=session.initial_orders[0], # Blue aperture
        template=sol_blue,
        wl_region=[4750, 4850])

    session.notes += \
        "Measured v_rad = {0:.1f} +/- {1:.1f} km/s\n".format(v_rad, v_err)

    # Shift the initial orders by -v_rad
    session.initial_orders = [spectrum.doppler_correct(-v_rad) \
        for spectrum in session.initial_orders]

    # Normalise all apertures and stitch them
    session.normalisation.normalise_apertures()
    session.normalisation.stitch_apertures()

    # We have already placed these spectra at rest
    session.doppler.set_already_at_rest()

    # Measure equivalent widths
    session.measure_equivalent_widths(session.balance_transitions)

    # Ignore transitions with extreme equivalent width measurements
    for atomic_transition in session.measurements:
        if not (100 >= atomic_transition.measured_equivalent_width >= 5): # milliAngstroms
            atomic_transition.is_acceptable = False

    # Solve for stellar parameters
    session.solve_stellar_parameters(
        max_attempts=10,
        total_tolerance=4e-6,
        individual_tolerances=[1e-3, 1e-3, 1e-3, 1e-3],
        outlier_rejection_total_tolerance=5e-3,
        outlier_limits={26.0: 1.5}
        )

    session.save_as(os.path.join(os.path.dirname(__file__),
        "benchmark-stars", "analysed", "{0}.smh".format(star)),
        clobber=True)