Skip to content

Add basic Turbospectrum support #268

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
andycasey wants to merge 14 commits into
base: master
Choose a base branch
from
Open

Add basic Turbospectrum support #268

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
efd3d82
initial
andycasey Sep 29, 2017
ef29a4b
sandboxing
andycasey Sep 29, 2017
b9c640d
add example etc
andycasey Sep 29, 2017
c0f9b8f
move photosphere pickler to scripts/
andycasey Sep 29, 2017
c5437e8
interpolate marcs models while retaining Prad, radius, and write MARC…
andycasey Sep 29, 2017
d69f7ca
hacky but worky
andycasey Sep 29, 2017
1043a51
write line list in turbospectrum format, even if we dont have all the…
andycasey Sep 29, 2017
09b036c
hack
andycasey Sep 29, 2017
2e0151e
add functionality to calculate equivalent widths for transitions usin…
andycasey Sep 29, 2017
8f1e68b
add function to calculate abundances from equivalent widths
andycasey Sep 29, 2017
efb7926
fix indexing problem that TS causes
andycasey Sep 29, 2017
c54cd82
Add defaults for turbospectrum
andycasey Sep 29, 2017
2513d5e
make equivalent_widths and abundances_ return same format as moog.abu…
andycasey Sep 29, 2017
fd6cc86
update zenodo reference to corrected marcs atmospheres
andycasey Sep 29, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 1 addition & 5 deletions smh/photospheres/pickler.py → scripts/photosphere_pickler.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,16 +8,12 @@
__author__ = "Andy Casey <arc@ast.cam.ac.uk>"

import cPickle as pickle
import gzip
import os
import sys
from glob import glob

import numpy as np

import marcs
import castelli_kurucz
import stagger
from smh.photospheres import (marcs, castelli_kurucz)


def pickle_photospheres(photosphere_filenames, kind, meta=None):
Expand Down
6 changes: 3 additions & 3 deletions setup.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,9 +31,9 @@ def read(filename):
# Castelli & Kurucz (2004)
("https://zenodo.org/record/14964/files/castelli-kurucz-2004.pkl",
"smh/photospheres/castelli-kurucz-2004.pkl"),
# MARCS (2008)
("https://zenodo.org/record/14964/files/marcs-2011-standard.pkl",
"smh/photospheres/marcs-2011-standard.pkl"),
# MARCS (2011)
("https://zenodo.org/record/999271/files/marcs-2011_m1.0_t02_st.pkl",
"smh/photospheres/marcs-2011_m1.0_t02_st.pkl"),
# Stagger-Grid <3D> (2013)
("https://zenodo.org/record/15077/files/stagger-2013-optical.pkl",
"smh/photospheres/stagger-2013-optical.pkl"),
Expand Down
38 changes: 38 additions & 0 deletions smh/linelists.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -718,6 +718,42 @@ def write_moog(self,filename):
EW = space
f.write(fmt.format(line['wavelength'],line['species'],line['expot'],line['loggf'],C6,D0,EW,line['comments'])+"\n")



def write_turbospectrum(self, filename):

output = ""

logger.warn("Using default (wrong) physics due to incomplete line list")

# Need to separate by species first.
for species in np.sort(np.unique(self["species"])):

match = self["species"] == species

ion = 1 + int(10 * (species - int(species)))
N = sum(match)

output += "'{0:20.6f}' {1:4.0f} {2:9.0f}\n".format(species, ion, N)
output += "'{0:7s}'\n".format(species_to_element(species))

for transition in self[match]:

row = dict(lower_orbital_type='x', upper_orbital_type='x',
fdamp=-7.750, upper_g=1.0, rad=1e8, equivalent_width_err=0)
row.update(dict(zip(transition.dtype.names, transition.data)))

if not np.isfinite(row["equivalent_width"]):
row["equivalent_width"] = 0

output += "{wavelength:10.3f} {expot:9.5f} {loggf:6.3f} {fdamp:8.3f} {upper_g:6.1f} {rad:9.2e} '{lower_orbital_type:s}' '{upper_orbital_type:s}' {equivalent_width:5.1f} {equivalent_width_err:6.1f} ''\n".format(**row)

with open(filename, "w") as fp:
fp.write(output)

return None


def write_latex(self,filename,sortby=['species','wavelength'],
write_cols = ['wavelength','element','expot','loggf']):
new_table = self.copy()
Expand All @@ -732,3 +768,5 @@ def _moog_identifier(*args, **kwargs):
registry.register_reader("moog", LineList, LineList.read_moog)
registry.register_identifier("moog", LineList, _moog_identifier)

registry.register_writer("turbospectrum", LineList, LineList.write_turbospectrum)

66 changes: 39 additions & 27 deletions smh/photospheres/interpolator.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,41 +102,56 @@ def __call__(self, *args, **kwargs):
return self.interpolate(*args, **kwargs)


def _return_photosphere(self, stellar_parameters, quantities):
"""
Prepare the interpolated photospheric quantities (with correct columns,
units, metadata, etc).
"""
@property
def stellar_parameters_grid(self):
return self.stellar_parameters.view(float).reshape(
len(self.stellar_parameters), -1)

meta = self.meta.copy()
meta["common_optical_scale"] = self.opacity_scale
meta["stellar_parameters"] = \

def _return_photosphere(self, stellar_parameters, quantities, meta=None):

_meta = self.meta.copy()
_meta["common_optical_scale"] = self.opacity_scale
_meta["stellar_parameters"] = \
dict(zip(self.stellar_parameters.dtype.names, stellar_parameters))
units = meta.pop("photospheric_units", None)

photosphere = Photosphere(data=quantities, meta=meta,
_meta.update(meta or {})
units = _meta.pop("photospheric_units", None)

photosphere = Photosphere(data=quantities, meta=_meta,
names=self.photospheric_quantities)
if units is not None:
for name, unit in zip(self.photospheric_quantities, units):
photosphere[name].unit = unit
return photosphere


def _prepare_photosphere(self, stellar_parameters, quantities,
neighbour_indices):
"""
Prepare the interpolated photospheric quantities (with correct columns,
units, metadata, etc).
"""

return self._return_photosphere(stellar_parameters, quantities)



def nearest_neighbours(self, point, n):
"""
Return the indices of the n nearest neighbours to the point.
"""

stellar_parameters = _recarray_to_array(self.stellar_parameters)
distances = np.sum(((point - stellar_parameters) \
/ np.ptp(stellar_parameters, axis=0))**2, axis=1)
spg = self.stellar_parameters_grid
distances = np.nansum(((point - spg) / np.ptp(spg, axis=0))**2, axis=1)
return distances.argsort()[:n]


def nearest(self, *point):
logger.warn("Living dangerously!")
return self._return_photosphere(point,
self.photospheres[self.nearest_neighbours(point, 1)[0]])
neighbour = self.nearest_neighbours(point, 1)[0]
return self._prepare_photosphere(point, self.photospheres[neighbour],
neighbour)


def _test_interpolator(self, *point):
Expand All @@ -161,34 +176,33 @@ def interpolate(self, *point, **kwargs):

__ignore_nearest = kwargs.pop("__ignore_nearest", False)

grid = self.stellar_parameters.view(float).reshape(
len(self.stellar_parameters), -1)
grid_index = np.all(grid == point, axis=1)
stellar_parameters_grid = self.stellar_parameters_grid
grid_index = np.all(stellar_parameters_grid == point, axis=1)
if np.any(grid_index) and not __ignore_nearest:
grid_index = np.where(grid_index)[0][0]
return self._return_photosphere(point, self.photospheres[grid_index])
return self._prepare_photosphere(point, self.photospheres[grid_index],
grid_index)

# Work out what the optical depth points will be in our (to-be)-
# interpolated photosphere.
if __ignore_nearest:
neighbours = self.nearest_neighbours(point, self.neighbours + 1)[1:]
else:
neighbours = self.nearest_neighbours(point, self.neighbours)
stellar_parameters = _recarray_to_array(self.stellar_parameters)

# Shapes required for griddata:
# points: (N, ndim)
# values: shape (N, )
# xi: shape (M, ndim)

# Protect Qhull from columns with a single value.
cols = _protect_qhull(stellar_parameters[neighbours])
cols = _protect_qhull(stellar_parameters_grid[neighbours])
shape = [self.neighbours] + list(self.photospheres.shape[1:])
if self.opacity_scale is not None:
opacity_index = self.photospheric_quantities.index(self.opacity_scale)
kwds = {
"xi": point[cols].reshape(1, len(cols)),
"points": stellar_parameters[neighbours][:, cols],
"points": stellar_parameters_grid[neighbours][:, cols],
"values": self.photospheres[neighbours, :, opacity_index],
"method": self.method,
"rescale": self.rescale
Expand Down Expand Up @@ -233,7 +247,7 @@ def interpolate(self, *point, **kwargs):
# Now interpolate the photospheric quantities.
kwds = {
"xi": point[cols].reshape(1, len(cols)),
"points": stellar_parameters[neighbours][:, cols],
"points": stellar_parameters_grid[neighbours][:, cols],
"values": neighbour_quantities,
"method": self.method,
"rescale": self.rescale
Expand All @@ -257,7 +271,8 @@ def interpolate(self, *point, **kwargs):
interpolated_quantities[:, index] = \
10**interpolated_quantities[:, index]

return self._return_photosphere(point, interpolated_quantities)
return self._prepare_photosphere(point, interpolated_quantities,
neighbours)


def resample_photosphere(opacities, photosphere, opacity_index):
Expand All @@ -282,9 +297,6 @@ def resample_photosphere(opacities, photosphere, opacity_index):
return resampled_photosphere


def _recarray_to_array(a, dtype=float):
return a.view(dtype).reshape(len(a), -1)

def _protect_qhull(a):
return np.where([np.unique(a[:, i]).size > 1 for i in range(a.shape[1])])[0]

Expand Down
51 changes: 43 additions & 8 deletions smh/photospheres/marcs.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,18 +35,28 @@ def __init__(self, **kwargs):
of 1 km/s in plane-parallel models and 2 km/s in spherical models.

"""
return super(self.__class__, self).__init__("marcs-2011-standard.pkl",
return super(self.__class__, self).__init__("marcs-2011_m1.0_t02_st.pkl",
**kwargs)


@property
def stellar_parameters_grid(self):
"""
Overload this function so that the radii are considered 0 or 1 for the
purposes of interpolation.
"""

grid = self.stellar_parameters.view(float).reshape(
len(self.stellar_parameters), -1).copy()
grid[:, -1] = np.clip(grid[:, -1], 0, 1)
return grid


def _spherical_or_plane_parallel(self, *point):

point = list(point) + [0.5] # equi-spaced from plane-parallel/spherical
neighbours = self.nearest_neighbours(point, 8) # 8 = 2**3

sph_or_pp = self.stellar_parameters.view(float).reshape(
len(self.stellar_parameters), -1)[:, -1]
return np.round(np.median(sph_or_pp[neighbours]))
return np.round(np.median(self.stellar_parameters_grid[neighbours, -1]))


def interpolate(self, *point, **kwargs):
Expand Down Expand Up @@ -75,6 +85,19 @@ def interpolate(self, *point, **kwargs):



def _prepare_photosphere(self, stellar_parameters, quantities,
neighbour_indices):
"""
Prepare the interpolated photospheric quantities (with correct columns,
units, metadata, etc).
"""

radius = np.mean(self.stellar_parameters["radius"][neighbour_indices])
return self._return_photosphere(
stellar_parameters, quantities, meta=dict(radius=radius))




def parse_filename(filename, full_output=False):
"""
Expand All @@ -85,17 +108,29 @@ def parse_filename(filename, full_output=False):
teff = basename[1:5]
logg = basename.split("_")[1][1:]
feh = basename.split("_")[5][1:]
parameters = map(float, [teff, logg, feh, int(basename[0].lower() == "s")])

is_spherical = int(basename[0].lower() == "s")
if is_spherical:
# Need to open the file and get the radius from line 8
f = gzip.open if filename.lower().endswith(".gz") else open
with f(filename, "r") as fp:
content = fp.readlines()
radius = content[7].split()[0]

else:
radius = 0

parameters = np.array([teff, logg, feh, radius]).astype(float)

if full_output:
names = ("effective_temperature", "surface_gravity", "metallicity",
"is_spherical?")
"radius")
return (parameters, names)
return parameters


def parse_photospheric_structure(filename, ndepth=56, line=25,
columns=("lgTau5", "Depth", "T", "Pe", "Pg"), full_output=False):
columns=("lgTau5", "Depth", "T", "Pe", "Pg", "Prad"), full_output=False):
"""
Parse the photospheric structure (optical depths, temperatures, electron
and gas pressures) from the filename provided.
Expand Down
64 changes: 54 additions & 10 deletions smh/photospheres/photosphere.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,8 +6,6 @@

from __future__ import division, absolute_import, print_function

__author__ = "Andy Casey <arc@ast.cam.ac.uk>"

import logging
from textwrap import dedent

Expand All @@ -18,13 +16,56 @@
logger = logging.getLogger(__name__)

class Photosphere(astropy.table.Table):
""" A model photosphere object. """
pass


def _turbospectrum_writer(photosphere, filename, **kwargs):
"""
A model photosphere object.
Writes a :class:`photospheres.Photosphere` to file in a format that can be
read by Turbospectrum.

:param photosphere:
The photosphere.

:param filename:
The filename to write the photosphere to.
"""
pass

if photosphere.meta["kind"] != "marcs":
raise ValueError("only marcs photospheres supported with turbospectrum")

output = ""
radius = photosphere.meta["radius"]
if radius > 0:
output += \
"spherical model\n"\
" 1.0 Mass [Msun]\n"\
" {:.4e} Radius [cm] At Tau(Rosseland)=1.0\n".format(radius)

else:
output += \
"plane-parallel model\n"\
" 0.0 No mass for plane-parallel models\n"\
" 1.0000E+00 1 cm radius for plane-parallel models\n"

output += " {:.0f} Number of depth points\n".format(len(photosphere))
output += "Model structure\n"
output += " k lgTauR lgTau5 Depth T Pe Pg Prad Pturb\n"

lgTauR = -5
for i, layer in enumerate(photosphere):
lgTauR = -5 + i * 0.20
output += "{:3.0f} {:5.2f} {:7.4f} {:10.3e} {:7.1f} {:10.3e} {:10.3e} {:10.3e} {:10.3e}\n"\
.format(i + 1, lgTauR, layer["lgTau5"], layer["Depth"], layer["T"],
layer["Pe"], layer["Pg"], layer["Prad"], 0)

with open(filename, "w") as fp:
fp.write(output)

return None


# MOOG writer and identifier.
def _moog_writer(photosphere, filename, **kwargs):
"""
Writes an :class:`photospheres.photosphere` to file in a MOOG-friendly
Expand Down Expand Up @@ -105,9 +146,12 @@ def _get_xi():
return None


def _moog_identifier(*args, **kwargs):
return isinstance(args[0], basestring) and args[0].lower().endswith(".moog")

# Register the MOOG writer and identifier
# Register writers.
astropy.io.registry.register_writer("moog", Photosphere, _moog_writer)
astropy.io.registry.register_identifier("moog", Photosphere, _moog_identifier)
astropy.io.registry.register_writer("turbospectrum", Photosphere, _turbospectrum_writer)

# Register identifiers.
astropy.io.registry.register_identifier("moog", Photosphere,
lambda *a, **k: "{0}".format(a[0]).lower().endswith(".moog"))
astropy.io.registry.register_identifier("turbospectrum", Photosphere,
lambda *a, **k: "{0}".format(a[0]).lower().endswith(".turbospectrum"))
Loading