Switch to normalized oot, and check_call stsp runs

rms/lightcurve.py

@@ -307,7 +307,7 @@ def mask_out_of_transit(self, star, oot_duration_fraction=0.25,
         near_transit = ((phased < star.planet.duration*(0.5 + oot_duration_fraction)) |
                         (phased > star.planet.per - star.planet.duration*(0.5 + oot_duration_fraction)))
         if flip:
-            near_transit = ~near_transit
+            near_transit = np.logical_not(near_transit)
         sort_by_time = np.argsort(self.times[near_transit].jd)
         return dict(times=self.times[near_transit][sort_by_time],
                     fluxes=self.fluxes[near_transit][sort_by_time],

rms/planet.py

@@ -26,6 +26,29 @@
     "inc": 89.34708509841991
 }
 
+hat7_params = {
+    "a": 4.1502,
+    "fp": None,
+    "b": 0.50021254073449306,
+    "ecc": 0.0,
+    "inc_stellar": 0,  # Lund 2014
+    "rho_star": 0.18875,  # Lund 2014
+    "limb_dark": "quadratic",
+    "per": 2.204737,
+    "per_rot": 13.0, # Lund 2014
+    "lam": 155,  # Albrecht et al 2012
+    "t_secondary": None,
+    "t0": 2454954.357463,
+    "w": 0.0,
+    "duration": 0.12813004872191647,
+    "rp": 0.058330305324663184,
+    "u": [
+        0.3525,
+        0.168
+    ],
+    "inc": 83.111
+}
+
 
 class Planet(object):
     """
@@ -95,3 +118,7 @@ def __init__(self, per=None, inc=None, a=None, t0=None, rp=None, ecc=None,
     @classmethod
     def from_hat11(cls):
         return cls(**hat11_params)
+
+    @classmethod
+    def from_hat7(cls):
+        return cls(**hat7_params)

rms/spot.py

@@ -34,8 +34,8 @@ def at_random_position(cls, radius):
         return cls(lat, lon, radius)
 
     @classmethod
-    def from_sunspot_distribution(cls):
-        lat = draw_random_sunspot_latitudes(n=1)[0]
+    def from_sunspot_distribution(cls, mean_latitude=15):
+        lat = draw_random_sunspot_latitudes(n=1, mean_latitude=mean_latitude)[0]
         lon = 2*np.pi * np.random.rand() * u.rad
         radius = draw_random_sunspot_radii(n=1)[0]
         return cls(lat, lon, radius)
@@ -143,7 +143,7 @@ def generate_random_coord(n=1):
     return result
 
 
-def sunspot_distribution(latitude):
+def sunspot_distribution(latitude, mean_latitude=15):
     """
     Approximate un-normalized probability distribution of sunspots at
     ``latitude`` near activity maximum on the Sun.
@@ -153,15 +153,18 @@ def sunspot_distribution(latitude):
     latitude : `~numpy.ndarray`
         Latitude
 
+    mean_latitude : float
+        Mean active latitude in degrees
+
     Returns
     -------
     p : `~numpy.ndarray
         Probability (un-normalized)
     """
-    return np.exp(-0.5 * (abs(latitude) - 15)**2 / 6**2)
+    return np.exp(-0.5 * (abs(latitude) - mean_latitude)**2 / 6**2)
 
 
-def sunspot_latitude_inverse_transform(x):
+def sunspot_latitude_inverse_transform(x, mean_latitude=15):
     """
     Use inverse transform sampling to randomly draw spot latitudes from the
     sunspot latitude distribution, for a uniform random variate ``x`` on the
@@ -177,13 +180,13 @@ def sunspot_latitude_inverse_transform(x):
     lat : `~astropy.units.Quantity`
         Latitude of a sunspot drawn from the sunspot latitude distribution.
     """
-    lats = np.linspace(-60, 60, 1000)
-    prob = np.cumsum(sunspot_distribution(lats))
+    lats = np.linspace(-85, 85, 1000)
+    prob = np.cumsum(sunspot_distribution(lats, mean_latitude=mean_latitude))
     prob /= np.max(prob)
     return np.interp(x, prob, lats) * u.deg
 
 
-def draw_random_sunspot_latitudes(n):
+def draw_random_sunspot_latitudes(n, mean_latitude=15):
     """
     Draw one or more random samples from the sunspot latitude distribution.
 
@@ -197,7 +200,8 @@ def draw_random_sunspot_latitudes(n):
     lat : `~astropy.units.Quantity`
         Latitude of a sunspot drawn from the sunspot latitude distribution.
     """
-    return sunspot_latitude_inverse_transform(np.random.rand(n))
+    return sunspot_latitude_inverse_transform(np.random.rand(n),
+                                              mean_latitude=mean_latitude)
 
 
 def sunspot_umbral_area_distribution(log_area_uhem):

rms/stsp.py

@@ -9,6 +9,7 @@
 
 import numpy as np
 from astropy.io import ascii
+from astropy.io.ascii import InconsistentTableError
 
 from .lightcurve import LightCurve
 from .exceptions import OverlappingSpotsWarning, STSPMemoryWarning
@@ -47,7 +48,7 @@
 {start_time:2.10f}		; start time to start fitting the light curve
 {lc_duration:2.10f}		; duration of light curve to fit (days)
 {real_max:2.10f}		; real maximum of light curve data (corrected for noise), 0 -> use downfrommax
-0						; is light curve flattened (to zero) outside of transits?
+1						; is light curve flattened (to zero) outside of transits?
 #ACTION
 l						; l= generate light curve from parameters
 {spot_params}
@@ -67,17 +68,17 @@ def quadratic_to_nonlinear(u1, u2):
     return a1, a2, a3, a4
 
 
-def _spot_obj_to_params(spot):
+def _spot_obj_to_params(spot, quiet=False):
 
     if hasattr(spot, '__len__'):
-        validated_spot_list = find_overlapping_spots(spot)
+        validated_spot_list = find_overlapping_spots(spot, quiet=quiet)
         return np.concatenate([[s.r, s.theta, s.phi]
                                for s in validated_spot_list])
     else:
         return np.array([spot.r, spot.theta, spot.phi])
 
 
-def find_overlapping_spots(spot_list, tolerance=1.01):
+def find_overlapping_spots(spot_list, tolerance=1.05, quiet=False):
     """
     Find overlapping spots in a list of spot objects.
 
@@ -110,7 +111,7 @@ def find_overlapping_spots(spot_list, tolerance=1.01):
                         if i in save_these_spot_indices]
     toss_these_spots = [spot for i, spot in enumerate(spot_list)
                         if i in toss_these_spot_indices]
-    if len(spots_with_overlap) > 0:
+    if len(spots_with_overlap) > 0 and not quiet:
         warning_message = ('Some spots were overlapping. Tossing one of the two'
                            ' overlapping spots. \n\nSpots tossed:\n\n' +
                            '\n'.join(map(str, toss_these_spots)))
@@ -170,7 +171,8 @@ class STSP(object):
     """
     Context manager for working with STSP
     """
-    def __init__(self, times, star, spot, outdir=None, keep_dir=False):
+    def __init__(self, times, star, spot, outdir=None, keep_dir=False,
+                 quiet=False):
         """
         Parameters
         ----------
@@ -185,7 +187,7 @@ def __init__(self, times, star, spot, outdir=None, keep_dir=False):
         """
         self.times = times
         self.star = star
-        self.spot_params = _spot_obj_to_params(spot)
+        self.spot_params = _spot_obj_to_params(spot, quiet=quiet)
         self.spot_contrast = self.star.spot_contrast
 
         current_time = datetime.datetime.now().strftime("%Y%m%d_%H%M%S_%f")
@@ -281,22 +283,26 @@ def generate_lightcurve(self, n_ld_rings=40, stsp_exec=None):
             in_file.write(in_file_text)
 
         try:
-            stdout = subprocess.check_output([stsp_exec, 'test.in'], cwd=self.outdir)
+            stdout = subprocess.check_call([stsp_exec, 'test.in'], cwd=self.outdir)
         except subprocess.CalledProcessError as err:
-            pass  # print("Failed. Error:", err.output, err.stderr, err.stdout)
+            print("Failed. Error:", err.output, err.stderr, err.stdout)
 
-        #time.sleep(0.01)
         # Read the outputs
         if os.stat(os.path.join(self.outdir, 'test_lcout.txt')).st_size == 0:
             stsp_times = self.times.jd
-            stsp_fluxes = np.ones_like(self.times)
-            stsp_flag = 0 * np.ones_like(self.times)
+            stsp_fluxes = np.ones(len(self.times))
+            stsp_flag = 0 * np.ones(len(self.times))
 
         else:
-            tbl = ascii.read(os.path.join(self.outdir, 'test_lcout.txt'),
-                             format='fast_no_header')
-            stsp_times, stsp_fluxes, stsp_flag = tbl['col1'], tbl['col4'], tbl['col5']
-        return LightCurve(times=stsp_times, fluxes=stsp_fluxes.data, quarters=stsp_flag)
+            try:
+                tbl = ascii.read(os.path.join(self.outdir, 'test_lcout.txt'),
+                                 format='fast_no_header')
+                stsp_times, stsp_fluxes, stsp_flag = tbl['col1'], tbl['col4'].data, tbl['col5']
+            except InconsistentTableError:
+                stsp_times = self.times.jd
+                stsp_fluxes = np.ones(len(self.times)) * np.nan
+                stsp_flag = 0 * np.ones(len(self.times))
+        return LightCurve(times=stsp_times, fluxes=stsp_fluxes, quarters=stsp_flag)
 
 
 def _spot_params_to_string(spot_params):