From efd3d8254ec209178149d47c528770c3bc551a82 Mon Sep 17 00:00:00 2001 From: Andy Casey Date: Fri, 29 Sep 2017 12:48:16 +1000 Subject: [PATCH 01/14] initial --- smh/radiative_transfer/__init__.py | 3 ++- smh/radiative_transfer/turbospectrum/__init__.py | 12 ++++++++++++ 2 files changed, 14 insertions(+), 1 deletion(-) create mode 100644 smh/radiative_transfer/turbospectrum/__init__.py diff --git a/smh/radiative_transfer/__init__.py b/smh/radiative_transfer/__init__.py index 86ac8f43..81344579 100644 --- a/smh/radiative_transfer/__init__.py +++ b/smh/radiative_transfer/__init__.py @@ -20,4 +20,5 @@ def abundance_cog(photosphere, transitions, **kwargs) """ -from .moog import * \ No newline at end of file +from .moog import * +from .turbospectrum import * diff --git a/smh/radiative_transfer/turbospectrum/__init__.py b/smh/radiative_transfer/turbospectrum/__init__.py new file mode 100644 index 00000000..7d44cfb5 --- /dev/null +++ b/smh/radiative_transfer/turbospectrum/__init__.py @@ -0,0 +1,12 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +""" A standardized interface to the MOOG(SILENT) radiative transfer code. """ + +from __future__ import (division, print_function, absolute_import, + unicode_literals) + +#__all__ = ["abundance_cog", "synthesize", "RTError"] + +# See stackoverflow.com/questions/19913653/no-unicode-in-all-for-a-packages-init +#__all__ = [_.encode("ascii") for _ in __all__] From ef29a4bfdfb8ba33da810f4e7607eae71b28c75d Mon Sep 17 00:00:00 2001 From: Andy Casey Date: Fri, 29 Sep 2017 13:06:28 +1000 Subject: [PATCH 02/14] sandboxing --- sandbox_turbospectrum.py | 7 ++++ .../turbospectrum/synthesis.py | 41 +++++++++++++++++++ 2 files changed, 48 insertions(+) create mode 100644 sandbox_turbospectrum.py create mode 100644 smh/radiative_transfer/turbospectrum/synthesis.py diff --git a/sandbox_turbospectrum.py b/sandbox_turbospectrum.py new file mode 100644 index 00000000..85932821 --- /dev/null +++ b/sandbox_turbospectrum.py @@ -0,0 +1,7 @@ + +from smh.radiative_transfer import turbospectrum as ts +from smh.photospheres.marcs import Interpolator + +marcs = Interpolator() + +solar_photosphere = marcs(5777, 4.4, 0) diff --git a/smh/radiative_transfer/turbospectrum/synthesis.py b/smh/radiative_transfer/turbospectrum/synthesis.py new file mode 100644 index 00000000..cba95530 --- /dev/null +++ b/smh/radiative_transfer/turbospectrum/synthesis.py @@ -0,0 +1,41 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +""" Synthesis functionality using the Turbospectrum radiative transfer code. """ + +from __future__ import (division, print_function, absolute_import, + unicode_literals) + +import logging +import numpy as np +import yaml +from pkg_resources import resource_stream + +logger = logging.getLogger(__name__) + +# Load defaults. +with resource_stream(__name__, "defaults.yaml") as fp: + _turbospectrum_defaults = yaml.load(fp) + + + +def synthesize(photosphere, transitions, abundances=None, isotopes=None, + verbose=False, twd=None, **kwargs): + """ + Sythesize a stellar spectrum given the model photosphere and transitions + provided. + + """ + + raise a + # Get regions to synthesise. + + + # Update abundances + + # Calculate opacities + + # Synthesise + + + From b9c640dd0f5d1743df93afe0a90e312694f00d1e Mon Sep 17 00:00:00 2001 From: Andy Casey Date: Fri, 29 Sep 2017 13:48:15 +1000 Subject: [PATCH 03/14] add example etc --- sandbox_turbospectrum.py | 11 ++- .../turbospectrum/__init__.py | 2 + .../turbospectrum/babsma_lu.in | 15 ++++ .../turbospectrum/synthesis.py | 68 +++++++++++++++++++ smh/radiative_transfer/turbospectrum/utils.py | 24 +++++++ 5 files changed, 119 insertions(+), 1 deletion(-) create mode 100644 smh/radiative_transfer/turbospectrum/babsma_lu.in create mode 100644 smh/radiative_transfer/turbospectrum/utils.py diff --git a/sandbox_turbospectrum.py b/sandbox_turbospectrum.py index 85932821..673c2de9 100644 --- a/sandbox_turbospectrum.py +++ b/sandbox_turbospectrum.py @@ -1,7 +1,16 @@ from smh.radiative_transfer import turbospectrum as ts from smh.photospheres.marcs import Interpolator +from smh.linelists import LineList + + -marcs = Interpolator() +marcs = Interpolator() solar_photosphere = marcs(5777, 4.4, 0) +linelist = LineList.read("atomic.list") + +# Just restrict to a few lines. +linelist = linelist[500:505] + +ts.synthesize(solar_photosphere, linelist, None) diff --git a/smh/radiative_transfer/turbospectrum/__init__.py b/smh/radiative_transfer/turbospectrum/__init__.py index 7d44cfb5..fb3ea98e 100644 --- a/smh/radiative_transfer/turbospectrum/__init__.py +++ b/smh/radiative_transfer/turbospectrum/__init__.py @@ -10,3 +10,5 @@ # See stackoverflow.com/questions/19913653/no-unicode-in-all-for-a-packages-init #__all__ = [_.encode("ascii") for _ in __all__] + +from .synthesis import synthesize \ No newline at end of file diff --git a/smh/radiative_transfer/turbospectrum/babsma_lu.in b/smh/radiative_transfer/turbospectrum/babsma_lu.in new file mode 100644 index 00000000..7401e308 --- /dev/null +++ b/smh/radiative_transfer/turbospectrum/babsma_lu.in @@ -0,0 +1,15 @@ +'LAMBDA_MIN:' '{dispersion_min:.2f}' +'LAMBDA_MAX:' '{dispersion_max:.2f}' +'LAMBDA_STEP:' '{dispersion_delta:.2f}' +'MODELINPUT:' '{photosphere_path}' +'MARCS-FILE:' '.true.' +'MODELOPAC:' '{opacities_path}' +'METALLICITY:' '0.00' +'ALPHA/Fe :' '0.00' +'HELIUM :' '0.00' +'R-PROCESS :' '0.00' +'S-PROCESS :' '0.00' +'INDIVIDUAL ABUNDANCES:' '{num_abundances}' +{formatted_abundances} +'XIFIX:' 'T' +{microturbulence:.2f} diff --git a/smh/radiative_transfer/turbospectrum/synthesis.py b/smh/radiative_transfer/turbospectrum/synthesis.py index cba95530..b60590ba 100644 --- a/smh/radiative_transfer/turbospectrum/synthesis.py +++ b/smh/radiative_transfer/turbospectrum/synthesis.py @@ -11,6 +11,8 @@ import yaml from pkg_resources import resource_stream +from . import utils + logger = logging.getLogger(__name__) # Load defaults. @@ -19,18 +21,84 @@ + + def synthesize(photosphere, transitions, abundances=None, isotopes=None, verbose=False, twd=None, **kwargs): """ Sythesize a stellar spectrum given the model photosphere and transitions provided. + Optional keyword arguments include: + dispersion_min + dispersion_max + dispersion_delta + opacity_contribution + """ + + kwds = _turbospectrum_defaults.copy() + kwds.update(kwargs) + + # Set default regions to synthesize. + kwds.setdefault("dispersion_min", + min(transitions["wavelength"]) - kwds["opacity_contribution"]) + kwds.setdefault("dispersion_max", + max(transitions["wavelength"]) + kwds["opacity_contribution"] \ + + kwds["dispersion_delta"]) + + + # TODO put in abundances. + if kwds["dispersion_max"] - kwds["dispersion_min"] > 1000.0: + raise NotImplementedError("turbospectrum fails for large synth, " + "and chunking not implemented yet") + + path = utils.twd_path(twd=twd, **kwargs) + + + # Calculate opacities. + babsma_lu_kwds = kwds.copy() + babsma_lu_kwds["opacities_path"] = path("opacities.in") + babsma_lu_kwds["photosphere_path"] = path("photosphere.in") + + # Write atmosphere file. + photosphere.write(babsma_lu_kwds["photosphere_path"], format="turbospectrum") + + # Get babsma_lu template and format it. + with resource_stream(__name__, "babsma_lu.in") as fp: + babsma_lu_template = fp.read() + + # XI + # num abundances + # formatted abundances. + + # If no XI value, then estimate it. + xi = photosphere.meta.get("microturbulence", None) + if xi is None: + xi = 1.0 # [km/s] + logger.warn("No microturbulence given in photosphere. Assuming 1 km/s.") + + # TODO microturbulence may have been specified by kwargs.... + if "microturbulence" in babsma_lu_kwds: + raise a + + babsma_lu_kwds["microturbulence"] = xi + + # Set some default values??? + + babsma_lu_contents = babsma_lu_template.format(**babsma_lu_kwds) + + + + + + raise a # Get regions to synthesise. + # Update abundances # Calculate opacities diff --git a/smh/radiative_transfer/turbospectrum/utils.py b/smh/radiative_transfer/turbospectrum/utils.py new file mode 100644 index 00000000..609c6892 --- /dev/null +++ b/smh/radiative_transfer/turbospectrum/utils.py @@ -0,0 +1,24 @@ + + +import logging +import os +from smh.utils import mkdtemp + +logger = logging.getLogger(__name__) + +def twd_path(twd=None,**kwargs): + """ + Create a temporary working directory and return a function that will format + basenames from that temporary working directory. + """ + + if twd is None: + kwds = {} + kwds["dir"] = kwargs.get("dir", "/tmp/") + kwds["prefix"] = kwargs.get("prefix", "smh-") + kwds["suffix"] = kwargs.get("suffix", "") + #kwds = kwargs.copy() + #kwds.setdefault("dir", "/tmp/") + #kwds.setdefault("prefix", "smh-") + twd = mkdtemp(**kwds) + return lambda filename: os.path.join(twd, filename) \ No newline at end of file From c0f9b8faa42bcabad7cc627b72866f5e8f12f112 Mon Sep 17 00:00:00 2001 From: Andy Casey Date: Fri, 29 Sep 2017 14:51:35 +1000 Subject: [PATCH 04/14] move photosphere pickler to scripts/ --- .../pickler.py => scripts/photosphere_pickler.py | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) rename smh/photospheres/pickler.py => scripts/photosphere_pickler.py (97%) diff --git a/smh/photospheres/pickler.py b/scripts/photosphere_pickler.py similarity index 97% rename from smh/photospheres/pickler.py rename to scripts/photosphere_pickler.py index 965cbbd5..a635c6f6 100644 --- a/smh/photospheres/pickler.py +++ b/scripts/photosphere_pickler.py @@ -8,16 +8,12 @@ __author__ = "Andy Casey " 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): From c5437e8f9c830f7565b0bda2f9536f1d508c410c Mon Sep 17 00:00:00 2001 From: Andy Casey Date: Fri, 29 Sep 2017 16:26:56 +1000 Subject: [PATCH 05/14] interpolate marcs models while retaining Prad, radius, and write MARCS photospheres to turbospectrum format --- smh/photospheres/interpolator.py | 66 +++++++++++++++++++------------- smh/photospheres/marcs.py | 51 ++++++++++++++++++++---- smh/photospheres/photosphere.py | 64 ++++++++++++++++++++++++++----- 3 files changed, 136 insertions(+), 45 deletions(-) diff --git a/smh/photospheres/interpolator.py b/smh/photospheres/interpolator.py index 19ba71ea..6f1bafa8 100644 --- a/smh/photospheres/interpolator.py +++ b/smh/photospheres/interpolator.py @@ -102,19 +102,23 @@ 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): @@ -122,21 +126,32 @@ def _return_photosphere(self, stellar_parameters, quantities): 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): @@ -161,12 +176,12 @@ 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. @@ -174,7 +189,6 @@ def interpolate(self, *point, **kwargs): 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) @@ -182,13 +196,13 @@ def interpolate(self, *point, **kwargs): # 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 @@ -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 @@ -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): @@ -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] diff --git a/smh/photospheres/marcs.py b/smh/photospheres/marcs.py index e3fe20e2..d8c1f8e0 100644 --- a/smh/photospheres/marcs.py +++ b/smh/photospheres/marcs.py @@ -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): @@ -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): """ @@ -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. diff --git a/smh/photospheres/photosphere.py b/smh/photospheres/photosphere.py index 7bc466b9..9feb13a0 100644 --- a/smh/photospheres/photosphere.py +++ b/smh/photospheres/photosphere.py @@ -6,8 +6,6 @@ from __future__ import division, absolute_import, print_function -__author__ = "Andy Casey " - import logging from textwrap import dedent @@ -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 @@ -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")) From d69f7cac7355cd0a712a9d65ded711b71b65a479 Mon Sep 17 00:00:00 2001 From: Andy Casey Date: Fri, 29 Sep 2017 16:53:42 +1000 Subject: [PATCH 06/14] hacky but worky --- .../turbospectrum/bsyn_lu.in | 24 +++++++ .../turbospectrum/synthesis.py | 65 +++++++++++++++++++ 2 files changed, 89 insertions(+) create mode 100644 smh/radiative_transfer/turbospectrum/bsyn_lu.in diff --git a/smh/radiative_transfer/turbospectrum/bsyn_lu.in b/smh/radiative_transfer/turbospectrum/bsyn_lu.in new file mode 100644 index 00000000..bea70714 --- /dev/null +++ b/smh/radiative_transfer/turbospectrum/bsyn_lu.in @@ -0,0 +1,24 @@ +'LAMBDA_MIN:' '{dispersion_min:.2f}' +'LAMBDA_MAX:' '{dispersion_max:.2f}' +'LAMBDA_STEP:' '{dispersion_delta:.2f}' +'INTENSITY/FLUX:' 'Flux' +'COS(THETA):' '1.00' +'ABFIND :' '.false.' +'MODELOPAC:' '{opacities_path}' +'RESULTFILE': '{spectrum_path}' +'METALLICITY:' '0.00' +'ALPHA/Fe :' '0.00' +'HELIUM :' '0.00' +'R-PROCESS :' '0.00' +'S-PROCESS :' '0.00' +'INDIVIDUAL ABUNDANCES:' '{num_abundances}' +{formatted_abundances} +'ISOTOPES:' '{num_isotopes}' +{formatted_isotopes} +'NFILES:' '{n_files}' +{formatted_n_files} +'SPHERICAL:' '{is_spherical}' + 30 + 300.00 + 15 + 1.30 diff --git a/smh/radiative_transfer/turbospectrum/synthesis.py b/smh/radiative_transfer/turbospectrum/synthesis.py index b60590ba..cb07d49e 100644 --- a/smh/radiative_transfer/turbospectrum/synthesis.py +++ b/smh/radiative_transfer/turbospectrum/synthesis.py @@ -8,6 +8,8 @@ import logging import numpy as np +import os +import subprocess import yaml from pkg_resources import resource_stream @@ -69,6 +71,20 @@ def synthesize(photosphere, transitions, abundances=None, isotopes=None, with resource_stream(__name__, "babsma_lu.in") as fp: babsma_lu_template = fp.read() + if abundances is None: + abundances = {} + + # TODO Put the abundances on the right scale???? + + formatted_abundances = "\n".join(["{0} {1:.2f}".format(species, abundance) \ + for species, abundance in abundances.items()]) + + babsma_lu_kwds.update( + num_abundances=len(abundances), + formatted_abundances=formatted_abundances) + + + # XI # num abundances # formatted abundances. @@ -90,6 +106,55 @@ def synthesize(photosphere, transitions, abundances=None, isotopes=None, babsma_lu_contents = babsma_lu_template.format(**babsma_lu_kwds) + # Calcualte opacities. + command = kwds["turbospectrum_babsma_lu"] + + os.symlink(kwds["turbospectrum_data_dir"], path("DATA")) + proc = subprocess.Popen(command.split(), stdin=subprocess.PIPE, + stdout=subprocess.PIPE, stderr=subprocess.PIPE, + cwd=path("")) + + out, err = proc.communicate(input=babsma_lu_contents) + errcode = proc.returncode + + + # Calculate the spectrum. + os.symlink(kwds["turbospectrum_molecules_dir"], path("molecules")) + command = kwds["turbospectrum_bsyn_lu"] + + with resource_stream(__name__, "bsyn_lu.in") as fp: + bsyn_lu_template = fp.read() + + if isotopes is None: + isotopes = {} + formatted_isotopes = "\n".join(["{0} {1:.2f}".format(mass_number, relative_abundance) \ + for mass_number, relative_abundance in isotopes.items()]) + + + babsma_lu_kwds.update(spectrum_path=path("spectrum.out"), + num_isotopes=len(isotopes), formatted_isotopes=formatted_isotopes, + is_spherical=["F", "T"][photosphere.meta["radius"] > 0]) + + transitions_basename = "transitions.in" + transitions.write(path(transitions_basename), format="turbospectrum") + + # TODO: MOlecule files. + + n_files = ["DATA/Hlinedata", transitions_basename] + + babsma_lu_kwds.update(formatted_n_files="\n".join(n_files), + n_files=len(n_files)) + + + bsyn_lu_contents = bsyn_lu_template.format(**babsma_lu_kwds) + command = kwds["turbospectrum_bsyn_lu"].split() + + + proc_synth = subprocess.Popen(command, stdin=subprocess.PIPE, + stdout=subprocess.PIPE, stderr=subprocess.PIPE) + + out_synth, err_synth = proc_synth.communicate(input=bsyn_lu_contents) + errcode_synth = proc_synth.returncode From 1043a51e5baa26d5a119dfc78e777307fd5f2320 Mon Sep 17 00:00:00 2001 From: Andy Casey Date: Fri, 29 Sep 2017 18:04:27 +1000 Subject: [PATCH 07/14] write line list in turbospectrum format, even if we dont have all the transition properties --- smh/linelists.py | 38 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/smh/linelists.py b/smh/linelists.py index 70caae31..c118ce1e 100644 --- a/smh/linelists.py +++ b/smh/linelists.py @@ -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() @@ -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) + From 09b036cd78f7e2133c858d6f526e52726522dd99 Mon Sep 17 00:00:00 2001 From: Andy Casey Date: Fri, 29 Sep 2017 18:05:18 +1000 Subject: [PATCH 08/14] hack --- ni-test.list | 5 +++++ sandbox_turbospectrum.py | 6 +++--- smh/radiative_transfer/turbospectrum/bsyn_lu.in | 2 +- smh/radiative_transfer/turbospectrum/synthesis.py | 7 +++++-- 4 files changed, 14 insertions(+), 6 deletions(-) create mode 100644 ni-test.list diff --git a/ni-test.list b/ni-test.list new file mode 100644 index 00000000..d57f352e --- /dev/null +++ b/ni-test.list @@ -0,0 +1,5 @@ +5155.126 28.0 3.89800 -0.650 +5155.764 28.0 3.89800 0.074 +5157.707 28.0 3.70600 -1.264 +5157.980 28.0 3.60600 -1.510 +5162.916 28.0 3.83300 -3.200 \ No newline at end of file diff --git a/sandbox_turbospectrum.py b/sandbox_turbospectrum.py index 673c2de9..2ba710f1 100644 --- a/sandbox_turbospectrum.py +++ b/sandbox_turbospectrum.py @@ -8,9 +8,9 @@ marcs = Interpolator() solar_photosphere = marcs(5777, 4.4, 0) -linelist = LineList.read("atomic.list") +linelist = LineList.read("ni-test.list") # Just restrict to a few lines. -linelist = linelist[500:505] +#linelist = linelist[500:505] -ts.synthesize(solar_photosphere, linelist, None) +ts.synthesize(solar_photosphere, linelist, None, dispersion_min=5150, dispersion_max=5165) diff --git a/smh/radiative_transfer/turbospectrum/bsyn_lu.in b/smh/radiative_transfer/turbospectrum/bsyn_lu.in index bea70714..dbd489bb 100644 --- a/smh/radiative_transfer/turbospectrum/bsyn_lu.in +++ b/smh/radiative_transfer/turbospectrum/bsyn_lu.in @@ -5,7 +5,7 @@ 'COS(THETA):' '1.00' 'ABFIND :' '.false.' 'MODELOPAC:' '{opacities_path}' -'RESULTFILE': '{spectrum_path}' +'RESULTFILE:' '{spectrum_path}' 'METALLICITY:' '0.00' 'ALPHA/Fe :' '0.00' 'HELIUM :' '0.00' diff --git a/smh/radiative_transfer/turbospectrum/synthesis.py b/smh/radiative_transfer/turbospectrum/synthesis.py index cb07d49e..01b0274b 100644 --- a/smh/radiative_transfer/turbospectrum/synthesis.py +++ b/smh/radiative_transfer/turbospectrum/synthesis.py @@ -138,7 +138,6 @@ def synthesize(photosphere, transitions, abundances=None, isotopes=None, transitions_basename = "transitions.in" transitions.write(path(transitions_basename), format="turbospectrum") - # TODO: MOlecule files. n_files = ["DATA/Hlinedata", transitions_basename] @@ -151,13 +150,17 @@ def synthesize(photosphere, transitions, abundances=None, isotopes=None, proc_synth = subprocess.Popen(command, stdin=subprocess.PIPE, - stdout=subprocess.PIPE, stderr=subprocess.PIPE) + stdout=subprocess.PIPE, stderr=subprocess.PIPE, cwd=path("")) out_synth, err_synth = proc_synth.communicate(input=bsyn_lu_contents) errcode_synth = proc_synth.returncode + import matplotlib.pyplot as plt + spectrum = np.loadtxt(path("spectrum.out")) + fig, ax = plt.subplots() + ax.plot(spectrum[:,0], spectrum[:,1]) raise a # Get regions to synthesise. From 2e0151e291538ef2d8f5dbd3ecf95274325bcaa1 Mon Sep 17 00:00:00 2001 From: Andy Casey Date: Fri, 29 Sep 2017 20:42:04 +1000 Subject: [PATCH 09/14] add functionality to calculate equivalent widths for transitions using turbospectrum --- sandbox_turbospectrum.py | 34 +++- .../turbospectrum/__init__.py | 6 +- .../turbospectrum/bsyn_lu.in | 2 +- .../turbospectrum/equivalent_widths.py | 115 ++++++++++++ .../turbospectrum/eqwidt_lu.in | 24 +++ .../turbospectrum/synthesis.py | 173 ++++++++---------- 6 files changed, 246 insertions(+), 108 deletions(-) create mode 100644 smh/radiative_transfer/turbospectrum/equivalent_widths.py create mode 100644 smh/radiative_transfer/turbospectrum/eqwidt_lu.in diff --git a/sandbox_turbospectrum.py b/sandbox_turbospectrum.py index 2ba710f1..1ca35620 100644 --- a/sandbox_turbospectrum.py +++ b/sandbox_turbospectrum.py @@ -1,16 +1,40 @@ -from smh.radiative_transfer import turbospectrum as ts +import numpy as np + +from smh.radiative_transfer import (moog, turbospectrum as ts) from smh.photospheres.marcs import Interpolator from smh.linelists import LineList - +from time import time marcs = Interpolator() solar_photosphere = marcs(5777, 4.4, 0) -linelist = LineList.read("ni-test.list") +linelist = LineList.read("atomic.list") # Just restrict to a few lines. -#linelist = linelist[500:505] +linelist = linelist[500:505] + +t = [time()] +ts_result = ts.synthesize(solar_photosphere, linelist) + +t.append(time()) +moog_result = moog.synthesize(solar_photosphere, linelist) +t.append(time()) + +t_ts, t_moog = np.diff(t) + +ts_dispersion, ts_flux = ts_result[0][:2] +moog_dispersion, moog_flux = moog_result[0][:2] + + +import matplotlib.pyplot as plt +fig, ax = plt.subplots() + +ax.plot(ts_dispersion, ts_flux) +ax.plot(moog_dispersion, moog_flux) + +linelist["equivalent_width"] = np.array([ + 60, 0.1, 260, 90, 150]) -ts.synthesize(solar_photosphere, linelist, None, dispersion_min=5150, dispersion_max=5165) +result = ts.abundance_cog(solar_photosphere, linelist) \ No newline at end of file diff --git a/smh/radiative_transfer/turbospectrum/__init__.py b/smh/radiative_transfer/turbospectrum/__init__.py index fb3ea98e..316f2a1d 100644 --- a/smh/radiative_transfer/turbospectrum/__init__.py +++ b/smh/radiative_transfer/turbospectrum/__init__.py @@ -6,9 +6,13 @@ from __future__ import (division, print_function, absolute_import, unicode_literals) + #__all__ = ["abundance_cog", "synthesize", "RTError"] # See stackoverflow.com/questions/19913653/no-unicode-in-all-for-a-packages-init #__all__ = [_.encode("ascii") for _ in __all__] -from .synthesis import synthesize \ No newline at end of file +from .synthesis import synthesize +from .equivalent_widths import equivalent_widths +#from .cog import abundance_cog + diff --git a/smh/radiative_transfer/turbospectrum/bsyn_lu.in b/smh/radiative_transfer/turbospectrum/bsyn_lu.in index dbd489bb..7a718bda 100644 --- a/smh/radiative_transfer/turbospectrum/bsyn_lu.in +++ b/smh/radiative_transfer/turbospectrum/bsyn_lu.in @@ -5,7 +5,7 @@ 'COS(THETA):' '1.00' 'ABFIND :' '.false.' 'MODELOPAC:' '{opacities_path}' -'RESULTFILE:' '{spectrum_path}' +'RESULTFILE:' '{result_path}' 'METALLICITY:' '0.00' 'ALPHA/Fe :' '0.00' 'HELIUM :' '0.00' diff --git a/smh/radiative_transfer/turbospectrum/equivalent_widths.py b/smh/radiative_transfer/turbospectrum/equivalent_widths.py new file mode 100644 index 00000000..236a35ff --- /dev/null +++ b/smh/radiative_transfer/turbospectrum/equivalent_widths.py @@ -0,0 +1,115 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +""" +Functionality to derive equivalent widths from a transition's curve-of-growth using the +Turbospectrum radiative transfer code. +""" + +from __future__ import (division, print_function, absolute_import, + unicode_literals) + +import logging +import numpy as np +import os +import subprocess +import yaml +from pkg_resources import resource_stream + +from . import utils + +logger = logging.getLogger(__name__) + + +with resource_stream(__name__, "defaults.yaml") as fp: + _turbospectrum_defaults = yaml.load(fp) + + +def equivalent_widths(photosphere, transitions, full_output=False, verbose=False, + twd=None, **kwargs): + """ + Calculate equivalent widths from a curve-of-growth. + + :param photosphere: + A formatted photosphere. + + :param transitions: + A list of atomic transitions with measured equivalent widths. + + :param verbose: [optional] + Specify verbose flags. + """ + + kwds = _turbospectrum_defaults.copy() + kwds.update(kwargs) + + # Set default regions to synthesize. + kwds.setdefault("dispersion_min", + min(transitions["wavelength"]) - kwds["opacity_contribution"]) + kwds.setdefault("dispersion_max", + max(transitions["wavelength"]) + kwds["opacity_contribution"] \ + + kwds["dispersion_delta"]) + + # Update keywords. + xi = photosphere.meta.get("microturbulence", None) or 1.0 + kwds.update( + microturbulence=xi, + is_spherical=["F", "T"][photosphere.meta["radius"] > 0], + num_isotopes=0, formatted_isotopes="", + num_abundances=0, formatted_abundances="" + ) + + path = utils.twd_path(twd=twd, **kwargs) + + # Write atmosphere file. + photosphere.write(path(kwds["photosphere_path"]), format="turbospectrum") + + # Requires environment variable for the Turbospectrum data path. + os.symlink(os.environ.get("TURBODATA"), path("DATA")) + + # Calculate opacities. + op_proc = subprocess.Popen(["babsma_lu"], stdin=subprocess.PIPE, + stdout=subprocess.PIPE, stderr=subprocess.PIPE, + cwd=path("")) + + with resource_stream(__name__, "babsma_lu.in") as fp: + babsma_contents = fp.read() + + op_out, op_err = op_proc.communicate(input=babsma_contents.format(**kwds)) + + if verbose: + logger.info("Output from opacitices:\n{}".format(op_out)) + logger.warn("Error from opacities:\n{}".format(op_err)) + + if op_proc.returncode: + logging.exception( + "Exception when calculating opacities in Turbospectrum: {}".format( + op_err)) + raise ValueError("exception when calculating Turbospectrum opacities") + + transitions.write(path(kwds["transitions_path"]), format="turbospectrum") + + # TODO: Molecules + n_files = [kwds["transitions_path"]] + kwds.update(formatted_n_files="\n".join(n_files), n_files=len(n_files)) + + # Calculate the abundances. + proc = subprocess.Popen(["eqwidt_lu"], stdin=subprocess.PIPE, + stdout=subprocess.PIPE, stderr=subprocess.PIPE, cwd=path("")) + + with resource_stream(__name__, "eqwidt_lu.in") as fp: + bsyn_contents = fp.read() + + out, err = proc.communicate(input=bsyn_contents.format(**kwds)) + if proc.returncode: + logging.exception("Exception when calculating EWs in Turbospectrum:"\ + "\n{}".format(err)) + raise ValueError("exception when calculating EWs with Turbospectrum") + + transitions = transitions.copy() + transitions["equivalent_width"] = np.loadtxt( + path(kwds["result_path"]), usecols=(5, )).T + + return transitions + + diff --git a/smh/radiative_transfer/turbospectrum/eqwidt_lu.in b/smh/radiative_transfer/turbospectrum/eqwidt_lu.in new file mode 100644 index 00000000..7a718bda --- /dev/null +++ b/smh/radiative_transfer/turbospectrum/eqwidt_lu.in @@ -0,0 +1,24 @@ +'LAMBDA_MIN:' '{dispersion_min:.2f}' +'LAMBDA_MAX:' '{dispersion_max:.2f}' +'LAMBDA_STEP:' '{dispersion_delta:.2f}' +'INTENSITY/FLUX:' 'Flux' +'COS(THETA):' '1.00' +'ABFIND :' '.false.' +'MODELOPAC:' '{opacities_path}' +'RESULTFILE:' '{result_path}' +'METALLICITY:' '0.00' +'ALPHA/Fe :' '0.00' +'HELIUM :' '0.00' +'R-PROCESS :' '0.00' +'S-PROCESS :' '0.00' +'INDIVIDUAL ABUNDANCES:' '{num_abundances}' +{formatted_abundances} +'ISOTOPES:' '{num_isotopes}' +{formatted_isotopes} +'NFILES:' '{n_files}' +{formatted_n_files} +'SPHERICAL:' '{is_spherical}' + 30 + 300.00 + 15 + 1.30 diff --git a/smh/radiative_transfer/turbospectrum/synthesis.py b/smh/radiative_transfer/turbospectrum/synthesis.py index 01b0274b..039c048c 100644 --- a/smh/radiative_transfer/turbospectrum/synthesis.py +++ b/smh/radiative_transfer/turbospectrum/synthesis.py @@ -22,6 +22,24 @@ _turbospectrum_defaults = yaml.load(fp) +def _is_supported(): + + # Must have these things on path: + required_executables_on_path = ("babsma_lu", "bsyn_lu") + for each in required_executables_on_path: + try: + subprocess.check_output("which {}".format(each), shell=True) + + except subprocess.CalledProcessError: + raise OSError("cannot find executable {} on path".format(each)) + + # and environment variable $TURBODATA + if os.environ.get("TURBODATA") is None \ + or not os.path.exists(os.environ.get("TURBODATA")): + raise OSError("the TURBODATA environment variable must point to the "\ + "DATA directory for turbospectrum") + + return True @@ -39,6 +57,18 @@ def synthesize(photosphere, transitions, abundances=None, isotopes=None, """ + if abundances is None: + abundances = {} + + # TODO Put the abundances on the right scale???? + formatted_abundances = "\n".join(["{0} {1:.2f}".format(species, abundance) \ + for species, abundance in abundances.items()]) + + if isotopes is None: + isotopes = {} + + formatted_isotopes = "\n".join(["{0} {1:.2f}".format(mass_number, relative_abundance) \ + for mass_number, relative_abundance in isotopes.items()]) kwds = _turbospectrum_defaults.copy() kwds.update(kwargs) @@ -50,128 +80,69 @@ def synthesize(photosphere, transitions, abundances=None, isotopes=None, max(transitions["wavelength"]) + kwds["opacity_contribution"] \ + kwds["dispersion_delta"]) - # TODO put in abundances. if kwds["dispersion_max"] - kwds["dispersion_min"] > 1000.0: raise NotImplementedError("turbospectrum fails for large synth, " "and chunking not implemented yet") - path = utils.twd_path(twd=twd, **kwargs) - - - # Calculate opacities. - babsma_lu_kwds = kwds.copy() - babsma_lu_kwds["opacities_path"] = path("opacities.in") - babsma_lu_kwds["photosphere_path"] = path("photosphere.in") - - # Write atmosphere file. - photosphere.write(babsma_lu_kwds["photosphere_path"], format="turbospectrum") - - # Get babsma_lu template and format it. - with resource_stream(__name__, "babsma_lu.in") as fp: - babsma_lu_template = fp.read() - - if abundances is None: - abundances = {} - - # TODO Put the abundances on the right scale???? - - formatted_abundances = "\n".join(["{0} {1:.2f}".format(species, abundance) \ - for species, abundance in abundances.items()]) - - babsma_lu_kwds.update( + # Update keywords. + xi = photosphere.meta.get("microturbulence", None) or 1.0 + kwds.update( + microturbulence=xi, + is_spherical=["F", "T"][photosphere.meta["radius"] > 0], + num_isotopes=len(isotopes), + formatted_isotopes=formatted_isotopes, num_abundances=len(abundances), - formatted_abundances=formatted_abundances) - + formatted_abundances=formatted_abundances, + ) + path = utils.twd_path(twd=twd, **kwargs) - # XI - # num abundances - # formatted abundances. - - # If no XI value, then estimate it. - xi = photosphere.meta.get("microturbulence", None) - if xi is None: - xi = 1.0 # [km/s] - logger.warn("No microturbulence given in photosphere. Assuming 1 km/s.") - - # TODO microturbulence may have been specified by kwargs.... - if "microturbulence" in babsma_lu_kwds: - raise a - - babsma_lu_kwds["microturbulence"] = xi - - # Set some default values??? - - babsma_lu_contents = babsma_lu_template.format(**babsma_lu_kwds) - + # Write atmosphere file. + photosphere.write(path(kwds["photosphere_path"]), format="turbospectrum") - # Calcualte opacities. - command = kwds["turbospectrum_babsma_lu"] + # Requires environment variable for the Turbospectrum data path. + os.symlink(os.environ.get("TURBODATA"), path("DATA")) - os.symlink(kwds["turbospectrum_data_dir"], path("DATA")) - proc = subprocess.Popen(command.split(), stdin=subprocess.PIPE, + # Calculate opacities. + op_proc = subprocess.Popen(["babsma_lu"], stdin=subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.PIPE, cwd=path("")) - out, err = proc.communicate(input=babsma_lu_contents) - errcode = proc.returncode - - - # Calculate the spectrum. - os.symlink(kwds["turbospectrum_molecules_dir"], path("molecules")) - command = kwds["turbospectrum_bsyn_lu"] - - with resource_stream(__name__, "bsyn_lu.in") as fp: - bsyn_lu_template = fp.read() - - if isotopes is None: - isotopes = {} - formatted_isotopes = "\n".join(["{0} {1:.2f}".format(mass_number, relative_abundance) \ - for mass_number, relative_abundance in isotopes.items()]) - - - babsma_lu_kwds.update(spectrum_path=path("spectrum.out"), - num_isotopes=len(isotopes), formatted_isotopes=formatted_isotopes, - is_spherical=["F", "T"][photosphere.meta["radius"] > 0]) - - transitions_basename = "transitions.in" - transitions.write(path(transitions_basename), format="turbospectrum") - - - n_files = ["DATA/Hlinedata", transitions_basename] + with resource_stream(__name__, "babsma_lu.in") as fp: + babsma_contents = fp.read() - babsma_lu_kwds.update(formatted_n_files="\n".join(n_files), - n_files=len(n_files)) + op_out, op_err = op_proc.communicate(input=babsma_contents.format(**kwds)) + if op_proc.returncode: + logging.exception( + "Exception when calculating opacities in Turbospectrum: {}".format( + op_err)) + raise ValueError("exception when calculating Turbospectrum opacities") - bsyn_lu_contents = bsyn_lu_template.format(**babsma_lu_kwds) - command = kwds["turbospectrum_bsyn_lu"].split() + # Calculate the spectrum. + transitions.write(path(kwds["transitions_path"]), format="turbospectrum") + + # TODO: Molecules + n_files = ["DATA/Hlinedata", kwds["transitions_path"]] + kwds.update(formatted_n_files="\n".join(n_files), n_files=len(n_files)) - proc_synth = subprocess.Popen(command, stdin=subprocess.PIPE, + # Calculate the spectrum. + proc_synth = subprocess.Popen(["bsyn_lu"], stdin=subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.PIPE, cwd=path("")) - out_synth, err_synth = proc_synth.communicate(input=bsyn_lu_contents) - errcode_synth = proc_synth.returncode - - - import matplotlib.pyplot as plt - - spectrum = np.loadtxt(path("spectrum.out")) - fig, ax = plt.subplots() - ax.plot(spectrum[:,0], spectrum[:,1]) - - raise a - # Get regions to synthesise. - - - - # Update abundances - - # Calculate opacities + with resource_stream(__name__, "bsyn_lu.in") as fp: + bsyn_contents = fp.read() - # Synthesise + synth_out, synth_err = proc_synth.communicate(input=bsyn_contents.format(**kwds)) + if proc_synth.returncode: + logging.exception("Exception when calculating spectrum in Turbospectrum:"\ + "\n{}".format(synth_err)) + raise ValueError("exception when calculating spectrum with Turbospectrum") + dispersion, flux = np.loadtxt(path(kwds["result_path"]), usecols=(0, 1)).T + meta = {} + return [(dispersion, flux, meta)] From 8f1e68b0631284c95bc5fe42f4f27cba2f5c2350 Mon Sep 17 00:00:00 2001 From: Andy Casey Date: Fri, 29 Sep 2017 20:50:40 +1000 Subject: [PATCH 10/14] add function to calculate abundances from equivalent widths --- .../turbospectrum/__init__.py | 2 +- smh/radiative_transfer/turbospectrum/cog.py | 116 ++++++++++++++++++ .../turbospectrum/equivalent_widths.py | 5 +- .../turbospectrum/eqwidt_lu.in | 2 +- 4 files changed, 122 insertions(+), 3 deletions(-) create mode 100644 smh/radiative_transfer/turbospectrum/cog.py diff --git a/smh/radiative_transfer/turbospectrum/__init__.py b/smh/radiative_transfer/turbospectrum/__init__.py index 316f2a1d..9ec28c32 100644 --- a/smh/radiative_transfer/turbospectrum/__init__.py +++ b/smh/radiative_transfer/turbospectrum/__init__.py @@ -14,5 +14,5 @@ from .synthesis import synthesize from .equivalent_widths import equivalent_widths -#from .cog import abundance_cog +from .cog import abundance_cog diff --git a/smh/radiative_transfer/turbospectrum/cog.py b/smh/radiative_transfer/turbospectrum/cog.py new file mode 100644 index 00000000..0eb0446a --- /dev/null +++ b/smh/radiative_transfer/turbospectrum/cog.py @@ -0,0 +1,116 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +""" +Functionality to derive abundances from a transition's curve-of-growth using the +Turbospectrum radiative transfer code. +""" + +from __future__ import (division, print_function, absolute_import, + unicode_literals) + +import logging +import numpy as np +import os +import subprocess +import yaml +from pkg_resources import resource_stream + +from . import utils + +logger = logging.getLogger(__name__) + + +with resource_stream(__name__, "defaults.yaml") as fp: + _turbospectrum_defaults = yaml.load(fp) + + +def abundance_cog(photosphere, transitions, full_output=False, verbose=False, + twd=None, **kwargs): + """ + Calculate atomic line abundances by interpolating the measured + equivalent width from the curve-of-growth. + This wraps the MOOG `abfind` driver. + + :param photosphere: + A formatted photosphere. + + :param transitions: + A list of atomic transitions with measured equivalent widths. + + :param verbose: [optional] + Specify verbose flags to MOOG. This is primarily used for debugging. + """ + + kwds = _turbospectrum_defaults.copy() + kwds.update(kwargs) + + # Set default regions to synthesize. + kwds.setdefault("dispersion_min", + min(transitions["wavelength"]) - kwds["opacity_contribution"]) + kwds.setdefault("dispersion_max", + max(transitions["wavelength"]) + kwds["opacity_contribution"] \ + + kwds["dispersion_delta"]) + + # Update keywords. + xi = photosphere.meta.get("microturbulence", None) or 1.0 + kwds.update( + microturbulence=xi, + is_spherical=["F", "T"][photosphere.meta["radius"] > 0], + num_isotopes=0, formatted_isotopes="", + num_abundances=0, formatted_abundances="", + abfind="true" + ) + + path = utils.twd_path(twd=twd, **kwargs) + + # Write atmosphere file. + photosphere.write(path(kwds["photosphere_path"]), format="turbospectrum") + + # Requires environment variable for the Turbospectrum data path. + os.symlink(os.environ.get("TURBODATA"), path("DATA")) + + # Calculate opacities. + op_proc = subprocess.Popen(["babsma_lu"], stdin=subprocess.PIPE, + stdout=subprocess.PIPE, stderr=subprocess.PIPE, + cwd=path("")) + + with resource_stream(__name__, "babsma_lu.in") as fp: + babsma_contents = fp.read() + + op_out, op_err = op_proc.communicate(input=babsma_contents.format(**kwds)) + + if op_proc.returncode: + logging.exception( + "Exception when calculating opacities in Turbospectrum: {}".format( + op_err)) + raise ValueError("exception when calculating Turbospectrum opacities") + + transitions.write(path(kwds["transitions_path"]), format="turbospectrum") + + # TODO: Molecules + n_files = [kwds["transitions_path"]] + kwds.update(formatted_n_files="\n".join(n_files), n_files=len(n_files)) + + # Calculate the abundances. + proc = subprocess.Popen(["eqwidt_lu"], stdin=subprocess.PIPE, + stdout=subprocess.PIPE, stderr=subprocess.PIPE, cwd=path("")) + + with resource_stream(__name__, "eqwidt_lu.in") as fp: + bsyn_contents = fp.read() + + out, err = proc.communicate(input=bsyn_contents.format(**kwds)) + if proc.returncode: + logging.exception("Exception when calculating spectrum in Turbospectrum:"\ + "\n{}".format(err)) + raise ValueError("exception when calculating spectrum with Turbospectrum") + + # Order the transitions + transitions = transitions.copy() + transitions.sort(["species", "wavelength"]) + + transitions["abundance"] = np.loadtxt( + path(kwds["result_path"]), usecols=(11, ), skiprows=2) + + return transitions + diff --git a/smh/radiative_transfer/turbospectrum/equivalent_widths.py b/smh/radiative_transfer/turbospectrum/equivalent_widths.py index 236a35ff..ca4ebf0f 100644 --- a/smh/radiative_transfer/turbospectrum/equivalent_widths.py +++ b/smh/radiative_transfer/turbospectrum/equivalent_widths.py @@ -56,7 +56,8 @@ def equivalent_widths(photosphere, transitions, full_output=False, verbose=False microturbulence=xi, is_spherical=["F", "T"][photosphere.meta["radius"] > 0], num_isotopes=0, formatted_isotopes="", - num_abundances=0, formatted_abundances="" + num_abundances=0, formatted_abundances="", + abfind="false" ) path = utils.twd_path(twd=twd, **kwargs) @@ -107,6 +108,8 @@ def equivalent_widths(photosphere, transitions, full_output=False, verbose=False raise ValueError("exception when calculating EWs with Turbospectrum") transitions = transitions.copy() + transitions.sort(["species", "wavelength"]) + transitions["equivalent_width"] = np.loadtxt( path(kwds["result_path"]), usecols=(5, )).T diff --git a/smh/radiative_transfer/turbospectrum/eqwidt_lu.in b/smh/radiative_transfer/turbospectrum/eqwidt_lu.in index 7a718bda..2a878a7a 100644 --- a/smh/radiative_transfer/turbospectrum/eqwidt_lu.in +++ b/smh/radiative_transfer/turbospectrum/eqwidt_lu.in @@ -3,7 +3,7 @@ 'LAMBDA_STEP:' '{dispersion_delta:.2f}' 'INTENSITY/FLUX:' 'Flux' 'COS(THETA):' '1.00' -'ABFIND :' '.false.' +'ABFIND :' '.{abfind}.' 'MODELOPAC:' '{opacities_path}' 'RESULTFILE:' '{result_path}' 'METALLICITY:' '0.00' From efb7926e0a3fd2eaf69b711a0a0b4c4956c8968d Mon Sep 17 00:00:00 2001 From: Andy Casey Date: Fri, 29 Sep 2017 21:01:56 +1000 Subject: [PATCH 11/14] fix indexing problem that TS causes --- smh/radiative_transfer/turbospectrum/cog.py | 9 +++++++-- .../turbospectrum/equivalent_widths.py | 8 ++++++-- 2 files changed, 13 insertions(+), 4 deletions(-) diff --git a/smh/radiative_transfer/turbospectrum/cog.py b/smh/radiative_transfer/turbospectrum/cog.py index 0eb0446a..e6699625 100644 --- a/smh/radiative_transfer/turbospectrum/cog.py +++ b/smh/radiative_transfer/turbospectrum/cog.py @@ -109,8 +109,13 @@ def abundance_cog(photosphere, transitions, full_output=False, verbose=False, transitions = transitions.copy() transitions.sort(["species", "wavelength"]) - transitions["abundance"] = np.loadtxt( - path(kwds["result_path"]), usecols=(11, ), skiprows=2) + transitions["abundance"] = np.nan * np.ones(len(transitions)) + + wavelengths, abundances = np.loadtxt( + path(kwds["result_path"]), usecols=(2, 11, ), skiprows=2).T + + indices = np.where(np.in1d(transitions["wavelength"], wavelengths))[0] + transitions["abundance"][indices] = abundances return transitions diff --git a/smh/radiative_transfer/turbospectrum/equivalent_widths.py b/smh/radiative_transfer/turbospectrum/equivalent_widths.py index ca4ebf0f..4c11676d 100644 --- a/smh/radiative_transfer/turbospectrum/equivalent_widths.py +++ b/smh/radiative_transfer/turbospectrum/equivalent_widths.py @@ -109,9 +109,13 @@ def equivalent_widths(photosphere, transitions, full_output=False, verbose=False transitions = transitions.copy() transitions.sort(["species", "wavelength"]) + transitions["equivalent_width"] = np.nan * np.ones(len(transitions)) - transitions["equivalent_width"] = np.loadtxt( - path(kwds["result_path"]), usecols=(5, )).T + wavelengths, equivalent_widths = np.loadtxt( + path(kwds["result_path"]), usecols=(2, 5, )).T + + indices = np.where(np.in1d(transitions["wavelength"], wavelengths))[0] + transitions["equivalent_width"][indices] = equivalent_widths return transitions From c54cd822f85862675ba995ad450d05249d48a5f7 Mon Sep 17 00:00:00 2001 From: Andy Casey Date: Fri, 29 Sep 2017 21:07:50 +1000 Subject: [PATCH 12/14] Add defaults for turbospectrum --- smh/radiative_transfer/turbospectrum/defaults.yaml | 10 ++++++++++ 1 file changed, 10 insertions(+) create mode 100644 smh/radiative_transfer/turbospectrum/defaults.yaml diff --git a/smh/radiative_transfer/turbospectrum/defaults.yaml b/smh/radiative_transfer/turbospectrum/defaults.yaml new file mode 100644 index 00000000..400ac66e --- /dev/null +++ b/smh/radiative_transfer/turbospectrum/defaults.yaml @@ -0,0 +1,10 @@ +# The range in dispersion units (Angstroms) that each transition contributes. +opacity_contribution: 2.0 + +# The step size in dispersion points. +dispersion_delta: 0.01 + +photosphere_path: photosphere.in +opacities_path: opacities.in +transitions_path: transitions.in +result_path: result.out From 2513d5e2f68fc482c104dfe2cbedf4fdf15bc4b7 Mon Sep 17 00:00:00 2001 From: Andy Casey Date: Fri, 29 Sep 2017 21:11:56 +1000 Subject: [PATCH 13/14] make equivalent_widths and abundances_ return same format as moog.abundances --- smh/radiative_transfer/turbospectrum/cog.py | 7 +++---- .../turbospectrum/equivalent_widths.py | 11 ++++++----- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/smh/radiative_transfer/turbospectrum/cog.py b/smh/radiative_transfer/turbospectrum/cog.py index e6699625..146b7b20 100644 --- a/smh/radiative_transfer/turbospectrum/cog.py +++ b/smh/radiative_transfer/turbospectrum/cog.py @@ -109,13 +109,12 @@ def abundance_cog(photosphere, transitions, full_output=False, verbose=False, transitions = transitions.copy() transitions.sort(["species", "wavelength"]) - transitions["abundance"] = np.nan * np.ones(len(transitions)) - wavelengths, abundances = np.loadtxt( path(kwds["result_path"]), usecols=(2, 11, ), skiprows=2).T indices = np.where(np.in1d(transitions["wavelength"], wavelengths))[0] - transitions["abundance"][indices] = abundances - return transitions + _abundances = np.nan * np.ones(len(transitions)) + _abundances[indices] = abundances + return _abundances diff --git a/smh/radiative_transfer/turbospectrum/equivalent_widths.py b/smh/radiative_transfer/turbospectrum/equivalent_widths.py index 4c11676d..3ced3d32 100644 --- a/smh/radiative_transfer/turbospectrum/equivalent_widths.py +++ b/smh/radiative_transfer/turbospectrum/equivalent_widths.py @@ -109,14 +109,15 @@ def equivalent_widths(photosphere, transitions, full_output=False, verbose=False transitions = transitions.copy() transitions.sort(["species", "wavelength"]) - transitions["equivalent_width"] = np.nan * np.ones(len(transitions)) - - wavelengths, equivalent_widths = np.loadtxt( + + wavelengths, ews = np.loadtxt( path(kwds["result_path"]), usecols=(2, 5, )).T indices = np.where(np.in1d(transitions["wavelength"], wavelengths))[0] - transitions["equivalent_width"][indices] = equivalent_widths - return transitions + equivalent_widths = np.nan * np.ones(len(transitions)) + equivalent_widths[indices] = ews + + return equivalent_widths From fd6cc867787a0981991931ac025413a2587b3dd0 Mon Sep 17 00:00:00 2001 From: Andy Casey Date: Fri, 29 Sep 2017 21:23:26 +1000 Subject: [PATCH 14/14] update zenodo reference to corrected marcs atmospheres --- ni-test.list | 5 ----- sandbox_turbospectrum.py | 40 ---------------------------------------- setup.py | 6 +++--- 3 files changed, 3 insertions(+), 48 deletions(-) delete mode 100644 ni-test.list delete mode 100644 sandbox_turbospectrum.py diff --git a/ni-test.list b/ni-test.list deleted file mode 100644 index d57f352e..00000000 --- a/ni-test.list +++ /dev/null @@ -1,5 +0,0 @@ -5155.126 28.0 3.89800 -0.650 -5155.764 28.0 3.89800 0.074 -5157.707 28.0 3.70600 -1.264 -5157.980 28.0 3.60600 -1.510 -5162.916 28.0 3.83300 -3.200 \ No newline at end of file diff --git a/sandbox_turbospectrum.py b/sandbox_turbospectrum.py deleted file mode 100644 index 1ca35620..00000000 --- a/sandbox_turbospectrum.py +++ /dev/null @@ -1,40 +0,0 @@ - -import numpy as np - -from smh.radiative_transfer import (moog, turbospectrum as ts) -from smh.photospheres.marcs import Interpolator -from smh.linelists import LineList - -from time import time - - -marcs = Interpolator() -solar_photosphere = marcs(5777, 4.4, 0) -linelist = LineList.read("atomic.list") - -# Just restrict to a few lines. -linelist = linelist[500:505] - -t = [time()] -ts_result = ts.synthesize(solar_photosphere, linelist) - -t.append(time()) -moog_result = moog.synthesize(solar_photosphere, linelist) -t.append(time()) - -t_ts, t_moog = np.diff(t) - -ts_dispersion, ts_flux = ts_result[0][:2] -moog_dispersion, moog_flux = moog_result[0][:2] - - -import matplotlib.pyplot as plt -fig, ax = plt.subplots() - -ax.plot(ts_dispersion, ts_flux) -ax.plot(moog_dispersion, moog_flux) - -linelist["equivalent_width"] = np.array([ - 60, 0.1, 260, 90, 150]) - -result = ts.abundance_cog(solar_photosphere, linelist) \ No newline at end of file diff --git a/setup.py b/setup.py index bae2ee90..e3ff2646 100644 --- a/setup.py +++ b/setup.py @@ -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"),