coltess

12@dataclass
13class StarData:
14    """
15    This class stores the star's identifying information (name, coordinates,
16    Gaia ID) and, once analysis is performed, its light-curve data (times
17    and fluxes).
18
19    The object is intentionally mutable: analysis functions are expected to
20    populate photometry fields after creation.
21    """
22
23    # Identity / catalog info
24    name: str
25    ra: float
26    dec: float
27    gaia_id: Optional[str] = None
28
29    # Photometry results (optional)
30    times: Optional[np.ndarray] = field(default=None, repr=False)
31    fluxes: Optional[np.ndarray] = field(default=None, repr=False)
32    flux_errors: Optional[np.ndarray] = field(default=None, repr=False)
33
34    # Metadata
35    sector: Optional[int] = None
36
37    def has_photometry(self) -> bool:
38        """Return True if the star already has light-curve data loaded."""
39        return self.times is not None and self.fluxes is not None

 62class TessPhotometry:
 63    """
 64    Aperture photometry pipeline for TESS Full Frame Images (FFIs).
 65
 66    This class provides utilities to load a Gaia-based source catalog,
 67    identify catalog sources falling within a TESS image footprint,
 68    and perform aperture photometry with local background subtraction.
 69
 70    The class is designed to be safe for batch and parallel processing
 71    of large TESS datasets.
 72    """
 73
 74    def __init__(
 75        self,
 76        aperture_radius: int = 10,
 77        annulus_inner: int = 12,
 78        annulus_outer: int = 14,
 79        zeropoint: float = 20.4402281476,
 80    ):
 81        """
 82        Initialize the photometry configuration.
 83
 84        Parameters
 85        ----------
 86        aperture_radius : int, optional
 87            Radius of the circular photometric aperture in pixels.
 88        annulus_inner : int, optional
 89            Inner radius of the background annulus in pixels.
 90        annulus_outer : int, optional
 91            Outer radius of the background annulus in pixels.
 92        zeropoint : float, optional
 93            Photometric zeropoint used to convert fluxes to instrumental
 94            magnitudes.
 95
 96        Notes
 97        -----
 98        The gain value (e-/ADU) is fixed to the nominal TESS FFI gain.
 99        """
100        self.aperture_radius = aperture_radius
101        self.annulus_inner = annulus_inner
102        self.annulus_outer = annulus_outer
103        self.zeropoint = zeropoint
104        self.epadu = 5.22  # TESS gain (e-/ADU); overridden for e-/s data
105
106    def load_catalog(self, catalog_file: str):
107        """
108        Load a Gaia source catalog from a CSV file.
109
110        Parameters
111        ----------
112        catalog_file : str
113            Path to a CSV file containing at least `ra`, `dec`, and `source_id`
114            columns in degrees (ICRS).
115
116        Returns
117        -------
118        list of tuple
119            List of `(ra, dec, source_id)` entries for all catalog sources.
120        """
121        return _load_catalog_cached(catalog_file)
122
123    def process_fits(self, fits_path: str, catalog: List[tuple]) -> Optional[Table]:
124        """
125        Process a single TESS FITS image and perform photometry on catalog sources.
126
127        Parameters
128        ----------
129        fits_path : str
130            Path to the TESS FFI FITS file.
131        catalog : list of tuple
132            Catalog entries as `(ra, dec, source_id)`.
133
134        Returns
135        -------
136        astropy.table.Table or None
137            Table containing photometric measurements and metadata for all
138            detected catalog sources in the frame, or None if no sources fall
139            within the image footprint.
140        """
141        with fits.open(fits_path) as hdul:
142            image = hdul[1].data
143            header = hdul[1].header
144            wcs = WCS(header)
145
146            # Determine effective gain for Poisson noise term
147            bunit = header.get("BUNIT", "").strip().lower()
148            if "e-/s" in bunit or "e- /s" in bunit:
149                self.epadu = header.get("EXPTIME", header.get("TINT", 200.0))
150            else:
151                self.epadu = 5.22
152
153        # Filter objects in frame
154        objects_in_frame = []
155        for ra, dec, obj_id in catalog:
156            try:
157                x, y = SkyCoord(ra, dec, unit=u.deg).to_pixel(wcs)
158                if 0 <= x < image.shape[1] and 0 <= y < image.shape[0]:
159                    objects_in_frame.append((ra, dec, obj_id, x, y))
160            except Exception:
161                continue
162
163        if not objects_in_frame:
164            # print("no object in frame")
165            return None
166
167        # Extract positions
168        positions = [(pos[3], pos[4]) for pos in objects_in_frame]
169
170        # Perform photometry
171        result_table = self._perform_aperture_photometry(image, positions)
172
173        # Add metadata
174        result_table["RA"] = [pos[0] for pos in objects_in_frame]
175        result_table["DEC"] = [pos[1] for pos in objects_in_frame]
176        result_table["ID"] = [pos[2] for pos in objects_in_frame]
177        result_table["DATE-OBS"] = header.get("DATE-OBS", "")
178
179        return result_table
180
181    def _perform_aperture_photometry(
182        self, image: np.ndarray, positions: List[tuple]
183    ) -> Table:
184        """
185        Perform aperture photometry at specified pixel positions.
186
187        Parameters
188        ----------
189        image : numpy.ndarray
190            2D image array extracted from the FITS file.
191        positions : list of tuple
192            Initial `(x, y)` pixel coordinates for photometry.
193
194        Returns
195        -------
196        astropy.table.Table
197            Table containing fluxes, magnitudes, and magnitude uncertainties
198            for all successfully measured sources.
199        """
200        # Centroid refinement
201        x_init, y_init = zip(*positions)
202        x_cent, y_cent = centroid_sources(
203            image, x_init, y_init, box_size=3, centroid_func=centroid_com
204        )
205
206        # Remove failed centroids
207        valid = ~np.isnan(x_cent)
208        positions = [(x_cent[i], y_cent[i]) for i in range(len(x_cent)) if valid[i]]
209
210        # Aperture definitions
211        aperture = CircularAperture(positions, r=self.aperture_radius)
212        annulus = CircularAnnulus(
213            positions, r_in=self.annulus_inner, r_out=self.annulus_outer
214        )
215
216        # Photometry
217        phot_table = aperture_photometry(image, [aperture, annulus])
218
219        # Background subtraction
220        bkg_mean = phot_table["aperture_sum_1"] / annulus.area
221        bkg_sum = bkg_mean * aperture.area
222        final_flux = phot_table["aperture_sum_0"] - bkg_sum
223
224        # Calculate magnitudes
225        magnitudes = self.zeropoint - 2.5 * np.log10(np.abs(final_flux))
226
227        # Error estimation
228        _, _, std = sigma_clipped_stats(image, sigma=3.0)
229
230        # Flux uncertainty (in electrons)
231        flux_uncertainty = np.sqrt(
232            np.abs(final_flux) / self.epadu
233            + aperture.area * std**2
234            + aperture.area * std**2 / annulus.area
235        )
236
237        # Magnitude uncertainty
238        mag_error = 1.0857 * flux_uncertainty / np.abs(final_flux)
239
240        # Create output table
241        result = Table()
242        result["flux"] = final_flux
243        result["mag"] = magnitudes
244        result["mag_err"] = mag_error
245        result["flux_err"] = flux_uncertainty
246
247        return result
248
249    def process_image(
250        self,
251        fits_file: str,
252        catalog_file: str,
253        target_star: StarData,
254        max_sep_arcsec: float = 0.5,
255    ) -> Optional[Table]:
256        """
257        Process a single FITS file and extract photometry for a target source.
258
259        This method filters all detected sources in the frame and keeps only
260        the source below a maximum angular separation threshold.
261
262        Parameters
263        ----------
264        fits_file : str
265            Path to a single TESS FITS image.
266        catalog_file : str
267            Path to the Gaia catalog CSV file.
268        target_star : StarData
269            Target star information.
270        max_sep_arcsec : float, optional
271            Maximum angular separation threshold.
272
273        Returns
274        -------
275        astropy.table.Table or None
276            Table containing the photometric measurements of the target
277            source in the frame, or None if the target was not detected.
278
279        Notes
280        -----
281        This method is safe for parallel execution and is intended to be used
282        inside multiprocessing workers.
283        """
284
285        target_ra = target_star.ra
286        target_dec = target_star.dec
287
288        catalog = self.load_catalog(catalog_file)
289
290        target_coord = SkyCoord(target_ra, target_dec, unit=u.deg)
291
292        try:
293            result = self.process_fits(fits_file, catalog)
294            if result is None or len(result) == 0:
295                return None
296
297            coords = SkyCoord(result["RA"], result["DEC"], unit=u.deg)
298
299            separations = target_coord.separation(coords)
300            idx = np.argmin(separations)
301
302            # Reject frame if the star is not detected
303            if separations[idx].arcsec > max_sep_arcsec:
304                print(
305                    f"[PID {os.getpid()}] "
306                    f'Target not found (sep={separations[idx].arcsec:.2f}") '
307                    f"in {os.path.basename(fits_file)}",
308                    flush=True,
309                )
310                return None
311
312            # Keep only if it contains the star
313            lambda_tau_row = result[idx : idx + 1]
314
315            return lambda_tau_row
316        except Exception as e:
317            print(f"[PID {os.getpid()}] Error processing {fits_file}: {e}")
318            return None

321def analyze_image(
322    fits_file: str,
323    catalog_file: str,
324    target_star: StarData,
325    max_sep_arcsec: float = 0.5,
326) -> Optional[Table]:
327    """
328    Analyze a single local FITS image and return the target's photometry row.
329
330    Convenience wrapper around ``TessPhotometry().process_image`` for
331    users that do not need to configure the photometry class.
332
333    Parameters
334    ----------
335    fits_file : str
336        Path to a single TESS FITS image.
337    catalog_file : str
338        Path to the Gaia catalog CSV file.
339    target_star : StarData
340        Target star information.
341    max_sep_arcsec : float, optional
342        Maximum angular separation threshold.
343
344    Returns
345    -------
346    astropy.table.Table or None
347        Table containing the photometric measurements of the target
348        source in the frame, or None if the target was not detected.
349    """
350    return TessPhotometry().process_image(
351        fits_file,
352        catalog_file,
353        target_star,
354        max_sep_arcsec=max_sep_arcsec,
355    )

111def create_catalog(
112    star_name: str,
113    radius_arcmin: float = 10.0,
114    output_file: str = "gaia_catalog.csv",
115    simbad_server: str | None = None,
116) -> StarData:
117    """
118    Create and save a Gaia source catalog centered on a target star.
119
120    Parameters
121    ----------
122    star_name : str
123        Object name resolvable by SIMBAD.
124    radius_arcmin : float, optional
125        Search radius in arcminutes.
126    output_file : str, optional
127        Output CSV filename.
128    simbad_server : str, optional
129        SIMBAD server to use. Defaults to the astroquery default
130        (``simbad.cds.unistra.fr``); use ``"simbad.harvard.edu"`` for the
131        Harvard mirror when the default is unreachable.
132
133    Returns
134    -------
135    star : StarData
136        StarData with name, right acension, declination and gaiadr3 id.
137    """
138
139    star = get_star(star_name, simbad_server=simbad_server)
140    catalog = query_gaia_catalog(star, radius_arcmin)
141
142    catalog.to_csv(output_file, index=False)
143    print(f"Catalog saved to {output_file} with {len(catalog)} stars")
144
145    return star

20def get_star(name: str, simbad_server: str | None = None) -> StarData:
21    """
22    Resolve an object name using SIMBAD and return a
23    StarData with its sky coordinates and id.
24
25    Parameters
26    ----------
27    name : str
28        Object identifier resolvable by SIMBAD.
29    simbad_server : str, optional
30        SIMBAD server to use. Defaults to the astroquery default
31        (``simbad.cds.unistra.fr``); use ``"simbad.harvard.edu"`` for the
32        Harvard mirror when the default is unreachable.
33
34    Returns
35    -------
36    star : StarData
37        Returns the found right ascension and declination, can
38        be accessed as star.ra and star.dec respectively.
39
40    """
41    if simbad_server is not None:
42        Simbad.server = simbad_server
43
44    Simbad.add_votable_fields("ids")
45    result = Simbad.query_object(name)
46
47    if result is None:
48        raise ValueError(f"Object '{name}' not found in SIMBAD")
49
50    ra = result["ra"][0]
51    dec = result["dec"][0]
52
53    # Get Gaia ID
54    gaia_id = None
55    for id_str in result["ids"][0].split("|"):
56        if "Gaia DR3" in id_str:
57            gaia_id = id_str.replace("Gaia DR3 ", "").strip()
58            break
59
60    star = StarData(
61        name=name,
62        ra=ra,
63        dec=dec,
64        gaia_id=gaia_id,
65    )
66
67    return star

 70def query_gaia_catalog(star: StarData, radius_arcmin: float = 1.0) -> pd.DataFrame:
 71    """
 72    Query the Gaia DR3 catalog around a given sky position.
 73
 74    Parameters
 75    ----------
 76    star : StarData
 77        Target star.
 78    radius_arcmin : float, optional
 79        Search radius in arcminutes.
 80
 81    Returns
 82    -------
 83    pandas.DataFrame
 84        Table containing Gaia sources within the search radius.
 85    """
 86    ra = star.ra
 87    dec = star.dec
 88
 89    radius_deg = radius_arcmin / 60.0
 90
 91    # 1. Format the coordinates correctly using astropy
 92    coord = SkyCoord(ra=ra, dec=dec, unit=(u.deg, u.deg), frame="icrs")
 93
 94    # 2. Use the native cone search
 95    job = Gaia.cone_search(coordinate=coord, radius=radius_deg * u.deg)
 96    results = job.get_results()
 97
 98    # query = f"""
 99    # SELECT ra, dec, source_id
100    # FROM gaiadr3.gaia_source
101    # WHERE CONTAINS(POINT('ICRS', ra, dec),
102    #                CIRCLE('ICRS', {ra}, {dec}, {radius_deg})) = 1
103    # """
104    #
105    # job = Gaia.launch_job_async(query)
106    # results = job.get_results()
107
108    return results.to_pandas()

20def get_tess_sectors(target_star: StarData) -> pd.DataFrame:
21    """
22    Return the list of TESS observing sectors covering a sky position.
23
24    Parameters
25    ----------
26    target_star : StarData
27        StarData with right acension and declination.
28
29    Returns
30    -------
31    pandas.DataFrame
32        Table describing available TESS sectors for the target position.
33    """
34
35    ra = target_star.ra
36    dec = target_star.dec
37
38    coord = SkyCoord(ra=ra, dec=dec, unit=(u.deg, u.deg))
39    sectors = Tesscut.get_sectors(coordinates=coord)
40    return sectors

43def download_tess_sector_script(sector: int) -> str:
44    """
45    Download the official shell script for a TESS sector from:
46    https://archive.stsci.edu/missions/tess/download_scripts/sector/
47
48    The script contains curl commands for downloading all FFIs
49    associated with the specified sector.
50
51    Parameters
52    ----------
53    sector : int
54        TESS sector number.
55
56    Returns
57    -------
58    str
59        Path to the downloaded shell script.
60    """
61    script_name = f"tesscurl_sector_{sector}_ffic.sh"
62
63    if not os.path.exists(script_name):
64        url = f"https://archive.stsci.edu/missions/tess/download_scripts/sector/{script_name}"
65        response = requests.get(url)
66        response.raise_for_status()
67
68        with open(script_name, "wb") as f:
69            f.write(response.content)
70
71    return script_name

 74def download_tess_image(shell_command: str, output_dir: str) -> str:
 75    """
 76    Download a single TESS FFI using the passed shell command.
 77
 78    Parameters
 79    ----------
 80    shell_command : str
 81        Shell command used to download the FITS file.
 82    output_dir : str
 83        Directory where the FITS file will be written.
 84
 85    Returns
 86    -------
 87    str
 88        Path to the downloaded FITS file.
 89    """
 90    filename = os.path.basename(shell_command.split()[-1])
 91    output_path = os.path.join(output_dir, filename)
 92
 93    cmd = f"{shell_command} -o {output_path}"
 94    # print(cmd)
 95    subprocess.run(
 96        cmd,
 97        shell=True,
 98        stdout=subprocess.DEVNULL,
 99        stderr=subprocess.DEVNULL,
100        cwd=output_dir,
101    )
102
103    return output_path

106def download_tess_images(
107    script_path: str, start_idx: int, num_images: int, output_dir: str = "./tess_images"
108) -> None:
109    """
110    Download a subset of TESS FFIs listed in a sector shell script.
111
112    Parameters
113    ----------
114    script_path : str
115        Path to the TESS download shell script.
116    start_idx : int
117        Index of the first line to process.
118    num_images : int
119        Number of images to download.
120    output_dir : str, optional
121        Directory where FITS files will be written.
122
123    Notes
124    -----
125    This function executes shell commands and suppresses stdout/stderr.
126    It is safe to call concurrently for non-overlapping index ranges.
127    """
128    os.makedirs(output_dir, exist_ok=True)
129
130    with open(script_path, "r") as f:
131        lines = [line.strip() for line in f]
132
133    # Get the requested range
134    download_lines = lines[start_idx : start_idx + num_images]
135
136    for i, line in enumerate(download_lines, 1):
137        download_tess_image(line, output_dir)

 24def load_photometry_data(
 25    csv_path: str, target_star: StarData, max_sep_arcsec: float = 0.5
 26) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
 27    """
 28    Load photometry data and extract a target light curve.
 29
 30    ``csv_path`` may be either a single CSV file containing one row per
 31    frame (appended photometry output) or a directory of per-frame CSV
 32    files. Sources below a maximum angular separation are selected.
 33
 34    Parameters
 35    ----------
 36    csv_path : str
 37        Path to a combined photometry CSV file or a directory of
 38        per-frame CSV files.
 39    target_star : StarData
 40        Target star with RA/Dec coordinates.
 41    max_sep_arcsec : float, optional
 42        Maximum allowed separation for a valid detection.
 43
 44    Returns
 45    -------
 46    times : numpy.ndarray
 47        Observation times in Julian Date, sorted chronologically.
 48    fluxes : numpy.ndarray
 49        Measured fluxes corresponding to the target.
 50    flux_errors : numpy.ndarray
 51        1-sigma flux uncertainties corresponding to the target.
 52
 53    Raises
 54    ------
 55    RuntimeError
 56        If no photometry data is found or the target is not detected
 57        in any frame.
 58    """
 59
 60    target_ra = target_star.ra
 61    target_dec = target_star.dec
 62
 63    target_coord = SkyCoord(target_ra, target_dec, unit=u.deg)
 64
 65    times = []
 66    fluxes = []
 67    flux_errors = []
 68
 69    if os.path.isdir(csv_path):
 70        csv_files = sorted(Path(csv_path).glob("*.csv"))
 71
 72        if not csv_files:
 73            raise RuntimeError("No CSV files found.")
 74
 75        for csv_file in csv_files:
 76            df = pd.read_csv(csv_file)
 77
 78            coords = SkyCoord(df["RA"].values, df["DEC"].values, unit=u.deg)
 79
 80            seps = target_coord.separation(coords).arcsec
 81            idx = np.argmin(seps)
 82
 83            if seps[idx] > max_sep_arcsec:
 84                continue
 85
 86            flux = df.loc[idx, "flux"]
 87            flux_err = df.loc[idx, "flux_err"]
 88            date_obs = df.loc[idx, "DATE-OBS"]
 89
 90            jd = Time(date_obs, format="isot", scale="utc").jd
 91
 92            fluxes.append(flux)
 93            flux_errors.append(flux_err)
 94            times.append(jd)
 95
 96    else:
 97        if not os.path.exists(csv_path):
 98            raise RuntimeError(f"Photometry file not found: {csv_path}")
 99
100        df = pd.read_csv(csv_path)
101
102        coords = SkyCoord(df["RA"].values, df["DEC"].values, unit=u.deg)
103
104        seps = target_coord.separation(coords).arcsec
105
106        for i in range(len(df)):
107            if seps[i] > max_sep_arcsec:
108                continue
109
110            flux = df.loc[i, "flux"]
111            flux_err = df.loc[i, "flux_err"]
112            date_obs = df.loc[i, "DATE-OBS"]
113
114            jd = Time(date_obs, format="isot", scale="utc").jd
115
116            fluxes.append(flux)
117            flux_errors.append(flux_err)
118            times.append(jd)
119
120    if not times:
121        raise RuntimeError("Target not found in any frame.")
122
123    times_arr = np.array(times)
124    fluxes_arr = np.array(fluxes)
125    flux_errors_arr = np.array(flux_errors)
126
127    order = np.argsort(times_arr)
128
129    return times_arr[order], fluxes_arr[order], flux_errors_arr[order]

132def compute_periodogram(
133    times: np.ndarray,
134    fluxes: np.ndarray,
135    flux_errors: Optional[np.ndarray] = None,
136    min_period: float = 0.1,
137    max_period: float = 10.0,
138    oversampling: float = 5.0,
139) -> dict:
140    """Compute Lomb-Scargle periodogram with weighted fitting and FAP.
141
142    Parameters
143    ----------
144    times : np.ndarray
145        Observation times (e.g., Julian Date).
146    fluxes : np.ndarray
147        Measured fluxes.
148    flux_errors : np.ndarray, optional
149        1-sigma flux uncertainties. If provided, passed to LombScargle
150        as `dy` for weighted periodogram (generalized LS).
151    min_period : float, optional
152        Minimum period to search (days).
153    max_period : float, optional
154        Maximum period to search (days).
155    oversampling : float, optional
156        Oversampling factor relative to the frequency resolution 1/T.
157        Default 5 follows VanderPlas (2018) recommendation.
158
159    Returns
160    -------
161    dict
162        Dictionary containing:
163        - 'periods': array of periods (days)
164        - 'power': array of periodogram power
165        - 'frequency': array of frequencies (1/days)
166        - 'peak_indices': list of peak indices sorted by power
167        - 'primary_period': best period or None
168        - 'primary_period_uncertainty': 1-sigma uncertainty or None
169        - 'primary_fap': False alarm probability (Baluev 2008) or None
170        - 'secondary_period': second best period or None
171    """
172    times = np.asarray(times)
173    fluxes = np.asarray(fluxes)
174
175    mask = ~np.isnan(times) & ~np.isnan(fluxes)
176    if flux_errors is not None:
177        flux_errors = np.asarray(flux_errors)
178        mask &= ~np.isnan(flux_errors)
179
180    times = times[mask]
181    fluxes = fluxes[mask]
182    if flux_errors is not None:
183        flux_errors = flux_errors[mask]
184
185    if len(times) < 3:
186        return {
187            "periods": np.array([]),
188            "power": np.array([]),
189            "frequency": np.array([]),
190            "peak_indices": [],
191            "primary_period": None,
192            "primary_period_uncertainty": None,
193            "primary_fap": None,
194            "secondary_period": None,
195        }
196
197    baseline = times.max() - times.min()
198    df = 1.0 / baseline / oversampling
199    f_min = 1.0 / max_period
200    f_max = 1.0 / min_period
201    frequency = np.arange(f_min, f_max, df)
202    periods = 1.0 / frequency
203
204    ls = LombScargle(times, fluxes, dy=flux_errors, normalization="standard")
205    power = ls.power(frequency)
206
207    # Peak finding with prominence and distance thresholds
208    power_std = np.std(power)
209    min_distance = max(1, int(1.0 / df))  # at least one independent resolution element
210    peaks, _ = find_peaks(power, prominence=3.0 * power_std, distance=min_distance)
211    sorted_peaks = sorted(peaks, key=lambda x: power[x], reverse=True)
212
213    results = {
214        "periods": periods,
215        "power": power,
216        "frequency": frequency,
217        "peak_indices": sorted_peaks,
218        "primary_period": None,
219        "primary_period_uncertainty": None,
220        "primary_fap": None,
221        "secondary_period": None,
222    }
223
224    if len(sorted_peaks) >= 1:
225        primary_idx = sorted_peaks[0]
226        results["primary_period"] = float(periods[primary_idx])
227        results["primary_period_uncertainty"] = estimate_period_uncertainty(
228            frequency, power, primary_idx, baseline
229        )
230        # Baluev (2008) analytic FAP
231        fap = ls.false_alarm_probability(
232            power[primary_idx],
233            method="baluev",
234            minimum_frequency=frequency.min(),
235            maximum_frequency=frequency.max(),
236        )
237        results["primary_fap"] = float(fap)
238
239    if len(sorted_peaks) >= 2:
240        results["secondary_period"] = float(periods[sorted_peaks[1]])
241
242    return results

 37def process_images_parallel(
 38    script_file: str,
 39    catalog_file: str,
 40    star: StarData,
 41    output_file: str | None = None,
 42    start_idx: int | None = None,
 43    max_workers: int | None = None,
 44    keep_images_dir: str | None = None,
 45):
 46    """
 47    Download and process TESS images in parallel for a target star.
 48
 49    This function coordinates the parallel execution of photometry over
 50    a list of TESS image URLs or commands contained in a script file.
 51    Each worker downloads exactly one FITS file, performs photometry,
 52    appends one row per detection to a single output CSV file, and
 53    deletes temporary files.
 54
 55    Parameters
 56    ----------
 57    script_file : str
 58        Path to a TESS download script (one image per line).
 59    catalog_file : str
 60        Path to a Gaia catalog CSV.
 61    star : StarData
 62        Target star information.
 63    output_file : str or None, optional
 64        Path of the single CSV file that photometry rows are appended to.
 65        Defaults to ``<star name>_<sector>.csv``.
 66    start_idx : int or None, optional
 67        Line index in the script file from which to start processing.
 68        If None (default) and the output file exists with data,
 69        processing resumes automatically after the last processed image.
 70        Passing an explicit value overrides the automatic resume.
 71    max_workers : int or None, optional
 72        Number of parallel worker processes. Defaults to the number of
 73        available CPU cores.
 74    keep_images_dir : str or None, optional
 75        Directory where FITS images in which the star was detected are
 76        saved. If None, no images are kept.
 77
 78    Notes
 79    -----
 80    - Each FITS file is handled independently.
 81    - Rows are appended to the output file under an exclusive file lock,
 82      so concurrent workers never corrupt it. Duplicate rows (same curl
 83      command) are skipped.
 84    - Temporary FITS files are stored in per-worker directories and
 85      deleted after processing.
 86    - Pressing ``Ctrl+C`` terminates all workers immediately and exits
 87      with status code 130.
 88    """
 89    if max_workers is None:
 90        max_workers = mp.cpu_count()
 91
 92    sector = _sector_from_script(script_file)
 93    if sector is not None:
 94        star.sector = sector
 95
 96    if output_file is None:
 97        safe_name = star.name.replace(" ", "_")
 98        if sector is not None:
 99            output_file = f"{safe_name}_{sector}.csv"
100        else:
101            output_file = f"{safe_name}.csv"
102
103    output_dir = os.path.dirname(output_file) or "."
104    os.makedirs(output_dir, exist_ok=True)
105
106    with open(script_file) as f:
107        lines = f.readlines()
108
109    n_images = len(lines)
110
111    if start_idx is None:
112        start_idx = _resume_start_idx(lines, output_file)
113
114    print(f"Processing images {start_idx} -> {n_images - 1}")
115    print(f"Using {max_workers} workers")
116
117    indices = list(range(start_idx, n_images))
118
119    worker = partial(
120        worker_process_fits,
121        script_file,
122        catalog_file=catalog_file,
123        star=star,
124        output_file=output_file,
125        sector=sector,
126        keep_images_dir=keep_images_dir,
127    )
128
129    _run_pool(worker, indices, max_workers)

395def process_local_images_parallel(
396    fits_dir: str,
397    catalog_file: str,
398    star: StarData,
399    output_file: str | None = None,
400    max_workers: int | None = None,
401    keep_images_dir: str | None = None,
402    pattern: str = "*.fits",
403):
404    """
405    Analyze FITS images already stored locally, in parallel.
406
407    This function performs photometry over all FITS files in a directory
408    without downloading anything. One row per detection is appended to a
409    single output CSV file, with the FITS header ``SECTOR`` and the local
410    file path (in the ``CURL`` column) as metadata.
411
412    Parameters
413    ----------
414    fits_dir : str
415        Directory containing the FITS images.
416    catalog_file : str
417        Path to a Gaia catalog CSV.
418    star : StarData
419        Target star information.
420    output_file : str or None, optional
421        Path of the single CSV file that photometry rows are appended to.
422        Defaults to ``<star name>.csv``.
423    max_workers : int or None, optional
424        Number of parallel worker processes. Defaults to the number of
425        available CPU cores.
426    keep_images_dir : str or None, optional
427        Directory where FITS images in which the star was detected are
428        copied. If None, no images are kept.
429    pattern : str, optional
430        Glob pattern used to select files inside ``fits_dir``. Defaults
431        to ``"*.fits"``.
432
433    Notes
434    -----
435    - Images whose path is already stored in the output file are skipped,
436      so re-running resumes automatically. The deduplication keys on the
437      exact path string; different path spellings for the same file will
438      not be recognized as duplicates.
439    - Rows are appended under an exclusive file lock.
440    - Pressing ``Ctrl+C`` terminates all workers immediately and exits
441      with status code 130.
442    """
443    if max_workers is None:
444        max_workers = mp.cpu_count()
445
446    fits_files = sorted(glob.glob(os.path.join(fits_dir, pattern)))
447
448    if not fits_files:
449        raise FileNotFoundError(
450            f"No FITS files found in {fits_dir} (pattern {pattern!r})"
451        )
452
453    if output_file is None:
454        output_file = f"{star.name.replace(' ', '_')}.csv"
455
456    output_path_dir = os.path.dirname(output_file) or "."
457    os.makedirs(output_path_dir, exist_ok=True)
458
459    processed = _processed_sources(output_file)
460    if processed:
461        already = set(processed)
462        remaining = [f for f in fits_files if f not in already]
463        print(f"Skipping {len(fits_files) - len(remaining)} already-processed images")
464        fits_files = remaining
465
466    print(f"Processing {len(fits_files)} local images")
467    print(f"Using {max_workers} workers")
468
469    worker = partial(
470        worker_process_local_fits,
471        catalog_file=catalog_file,
472        star=star,
473        output_file=output_file,
474        keep_images_dir=keep_images_dir,
475    )
476
477    _run_pool(worker, fits_files, max_workers)