B q\@sdZddlmZddlZddlmZmZmZddl m Z m Z ddl m Z ddlmZdd lmZmZd d d d ddddddddddddddddddd d!d"d#gZd$ejZeeejjZd%ed%ed%ZGd&ddeZGd'ddeZ Gd(ddeZ!Gd)ddeZ"Gd*ddeZ#Gd+ddeZ$Gd,ddeZ%Gd-ddeZ&Gd.ddeZ'Gd/d0d0eZ(Gd1ddeZ)Gd2d#d#eZ*Gd3ddeZ+Gd4ddeZ,Gd5ddeZ-Gd6ddeZ.Gd7d"d"eZ/Gd8d9d9eZ0Gd:d;d;eZ1Gdd d eZ4Gd?d!d!eZ5Gd@ddeZ6GdAddeZ7GdBd d eZ8GdCd d eZ9GdDd d eZ:GdEddeZ;dS)FzMathematical models.) OrderedDictN)Fittable1DModelFittable2DModelModelDefinitionError) ParameterInputParameterError)ellipse_extent)units)Quantity UnitsError AiryDisk2DMoffat1DMoffat2DBox1DBox2DConst1DConst2D Ellipse2DDisk2D Gaussian1D Gaussian2DLinear1D Lorentz1D MexicanHat1D MexicanHat2DRedshiftScaleFactorScaleMultiplySersic1DSersic2DShiftSine1D Trapezoid1DTrapezoidDisk2DRing2DVoigt1Dg@c@sveZdZdZeddZeddZededfdZddd Z e d d Z e d d Z e ddZe ddZddZdS)ra One dimensional Gaussian model. Parameters ---------- amplitude : float Amplitude of the Gaussian. mean : float Mean of the Gaussian. stddev : float Standard deviation of the Gaussian. Notes ----- Model formula: .. math:: f(x) = A e^{- \frac{\left(x - x_{0}\right)^{2}}{2 \sigma^{2}}} Examples -------- >>> from astropy.modeling import models >>> def tie_center(model): ... mean = 50 * model.stddev ... return mean >>> tied_parameters = {'mean': tie_center} Specify that 'mean' is a tied parameter in one of two ways: >>> g1 = models.Gaussian1D(amplitude=10, mean=5, stddev=.3, ... tied=tied_parameters) or >>> g1 = models.Gaussian1D(amplitude=10, mean=5, stddev=.3) >>> g1.mean.tied False >>> g1.mean.tied = tie_center >>> g1.mean.tied Fixed parameters: >>> g1 = models.Gaussian1D(amplitude=10, mean=5, stddev=.3, ... fixed={'stddev': True}) >>> g1.stddev.fixed True or >>> g1 = models.Gaussian1D(amplitude=10, mean=5, stddev=.3) >>> g1.stddev.fixed False >>> g1.stddev.fixed = True >>> g1.stddev.fixed True .. plot:: :include-source: import numpy as np import matplotlib.pyplot as plt from astropy.modeling.models import Gaussian1D plt.figure() s1 = Gaussian1D() r = np.arange(-5, 5, .01) for factor in range(1, 4): s1.amplitude = factor plt.plot(r, s1(r), color=str(0.25 * factor), lw=2) plt.axis([-5, 5, -1, 4]) plt.show() See Also -------- Gaussian2D, Box1D, Moffat1D, Lorentz1D r)defaultrN)r(bounds@cCs |j}||j}||||fS)a Tuple defining the default ``bounding_box`` limits, ``(x_low, x_high)`` Parameters ---------- factor : float The multiple of `stddev` used to define the limits. The default is 5.5, corresponding to a relative error < 1e-7. Examples -------- >>> from astropy.modeling.models import Gaussian1D >>> model = Gaussian1D(mean=0, stddev=2) >>> model.bounding_box (-11.0, 11.0) This range can be set directly (see: `Model.bounding_box `) or by using a different factor, like: >>> model.bounding_box = model.bounding_box(factor=2) >>> model.bounding_box (-4.0, 4.0) )meanstddev)selffactorx0dxr1Alib/python3.7/site-packages/astropy/modeling/functional_models.py bounding_boxzs zGaussian1D.bounding_boxcCs |jtS)z$Gaussian full width at half maximum.)r,GAUSSIAN_SIGMA_TO_FWHM)r-r1r1r2fwhmszGaussian1D.fwhmcCs"|td||d|dS)z, Gaussian1D model function. gr')npexp)x amplituder+r,r1r1r2evaluateszGaussian1D.evaluatecCs\td|d||d}|||||d}||||d|d}|||gS)z8 Gaussian1D model function derivatives. gr')r6r7)r8r9r+r, d_amplitudeZd_meanZd_stddevr1r1r2 fit_derivszGaussian1D.fit_derivcCs |jjdkrdSd|jjiSdS)Nr8)r+unit)r-r1r1r2 input_unitss zGaussian1D.input_unitscCs&td|dfd|dfd|dfgS)Nr+r8r,r9y)r)r- inputs_unit outputs_unitr1r1r2_parameter_units_for_data_unitss  z*Gaussian1D._parameter_units_for_data_units)r*)__name__ __module__ __qualname____doc__rr9r+ FLOAT_EPSILONr,r3propertyr5 staticmethodr:r=r?rCr1r1r1r2r!sP     cseZdZdZeddZeddZeddZeddZeddZ eddZ ej ej ej ddddffdd Z e d d Ze d d ZdddZeddZeddZe ddZddZZS)ral Two dimensional Gaussian model. Parameters ---------- amplitude : float Amplitude of the Gaussian. x_mean : float Mean of the Gaussian in x. y_mean : float Mean of the Gaussian in y. x_stddev : float or None Standard deviation of the Gaussian in x before rotating by theta. Must be None if a covariance matrix (``cov_matrix``) is provided. If no ``cov_matrix`` is given, ``None`` means the default value (1). y_stddev : float or None Standard deviation of the Gaussian in y before rotating by theta. Must be None if a covariance matrix (``cov_matrix``) is provided. If no ``cov_matrix`` is given, ``None`` means the default value (1). theta : float, optional Rotation angle in radians. The rotation angle increases counterclockwise. Must be None if a covariance matrix (``cov_matrix``) is provided. If no ``cov_matrix`` is given, ``None`` means the default value (0). cov_matrix : ndarray, optional A 2x2 covariance matrix. If specified, overrides the ``x_stddev``, ``y_stddev``, and ``theta`` defaults. Notes ----- Model formula: .. math:: f(x, y) = A e^{-a\left(x - x_{0}\right)^{2} -b\left(x - x_{0}\right) \left(y - y_{0}\right) -c\left(y - y_{0}\right)^{2}} Using the following definitions: .. math:: a = \left(\frac{\cos^{2}{\left (\theta \right )}}{2 \sigma_{x}^{2}} + \frac{\sin^{2}{\left (\theta \right )}}{2 \sigma_{y}^{2}}\right) b = \left(\frac{\sin{\left (2 \theta \right )}}{2 \sigma_{x}^{2}} - \frac{\sin{\left (2 \theta \right )}}{2 \sigma_{y}^{2}}\right) c = \left(\frac{\sin^{2}{\left (\theta \right )}}{2 \sigma_{x}^{2}} + \frac{\cos^{2}{\left (\theta \right )}}{2 \sigma_{y}^{2}}\right) If using a ``cov_matrix``, the model is of the form: .. math:: f(x, y) = A e^{-0.5 \left(\vec{x} - \vec{x}_{0}\right)^{T} \Sigma^{-1} \left(\vec{x} - \vec{x}_{0}\right)} where :math:`\vec{x} = [x, y]`, :math:`\vec{x}_{0} = [x_{0}, y_{0}]`, and :math:`\Sigma` is the covariance matrix: .. math:: \Sigma = \left(\begin{array}{ccc} \sigma_x^2 & \rho \sigma_x \sigma_y \\ \rho \sigma_x \sigma_y & \sigma_y^2 \end{array}\right) :math:`\rho` is the correlation between ``x`` and ``y``, which should be between -1 and +1. Positive correlation corresponds to a ``theta`` in the range 0 to 90 degrees. Negative correlation corresponds to a ``theta`` in the range of 0 to -90 degrees. See [1]_ for more details about the 2D Gaussian function. See Also -------- Gaussian1D, Box2D, Moffat2D References ---------- .. [1] https://en.wikipedia.org/wiki/Gaussian_function r)r(rgNc s|dkr@|dkr|jjj}|dkr,|jjj}|dkr|jjj}n|dk sX|dk sX|dk rbtdn^t|}|jdkr~t dtj |\} } t | \}}| dddf} t | d| d}|di|ddtdf|ddtdftjf||||||d |dS) Nz3Cannot specify both cov_matrix and x/y_stddev/theta)r'r'zCovariance matrix must be 2x2rrr)x_stddevy_stddev)r9x_meany_meanrKrLtheta) __class__rKr(rLrOrr6arrayshape ValueErrorZlinalgZeigsqrtZarctan2 setdefaultrHsuper__init__) r-r9rMrNrKrLrOZ cov_matrixkwargsZeig_valsZeig_vecsZy_vec)rPr1r2rWs,       zGaussian2D.__init__cCs |jtS)z)Gaussian full width at half maximum in X.)rKr4)r-r1r1r2x_fwhm>szGaussian2D.x_fwhmcCs |jtS)z)Gaussian full width at half maximum in Y.)rLr4)r-r1r1r2y_fwhmCszGaussian2D.y_fwhm@cCsT||j}||j}|jj}t|||\}}|j||j|f|j||j|ffS)a Tuple defining the default ``bounding_box`` limits in each dimension, ``((y_low, y_high), (x_low, x_high))`` The default offset from the mean is 5.5-sigma, corresponding to a relative error < 1e-7. The limits are adjusted for rotation. Parameters ---------- factor : float, optional The multiple of `x_stddev` and `y_stddev` used to define the limits. The default is 5.5. Examples -------- >>> from astropy.modeling.models import Gaussian2D >>> model = Gaussian2D(x_mean=0, y_mean=0, x_stddev=1, y_stddev=2) >>> model.bounding_box ((-11.0, 11.0), (-5.5, 5.5)) This range can be set directly (see: `Model.bounding_box `) or by using a different factor like: >>> model.bounding_box = model.bounding_box(factor=2) >>> model.bounding_box ((-4.0, 4.0), (-2.0, 2.0)) )rKrLrOvaluer rNrM)r-r.abrOr0dyr1r1r2r3Hs   zGaussian2D.bounding_boxcCst|d}t|d} td|} |d} |d} ||} ||}d|| | | }d| | | | }d| | || }|t|| d|| |||d S)z!Two dimensional Gaussian functionr'g@g?)r6cossinr7)r8r@r9rMrNrKrLrOcost2sint2sin2txstd2ystd2xdiffydiffr]r^cr1r1r2r:nszGaussian2D.evaluatec)Cst|}t|} t|d} t|d} td|} td|} |d}|d}|d}|d}||}||}|d}|d}d| || |}d| || |}d| || |}|t||||||| }| |d|d|}| |}| |}| || |}| |}| |}| } | |}!| |}"||}#|d||||}$|||d||}%|||||||!| }&|||||||"| }'||||||| | }(|#|$|%|&|'|(gS)zGTwo dimensional Gaussian function derivative with respect to parametersr'g@r;g?g?)r6r`rar7))r8r@r9rMrNrKrLrOcostsintrbrcZcos2trdrerfZxstd3Zystd3rgrhZxdiff2Zydiff2r]r^rigZ da_dthetaZ da_dx_stddevZ da_dy_stddevZ db_dthetaZ db_dx_stddevZ db_dy_stddevZ dc_dthetaZ dc_dx_stddevZ dc_dy_stddevZdg_dAZ dg_dx_meanZ dg_dy_meanZ dg_dx_stddevZ dg_dy_stddevZ dg_dthetar1r1r2r=sT        zGaussian2D.fit_derivcCs2|jjdkr|jjdkrdS|jj|jjdSdS)N)r8r@)rMr>rN)r-r1r1r2r?szGaussian2D.input_unitsc CsZ|d|dkrtdtd|dfd|dfd|dfd|dfdtjfd |d fgS) Nr8r@z(Units of 'x' and 'y' inputs should matchrMrNrKrLrOr9z)r rurad)r-rArBr1r1r2rCs    z*Gaussian2D._parameter_units_for_data_units)r[)rDrErFrGrr9rMrNrKrLrOr(rWrIrYrZr3rJr:r=r?rC __classcell__r1r1)rPr2rs"M      (   &  / c@sbeZdZdZdZdZeddZdZe ddZ e dd Z e d d Z e d d Ze ddZdS)r!zv Shift a coordinate. Parameters ---------- offset : float Offset to add to a coordinate. )r8r)r(TcCs |jjdkrdSd|jjiSdS)Nr8)offsetr>)r-r1r1r2r?s zShift.input_unitscCs|}|jd9_|S)z,One dimensional inverse Shift model function)copyrq)r-invr1r1r2inversesz Shift.inversecCs||S)z$One dimensional Shift model functionr1)r8rqr1r1r2r:szShift.evaluatecCs|S)z=Evaluate the implicit term (x) of one dimensional Shift modelr1)r8r1r1r2sum_of_implicit_termsszShift.sum_of_implicit_termscGst|}|gS)z@One dimensional Shift model derivative with respect to parameter)r6 ones_like)r8paramsZd_offsetr1r1r2r=s zShift.fit_derivN)rDrErFrGinputsoutputsrrqlinearrIr?rurJr:rvr=r1r1r1r2r!s     c@sbeZdZdZdZdZeddZdZdZ dZ dZ e ddZ e dd Zed d Zed d ZdS)ra  Multiply a model by a dimensionless factor. Parameters ---------- factor : float Factor by which to scale a coordinate. Notes ----- If ``factor`` is a `~astropy.units.Quantity` then the units will be stripped before the scaling operation. )r8r)r(TcCs |jjdkrdSd|jjiSdS)Nr8)r.r>)r-r1r1r2r?s zScale.input_unitscCs|}d|j|_|S)z,One dimensional inverse Scale model functionr)rsr.)r-rtr1r1r2rus z Scale.inversecCst|tjr|j}||S)z$One dimensional Scale model function) isinstancernr r\)r8r.r1r1r2r:s zScale.evaluatecGs |}|gS)z@One dimensional Scale model derivative with respect to parameterr1)r8rxd_factorr1r1r2r=%szScale.fit_derivN)rDrErFrGryrzrr.r{fittableZ_input_units_strictZ _input_units_allow_dimensionlessrIr?rurJr:r=r1r1r1r2rs    c@sNeZdZdZdZdZeddZdZdZ e ddZ e dd Z e d d Zd S) rz Multiply a model by a quantity or number. Parameters ---------- factor : float Factor by which to multiply a coordinate. )r8r)r(TcCs|}d|j|_|S)z/One dimensional inverse multiply model functionr)rsr.)r-rtr1r1r2ru>s zMultiply.inversecCs||S)z'One dimensional multiply model functionr1)r8r.r1r1r2r:EszMultiply.evaluatecGs |}|gS)zCOne dimensional multiply model derivative with respect to parameterr1)r8rxr}r1r1r2r=JszMultiply.fit_derivN)rDrErFrGryrzrr.r{r~rIrurJr:r=r1r1r1r2r-s   c@s@eZdZdZedddZeddZeddZe d d Z d S) rz One dimensional redshift scale factor model. Parameters ---------- z : float Redshift value. Notes ----- Model formula: .. math:: f(x) = x (1 + z) Zredshiftr)Z descriptionr(cCs d||S)z2One dimensional RedshiftScaleFactor model functionrr1)r8rmr1r1r2r:dszRedshiftScaleFactor.evaluatecCs |}|gS)z4One dimensional RedshiftScaleFactor model derivativer1)r8rmZd_zr1r1r2r=jszRedshiftScaleFactor.fit_derivcCs |}dd|jd|_|S)z!Inverse RedshiftScaleFactor modelg?)rsrm)r-rtr1r1r2ruqszRedshiftScaleFactor.inverseN) rDrErFrGrrmrJr:r=rIrur1r1r1r2rRs    c@sReZdZdZeddZeddZeddZdZe ddZ e dd Z d d Z dS) ra One dimensional Sersic surface brightness profile. Parameters ---------- amplitude : float Surface brightness at r_eff. r_eff : float Effective (half-light) radius n : float Sersic Index. See Also -------- Gaussian1D, Moffat1D, Lorentz1D Notes ----- Model formula: .. math:: I(r)=I_e\exp\left\{-b_n\left[\left(\frac{r}{r_{e}}\right)^{(1/n)}-1\right]\right\} The constant :math:`b_n` is defined such that :math:`r_e` contains half the total luminosity, and can be solved for numerically. .. math:: \Gamma(2n) = 2\gamma (b_n,2n) Examples -------- .. plot:: :include-source: import numpy as np from astropy.modeling.models import Sersic1D import matplotlib.pyplot as plt plt.figure() plt.subplot(111, xscale='log', yscale='log') s1 = Sersic1D(amplitude=1, r_eff=5) r=np.arange(0, 100, .01) for n in range(1, 10): s1.n = n plt.plot(r, s1(r), color=str(float(n) / 15)) plt.axis([1e-1, 30, 1e-2, 1e3]) plt.xlabel('log Radius') plt.ylabel('log Surface Brightness') plt.text(.25, 1.5, 'n=1') plt.text(.25, 300, 'n=10') plt.xticks([]) plt.yticks([]) plt.show() References ---------- .. [1] http://ned.ipac.caltech.edu/level5/March05/Graham/Graham2.html r)r(NcCsn|jdkr>yddlm}||_Wntk r<tdYnX|t|d|d ||d|dS)z(One dimensional Sersic profile function.Nr) gammaincinvz%Sersic1D model requires scipy > 0.11.r'g?r) _gammaincinv scipy.specialrrS ImportErrorr6r7)clsrr9r_effnrr1r1r2r:s   zSersic1D.evaluatecCs |jjdkrdSd|jjiSdS)Nr8)rr>)r-r1r1r2r?s zSersic1D.input_unitscCstd|dfd|dfgS)Nrr8r9r@)r)r-rArBr1r1r2rCs z(Sersic1D._parameter_units_for_data_units)rDrErFrGrr9rrr classmethodr:rIr?rCr1r1r1r2rzs>     c@sZeZdZdZeddZeddZeddZeddZ eddZ e d d Z d d Z d S)r"aL One dimensional Sine model. Parameters ---------- amplitude : float Oscillation amplitude frequency : float Oscillation frequency phase : float Oscillation phase See Also -------- Const1D, Linear1D Notes ----- Model formula: .. math:: f(x) = A \sin(2 \pi f x + 2 \pi p) Examples -------- .. plot:: :include-source: import numpy as np import matplotlib.pyplot as plt from astropy.modeling.models import Sine1D plt.figure() s1 = Sine1D(amplitude=1, frequency=.25) r=np.arange(0, 10, .01) for amplitude in range(1,4): s1.amplitude = amplitude plt.plot(r, s1(r), color=str(0.25 * amplitude), lw=2) plt.axis([0, 10, -5, 5]) plt.show() r)r(rcCs.t|||}t|tr |j}|t|S)z#One dimensional Sine model function)TWOPIr|r r\r6ra)r8r9 frequencyphaseZargumentr1r1r2r: s zSine1D.evaluatecCsltt||t|}t||tt||t|}t|tt||t|}|||gS)z%One dimensional Sine model derivative)r6rarr`)r8r9rrr<Z d_frequencyZd_phaser1r1r2r=s  zSine1D.fit_derivcCs$|jjdkrdSdd|jjiSdS)Nr8g?)rr>)r-r1r1r2r?#s zSine1D.input_unitscCs td|ddfd|dfgS)Nrr8rrr9r@)r)r-rArBr1r1r2rC*sz&Sine1D._parameter_units_for_data_unitsN)rDrErFrGrr9rrrJr:r=rIr?rCr1r1r1r2r"s,    c@s`eZdZdZeddZeddZdZeddZ edd Z e d d Z e d d Z ddZdS)ra( One dimensional Line model. Parameters ---------- slope : float Slope of the straight line intercept : float Intercept of the straight line See Also -------- Const1D Notes ----- Model formula: .. math:: f(x) = a x + b r)r(rTcCs |||S)z#One dimensional Line model functionr1)r8slope interceptr1r1r2r:JszLinear1D.evaluatecCs|}t|}||gS)z@One dimensional Line model derivative with respect to parameters)r6rw)r8rrZd_slope d_interceptr1r1r2r=Ps zLinear1D.fit_derivcCs&|jd}|j |j}|j||dS)Nrr)rr)rrrP)r-Z new_slopeZ new_interceptr1r1r2ruXs zLinear1D.inversecCs4|jjdkr|jjdkrdSd|jj|jjiSdS)Nr8)rr>r)r-r1r1r2r?^szLinear1D.input_unitscCs$td|dfd|d|dfgS)Nrr@rr8)r)r-rArBr1r1r2rCes z(Linear1D._parameter_units_for_data_unitsN)rDrErFrGrrrr{rJr:r=rIrur?rCr1r1r1r2r/s      c@sJeZdZdZeddZeddZeddZdZe ddZ e dd Z d S) Planar2Da Two dimensional Plane model. Parameters ---------- slope_x : float Slope of the straight line in X slope_y : float Slope of the straight line in Y intercept : float Z-intercept of the straight line See Also -------- Linear1D, Polynomial2D Notes ----- Model formula: .. math:: f(x, y) = a x + b y + c r)r(rTcCs|||||S)z$Two dimensional Plane model functionr1)r8r@slope_xslope_yrr1r1r2r:szPlanar2D.evaluatecCs|}|}t|}|||gS)zATwo dimensional Plane model derivative with respect to parameters)r6rw)r8r@rrrZ d_slope_xZ d_slope_yrr1r1r2r=s zPlanar2D.fit_derivN) rDrErFrGrrrrr{rJr:r=r1r1r1r2rjs    rc@sdeZdZdZeddZeddZeddZeddZ eddZ dd d Z e d d Z ddZdS)ra_ One dimensional Lorentzian model. Parameters ---------- amplitude : float Peak value x_0 : float Position of the peak fwhm : float Full width at half maximum See Also -------- Gaussian1D, Box1D, MexicanHat1D Notes ----- Model formula: .. math:: f(x) = \frac{A \gamma^{2}}{\gamma^{2} + \left(x - x_{0}\right)^{2}} Examples -------- .. plot:: :include-source: import numpy as np import matplotlib.pyplot as plt from astropy.modeling.models import Lorentz1D plt.figure() s1 = Lorentz1D() r = np.arange(-5, 5, .01) for factor in range(1, 4): s1.amplitude = factor plt.plot(r, s1(r), color=str(0.25 * factor), lw=2) plt.axis([-5, 5, -1, 4]) plt.show() r)r(rcCs(||dd||d|ddS)z)One dimensional Lorentzian model functiong@r'r1)r8r9x_0r5r1r1r2r:szLorentz1D.evaluatecCsj|d|d||d}||d|d||d||d}d|||d|}|||gS)zFOne dimensional Lorentzian model derivative with respect to parametersr'rr1)r8r9rr5r<d_x_0Zd_fwhmr1r1r2r=s zLorentz1D.fit_derivcCs |j}||j}||||fS)arTuple defining the default ``bounding_box`` limits, ``(x_low, x_high)``. Parameters ---------- factor : float The multiple of FWHM used to define the limits. Default is chosen to include most (99%) of the area under the curve, while still showing the central feature of interest. )rr5)r-r.r/r0r1r1r2r3s  zLorentz1D.bounding_boxcCs |jjdkrdSd|jjiSdS)Nr8)rr>)r-r1r1r2r?s zLorentz1D.input_unitscCs&td|dfd|dfd|dfgS)Nrr8r5r9r@)r)r-rArBr1r1r2rCs  z)Lorentz1D._parameter_units_for_data_unitsN)r)rDrErFrGrr9rr5rJr:r=r3rIr?rCr1r1r1r2rs-      c @seZdZdZeddZeddZedejdZ ee ddZ e ddddgdd d d gd d d d gddddggZ eddZeddZeddZddZdS)r&a! One dimensional model for the Voigt profile. Parameters ---------- x_0 : float Position of the peak amplitude_L : float The Lorentzian amplitude fwhm_L : float The Lorentzian full width at half maximum fwhm_G : float The Gaussian full width at half maximum See Also -------- Gaussian1D, Lorentz1D Notes ----- Algorithm for the computation taken from McLean, A. B., Mitchell, C. E. J. & Swanston, D. M. Implementation of an efficient analytical approximation to the Voigt function for photoemission lineshape analysis. Journal of Electron Spectroscopy and Related Phenomena 69, 125-132 (1994) Examples -------- .. plot:: :include-source: import numpy as np from astropy.modeling.models import Voigt1D import matplotlib.pyplot as plt plt.figure() x = np.arange(0, 10, 0.01) v1 = Voigt1D(x_0=5, amplitude_L=10, fwhm_L=0.5, fwhm_G=0.9) plt.plot(x, v1(x)) plt.show() r)r(rr'gq= ףpgQIg??g~:?g?g~:ؿgX9vӿg.1?g/$?g~k g/$g~k ?cCs|j\}}}} ttd} ||d| |} t| dtjf} || |} t| dtjf} tj|| || | || |d| |ddd} ||ttj| || S)Nr'.rr)axis)_abcdr6rTlog atleast_1dnewaxissumpi)rr8r amplitude_Lfwhm_Lfwhm_GABCDsqrt_ln2XYVr1r1r2r:2s :zVoigt1D.evaluatec Cs|j\}}}} ttd} ||d| |} t| ddtjf} || |} t| ddtjf} ||ttj| |} || || | |}| |d| |d}tj||dd}tj| |d| ||t|dd}tj||d| ||t|dd}| |d| || ||| |||| || || |d||||||g}|S)Nr'rr)r) rr6rTrrrrrZsquare)rr8rrrrrrrrrrrZconstantalphaZbetarZdVdxZdVdyZdydar1r1r2r=@s" ,, 0zVoigt1D.fit_derivcCs |jjdkrdSd|jjiSdS)Nr8)rr>)r-r1r1r2r?Ws zVoigt1D.input_unitscCs0td|dfd|dfd|dfd|dfgS)Nrr8rrrr@)r)r-rArBr1r1r2rC^s   z'Voigt1D._parameter_units_for_data_unitsN)rDrErFrGrrrr6rrrrrQrrr:r=rIr?rCr1r1r1r2r&s)        c@sJeZdZdZeddZdZeddZeddZ e d d Z d d Z d S)ra One dimensional Constant model. Parameters ---------- amplitude : float Value of the constant function See Also -------- Const2D Notes ----- Model formula: .. math:: f(x) = A Examples -------- .. plot:: :include-source: import numpy as np import matplotlib.pyplot as plt from astropy.modeling.models import Const1D plt.figure() s1 = Const1D() r = np.arange(-5, 5, .01) for factor in range(1, 4): s1.amplitude = factor plt.plot(r, s1(r), color=str(0.25 * factor), lw=2) plt.axis([-5, 5, -1, 4]) plt.show() r)r(TcCs\|jdkr(tj|dd}||n|tj|dd}t|trTt||jddS|SdS)z'One dimensional Constant model functionrF)subok)r>rsN) sizer6 empty_likefillitemrwr|r r>)r8r9r1r1r2r:s  zConst1D.evaluatecCst|}|gS)zDOne dimensional Constant model derivative with respect to parameters)r6rw)r8r9r<r1r1r2r=s zConst1D.fit_derivcCsdS)Nr1)r-r1r1r2r?szConst1D.input_unitscCstd|dfgS)Nr9r@)r)r-rArBr1r1r2rCsz'Const1D._parameter_units_for_data_unitsN) rDrErFrGrr9r{rJr:r=rIr?rCr1r1r1r2res'    c@s>eZdZdZeddZdZeddZe ddZ d d Z d S) rz Two dimensional Constant model. Parameters ---------- amplitude : float Value of the constant function See Also -------- Const1D Notes ----- Model formula: .. math:: f(x, y) = A r)r(TcCs\|jdkr(tj|dd}||n|tj|dd}t|trTt||jddS|SdS)z'Two dimensional Constant model functionrF)r)r>rsN) rr6rrrrwr|r r>)r8r@r9r1r1r2r:s  zConst2D.evaluatecCsdS)Nr1)r-r1r1r2r?szConst2D.input_unitscCstd|dfgS)Nr9rm)r)r-rArBr1r1r2rCsz'Const2D._parameter_units_for_data_unitsN) rDrErFrGrr9r{rJr:rIr?rCr1r1r1r2rs    c@sxeZdZdZeddZeddZeddZeddZeddZ eddZ e ddZ e ddZe d d Zd d Zd S)raA A 2D Ellipse model. Parameters ---------- amplitude : float Value of the ellipse. x_0 : float x position of the center of the disk. y_0 : float y position of the center of the disk. a : float The length of the semimajor axis. b : float The length of the semiminor axis. theta : float The rotation angle in radians of the semimajor axis. The rotation angle increases counterclockwise from the positive x axis. See Also -------- Disk2D, Box2D Notes ----- Model formula: .. math:: f(x, y) = \left \{ \begin{array}{ll} \mathrm{amplitude} & : \left[\frac{(x - x_0) \cos \theta + (y - y_0) \sin \theta}{a}\right]^2 + \left[\frac{-(x - x_0) \sin \theta + (y - y_0) \cos \theta}{b}\right]^2 \leq 1 \\ 0 & : \mathrm{otherwise} \end{array} \right. Examples -------- .. plot:: :include-source: import numpy as np from astropy.modeling.models import Ellipse2D from astropy.coordinates import Angle import matplotlib.pyplot as plt import matplotlib.patches as mpatches x0, y0 = 25, 25 a, b = 20, 10 theta = Angle(30, 'deg') e = Ellipse2D(amplitude=100., x_0=x0, y_0=y0, a=a, b=b, theta=theta.radian) y, x = np.mgrid[0:50, 0:50] fig, ax = plt.subplots(1, 1) ax.imshow(e(x, y), origin='lower', interpolation='none', cmap='Greys_r') e2 = mpatches.Ellipse((x0, y0), 2*a, 2*b, theta.degree, edgecolor='red', facecolor='none') ax.add_patch(e2) plt.show() r)r(rcCs||}||} t|} t|} || | | } ||  | | } | |d| |ddk}t|g|g}t|trt||jddS|SdS)z'Two dimensional Ellipse model function.r'g?F)r>rsN)r6r`raselectr|r r>)r8r@r9ry_0r]r^rOZxxZyyrjrkZ numerator1Z numerator2Z in_ellipseresultr1r1r2r:0s   zEllipse2D.evaluatecCsL|j}|j}|jj}t|||\}}|j||j|f|j||j|ffS)zu Tuple defining the default ``bounding_box`` limits. ``((y_low, y_high), (x_low, x_high))`` )r]r^rOr\r rr)r-r]r^rOr0r_r1r1r2r3Bs zEllipse2D.bounding_boxcCs&|jjdkrdS|jj|jjdSdS)N)r8r@)rr>r)r-r1r1r2r?Rs zEllipse2D.input_unitsc CsZ|d|dkrtdtd|dfd|dfd|dfd|dfdtjfd |d fgS) Nr8r@z(Units of 'x' and 'y' inputs should matchrrr]r^rOr9rm)r rrnro)r-rArBr1r1r2rCZs    z)Ellipse2D._parameter_units_for_data_unitsN)rDrErFrGrr9rrr]r^rOrJr:rIr3r?rCr1r1r1r2rsD         c@sdeZdZdZeddZeddZeddZeddZe ddZ e ddZ e d d Z d d Zd S)raf Two dimensional radial symmetric Disk model. Parameters ---------- amplitude : float Value of the disk function x_0 : float x position center of the disk y_0 : float y position center of the disk R_0 : float Radius of the disk See Also -------- Box2D, TrapezoidDisk2D Notes ----- Model formula: .. math:: f(r) = \left \{ \begin{array}{ll} A & : r \leq R_0 \\ 0 & : r > R_0 \end{array} \right. r)r(rcCsR||d||d}t||dkg|g}t|trJt||jddS|SdS)z#Two dimensional Disk model functionr'F)r>rsN)r6rr|r r>)r8r@r9rrR_0rrrr1r1r2r:s  zDisk2D.evaluatecCs0|j|j|j|jf|j|j|j|jffS)zu Tuple defining the default ``bounding_box`` limits. ``((y_low, y_high), (x_low, x_high))`` )rrr)r-r1r1r2r3szDisk2D.bounding_boxcCs2|jjdkr|jjdkrdS|jj|jjdSdS)N)r8r@)rr>r)r-r1r1r2r?szDisk2D.input_unitscCsH|d|dkrtdtd|dfd|dfd|dfd|dfgS) Nr8r@z(Units of 'x' and 'y' inputs should matchrrrr9rm)r r)r-rArBr1r1r2rCs    z&Disk2D._parameter_units_for_data_unitsN)rDrErFrGrr9rrrrJr:rIr3r?rCr1r1r1r2rhs     cseZdZdZeddZeddZeddZeddZeddZ ej ej ej ej e j dffdd Z e dd Z ed d Zed d ZddZZS)r%a7 Two dimensional radial symmetric Ring model. Parameters ---------- amplitude : float Value of the disk function x_0 : float x position center of the disk y_0 : float y position center of the disk r_in : float Inner radius of the ring width : float Width of the ring. r_out : float Outer Radius of the ring. Can be specified instead of width. See Also -------- Disk2D, TrapezoidDisk2D Notes ----- Model formula: .. math:: f(r) = \left \{ \begin{array}{ll} A & : r_{in} \leq r \leq r_{out} \\ 0 & : \text{else} \end{array} \right. Where :math:`r_{out} = r_{in} + r_{width}`. r)r(rNc sX|dk r&||jjkrtd||}n|dkr6|jj}tjf|||||d|dS)Nz6Cannot specify both width and outer radius separately.)r9rrr_inwidth)rr(rrVrW)r-r9rrrrZr_outrX)rPr1r2rWs  zRing2D.__init__c Csj||d||d}t||dk|||dk}t|g|g} t|trbt| |jddS| SdS)z$Two dimensional Ring model function.r'F)r>rsN)r6 logical_andrr|r r>) r8r@r9rrrrrZr_rangerr1r1r2r:s   zRing2D.evaluatecCs4|j|j}|j||j|f|j||j|ffS)zn Tuple defining the default ``bounding_box``. ``((y_low, y_high), (x_low, x_high))`` )rrrr)r-drr1r1r2r3s zRing2D.bounding_boxcCs&|jjdkrdS|jj|jjdSdS)N)r8r@)rr>r)r-r1r1r2r?s zRing2D.input_unitscCsR|d|dkrtdtd|dfd|dfd|dfd|dfd|d fgS) Nr8r@z(Units of 'x' and 'y' inputs should matchrrrrr9rm)r r)r-rArBr1r1r2rCs    z&Ring2D._parameter_units_for_data_units)rDrErFrGrr9rrrrr(rWrJr:rIr3r?rCrpr1r1)rPr2r%s%       c@seZdZdZddZdS)Delta1Dz%One dimensional Dirac delta function.cCs tddS)NzNot implemented)r)r-r1r1r2rW(szDelta1D.__init__N)rDrErFrGrWr1r1r1r2r%src@seZdZdZddZdS)Delta2Dz%Two dimensional Dirac delta function.cCs tddS)NzNot implemented)r)r-r1r1r2rW/szDelta2D.__init__N)rDrErFrGrWr1r1r1r2r,src@sZeZdZdZeddZeddZeddZeddZ e ddZ e d d Z d d Z d S)ra One dimensional Box model. Parameters ---------- amplitude : float Amplitude A x_0 : float Position of the center of the box function width : float Width of the box See Also -------- Box2D, TrapezoidDisk2D Notes ----- Model formula: .. math:: f(x) = \left \{ \begin{array}{ll} A & : x_0 - w/2 \leq x \leq x_0 + w/2 \\ 0 & : \text{else} \end{array} \right. Examples -------- .. plot:: :include-source: import numpy as np import matplotlib.pyplot as plt from astropy.modeling.models import Box1D plt.figure() s1 = Box1D() r = np.arange(-5, 5, .01) for factor in range(1, 4): s1.amplitude = factor s1.width = factor plt.plot(r, s1(r), color=str(0.25 * factor), lw=2) plt.axis([-5, 5, -1, 4]) plt.show() r)r(rcCsXt|||dk|||dk}t|g|gd}t|trPt||jddS|SdS)z"One dimensional Box model functiong@rF)r>rsN)r6rrr|r r>)r8r9rrZinsiderr1r1r2r:ls $ zBox1D.evaluatecCs|jd}|j||j|fS)zc Tuple defining the default ``bounding_box`` limits. ``(x_low, x_high))`` r')rr)r-r0r1r1r2r3xs zBox1D.bounding_boxcCs |jjdkrdSd|jjiSdS)Nr8)rr>)r-r1r1r2r?s zBox1D.input_unitscCs&td|dfd|dfd|dfgS)Nrr8rr9r@)r)r-rArBr1r1r2rCs  z%Box1D._parameter_units_for_data_unitsN)rDrErFrGrr9rrrJr:rIr3r?rCr1r1r1r2r3s3    c@sneZdZdZeddZeddZeddZeddZeddZ e ddZ e ddZ e d d Zd d Zd S)ra Two dimensional Box model. Parameters ---------- amplitude : float Amplitude A x_0 : float x position of the center of the box function x_width : float Width in x direction of the box y_0 : float y position of the center of the box function y_width : float Width in y direction of the box See Also -------- Box1D, Gaussian2D, Moffat2D Notes ----- Model formula: .. math:: f(x, y) = \left \{ \begin{array}{ll} A : & x_0 - w_x/2 \leq x \leq x_0 + w_x/2 \text{ and} \\ & y_0 - w_y/2 \leq y \leq y_0 + w_y/2 \\ 0 : & \text{else} \end{array} \right. r)r(rc Cst|||dk|||dk}t|||dk|||dk}tt||g|gd} t|tr|t| |jddS| SdS)z"Two dimensional Box model functiong@rF)r>rsN)r6rrr|r r>) r8r@r9rrx_widthy_widthZx_rangeZy_rangerr1r1r2r:s zBox2D.evaluatecCs<|jd}|jd}|j||j|f|j||j|ffS)zn Tuple defining the default ``bounding_box``. ``((y_low, y_high), (x_low, x_high))`` r')rrrr)r-r0r_r1r1r2r3s  zBox2D.bounding_boxcCs&|jjdkrdS|jj|jjdSdS)N)r8r@)rr>r)r-r1r1r2r?s zBox2D.input_unitscCs:td|dfd|dfd|dfd|dfd|dfgS) Nrr8rr@rrr9rm)r)r-rArBr1r1r2rCs     z%Box2D._parameter_units_for_data_unitsN)rDrErFrGrr9rrrrrJr:rIr3r?rCr1r1r1r2rs#        c@sdeZdZdZeddZeddZeddZeddZe ddZ e ddZ e d d Z d d Zd S)r#ae One dimensional Trapezoid model. Parameters ---------- amplitude : float Amplitude of the trapezoid x_0 : float Center position of the trapezoid width : float Width of the constant part of the trapezoid. slope : float Slope of the tails of the trapezoid See Also -------- Box1D, Gaussian1D, Moffat1D Examples -------- .. plot:: :include-source: import numpy as np import matplotlib.pyplot as plt from astropy.modeling.models import Trapezoid1D plt.figure() s1 = Trapezoid1D() r = np.arange(-5, 5, .01) for factor in range(1, 4): s1.amplitude = factor s1.width = factor plt.plot(r, s1(r), color=str(0.25 * factor), lw=2) plt.axis([-5, 5, -1, 4]) plt.show() r)r(rcCs||d}||d}|||}|||}t||k||k} t||k||k} t||k||k} |||} |} |||}t| | | g| | |g}t|trt||jddS|SdS)z(One dimensional Trapezoid model functiong@F)r>rsN)r6rrr|r r>)r8r9rrrZx2Zx3Zx1Zx4Zrange_aZrange_bZrange_cZval_aZval_bZval_crr1r1r2r:s       zTrapezoid1D.evaluatecCs*|jd|j|j}|j||j|fS)zc Tuple defining the default ``bounding_box`` limits. ``(x_low, x_high))`` r')rr9rr)r-r0r1r1r2r32szTrapezoid1D.bounding_boxcCs |jjdkrdSd|jjiSdS)Nr8)rr>)r-r1r1r2r?>s zTrapezoid1D.input_unitscCs8td|dfd|dfd|d|dfd|dfgS)Nrr8rrr@r9)r)r-rArBr1r1r2rCEs  z+Trapezoid1D._parameter_units_for_data_unitsN)rDrErFrGrr9rrrrJr:rIr3r?rCr1r1r1r2r#s(      c@sneZdZdZeddZeddZeddZeddZeddZ e ddZ e ddZ e d d Zd d Zd S)r$a Two dimensional circular Trapezoid model. Parameters ---------- amplitude : float Amplitude of the trapezoid x_0 : float x position of the center of the trapezoid y_0 : float y position of the center of the trapezoid R_0 : float Radius of the constant part of the trapezoid. slope : float Slope of the tails of the trapezoid in x direction. See Also -------- Disk2D, Box2D r)r(rc Cst||d||d}||k}t||k||||k} |} ||||} t|| g| | g} t|trt| |jddS| SdS)z-Two dimensional Trapezoid Disk model functionr'F)r>rsN)r6rTrrr|r r>) r8r@r9rrrrrZrange_1Zrange_2Zval_1Zval_2rr1r1r2r:hs zTrapezoidDisk2D.evaluatecCs:|j|j|j}|j||j|f|j||j|ffS)zn Tuple defining the default ``bounding_box``. ``((y_low, y_high), (x_low, x_high))`` )rr9rrr)r-rr1r1r2r3xszTrapezoidDisk2D.bounding_boxcCs2|jjdkr|jjdkrdS|jj|jjdSdS)N)r8r@)rr>r)r-r1r1r2r?szTrapezoidDisk2D.input_unitscCsZ|d|dkrtdtd|dfd|dfd|dfd|d|dfd |dfgS) Nr8r@z(Units of 'x' and 'y' inputs should matchrrrrrmr9)r r)r-rArBr1r1r2rCs   z/TrapezoidDisk2D._parameter_units_for_data_unitsN)rDrErFrGrr9rrrrrJr:rIr3r?rCr1r1r1r2r$Ls       c@sXeZdZdZeddZeddZeddZeddZ ddd Z e d d Z d d Z dS)ra One dimensional Mexican Hat model. Parameters ---------- amplitude : float Amplitude x_0 : float Position of the peak sigma : float Width of the Mexican hat See Also -------- MexicanHat2D, Box1D, Gaussian1D, Trapezoid1D Notes ----- Model formula: .. math:: f(x) = {A \left(1 - \frac{\left(x - x_{0}\right)^{2}}{\sigma^{2}}\right) e^{- \frac{\left(x - x_{0}\right)^{2}}{2 \sigma^{2}}}} Examples -------- .. plot:: :include-source: import numpy as np import matplotlib.pyplot as plt from astropy.modeling.models import MexicanHat1D plt.figure() s1 = MexicanHat1D() r = np.arange(-5, 5, .01) for factor in range(1, 4): s1.amplitude = factor s1.width = factor plt.plot(r, s1(r), color=str(0.25 * factor), lw=2) plt.axis([-5, 5, -2, 4]) plt.show() r)r(rcCs4||dd|d}|dd|t| S)z*One dimensional Mexican Hat model functionr'r)r6r7)r8r9rsigmaZxx_wwr1r1r2r:szMexicanHat1D.evaluate$@cCs |j}||j}||||fS)zTuple defining the default ``bounding_box`` limits, ``(x_low, x_high)``. Parameters ---------- factor : float The multiple of sigma used to define the limits. )rr)r-r.r/r0r1r1r2r3s  zMexicanHat1D.bounding_boxcCs |jjdkrdSd|jjiSdS)Nr8)rr>)r-r1r1r2r?s zMexicanHat1D.input_unitscCs&td|dfd|dfd|dfgS)Nrr8rr9r@)r)r-rArBr1r1r2rCs  z,MexicanHat1D._parameter_units_for_data_unitsN)r)rDrErFrGrr9rrrJr:r3rIr?rCr1r1r1r2rs/      c@sXeZdZdZeddZeddZeddZeddZe ddZ e ddZ d d Z d S) raY Two dimensional symmetric Mexican Hat model. Parameters ---------- amplitude : float Amplitude x_0 : float x position of the peak y_0 : float y position of the peak sigma : float Width of the Mexican hat See Also -------- MexicanHat1D, Gaussian2D Notes ----- Model formula: .. math:: f(x, y) = A \left(1 - \frac{\left(x - x_{0}\right)^{2} + \left(y - y_{0}\right)^{2}}{\sigma^{2}}\right) e^{\frac{- \left(x - x_{0}\right)^{2} - \left(y - y_{0}\right)^{2}}{2 \sigma^{2}}} r)r(rcCs<||d||dd|d}|d|t| S)z*Two dimensional Mexican Hat model functionr'r)r6r7)r8r@r9rrrZrr_wwr1r1r2r:s$zMexicanHat2D.evaluatecCs&|jjdkrdS|jj|jjdSdS)N)r8r@)rr>r)r-r1r1r2r?s zMexicanHat2D.input_unitscCsH|d|dkrtdtd|dfd|dfd|dfd|dfgS) Nr8r@z(Units of 'x' and 'y' inputs should matchrrrr9rm)r r)r-rArBr1r1r2rC%s    z,MexicanHat2D._parameter_units_for_data_unitsN)rDrErFrGrr9rrrrJr:rIr?rCr1r1r1r2rs      c@s`eZdZdZeddZeddZeddZeddZdZ dZ e ddZ e dd Zd d ZdS) r a Two dimensional Airy disk model. Parameters ---------- amplitude : float Amplitude of the Airy function. x_0 : float x position of the maximum of the Airy function. y_0 : float y position of the maximum of the Airy function. radius : float The radius of the Airy disk (radius of the first zero). See Also -------- Box2D, TrapezoidDisk2D, Gaussian2D Notes ----- Model formula: .. math:: f(r) = A \left[\frac{2 J_1(\frac{\pi r}{R/R_z})}{\frac{\pi r}{R/R_z}}\right]^2 Where :math:`J_1` is the first order Bessel function of the first kind, :math:`r` is radial distance from the maximum of the Airy function (:math:`r = \sqrt{(x - x_0)^2 + (y - y_0)^2}`), :math:`R` is the input ``radius`` parameter, and :math:`R_z = 1.2196698912665045`). For an optical system, the radius of the first zero represents the limiting angular resolution and is approximately 1.22 * lambda / D, where lambda is the wavelength of the light and D is the diameter of the aperture. See [1]_ for more details about the Airy disk. References ---------- .. [1] https://en.wikipedia.org/wiki/Airy_disk r)r(rNc Cs|jdkrXy0ddlm}m}|dddtj|_||_Wntk rVtdYnXt ||d||d||j} t | t r| t j} t| j} tj| | dk} d|| | d| | dk<t |t rt | t jdd } | |9} | S) z#Two dimensional Airy model functionNr)j1jn_zerosrz'AiryDisk2D model requires scipy > 0.11.r'g@F)rs)_rzrrrr6r_j1rSrrTr|r Zto_valuernZdimensionless_unscaledZonesrR) rr8r@r9rrradiusrrrrmZrtr1r1r2r:cs"  (    zAiryDisk2D.evaluatecCs&|jjdkrdS|jj|jjdSdS)N)r8r@)rr>r)r-r1r1r2r?s zAiryDisk2D.input_unitscCsH|d|dkrtdtd|dfd|dfd|dfd|dfgS) Nr8r@z(Units of 'x' and 'y' inputs should matchrrrr9rm)r r)r-rArBr1r1r2rCs    z*AiryDisk2D._parameter_units_for_data_units)rDrErFrGrr9rrrrrrr:rIr?rCr1r1r1r2r 1s)      c@speZdZdZeddZeddZeddZeddZe ddZ e ddZ e d d Z e d d Zd dZdS)ra One dimensional Moffat model. Parameters ---------- amplitude : float Amplitude of the model. x_0 : float x position of the maximum of the Moffat model. gamma : float Core width of the Moffat model. alpha : float Power index of the Moffat model. See Also -------- Gaussian1D, Box1D Notes ----- Model formula: .. math:: f(x) = A \left(1 + \frac{\left(x - x_{0}\right)^{2}}{\gamma^{2}}\right)^{- \alpha} Examples -------- .. plot:: :include-source: import numpy as np import matplotlib.pyplot as plt from astropy.modeling.models import Moffat1D plt.figure() s1 = Moffat1D() r = np.arange(-5, 5, .01) for factor in range(1, 4): s1.amplitude = factor s1.width = factor plt.plot(r, s1(r), color=str(0.25 * factor), lw=2) plt.axis([-5, 5, -1, 4]) plt.show() r)r(rcCs"d|jtdd|jdS)z Moffat full width at half maximum. Derivation of the formula is available in `this notebook by Yoonsoo Bach `_. g@g?)gammar6rTr)r-r1r1r2r5sz Moffat1D.fwhmcCs|d|||d| S)z%One dimensional Moffat model functionrr'r1)r8r9rrrr1r1r2r:szMoffat1D.evaluatec Csd||d|d}|| }d|||||||d}d||||d|||d}| |t|} |||| gS)zBOne dimensional Moffat model derivative with respect to parametersrr'r;)r6r) r8r9rrrZfacd_Ard_gammad_alphar1r1r2r=s $zMoffat1D.fit_derivcCs |jjdkrdSd|jjiSdS)Nr8)rr>)r-r1r1r2r?s zMoffat1D.input_unitscCs&td|dfd|dfd|dfgS)Nrr8rr9r@)r)r-rArBr1r1r2rCs  z(Moffat1D._parameter_units_for_data_unitsN)rDrErFrGrr9rrrrIr5rJr:r=r?rCr1r1r1r2rs0      c@szeZdZdZeddZeddZeddZeddZeddZ e ddZ e ddZ e d d Ze d d Zd dZdS)rak Two dimensional Moffat model. Parameters ---------- amplitude : float Amplitude of the model. x_0 : float x position of the maximum of the Moffat model. y_0 : float y position of the maximum of the Moffat model. gamma : float Core width of the Moffat model. alpha : float Power index of the Moffat model. See Also -------- Gaussian2D, Box2D Notes ----- Model formula: .. math:: f(x, y) = A \left(1 + \frac{\left(x - x_{0}\right)^{2} + \left(y - y_{0}\right)^{2}}{\gamma^{2}}\right)^{- \alpha} r)r(rcCs"d|jtdd|jdS)z Moffat full width at half maximum. Derivation of the formula is available in `this notebook by Yoonsoo Bach `_. g@g?)rr6rTr)r-r1r1r2r5 sz Moffat2D.fwhmcCs2||d||d|d}|d|| S)z%Two dimensional Moffat model functionr'rr1)r8r@r9rrrrrr_ggr1r1r2r:# s zMoffat2D.evaluatec Cs||d||d|d}d|| }d||||||dd|} d||||||dd|} | |td|} d|||||dd|} || | | | gS)zBTwo dimensional Moffat model derivative with respect to parametersr'rr;)r6r) r8r@r9rrrrrrrZd_y_0rrr1r1r2r=* s zMoffat2D.fit_derivcCs&|jjdkrdS|jj|jjdSdS)N)r8r@)rr>r)r-r1r1r2r?9 s zMoffat2D.input_unitscCsH|d|dkrtdtd|dfd|dfd|dfd|dfgS) Nr8r@z(Units of 'x' and 'y' inputs should matchrrrr9rm)r r)r-rArBr1r1r2rCA s    z(Moffat2D._parameter_units_for_data_unitsN)rDrErFrGrr9rrrrrIr5rJr:r=r?rCr1r1r1r2rs        c@szeZdZdZeddZeddZeddZeddZeddZ eddZ eddZ dZ e ddZed d Zd d ZdS) r a Two dimensional Sersic surface brightness profile. Parameters ---------- amplitude : float Surface brightness at r_eff. r_eff : float Effective (half-light) radius n : float Sersic Index. x_0 : float, optional x position of the center. y_0 : float, optional y position of the center. ellip : float, optional Ellipticity. theta : float, optional Rotation angle in radians, counterclockwise from the positive x-axis. See Also -------- Gaussian2D, Moffat2D Notes ----- Model formula: .. math:: I(x,y) = I(r) = I_e\exp\left\{-b_n\left[\left(\frac{r}{r_{e}}\right)^{(1/n)}-1\right]\right\} The constant :math:`b_n` is defined such that :math:`r_e` contains half the total luminosity, and can be solved for numerically. .. math:: \Gamma(2n) = 2\gamma (b_n,2n) Examples -------- .. plot:: :include-source: import numpy as np from astropy.modeling.models import Sersic2D import matplotlib.pyplot as plt x,y = np.meshgrid(np.arange(100), np.arange(100)) mod = Sersic2D(amplitude = 1, r_eff = 25, n=4, x_0=50, y_0=50, ellip=.5, theta=-1) img = mod(x, y) log_img = np.log10(img) plt.figure() plt.imshow(log_img, origin='lower', interpolation='nearest', vmin=-1, vmax=2) plt.xlabel('x') plt.ylabel('y') cbar = plt.colorbar() cbar.set_label('Log Brightness', rotation=270, labelpad=25) cbar.set_ticks([-1, 0, 1, 2], update_ticks=True) plt.show() References ---------- .. [1] http://ned.ipac.caltech.edu/level5/March05/Graham/Graham2.html r)r(rrNc Cs|jdkr>yddlm} | |_Wntk r<tdYnX|d|d} |d||} } t| t| }}||||||}|| ||||}t|| d|| d}|t | |d|dS) z(Two dimensional Sersic profile function.Nr)rz%Sersic2D model requires scipy > 0.11.g@g?rr') rrrrSrr6r`rarTr7)rr8r@r9rrrrelliprOrZbnr]r^Z cos_thetaZ sin_thetaZx_majZx_minrmr1r1r2r: s   zSersic2D.evaluatecCs&|jjdkrdS|jj|jjdSdS)N)r8r@)rr>r)r-r1r1r2r? s zSersic2D.input_unitscCsP|d|dkrtdtd|dfd|dfd|dfdtjfd|d fgS) Nr8r@z(Units of 'x' and 'y' inputs should matchrrrrOr9rm)r rrnro)r-rArBr1r1r2rC s   z(Sersic2D._parameter_units_for_data_units)rDrErFrGrr9rrrrrrOrrr:rIr?rCr1r1r1r2r M sG         )sb      09%(_V;/ciM1Ql^YbNX?e_X