PSOAP

Precision Spectroscopic Orbits A-Parametrically (pronounced ‘soap’)

PSOAP is a package to model astronomical spectra nonparametrically, for the purposes of disentangling spectra (as in double-lined spectroscopic binaries) or simply determining orbits (in the case of (single-lined spectroscopic binaries or exoplanet hosts). For more information about the mathematical framework underlying PSOAP, please see our paper

Disentangling Time Series Spectra with Gaussian Processes: Applications to Radial Velocity Analysis, Czekala et al., 2017ApJ...840...49C

PSOAP is actively developed on github here

Contents:

Getting Started

Installation

PSOAP requires a few packages standard to the scientific Python ecosystem. It is written and tested for currently maintained Python releases (Python 3.4+ as of Nov 2017); it has not been tested on the Python 2.x series and I do not anticipate it to work on it either.

• numpy
• scipy
• cython
• astropy
• h5py
• matplotlib
• celerite (optional)

All of these packages can be installed via an Anaconda Python installation or your normal means of managing your Python packages. Once you have installed them, clone the PSOAP package from the github repository

$ git clone https://github.com/iancze/PSOAP.git
$ cd PSOAP

and change to the top level PSOAP directory. Build the package via

$ python setup.py install

Which should build the cython extensions (used for faster matrix evaluations) and install the system scripts to your shell PATH. To check that you’ve got everything installed properly, try running from your shell

$ psoap-initialize --check
PSOAP successfully installed and linked.
Using Python Version 3.6.3 |Anaconda custom (64-bit)| (default, Nov  3 2017, 19:19:16)

If this doesn’t work, try double-checking the output from your install process to see if any errors popped up. If you are unable to fix these issues via the normal means of debugging python installs, please raise an issue with specifics about your system.

PSOAP has preliminary support for using the celerite package, which implements fast, one dimensional Gaussian processes which are used when fitting a single stationary star, or a single-lined spectroscopic binary or triple. You can optionally install this package following the link above. Unfortunately, this speedup is not available when fitting double or triple-lined spectroscopic binaries, though there may exist approximations which make this possible in the future.

Testing

If you really want to make sure everything works on your system, you can run the test suite by installing the pytest package, changing to the directory where you cloned the repository, and then running

$ py.test -v

If any of these tests fail, please report them by raising an issue with specifics about your system.

Citing

If you use our paper, code, or a derivative of it in your research, we would really appreciate a citation to Czekala et al. 2017

@ARTICLE{2017ApJ...840...49C,
    author = {{Czekala}, I. and {Mandel}, K.~S. and {Andrews}, S.~M. and {Dittmann}, J.~A. and
    {Ghosh}, S.~K. and {Montet}, B.~T. and {Newton}, E.~R.},
    title = "{Disentangling Time-series Spectra with Gaussian Processes: Applications to Radial Velocity Analysis}",
    journal = {\apj},
    archivePrefix = "arXiv",
    eprint = {1702.05652},
    primaryClass = "astro-ph.SR",
    keywords = {binaries: spectroscopic, celestial mechanics, stars: fundamental parameters, stars: individual: LP661-13, techniques: radial velocities, techniques: spectroscopic},
     year = 2017,
    month = may,
    volume = 840,
      eid = {49},
    pages = {49},
      doi = {10.3847/1538-4357/aa6aab},
    adsurl = {http://adsabs.harvard.edu/abs/2017ApJ...840...49C},
    adsnote = {Provided by the SAO/NASA Astrophysics Data System}
}

Copyright Ian Czekala and collaborators 2016-17.

Because the PSOAP package is still under occasional development, it may be wise to install the package in ‘development mode’ via the following commands:

$ pip install -e .
$ python setup.py build_ext --inplace

If you installed the package via development mode, then it in the case that the package is upgraded, it is easy to upgrade your local copy by simply pulling down the latest changes and rerunning the build script:

$ git pull
$ python setup.py build_ext --inplace

Scripts

The following scripts are made available on your command-line by the PSOAP package.

Initialization

!psoap-initialize --help

usage: psoap-initialize [-h] [--check] [--model {SB1,SB2,ST3}]

Initialize a new directory to do inference.

optional arguments:
  -h, --help            show this help message and exit
  --check               To help folks check whether the package was installed
                        properly.
  --model {SB1,SB2,ST3}
                        Which type of model to use, SB1, SB2, ST1, or SB3.

!psoap-generate-chunks --help

/usr/bin/sh: psoap-generate-chunks: command not found

!psoap-process-chunks --help

/usr/bin/sh: psoap-process-chunks: command not found

!psoap-generate-masks --help

/usr/bin/sh: psoap-generate-masks: command not found

!psoap-process-masks --help

/usr/bin/sh: psoap-process-masks: command not found

Sampling

The following scripts are used in sampling the posterior distribution.

!psoap-sample --help

You need to copy a config.yaml file to this directory, and then edit the values to your particular case.
Traceback (most recent call last):
  File "/home/ian/.build/anaconda/bin/psoap-sample", line 11, in <module>
    load_entry_point('psoap==0.0.1', 'console_scripts', 'psoap-sample')()
  File "/home/ian/.build/anaconda/lib/python3.6/site-packages/pkg_resources/__init__.py", line 570, in load_entry_point
    return get_distribution(dist).load_entry_point(group, name)
  File "/home/ian/.build/anaconda/lib/python3.6/site-packages/pkg_resources/__init__.py", line 2751, in load_entry_point
    return ep.load()
  File "/home/ian/.build/anaconda/lib/python3.6/site-packages/pkg_resources/__init__.py", line 2405, in load
    return self.resolve()
  File "/home/ian/.build/anaconda/lib/python3.6/site-packages/pkg_resources/__init__.py", line 2411, in resolve
    module = __import__(self.module_name, fromlist=['__name__'], level=0)
  File "/home/ian/.build/anaconda/lib/python3.6/site-packages/psoap-0.0.1-py3.6-linux-x86_64.egg/psoap/sample.py", line 33, in <module>
    f = open("config.yaml")
FileNotFoundError: [Errno 2] No such file or directory: 'config.yaml'

!psoap-sample-parallel --help

usage: psoap-sample-parallel [-h] [--debug] run_index

Sample the distribution across multiple chunks.

positional arguments:
  run_index   Which output subdirectory to save this particular run, in the
              case you may be running multiple concurrently.

optional arguments:
  -h, --help  show this help message and exit
  --debug     Print out debug commands to log.log

Configuration

PSOAP Spectra File Format

PSOAP uses HDF5 files to store the dataset of echelle spectra in a commonly used binary format. An example of spectra in this format is provided via the LP661-13 dataset, which you can download here. Since PSOAP is primarily designed to be used with high resolution spectra, which are commonly acquired with echelle spectrographs, it presumes an echelle-like format. If you have a regular spectrum, then you would just treat your dataset as an echelle with only one order.

Assumptions of this format (and of PSOAP):

• all spectra are taken with the same instrument
• all spectra have been corrected to the barycentric frame
• all orders of the echelle have the same number of pixels

The limitation on the same number of pixels per order is not necessarily a strict constraint. Let n_pix be the number of pixels in each echelle order, n_orders be the number of echelle orders of the spectrograph, and n_epochs be the number of epochs of data you have. If you have telescope data that does not have an equal number of n_pix for each order, you can choose the largest order as defining n_pix and then mask (see below) the blank regions in the shorter orders. Future versions of PSOAP may have the ability to incorporate data from multiple instruments, but for the foreseeable future it should be assumed that all data must be acquired from the same spectrograph (with the same spectral resolution/line spread function).

The HDF5 file has the following top-level entries as datasets:

• BCV: a length n_epochs array that provides the barycentric correction that was applied for each epoch. This value is only stored here for telluric/debugging purposes.
• JD: a length n_epochs array that provides the barycentric Julian Date for each epoch.
• `wl: a size (n_epochs, n_orders, n_pix) array that contains the pixel wavelengths of the spectrum, already corrected to the barycentric frame.
• fl: a size (n_epochs, n_orders, n_pix) array that contains the pixel fluxes of the spectrum.
• sigma: a size (n_epochs, n_orders, n_pix) array that contains the pixel flux uncertainties of the spectrum.

If you write scripts for converting spectra from your telescope, please consider adding them to this package via a pull request (into a scripts/your-telescope-name/ directory) so that they may be useful for other users with similar data. For example, there are already some scripts for converting data taken with the TRES spectrograph into the HDF5 format here.

Chunking the dataset

Masking

To start a new project, first decide what type of model you might like to use and navigate to a fresh directory. For example, let’s choose an SB2

Then, within your new directory

$ psoap-initialize --model SB2

This will copy the appropriate scripts to your directory. You are encouraged to use your favorite text editor and inspect the contents of the config.yaml file, which will be used frequently by this package.

One of the first things you should notice is the field data_file: data.hdf5.

Models

As PSOAP has grown, there a number of different possible models available to fit a number of different astrophysical systems.

Gravitational bodies

First, we use the following terminology to describe the number of gravitationally significant bodies in a system, meaning those bodies that must be factored into an orbital model. A single star would be a “spectroscopic single,” or SS, a binary star or a single star with an exoplanet would both be called a “spectroscopic binary,” or SB (hecklers—please excuse the terminology or suggest an improvement), a hierarchical triple star would be called a “spectroscopic triple,” or ST, and a hierarchical quadruple star (or a double binary) would be called a “spectroscopic quadruple,” or SQ. More complicated orbital hierarchies can be implemented, but we will need to think more critically about how to name these.

Spectroscopic bodies

Second, we use a number to describe the number of spectroscopically significant bodies in a system. In this case, a single star is 1, and a star with an exoplanet is also 1. A double-lined spectroscopic binary would be 2, and a triple-lined spectroscopic triple would be 3. However, a single-lined spectroscopic triple would be 1.

The final model specification is the concatenation of both of these. So a double lined spectroscopic binary is SB2, while a single-lined spectroscopic triple is ST1. Fortunately, the 1 models permit the usage of the extremely fast celerite framework by [Foreman-Mackey et al. 2017](http://adsabs.harvard.edu/abs/2017arXiv170309710F). However, if the system is a multiple star and there is in fact significant flux contribution from the other components, using a single-lined model can deliver biased results.

Specifying orbital models

Although one would think that we would only need the number of gravitational bodies to specify an orbit (and this is true), we do usually want to know how many spectroscopic bodies there are as well, since there are a different number of orbital parameters that we can constrain based upon how many spectroscopic bodies we see. Therefore orbit.py has orbital models for a wide range of models (e.g., ST1, SB2, etc...).

When this framework grows to include telluric models, or veiling models, these will be denoted with a +. For example, a single star with a telluric model would be denoted by SS1+T. Unfortunately these additional models are not able to use the fast celerite framework.

Note

This tutorial was generated from an IPython notebook that can be downloaded here.

Tutorial

This tutorial works through the example of fitting the LP661-13 dataset which appeared in Czekala et al. 2017. The spectra were originally acquired by Dittmann et al. 2017, and can be downloaded in HDF5 format **here**.

If you are looking to use data from a different telescope, you will need to process these spectra into a format like this HDF5 file. Some additional notes on how to do this are in processing your spectra to an HDF5 file.

This tutorial assumes that you have already followed the installation instructions.

Visualizing the dataset

Before we do any analysis with PSOAP, it’s a good idea to plot up all of your data. That way, we can see if there are any regions of the spectrum we may want to pay special attention to

!psoap_hdf5_exploder.py --help

usage: psoap_hdf5_exploder.py [-h] [--orders [ORDERS [ORDERS ...]]] [--SNR]
                              [--topo]

Make summary plots for a full HDF5 dataset.

optional arguments:
  -h, --help            show this help message and exit
  --orders [ORDERS [ORDERS ...]]
                        Which orders to plot. By default, all orders are
                        plotted. Can add more than one order in a spaced list,
                        e.g., --orders 22 23 24 but not --orders=22,23,24
  --SNR                 Plot spectra in order of highest SNR first, instead of
                        by date. Default is by date.
  --topo                Plot spectra in topocentric frame instead of
                        barycentric frame. Default is barycentric frame.

This will produce a bunch of plots in a newly-created plots directory.

Creating a configuration file

PSOAP generally relies upon a configuration text file for many of the project-specific settings. To create one from scratch, use the psoap_initialize.py command

!psoap_initialize.py --help

Using Python Version 3.6.1 |Anaconda 4.4.0 (64-bit)| (default, May 11 2017, 13:09:58)
[GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]
usage: psoap_initialize.py [-h] [--check] [--model {SB1,SB2,ST3}]

Initialize a new directory to do inference.

optional arguments:
  -h, --help            show this help message and exit
  --check               To help folks check whether the package was installed
                        properly.
  --model {SB1,SB2,ST3}
                        Which type of model to use, SB1, SB2, ST1, or SB3.

For this project, we’ll do

!psoap_initialize.py --model SB2

Open up the new config.yaml file in your directory with your favorite text editor, and familiarize yourself with the settings. For more information, check out :ref:configuration.

Setting up the chunks file

Because Gaussian processes are generally very computationally intensive, we’ll need to split the spectrum up into chunks so that it can be processed in parallel. The easiest way to get started is with psoap_generate_chunks.py

!psoap_generate_chunks.py --help

usage: psoap_generate_chunks.py [-h] [--pixels PIXELS] [--overlap OVERLAP]
                                [--start START] [--end END]

Auto-generate comprehensive chunks.dat file, which can be later edited by
hand.

optional arguments:
  -h, --help         show this help message and exit
  --pixels PIXELS    Roughly how many pixels should we keep in each chunk?
  --overlap OVERLAP  How many pixels of overlap to aim for.
  --start START      Starting wavelength.
  --end END          Ending wavelength.

Try running this command with the default values, and then open up the chunks.dat file that now exists in your local directory. You can try playing around with the specific values, but if you want to regenerate the file, you’ll need to delete the existing chunks.dat file from the directory first. To make things go quickly for this tutorial, we’re only going to use a limited section of the spectrum. Therefore, we’re going to open up chunks.dat and delete the chunks blueward of XX AA and redward of AA, leaving only 3 actual chunks. If you were doing this for real, you could choose your chunks more wisely. The inference procedure is set up so that it’s one chunk per CPU core, so generally feel free to use as many chunks as you have CPU cores, since there is no additional time penalty.

Orbit Routines

The psoap.orbit module is the backbone for the main PSOAP tasks of computing radial velocities to shift the wavelength vectors. For a description of the various radial velocity orbital models, see Models. There is also an additional module included within the PSOAP packaged, psoap.orbit_astrometry, which is not used for the main Gaussian process routines, but includes functions for jointly modeling radial velocity measurements and relative astrometric measurements.

psoap.orbit

Orbital conventions derived according to

\[\begin{split}v_{r,1} = K_1 (\cos (\omega_1 + f) + e \cos \omega_1) \\ v_{r,2} = K_2 (\cos (\omega_2 + f) + e \cos \omega_2)\end{split}\]

where a positive \(v_r\) indicates a redshifted source and \(\omega_2 = \omega_1 + \pi\). In an RV-only paper, typically authors report \(\omega = \omega_1\).

In general, quantities like vA, vB, and vC refer to radial velocities relative to the heliocentric frame of the earth, i.e., quantities to directly compare to actual radial velocity measurements.

(Private) quantities like v1, v2, v1_in, and v1_out have more subtle meanings—these are generally partial velocities used in the process of constructing vA, etc...

class psoap.orbit.SB1(K, e, omega, P, T0, gamma, obs_dates=None, **kwargs)

Orbital model for a single-line spectroscopic binary.

Parameters:

• K (float) – velocity semi-amplitude
• e (float) – eccentricity (must be between [0.0, 1.0))
• omega (float) – argument of periastron for the primary star (degrees)
• P (float) – period (days)
• T0 (float) – epoch of periastron (JD)
• gamma (float) – systemic velocity (km/s)
• obs_dates (1D np.array) – dates of observation (JD)

get_velocities(dates=None)

Calculate the velocities of the primary (vA) for all dates provided.

Parameters:dates (optional) – if provided, calculate velocities at this new vector of dates, rather than the one provided when the object was initialized. Returns:A (1, nvels) shape array of the primary velocities Return type:np.array

class psoap.orbit.SB2(q, K, e, omega, P, T0, gamma, obs_dates=None, **kwargs)

Orbital model for a double-line spectroscopic binary.

Parameters:

• q (float) – mass ratio M_B / M_A
• K (float) – velocity semi-amplitude
• e (float) – eccentricity (must be between [0.0, 1.0))
• omega (float) – argument of periastron for primary star (degrees)
• P (float) – period (days)
• T0 (float) – epoch of periastron (JD)
• gamma (float) – systemic velocity (km/s)
• obs_dates (1D np.array) – dates of observation (JD)

get_velocities(dates=None)

Calculate the velocities of the primary (vA) and secondary (vB) for all dates provided.

Parameters:dates (optional) – if provided, calculate velocities at this new vector of dates, rather than the one provided when the object was initialized. Returns:A (2, nvels) shape array of the primary and secondary velocities Return type:np.array

class psoap.orbit.ST1(K_in, e_in, omega_in, P_in, T0_in, K_out, e_out, omega_out, P_out, T0_out, gamma, obs_dates=None, **kwargs)

Orbital model for a single-line hierarchical spectroscopic triple.

Parameters:

• K_in (float) – velocity semi-amplitude for inner orbit
• e_in (float) – eccentricity for inner orbit (must be between [0.0, 1.0))
• omega_in (float) – argument of periastron for inner orbit (primary star) (degrees)
• P_in (float) – inner period (days)
• T0_in (float) – epoch of periastron for inner orbit (JD)
• K_out (float) – velocity semi-amplitude for outer orbit
• e_out (float) – eccentricity for outer orbit (must be between [0.0, 1.0))
• omega_out (float) – argument of periastron for outer orbit (primary stars) (degrees)
• P_out (float) – outer period (days)
• T0_out (float) – epoch of periastron for outer orbit (JD)
• gamma (float) – systemic velocity (km/s)
• obs_dates (1D np.array) – dates of observation (JD)

get_component_velocities(dates=None)

Routine to grab v1_inner and v1_outer.

Parameters:dates (optional) – if provided, calculate velocities at this new vector of dates, rather than the one provided when the object was initialized. Returns:A (1, nvels) shape array of the primary velocities for the outer motion Return type:np.array

get_velocities(dates=None)

Calculate the velocities of the primary (vA) for all dates provided.

Parameters:dates (optional) – if provided, calculate velocities at this new vector of dates, rather than the one provided when the object was initialized. Returns:A (1, nvels) shape array of the primary velocities Return type:np.array

class psoap.orbit.ST2(q_in, K_in, e_in, omega_in, P_in, T0_in, K_out, e_out, omega_out, P_out, T0_out, gamma, obs_dates=None, **kwargs)

Orbital model for a double-line (inner orbit) hierarchical spectroscopic triple.

Parameters:

• q_in (float) – mass ratio M_B / M_A
• K_in (float) – velocity semi-amplitude for inner orbit
• e_in (float) – eccentricity for inner orbit (must be between [0.0, 1.0))
• omega_in (float) – argument of periastron for inner orbit (degrees)
• P_in (float) – inner period (days)
• T0_in (float) – epoch of periastron for inner orbit (JD)
• K_out (float) – velocity semi-amplitude for outer orbit
• e_out (float) – eccentricity for outer orbit (must be between [0.0, 1.0))
• omega_out (float) – argument of periastron for outer orbit (degrees)
• P_out (float) – outer period (days)
• T0_out (float) – epoch of periastron for outer orbit (JD)
• gamma (float) – systemic velocity (km/s)
• obs_dates (1D np.array) – dates of observation (JD)

get_component_velocities(dates=None)

Routine to grab v1_inner and v1_outer.

Parameters:dates (optional) – if provided, calculate velocities at this new vector of dates, rather than the one provided when the object was initialized. Returns:A (1, nvels) shape array of the primary velocities for the outer motion Return type:np.array

get_velocities(dates=None)

Calculate the velocities of the primary (vA) and secondary (vB) for all dates provided.

Parameters:dates (optional) – if provided, calculate velocities at this new vector of dates, rather than the one provided when the object was initialized. Returns:A (2, nvels) shape array of the primary and secondary velocities Return type:np.array

class psoap.orbit.ST3(q_in, K_in, e_in, omega_in, P_in, T0_in, q_out, K_out, e_out, omega_out, P_out, T0_out, gamma, obs_dates=None, **kwargs)

Orbital model for a double-line (inner orbit) hierarchical spectroscopic triple.

Parameters:

• q_in (float) – mass ratio M_B / M_A
• K_in (float) – velocity semi-amplitude for inner orbit
• e_in (float) – eccentricity for inner orbit (must be between [0.0, 1.0))
• omega_in (float) – argument of periastron for inner orbit (degrees)
• P_in (float) – inner period (days)
• T0_in (float) – epoch of periastron for inner orbit (JD)
• q_out (float) – mass ratio M_C / (M_A + M_B)
• K_out (float) – velocity semi-amplitude for outer orbit
• e_out (float) – eccentricity for outer orbit (must be between [0.0, 1.0))
• omega_out (float) – argument of periastron for outer orbit (degrees)
• P_out (float) – outer period (days)
• T0_out (float) – epoch of periastron for outer orbit (JD)
• gamma (float) – systemic velocity (km/s)
• obs_dates (1D np.array) – dates of observation (JD)

get_component_velocities(dates=None)

Routine to grab v1_inner and v1_outer.

Parameters:dates (optional) – if provided, calculate velocities at this new vector of dates, rather than the one provided when the object was initialized. Returns:A (1, nvels) shape array of the primary velocities for the outer motion Return type:np.array

get_velocities(dates=None)

Calculate the velocities of the primary (vA), secondary (vB), and tertiary (vC) for all dates provided.

Parameters:dates (optional) – if provided, calculate velocities at this new vector of dates, rather than the one provided when the object was initialized. Returns:A (3, nvels) shape array of the primary, secondary, and tertiary velocities Return type:np.array

psoap.orbit_astrometry

This is an alternate parameterization suited for joint astrometric and radial velocity orbital modeling.

Orbital conventions derived according to

\[\begin{split}\Delta \delta = X = \rho \cos \theta \\ \Delta \alpha \cos \delta = Y = \rho \sin \theta\end{split}\]

The positive X-axis points North, and the positive Y-axis points East. Rotation angle \(\theta\) is measured in degrees East of North, i.e.

\[\theta = \arctan(Y / X).\]

In our implementation, we use np.atan2(Y/X) to resolve the quadrant ambiguity.

The ascending node is defined as the point where the secondary crosses the plane of the sky receeding from the observer. It is measured in degrees East of North.

Inclination is defined such that \(0 < i < \pi/2\) yield prograde orbits (counter-clockwise, such that \(\theta\) increases), and \(\pi/2 < i < \pi\) yield retrograde orbits (clockwise, such that \(\theta\) decreases).

These lead to the following equations of motion

\[\begin{split}X = r (\cos \Omega \cos(\omega + f) - \sin(\Omega) \sin(\omega + f) \cos(i)) \\ Y = r (\sin \Omega \cos(\omega + f) + \cos(\Omega) \sin(\omega + f) \cos(i)) \\ Z = - r \sin(\omega + f) \sin(i)\end{split}\]

and

\[\begin{split}v_{r,1} = K_1 (\cos (\omega_1 + f) + e \cos \omega_1) \\ v_{r,2} = K_2 (\cos (\omega_2 + f) + e \cos \omega_2).\end{split}\]

class psoap.orbit_astrometry.Binary(a, e, i, omega, Omega, T0, M_tot, M_2, gamma, obs_dates=None, **kwargs)

Binary orbital model that can deliver absolute astrometric position, relative astrometric position (B relative to A), and radial velocities for A and B.

Parameters:

• a (float) – semi-major axis [AU]
• e (float) – eccentricity (must be between [0.0, 1.0))
• i (float) – inclination [deg]
• omega (float) – argument of periastron of the primary, i.e. \(\omega_1\) [degrees]
• Omega (float) – position angle of the ascending node (going into sky) [deg] east of north
• T0 (float) – epoch of periastron passage [JD]
• M_tot (float) – sum of the masses \(M_\mathrm{A} + M_\mathrm{B}\) [\(M_\odot\)]
• M_2 (float) – mass of B [\(M_\odot\)]
• gamma (float) – systemic velocity (km/s)
• obs_dates (1D np.array) – dates of observation (JD)

get_full_orbit(dates=None)

Deliver the full set of astrometric and radial velocity quantities, namely the radial velocities vA, vB, the position of A and B relative to the center of mass in the plane of the sky (XY_A and XY_B, respectively), the position of B relative to the position of A in the plane of the sky (XY_AB), the position of A and B in the plane of the orbit (xy_A and xy_B, respectively), and the position of B relative to the position of A in the plane of the orbit (xy_AB), for all dates provided. All positions are given in units of AU, and so must be converted to arcseconds after assuming a distance to the system.

Parameters:dates (optional) – if provided, calculate quantities at this new vector of dates, rather than the one provided when the object was initialized. Returns:A dictionary with items of "vAs", "vBs", "XYZ_As", "XYZ_Bs", "XYZ_ABs", "xy_As", "xy_Bs", "xy_ABs" Return type:dict

get_orbit(dates=None)

Deliver only the main quantities useful for performing a joint astrometric + RV fit to real data, namely the radial velocities vA, vB, the relative offsets \(\rho\), and relative position angles \(\theta\), for all dates provided. Relative offsets are provided in AU, and so must be converted to arcseconds after assuming a distance to the system. Relative position angles are given in degrees east of north.

Parameters:dates (optional) – if provided, calculate quantities at this new vector of dates, rather than the one provided when the object was initialized. Returns:A dictionary with items "vAs", "vBs", "rhos", "thetas". Return type:dict

class psoap.orbit_astrometry.Triple(a_in, e_in, i_in, omega_in, Omega_in, T0_in, a_out, e_out, i_out, omega_out, Omega_out, T0_out, M_1, M_2, M_3, gamma, obs_dates=None, **kwargs)

Triple orbital model that can deliver absolute astrometric position, relative astrometric position (B relative to A, and C relative to A), and radial velocities for A, B, and C.

Parameters:

• a_in (float) – semi-major axis for inner orbit [AU]
• e_in (float) – eccentricity for inner orbit (must be between [0.0, 1.0))
• i_in (float) – inclination for inner orbit [deg]
• omega_in (float) – argument of periastron for inner orbit, for the primary (\(\omega_1\)) [degrees]
• Omega_in (float) – position angle of the ascending node [deg] east of north for inner orbit
• T0_in (float) – epoch of periastron passage for inner orbit [JD]
• a_out (float) – semi-major axis for outer orbit [AU]
• e_out (float) – eccentricity for outer orbit (must be between [0.0, 1.0))
• i_out (float) – inclination for outer orbit [deg]
• omega_out (float) – argument of periastron for outer orbit, for the “primary” (\(\omega_{12}\)) [degrees]
• Omega_out (float) – position angle of the ascending node [deg] east of north for outer orbit
• T0_out (float) – epoch of periastron passage for outer orbit [JD]
• M_1 (float) – mass of A [\(M_\odot\)]
• M_2 (float) – mass of B [\(M_\odot\)]
• M_3 (float) – mass of C [\(M_\odot\)]
• gamma (float) – systemic velocity (km/s)
• obs_dates (1D np.array) – dates of observation (JD)

get_full_orbit(dates=None)

Deliver the full set of astrometric and radial velocity quantities, namely the radial velocities \(v_{r,A}\), \(v_{r,B}\), \(v_{r,C}\), the position of A, B, and C relative to the center of mass in the plane of the sky (XYZ_A, XYZ_B, and XYZ_C, respectively), the absolute position of the center of mass of (AB), (XYZ_AB), the position of A relative to the center of mass of AB (XYZ_A_loc), the position of B relative to the center of mass of (AB) (XYZ_B_loc), the absolute position of C in the plane of the orbit (xy_C), the absolute positon of the center of mass of AB in the plane of the orbit (xy_AB), the position of A in the plane of the orbit, relative to the center of mass of AB (xy_A_locs), and the position of B in the plane of the orbit, relative to the center of mass of AB (xy_B_locs), for all dates provided. All positions are given in units of AU, and so must be converted to arcseconds after assuming a distance to the system.

Parameters:dates (optional) – if provided, calculate quantities at this new vector of dates, rather than the one provided when the object was initialized. Returns:A dictionary with keys vAs, vBs, vCs, XYZ_As, XYZ_Bs, XYZ_Cs, XYZ_ABs, XYZ_A_locs, XYZ_B_locs, xy_Cs, xy_ABs, xy_A_locs, xy_B_locs Return type:dict

get_orbit(dates=None)

Deliver only the main quantities useful for performing a joint astrometric + RV fit to real data, namely the radial velocities \(v_{r,A}\), \(v_{r,B}\), \(v_{r,C}\), the relative offsets of B to A \(\rho_\mathrm{AB}\), and relative position angles \(\theta_\mathrm{AB}\), and the relative offsets of C to A \(\rho_\mathrm{AC}\) and \(\theta_\mathrm{AC}\) for all dates provided. Relative offsets are provided in AU, and so must be converted to arcseconds after assuming a distance to the system. Relative position angles are given in degrees east of north.

Parameters:dates (optional) – if provided, calculate quantities at this new vector of dates, rather than the one provided when the object was initialized. Returns:A dictionary with keys "vAs", "vBs", "vCs", "rho_ABs", "theta_ABs", "rho_ACs", "theta_ACs" Return type:dict

Data

Data formats

PSOAP relies upon chunks of data. When working with real data, there are a few things to keep in mind.

First, it may so happen that certain pixels may need to be masked, for example due to cosmic ray hits. This means that actual data chunks will probably have an un-equal number of pixels per epoch. This is OK.

So, all chunks will be generated with their full complement of data, but when executing any inference routines, the masks will be applied to the data.

All data and chunks are stored in an HDF5 format.

Data module

class psoap.data.Chunk(wl, fl, sigma, date, mask=None)

Hold a chunk of data. Each chunk is shape (n_epochs, n_pix) and has components wl, fl, sigma, date, and mask (all the same length).

apply_mask()

Apply the mask to all of the attributes, so now we return 1D arrays.

date = None

date vector

date1D = None

data vector of length n_epochs

fl = None

flux vector

lwl = None

natural log of the wavelength vector

classmethod open(order, wl0, wl1, limit=100, prefix='')

Load a spectrum from a directory link pointing to HDF5 output. :param fname: HDF5 file containing files on disk.

sigma = None

measurement uncertainty vector

wl = None

wavelength vector

class psoap.data.Spectrum(fname)

Data structure for the raw spectra, stored in an HDF5 file.

This is the main datastructure used to interact with your dataset. The key is getting your spectra into an HDF5 format first.

Parameters:fname (string) – location of the HDF5 file. Returns:the instantiated Spectrum object. Return type:Spectrum

sort_by_SN(order=22)

Sort the dataset in order of decreasing signal to noise. This is designed to make it easy to limit the analysis to the highest SNR epochs, if you wish to speed things up.

Parameters:order (int) – the order to calculate the signal-to-noise. By default, the TRES Mg b order is chosen, which is generally a good order for TRES data. If you are using data from a different telescope, you will likely need to adjust this value.

psoap.data.lredshift(lwl, v)

Redshift a vector of wavelengths that are already in log-lamba (natural log). A positive velocity corresponds to a lengthening (increase) of the wavelengths in the array.

Parameters:

• wl (np.array, arbitrary shape) – the input ln(wavelengths).
• velocity (float) – the velocity by which to redshift the wavelengths

Returns:

A redshifted version of the wavelength vector

Return type:

np.array

psoap.data.redshift(wl, v)

Redshift a vector of wavelengths. A positive velocity corresponds to a lengthening (increase) of the wavelengths in the array.

Parameters:

• wl (np.array, arbitrary shape) – the input wavelengths
• velocity (float) – the velocity by which to redshift the wavelengths

Returns:

A redshifted version of the wavelength vector

Return type:

np.array

psoap.data.replicate_wls(lwls, velocities, mask)

Using the set of velocities calculated from an orbit, copy and blue-shift the input ln(wavelengths), so that they correspond to the rest-frame wavelengths of the individual components. This routine is primarily for producing replicated ln-wavelength vectors ready to feed to the GP routines.

Parameters:

• lwls (1D np.array with length (n_epochs * n_good_pixels)) – this dataproduct is the 1D representation of the natural log of the (masked) input wavelength vectors. The masking process naturally makes it 1D.
• velocities (2D np.array with shape (n_components, n_epochs)) – a set of velocities determined from an orbital model.
• mask – the np.bool mask used to select the good datapoints. It is necessary for properly replicating the velocities to the right epoch.

Returns:

A 2D (n_components, n_epochs * n_good_pixels) shape array of the wavelength vectors blue-shifted according to the velocities. This means that for each component, the arrays are flattened into 1D vectors.

Return type:

np.array

Utils module

psoap.utils.convert_dict(model, fix_params, **kwargs)

Used to turn a dictionary of parameter values (from config.yaml) directly into a parameter type. Generally used for synthesis and plotting command line scripts.

psoap.utils.convert_vector(p, model, fix_params, **kwargs)

Unroll a vector of parameter values into a parameter type, using knowledge about which model we are fitting, the parameters we are fixing, and the default values of those parameters.

Parameters:

• p (np.float) – 1D input array of only a subset of parameter values.
• model (str) – “SB1”, “SB2”, etc..
• fix_params (list of str) – names of parameters that will be fixed
• **kwargs – input for {param_name: default_value} pairs

Returns:

a 2-tuple of the full vectors for the orbital parameters, and the GP parameters, augmented with the previously missing values.

Return type:

(np.float, np.float)

psoap.utils.gelman_rubin(samplelist)

Given a list of flatchains from separate runs (that already have burn in cut and have been trimmed, if desired), compute the Gelman-Rubin statistics in Bayesian Data Analysis 3, pg 284. If you want to compute this for fewer parameters, then slice the list before feeding it in.

psoap.utils.get_labels(model, fix_params)

Collect the labels for each model, so that we can plot.

Covariance Routines

In the original paper, the kernel was specified using the velocity distance between two wavelengths (\(\lambda_i\), \(\lambda_j\)). In subsequent versions, we choose to simply use the absolute magnitude of the distance between the natural log of the wavelengths, which is functionally equivalent to the former definition, but makes the computation easier. We let

\[\zeta_i = \ln \lambda_i\]

\[r_{ij} = |\zeta_i - \zeta_j|\]

Of course, when we are fitting for two spectral components, f and g, then there will be two sets of input wavelength vectors and thus two distances

\[\begin{split}r_{f,ij} = |\zeta_{f,i} - \zeta_{f,j}| \\ r_{g,ij} = |\zeta_{g,i} - \zeta_{g,j}|\end{split}\]

With two stars, to evaluate the kernel for a specific element in the matrix actually requires four inputs

\[k(\zeta_{f,i}, \zeta_{f,j}, \zeta_{g,i}, \zeta_{g,j} | a_f, l_f, a_g, l_g) = a_f^2 \exp \left ( - \frac{|\zeta_{f,i} - \zeta_{f,j}|^2}{2 l_f^2} \right) + a_g^2 \exp \left ( - \frac{|\zeta_{g,i} - \zeta_{g,j}|^2}{2 l_g^2} \right)\]

psoap.covariance.cycle_calibration(wl, fl, sigma, amp_f, l_f, ncycles, order=1, limit_array=3, mu_GP=1.0, soften=1.0)

Given a chunk of spectra, cycle n_cycles amongst all spectra and return the spectra with inferred calibration adjustments.

order : what order of Chebyshev polynomials to use. 1st order = line.

Only use limit_array number of spectra to save memory.

psoap.covariance.cycle_calibration_chunk(chunk, amp_f, l_f, n_cycles, order=1, limit_array=3, mu_GP=1.0, soften=1.0)

Do the calibration on a chunk at at time, incorporating the masks.

psoap.covariance.lnlike_f(V11, wl_f, fl, sigma, amp_f, l_f, mu_GP=1.0)

Calculate the log-likelihood for a single-lined spectrum.

This function takes a pre-allocated array and fills out the covariance matrices and evaluates the likelihood function for a single-lined spectrum, assuming a squared-exponential kernel (does not celerite).

Parameters:

• V11 (numpy 2D array) – Description of arg1
• wl_f (numpy 1D array) – Description of arg2
• fl (numpy 1D array) – ae
• amp_f (float) – amplitude of GP
• l_f (float) – length scale of GP
• mu_GP (float) – mean of GP

Returns:

The log-likelihood value

Return type:

float

psoap.covariance.lnlike_f_g(V11, wl_f, wl_g, fl, sigma, amp_f, l_f, amp_g, l_g, mu_GP=1.0)

V11 is a matrix to be allocated. wl_known, fl_known, and sigma_known are flattened 1D arrays.

psoap.covariance.lnlike_f_g_george(lwl_f, lwl_g, fl, sigma, amp_f, l_f, amp_g, l_g, mu_GP=1.0)

Evaluate the joint likelihood for f and g using George.

psoap.covariance.lnlike_f_g_h(V11, wl_f, wl_g, wl_h, fl, sigma, amp_f, l_f, amp_g, l_g, amp_h, l_h, mu_GP=1.0)

V11 is a matrix to be allocated. wl_known, fl_known, and sigma_known are flattened 1D arrays.

psoap.covariance.optimize_GP_f(wl_known, fl_known, sigma_known, amp_f, l_f, mu_GP=1.0)

Optimize the GP hyperparameters for the given slice of data. Amp and lv are starting guesses.

psoap.covariance.optimize_calibration(lwl0, lwl1, lwl_cal, fl_cal, fl_fixed, A, B, C, order=1, mu_GP=1.0)

Determine the calibration parameters for this epoch of observations. This is a more general method than psoap.covariance.optimize_calibration_static(), since it allows arbitrary covariance matrices, which should be used when there is orbital motion. Assumes that covariance matrices are appropriately filled out.

Parameters:

• lwl0 (float) – left side evaluation point for Chebyshev
• lwl1 (float) – right side evaluation point for Chebyshev
• lwl_cal (np.array) – the wavelengths corresponding to the epoch we want to calibrate
• fl_cal (np.array) – the fluxes corresponding to the epoch we want to calibrate
• fl_fixed (np.array) – the remaining epochs of data to calibrate in reference to.
• A (2D np.array) – matrix_functions.fill_V11_f(A, lwl_cal, amp, l_f) with sigma_cal already added to the diagonal
• B (2D np.array) – matrix_functions.fill_V11_f(B, lwl_fixed, amp, l_f) with sigma_fixed already added to the diagonal
• C (2D np.array) – matrix_functions.fill_V12_f(C, lwl_cal, lwl_fixed, amp, l_f) cross matrix (with no sigma added, since these are independent measurements).
• order (int) – the degree polynomial to use. order = 1 is a line, order = 2 is a line + parabola

Returns:

a tuple of two data products. The first is the fl_cal vector, now calibrated. The second is the array of the Chebyshev coefficients, in case one wants to re-evaluate the calibration polynomials.

Return type:

(np.array, np.array)

psoap.covariance.optimize_calibration_ST1(lwl0, lwl1, lwl_cal, fl_cal, fl_fixed, gp, A, C, mu_GP=1.0, order=1)

Determine the calibration parameters for this epoch of observations.

lwl0, lwl1: set the points for the Chebyshev.

This is a more general method than optimize_calibration_static, since it allows arbitrary covariance matrices, which should be used when there is orbital motion.

lwl_cal: the wavelengths corresponding to the epoch we want to calibrate fl_cal: the fluxes corresponding to the epoch we want to calibrate

fl_fixed: the remaining epochs of data to calibrate in reference to.

gp: the celerite GP

order: the degree polynomial to use. order = 1 is a line, order = 2 is a line + parabola

Assumes that covariance matrices are appropriately filled out.

psoap.covariance.optimize_calibration_static(wl0, wl1, wl_cal, fl_cal, sigma_cal, wl_fixed, fl_fixed, sigma_fixed, amp, l_f, order=1, mu_GP=1.0)

Determine the calibration parameters for this epoch of observations. Assumes all wl, fl arrays are 1D, and that the relative velocities between all epochs are zero.

Parameters:

• wl0 (float) – left wl point to evaluate the Chebyshev
• wl1 (float) – right wl point to evaluate the Chebyshev
• wl_cal (np.array) – the wavelengths of the epoch to calibrate
• fl_cal (np.array) – the fluxes of the epoch to calibrate
• sigma_cal (np.array) – the sigmas of the epoch to calibrate
• wl_fixed (np.array) – the 1D (flattened) array of the reference wavelengths
• fl_fixed (np.array) – the 1D (flattened) array of the reference fluxes
• sigma_fixed (np.array) – the 1D (flattened) array of the reference sigmas
• amp (float) – the GP amplitude
• l_f (float) – the GP length
• order (int) – the Chebyshev order to use
• mu_GP (optional) – the mean of the GP to assume.

Returns:

a tuple of two data products. The first is the fl_cal vector, now calibrated. The second is the array of the Chebyshev coefficients, in case one wants to re-evaluate the calibration polynomials.

Return type:

(np.array, np.array)

psoap.covariance.optimize_epoch_velocity_f(lwl_epoch, fl_epoch, sigma_epoch, lwl_fixed, fl_fixed, sigma_fixed, gp)

Optimize the wavelengths of the chosen epoch relative to the fixed wavelengths. Identify the velocity required to redshift the chosen epoch.

psoap.covariance.predict_f(lwl_known, fl_known, sigma_known, lwl_predict, amp_f, l_f, mu_GP=1.0)

wl_known are known wavelengths. wl_predict are the prediction wavelengths. Assumes all inputs are 1D arrays.

psoap.covariance.predict_f_g(lwl_f, lwl_g, fl_fg, sigma_fg, lwl_f_predict, lwl_g_predict, mu_f, amp_f, l_f, mu_g, amp_g, l_g, get_Sigma=True)

Given that f + g is the flux that we’re modeling, jointly predict the components.

psoap.covariance.predict_f_g_h(lwl_f, lwl_g, lwl_h, fl_fgh, sigma_fgh, lwl_f_predict, lwl_g_predict, lwl_h_predict, mu_f, mu_g, mu_h, amp_f, l_f, amp_g, l_g, amp_h, l_h)

Given that f + g + h is the flux that we’re modeling, jointly predict the components.

psoap.covariance.predict_f_g_h_sum(lwl_f, lwl_g, lwl_h, fl_fgh, sigma_fgh, lwl_f_predict, lwl_g_predict, lwl_h_predict, mu_fgh, amp_f, l_f, amp_g, l_g, amp_h, l_h)

Given that f + g + h is the flux that we’re modeling, predict the joint sum.

psoap.covariance.predict_python(wl_known, fl_known, sigma_known, wl_predict, amp_f, l_f, mu_GP=1.0)

wl_known are known wavelengths. wl_predict are the prediction wavelengths.

This is the covariance module.

Sampling the Parameters

Options to use SB2 and SB1. Show an example using the fast celerite GPs for single lined systems, including options to determine radial velocity on a per-epoch basis. Also strong warnings about contamination from secondary light.

Reconstructing Spectra

Once you have determined appropriate parameters for your orbit and Gaussian processes, then you can reconstruct the orbit.

Changelog

v0.1.1

standardized different model selections

Now models can be specified as combinations of the number of gravitationally significant bodies, and the number of spectroscopically significant bodies. See models.md for more information. This has an impact for the way orbit.py is used.

Calibration optimization

For now, only implemented for the ST3 model.

log-lambda input coordinates

Using log-lambda input coordinates instead of lambda enables a small but not insignificant speedup in the calculation of the kernel spacing, since the calculation of velocity spacing becomes a simple subtraction.

The following modules and scripts have been converted

• psoap_sample_parallel.py
• covariance.py
• matrix_functions.pyx
• psoap_retrieve_SB2.py
• psoap_retrieve_ST3.py

v0.1.0

Beta Release corresponding to code used for paper on the arXiv.

Indices and tables

• Index
• Module Index
• Search Page