Added option to save photometry of fitted BG stars

pipe/pipe_param.py

@@ -142,7 +142,7 @@ def __init__(self, name, visit, version=None, outdir=None, sa_range=None,
         # Save switches
         self.save_mask_cube = True   # Save mask used to filter out bad
                                      # data (as fits file)
-        self.save_bg_mask_cube = True  # If BG star mask defined, save is as fits
+        self.save_bg_mask_cube = True  # If BG star mask defined, save it as fits
         self.save_resid_cube = True  # Save cube of residuals (as fits file)
         self.save_bg_cube = False    # Save cube of residuals with bg stars (as fits file)
         self.save_bg_models = False  # Save model of background, incl stars, smearing, static
@@ -151,7 +151,8 @@ def __init__(self, name, visit, version=None, outdir=None, sa_range=None,
         self.save_psf_list = True    # Save list of filenames of PSFs used
         self.save_motion_mat = False # Save fitted motion blur matrix
         self.save_noise_cubes = False # Save estimated noise (raw/PSF/empiric) as fits cubes
-        self.save_gain = False       # Save  estimated gain table (with columns MJD, gain)
+        self.save_gain = False       # Save estimated gain table (with columns MJD, gain)
+        self.save_bg_star_phot = True # If BG stars are fitted, save their photometry
         self.save_astrometry = False # For binaries, saves text file with separation
 
         # Extraction parameters

pipe/psf_phot.py

@@ -17,15 +17,13 @@
 import os
 import pickle
 import numpy as np
-from scipy.ndimage import shift
 from astropy.io import fits
 
 from .analyse import psf_phot_cube
 from .cent import (
     flux as cent_flux
 )
 from .multi_cent import (
-    psf as multi_cent_psf,
     deconvolve as multi_cent_deconvolve,
     binary_deconvolve as multi_cent_binary_deconvolve,
     binary_psf as multi_cent_binary_psf,
@@ -42,7 +40,7 @@
 from .read import (
     imagette_offset, raw_datacube as read_raw_datacube, attitude, gain as read_gain,
     bias_ron_adu, thermFront_2, mjd2bjd, nonlinear, flatfield, starcat,
-    save_eigen_fits, save_binary_eigen_fits, sub_image_indices,
+    save_eigen_fits, save_bg_star_phot_fits, save_binary_eigen_fits, sub_image_indices,
     dark as read_dark, bad as read_bad, PSFs as load_PSFs, read_psf_filenames,
     save_psf_filenames
 )
@@ -1813,6 +1811,18 @@ def refine_star_bg_sa(self):
                                   krn_scl=self.pps.motion_step,
                                   krn_rad=self.pps.motion_nsteps,
                                   nthreads=self.pps.nthreads)
+        if len(starids) > 0 and self.pps.save_bg_star_phot:
+            gaiaID = [self.starcat.gaiaID[starid] for starid in starids]
+            fluxes = np.zeros((len(self.sa_workcat), len(starids)))
+            for n in range(len(self.sa_workcat)):
+                fluxes[n] = self.sa_workcat[n].fscale[starids]
+            filename = os.path.join(self.pps.outdir, 'bg_star_photometry_sa.fits')
+            save_bg_star_phot_fits(filename,
+                                self.sa_att[:,0],
+                                self.mjd2bjd(self.sa_att[:,0]),
+                                fluxes,
+                                gaiaID,
+                                self.sa_hdr)
 
 
     def refine_star_bg_im(self):
@@ -1837,9 +1847,18 @@ def refine_star_bg_im(self):
                                   krn_scl=self.pps.motion_step,
                                   krn_rad=self.pps.motion_nsteps,
                                   nthreads=self.pps.nthreads)
-
-
-
+        if len(starids) > 0 and self.pps.save_bg_star_phot:
+            gaiaID = [self.starcat.gaiaID[starid] for starid in starids]
+            fluxes = np.zeros((len(self.im_workcat), len(starids)))
+            for n in range(len(self.im_workcat)):
+                fluxes[n] = self.im_workcat[n].fscale[starids]
+            filename = os.path.join(self.pps.outdir, 'bg_star_photometry_im.fits')
+            save_bg_star_phot_fits(filename,
+                                self.im_att[:,0],
+                                self.mjd2bjd(self.im_att[:,0]),
+                                fluxes,
+                                gaiaID,
+                                self.im_hdr)
 
 
     def has_source(self, data, mask=None, clip=5, niter=10):

pipe/read.py

@@ -340,6 +340,25 @@ def imagette_offset(filename, frame_range=None):
     # raise Exception('[imagette_offset] Error: {:s} not found'.format(filename))
 
 
+
+def save_bg_star_phot_fits(filename, t, bjd, fluxes, gaia_IDs, header):
+    """Save lightcurve data of fitted background stars, as defined by
+    arguments, to fits table in binary format. The stars are identified by 
+    their Gaia ID numbers in the comment for respective fits-column
+    """
+    c = []
+    c.append(fits.Column(name='MJD_TIME', format='D', unit='day', array=t))
+    c.append(fits.Column(name='BJD_TIME', format='D', unit='day', array=bjd))
+    for n in range(fluxes.shape[1]):
+        c.append(fits.Column(name='f{:d}'.format(n), format='D', unit='electrons',
+                             array=fluxes[:, n]))
+    tab = fits.BinTableHDU.from_columns(c, header=header)
+    for n in range(fluxes.shape[1]):
+        key = f'TTYPE{n+3}'
+        tab.header[key] = (tab.header[key], f'Gaia {gaia_IDs[n]}')
+    tab.writeto(filename, overwrite=True, checksum=True)
+
+
 def save_eigen_fits(filename, t, bjd, sc, err, bg, roll, xc, yc, flag,
                     w, thermFront_2, header):
     """Save lightcurve data as defined by arguments to fits table in binary

pipe/syntstar.py

@@ -45,7 +45,7 @@ def __init__(self, starcatfile, psf_lib, maxrad=None,
         self.psf_mod_Teff = np.array([3000.0, 4000.0, 5000.0, 6000.0, 8000.0, 10000.0])
         self.default_psf_id = 2
         self.psf_lib = psf_lib
-        self.xpos, self.ypos, self.fscale, self.Teff = \
+        self.xpos, self.ypos, self.fscale, self.Teff, self.gaiaID = \
             self.read_starcat(starcatfile, maxrad=maxrad, fscalemin=fscalemin)
         self.catsize = len(self.fscale)
         self.psf_ids, self.psfs = self.assign_psf()
@@ -94,9 +94,10 @@ def read_starcat(self, starcatfile, maxrad=None, fscalemin=1e-5):
                    np.cos(np.deg2rad(cat['DEC'][0])) * 3600.0 / self.pxl_scl)
             dy = ((cat['DEC'][:N]-cat['DEC'][0]) * 3600.0 / self.pxl_scl)
             Teff = cat['T_EFF'][:N]
+            gaiaID = cat['ID'][:N]
 
             sel = fscale > fscalemin
-            return dx[sel], dy[sel], fscale[sel], Teff[sel]
+            return dx[sel], dy[sel], fscale[sel], Teff[sel], gaiaID[sel]
 
 
     def rotate_cat(self, rolldeg, maxrad=None):