StarDICE

config/stardice.ini

@@ -12,17 +12,17 @@ SPECTRACTOR_COMPUTE_ROTATION_ANGLE = disperser
 # parameter to prevent from background subtraction
 SPECTRACTOR_BACKGROUND_SUBTRACTION = True
 # deconvolve spectrogram with simple 2D PSF analysis: False, True
-SPECTRACTOR_DECONVOLUTION_PSF2D = False
+SPECTRACTOR_DECONVOLUTION_PSF2D = True
 # deconvolve spectrogram with full forward model: False, True
-SPECTRACTOR_DECONVOLUTION_FFM = False
+SPECTRACTOR_DECONVOLUTION_FFM = True
 # value of sigma clip parameter for the spectractor deconvolution process PSF2D and FFM
 SPECTRACTOR_DECONVOLUTION_SIGMA_CLIP = 100
 # library to compute atmospheric transmission: none, libradtran, getobsatmo
 SPECTRACTOR_ATMOSPHERE_SIM = getobsatmo
 
 [instrument]
 # instrument name
-OBS_NAME = STARDICE
+OBS_NAME = OHP 
 # observatory altitude in km
 OBS_ALTITUDE = 0.650
 # observatory latitude
@@ -108,11 +108,11 @@ PIXWIDTH_BOXSIZE = 40
 # the PSF model: Gauss, Moffat or MoffatGauss
 PSF_TYPE = Moffat
 # the order of the polynomials to model wavelength dependence of the PSF shape parameters
-PSF_POLY_ORDER = 2
+PSF_POLY_ORDER = 4
 # regularisation parameter for the chisq minimisation to extract the spectrum
 PSF_FIT_REG_PARAM = 0.1
 # step size in pixels for the first transverse PSF1D fit
-PSF_PIXEL_STEP_TRANSVERSE_FIT = 10
+PSF_PIXEL_STEP_TRANSVERSE_FIT = 1
 # PSF is not evaluated outside a region larger than max(PIXWIDTH_SIGNAL, PSF_FWHM_CLIP*fwhm) pixels
 PSF_FWHM_CLIP = 2
 

runFitter.py

@@ -18,7 +18,7 @@
     parser.add_argument("-o", "--output_directory", dest="output_directory", default="outputs/",
                         help="Write results in given output directory (default: ./outputs/).")
     parser.add_argument("-c", "--config", dest="config", default="",
-                        help="config file to be given for spectra extracted with Spectractor<2.4. (default: ''.")
+                        help="config file to be given for spectra extracted with Spectractor<2.4. (default: ''.)")
     args = parser.parse_args()
 
     parameters.VERBOSE = args.verbose

spectractor/astrometry.py

@@ -323,7 +323,7 @@ def __init__(self, image, wcs_file_name="", gaia_file_name="", output_directory=
         self.gaia_radec_positions_after_pm = None
         if os.path.isfile(self.gaia_file_name):
             self.my_logger.info(f"\n\tLoad Gaia catalog from {self.gaia_file_name}.")
-            self.gaia_catalog = ascii.read(self.gaia_file_name, format="ecsv")
+            self.gaia_catalog = Table(ascii.read(self.gaia_file_name, format="ecsv"), masked=True)
             self.gaia_radec_positions_after_pm = get_gaia_coords_after_proper_motion(self.gaia_catalog, self.image.date_obs)
         self.sources = None
         self.sources_radec_positions = None
@@ -1095,9 +1095,16 @@ def run_simple_astrometry(self, extent=None, sources=None):
                 self.write_sources()
                 solve_field_input = self.sources_file_name
             elif self.source_extractor == "astrometrynet":
+                # use data without spectrums if it exists
+                if "data_wo_spectrums" in self.image.__dict__.keys():
+                    if getattr(self.image, "data_wo_spectrums") is not None:
+                        data = getattr(self.image, "data_wo_spectrums").copy()
+                    else:
+                        data = self.image.data
+
                 self.my_logger.info(f"\n\tSource extraction directly with solve-field.")
                 # must write a temporary image file with Spectractor flips and rotations
-                fits.writeto(tmp_image_file_name, self.image.data, header=self.image.header, overwrite=True)
+                fits.writeto(tmp_image_file_name, data, header=self.image.header, overwrite=True)
                 solve_field_input = tmp_image_file_name
             else:
                 raise ValueError(f"Got {self.source_extractor=}. Must be either 'iraf' or 'astrometrynet' "
@@ -1121,7 +1128,7 @@ def run_simple_astrometry(self, extent=None, sources=None):
         try:
             subprocess.run(command, shell=True, check=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, encoding='ascii')
         except subprocess.CalledProcessError as e:  # pragma: nocover
-            self.my_logger.warning(f"\n\tAstrometry command:\n{command}")
+            self.my_logger.info(f"\n\tAstrometry command:\n{command}")
             self.my_logger.error(f"\n\t{e.stderr}")
             sys.exit()
         if os.path.isfile(tmp_image_file_name):
@@ -1218,7 +1225,7 @@ def run_gaia_astrometry(self, min_range=3 * u.arcsec, max_range=5 * u.arcmin, ma
 
         # update coordinates with proper motion data
         self.my_logger.info(f"\n\tUpdate object coordinates with proper motion at time={self.image.date_obs}.")
-        self.image.target_radec_position_after_pm = self.image.target.get_radec_position_after_pm(self.image.date_obs)
+        self.image.target_radec_position_after_pm = self.image.target.get_radec_position_after_pm(self.image.target.simbad_table, date_obs=self.image.date_obs)
         self.gaia_radec_positions_after_pm = get_gaia_coords_after_proper_motion(self.gaia_catalog, self.image.date_obs)
         if parameters.DEBUG:
             self.plot_sources_and_gaia_catalog(sources=self.sources, gaia_coord=self.gaia_radec_positions_after_pm,
@@ -1407,8 +1414,8 @@ def run_full_astrometry(self, extent=None, maxiter=20, min_range=3 * u.arcsec, m
                 gaia_residuals_quad_sum = np.sum(np.sqrt(np.sum(gaia_residuals ** 2, axis=1)))
                 if parameters.DEBUG:
                     self.plot_sources_and_gaia_catalog(sources=self.sources, gaia_coord=self.gaia_matches, margin=20,
-                                                       quad=self.quad_stars_pixel_positions,
-                                                       label=self.image.target.label)
+                                                        quad=self.quad_stars_pixel_positions,
+                                                        label=self.image.target.label)
                     self.plot_astrometry_shifts(vmax=3)
                     self.plot_quad_stars()
                 target_x, target_y = self.get_target_pixel_position()
@@ -1419,7 +1426,7 @@ def run_full_astrometry(self, extent=None, maxiter=20, min_range=3 * u.arcsec, m
                                        f"Stop the loop here and look for best solution.")
                 k -= 1
                 break
-        t.pprint_all()
+        #t.pprint_all()
         if len(t) == 0:
             raise IndexError(f"Astrometry has failed at every iteration, empty table {t=}.")
         best_iter = int(np.argmin(t["gaia_residuals_quad_sum"]))

spectractor/extractor/background.py

@@ -275,7 +275,8 @@ def extract_spectrogram_background_sextractor(data, err, ws=(20, 30), mask_signa
     bkg = Background2D(data, (parameters.PIXWIDTH_BOXSIZE, parameters.PIXWIDTH_BOXSIZE),
                        filter_size=(filter_size, filter_size),
                        sigma_clip=sigma_clip, bkg_estimator=bkg_estimator,
-                       mask=mask)
+                       mask=mask,
+                       exclude_percentile=50)
     # reset at zero the edges
     bkg.background[data == 0] = 0
     bgd_model_func_interp = RegularGridInterpolator((np.arange(Nx), np.arange(Ny)), bkg.background.T, method='linear',

spectractor/extractor/chromaticpsf.py

@@ -1340,11 +1340,14 @@ def from_profile_params_to_shape_params(self, profile_params):
         # pixels.shape = (2, Ny, Nx): self.pixels[1<-y, :, 0<-first pixel value column]
         # TODO: account for rotation ad projection effects is PSF is not round
         pixel_eval = np.arange(self.pixels[1, 0, 0], self.pixels[1, -1, 0], 0.5, dtype=np.float32)
+        center = (np.max(pixel_eval) - np.min(pixel_eval)) / 2
         for ix, x in enumerate(pixel_x):
-            p = profile_params[x, :]
-            # compute FWHM transverse to dispersion axis (assuming revolution symmetry of the PSF)
+            # centering the PSF on the pixel grid to compute FWHM
+            p = np.copy(profile_params[x, :])
+            p[2] = center
+            # compute FWHM transverse in 1D (assuming revolution symmetry of the PSF)
             out = self.psf.evaluate(pixel_eval, values=p)
-            fwhms[ix] = compute_fwhm(pixel_eval, out, center=p[2], minimum=0, epsilon=1e-2)
+            fwhms[ix] = compute_fwhm(pixel_eval, out, center=center, minimum=0, epsilon=1e-2)
         # clean fwhm bad points
         mask = np.logical_and(fwhms > 1, fwhms < self.Ny // 2)  # more than 1 pixel or less than window
         self.table['fwhm'] = interp1d(pixel_x[mask], fwhms[mask], kind="linear",
@@ -1854,7 +1857,10 @@ def fit_chromatic_psf(self, data, mask=None, bgd_model_func=None, data_errors=No
             raise ValueError(f"mode argument must be '1D' or '2D'. Got {mode=}.")
         if amplitude_priors_method == "psf1d":
             w_reg = RegFitWorkspace(w, opt_reg=parameters.PSF_FIT_REG_PARAM, verbose=verbose)
-            w_reg.run_regularisation(Ndof=w.trace_r)
+            weighted_mean_fwhm = np.average(self.table['fwhm'], weights=self.table['amplitude'])
+            self.my_logger.info(f"\n\tMean FWHM: {weighted_mean_fwhm} pixels (weighted with spectrum amplitude)"
+                                f"\n\tExpected Ndof: {self.Nx / weighted_mean_fwhm}")
+            w_reg.run_regularisation(Ndof=self.Nx / weighted_mean_fwhm)
             w.reg = np.copy(w_reg.opt_reg)
             w.trace_r = np.trace(w_reg.resolution)
             self.opt_reg = w_reg.opt_reg