blscint API reference

blscint.ne2001 module

blscint.sample_t_d module

blscint.sample_t_d.mcmillan_rho_bulge(R, z)

blscint.sample_t_d.mcmillan_rho_thin(R, z)

blscint.sample_t_d.mcmillan_rho_thick(R, z)

blscint.sample_t_d.mcmillan_rho_tot(R, z)

blscint.sample_t_d.carroll_ostlie_n(R, z)

blscint.sample_t_d.coverage(t_ds, start=None, stop=None)

Get fractional coverage of scintillation timescales for a specified range.

blscint.sample_t_d.central(t_ds, p)

Get central fraction p of array.

blscint.sample_t_d.visualize_mc(t_ds)

blscint.sample_t_d.transition_freqs(l, b, d=(0.001, 20), d_steps=1000)

Get transition from strong to weak scattering.

blscint.sample_t_d.min_d_ss(l, b, d=(0.001, 20), f=(4, 8), delta_d=0.01)

Get min distance for strong scattering.

class blscint.sample_t_d.NESampler(l, b, d=(0.01, 20), delta_d=0.01)

Bases: object

Class for sampling scintillation estimates from NE2001 electron density model.

save_pickle(filename)

Save entire sampler (including NE2001 calculations) as a pickled file (.pickle).

classmethod load_pickle(filename)

Load sampler object from a pickled file (.pickle), created with NESampler.save_pickle.

tau_d(nu_d, C1=1.16)

Temporal broadening. C1 is 1.53 for a medium with a square-law structure function. NE2001 paper.

nu_sb(t_d, C2=2.02)

Spectral broadening.

sample(n=1000, f=(4, 8), v=(5, 100), d=None, d_sampling_type='mcmillan', weight_by_flux=False, galcen_distance=8.5, scint_regime='moderate', verbose=False)

Sample frequencies, transverse velocities, and scale base model values appropriately.

Sampling type can be: ‘uniform’, ‘mcmillan’, ‘carrollostlie’. McMillan model uses 8.21 kpc to the galactic center, but NE2001 uses 8.5 kpc.

Set weight_by_flux=True to account for inverse square law.

class blscint.sample_t_d.NEData(t_d, nu_d, d=array([], dtype=float64), f=array([], dtype=float64), v=array([], dtype=float64))

Bases: object

property tau_d

property nu_sb

property data_dict

mask(mask)

report()

diagnostic_plot()

plot(num_bins=25, masks=None, labels=None)

Plot 2x2 grid of histograms of sampled scintillation properties.

Optionally, specify boolean masks and labels to plot sub-histograms within each figure, to get a sense of the sampling breakdown.

plot_t_d()

blscint.gen_arta module

Generate scintillated signals matching Gaussian pulse profiles and exponential intensity distributions using the autogregressive to anything (ARTA) algorithm.

Cario & Nelson 1996: https://www.sciencedirect.com/science/article/pii/016763779600017X

Scintillation on narrowband signals references:

Cordes & Lazio 1991: http://articles.adsabs.harvard.edu/pdf/1991ApJ…376..123C

Cordes, Lazio, & Sagan 1997: https://iopscience.iop.org/article/10.1086/304620/pdf

blscint.gen_arta.get_rho(t_d, dt, p, pow=1.6666666666666667)

Get target autocorrelation guesses for ARTA with scintillation timescale t_d and time resolution dt, up to lag p. Modeled as a power law with exponent 5/3 or 2.

Parameters:

• t_d (float) – Scintillation timescale (s)

• dt (float) – Time resolution (s)

• p (int) – Number of lags to calculate

• pow (float, optional) – Exponent for ACF fit, either 5/3 or 2 (arising from phase structure function)

Returns:

r – Array of autocorrelations, starting with lag 1

Return type:

np.ndarray

blscint.gen_arta.psi(r)

Return covariance matrix for initial multivariate normal distribution.

Parameters:

r (np.ndarray) – Array of autocorrelation guesses, starting with lag 1

Returns:

cov_mat – Covariance matrix

Return type:

np.ndarray

blscint.gen_arta.build_Z(r, T)

Build full baseline Z array.

Parameters:

• r (np.ndarray) – Array of autocorrelation guesses, starting with lag 1

• T (int) – Final length of array Z, should be greater than p

Returns:

Z – Array of Z values, as in ARTA

Return type:

np.ndarray

blscint.gen_arta.inv_exp_cdf(x, rate=1)

Inverse exponential distribution CDF.

Parameters:

• x (float, np.ndarray) – Input value(s) from [0, 1)

• rate (float) – Rate parameter lambda

Returns:

v – Output of inverse CDF

Return type:

float, np.ndarray

blscint.gen_arta.inv_levy_cdf(x, loc=0, scale=1)

blscint.gen_arta.get_Y(Z, dist='exp')

Get final values specific to an overall exponential distribution, normalized to mean of 1.

Parameters:

Z (np.ndarray) – Array of Z values, as in ARTA

Returns:

Y – Final synthetic scintillated time series (Y values)

Return type:

np.ndarray

blscint.gen_arta.get_ts_arta(t_d, dt, num_samples, p=2, pow=1.6666666666666667, dist='exp')

Produce time series data via an ARTA process.

Parameters:

• t_d (float) – Scintillation timescale (s)

• dt (float) – Time resolution (s)

• num_samples (int) – Number of synthetic samples to produce

• p (int, optional) – Number of lags to calculate

• pow (float, optional) – Exponent for ACF fit, either 5/3 or 2 (arising from phase structure function)

Returns:

Y – Final synthetic scintillated time series (Y values)

Return type:

np.ndarray

blscint.gen_pdf module

Generate scintillated signals matching Gaussian pulse profiles and exponential intensity distributions using theoretical PDFs.

Scintillation on narrowband signals references:

Cordes & Lazio 1991: http://articles.adsabs.harvard.edu/pdf/1991ApJ…376..123C

Cordes, Lazio, & Sagan 1997: https://iopscience.iop.org/article/10.1086/304620/pdf

blscint.gen_pdf.find_nearest(arr, val)

Return index of closest value.

Parameters:

• arr (np.ndarray) – Input array

• val (float) – Target value

Returns:

idx – Closest index

Return type:

int

blscint.gen_pdf.bpdf(g1, g2, ac)

Calculate joint probability of g1, g2 separated by a temporal auto-correlation value of ac, according to Eq. B12 of Cordes, Lazio, & Sagan 1997.

Parameters:

• g1 (float) – Signal gain 1

• g2 (float) – Signal gain 2

• ac (float) – Autocorrelation value in [0, 1)

Returns:

f_2g – Joint probability

Return type:

float

blscint.gen_pdf.get_ts_pdf(t_d, dt, num_samples, max_g=5, steps=1000)

Produce time series data via bivariate pdf for the gain. With a maximum gain max_g and number of potential gain levels within [0, max_g].

Parameters:

• t_d (float) – Scintillation timescale (s)

• dt (float) – Time resolution (s)

• num_samples (int) – Number of synthetic samples to produce

• max_g (float, optional) – Maximum possible gain (for computation)

• steps (int, optional) – Number of possible gain levels within [0, max_g] (for computation)

Returns:

Y – Final synthetic scintillated time series (Y values)

Return type:

np.ndarray

blscint.analyze_obs module

blscint.dataframe module

blscint.dataframe.make_dataframe(dat_file)

Create Pandas dataframe from TurboSETI output.

Parameters:

dat_file (str) – TurboSETI .dat filename to format as a dataframe

Returns:

dataframe – Dataframe with extracted TurboSETI parameters

Return type:

DataFrame

blscint.dataframe.get_frame_params(fn)

Get frame resolution from a spectrogram file without loading the actual data.

Parameters:

fn (str) – .fil or .h5 filename

Returns:

params – Dictionary with tchans, df, dt

Return type:

dict

blscint.dataframe.centered_frame(fn, drift_rate, center_freq, fchans=256, frame_params=None)

Here, center_freq is in MHz.

blscint.dataframe.turbo_centered_frame(i, dataframe, fn=None, fchans=256, tchans=64, df=2.7939677238464355, dt=9.305762474666658, **kwargs)

Create Frame centered at a target signal from a TurboSETI-created .dat file. Does not remove drift – use Frame.dedrift separately to do so.

Parameters:

• i (int) – Signal index

• dataframe (DataFrame) – Pandas dataframe with TurboSETI parameters

• fn (str, optional) – Filename of datafile (unless filename is also in the dataframe, under ‘fn’)

• fchans (int, optional) – Number of frequency bins in target frame

• tchans (int, optional) – Number of time samples in data

• df (float, optional) – Frequency resolution (Hz)

• dt (float, optional) – Time resolution (Hz)

Returns:

frame – Generated frame

Return type:

stg.Frame

blscint.bounds module

blscint.bounds.plot_bounds(frame, l, r, use_db=False, cb=True, lw=2)

Plot frame data with overlaid bounding boxes.

Parameters:

• frame (stg.Frame) – Analysis frame

• l (int) – Left bound

• r (int) – Right bound

• use_db (bool, optional) – Option to convert intensities to dB

• cb (bool, optional) – Option to display colorbar

• lw (int, optional) – Line width in matplotlib

blscint.bounds.polyfit_bounds(spec, deg=1, snr_threshold=10)

Bounding box set by a polynomial fit to the background. Edges are set by where the fit intersects the data on either side of the central peak.

If the signal is broad compared to the spectrum or frame, the background fit will highly intersect or even model out the actual signal, hiding any peaks.

Degrees to try: 1, 7.

Parameters:

• spec (ndarray) – Intensity spectra

• deg (int, optional) – Degree of polynomial fit

• snr_threshold (int, float, optional) – Threshold for peak detection, in units of noise standard deviations

Returns:

• l (int) – Left bound

• r (int) – Right bound

• metadata (dict) – Dictionary with metadata related to peak-finding. Contains polynomial fit object, number of peaks, and detected peak information.

blscint.bounds.threshold_bounds(spec, half_width=3)

Create bounds based on integrated intensity on either side of the central peak. Threshold is set by ideal Gaussian profile, in units of standard deviations (sigma).

Parameters:

• spec (ndarray) – Intensity spectra

• half_width (float, optional) – Assuming a Gaussian signal profile, the half_width determines where to set the bounds, in units of sigma from the peak.

Returns:

• l (int) – Left bound

• r (int) – Right bound

• metadata (dict) – Dictionary with metadata. Contains noise mean and spectra maximum, which are used to normalize spec to the spectra maximum.

blscint.bounds.threshold_baseline_bounds(spec, p=0.01)

Create bounds based on integrated intensity on either side of the central peak, as a fraction of the peak integrated intensity. Uses a 1D fit to the noise baseline.

Parameters:

• spec (ndarray) – Intensity spectra

• p (float, optional) – Fraction of peak, used to set left and right bounds

Returns:

• l (int) – Left bound

• r (int) – Right bound

• metadata (dict) – Dictionary with metadata. Contains noise mean and spectra maximum, which are used to normalize spec to the spectra maximum.

blscint.bounds.snr_bounds(spec, snr=5)

Create bounds based on intensity attentuation on either side of the central peak, as an SNR threshold above the noise background. Uses a 1D fit to the noise baseline.

Parameters:

• spec (ndarray) – Intensity spectra

• snr (float, optional) – SNR threshold to set left and right bounds

Returns:

• l (int) – Left bound

• r (int) – Right bound

• metadata (dict) – Dictionary with metadata. Contains noise mean and spectra maximum, which are used to normalize spec to the spectra maximum.

blscint.bounds.gaussian_bounds(spec, half_width=3, peak_guess=None)

Create bounds based on a Gaussian fit to the central peak.

Parameters:

• spec (ndarray) – Intensity spectra

• half_width (float, optional) – Assuming a Gaussian signal profile, the half_width determines where to set the bounds, in units of sigma from the peak.

• peak_guess (int, optional) – Guess for peak index. Should normally be the center of the spectrum, but can have strange side effects.

Returns:

• l (int) – Left bound

• r (int) – Right bound

• metadata (dict) – Dictionary with metadata related to peak-finding. Contains Gaussian fit information.

blscint.bounds.clipped_bounds(spec, min_empty_bins=2)

Run sigma clip on spectrum, and find central peak.

This method reduces the spectrum to a mask based on sigma_clip, so it can be unreliable if the signal is broad compared to the size of the spectrum or frame.

Parameters:

• spec (ndarray) – Intensity spectra

• min_empty_bins (int, optional) – Minimum number of adjacent “empty”, or non-clipped, bins required to mark either bound with respect to the central peak.

Returns:

• l (int) – Left bound

• r (int) – Right bound

• metadata (dict) – Dictionary with metadata related to peak-finding. Contains detected peak information.

blscint.bounds.clipped_2D_bounds(frame, min_empty_bins=2, min_clipped=1, peak_prominence=4)

Run sigma clip on 2D data array, and find central peak above a certain number of clipped values along the time axis, per frequency bin.

This will return an IndexError if the signal passes outside of the frame (i.e. for very wide signals).

Note that this function accepts a frame instead of a 1D numpy array spectrum.

Parameters:

• frame (setigen.Frame) – Frame of intensity data

• min_empty_bins (int, optional) – Minimum number of adjacent “empty”, or non-clipped, bins required to mark either bound with respect to the central peak.

• min_clipped (int, optional) – Minimum number of pixels clipped by sigma_clip in the time direction in order for a frequency bin to be considered part of the signal.

• peak_prominence (int, optional) – Prominence in units of clipped pixels, for use in peak finding.

Returns:

• l (int) – Left bound

• r (int) – Right bound

• metadata (dict) – Dictionary with metadata related to peak-finding. Contains detected peak information.

blscint.bounds.boxcar_bounds(spec, window_sizes=None)

Bounding box set by boxcar filter. Simplified version of matched filtering. Default window sizes: powers of 2.

Parameters:

• spec (ndarray) – Intensity spectra

• window_sizes (list, array-like, optional) – List of boxcar window sizes to try. If None, automatically uses powers of 2 as window sizes.

Returns:

• l (int) – Left bound

• r (int) – Right bound

• metadata (dict) – Dictionary with metadata related to peak-finding. Contains peak SNR corresponding to the “best” boxcar filter.

blscint.frame_processing module

blscint.frame_processing.t_norm_frame(frame, as_data=None, divide_std=False)

Normalize frame by subtracting out noise background, along time axis. Additional option to divide out by the standard deviation of each individual spectrum.

Parameters:

• frame (stg.Frame) – Raw spectrogram frame

• as_data (stg.Frame, optional) – Use alternate frame to compute noise stats. If desired, use a more isolated region of time-frequency space for cleaner computation.

• divide_std (bool, optional) – Normalize each spectrum by dividing by its standard deviation

Returns:

n_frame – Normalized frame

Return type:

stg.Frame

blscint.frame_processing.dedrift_frame(frame, drift_rate=None)

blscint.ts_statistics module

blscint.ts_statistics.autocorr(ts, remove_spike=False)

Calculate full autocorrelation, normalizing time series to zero mean and unit variance.

blscint.ts_statistics.acf(ts, remove_spike=False)

blscint.ts_statistics.get_stats(ts, pow=1.6666666666666667, use_triangle=True)

Calculate statistics based on normalized time series (to mean 1).

blscint.ts_statistics.triangle(x, L)

blscint.ts_statistics.scint_acf(x, t_d, pow=2)

pow is 2 for square-law; 5/3 for Kolmogorov.

blscint.ts_statistics.noisy_scint_acf(x, t_d, A, W, pow=2, use_triangle=True)

pow is 2 for square-law; 5/3 for Kolmogorov. use_triangle weights the acf model by the triangular function for the acf calculation.

blscint.ts_statistics.noisy_scint_acf_gen(pow=2, use_triangle=True)

pow is 2 for square-law; 5/3 for Kolmogorov.

blscint.ts_statistics.fit_acf(acf, pow=2, use_triangle=True, remove_spike=False)

Routine to fit ideal ACF shapes to empirical autocorrelations. pow is 2 for square-law; 5/3 for Kolmogorov.

blscint.ts_statistics.ts_plots(ts, xlim=None, bins=None)

Plot time series, autocorrelation, and histogram.

blscint.ts_statistics.ts_stat_plots(ts_arr, t_d=None, dt=None)

Plot relevant statistics over a list of time series arrays. Plot in reference to a given scintillation timescale and time resolution.

blscint.ts_statistics.ts_ac_plot(ts_arr, t_d, dt, p=2)

Plot autocorrelations, mean += std dev at each lag over a list of time series arrays. Plot in reference to scintillation timescale and time resolution, up to lag p.

blscint.factors module

Multiplicative factors to take a standard deviation (sigma) value to various measures such as full width at half maximum (FWHM).