HARPSpol Reduction Pipeline

Reducing all HARPS polarimetry data from the ESO archive using PyReduce. The goal is a consistent, polarimetry-focussed reduction of all public data for publication on PolarBase.

Scope (v1): circular polarimetry only, and only sequences of exactly four exposures. Both restrictions are enforced at the archive-query step (01_query_archive.py) so rejected data never enters the pipeline:

• ECHELLE,CIRPOL only. Linear polarimetry needs an inverted beam ratio and a Q-vs-U branch that does not exist yet; running it through the circular formula would silently produce wrong Stokes parameters. LINPOL is deferred to v2. 04_demodulate.py also refuses LINPOL outright, in case such data is already on disk from an earlier query.
• Exactly 4 exposures per TPL START. The ratio method consumes one 4-exposure QWP cycle. Longer templates (most often 8, i.e. two cycles filed under one TPL START) would need splitting into cycles first, which the demodulator does not do. They are rejected and counted at query time rather than being dropped without a record further down.

Requirements

• uv (manages Python 3.13+ and all dependencies)
• GNU parallel (for batch processing)
• An ESO archive account. config.py defaults to alavail; override without editing it by setting ESO_USERNAME in the environment.

All Python dependencies are managed automatically by uv via pyproject.toml.

Important version notes:

• astroquery >= 0.4.11 is required. Versions <= 0.4.7 are broken due to ESO archive API changes (July 2025). The login API also changed to keyword-only arguments in 0.4.8+.
• pyreduce-astro >= 0.9. The 0.9 series replaced the extraction C backend with slitdec; validated against 0.8a4 on HD 120411 (2026-03-18): intensity spectra agree to <1% per order, and the blue-channel Null improves from systematics-dominated (|N|/σ ≈ 3.2) to noise-consistent (0.84). Blue-channel products extracted with 0.8a4 should be considered suspect and re-reduced.
• 0.9 proper (2026-07-31) also made the scatter step usable — see "Scattered-light subtraction" below. Products reduced with 0.9b1 carry no scattered-light correction.

Tested with: astroquery 0.4.11, pyreduce 0.9, barycorrpy 0.4.4, numpy 2.4, scipy 1.17, astropy 7.2.

Reproducibility

uv.lock is committed. It pins the exact resolved version of all 63 dependencies, so uv sync reconstructs the environment a given reduction was produced with — the versions above are prose, the lockfile is the record.

One dependency it cannot pin: PyReduce itself. It is an editable install from a local path (pyreduce-astro = { path = "../PyReduce", editable = true }), so the lockfile records the path, not a version. Since the extraction backend is the single largest influence on the products, the PyReduce commit has to be recorded separately when a reduction is published:

git -C ../PyReduce describe --tags --always --dirty

At the time of writing that is v0.9-10-g0269ec4. Quote that string, not just "0.9" — the 0.9 series replaced the extraction C backend, and intermediate commits have changed HARPSPOL behaviour (e.g. the wavecal linelist order convention, and whether the scatter step masks anything).

Setup

git clone ssh://git@codeberg.org/verveine/harpspol.git
cd harpspol
uv sync

Tests

uv run --group dev pytest          # ~2 s, no network, no data required

The suite covers the numerical core on synthetic inputs — the continuum fit, the physical-order identification, the flat-blaze row layout, the sequence QC and BERV handling, plus an end-to-end demodulation of a synthetic 4-exposure sequence. These are the functions that have each been got wrong at least once (see TODO.md), so the tests are written to fail if the specific historical bug is reintroduced, not merely to exercise the code. Each was checked by re-introducing the bug and confirming the suite goes red.

PyReduce is expected at ../PyReduce (editable install via pyproject.toml). Adjust the path in [tool.uv.sources] if your checkout is elsewhere.

Pipeline Overview

Script	Purpose
01_query_archive.py	Query ESO archive for HARPSpol science observations (CIRPOL, 4-exposure sequences)
02_download.py	Download science and calibration FITS files
03_reduce.py	Run PyReduce (calibrations and/or science extraction)
04_demodulate.py	Polarimetric demodulation (4-exposure Stokes V/I)
05_combine_channels.py	Merge blue + red demodulated spectra
filter_target.py	Subset a datasets.ecsv to a single star (single-target runs)
check_wavecal.py	Check wavelength-calibration consistency across a target's dates
plot_stokes.py	Quick-look I/V/Null plot of a demodulated file
config.py	Shared configuration (paths, channels, constants)
run_all.sh	Batch runner with GNU parallel

Quick Start

1. Query the ESO archive

# All HARPSpol data from 2012
uv run python 01_query_archive.py --start 2012-01-01 --end 2013-01-01

# Filter by ESO program ID
uv run python 01_query_archive.py --start 2012-01-01 --end 2013-01-01 --program 187.D-0917

Outputs datasets.ecsv (QC-approved science files) and calib_nights.ecsv (calibration nights). The run prints how many sequences were rejected and with what exposure counts — worth recording, since that is the pipeline's completeness figure.

2. Download data

uv run python 02_download.py --all                # science + calibrations
uv run python 02_download.py --science             # science only
uv run python 02_download.py --calibs              # calibrations only
uv run python 02_download.py --calibs --night 2012-07-14  # single night

Files are organized as:

• data/science/{TARGET}_{TPL_START}/raw/ — 4 raw science FITS per observation
• data/calibs/{YYYY-MM-DD}/raw/ — bias, flat, ThAr for each night

3. Reduce calibrations

uv run python 03_reduce.py calibs --night 2012-07-14

Runs: bias, flat, trace, curvature, norm_flat, wavecal_master, wavecal. Output goes to data/calibs/{night}/reduced/.

Add --scatter to also run the scatter step and produce scatter.npz (see "Scattered-light subtraction" below). The full set a night must produce before science extraction will use it is CALIB_PRODUCTS in config.py — bias, flat, traces and flat_norm per channel, plus wavecal_master and linelist per channel and per beam (16 files, or 18 with --scatter). 03_reduce.night_has_all_calibs checks exactly that set, and required_calib_names(with_scatter=) selects which of the two it is.

Scattered-light subtraction

Scattered light is unpolarised, so leaving it in dilutes V/I — a multiplicative bias that always underestimates the field. The step is enabled with --scatter on both calibs and science (or once on run_all.sh): the calibration run produces scatter.npz, and the science run has to be told to require and link it.

It was off by default until 2026-08-05 because it destroyed the products, which turned out to be two independent defects in PyReduce rather than a tuning problem:

1. The model was fitted on the LAMP,LAMP,TUN flat and subtracted from science frames unscaled. Fixed upstream in 0.9 — consumers now re-estimate on the calibrated frame they are about to extract.
2. A fractional extraction_height resolved to a zero-height aperture, so no trace was masked and the polynomial was fitted to the order flux itself. Our trace files store each trace twice (grouped + raw_traces, identical polynomials for HARPSPOL), which drove the median order spacing to 0. Fixed in ivh/PyReduce#40.

With both fixed, the fitted background sits within 20% of the frame's own inter-order floor (1.18x blue, 1.02x red at extraction_height: 0.7), and on γ Equ the correction removes 2.9% of the blue and 1.4% of the red median intensity while leaving the V-to-null rms ratio unchanged.

4. Extract science spectra

uv run python 03_reduce.py science --obs HD-96446_2012-07-14T22-55-44
uv run python 03_reduce.py science --obs HD-96446_2012-07-14T22-55-44 --berv  # barycentric WAVE

Automatically finds the nearest calibration night (nearest complete one within --max-calib-days, default 30), symlinks products, and runs optimal extraction. Output goes to data/science/{obs}/reduced/.

The extracted WAVE is topocentric, so the same stellar lines sit at different wavelengths across dates (Earth's barycentric motion, up to ~0.45 Å over a year). With --berv, WAVE is barycentric-corrected in place (WAVE * (1 + BERV/c), BERV via barycorrpy) and a BERVCORR header flag is set. 04_demodulate reads that flag and skips its own BERV correction, so the demodulated products are not double-corrected. The correction is idempotent (re-running is a no-op) and reversible (topocentric WAVE = WAVE / (1 + BERV/c), BERV is in the header).

One BERV per sequence, evaluated at mid-sequence. BERV is computed once, at the midpoint of the 4-exposure cycle (start of the first exposure to the end of the last), and the same value is applied to all four exposures. Both parts matter:

• The demodulator assumes the four exposures share one wavelength grid — it takes the grid from the first exposure and never re-interpolates the lower beams. A per-exposure BERV would leave them on grids differing by ~20 m/s across a cycle.
• BERV drifts by ~20 m/s over a cycle, so referencing the first exposure would put the wavelength frame ~10 m/s away from the DEMOD JD_UTC MIDSEQ epoch the product is stamped with.

The epoch used is recorded as BERVJD. A file already carrying BERVCORR with a different BERV (e.g. reduced by an earlier per-exposure version) is re-referenced by (1 + BERV_new/c) / (1 + BERV_old/c) rather than skipped, so existing products converge to the sequence frame without a full re-reduction.

5. Demodulate

uv run python 04_demodulate.py --obs HD-96446_2012-07-14T22-55-44
uv run python 04_demodulate.py --obs HD-96446_2012-07-14T22-55-44 --plot  # with diagnostics
uv run python 04_demodulate.py --obs HD-96446_2012-07-14T22-55-44 --channel BLUE  # single channel
uv run python 04_demodulate.py --obs HD-96446_2012-07-14T22-55-44 --debug   # re-raise instead of skipping

Computes Stokes V/I, null spectrum, continuum normalization, and barycentric correction (one BERV per sequence, at mid-sequence — see above). Output: data/science/{obs}/demodulated/demodulated_{blue,red}.fits. Linear-polarimetry sequences are refused with a message rather than demodulated.

6. Combine channels

uv run python 05_combine_channels.py --obs HD-96446_2012-07-14T22-55-44

Concatenates blue and red into data/science/{obs}/demodulated/, named after the first science exposure — e.g. HARPS.2013-08-02T...658_demodulated.fits.

Batch Processing

Run everything

./run_all.sh                       # all steps, 4 parallel jobs
./run_all.sh -j 8                  # 8 parallel jobs
./run_all.sh --program 187.D-0917  # filter by program
./run_all.sh --berv                # barycentric-correct WAVE during extraction
./run_all.sh --max-calib-days 45   # widen the calibration search window
./run_all.sh --help                # full usage

--berv and --max-calib-days are forwarded to 03_reduce.py science, so the barycentric products described above are reachable from the batch runner and not only from a single-observation invocation.

A step with no work to do prints a message and is skipped; it does not stop the run, so --step science,demod,combine on a fresh tree still reaches every stage.

Run specific steps

./run_all.sh --step calibs                  # calibrations only
./run_all.sh --step science                 # science extraction only
./run_all.sh --step science,demod,combine   # skip calibs
./run_all.sh --step demod -j 12             # re-run demodulation, 12 jobs

Available steps: calibs, science, demod, combine (or all).

Process all observations for a program

# Demodulate + combine all observations from program 187.D-0917
./run_all.sh --step demod,combine --program 187.D-0917 -j 8

Process a batch with the individual scripts

# All observations, parallel science extraction
ls -d data/science/*/raw | xargs -I{} dirname {} | xargs -n1 basename | \
    parallel -j4 uv run python 03_reduce.py science --obs {}

# Demodulate all reduced observations
uv run python 04_demodulate.py --all
uv run python 04_demodulate.py --all --program 187.D-0917

# Combine all demodulated observations
uv run python 05_combine_channels.py --all
uv run python 05_combine_channels.py --all --program 187.D-0917

Program ID Filtering

All scripts support --program to filter by ESO program ID (substring match):

uv run python 01_query_archive.py --program 187.D-0917   # filter archive query
uv run python 02_download.py --all --program 187.D-0917   # download only this program
uv run python 04_demodulate.py --all --program 187.D-0917 # demodulate only this program
./run_all.sh --program 187.D-0917                         # full pipeline for one program

For 03_reduce.py, the --program flag works with the all subcommand:

uv run python 03_reduce.py all --program 187.D-0917

Single-Target Reduction

The query and download scripts filter by date and program, not by star, so a single-target run goes through filter_target.py, which subsets a datasets.ecsv to one object. Archive OBJECT headers are inconsistent (observer-set), so it matches a list of normalized aliases plus distinctive substring tokens (e.g. an HD number). Example for eps Ind (HD 209100), observed 2013–2019:

# 1. Query the whole archive over the target's era (all stars in the window)
uv run python 01_query_archive.py --start 2013-01-01 --end 2020-01-01

# 2. Subset datasets.ecsv to eps Ind -> eps_ind.ecsv (prints matched names,
#    sequence count, programs, and the calibration nights needed)
uv run python filter_target.py \
    --aliases "HD-209100,eps Ind,epsInd,HD209100,GJ845,HIP108870" \
    --tokens 209100 -o eps_ind.ecsv

# 3. Download science + its calibration nights, then reduce/demodulate/combine
uv run python 02_download.py --all --datasets eps_ind.ecsv
uv run python 03_reduce.py all
uv run python 04_demodulate.py --all
uv run python 05_combine_channels.py --all

If nothing matches, filter_target.py prints every distinct OBJECT name in the input (with its normalized form) so you can spot the archive's spelling.

Directory Structure

harpspol/
├── config.py                    # shared configuration
├── 01_query_archive.py          # ESO archive query
├── 02_download.py               # data download
├── 03_reduce.py                 # PyReduce calibration + science
├── 04_demodulate.py             # polarimetric demodulation
├── 05_combine_channels.py       # merge blue + red channels
├── run_all.sh                   # batch runner
├── pyproject.toml               # uv project (dependencies)
├── datasets.ecsv                # (generated) science file list
├── calib_nights.ecsv            # (generated) calibration nights
│
├── data/
│   ├── calibs/
│   │   └── {YYYY-MM-DD}/
│   │       ├── raw/             # raw calibration FITS
│   │       └── reduced/         # PyReduce output (bias, flat, traces, wavecal)
│   │
│   └── science/
│       └── {TARGET}_{TPL_START}/
│           ├── raw/             # 4 raw science FITS
│           ├── reduced/         # symlinked calibs + extracted science spectra
│           └── demodulated/     # demodulated Stokes spectra
│
├── doc/                         # paper catalog scripts & data
└── paper/                       # manuscript

Running against another tree

The layout above is the default and needs no configuration. It is not the only one: the deployed PolarBase reduction files its data by year and month, and keeps raw and reduced products under separate roots. config.py resolves both through a Layout object, selected by environment variable, and every path in the pipeline goes through it.

	HARPSPOL_LAYOUT=flat (default)	HARPSPOL_LAYOUT=dated
calibrations	{base}/calibs/{night}/{raw,reduced}	{root}/{YYYY}/{YYYY-MM}/calibration/{night}
science	{base}/science/{TARGET}_{TPL_START}/{raw,reduced,demodulated}	{root}/{YYYY}/{YYYY-MM}/science/{TPL_START}
raw vs reduced	side by side under one root	separate roots (IN and REDUCTION)
target name	parsed from the directory name	read from the OBJECT header

# the default tree, rooted somewhere other than ./data
HARPSPOL_BASE_DIR=/scratch/harpspol ./run_all.sh

# the deployed tree
export HARPSPOL_LAYOUT=dated
export HARPSPOL_RAW_ROOT=$BAY/IN/harpspol
export HARPSPOL_REDUCED_ROOT=$BAY/REDUCTION/harpspol
# or just BAY_VOLUME, from which both roots are derived
./run_all.sh

Two subtleties the dated layout has to get right. Science directories are filed by the calendar date of TPL START, which is how the download step places them, while calibrations are matched on the observing night — a sequence at 02:00 on the 1st of a month lives under that month but belongs to the night before. And because the directory name carries no target, the target is read from OBJECT; where that fails the extraction matches any target, which is safe because the directory holds exactly one sequence.

Demodulated Output Format

The output FITS files contain a binary table with these columns:

Column	Unit	Description
WAVE	nm	Wavelength (barycentric-corrected)
I	ADU	Stokes I (total intensity)
I/Ic		Continuum-normalized intensity
ERR_I	ADU	Error on I
ERR_I/Ic		Error on I/Ic
Ic	ADU	Continuum
STOKES		Stokes V/I (or Q/I, U/I)
ERR_STOKES		Error on Stokes parameter
STOKES/Ic		Stokes parameter scaled by normalized intensity
ERR_STOKES/Ic		Error on scaled Stokes
NULL		Null spectrum (diagnostic)
ERR_NULL		Error on null
NULL/Ic		Null scaled by normalized intensity
ERR_NULL/Ic		Error on scaled null
ORDER		Physical echelle order label (e.g., blue161, red089)

Provenance keywords

The primary header records how the product was made. The ones worth knowing:

Keyword	Meaning
DEMOD_VERS_ID	Demodulator version (6.7 at the time of writing)
DEMOD REDUCTION WHEN	Date of the run that produced the file
DEMOD QWP PAIRING	Which exposure split was used, e.g. 1+4 vs 2+3 (see step 5)
DEMOD QWP ANG1..4	The retarder angles it was derived from
DEMOD STOKES PARAM	V — circular only in v1
DEMOD ARCFILE1..4	The four raw exposures that went in
DEMOD BLAZE SOURCE	flat field or the science-envelope fallback
DEMOD JD_UTC MIDSEQ	Mid-sequence epoch the product is stamped with
DEMOD BERV	Barycentric velocity applied, m/s, at that epoch
DEMOD ORDER MIN/MAX	Physical echelle order range
DEMOD SNRMAX BLUE/RED/FULLSPECTRUM	Peak SNR
NOCONTN, NOCONT	Count and list of orders with no usable continuum

DEMOD QWP PAIRING is the one to tally across a batch: it is how the fraction of ABBA sequences in a programme (or in the archive) gets established.

Orders with no usable continuum

Every order that was demodulated is kept in the table. But an order can have no continuum to normalize against — in the far blue of a cool star there is essentially no flux, and the summed intensity scatters about zero. There is no honest value for Ic there, so Ic and the four /Ic columns (I/Ic, ERR_I/Ic, STOKES/Ic, ERR_NULL/Ic …) are set to NaN. WAVE, I, ERR_I, STOKES, ERR_STOKES, NULL and ERR_NULL are unaffected.

Two header keywords record this:

Keyword	Meaning
NOCONTN	number of orders with no usable continuum (0 for most targets)
NOCONT	comma-separated list of those order labels (absent if NOCONTN = 0)

The test is on the data, not on a fixed order list: it fires only when an order has no positive flux at all (95th percentile and median of the summed intensity both ≤ 0). Hot stars have real flux in the far blue, so their blue orders never reach this branch. Verified on HD 120411 and γ Equ — NOCONTN = 0 in both channels, and all science columns bit-identical to the v6.6 products. For eps Ind (K5V) it flags blue154–blue161.

Before v6.7 these orders were emitted with a placeholder continuum — the 95th percentile of the noise, or literally Ic = 1 ADU when even that was negative — which turned raw ADU into the normalized column and produced I/Ic values of several thousand, with the wrong sign. Products written before v6.7 should be re-demodulated if they contain cool stars.

How It Works

This pipeline replaces the previous 11-step workflow that required dual-path calibration (1-beam for wavecal, 2-beam for traces) and manual copying of products between directories. PyReduce now handles dual-beam instruments natively via fibers_per_order: 2 in the instrument config, which auto-pairs the two beams from the polarising beam splitter.

The demodulation algorithm (ratio method for 4 circular polarimetry exposures):

1. Extract upper/lower beam spectra from each exposure using the GROUP and M columns. PyReduce numbers traces sequentially, so M is not the diffraction order; the physical echelle order is recovered from the central wavelengths of consecutive traces via m·λ(m) = (m+1)·λ(m+1) and used for the ORDER label.

2. Divide out the blaze. The instrumental blaze comes from the flat field (flat_norm.npz), per beam — upper and lower are each divided by their own flat blaze. If the flat products are missing, fall back to an envelope estimated from the science data itself (rolling 95th-percentile filter plus light Gaussian smoothing). The fallback is second choice by design: an estimator built from the star's own spectrum cannot separate the blaze from broad stellar structure.

3. Cross-correlate the beams for a flux scale and a wavelength offset (χ² fit, offset bounded to ±2 Å), then interpolate the upper beam onto the lower beam's wavelength grid.

4. Scale each sub-exposure to the sequence median. These factors cancel out of the double ratio below by construction, so they do not bias V/I.

5. Pair the exposures by retarder state. The ratio method puts the two exposures sharing a quarter-wave-plate state over the two sharing the other, and which exposures those are is not fixed — the archive holds at least two QWP orderings:

   QWP angles	states (mod 180)	pairing
   45, 135, 225, 315	A B A B	1+3 against 2+4
   45, 135, 315, 225	A B B A	1+4 against 2+3

   qwp_pairing derives the split per sequence from ESO INS RET25 POS (a QWP is periodic in 180°, so the state is the angle modulo 180). The keyword is required: a sequence whose angles are missing, or which do not split two-and-two, is refused rather than guessed at. The ordering tracks the observing programme, not the date. A hardcoded pairing is wrong for one of the two orderings and largely cancels the polarisation — measured on γ Equ, deliberately mis-pairing drops rms |V/I| by a factor of 2.8. The split actually used is recorded as DEMOD QWP PAIRING.

6. With r_i = u_i / d_i and the two state groups (a1, a2) and (b1, b2) from step 5:

   • R = (r_a1 · r_a2) / (r_b1 · r_b2), then V/I = (R^{1/4} - 1) / (R^{1/4} + 1).
   • RN = (r_a1 · r_b1) / (r_a2 · r_b2) — one exposure of each state on either side, so a real signal cancels and the null is left showing the noise floor.

7. Propagate errors following Bagnulo et al. 2009 (PASP 121, 993) Eq. A10, with R^{1/2N} = R^{1/4} for N = 2 pairs.

8. Fit the continuum on the summed 8-beam intensity: an iterative degree-2 polynomial with asymmetric clipping, where the clip-band RMS is computed once from the initial smoothed fit and held fixed. Orders dominated by a broad line (Balmer series, Ca II H&K — FLAT_CONTINUUM_ORDERS) have no continuum window and get a degree-0 fit instead. Noise-dominated orders fall back to a flat robust level, and orders with no positive flux at all get Ic = NaN (see above).

9. Trim 15 pixels from each order edge, apply the barycentric correction (one BERV per sequence, at the mid-sequence epoch), and convert wavelengths from Å to nm.

Contact

Alexis Lavail