PSF and image simulation (wcc_etc.psfsim)

wcc_etc.psfsim.DEFOCUS_REF_WAVELENGTH_M = 5e-07

Reference wavelength of the bundled Zemax defocus runs (from their filenames and Zemax configuration: monochromatic 500 nm Huygens PSFs).

wcc_etc.psfsim.DEFOCUS_REF_FNUM = 15.0

f-number of the Lazuli WCC design the defocus products were computed for (data/config/lazuli.toml).

class wcc_etc.psfsim.DetectorPSFContext(npix, pixel_size_um, plate_scale_mas, wavelength_m, diameter_m, fnum, jitter_sigma_mas=0.0, center=None, oversample=11)[source]

Bases: object

Detector + optics parameters a PSF source needs to render onto the grid.

wavelength_m is the source-weighted effective wavelength of the bandpass times the source SED (wcc_etc.Simulation.effective_wavelength), not the filter pivot wavelength: the diffraction scale is linear in wavelength, so an M5V and an O5V observed through the same broad band do not share a PSF (issue #65).

Parameters:
npix: int
pixel_size_um: float
plate_scale_mas: float
wavelength_m: float
diameter_m: float
fnum: float
jitter_sigma_mas: float = 0.0
center: tuple | None = None
oversample: int = 11
wcc_etc.psfsim.normalize_psf(psf)[source]

Clip negatives and normalize a 2D PSF so it sums to 1.

wcc_etc.psfsim.grid_center(npix)[source]

Nominal (cx, cy) of a PSF rendered with ctx.center=None.

The integer pixel (npix-1)//2 on both axes: the geometric centre of an odd grid, and the pixel just below it on an even grid (where the Airy renderer’s odd stamp lands when cropped down). It is also the origin of an fftconvolve(..., mode="same") kernel, so Simulation.extended_rate_image convolving a host placed here with a PSF leaves the host on the source.

wcc_etc.psfsim.crop_start(n, npix)[source]

First index of an n-long axis kept by center_crop_or_pad to npix.

The source centre (n-1)/2 lands at (n-1)/2 - crop_start(n, npix): grid_center(npix) when n is odd, half a pixel off it when even.

wcc_etc.psfsim.center_crop_or_pad(img, npix, fill=0.0)[source]

Center-crop or zero-pad a 2D array to (npix, npix), preserving the center.

wcc_etc.psfsim.recenter(psf, center, origin=None)[source]

Sub-pixel shift a PSF centred at origin so its centre lands at (cx, cy).

origin defaults to grid_center(npix); pass the actual centre when the array was cropped from an even-sized stamp (see crop_start).

class wcc_etc.psfsim.PSFSource[source]

Bases: object

Base class: render a normalized (sum=1) PSF onto a DetectorPSFContext.

render(ctx)[source]
cache_key()[source]

Hashable key identifying this PSF source’s parameters (not the rendered output). The render also depends on the DetectorPSFContext (wavelength, optics, jitter, npix, oversample), so a render cache must combine this key with the context and/or be invalidated when the simulation state changes — see Simulation._image_render_bundle.

class wcc_etc.psfsim.AiryPSF[source]

Bases: PSFSource

Diffraction-limited Airy PSF rendered on the detector grid (default).

Monochromatic, at ctx.wavelength_m – which the ETC sets to the source-weighted effective wavelength, so the PSF width tracks the colour of the source. For a coadd across the band (which also washes out the monochromatic ring contrast) use PolychromaticPSF or wcc_etc.Simulation.polychromatic_psf().

render(ctx)[source]
class wcc_etc.psfsim.PolychromaticPSF(wavelengths, weights, base=None)[source]

Bases: PSFSource

Photon-weighted coadd of a base PSF rendered at several wavelengths.

Rendering at a single wcc_etc.Simulation.effective_wavelength already removes the colour bias in the PSF width, which is the dominant error (issue #65). This class goes one step further: it renders the base PSF at n_sub sub-band wavelengths and coadds them with the sub-bands’ photon weights, so the diffraction rings are washed out the way a real broadband PSF’s are rather than staying at full monochromatic contrast.

The coadd’s own wavelengths override ctx.wavelength_m; the context still supplies the optics, grid, jitter and centering.

Parameters:
  • wavelengths (Quantity or array_like) – Sub-band wavelengths. A bare array is read as metres.

  • weights (array_like) – Relative photon weight per sub-band; normalized internally to sum to 1.

  • base (PSFSource, optional) – PSF to render at each wavelength. Default AiryPSF.

See also

wcc_etc.Simulation.polychromatic_psf

build one for a simulation.

wcc_etc.spectral.photon_weighted_subbands

chooses the wavelengths and weights.

classmethod from_bandpass(bandpass, spectrum=None, n_sub=7, base=None)[source]

Build the coadd from a bandpass and a source spectrum.

Parameters:
  • bandpass (SpectralElement) – Total throughput.

  • spectrum (SourceSpectrum, optional) – Source spectrum; None weights by throughput alone. The normalization is irrelevant, so an un-normalized spectrum is fine.

  • n_sub (int, optional) – Number of equal-photon-weight sub-bands. Default 7.

  • base (PSFSource, optional) – PSF rendered at each sub-band wavelength. Default AiryPSF.

Return type:

PolychromaticPSF

render(ctx)[source]

Coadd the base PSF rendered at each sub-band wavelength.

cache_key()[source]

Hashable key identifying this PSF source’s parameters (not the rendered output). The render also depends on the DetectorPSFContext (wavelength, optics, jitter, npix, oversample), so a render cache must combine this key with the context and/or be invalidated when the simulation state changes — see Simulation._image_render_bundle.

wcc_etc.psfsim.load_huygens_psf(path, encoding='utf-16')[source]

Load a Zemax Huygens PSF text file into a 2D float array of intensities.

Zemax writes the top row (largest +Y) first. The rows are flipped on load so that, as for a FITS image, the row index increases with +Y and the array can be displayed with origin="lower"; the chief ray from Center point in the header is then at 0-based (x - 1, y - 1).

wcc_etc.psfsim.parse_huygens_header(path, encoding='utf-16')[source]

Read the # comment header of a Zemax Huygens PSF export.

Zemax puts the things that make the array interpretable – the detector data spacing, the defocus, the wavelength range, the Strehl ratio – in comment lines above the numbers. What it does not put there is the reference wavelength of the run or the f-number of the design; those come from the Zemax configuration and are supplied to huygens_txt_to_fits().

Parameters:
  • path (str) – Path to the Zemax .txt export.

  • encoding (str, optional) – Text encoding. Zemax writes UTF-16 by default.

Returns:

Any of src_um_per_pix, defocus_waves, strehl, wl_min_um, wl_max_um, grid_size, data_area_um that could be parsed, plus chief_ray_pix: the 0-based (x, y) of Zemax’s Center point in the array load_huygens_psf() returns.

Return type:

dict

wcc_etc.psfsim.huygens_txt_to_fits(txt_path, fits_path, ref_wavelength_m=5e-07, ref_fnum=15.0, overwrite=False, encoding='utf-16')[source]

Convert a Zemax Huygens PSF export into a self-describing FITS product.

A bare intensity array is not enough to render a PSF onto a detector: you also need to know the sampling it was computed at, and the wavelength and f-number it was computed for, or a renderer has no choice but to ignore the optics it is asked about. This writes all of that into the header (PIXSCALE, WAVELEN, FNUM, DEFOCUSW) so DefocusPSF can scale rather than guess.

Parameters:
  • txt_path (str) – Zemax Huygens PSF .txt export.

  • fits_path (str) – Output FITS path.

  • ref_wavelength_m (float, optional) – Wavelength the Zemax run was computed at, in metres. Not recoverable from the export, so it must be supplied; defaults to the 500 nm the bundled WCC products were run at.

  • ref_fnum (float, optional) – f-number of the design the run was computed for. Defaults to the Lazuli WCC f/15.

  • overwrite (bool, optional) – Overwrite an existing output file. Default False.

  • encoding (str, optional) – Encoding of the .txt export. Default ‘utf-16’.

Returns:

The path written.

Return type:

str

wcc_etc.psfsim.load_psf_fits(path)[source]

Load a FITS PSF product into (data, meta).

Parameters:

path (str) – Path to a FITS file whose primary HDU holds a 2D PSF.

Returns:

  • data (ndarray) – The 2D intensity array.

  • meta (dict) – src_um_per_pix, ref_wavelength_m and ref_fnum where the header supplies PIXSCALE, WAVELEN (nm) and FNUM, plus defocus_waves from DEFOCUSW, and chief_ray_pix (0-based (x, y)) from CRPIX1/CRPIX2. Missing keywords are simply absent.

wcc_etc.psfsim.WAVELENGTH_SCALING_MODES = ('despace', 'waves', 'none')

Scaling laws _ResampledPSF can use to map its sampled array onto the optics of a DetectorPSFContext. See the class docstring.

class wcc_etc.psfsim.DefocusPSF(path, src_um_per_pix=None, ref_wavelength_m=None, ref_fnum=None, wavelength_scaling='despace', encoding='utf-16')[source]

Bases: _ResampledPSF

A defocused PSF loaded from a bundled FITS product or a Zemax .txt export.

A FITS product supplies its own sampling and reference optics from its header (PIXSCALE, WAVELEN, FNUM), so the render can honour the optics it is handed; see _ResampledPSF for the scaling laws. A raw .txt export only carries the sampling, so scaling it needs ref_wavelength_m and ref_fnum passed in.

Parameters:
  • path (str) – Path to a .fits PSF product or a Zemax Huygens .txt export.

  • src_um_per_pix (float, optional) – Override the sampling read from the file.

  • ref_wavelength_m (float, optional) – Override the reference optics read from the file.

  • ref_fnum (float, optional) – Override the reference optics read from the file.

  • wavelength_scaling (str, optional) – See _ResampledPSF. Default "despace".

  • encoding (str, optional) – Encoding of a .txt export. Default ‘utf-16’.

cache_key()[source]

Hashable key identifying this PSF source’s parameters (not the rendered output). The render also depends on the DetectorPSFContext (wavelength, optics, jitter, npix, oversample), so a render cache must combine this key with the context and/or be invalidated when the simulation state changes — see Simulation._image_render_bundle.

class wcc_etc.psfsim.CustomPSF(source, src_um_per_pix, ref_wavelength_m=None, ref_fnum=None, wavelength_scaling='despace', encoding='utf-16', chief_ray_pix=None)[source]

Bases: _ResampledPSF

A custom PSF from an ndarray, a FITS product, or a Huygens .txt export.

Declare ref_wavelength_m and ref_fnum (or read them from a FITS header) if the array should track the optics it is rendered onto; otherwise pass wavelength_scaling="none" to render it at its native sampling. See _ResampledPSF.

wcc_etc.psfsim.saturation_mask_from_image_e(sensor, image_e)[source]

Boolean mask of pixels at/over the sensor’s per-frame saturation ceiling.

Parameters:
  • sensor (Sensor) – Provides saturation_e = min(full well, (adc_max - bias) * gain).

  • image_e (ndarray) – Per-pixel charge in electrons (per frame). Pass a clean image for a planning answer; a noisy realization (ImageSimulator.simulate with add_noise=True) can flag pixels the clean prediction does not.

Returns:

True where image_e >= sensor.saturation_e.

Return type:

ndarray of bool

class wcc_etc.psfsim.SimulatedImage(image_e, image_clean, saturation_mask, gain, bias_level, npix, pixel_scale_mas, psf)[source]

Bases: object

Result of an ImageSimulator.simulate() call.

Parameters:
image_e: ndarray
image_clean: ndarray
saturation_mask: ndarray
gain: float
bias_level: float
npix: int
pixel_scale_mas: float
psf: PSFSource
to_adu()[source]

Electrons -> ADU via gain, plus the bias level.

to_fitsimg()[source]

Wrap the electron image in a FitsImg for photometry/plotting.

plot_image(backend='mpl', **kwargs)[source]

Plot this image (single panel). backend=’mpl’ or ‘bokeh’.

plot_image_row(backend='mpl', **kwargs)[source]

Plot the 3-panel row (PSF+noise, PSF, saturation mask).

plot_radial(backend='mpl', **kwargs)[source]

Plot the azimuthally-averaged radial profile.

plot_encircled_energy(backend='mpl', **kwargs)[source]

Plot the encircled-energy curve (marks the 90% EE radius by default; pass ee_target=None to disable or another fraction to override).

plot(backend='mpl', **kwargs)

Plot this image (single panel). backend=’mpl’ or ‘bokeh’.

class wcc_etc.psfsim.ImageSimulator(simulation, npix=300, oversample=11)[source]

Bases: object

Render a point source onto a detector grid with noise, driven by an ETC Simulation.

classmethod from_sensor_and_scene(sensor, scene, npix=300, oversample=11)[source]

Build an ImageSimulator from a sensor name (e.g. ‘sony:r’) and a Scene.

classmethod from_sensorfilter(sensorfilter, scene, npix=300, oversample=11)[source]

Build an ImageSimulator from a sensorfilter label (e.g. ‘zwo:r+1’).

The PSF is chosen from the sensor’s focus_level in sensor_info and exposed as default_psf. simulate(), get_image_snr(), and get_image_exptime_for_snr() all use it when no psf= argument is given.

property plate_scale_mas

Detector plate scale in mas/pixel, from the sensor + telescope.

property default_psf

The PSF used when simulate() is called without psf=.

render_psf(psf=None, jitter_sigma_mas=None, npix=None, center=None)[source]

Render a PSF onto the detector grid, normalized to unit sum.

This is the bare PSF with no source flux, background, or noise — useful for comparing PSF shapes, peak-pixel fractions, radial profiles, and encircled energy across focus levels or jitter values. Pair it with plate_scale_mas to plot in mas.

Parameters:
  • psf (PSFSource, optional) – PSF model to render. Defaults to default_psf.

  • jitter_sigma_mas (float, optional) – Override the telescope jitter (mas). Defaults to the telescope value; pass 0 to see the un-blurred diffraction limit.

  • npix (int, optional) – Override the (square) grid size for this render only.

  • center (tuple, optional) – Sub-pixel (cx, cy) center. Defaults to grid_center(npix), the integer pixel (npix-1)//2.

Returns:

The (npix, npix) normalized PSF image (sum = 1).

Return type:

numpy.ndarray

plot_psf(psf=None, **kwargs)[source]

Plot the bare normalized PSF; see plotting.plot_psf.

simulate(time=None, psf=None, jitter_sigma_mas=None, center=None, add_noise=True, seed=None)[source]

Simulate a detector image for the given exposure time and PSF.

Parameters:
  • time (float or Quantity, optional) – Exposure time (seconds if a bare float). Defaults to the Simulation’s meta[‘time’].

  • psf (PSFSource, optional) – PSF model to render. Defaults to default_psf — diffraction-limited unless the simulation came from from_sensorfilter, in which case it is the PSF for that filter’s focus level.

  • jitter_sigma_mas (float, optional) – Override the telescope jitter (mas). Defaults to the telescope value.

  • center (tuple, optional) – Sub-pixel (cx, cy) center for the PSF (and any Sersic host attached to it). Defaults to grid_center(npix), the integer pixel (npix-1)//2.

  • add_noise (bool, optional) – If True, apply Poisson shot + dark noise and Gaussian read noise. If False, return the noiseless electron image. Default True.

  • seed (int, optional) – Seed for the random generator, for reproducible noise.

Returns:

Holds the electron image, the noiseless image, and a saturation mask (computed from the returned image_e against adc_max / well).

Return type:

SimulatedImage

wcc_etc.psfsim.howell_center(postage_stamp)[source]

Howell centroiding, from Howell’s Handbook of CCD astronomy

INPUT:

postage_stamp - A 2d numpy array to do the centroiding

OUTPUT:

x and y center of the numpy array

NOTES: Many thanks to Thomas Beatty and the MINERVAphot.py pipeline for this method see here: https://github.com/TGBeatty/MINERVAphot/blob/master/MINERVAphot.py

wcc_etc.psfsim.psf_center(image)[source]

Measured center (xc, yc) of a PSF image for radial reducers.

Uses Howell centroiding so radii are measured from the actual light distribution rather than a fixed grid index. This matters because render_detector_psf forces an odd grid (peak on the exact center pixel) and center_crop_or_pad then brings an even target grid back down, landing the centroid on an integer pixel — a full pixel away from n//2 and half a pixel from the geometric center (n-1)/2. It also tracks the centroid of asymmetric defocused (DefocusPSF) bands, which are not guaranteed symmetric about the grid center.

Falls back to the geometric center ((nx-1)/2, (ny-1)/2) for an empty or degenerate image where the centroid is undefined (zero/negative total flux or a non-finite centroid, e.g. a flat field).

wcc_etc.psfsim.apply_jitter(data, jitter_mas, pixel_scale)[source]

Apply jitter to the input data.

INPUT:

data: 2d data jitter_mas: jitter in mas pixel_scale: mas/pixel

wcc_etc.psfsim.tiff_to_fits(tiff_file, output_dir)[source]

Convert a multi-frame TIFF file to a series of FITS files.

wcc_etc.psfsim.create_master_flat_from_fits(directory, output_filename)[source]
wcc_etc.psfsim.apply_nl_scaling(df_nl, data, how='makenonlinear', scale=1)[source]

Apply non-linear scaling to the input data using the provided scaling factors.

Parameters:

df_nl (pd.DataFrame): DataFrame containing ‘mean_value’ and ‘scaling_factor’ columns. data (np.ndarray): Input data array to be scaled. how (str): ‘makenonlinear’ to apply non-linearity, ‘correctnonlinear’ to reverse it.

RETURNS:

data_scaled (np.ndarray): Scaled data array.

EXAMPLE:

NOTES:
  • only applies scaling to positive values; zero or negative values are unchanged.

  • Should only be applied after bias subtraction.

class wcc_etc.psfsim.FitsImgList(data_list, center_list=None, **kwargs)[source]

Bases: object

aperture_photometry(**kwargs)[source]

Perform aperture photometry on all images in the list.

Parameters:

**kwargs (dict) – Keyword arguments to pass to the aperture_photometry method of FitsImg.

Returns:

results_list – List of dictionaries containing photometry results for each image.

Return type:

list of dict

plot_photometry(axes=None)[source]

Plot the photometry results stored in self.df_phot.

class wcc_etc.psfsim.FitsImg(filename=None, data=None, header=None, center=None, dark_current_rate=0.0, exp_time=1.0, read_noise_rms=0.0, imgnumber=0)[source]

Bases: object

crop(x, y, w, h)[source]

Crop to a box centered at (x,y), of size w x h

cropcenter(w, h, points=False)[source]

Returns an image array around the center of an image array.

cropcentroid(w, h)[source]

Crop to a box centered at the centroid, of size w x h

get_centroid(plot_cross=False, ax=None, plot_lines=False)[source]

Find centroid using Howell centroiding.

See phothelp for the method

get_centroid_line_cut(line='X', plot=False, ax=None)[source]

Horizontal (line='X') or vertical (line='Y') cut through the centroid.

INPUT:

line - ‘X’ for the row through the centroid, ‘Y’ for the column plot - draw the cut on ax (or a new axis)

OUTPUT:

cut - 1D array along the requested direction

plot(stretch='hist', cmap='gray', origin='lower', ax=None, colorbar=False, title='', vmin=None, vmax=None, dpi=200)[source]
get_radial_profile(rmax=None, plot=False, z=2.0, return_hwzm=False, ax=None, xcen=None, ycen=None, annulus_width=1, subtract_min=False)[source]

Plot radial profile

aperture_photometry(r_ap=3.0, r_in=6.0, r_out=8.0, gain=1.0, center=None, bkg_sigma_clip=3.0, bkg_maxiters=5, plot=True, ax=None, verbose=True, vmin=None, vmax=None, cmap='viridis', origin='lower', stretch='hist', colorbar=True)[source]

Perform circular aperture photometry on the current self.data image using photutils.

Parameters:
  • r_ap (float) – Aperture radius in pixels.

  • r_in (float) – Inner and outer radii for background annulus in pixels. If None or invalid, the image median is used as a background estimate.

  • r_out (float) – Inner and outer radii for background annulus in pixels. If None or invalid, the image median is used as a background estimate.

  • gain (float) – e/ADU (set to 1.0 if self.data is already in electrons).

  • center (None or (y,x)) – If provided, use this center; otherwise the method will try self.center, then the peak pixel, then image center.

  • bkg_sigma_clip (float) – Sigma for sigma-clipped background estimation in the annulus.

  • bkg_maxiters (int) – Max iterations for sigma clipping.

Returns:

Dictionary containing ap_sum, bkg_mean_per_pix, bkg_std_per_pix, bkg_sum, net_flux, flux_err, snr, ap_area, ann_area, ann_pixels_used, center

Return type:

dict

wcc_etc.psfsim.solve_time_for_snr(snr, A, B, C)[source]

Solve SNR = A*t / sqrt(B*t + C) for the positive root t.

A, B, C may be scalars or broadcastable arrays (A = signal rate, B = variance rate, C = constant read-noise variance). Returns t in the same shape (a float if all inputs are scalar). Entries with A <= 0 return +inf (no solution). snr must be finite and > 0: a negative target would be squared into a positive time and NaN would propagate silently (#86).

wcc_etc.psfsim.RADIUS_TIE_TOL_PIX = 1e-06

Pixels whose centre radii agree within this (pixels) form one ring.

wcc_etc.psfsim.diffuse_enclosed(diffuse_per_pix, order, n_pix)[source]

Enclosed diffuse charge per radius: a scalar is uniform (per_pix * n_pix); a 2D image (e.g. a Sersic host) is summed in the same radial order.

wcc_etc.psfsim.aperture_snr_radial(psf_norm, plate_scale_mas, source_e_total, diffuse_per_pix, dark_per_pix, read_noise, contaminant_e_total=0.0)[source]

SNR as a function of circular-aperture radius for a rendered PSF.

The PSF (sum=1) sets how much source light falls inside each radius; the per-pixel diffuse (sky+host), dark, and read-noise terms set the background noise that grows with the number of aperture pixels.

Parameters:
  • psf_norm (ndarray) – Normalized (sum=1) PSF on the detector grid.

  • plate_scale_mas (float) – Detector plate scale, mas/pixel (to report radii in mas).

  • source_e_total (float) – Total source electrons (all of the PSF, before aperture clipping).

  • diffuse_per_pix (float or ndarray) – Per-pixel diffuse (sky / resolved host) electrons; a 2D image on the PSF grid for a spatially varying host (see Simulation.extended_rate_image).

  • dark_per_pix (float) – Per-pixel dark-current electrons.

  • read_noise (float) – Read noise (electrons rms per pixel).

  • contaminant_e_total (float, optional) – Total electrons from unresolved contaminants co-located with the source (e.g. an unresolved host). Distributed by the same PSF as the source, so it adds contaminant_e_total * enclosed_fraction of shot-noise variance inside each aperture, and contributes no signal. Default 0.

Returns:

‘r_mas’, ‘enclosed_fraction’, ‘n_pix’, ‘signal_e’, ‘noise_e’, ‘snr’.

Return type:

dict of ndarrays, sorted by ascending radius

wcc_etc.psfsim.select_aperture(profile, r_aper_mas=None, ee_frac=None, optimize=False)[source]

Index into a radial aperture profile for the chosen aperture mode.

Precedence: optimize > explicit r_aper_mas > ee_frac. optimize maximizes profile['snr'] (forward) or, if present, minimizes profile['time_s'] (inverse). Every profile entry is a complete ring (see _radial_cumulative), so the selected radius reproduces its aperture.

r_aper_mas must be finite, > 0 and inside the render grid (a radius below the innermost ring selects that ring); ee_frac must be finite in (0, 1]. Out-of-range requests raise rather than being clipped silently.

wcc_etc.psfsim.aperture_time_for_snr(psf_norm, plate_scale_mas, source_rate_total, diffuse_rate_per_pix, dark_rate_per_pix, read_noise, n_reads=1, snr=None, r_aper_mas=None, ee_frac=None, optimize=False, contaminant_rate_total=0.0)[source]

Exposure time (s) to reach snr for a rendered PSF, per aperture mode.

Rates are per second (the time dependence is solved for analytically). The read-noise variance is incurred n_reads times. Aperture precedence matches select_aperture: optimize (fastest radius) > r_aper_mas > ee_frac.

Returns {‘time_s’, ‘snr’, ‘r_aper_mas’, ‘enclosed_fraction’, ‘n_pix’}. Raises ValueError if the target is unreachable (no source flux inside the chosen aperture) instead of returning an infinite time.

wcc_etc.psfsim.calc_hwzm(x, y, z=20)[source]

Calculates the HWHM at the Z-th maximum for a given dataset, by finding the roots of splines.

INPUTS:

x - x input array y - y input array

OUTPUT:

HWZM The Half Width at Z-th Max of the data

EXAMPLE:

wcc_etc.psfsim.calc_hwhm(x, y)[source]

Calculates the HWHM for a given dataset, by finding the roots of splines.

INPUTS:

x - x input array y - y input array

OUTPUT:

HWZM The Half Width at Z-th Max of the data

EXAMPLE: