Merge branch 'develop'

CHANGELOG.rst

@@ -0,0 +1,29 @@
+CHANGELOG
+=========
+
+* All notable changes to this project will be documented in this file.
+* This project adheres to `Semantic Versioning <http://semver.org/>`_.
+* This change log is adapted from to `Keep a CHANGELOG <http://keepachangelog.com/>`_.
+
+[v0.1.3] - UNRELEASED
+---------------------
+Added
+^^^^^
+* `Issue 20 <https://github.com/ccd-utexas/binstarsolver/issues/20>`_: Added CHANGELOG.rst
+
+Fixed
+^^^^^
+* `Issue 19 <https://github.com/ccd-utexas/binstarsolver/issues/19>`_: Solver for inclination must give inclination < 90 deg.
+* `Issue 21 <https://github.com/ccd-utexas/binstarsolver/issues/21>`_: Solver for inclination did not converge for special cases.
+
+[v0.1.2] - 20150123T172000
+--------------------------
+Added
+^^^^^
+* Initial release to PyPI.
+
+[v0.1.1] - 20150123T140000
+-----------------------
+Added
+^^^^^
+* Initial release.

CONTRIBUTING.rst

@@ -8,3 +8,4 @@ Thanks for your interest in contributing! Please use the following:
 * Docstring style: https://github.com/numpy/numpy/blob/master/doc/example.py
 * Coding style: https://google-styleguide.googlecode.com/svn/trunk/pyguide.html
 * Versioning: http://semver.org/
+* Change log: http://keepachangelog.com/

binstarsolver/__init__.py

@@ -1,10 +1,20 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
 """Calculate physical quantities of a binary star system from observed quantities.
 
-Docstrings are adapted from `numpy` convention:
-https://github.com/numpy/numpy/blob/master/doc/example.py
+Notes
+-----
+Docstrings are adapted from `numpy` convention [1]_
+
+References
+----------
+..[1] https://github.com/numpy/numpy/blob/master/doc/example.py
 
 """
 
+
+# Import standard packages.
 from __future__ import absolute_import, division, print_function
+# Import local packages.
 from . import main
 from . import utils

binstarsolver/main.py

@@ -1,8 +1,13 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
 """Top-level script for calculating physical quantities from observed quantities.
 
 """
 
+
+# Import standard packages.
 from __future__ import absolute_import, division, print_function
+# Import local packages.
 from . import utils
 
 

binstarsolver/utils.py

@@ -1,9 +1,15 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
 """Utilities for calculating physical quantities from observed quantities.
 
 """
 
+
+# Import standard packages.
 from __future__ import absolute_import, division, print_function
 import sys
+import warnings
+# Import installed packages.
 import numpy as np
 import scipy.optimize as sci_opt
 import scipy.constants as sci_con
@@ -259,10 +265,10 @@ def calc_radii_ratio_from_rads(radius_sep_s, radius_sep_g):
     return radii_ratio_rad
 
 
-def calc_incl_from_radii_ratios_phase_incl(radii_ratio_lt, phase_orb_ext, phase_orb_int,
-                                           incl_init=np.deg2rad(85.0), show_plots=False):
+def calc_incl_from_radii_ratios_phase_incl(
+    radii_ratio_lt, phase_orb_ext, phase_orb_int, tol=1e-4, maxiter=10, show_plots=False):
     """Calculate inclination angle by minimizing difference of ratios of stellar radii as calulated
-    from light levels and lightcurve events (tangencies) for various values of inclination.
+    from light levels and light curve events (tangencies) for various values of inclination.
     
     Parameters
     ----------
@@ -274,89 +280,141 @@ def calc_incl_from_radii_ratios_phase_incl(radii_ratio_lt, phase_orb_ext, phase_
         Orbital phase angle at external tangencies (e.g. begin ingress, end egress). Unit is radians.
     phase_orb_int : float
         Orbital phase angle at internal tangencies (e.g. end ingress, begin egress). Unit is radians.
-    incl_init : {numpy.deg2rad(85.0)}, float, optional
-        Initial guess for inclination angle in radians for minimization procedure of differences in radii ratios.
-        Radii ratio difference has a local maximum for `incl` = numpy.deg2rad(90).
+    tol : {1e-4}, float, optional
+        Maximum tolerance for difference in radii ratios at a self-consistent solution for inclination,
+        i.e. at solved inclination:
+        abs('radii ratio from light levels' - 'radii ratio from eclipse events') < tol.
+    maxiter : {10}, int, optional
+        Maximum number of iterations to perform when solving for inclination.
+        For `tol = 1e-4`, the solution typically converges within 5 iterations.
     show_plots : {True}, bool, optional
-        Create and show diagnostic plots of difference in independent radii_ratio values vs inclination angle.
+        Create and show diagnostic plots of difference in independent radii ratio values vs inclination angle.
         Use to check solution in case initial guess for inclination angle caused convergence to wrong solution
         for inclination angle.
     
     Returns
     -------
     incl : float
         Orbital inclination. Angle between line of sight and the axis of the orbit. Unit is radians.
+        If solution does not exist or does not converge, `incl = numpy.nan`
     
     See Also
     --------
     calc_sep_proj_from_incl_phase, calc_radius_sep_s_from_sep, calc_radius_sep_g_from_sep, calc_radii_ratio_from_light
     
     Notes
     -----
-    Note: Method uses calc_radii_ratio_from_light, which may not be valid for stars in different stages of evolution
-    or for radii ratios > 10 (e.g. a binary system with main sequence star and a red giant)
-    Compares two independent ways of calculating radii_ratio in order to infer inclination angle. 
-    Equations from section 7.3 of [1]_.
+    - Method uses calc_radii_ratio_from_light, which may not be valid for stars in different stages of evolution
+        or for radii ratios > 10 (e.g. a binary system with main sequence star and a red giant)
+    - Compares two independent ways of calculating radii_ratio in order to infer inclination angle. 
+    - Equations from section 7.3 of [1]_.
     
     References
     ----------
     .. [1] Budding, 2007, Introduction to Astronomical Photometry
 
     """
     # Make radii_ratio_rad and all dependencies functions of inclination.
+    # Note: stdout and stderr are delayed if called within a loop in an IPython Notebook.
     sep_proj_ext = lambda incl: calc_sep_proj_from_incl_phase(incl=incl, phase_orb=phase_orb_ext)
     sep_proj_int = lambda incl: calc_sep_proj_from_incl_phase(incl=incl, phase_orb=phase_orb_int)
-    radii_sep = lambda incl: calc_radii_sep_from_seps(sep_proj_ext=sep_proj_ext(incl=incl),
-                                                      sep_proj_int=sep_proj_int(incl=incl))
+    radii_sep = lambda incl: calc_radii_sep_from_seps(
+        sep_proj_ext=sep_proj_ext(incl=incl), sep_proj_int=sep_proj_int(incl=incl))
     radii_ratio_rad = lambda incl: calc_radii_ratio_from_rads(*radii_sep(incl=incl))
-    # Minimize difference between independent radii_ratio values.
-    diff_radii_ratios = lambda incl: abs(radii_ratio_lt - radii_ratio_rad(incl=incl))
-    result = sci_opt.minimize(fun=diff_radii_ratios, x0=incl_init)
-    incl = result['x'][0]
-    # Create and show diagnostic plot.
+    diff_radii_ratios = lambda incl: radii_ratio_lt - radii_ratio_rad(incl=incl)
+    fmt_parameters =  \
+      ("    radii_ratio_lt = {rrl}\n" +
+       "    phase_orb_ext  = {poe}\n" +
+       "    phase_orb_int  = {poi}\n" +
+       "    tol            = {tol}").format(
+           rrl=radii_ratio_lt, poe=phase_orb_ext, poi=phase_orb_int, tol=tol)
+    # Minimize difference between independent radii_ratio values to within a tolerance.
+    # Note: A naive implentation of scipy.optimize.minimize will not find the solution
+    # for some parameters due to the non-differentiability of `abs(diff_radii_ratios)`
+    # Note: diff_radii_ratios is monotonically decreasing with a zero
+    # at the solution for self-consistent inclination:
+    # - For incl < incl_soln, diff_radii_ratios > 0.
+    # - For incl > incl_soln, diff_radii_ratios < 0.
+    # Find the solution to within `tol` by recursively zooming in where
+    # `diff_radii_ratios` changes sign.
+    incls = np.deg2rad(np.linspace(start=0.0, stop=90.0, num=10, endpoint=True))
+    diffs = np.asarray(map(diff_radii_ratios, incls))
+    if not (np.all(np.diff(diffs) <= 0.0)):
+        raise AssertionError(
+            "Program error. `diff_radii_ratios` must be monotonically decreasing.")
+    if (diffs[0] > 0.0) and (diffs[-1] < 0.0):
+        itol = diffs[0] - diffs[-1]
+        inum = 0
+        while (itol > tol) and (inum < maxiter):
+            if inum > 0:
+                incls = \
+                    np.linspace(
+                        start=incls[idx_diff_least_pos], stop=incls[idx_diff_least_neg],
+                        num=10, endpoint=True)
+                diffs = np.asarray(map(diff_radii_ratios, incls))
+            idx_diff_least_pos = len(diffs[diffs > 0.0]) - 1
+            idx_diff_least_neg = -1 * len(diffs[diffs < 0.0])
+            incl = incls[idx_diff_least_pos]
+            itol = abs(diff_radii_ratios(incl=incl))
+            inum += 1
+        # Check exit condition and solution.
+        if (itol > tol) and (inum >= maxiter):
+            incl = np.nan
+            warnings.warn(
+                ("\n" +
+                "    Difference in radii ratios did not converge to within tolerance.\n" +
+                "    Input parameters:\n" +
+                fmt_parameters))
+        if radii_ratio_rad(incl=incl) < 0.1:
+            warnings.warn(
+                ("\n" +
+                "    From eclipse timing events, ratio of smaller star's radius\n" +
+                "    to greater star's radius is < 0.1. The radii ratio as calculated\n" +
+                "    from light levels may not be valid (e.g. for a binary system with\n" +
+                "    a main sequence star and a red giant).\n" +
+                "    VALID:\n" +
+                "    radii_ratio_rad = radius_s / radius_g from eclipse timings = {rtime}\n" +
+                "    MAYBE INVALID:\n" +
+                "    radii_ratio_lt  = radius_s / radius_g from light levels = {rlt}").format(
+                    rtime=radii_ratio_rad(incl=incl),
+                    rlt=radii_ratio_lt))
+    else:
+        incl = np.nan
+        warnings.warn(
+            ("\n" +
+             "    Inclination does not yield self-consistent solution for model.\n" +
+             "    Input parameters cannot be fit by model:\n" +
+             fmt_parameters))
+    # Create and show diagnostic plots.
     if show_plots:
-        incls_out = np.deg2rad(np.linspace(0, 90, num=1000))
+        incls_out = np.deg2rad(np.linspace(start=0, stop=90, num=100))
         diff_radii_ratios_out = map(diff_radii_ratios, incls_out)
         plt.plot(np.rad2deg(incls_out), diff_radii_ratios_out)
+        plt.axhline(0.0, color='black', linestyle='--')
         plt.title("Difference between independent radii ratio values\n" +
                   "vs inclination angle (zoomed out view)")
-        plt.ylabel("abs(radii_ratio_lt - radii_ratio_rad(incl))")
-        plt.xlabel("inclination angle (degrees)")
+        xlabel = "inclination angle (degrees)"
+        ylabel = "radii ratio from light levels - radii ratio from eclipse events"
+        plt.xlabel(xlabel)
+        plt.ylabel(ylabel)
         plt.show()
-        incls_in = np.deg2rad(np.linspace(np.rad2deg(incl) - 1.0, np.rad2deg(incl) + 1.0, num=1000))
-        diff_radii_ratios_in = map(diff_radii_ratios, incls_in)
-        plt.plot(np.rad2deg(incls_in), diff_radii_ratios_in)
-        plt.title("Difference between independent radii ratio values\n" +
-                  "vs inclination angle (zoomed in view)")
-        plt.ylabel("abs(radii_ratio_lt - radii_ratio_rad(incl))")
-        plt.xlabel("inclination angle (degrees)")
-        plt.show()
-    # Check solutions.
-    if radii_ratio_rad(incl=incl) < 0.1:
-        print(("WARNING: From eclipse timing events, ratio of smaller star's radius\n" +
-               "    to greater star's radius is < 0.1. The radii ratio as calculated\n" +
-               "    from light levels may not be valid (e.g. for a binary system with\n" +
-               "    a main sequence star and a red giant).\n" +
-               "    VALID:\n" +
-               "    radii_ratio_rad = radius_s / radius_g from eclipse timings = {rtime}\n" +
-               "    MAYBE INVALID:\n" +
-               "    radii_ratio_lt  = radius_s / radius_g from light levels = {rlt}").format(rtime=radii_ratio_rad(incl=incl),
-                                                                                             rlt=radii_ratio_lt),
-              file=sys.stderr)
-    if diff_radii_ratios(incl) < 1e-3:
-        print("INFO: Inclination yields self-consistent solution for model.")
-    else:
-        # Note: Warning message is delayed if called within a loop in an IPython Notebook.
-        print(("WARNING: Inclination does not yield self-consistent solution for model.\n" +
-               "    Input parameters cannot be fit by model:\n" +
-               "    radii_ratio_lt   = {rrl}\n" +
-               "    phase_orb_ext    = {poe}\n" +
-               "    phase_orb_int    = {poi}\n" +
-               "    incl_init        = {ii}").format(rrl=radii_ratio_lt,
-                                                     poe=phase_orb_ext,
-                                                     poi=phase_orb_int,
-                                                     ii=incl_init),
-              file=sys.stderr)
+        if incl is not np.nan:
+            incls_in = \
+                np.deg2rad(
+                    np.linspace(
+                        start=np.rad2deg(incl) - 1.0,
+                        stop=min(np.rad2deg(incl) + 1.0, 90.0),
+                        num=100))
+            diff_radii_ratios_in = map(diff_radii_ratios, incls_in)
+            plt.plot(np.rad2deg(incls_in), diff_radii_ratios_in)
+            plt.axhline(0.0, color='black', linestyle='--')
+            plt.title("Difference between independent radii ratio values\n" +
+                      "vs inclination angle (zoomed in view)")
+            plt.xlabel(xlabel)
+            plt.ylabel(ylabel)
+            plt.show()
+        else:
+            warnings.warn("\nNo inclination solution found.")
     return incl
 
 

setup.py

@@ -1,5 +1,5 @@
 #!/usr/bin/env python
-"""Adapted from [1]_, [2]_.
+"""Setup script adapted from [1]_, [2]_.
 
 References
 ----------
@@ -27,8 +27,8 @@
     description='Estimate physical quantities of a binary star system from observed quantities.',
     long_description=long_description,
     url='http://binstarsolver.readthedocs.org',
-    author='White Dwarf Research Group, Astronomy Dept, UT Austin',
-    author_email='harrold@astro.as.utexas.edu',
+    author='Samuel Harrold; White Dwarf Research Group, Astronomy Dept, UT Austin',
+    author_email='samuel.harrold@gmail.com',
     license='MIT',
     classifiers=[
         # Project maturity values: 3 - Alpha; 4 - Beta; 5 - Production/Stable

tests/test_main.py

@@ -1,12 +1,15 @@
-"""Tests for binstarsolver/main.py.
-
-Notes
------
-Tests are executed using pytest.
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""Pytest style tests for binstarsolver/main.py.
 
 """
 
+
+# Import standard packages.
 from __future__ import absolute_import, division, print_function
+import sys
+sys.path.insert(0, '.') # Test the code in this repository.
+# Import local packages.
 import binstarsolver as bss