# Source code for stingray.pulse.modeling

```import numpy as np
from astropy.modeling import models, fitting

__all__ = ["sinc_square_model", "sinc_square_deriv", "fit_sinc",
"fit_gaussian", "SincSquareModel"]

def sinc(x):
"""
Calculate a sinc function.

sinc(x)=sin(x)/x

Parameters
----------
x : array-like

Returns
-------
values : array-like
"""
values = np.sinc(x/np.pi)
return values

[docs]def sinc_square_model(x, amplitude=1., mean=0., width=1.):
"""
Calculate a sinc-squared function.

(sin(x)/x)**2

Parameters
----------
x: array-like

Other Parameters
----------
amplitude : float
the value for x=mean
mean : float
mean of the sinc function
width : float
width of the sinc function

Returns
-------
sqvalues : array-like
Return square of sinc function

Examples
--------
>>> sinc_square_model(0, amplitude=2.)
2.0
"""
sqvalues = amplitude * sinc((x-mean)/width) ** 2
return sqvalues

[docs]def sinc_square_deriv(x, amplitude=1., mean=0., width=1.):
"""
Calculate partial derivatives of sinc-squared.

Parameters
----------
x: array-like

Other Parameters
----------
amplitude : float
the value for x=mean
mean : float
mean of the sinc function
width : float
width of the sinc function

Returns
-------
d_amplitude : array-like
partial derivative of sinc-squared function
with respect to the amplitude
d_mean : array-like
partial derivative of sinc-squared function
with respect to the mean
d_width : array-like
partial derivative of sinc-squared function
with respect to the width

Examples
--------
>>> np.allclose(sinc_square_deriv(0, amplitude=2.), [1., 0., 0.])
True
"""
x_is_zero = x == mean

d_x = 2 * amplitude * \
sinc((x-mean)/width) * (
x * np.cos((x-mean)/width) -
np.sin((x - mean) / width)) / ((x - mean) / width) ** 2
d_x = np.asarray(d_x)
d_amplitude = sinc((x-mean)/width)**2
d_x[x_is_zero] = 0

d_mean = d_x*(-1/width)
d_width = d_x*(-(x-mean)/(width)**2)

return [d_amplitude, d_mean, d_width]

_SincSquareModel = models.custom_model(sinc_square_model,
fit_deriv=sinc_square_deriv)

[docs]class SincSquareModel(_SincSquareModel):
def __reduce__(cls):
members = dict(cls.__dict__)
return (type(cls), (), members)

[docs]def fit_sinc(x, y, amp=1.5, mean=0., width=1., tied={}, fixed={}, bounds={},
obs_length=None):
"""
Fit a sinc function to x,y values.

Parameters
----------
x : array-like
y : array-like

Other Parameters
----------------
amp : float
The initial value for the amplitude

mean : float
The initial value for the mean of the sinc

obs_length : float
The length of the observation. Default None. If it's defined, it
fixes width to 1/(pi*obs_length), as expected from epoch folding
periodograms

width : float
The initial value for the width of the sinc. Only valid if
obs_length is 0

tied : dict

fixed : dict

bounds : dict
Parameters to be passed to the [astropy models]_

Returns
-------
sincfit : function
The best-fit function, accepting x as input
and returning the best-fit model as output

References
----------
.. [astropy models] http://docs.astropy.org/en/stable/api/astropy.modeling.functional_models.Gaussian1D.html
"""
if obs_length is not None:
width = 1 / (np.pi * obs_length)
fixed["width"] = True

sinc_in = SincSquareModel(amplitude=amp, mean=mean, width=width, tied=tied,
fixed=fixed, bounds=bounds)
fit_s = fitting.LevMarLSQFitter()
sincfit = fit_s(sinc_in, x, y)
return sincfit

[docs]def fit_gaussian(x, y, amplitude=1.5, mean=0., stddev=2., tied={}, fixed={},
bounds={}):
"""
Fit a gaussian function to x,y values.

Parameters
----------
x : array-like
y : array-like

Other Parameters
----------------
amplitude : float
The initial value for the amplitude
mean : float
The initial value for the mean of the gaussian function
stddev : float
The initial value for the standard deviation of the gaussian function
tied : dict
fixed : dict
bounds : dict
Parameters to be passed to the [astropy models]_

Returns
-------
g : function
The best-fit function, accepting x as input
and returning the best-fit model as output
"""
g_in = models.Gaussian1D(amplitude=amplitude, mean=mean, stddev=stddev,
tied=tied, fixed=fixed, bounds=bounds)
fit_g = fitting.LevMarLSQFitter()
g = fit_g(g_in, x, y)
return g
```