Merge pull request #924 from rzellem/issue922

Added index constraint on centroid fit near image bounds and update to light curve plot

exotic/api/elca.py

@@ -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')