From 32f960b840122bd14f47440f69ceac44c773afa4 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 3 Mar 2026 01:58:27 +0000 Subject: [PATCH 1/4] Initial plan From b0676ede95af2cd65dec3af34d09eadbd29cd4d1 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 3 Mar 2026 02:10:58 +0000 Subject: [PATCH 2/4] Refactor: remove Python 2 compat, improve docstrings, update docs Co-authored-by: danielreardon <40694862+danielreardon@users.noreply.github.com> --- docs/source/dynspec.rst | 87 ++++++++++- docs/source/examples.rst | 267 ++++++++++++++++++++++++++++++++++ docs/source/index.rst | 2 + scintools/__init__.py | 2 - scintools/dynspec.py | 3 - scintools/scint_models.py | 162 ++++++++++++++++++++- scintools/scint_sim.py | 108 +++++++++++--- scintools/scint_utils.py | 297 ++++++++++++++++++++++++++++++++++---- 8 files changed, 866 insertions(+), 62 deletions(-) create mode 100644 docs/source/examples.rst diff --git a/docs/source/dynspec.rst b/docs/source/dynspec.rst index a86c13e..9bb3818 100644 --- a/docs/source/dynspec.rst +++ b/docs/source/dynspec.rst @@ -1,35 +1,106 @@ Dynspec class ============= -Note: documentation is currently being developed by a student. In the meantime, please email Daniel for any questions on usage: dreardon@swin.edu.au +The :class:`~scintools.dynspec.Dynspec` class is the primary interface for +loading, processing, and analysing pulsar dynamic spectra in scintools. -This will serve as a list of features available within the Dynspec class, in dynspec.py +For questions or support please email Daniel Reardon: dreardon@swin.edu.au + +Overview +-------- + +A dynamic spectrum is a two-dimensional array of flux density as a function of +time and radio frequency. The :class:`~scintools.dynspec.Dynspec` class +provides methods to: + +- Load dynamic spectra from ``psrflux``-format files or from + :class:`~scintools.scint_sim.Simulation` objects. +- Clean, crop, normalise, and resample data. +- Compute secondary spectra and fit parabolic arc curvatures. +- Compute autocorrelation functions (ACFs) and measure scintillation + timescales and bandwidths. +- Visualise all intermediate and final data products. + +See :doc:`examples` for complete, runnable worked examples. Importing dynamic spectra ------------------------- -A Dynspec object can be loaded from a `psrflux`-format dynamic spectrum, a simulation object, or ... +A :class:`~scintools.dynspec.Dynspec` object can be loaded from a +``psrflux``-format dynamic spectrum file:: -blah blah this is how you do things + from scintools.dynspec import Dynspec + dyn = Dynspec(filename='myobs.dynspec', verbose=True) +It can also be initialised directly from a +:class:`~scintools.scint_sim.Simulation` object:: + + from scintools.scint_sim import Simulation + from scintools.dynspec import Dynspec + sim = Simulation() + dyn = Dynspec(dyn=sim, verbose=True) Basic plotting -------------- -To view the dynamic spectrum, use +To view the dynamic spectrum use:: dyn.plot_dyn() +Additional keyword arguments control the colour scale, axis labels, and +whether to display or save the figure. + Processing ---------- -the dynamic spectrum can be cleaned, cropped, normalised, and resampled in several ways +The dynamic spectrum can be cleaned, cropped, normalised, and resampled in +several ways:: + dyn.trim_edges() # remove zero-padded edge channels/subints + dyn.zap() # interactive RFI flagging + dyn.refill() # interpolate over flagged data + dyn.crop_dyn(tmin, tmax, fmin, fmax) # crop in time and/or frequency + dyn.scale_dyn() # normalise each sub-integration Secondary spectra ----------------- -Calculation of the secondary spectrum. Includes pre-whitening and post-darkening by default .... +The secondary spectrum is computed via a 2-D Fourier transform of the dynamic +spectrum. Pre-whitening and post-darkening are applied by default:: + + dyn.calc_sspec() + dyn.plot_sspec() + +Arc curvature +------------- + +Parabolic arc curvatures are measured using +:meth:`~scintools.dynspec.Dynspec.fit_arc`, which searches for the arc +curvature :math:`\eta` that maximises power along the arc:: + + dyn.calc_sspec() + dyn.fit_arc(plot=True) + print(dyn.eta, dyn.etaerr) + +See :doc:`examples` for a complete arc-fitting workflow and guidance on +selecting input parameters. + +ACF analysis +------------ + +The autocorrelation function is computed with +:meth:`~scintools.dynspec.Dynspec.calc_acf` and fitted with +:meth:`~scintools.dynspec.Dynspec.fit_acf`:: + + dyn.calc_acf() + dyn.fit_acf() + print(dyn.tau, dyn.dnu) # scintillation timescale and bandwidth + +API reference +------------- +.. autoclass:: scintools.dynspec.Dynspec + :members: + :undoc-members: + :show-inheritance: -arcs are cool diff --git a/docs/source/examples.rst b/docs/source/examples.rst new file mode 100644 index 0000000..72a3ce3 --- /dev/null +++ b/docs/source/examples.rst @@ -0,0 +1,267 @@ +Examples +======== + +This page contains practical, runnable examples demonstrating common +scintools workflows. All examples assume that scintools and its dependencies +are installed (see the project `README `_). + +.. contents:: Contents + :local: + :depth: 2 + + +Loading and visualising a dynamic spectrum +------------------------------------------ + +**Loading from file** + +Dynamic spectra stored in ``psrflux`` format can be loaded directly:: + + from scintools.dynspec import Dynspec + + dyn = Dynspec(filename='J0437-4715.dynspec', verbose=True) + +The ``verbose=True`` flag (default) prints a summary of the observation +metadata (pulsar name, MJD, frequency, bandwidth, etc.) on load. + +**Basic plotting** + +Display the dynamic spectrum with the default colour scale:: + + dyn.plot_dyn() + +Save the figure to a file instead of displaying it:: + + dyn.plot_dyn(filename='dynspec.png', display=False) + +Customise the colour limits:: + + dyn.plot_dyn(vmin=0, vmax=5) + +**Loading from a simulation** + +A :class:`~scintools.scint_sim.Simulation` object can be passed directly:: + + from scintools.scint_sim import Simulation + from scintools.dynspec import Dynspec + + sim = Simulation(mb2=2, rf=1, ds=0.01, alpha=5/3) + dyn = Dynspec(dyn=sim, verbose=True) + dyn.plot_dyn() + + +Processing dynamic spectra +-------------------------- + +**Cleaning and normalisation** + +Remove short or bad sub-integrations and normalise each sub-integration to +remove bandpass structure:: + + dyn.trim_edges() # strip zero-padded boundary channels/sub-integrations + dyn.refill() # interpolate over flagged (zero) channels + dyn.scale_dyn() # normalise each sub-integration + +**Cropping** + +Crop to a sub-band and time range (units: MHz and minutes):: + + dyn.crop_dyn(tmin=0, tmax=30, fmin=1300, fmax=1500) + +**Resampling** + +Reduce frequency or time resolution by averaging:: + + dyn.zap_minmax() # flag outlier channels/sub-integrations + dyn.refill() + +**Common pitfalls** + +- Always call ``refill()`` after flagging to ensure that masked pixels do + not propagate as zeros into the FFT. +- Normalisation (``scale_dyn``) should be applied before computing the + secondary spectrum or ACF. + + +Computing secondary spectra +--------------------------- + +**Basic calculation** + +Compute the secondary spectrum using the default pre-whitening and +post-darkening:: + + dyn.calc_sspec() + dyn.plot_sspec() + +The secondary spectrum is stored in ``dyn.sspec`` with conjugate variables +``dyn.tdel`` (differential delay, μs) and ``dyn.fdop`` (differential Doppler +frequency, mHz). + +**Wavelength-scaled secondary spectrum** + +For a display that is independent of the observing frequency, use equal +steps in wavelength (``lamsteps=True``):: + + dyn.calc_sspec(lamsteps=True) + dyn.plot_sspec(lamsteps=True) + +**Interpreting results** + +Parabolic arcs in the secondary spectrum arise from a thin scattering screen. +The curvature :math:`\eta` (in units of s³) depends on the effective velocity +and the geometry of the scattering screen. See +:meth:`~scintools.dynspec.Dynspec.fit_arc` for how to measure :math:`\eta`. + + +Arc fitting workflow +-------------------- + +**Complete example** + +:: + + from scintools.dynspec import Dynspec + + dyn = Dynspec(filename='J0437-4715.dynspec') + dyn.trim_edges() + dyn.refill() + dyn.scale_dyn() + dyn.calc_sspec() + dyn.fit_arc(plot=True) + + print(f"Arc curvature eta = {dyn.eta:.4f} +/- {dyn.etaerr:.4f} s^3") + +**Parameter selection guide** + ++-------------------+------------------------------------------------------------+ +| Parameter | Guidance | ++===================+============================================================+ +| ``delmax`` | Maximum differential delay (μs) to include in the fit. | +| | Set to ~half the total delay range of your spectrum. | ++-------------------+------------------------------------------------------------+ +| ``numsteps`` | Number of trial curvatures; 10 000 is sufficient for most | +| | observations. | ++-------------------+------------------------------------------------------------+ +| ``startbin`` | Rows near zero-delay to exclude (avoids DC artefacts). | ++-------------------+------------------------------------------------------------+ +| ``cutmid`` | Columns near zero-Doppler to exclude (avoids DC artefacts).| ++-------------------+------------------------------------------------------------+ +| ``constraint`` | ``[eta_min, eta_max]`` to limit the search range if the | +| | approximate curvature is known a priori. | ++-------------------+------------------------------------------------------------+ +| ``nsmooth`` | Window length for the Savitzky–Golay smoothing filter | +| | applied to the summed power profile. Must be an odd | +| | integer ≥ 3. | ++-------------------+------------------------------------------------------------+ + +**Fitting each side of the arc separately** + +When the arc is clearly asymmetric (e.g. due to anisotropic scattering), +use ``asymm=True`` to fit the left and right wings independently:: + + dyn.fit_arc(asymm=True, plot=True) + print(dyn.eta_left, dyn.eta_right) + +**Troubleshooting** + +- *No clear arc visible*: try increasing ``startbin`` or adjusting + ``low_power_diff``/``high_power_diff`` to widen the fitting region. +- *Fit returns forward parabola*: the peak power sits at zero curvature. + Try increasing ``constraint[0]`` or reducing ``delmax``. +- *Error is very large*: the arc is broad or noisy. Try reducing + ``nsmooth`` or using ``log_parabola=True``. + + +ACF analysis +------------ + +**Computing the ACF** + +:: + + dyn.calc_acf() + dyn.plot_acf() + +The 2-D ACF is stored in ``dyn.acf`` with lag axes ``dyn.tau_samples`` +(time lags in minutes) and ``dyn.dnu_samples`` (frequency lags in MHz). + +**Fitting the ACF** + +Fit a 1-D model along the time and frequency axes simultaneously:: + + dyn.fit_acf() + print(f"Scintillation timescale tau = {dyn.tau:.1f} +/- {dyn.tauerr:.1f} s") + print(f"Decorrelation bandwidth dnu = {dyn.dnu:.3f} +/- {dyn.dnuerr:.3f} MHz") + +For a full 2-D fit that also measures the phase gradient (drift) in the +ACF:: + + dyn.fit_acf(method='2d') + print(dyn.acf_tilt, dyn.phasegrad) + +**Extracting scintillation parameters** + +Once ``tau`` and ``dnu`` are measured, derive the scintillation velocity:: + + from scintools.scint_utils import scint_velocity + viss = scint_velocity(params=None, dnu=dyn.dnu, tau=dyn.tau, + freq=dyn.freq) + print(f"Scintillation velocity Viss = {viss:.1f} km/s") + + +Simulations +----------- + +**Creating a simulated dynamic spectrum** + +Simulate a dynamic spectrum using the electromagnetic wave propagation code +of Coles et al. (2010):: + + from scintools.scint_sim import Simulation + from scintools.dynspec import Dynspec + + sim = Simulation(mb2=2, rf=1, ds=0.01, alpha=5/3, ar=1, psi=0, + ns=256, nf=256, dlam=0.25, freq=1400, dt=30) + dyn = Dynspec(dyn=sim, verbose=True) + dyn.plot_dyn() + +**Inspecting the theoretical arc curvature** + +The simulation computes the theoretical arc curvature automatically:: + + print(f"Theoretical eta = {sim.eta:.4f} s^3") + + dyn.calc_sspec() + dyn.fit_arc(plot=True) + print(f"Measured eta = {dyn.eta:.4f} +/- {dyn.etaerr:.4f} s^3") + +**Parameter studies** + +Run multiple simulations to study how arc curvature varies with screen +geometry:: + + import numpy as np + from scintools.scint_sim import Simulation + from scintools.dynspec import Dynspec + + screen_distances = np.linspace(0.1, 0.9, 9) + eta_values = [] + for s in screen_distances: + sim = Simulation(mb2=2, rf=1, ds=0.01, seed=42) + dyn = Dynspec(dyn=sim, verbose=False) + dyn.calc_sspec() + dyn.fit_arc(plot=False, display=False) + eta_values.append(dyn.eta) + + print("Screen distance vs measured curvature:") + for s, eta in zip(screen_distances, eta_values): + print(f" s={s:.1f} eta={eta:.4f}") + +**Common pitfalls** + +- The ``seed`` parameter controls the random screen realisation. Set it to + a fixed value for reproducible results. +- For lamsteps simulations, pass ``lamsteps=True`` to both + :class:`~scintools.scint_sim.Simulation` and + :class:`~scintools.dynspec.Dynspec`. diff --git a/docs/source/index.rst b/docs/source/index.rst index d3e0349..7bb5d6b 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -11,6 +11,8 @@ Welcome to the documentation for Scintools! :caption: Contents: dynspec + + examples simulation diff --git a/scintools/__init__.py b/scintools/__init__.py index 26fa278..8e2beb7 100644 --- a/scintools/__init__.py +++ b/scintools/__init__.py @@ -1,6 +1,4 @@ # -*- coding: utf-8 -*- -from __future__ import (absolute_import, division, - print_function, unicode_literals) __author__ = """Daniel John Reardon""" __email__ = 'dreardon@swin.edu.au' diff --git a/scintools/dynspec.py b/scintools/dynspec.py index 92311d4..afc5d6c 100755 --- a/scintools/dynspec.py +++ b/scintools/dynspec.py @@ -6,9 +6,6 @@ Dynamic spectrum class """ -from __future__ import (absolute_import, division, - print_function, unicode_literals) - import time import os from os.path import split diff --git a/scintools/scint_models.py b/scintools/scint_models.py index 23ec251..e2fe78f 100755 --- a/scintools/scint_models.py +++ b/scintools/scint_models.py @@ -19,8 +19,6 @@ Some functions use additional inputs """ -from __future__ import (absolute_import, division, - print_function, unicode_literals) import numpy as np from scintools.scint_sim import ACF from lmfit import Minimizer @@ -29,6 +27,53 @@ def fitter(model, params, args, mcmc=False, pos=None, nwalkers=100, steps=1000, burn=0.2, progress=True, workers=1, nan_policy='raise', max_nfev=None, thin=10, is_weighted=True): + """ + Fit a model to data using least-squares or MCMC via :mod:`lmfit`. + + Parameters + ---------- + model : callable + Residual function with signature ``model(params, *args)`` that + returns an array of weighted residuals ``(ydata - model) * weights``. + params : lmfit.Parameters + Initial parameter values and bounds. + args : tuple + Extra positional arguments passed to ``model`` after ``params``. + mcmc : bool, optional + If ``True``, run MCMC sampling with ``emcee`` instead of + least-squares. Default is ``False``. + pos : array_like or None, optional + Initial positions of the MCMC walkers. Passed to + :meth:`lmfit.Minimizer.emcee`. Default is ``None``. + nwalkers : int, optional + Number of MCMC ensemble walkers. Default is ``100``. + steps : int, optional + Total number of MCMC steps per walker. Default is ``1000``. + burn : float, optional + Fraction of steps to discard as burn-in (between 0 and 1). + Default is ``0.2``. + progress : bool, optional + Show a progress bar during MCMC sampling. Default is ``True``. + workers : int, optional + Number of parallel workers for MCMC. Default is ``1``. + nan_policy : str, optional + How to handle NaN residuals in least-squares mode. Passed to + :class:`lmfit.Minimizer`. Default is ``'raise'``. + max_nfev : int or None, optional + Maximum number of function evaluations for least-squares. Default + is ``None`` (no limit). + thin : int, optional + Thinning factor for MCMC chains. Default is ``10``. + is_weighted : bool, optional + Whether the residuals are weighted. Passed to + :meth:`lmfit.Minimizer.emcee`. Default is ``True``. + + Returns + ------- + results : lmfit.MinimizerResult + Fit results including best-fit parameter values, uncertainties, and + (for MCMC) the sampler flatchain. + """ # Do fit if mcmc: @@ -47,6 +92,24 @@ def fitter(model, params, args, mcmc=False, pos=None, nwalkers=100, def powerspectrum_model(params, xdata, ydata): + """ + Power-law model for the power spectrum with a white-noise floor. + + Parameters + ---------- + params : lmfit.Parameters + Must contain ``amp`` (power-law amplitude), ``wn`` (white-noise + level), and ``alpha`` (spectral index). + xdata : array_like + Frequency (or wavenumber) values. + ydata : array_like + Observed power spectrum values. + + Returns + ------- + residuals : numpy.ndarray + Array ``ydata - model``, where ``model = wn + amp * xdata**alpha``. + """ parvals = params.valuesdict() @@ -351,9 +414,50 @@ def arc_curvature(params, ydata, weights, true_anomaly, vearth_ra, vearth_dec, mjd=None, model_only=False, return_veff=False): """ - arc curvature model + Thin-screen arc curvature model for scintillation arcs. - ydata: arc curvature + Computes the expected parabolic arc curvature :math:`\\eta` from the + effective velocity at the scattering screen, optionally returning + weighted residuals for use with a least-squares fitter. + + Parameters + ---------- + params : lmfit.Parameters or dict + Must contain ``d`` (pulsar distance in kpc) and ``s`` (fractional + screen distance). May also contain ``zeta`` / ``vism_zeta`` (for + anisotropic screens) or ``vism_ra`` / ``vism_dec`` (ISM velocity + components in km/s). + ydata : array_like + Observed arc curvature values (1/(m·mHz²)). + weights : array_like or None + Weights for the residuals. Pass ``None`` for uniform weighting. + true_anomaly : array_like + Orbital true anomaly in radians at each epoch. + vearth_ra : array_like + Earth velocity in the RA direction (km/s) at each epoch. + vearth_dec : array_like + Earth velocity in the Dec direction (km/s) at each epoch. + mjd : array_like or None, optional + Barycentric MJD at each epoch. Required when ``OMDOT`` is in + ``params``. Default is ``None``. + model_only : bool, optional + If ``True``, return only the model curvature array rather than + weighted residuals. Default is ``False``. + return_veff : bool, optional + If ``True`` and ``model_only`` is ``True``, also return the RA and + Dec components of the effective velocity. Default is ``False``. + + Returns + ------- + residuals : numpy.ndarray + Weighted residuals ``(ydata - model) * weights``. Returned when + ``model_only=False``. + model : numpy.ndarray + Arc curvature model values. Returned when ``model_only=True`` and + ``return_veff=False``. + model, veff_ra, veff_dec : tuple + Arc curvature model and effective velocity components. Returned when + ``model_only=True`` and ``return_veff=True``. """ # ensure dimensionality of arrays makes sense @@ -504,8 +608,54 @@ def veff_thin_screen(params, ydata, weights, true_anomaly, def effective_velocity_annual(params, true_anomaly, vearth_ra, vearth_dec, mjd=None): """ - Effective velocity with annual and pulsar terms - Note: Does NOT include IISM velocity, but returns veff in IISM frame + Compute the annual and orbital effective velocity at the scattering screen. + + Accounts for Earth's orbital motion, the pulsar's proper motion, and (if + Keplerian parameters are present) the pulsar's orbital motion. The ISM + velocity is *not* included here; it must be subtracted by the calling + function. + + Parameters + ---------- + params : lmfit.Parameters or dict + May contain Keplerian orbital parameters (``A1``, ``PB``, ``ECC``, + ``OM``, ``KOM``, ``KIN``/``COSI``/``SINI``) and proper-motion + parameters (``PMRA``, ``PMDEC``). Must also contain ``s`` + (fractional screen distance, 0–1) and ``d`` (pulsar distance in + kpc). + true_anomaly : array_like + Orbital true anomaly in radians at each epoch. + vearth_ra : array_like + Earth velocity in the RA direction (km/s) at each epoch. + vearth_dec : array_like + Earth velocity in the Dec direction (km/s) at each epoch. + mjd : array_like or None, optional + Barycentric MJD at each epoch. Required when ``OMDOT`` is in + ``params``. Default is ``None``. + + Returns + ------- + veff_ra : numpy.ndarray + Effective velocity in the RA direction (km/s). + veff_dec : numpy.ndarray + Effective velocity in the Dec direction (km/s). + vp_ra : numpy.ndarray + Pulsar orbital velocity in the RA direction (km/s). + vp_dec : numpy.ndarray + Pulsar orbital velocity in the Dec direction (km/s). + + Notes + ----- + The effective velocity is defined as: + + .. math:: + + \\mathbf{v}_{\\mathrm{eff}} = s\\,\\mathbf{v}_{\\oplus} + + (1-s)\\,(\\mathbf{v}_{p} + \\boldsymbol{\\mu}_p) + + where :math:`s` is the fractional screen distance, :math:`\\mathbf{v}_{\\oplus}` + is Earth's velocity, and :math:`\\mathbf{v}_p + \\boldsymbol{\\mu}_p` is the + pulsar's velocity (orbital + proper motion). """ # Define some constants v_c = 299792.458 # km/s diff --git a/scintools/scint_sim.py b/scintools/scint_sim.py index 38ac9cb..a0dea73 100755 --- a/scintools/scint_sim.py +++ b/scintools/scint_sim.py @@ -6,9 +6,6 @@ Simulate scintillation. Based on original MATLAB code by Coles et al. (2010) """ -from __future__ import (absolute_import, division, - print_function, unicode_literals) - import numpy as np from numpy import random from numpy.random import randn @@ -28,20 +25,97 @@ def __init__(self, mb2=2, rf=1, ds=0.01, alpha=5/3, ar=1, psi=0, verbose=False, freq=1400, dt=30, mjd=60000, nsub=None, efield=False, noise=None): """ - Electromagnetic simulator based on original code by Coles et al. (2010) - - mb2: Max Born parameter for strength of scattering - rf: Fresnel scale - ds (or dx,dy): Spatial step sizes with respect to rf - alpha: Structure function exponent (Kolmogorov = 5/3) - ar: Anisotropy axial ratio - psi: Anisotropy orientation - inner: Inner scale w.r.t rf - should generally be smaller than ds - ns (or nx,ny): Number of spatial steps - nf: Number of frequency steps. - dlam: Fractional bandwidth relative to centre frequency - lamsteps: Boolean to choose whether steps in lambda or freq - seed: Seed number, or use "-1" to shuffle + Electromagnetic simulator based on original code by Coles et al. (2010). + + Generates a simulated dynamic spectrum by propagating a plane wave + through a thin, phase-modulating scattering screen described by a + power-law structure function. + + Parameters + ---------- + mb2 : float, optional + Maximum Born parameter, controlling the strength of scattering. + Default is 2. + rf : float, optional + Fresnel scale (in arbitrary spatial units). Default is 1. + ds : float, optional + Spatial step size relative to ``rf``. Used for both axes unless + ``dx`` / ``dy`` are set separately. Default is 0.01. + alpha : float, optional + Structure function exponent. Kolmogorov turbulence gives + ``alpha = 5/3``. Default is ``5/3``. + ar : float, optional + Anisotropy axial ratio of the scattering screen. Default is 1 + (isotropic). + psi : float, optional + Orientation angle of the anisotropy (degrees). Default is 0. + inner : float, optional + Inner scale of the turbulence relative to ``rf``. Should be + smaller than ``ds``. Default is 0.001. + ns : int, optional + Number of spatial steps along each axis (used when ``nx`` / + ``ny`` are not set separately). Default is 256. + nf : int, optional + Number of frequency steps in the dynamic spectrum. Default is + 256. + dlam : float, optional + Fractional bandwidth relative to the centre frequency. Default + is 0.25. + lamsteps : bool, optional + If ``True``, use equal steps in wavelength rather than + frequency. Default is ``False``. + seed : int or None, optional + Random-number seed for the screen realisation. Use ``-1`` to + shuffle (random seed). Default is ``None``. + nx : int or None, optional + Number of spatial steps in the x-direction. Overrides ``ns`` + when set. Default is ``None``. + ny : int or None, optional + Number of spatial steps in the y-direction. Overrides ``ns`` + when set. Default is ``None``. + dx : float or None, optional + Spatial step size in x relative to ``rf``. Overrides ``ds`` + when set. Default is ``None``. + dy : float or None, optional + Spatial step size in y relative to ``rf``. Overrides ``ds`` + when set. Default is ``None``. + plot : bool, optional + Plot all intermediate results after the simulation. Default is + ``False``. + verbose : bool, optional + Print progress messages during the simulation. Default is + ``False``. + freq : float, optional + Centre frequency of the observation in MHz. Default is 1400. + dt : float, optional + Time step between sub-integrations in seconds. Default is 30. + mjd : float, optional + Reference Modified Julian Date of the observation. Default is + 60000. + nsub : int or None, optional + Number of sub-integrations to keep. If ``None`` (default) all + available sub-integrations are used. + efield : bool, optional + If ``True``, use the real part of the electric field for the + dynamic spectrum instead of the intensity. Default is ``False``. + noise : float or None, optional + RMS noise level to add to the simulated intensity. If ``None`` + (default) no noise is added. + + Notes + ----- + After initialisation the object exposes the following key attributes: + + - ``dyn`` : 2-D intensity array of shape ``(nchan, nsub)``. + - ``freqs`` : Centre frequencies of each channel (MHz). + - ``times`` : Sub-integration mid-times (seconds from start). + - ``eta`` : Theoretical arc curvature in physical units (s³). + - ``betaeta`` : Theoretical arc curvature for wavelength-scaled + secondary spectrum. + + References + ---------- + Coles, W. A., et al. (2010), ApJ, 717, 1206. """ self.mb2 = mb2 diff --git a/scintools/scint_utils.py b/scintools/scint_utils.py index fc0ad50..7a7e635 100755 --- a/scintools/scint_utils.py +++ b/scintools/scint_utils.py @@ -6,9 +6,6 @@ Useful functions for scintools """ -from __future__ import (absolute_import, division, - print_function, unicode_literals) - import numpy as np import os import sys @@ -66,7 +63,27 @@ def clean_archive(archive, template=None, bandwagon=0.99, channel_threshold=5, def autocorr(arr): """ - Do a slow calculation of the 2d autocorrelation of an input masked array + Compute the 2D autocorrelation of a masked array. + + Uses a slow direct calculation that properly handles masked (missing) data. + + Parameters + ---------- + arr : numpy.ma.MaskedArray + 2D masked array to autocorrelate, shape (nr, nc). + + Returns + ------- + autocorr : numpy.ndarray + Normalised 2D autocorrelation of shape (2*nr, 2*nc). The zero-lag + peak is located at index (nr, nc) and the array is normalised so + that its maximum value is 1. + + Notes + ----- + This is a direct O(n^4) implementation intended for small arrays where + masked pixels must be excluded from each lag sum individually. For large + unmasked arrays, FFT-based methods are much faster. """ mean = np.ma.mean(arr) std = np.ma.std(arr) @@ -86,14 +103,42 @@ def autocorr(arr): def is_valid(array): """ - Returns boolean array of values that are finite an not nan + Return a boolean array of elements that are finite and not NaN. + + Parameters + ---------- + array : array_like + Input array to test. + + Returns + ------- + mask : numpy.ndarray of bool + Boolean array of the same shape as ``array``, where ``True`` indicates + that the corresponding element is both finite (not inf) and not NaN. + + Examples + -------- + >>> import numpy as np + >>> from scintools.scint_utils import is_valid + >>> is_valid(np.array([1.0, np.nan, np.inf, -np.inf, 2.0])) + array([ True, False, False, False, True]) """ return np.isfinite(array)*(~np.isnan(array)) def read_dynlist(file_path): """ - Reads list of dynamic spectra filenames from path + Read a list of dynamic spectrum filenames from a plain-text file. + + Parameters + ---------- + file_path : str + Path to the text file containing one filename per line. + + Returns + ------- + dynfiles : list of str + List of filenames read from ``file_path``. """ with open(file_path) as file: dynfiles = file.read().splitlines() @@ -102,7 +147,21 @@ def read_dynlist(file_path): def write_results(filename, dyn=None): """ - Appends dynamic spectrum information and parameters of interest to file + Append dynamic spectrum parameters to a CSV results file. + + Writes a header line on the first call (when the file is empty), then + appends one data row per call. The set of columns grows dynamically + depending on which optional attributes exist on ``dyn``. + + Parameters + ---------- + filename : str + Path to the output CSV file. Created if it does not exist. + dyn : Dynspec, optional + Dynamic spectrum object whose attributes are written. Expected to + have at minimum: ``name``, ``mjd``, ``freq``, ``bw``, ``tobs``, + ``dt``, ``df``. Optional attributes (``tau``, ``dnu``, ``eta``, …) + are included when present. """ header = "name,mjd,freq,bw,tobs,dt,df" @@ -204,7 +263,19 @@ def write_results(filename, dyn=None): def read_results(filename): """ - Reads a CSV results file written by write_results() + Read a CSV results file written by :func:`write_results`. + + Parameters + ---------- + filename : str + Path to the CSV file to read. + + Returns + ------- + param_dict : dict + Dictionary mapping each column header to a list of string values. + Use :func:`float_array_from_dict` to convert individual columns to + numeric arrays. """ csv_data = open(filename, 'r') @@ -233,7 +304,18 @@ def search_and_replace(filename, search, replace): def cov_to_corr(cov): """ - Calculate correlation matrix from covariance + Convert a covariance matrix to a correlation matrix. + + Parameters + ---------- + cov : numpy.ndarray, shape (n, n) + Symmetric positive-semidefinite covariance matrix. + + Returns + ------- + corr : numpy.ndarray, shape (n, n) + Correlation matrix derived from ``cov``. Elements where the + corresponding covariance is zero are set to zero. """ std = np.sqrt(np.diag(cov)) outer_std = np.outer(std, std) @@ -244,7 +326,24 @@ def cov_to_corr(cov): def float_array_from_dict(dictionary, key): """ - Convert an array stored in dictionary to a numpy array + Extract a numeric array from a string-valued results dictionary. + + Replaces ``'None'`` string entries with ``'nan'`` before converting, so + that missing values become ``numpy.nan`` rather than raising an error. + + Parameters + ---------- + dictionary : dict + Dictionary as returned by :func:`read_results`, where values are + lists of strings. + key : str + Column name to extract and convert. + + Returns + ------- + arr : numpy.ndarray of float + 1-D float array of the requested column, squeezed to remove + length-1 dimensions. """ ind = np.argwhere(np.array(dictionary[key]) == 'None').ravel() @@ -269,8 +368,22 @@ def save_fits(filename, dyn): def difference(x): """ - unlike np.diff, computes differences between centres of elements in x, - returns numpy array same size as x + Compute centre-to-centre differences for an array, returning the same size. + + Unlike :func:`numpy.diff`, this function returns an array of the same + length as the input by using one-sided differences at the boundaries. + + Parameters + ---------- + x : array_like + 1-D input array of values. + + Returns + ------- + dx : numpy.ndarray + Array of differences the same length as ``x``. Interior elements + are ``(x[i+1] - x[i-1]) / 2``; boundary elements use the adjacent + pair only. """ dx = [] for i in range(0, len(x)): @@ -285,8 +398,27 @@ def difference(x): def get_ssb_delay(mjds, raj, decj, message=True): """ - Get Romer delay to Solar System Barycentre (SSB) for correction of site - arrival times to barycentric. + Compute the Solar System Barycentre (SSB) Römer delay for each epoch. + + Calculates the light-travel-time delay between a ground-based observatory + and the SSB along the line of sight to the pulsar. + + Parameters + ---------- + mjds : array_like + Modified Julian Dates (MJD) of the observations (barycentric). + raj : str + Right ascension in ``HH:MM:SS.S`` format. + decj : str + Declination in ``±DD:MM:SS.S`` format. + message : bool, optional + If ``True`` (default), print a reminder that the returned delays + should be *added* to site arrival times. + + Returns + ------- + delays : numpy.ndarray + Römer delays in seconds, one per entry in ``mjds``. """ from astropy.constants import au, c @@ -348,8 +480,32 @@ def make_lsr(d, raj, decj, pmra, pmdec, vr=0): def get_earth_velocity(mjds, raj, decj, radial=False): """ - Calculates the component of Earth's velocity transverse to the line of - sight, in RA and DEC. Optionally returns the radial velocity + Compute the transverse Earth velocity projected onto the sky. + + Returns the component of Earth's barycentric velocity projected onto the + RA and Dec axes (and optionally the radial component) at each epoch. + + Parameters + ---------- + mjds : array_like + Modified Julian Dates (MJD) at which to evaluate the velocity. + raj : str + Right ascension in ``HH:MM:SS.S`` format. + decj : str + Declination in ``±DD:MM:SS.S`` format. + radial : bool, optional + If ``True``, also return the velocity component along the line of + sight. Default is ``False``. + + Returns + ------- + vearth_ra : numpy.ndarray + Earth velocity component in the RA direction (km/s). + vearth_dec : numpy.ndarray + Earth velocity component in the Dec direction (km/s). + vearth_radial : numpy.ndarray + Earth velocity along the line of sight (km/s). Only returned when + ``radial=True``. """ from astropy.time import Time @@ -397,7 +553,21 @@ def get_earth_velocity(mjds, raj, decj, radial=False): def read_par(parfile): """ - Reads a par file and return a dictionary of parameter names and values + Read a TEMPO/TEMPO2 ``.par`` file and return a parameter dictionary. + + Parameters + ---------- + parfile : str + Path to the pulsar parameter file. + + Returns + ------- + par : dict + Dictionary of parameter names to values. Numeric parameters are + stored as ``int`` or ``float``; string parameters remain as ``str``. + Fit errors are stored under ``_ERR`` and the numeric type + under ``_TYPE`` (``'d'`` for int, ``'f'`` or ``'e'`` for + float, ``'s'`` for string). """ par = {} @@ -452,7 +622,17 @@ def read_par(parfile): def mjd_to_year(mjd): """ - converts mjd to year + Convert a Modified Julian Date (MJD) to a Besselian year. + + Parameters + ---------- + mjd : float or array_like + Modified Julian Date(s) to convert. + + Returns + ------- + year : float or numpy.ndarray + Besselian year(s) corresponding to ``mjd``. """ t = Time(mjd, format='mjd') yrs = t.byear # observation year @@ -461,7 +641,20 @@ def mjd_to_year(mjd): def find_nearest(arr, val): """ - Returns the index of an array (arr) that is nearest to value (val) + Return the index of the element in ``arr`` closest to ``val``. + + Parameters + ---------- + arr : array_like + 1-D array to search. + val : float + Target value. + + Returns + ------- + ind : int + Index of the element in ``arr`` with the smallest absolute difference + from ``val``. """ arr = np.asarray(arr) ind = np.argmin(np.abs(arr - val)) @@ -479,10 +672,30 @@ def longest_run_of_zeros(arr): def pars_to_params(pars, params=None): """ - Converts a dictionary of par file parameters from read_par() to an - lmfit Parameters() class to use in models + Convert a par-file dictionary to an :class:`lmfit.Parameters` object. + + Numeric parameters are added with ``vary=False``. String parameters and + the special ``RAJ``/``RA`` and ``DECJ`` entries are handled separately. - By default, parameters are not varied + Parameters + ---------- + pars : dict + Parameter dictionary as returned by :func:`read_par`. + params : lmfit.Parameters, optional + Existing :class:`lmfit.Parameters` instance to append to. If + ``None`` (default), a new instance is created. + + Returns + ------- + params : lmfit.Parameters + Updated (or newly created) :class:`lmfit.Parameters` object + containing all numeric parameters from ``pars``. + + Notes + ----- + By default, all parameters are fixed (``vary=False``). Enable fitting + by setting individual parameter ``vary`` attributes after calling this + function. """ from lmfit import Parameters @@ -508,8 +721,24 @@ def pars_to_params(pars, params=None): def get_true_anomaly(mjds, pars): """ - Calculates true anomalies for an array of barycentric MJDs and a parameter - dictionary + Calculate orbital true anomalies from barycentric MJDs. + + Solves Kepler's equation to obtain the eccentric anomaly and converts to + true anomaly. Supports both T0/ECC and TASC/EPS1/EPS2 parameterisations. + + Parameters + ---------- + mjds : array_like + Barycentric Modified Julian Dates at which to evaluate the anomaly. + pars : dict + Pulsar timing parameter dictionary (from :func:`read_par`) containing + at minimum ``PB`` and either ``T0``+``ECC`` or ``TASC``+``EPS1``+ + ``EPS2``. Optional keys: ``PBDOT``. + + Returns + ------- + U : numpy.ndarray + True anomaly in radians, in the range ``[0, 2π)``. """ if 'TASC' in pars.keys(): @@ -556,8 +785,24 @@ def get_true_anomaly(mjds, pars): def get_binphase(mjds, pars): """ - Calculates binary phase for an array of barycentric MJDs and a parameter - dictionary + Calculate the orbital binary phase from barycentric MJDs. + + The binary phase is defined as the true anomaly plus the argument of + periastron, ``U + OM``. + + Parameters + ---------- + mjds : array_like + Barycentric Modified Julian Dates. + pars : dict + Pulsar timing parameter dictionary (from :func:`read_par`). If + ``TASC`` is present the orbit is assumed circular (``OM = 0``); + otherwise ``OM`` (and optionally ``OMDOT``) are used. + + Returns + ------- + phase : numpy.ndarray + Binary phase in radians. """ U = get_true_anomaly(mjds, pars) From dec5da05eac213d81eed640fdad4e2f889c6c5ce Mon Sep 17 00:00:00 2001 From: Daniel Reardon <40694862+danielreardon@users.noreply.github.com> Date: Wed, 4 Mar 2026 07:16:31 +1100 Subject: [PATCH 3/4] Correct parameter description for mb2 in scint_sim.py --- scintools/scint_sim.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scintools/scint_sim.py b/scintools/scint_sim.py index a0dea73..c6f519a 100755 --- a/scintools/scint_sim.py +++ b/scintools/scint_sim.py @@ -34,7 +34,7 @@ def __init__(self, mb2=2, rf=1, ds=0.01, alpha=5/3, ar=1, psi=0, Parameters ---------- mb2 : float, optional - Maximum Born parameter, controlling the strength of scattering. + Max Born parameter, controlling the strength of scattering. Default is 2. rf : float, optional Fresnel scale (in arbitrary spatial units). Default is 1. From efc62511483b6b5f146b94c65cede52c63b7c2f7 Mon Sep 17 00:00:00 2001 From: Daniel Reardon <40694862+danielreardon@users.noreply.github.com> Date: Wed, 4 Mar 2026 07:23:48 +1100 Subject: [PATCH 4/4] Update documentation for dynamic spectra processing --- docs/source/dynspec.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/source/dynspec.rst b/docs/source/dynspec.rst index 9bb3818..2d9eb45 100644 --- a/docs/source/dynspec.rst +++ b/docs/source/dynspec.rst @@ -21,7 +21,7 @@ provides methods to: timescales and bandwidths. - Visualise all intermediate and final data products. -See :doc:`examples` for complete, runnable worked examples. +See :doc:`examples` for runnable worked examples. Importing dynamic spectra ------------------------- @@ -60,7 +60,7 @@ several ways:: dyn.zap() # interactive RFI flagging dyn.refill() # interpolate over flagged data dyn.crop_dyn(tmin, tmax, fmin, fmax) # crop in time and/or frequency - dyn.scale_dyn() # normalise each sub-integration + dyn.scale_dyn() # resample the dynamic spectrum to a new scale Secondary spectra -----------------