# -*- coding: utf-8 -*-
import logging
import numpy as np
from scipy.interpolate import interp1d
from scipy.ndimage.filters import convolve
logger = logging.getLogger(__name__)
[docs]def apply_broadening(ipres, x_seg, y_seg, type="gauss", sme=None):
"""
Broaden a spectrum by instrument resolution, with a given broadening type
Parameters
----------
ipres : float
instrument resolution
x_seg : array (npoints,)
x values (wavelength) of the spectrum to broaden
y_seg : array (npoints,)
y values (intensities) of the spectrum to broaden
type : {str, None}, optional
broadening type to apply. Options are "gauss", "sinc", "table", None.
If None, will try to use sme.iptype to determine type. "table" requires
keyword sme to be passed as well. See functions the respective for details.
(default: "gauss")
sme : SME_Struct, optional
sme structure with instrument profile data, required only for
type="table" or type=None (default: None)
Raises
------
AttributeError
if type requires SME_Struct, but its missing
passed type not recognized
Returns
-------
y_seg : array (npoints,)
broadened intensity spectrum
"""
if sme is None and type in ["table", None]:
raise AttributeError(f"SME structure needs to be passed when using type={type}")
# Using the log-linear wavelength grid requires using the first point
# for specifying the width of the instrumental profile
hwhm = 0.5 * x_seg[0] / ipres if ipres > 0 else 0
if type is None:
type = sme.iptype
type = type.casefold()
if type == "table":
y_seg = tablebroad(x_seg, y_seg, sme.ip_x, sme.ip_y)
elif type == "gauss":
y_seg = gaussbroad(x_seg, y_seg, hwhm)
elif type == "sinc":
y_seg = sincbroad(x_seg, y_seg, hwhm)
else:
raise AttributeError(f"Unknown instrument profile type - {type}")
return y_seg
[docs]def tablebroad(w, s, xip, yip):
"""
Convolves a spectrum with an arbitrary instrumental profile.
Parameters
-------
w : array of size (n,)
wavelength scale of spectrum to be smoothed
s : array of size (n,)
spectrum to be smoothed
xip : array of size (m,)
x points of the instrument profile
yip : array of size (m,)
y points of the instrument profile
Returns
-------
sout: array[n]
the smoothed spectrum.
"""
"""
History
-------
22-May-92 JAV
Switched instrumental profile from multiple gaussians
to gaussian with power-law wings.
04-Aug-92 JAV
Renamed from ipsmo.pro# changed f/ procedure to function.
Switched f/ 10 to 15 Hamilton pixels in each wing.
20-Oct-92 JAV
Switched from gpfunc to ipfun (4 to 5 par).
23-Aug-94 JAV
Switched to explicitly passed IPs.
Oct-18 AW
Python version
"""
# Define sizes
dsdh = np.abs(np.min(np.diff(xip)))
nip = 2 * int(15 / dsdh) + 1 ## profile points
# Generate instrumental profile on model pixel scale.
x = (
np.arange(nip, dtype=float) - (nip - 1) / 2
) * dsdh # offset in Hamilton pixels
ip = interp1d(xip, yip, kind="cubic")(x)
# ip = bezier_interp(xip, yip, x) # spline onto new scale
ip = ip[::-1] # reverse for convolution
ip = ip / np.sum(ip) # ensure unit area
# Pad spectrum ends to minimize impact of Fourier ringing.
sout = convolve(s, ip, mode="nearest")
return sout # return convolved spectrum
[docs]def gaussbroad(w, s, hwhm):
"""
Smooths a spectrum by convolution with a gaussian of specified hwhm.
Parameters
-------
w : array[n]
wavelength scale of spectrum to be smoothed
s : array[n]
spectrum to be smoothed
hwhm : float
half width at half maximum of smoothing gaussian.
Returns
-------
sout: array[n]
the gaussian-smoothed spectrum.
"""
"""
History
--------
Dec-90 GB,GM
Rewrote with fourier convolution algorithm.
Jul-91 AL
Translated from ANA to IDL.
22-Sep-91 JAV
Relaxed constant dispersion check# vectorized, 50% faster.
05-Jul-92 JAV
Converted to function, handle nonpositive hwhm.
Oct-18 AW
Python version
"""
# Warn user if hwhm is negative.
if hwhm < 0:
logger.warning("Forcing negative smoothing width to zero.")
# Return input argument if half-width is nonpositive.
if hwhm <= 0:
return s # true: no broadening
# Calculate (uniform) dispersion.
nw = len(w) ## points in spectrum
wrange = w[-1] - w[0]
dw = wrange / (nw - 1) # wavelength change per pixel
# Make smoothing gaussian; extend to 4 sigma.
# 4.0 / sqrt(2.0 * alog(2.0)) = 3.3972872
# sqrt(alog(2.0)) = 0.83255461
# sqrt(alog(2.0) / pi) = 0.46971864
# (*1.0000632 to correct for >4 sigma wings)
if hwhm >= 5 * wrange:
return np.full(nw, np.sum(s) / nw)
## points in half gaussian
nhalf = int(3.3972872 * hwhm / dw)
## points in gaussian (odd!)
ng = 2 * nhalf + 1
# wavelength scale of gaussian
wg = dw * (np.arange(ng, dtype=float) - (ng - 1) / 2)
# convenient absisca
xg = (0.83255461 / hwhm) * wg
# unit area gaussian w / FWHM
gpro = (0.46974832 * dw / hwhm) * np.exp(-xg * xg)
gpro = gpro / np.sum(gpro)
# Pad spectrum ends to minimize impact of Fourier ringing.
sout = convolve(s, gpro, mode="nearest")
return sout
[docs]def sincbroad(w, s, hwhm):
"""
Smooths a spectrum by convolution with a sinc function of specified hwhm.
Parameters
------
w : array of size (n,)
wavelength scale of spectrum to be smoothed
s : array of size (n,)
spectrum to be smoothed
hwhm : float
half width at half maximum of smoothing gaussian.
Returns
-------
sout : array of size (n,)
the sinc-smoothed spectrum.
"""
"""
History
-------
Dec-90 GB,GM
Rewrote with fourier convolution algorithm.
Jul-91 AL
Translated from ANA to IDL.
22-Sep-91 JAV
Relaxed constant dispersion check# vectorized, 50% faster.
05-Jul-92 JAV
Converted to function, handle nonpositive hwhm.
14-Nov-93 JAV
Adapted from macbro.pro
23-Apr-93 JAV
Verified that convolution kernel has specified hwhm. For IR FTS
spectra: hwhm=0.0759 Angstroms, max change in profile is 0.4% of continuum.
Oct-18 AW
Python Version
"""
# Warn user if hwhm is negative.
if hwhm < 0:
logger.warning("Forcing negative smoothing width to zero.")
# Return input argument if half-width is nonpositive.
if hwhm <= 0:
return s # true: no broadening
# Calculate (uniform) dispersion.
nw = len(w) ## points in spectrum
dw = (w[-1] - w[0]) / (nw - 1) # wavelength change per pixel
# Make sinc function out to 20th zero-crossing on either side. Error due to
# ignoring additional lobes is less than 0.2% of continuum. Reducing extent
# to 10th zero-crossing doubles maximum error.
fwhm = 2.0 * hwhm # full width at half maximum
rperfw = 0.26525 # radians per fwhm of sinc
xrange = 20 * np.pi # 20th zero of sinc (radians)
wrange = xrange * fwhm * rperfw # 20th zero of sinc (wavelength)
nhalf = int(wrange / dw + 0.999) ## points in half sinc
nsinc = 2 * nhalf + 1 ## points in sinc (odd!)
wsinc = (np.arange(nsinc, dtype=float) - nhalf) * dw # absissca (wavelength)
xsinc = wsinc / (fwhm * rperfw) # absissca (radians)
xsinc[nhalf] = 1.0 # avoid divide by zero
sinc = np.sin(xsinc) / xsinc # calculate sinc
sinc[nhalf] = 1.0 # insert midpoint
xsinc[nhalf] = 0.0 # fix xsinc
sinc = sinc / np.sum(sinc) # normalize sinc
# Pad spectrum ends to minimize impact of Fourier ringing.
sout = convolve(s, sinc, mode="nearest")
return sout