Merge pull request #925 from rzellem/develop

v1.9.0: update to centroiding and lightcurve plotter

Author: rzellem

exotic/api/elca.py

@@ -123,7 +123,7 @@ def time_bin(time, flux, dt=1./(60*24)):
         if mask.sum() > 0:
             bflux[i] = np.nanmean(flux[mask])
             btime[i] = np.nanmean(time[mask])
-            bstds[i] = np.nanstd(flux[mask])/(1+mask.sum())**0.5
+            bstds[i] = np.nanstd(flux[mask])/(mask.sum()**0.5)
     zmask = (bflux==0) | (btime==0) | np.isnan(bflux) | np.isnan(btime)
     return btime[~zmask], bflux[~zmask], bstds[~zmask]
 
@@ -223,10 +223,12 @@ def lc2min_airmass(pars):
     def create_fit_variables(self):
         self.phase = get_phase(self.time, self.parameters['per'], self.parameters['tmid'])
         self.transit = transit(self.time, self.parameters)
+        self.time_upsample = np.linspace(min(self.time), max(self.time),1000)
+        self.transit_upsample = transit(self.time_upsample, self.parameters)
+        self.phase_upsample = get_phase(self.time_upsample, self.parameters['per'], self.parameters['tmid'])
         if self.mode == "ns":
             self.parameters['a1'], self.errors['a1'] = mc_a1(self.parameters.get('a2',0), self.errors.get('a2',1e-6),
                                                              self.transit, self.airmass, self.data)
-
         if np.ndim(self.airmass) == 2:
             detrended = self.data/self.transit
             self.wf = weightedflux(detrended, self.gw, self.nearest)
@@ -409,7 +411,9 @@ def plot_bestfit(self, title="", bin_dt=30./(60*24), zoom=False, phase=True):
             si = np.argsort(self.phase)
             bt2, bf2, bs = time_bin(self.phase[si]*self.parameters['per'], self.detrended[si], bin_dt)
             axs[0].errorbar(bt2/self.parameters['per'],bf2,yerr=bs,alpha=1,zorder=2,color='blue',ls='none',marker='s')
-            axs[0].plot(self.phase[si], self.transit[si], 'r-', zorder=3, label=lclabel)
+            #axs[0].plot(self.phase[si], self.transit[si], 'r-', zorder=3, label=lclabel)
+            sii = np.argsort(self.phase_upsample)
+            axs[0].plot(self.phase_upsample[sii], self.transit_upsample[sii], 'r-', zorder=3, label=lclabel)
             axs[0].set_xlim([min(self.phase), max(self.phase)])
             axs[0].set_xlabel("Phase ", fontsize=14)
         else:
@@ -421,8 +425,9 @@ def plot_bestfit(self, title="", bin_dt=30./(60*24), zoom=False, phase=True):
 
             bt, bf, bs = time_bin(self.time, self.detrended, bin_dt)
             si = np.argsort(self.time)
+            sii = np.argsort(self.time_upsample)
             axs[0].errorbar(bt,bf,yerr=bs,alpha=1,zorder=2,color='blue',ls='none',marker='s')
-            axs[0].plot(self.time[si], self.transit[si], 'r-', zorder=3, label=lclabel)
+            axs[0].plot(self.time_upsample[sii], self.transit_upsample[sii], 'r-', zorder=3, label=lclabel)
             axs[0].set_xlim([min(self.time), max(self.time)])
             axs[0].set_xlabel("Time [day]", fontsize=14)
 

exotic/exotic.py

@@ -900,28 +900,34 @@ def gaussian_psf(x, y, x0, y0, a, sigx, sigy, rot, b):
     return a * gausx * gausy + b
 
 
-def mesh_box(pos, box):
+def mesh_box(pos, box, maxx=0, maxy=0):
     pos = [int(np.round(pos[0])), int(np.round(pos[1]))]
-    x = np.arange(pos[0] - box, pos[0] + box + 1)
-    y = np.arange(pos[1] - box, pos[1] + box + 1)
+    if maxx:
+        x = np.arange(max(0,pos[0] - box), min(maxx, pos[0] + box + 1))
+    else:
+        x = np.arange(max(0,pos[0] - box), pos[0] + box + 1)
+    if maxy:
+        y = np.arange(max(0,pos[1] - box), min(maxy, pos[1] + box + 1))
+    else:
+        y = np.arange(max(0,pos[1] - box), pos[1] + box + 1)
     xv, yv = np.meshgrid(x, y)
     return xv.astype(int), yv.astype(int)
 
 
 # Method fits a 2D gaussian function that matches the star_psf to the star image and returns its pixel coordinates
 def fit_centroid(data, pos, psf_function=gaussian_psf, box=10):
     # get sub field in image
-    xv, yv = mesh_box(pos, box)
-
-    init = [np.nanmax(data[yv, xv]) - np.nanmin(data[yv, xv]), 1, 1, 0, np.nanmin(data[yv, xv])]
-
+    xv, yv = mesh_box(pos, box, maxx=data.shape[1], maxy=data.shape[0])
+    subarray = data[yv, xv]
+    init = [np.nanmax(subarray) - np.nanmin(subarray), 1, 1, 0, np.nanmin(subarray)]
+    
     # lower bound: [xc, yc, amp, sigx, sigy, rotation,  bg]
-    lo = [pos[0] - box * 0.5, pos[1] - box * 0.5, 0, 0.5, 0.5, -np.pi / 4, np.nanmin(data[yv, xv]) - 1]
-    up = [pos[0] + box * 0.5, pos[1] + box * 0.5, 1e7, 20, 20, np.pi / 4, np.nanmax(data[yv, xv]) + 1]
+    lo = [pos[0] - box * 0.5, pos[1] - box * 0.5, 0, 0.5, 0.5, -np.pi / 4, np.nanmin(subarray) - 1]
+    up = [pos[0] + box * 0.5, pos[1] + box * 0.5, 1e7, 20, 20, np.pi / 4, np.nanmax(subarray) + 1]
 
     def fcn2min(pars):
         model = psf_function(xv, yv, *pars)
-        return (data[yv, xv] - model).flatten()
+        return (subarray - model).flatten()
 
     try:
         res = least_squares(fcn2min, x0=[*pos, *init], bounds=[lo, up], jac='3-point', xtol=None, method='trf')

exotic/version.py

@@ -1 +1 @@
-__version__ = "1.8.4"
+__version__ = "1.9.0"

requirements.txt

@@ -7,7 +7,7 @@ astroalign~=2.4
 astropy>=4.3
 astroquery~=0.4
 barycorrpy~=0.4
-cython~=0.29.24 ; platform_system != "Windows"
+cython~=0.29.25 ; platform_system != "Windows"
 dynesty~=1.1 ; platform_system == "Windows"
 holoviews~=1.14
 LDTk~=1.7