1# Licensed under a 3-clause BSD style license - see LICENSE.rst 2 3"""Mathematical models.""" 4# pylint: disable=line-too-long, too-many-lines, too-many-arguments, invalid-name 5import numpy as np 6 7from astropy import units as u 8from astropy.units import Quantity, UnitsError 9from .core import (Fittable1DModel, Fittable2DModel) 10 11from .parameters import Parameter, InputParameterError 12from .utils import ellipse_extent 13 14 15__all__ = ['AiryDisk2D', 'Moffat1D', 'Moffat2D', 'Box1D', 'Box2D', 'Const1D', 16 'Const2D', 'Ellipse2D', 'Disk2D', 'Gaussian1D', 'Gaussian2D', 17 'Linear1D', 'Lorentz1D', 'RickerWavelet1D', 'RickerWavelet2D', 18 'RedshiftScaleFactor', 'Multiply', 'Planar2D', 'Scale', 19 'Sersic1D', 'Sersic2D', 'Shift', 20 'Sine1D', 'Cosine1D', 'Tangent1D', 21 'ArcSine1D', 'ArcCosine1D', 'ArcTangent1D', 22 'Trapezoid1D', 'TrapezoidDisk2D', 'Ring2D', 'Voigt1D', 23 'KingProjectedAnalytic1D', 'Exponential1D', 'Logarithmic1D'] 24 25TWOPI = 2 * np.pi 26FLOAT_EPSILON = float(np.finfo(np.float32).tiny) 27 28# Note that we define this here rather than using the value defined in 29# astropy.stats to avoid importing astropy.stats every time astropy.modeling 30# is loaded. 31GAUSSIAN_SIGMA_TO_FWHM = 2.0 * np.sqrt(2.0 * np.log(2.0)) 32 33 34class Gaussian1D(Fittable1DModel): 35 """ 36 One dimensional Gaussian model. 37 38 Parameters 39 ---------- 40 amplitude : float or `~astropy.units.Quantity`. 41 Amplitude (peak value) of the Gaussian - for a normalized profile 42 (integrating to 1), set amplitude = 1 / (stddev * np.sqrt(2 * np.pi)) 43 mean : float or `~astropy.units.Quantity`. 44 Mean of the Gaussian. 45 stddev : float or `~astropy.units.Quantity`. 46 Standard deviation of the Gaussian with FWHM = 2 * stddev * np.sqrt(2 * np.log(2)). 47 48 Notes 49 ----- 50 Either all or none of input ``x``, ``mean`` and ``stddev`` must be provided 51 consistently with compatible units or as unitless numbers. 52 53 Model formula: 54 55 .. math:: f(x) = A e^{- \\frac{\\left(x - x_{0}\\right)^{2}}{2 \\sigma^{2}}} 56 57 Examples 58 -------- 59 >>> from astropy.modeling import models 60 >>> def tie_center(model): 61 ... mean = 50 * model.stddev 62 ... return mean 63 >>> tied_parameters = {'mean': tie_center} 64 65 Specify that 'mean' is a tied parameter in one of two ways: 66 67 >>> g1 = models.Gaussian1D(amplitude=10, mean=5, stddev=.3, 68 ... tied=tied_parameters) 69 70 or 71 72 >>> g1 = models.Gaussian1D(amplitude=10, mean=5, stddev=.3) 73 >>> g1.mean.tied 74 False 75 >>> g1.mean.tied = tie_center 76 >>> g1.mean.tied 77 <function tie_center at 0x...> 78 79 Fixed parameters: 80 81 >>> g1 = models.Gaussian1D(amplitude=10, mean=5, stddev=.3, 82 ... fixed={'stddev': True}) 83 >>> g1.stddev.fixed 84 True 85 86 or 87 88 >>> g1 = models.Gaussian1D(amplitude=10, mean=5, stddev=.3) 89 >>> g1.stddev.fixed 90 False 91 >>> g1.stddev.fixed = True 92 >>> g1.stddev.fixed 93 True 94 95 .. plot:: 96 :include-source: 97 98 import numpy as np 99 import matplotlib.pyplot as plt 100 101 from astropy.modeling.models import Gaussian1D 102 103 plt.figure() 104 s1 = Gaussian1D() 105 r = np.arange(-5, 5, .01) 106 107 for factor in range(1, 4): 108 s1.amplitude = factor 109 plt.plot(r, s1(r), color=str(0.25 * factor), lw=2) 110 111 plt.axis([-5, 5, -1, 4]) 112 plt.show() 113 114 See Also 115 -------- 116 Gaussian2D, Box1D, Moffat1D, Lorentz1D 117 """ 118 119 amplitude = Parameter(default=1, description="Amplitude (peak value) of the Gaussian") 120 mean = Parameter(default=0, description="Position of peak (Gaussian)") 121 122 # Ensure stddev makes sense if its bounds are not explicitly set. 123 # stddev must be non-zero and positive. 124 stddev = Parameter(default=1, bounds=(FLOAT_EPSILON, None), description="Standard deviation of the Gaussian") 125 126 def bounding_box(self, factor=5.5): 127 """ 128 Tuple defining the default ``bounding_box`` limits, 129 ``(x_low, x_high)`` 130 131 Parameters 132 ---------- 133 factor : float 134 The multiple of `stddev` used to define the limits. 135 The default is 5.5, corresponding to a relative error < 1e-7. 136 137 Examples 138 -------- 139 >>> from astropy.modeling.models import Gaussian1D 140 >>> model = Gaussian1D(mean=0, stddev=2) 141 >>> model.bounding_box 142 (-11.0, 11.0) 143 144 This range can be set directly (see: `Model.bounding_box 145 <astropy.modeling.Model.bounding_box>`) or by using a different factor, 146 like: 147 148 >>> model.bounding_box = model.bounding_box(factor=2) 149 >>> model.bounding_box 150 (-4.0, 4.0) 151 """ 152 153 x0 = self.mean 154 dx = factor * self.stddev 155 156 return (x0 - dx, x0 + dx) 157 158 @property 159 def fwhm(self): 160 """Gaussian full width at half maximum.""" 161 return self.stddev * GAUSSIAN_SIGMA_TO_FWHM 162 163 @staticmethod 164 def evaluate(x, amplitude, mean, stddev): 165 """ 166 Gaussian1D model function. 167 """ 168 return amplitude * np.exp(- 0.5 * (x - mean) ** 2 / stddev ** 2) 169 170 @staticmethod 171 def fit_deriv(x, amplitude, mean, stddev): 172 """ 173 Gaussian1D model function derivatives. 174 """ 175 176 d_amplitude = np.exp(-0.5 / stddev ** 2 * (x - mean) ** 2) 177 d_mean = amplitude * d_amplitude * (x - mean) / stddev ** 2 178 d_stddev = amplitude * d_amplitude * (x - mean) ** 2 / stddev ** 3 179 return [d_amplitude, d_mean, d_stddev] 180 181 @property 182 def input_units(self): 183 if self.mean.unit is None: 184 return None 185 return {self.inputs[0]: self.mean.unit} 186 187 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 188 return {'mean': inputs_unit[self.inputs[0]], 189 'stddev': inputs_unit[self.inputs[0]], 190 'amplitude': outputs_unit[self.outputs[0]]} 191 192 193class Gaussian2D(Fittable2DModel): 194 r""" 195 Two dimensional Gaussian model. 196 197 Parameters 198 ---------- 199 amplitude : float or `~astropy.units.Quantity`. 200 Amplitude (peak value) of the Gaussian. 201 x_mean : float or `~astropy.units.Quantity`. 202 Mean of the Gaussian in x. 203 y_mean : float or `~astropy.units.Quantity`. 204 Mean of the Gaussian in y. 205 x_stddev : float or `~astropy.units.Quantity` or None. 206 Standard deviation of the Gaussian in x before rotating by theta. Must 207 be None if a covariance matrix (``cov_matrix``) is provided. If no 208 ``cov_matrix`` is given, ``None`` means the default value (1). 209 y_stddev : float or `~astropy.units.Quantity` or None. 210 Standard deviation of the Gaussian in y before rotating by theta. Must 211 be None if a covariance matrix (``cov_matrix``) is provided. If no 212 ``cov_matrix`` is given, ``None`` means the default value (1). 213 theta : float or `~astropy.units.Quantity`, optional. 214 Rotation angle (value in radians). The rotation angle increases 215 counterclockwise. Must be None if a covariance matrix (``cov_matrix``) 216 is provided. If no ``cov_matrix`` is given, ``None`` means the default 217 value (0). 218 cov_matrix : ndarray, optional 219 A 2x2 covariance matrix. If specified, overrides the ``x_stddev``, 220 ``y_stddev``, and ``theta`` defaults. 221 222 Notes 223 ----- 224 Either all or none of input ``x, y``, ``[x,y]_mean`` and ``[x,y]_stddev`` 225 must be provided consistently with compatible units or as unitless numbers. 226 227 Model formula: 228 229 .. math:: 230 231 f(x, y) = A e^{-a\left(x - x_{0}\right)^{2} -b\left(x - x_{0}\right) 232 \left(y - y_{0}\right) -c\left(y - y_{0}\right)^{2}} 233 234 Using the following definitions: 235 236 .. math:: 237 a = \left(\frac{\cos^{2}{\left (\theta \right )}}{2 \sigma_{x}^{2}} + 238 \frac{\sin^{2}{\left (\theta \right )}}{2 \sigma_{y}^{2}}\right) 239 240 b = \left(\frac{\sin{\left (2 \theta \right )}}{2 \sigma_{x}^{2}} - 241 \frac{\sin{\left (2 \theta \right )}}{2 \sigma_{y}^{2}}\right) 242 243 c = \left(\frac{\sin^{2}{\left (\theta \right )}}{2 \sigma_{x}^{2}} + 244 \frac{\cos^{2}{\left (\theta \right )}}{2 \sigma_{y}^{2}}\right) 245 246 If using a ``cov_matrix``, the model is of the form: 247 .. math:: 248 f(x, y) = A e^{-0.5 \left(\vec{x} - \vec{x}_{0}\right)^{T} \Sigma^{-1} \left(\vec{x} - \vec{x}_{0}\right)} 249 250 where :math:`\vec{x} = [x, y]`, :math:`\vec{x}_{0} = [x_{0}, y_{0}]`, 251 and :math:`\Sigma` is the covariance matrix: 252 253 .. math:: 254 \Sigma = \left(\begin{array}{ccc} 255 \sigma_x^2 & \rho \sigma_x \sigma_y \\ 256 \rho \sigma_x \sigma_y & \sigma_y^2 257 \end{array}\right) 258 259 :math:`\rho` is the correlation between ``x`` and ``y``, which should 260 be between -1 and +1. Positive correlation corresponds to a 261 ``theta`` in the range 0 to 90 degrees. Negative correlation 262 corresponds to a ``theta`` in the range of 0 to -90 degrees. 263 264 See [1]_ for more details about the 2D Gaussian function. 265 266 See Also 267 -------- 268 Gaussian1D, Box2D, Moffat2D 269 270 References 271 ---------- 272 .. [1] https://en.wikipedia.org/wiki/Gaussian_function 273 """ 274 275 amplitude = Parameter(default=1, description="Amplitude of the Gaussian") 276 x_mean = Parameter(default=0, description="Peak position (along x axis) of Gaussian") 277 y_mean = Parameter(default=0, description="Peak position (along y axis) of Gaussian") 278 x_stddev = Parameter(default=1, description="Standard deviation of the Gaussian (along x axis)") 279 y_stddev = Parameter(default=1, description="Standard deviation of the Gaussian (along y axis)") 280 theta = Parameter(default=0.0, description="Rotation angle [in radians] (Optional parameter)") 281 282 def __init__(self, amplitude=amplitude.default, x_mean=x_mean.default, 283 y_mean=y_mean.default, x_stddev=None, y_stddev=None, 284 theta=None, cov_matrix=None, **kwargs): 285 if cov_matrix is None: 286 if x_stddev is None: 287 x_stddev = self.__class__.x_stddev.default 288 if y_stddev is None: 289 y_stddev = self.__class__.y_stddev.default 290 if theta is None: 291 theta = self.__class__.theta.default 292 else: 293 if x_stddev is not None or y_stddev is not None or theta is not None: 294 raise InputParameterError("Cannot specify both cov_matrix and " 295 "x/y_stddev/theta") 296 # Compute principle coordinate system transformation 297 cov_matrix = np.array(cov_matrix) 298 299 if cov_matrix.shape != (2, 2): 300 raise ValueError("Covariance matrix must be 2x2") 301 302 eig_vals, eig_vecs = np.linalg.eig(cov_matrix) 303 x_stddev, y_stddev = np.sqrt(eig_vals) 304 y_vec = eig_vecs[:, 0] 305 theta = np.arctan2(y_vec[1], y_vec[0]) 306 307 # Ensure stddev makes sense if its bounds are not explicitly set. 308 # stddev must be non-zero and positive. 309 # TODO: Investigate why setting this in Parameter above causes 310 # convolution tests to hang. 311 kwargs.setdefault('bounds', {}) 312 kwargs['bounds'].setdefault('x_stddev', (FLOAT_EPSILON, None)) 313 kwargs['bounds'].setdefault('y_stddev', (FLOAT_EPSILON, None)) 314 315 super().__init__( 316 amplitude=amplitude, x_mean=x_mean, y_mean=y_mean, 317 x_stddev=x_stddev, y_stddev=y_stddev, theta=theta, **kwargs) 318 319 @property 320 def x_fwhm(self): 321 """Gaussian full width at half maximum in X.""" 322 return self.x_stddev * GAUSSIAN_SIGMA_TO_FWHM 323 324 @property 325 def y_fwhm(self): 326 """Gaussian full width at half maximum in Y.""" 327 return self.y_stddev * GAUSSIAN_SIGMA_TO_FWHM 328 329 def bounding_box(self, factor=5.5): 330 """ 331 Tuple defining the default ``bounding_box`` limits in each dimension, 332 ``((y_low, y_high), (x_low, x_high))`` 333 334 The default offset from the mean is 5.5-sigma, corresponding 335 to a relative error < 1e-7. The limits are adjusted for rotation. 336 337 Parameters 338 ---------- 339 factor : float, optional 340 The multiple of `x_stddev` and `y_stddev` used to define the limits. 341 The default is 5.5. 342 343 Examples 344 -------- 345 >>> from astropy.modeling.models import Gaussian2D 346 >>> model = Gaussian2D(x_mean=0, y_mean=0, x_stddev=1, y_stddev=2) 347 >>> model.bounding_box 348 ((-11.0, 11.0), (-5.5, 5.5)) 349 350 This range can be set directly (see: `Model.bounding_box 351 <astropy.modeling.Model.bounding_box>`) or by using a different factor 352 like: 353 354 >>> model.bounding_box = model.bounding_box(factor=2) 355 >>> model.bounding_box 356 ((-4.0, 4.0), (-2.0, 2.0)) 357 """ 358 359 a = factor * self.x_stddev 360 b = factor * self.y_stddev 361 theta = self.theta.value 362 dx, dy = ellipse_extent(a, b, theta) 363 364 return ((self.y_mean - dy, self.y_mean + dy), 365 (self.x_mean - dx, self.x_mean + dx)) 366 367 @staticmethod 368 def evaluate(x, y, amplitude, x_mean, y_mean, x_stddev, y_stddev, theta): 369 """Two dimensional Gaussian function""" 370 371 cost2 = np.cos(theta) ** 2 372 sint2 = np.sin(theta) ** 2 373 sin2t = np.sin(2. * theta) 374 xstd2 = x_stddev ** 2 375 ystd2 = y_stddev ** 2 376 xdiff = x - x_mean 377 ydiff = y - y_mean 378 a = 0.5 * ((cost2 / xstd2) + (sint2 / ystd2)) 379 b = 0.5 * ((sin2t / xstd2) - (sin2t / ystd2)) 380 c = 0.5 * ((sint2 / xstd2) + (cost2 / ystd2)) 381 return amplitude * np.exp(-((a * xdiff ** 2) + (b * xdiff * ydiff) + 382 (c * ydiff ** 2))) 383 384 @staticmethod 385 def fit_deriv(x, y, amplitude, x_mean, y_mean, x_stddev, y_stddev, theta): 386 """Two dimensional Gaussian function derivative with respect to parameters""" 387 388 cost = np.cos(theta) 389 sint = np.sin(theta) 390 cost2 = np.cos(theta) ** 2 391 sint2 = np.sin(theta) ** 2 392 cos2t = np.cos(2. * theta) 393 sin2t = np.sin(2. * theta) 394 xstd2 = x_stddev ** 2 395 ystd2 = y_stddev ** 2 396 xstd3 = x_stddev ** 3 397 ystd3 = y_stddev ** 3 398 xdiff = x - x_mean 399 ydiff = y - y_mean 400 xdiff2 = xdiff ** 2 401 ydiff2 = ydiff ** 2 402 a = 0.5 * ((cost2 / xstd2) + (sint2 / ystd2)) 403 b = 0.5 * ((sin2t / xstd2) - (sin2t / ystd2)) 404 c = 0.5 * ((sint2 / xstd2) + (cost2 / ystd2)) 405 g = amplitude * np.exp(-((a * xdiff2) + (b * xdiff * ydiff) + 406 (c * ydiff2))) 407 da_dtheta = (sint * cost * ((1. / ystd2) - (1. / xstd2))) 408 da_dx_stddev = -cost2 / xstd3 409 da_dy_stddev = -sint2 / ystd3 410 db_dtheta = (cos2t / xstd2) - (cos2t / ystd2) 411 db_dx_stddev = -sin2t / xstd3 412 db_dy_stddev = sin2t / ystd3 413 dc_dtheta = -da_dtheta 414 dc_dx_stddev = -sint2 / xstd3 415 dc_dy_stddev = -cost2 / ystd3 416 dg_dA = g / amplitude 417 dg_dx_mean = g * ((2. * a * xdiff) + (b * ydiff)) 418 dg_dy_mean = g * ((b * xdiff) + (2. * c * ydiff)) 419 dg_dx_stddev = g * (-(da_dx_stddev * xdiff2 + 420 db_dx_stddev * xdiff * ydiff + 421 dc_dx_stddev * ydiff2)) 422 dg_dy_stddev = g * (-(da_dy_stddev * xdiff2 + 423 db_dy_stddev * xdiff * ydiff + 424 dc_dy_stddev * ydiff2)) 425 dg_dtheta = g * (-(da_dtheta * xdiff2 + 426 db_dtheta * xdiff * ydiff + 427 dc_dtheta * ydiff2)) 428 return [dg_dA, dg_dx_mean, dg_dy_mean, dg_dx_stddev, dg_dy_stddev, 429 dg_dtheta] 430 431 @property 432 def input_units(self): 433 if self.x_mean.unit is None and self.y_mean.unit is None: 434 return None 435 return {self.inputs[0]: self.x_mean.unit, 436 self.inputs[1]: self.y_mean.unit} 437 438 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 439 # Note that here we need to make sure that x and y are in the same 440 # units otherwise this can lead to issues since rotation is not well 441 # defined. 442 if inputs_unit[self.inputs[0]] != inputs_unit[self.inputs[1]]: 443 raise UnitsError("Units of 'x' and 'y' inputs should match") 444 return {'x_mean': inputs_unit[self.inputs[0]], 445 'y_mean': inputs_unit[self.inputs[0]], 446 'x_stddev': inputs_unit[self.inputs[0]], 447 'y_stddev': inputs_unit[self.inputs[0]], 448 'theta': u.rad, 449 'amplitude': outputs_unit[self.outputs[0]]} 450 451 452class Shift(Fittable1DModel): 453 """ 454 Shift a coordinate. 455 456 Parameters 457 ---------- 458 offset : float 459 Offset to add to a coordinate. 460 """ 461 462 offset = Parameter(default=0, description="Offset to add to a model") 463 linear = True 464 465 _has_inverse_bounding_box = True 466 467 @property 468 def input_units(self): 469 if self.offset.unit is None: 470 return None 471 return {self.inputs[0]: self.offset.unit} 472 473 @property 474 def inverse(self): 475 """One dimensional inverse Shift model function""" 476 477 inv = self.copy() 478 inv.offset *= -1 479 480 try: 481 self.bounding_box 482 except NotImplementedError: 483 pass 484 else: 485 inv.bounding_box = tuple(self.evaluate(x, self.offset) for x in self.bounding_box) 486 487 return inv 488 489 @staticmethod 490 def evaluate(x, offset): 491 """One dimensional Shift model function""" 492 return x + offset 493 494 @staticmethod 495 def sum_of_implicit_terms(x): 496 """Evaluate the implicit term (x) of one dimensional Shift model""" 497 return x 498 499 @staticmethod 500 def fit_deriv(x, *params): 501 """One dimensional Shift model derivative with respect to parameter""" 502 503 d_offset = np.ones_like(x) 504 return [d_offset] 505 506 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 507 return {'offset': outputs_unit[self.outputs[0]]} 508 509 510class Scale(Fittable1DModel): 511 """ 512 Multiply a model by a dimensionless factor. 513 514 Parameters 515 ---------- 516 factor : float 517 Factor by which to scale a coordinate. 518 519 Notes 520 ----- 521 522 If ``factor`` is a `~astropy.units.Quantity` then the units will be 523 stripped before the scaling operation. 524 525 """ 526 527 factor = Parameter(default=1, description="Factor by which to scale a model") 528 linear = True 529 fittable = True 530 531 _input_units_strict = True 532 _input_units_allow_dimensionless = True 533 534 _has_inverse_bounding_box = True 535 536 @property 537 def input_units(self): 538 if self.factor.unit is None: 539 return None 540 return {self.inputs[0]: self.factor.unit} 541 542 @property 543 def inverse(self): 544 """One dimensional inverse Scale model function""" 545 inv = self.copy() 546 inv.factor = 1 / self.factor 547 548 try: 549 self.bounding_box 550 except NotImplementedError: 551 pass 552 else: 553 inv.bounding_box = tuple(self.evaluate(x, self.factor) for x in self.bounding_box.bounding_box()) 554 555 return inv 556 557 @staticmethod 558 def evaluate(x, factor): 559 """One dimensional Scale model function""" 560 if isinstance(factor, u.Quantity): 561 factor = factor.value 562 563 return factor * x 564 565 @staticmethod 566 def fit_deriv(x, *params): 567 """One dimensional Scale model derivative with respect to parameter""" 568 569 d_factor = x 570 return [d_factor] 571 572 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 573 return {'factor': outputs_unit[self.outputs[0]]} 574 575 576class Multiply(Fittable1DModel): 577 """ 578 Multiply a model by a quantity or number. 579 580 Parameters 581 ---------- 582 factor : float 583 Factor by which to multiply a coordinate. 584 """ 585 586 factor = Parameter(default=1, description="Factor by which to multiply a model") 587 linear = True 588 fittable = True 589 590 _has_inverse_bounding_box = True 591 592 @property 593 def inverse(self): 594 """One dimensional inverse multiply model function""" 595 inv = self.copy() 596 inv.factor = 1 / self.factor 597 598 try: 599 self.bounding_box 600 except NotImplementedError: 601 pass 602 else: 603 inv.bounding_box = tuple(self.evaluate(x, self.factor) for x in self.bounding_box.bounding_box()) 604 605 return inv 606 607 @staticmethod 608 def evaluate(x, factor): 609 """One dimensional multiply model function""" 610 return factor * x 611 612 @staticmethod 613 def fit_deriv(x, *params): 614 """One dimensional multiply model derivative with respect to parameter""" 615 616 d_factor = x 617 return [d_factor] 618 619 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 620 return {'factor': outputs_unit[self.outputs[0]]} 621 622 623class RedshiftScaleFactor(Fittable1DModel): 624 """ 625 One dimensional redshift scale factor model. 626 627 Parameters 628 ---------- 629 z : float 630 Redshift value. 631 632 Notes 633 ----- 634 Model formula: 635 636 .. math:: f(x) = x (1 + z) 637 """ 638 639 z = Parameter(description='Redshift', default=0) 640 641 _has_inverse_bounding_box = True 642 643 @staticmethod 644 def evaluate(x, z): 645 """One dimensional RedshiftScaleFactor model function""" 646 647 return (1 + z) * x 648 649 @staticmethod 650 def fit_deriv(x, z): 651 """One dimensional RedshiftScaleFactor model derivative""" 652 653 d_z = x 654 return [d_z] 655 656 @property 657 def inverse(self): 658 """Inverse RedshiftScaleFactor model""" 659 660 inv = self.copy() 661 inv.z = 1.0 / (1.0 + self.z) - 1.0 662 663 try: 664 self.bounding_box 665 except NotImplementedError: 666 pass 667 else: 668 inv.bounding_box = tuple(self.evaluate(x, self.z) for x in self.bounding_box.bounding_box()) 669 670 return inv 671 672 673class Sersic1D(Fittable1DModel): 674 r""" 675 One dimensional Sersic surface brightness profile. 676 677 Parameters 678 ---------- 679 amplitude : float 680 Surface brightness at r_eff. 681 r_eff : float 682 Effective (half-light) radius 683 n : float 684 Sersic Index. 685 686 See Also 687 -------- 688 Gaussian1D, Moffat1D, Lorentz1D 689 690 Notes 691 ----- 692 Model formula: 693 694 .. math:: 695 696 I(r)=I_e\exp\left\{-b_n\left[\left(\frac{r}{r_{e}}\right)^{(1/n)}-1\right]\right\} 697 698 The constant :math:`b_n` is defined such that :math:`r_e` contains half the total 699 luminosity, and can be solved for numerically. 700 701 .. math:: 702 703 \Gamma(2n) = 2\gamma (b_n,2n) 704 705 Examples 706 -------- 707 .. plot:: 708 :include-source: 709 710 import numpy as np 711 from astropy.modeling.models import Sersic1D 712 import matplotlib.pyplot as plt 713 714 plt.figure() 715 plt.subplot(111, xscale='log', yscale='log') 716 s1 = Sersic1D(amplitude=1, r_eff=5) 717 r=np.arange(0, 100, .01) 718 719 for n in range(1, 10): 720 s1.n = n 721 plt.plot(r, s1(r), color=str(float(n) / 15)) 722 723 plt.axis([1e-1, 30, 1e-2, 1e3]) 724 plt.xlabel('log Radius') 725 plt.ylabel('log Surface Brightness') 726 plt.text(.25, 1.5, 'n=1') 727 plt.text(.25, 300, 'n=10') 728 plt.xticks([]) 729 plt.yticks([]) 730 plt.show() 731 732 References 733 ---------- 734 .. [1] http://ned.ipac.caltech.edu/level5/March05/Graham/Graham2.html 735 """ 736 737 amplitude = Parameter(default=1, description="Surface brightness at r_eff") 738 r_eff = Parameter(default=1, description="Effective (half-light) radius") 739 n = Parameter(default=4, description="Sersic Index") 740 _gammaincinv = None 741 742 @classmethod 743 def evaluate(cls, r, amplitude, r_eff, n): 744 """One dimensional Sersic profile function.""" 745 746 if cls._gammaincinv is None: 747 from scipy.special import gammaincinv 748 cls._gammaincinv = gammaincinv 749 750 return (amplitude * np.exp( 751 -cls._gammaincinv(2 * n, 0.5) * ((r / r_eff) ** (1 / n) - 1))) 752 753 @property 754 def input_units(self): 755 if self.r_eff.unit is None: 756 return None 757 return {self.inputs[0]: self.r_eff.unit} 758 759 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 760 return {'r_eff': inputs_unit[self.inputs[0]], 761 'amplitude': outputs_unit[self.outputs[0]]} 762 763 764class _Trigonometric1D(Fittable1DModel): 765 """ 766 Base class for one dimensional trigonometric and inverse trigonometric models 767 768 Parameters 769 ---------- 770 amplitude : float 771 Oscillation amplitude 772 frequency : float 773 Oscillation frequency 774 phase : float 775 Oscillation phase 776 """ 777 778 amplitude = Parameter(default=1, description="Oscillation amplitude") 779 frequency = Parameter(default=1, description="Oscillation frequency") 780 phase = Parameter(default=0, description="Oscillation phase") 781 782 @property 783 def input_units(self): 784 if self.frequency.unit is None: 785 return None 786 return {self.inputs[0]: 1. / self.frequency.unit} 787 788 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 789 return {'frequency': inputs_unit[self.inputs[0]] ** -1, 790 'amplitude': outputs_unit[self.outputs[0]]} 791 792 793class Sine1D(_Trigonometric1D): 794 """ 795 One dimensional Sine model. 796 797 Parameters 798 ---------- 799 amplitude : float 800 Oscillation amplitude 801 frequency : float 802 Oscillation frequency 803 phase : float 804 Oscillation phase 805 806 See Also 807 -------- 808 ArcSine1D, Cosine1D, Tangent1D, Const1D, Linear1D 809 810 811 Notes 812 ----- 813 Model formula: 814 815 .. math:: f(x) = A \\sin(2 \\pi f x + 2 \\pi p) 816 817 Examples 818 -------- 819 .. plot:: 820 :include-source: 821 822 import numpy as np 823 import matplotlib.pyplot as plt 824 825 from astropy.modeling.models import Sine1D 826 827 plt.figure() 828 s1 = Sine1D(amplitude=1, frequency=.25) 829 r=np.arange(0, 10, .01) 830 831 for amplitude in range(1,4): 832 s1.amplitude = amplitude 833 plt.plot(r, s1(r), color=str(0.25 * amplitude), lw=2) 834 835 plt.axis([0, 10, -5, 5]) 836 plt.show() 837 """ 838 839 @staticmethod 840 def evaluate(x, amplitude, frequency, phase): 841 """One dimensional Sine model function""" 842 # Note: If frequency and x are quantities, they should normally have 843 # inverse units, so that argument ends up being dimensionless. However, 844 # np.sin of a dimensionless quantity will crash, so we remove the 845 # quantity-ness from argument in this case (another option would be to 846 # multiply by * u.rad but this would be slower overall). 847 argument = TWOPI * (frequency * x + phase) 848 if isinstance(argument, Quantity): 849 argument = argument.value 850 return amplitude * np.sin(argument) 851 852 @staticmethod 853 def fit_deriv(x, amplitude, frequency, phase): 854 """One dimensional Sine model derivative""" 855 856 d_amplitude = np.sin(TWOPI * frequency * x + TWOPI * phase) 857 d_frequency = (TWOPI * x * amplitude * 858 np.cos(TWOPI * frequency * x + TWOPI * phase)) 859 d_phase = (TWOPI * amplitude * 860 np.cos(TWOPI * frequency * x + TWOPI * phase)) 861 return [d_amplitude, d_frequency, d_phase] 862 863 @property 864 def inverse(self): 865 """One dimensional inverse of Sine""" 866 867 return ArcSine1D(amplitude=self.amplitude, frequency=self.frequency, phase=self.phase) 868 869 870class Cosine1D(_Trigonometric1D): 871 """ 872 One dimensional Cosine model. 873 874 Parameters 875 ---------- 876 amplitude : float 877 Oscillation amplitude 878 frequency : float 879 Oscillation frequency 880 phase : float 881 Oscillation phase 882 883 See Also 884 -------- 885 ArcCosine1D, Sine1D, Tangent1D, Const1D, Linear1D 886 887 888 Notes 889 ----- 890 Model formula: 891 892 .. math:: f(x) = A \\cos(2 \\pi f x + 2 \\pi p) 893 894 Examples 895 -------- 896 .. plot:: 897 :include-source: 898 899 import numpy as np 900 import matplotlib.pyplot as plt 901 902 from astropy.modeling.models import Cosine1D 903 904 plt.figure() 905 s1 = Cosine1D(amplitude=1, frequency=.25) 906 r=np.arange(0, 10, .01) 907 908 for amplitude in range(1,4): 909 s1.amplitude = amplitude 910 plt.plot(r, s1(r), color=str(0.25 * amplitude), lw=2) 911 912 plt.axis([0, 10, -5, 5]) 913 plt.show() 914 """ 915 916 @staticmethod 917 def evaluate(x, amplitude, frequency, phase): 918 """One dimensional Cosine model function""" 919 # Note: If frequency and x are quantities, they should normally have 920 # inverse units, so that argument ends up being dimensionless. However, 921 # np.sin of a dimensionless quantity will crash, so we remove the 922 # quantity-ness from argument in this case (another option would be to 923 # multiply by * u.rad but this would be slower overall). 924 argument = TWOPI * (frequency * x + phase) 925 if isinstance(argument, Quantity): 926 argument = argument.value 927 return amplitude * np.cos(argument) 928 929 @staticmethod 930 def fit_deriv(x, amplitude, frequency, phase): 931 """One dimensional Cosine model derivative""" 932 933 d_amplitude = np.cos(TWOPI * frequency * x + TWOPI * phase) 934 d_frequency = - (TWOPI * x * amplitude * 935 np.sin(TWOPI * frequency * x + TWOPI * phase)) 936 d_phase = - (TWOPI * amplitude * 937 np.sin(TWOPI * frequency * x + TWOPI * phase)) 938 return [d_amplitude, d_frequency, d_phase] 939 940 @property 941 def inverse(self): 942 """One dimensional inverse of Cosine""" 943 944 return ArcCosine1D(amplitude=self.amplitude, frequency=self.frequency, phase=self.phase) 945 946 947class Tangent1D(_Trigonometric1D): 948 """ 949 One dimensional Tangent model. 950 951 Parameters 952 ---------- 953 amplitude : float 954 Oscillation amplitude 955 frequency : float 956 Oscillation frequency 957 phase : float 958 Oscillation phase 959 960 See Also 961 -------- 962 Sine1D, Cosine1D, Const1D, Linear1D 963 964 965 Notes 966 ----- 967 Model formula: 968 969 .. math:: f(x) = A \\tan(2 \\pi f x + 2 \\pi p) 970 971 Note that the tangent function is undefined for inputs of the form 972 pi/2 + n*pi for all integers n. Thus thus the default bounding box 973 has been restricted to: 974 975 .. math:: [(-1/4 - p)/f, (1/4 - p)/f] 976 977 which is the smallest interval for the tangent function to be continuous 978 on. 979 980 Examples 981 -------- 982 .. plot:: 983 :include-source: 984 985 import numpy as np 986 import matplotlib.pyplot as plt 987 988 from astropy.modeling.models import Tangent1D 989 990 plt.figure() 991 s1 = Tangent1D(amplitude=1, frequency=.25) 992 r=np.arange(0, 10, .01) 993 994 for amplitude in range(1,4): 995 s1.amplitude = amplitude 996 plt.plot(r, s1(r), color=str(0.25 * amplitude), lw=2) 997 998 plt.axis([0, 10, -5, 5]) 999 plt.show() 1000 """ 1001 1002 @staticmethod 1003 def evaluate(x, amplitude, frequency, phase): 1004 """One dimensional Tangent model function""" 1005 # Note: If frequency and x are quantities, they should normally have 1006 # inverse units, so that argument ends up being dimensionless. However, 1007 # np.sin of a dimensionless quantity will crash, so we remove the 1008 # quantity-ness from argument in this case (another option would be to 1009 # multiply by * u.rad but this would be slower overall). 1010 argument = TWOPI * (frequency * x + phase) 1011 if isinstance(argument, Quantity): 1012 argument = argument.value 1013 return amplitude * np.tan(argument) 1014 1015 @staticmethod 1016 def fit_deriv(x, amplitude, frequency, phase): 1017 """One dimensional Tangent model derivative""" 1018 1019 sec = 1 / (np.cos(TWOPI * frequency * x + TWOPI * phase))**2 1020 1021 d_amplitude = np.tan(TWOPI * frequency * x + TWOPI * phase) 1022 d_frequency = TWOPI * x * amplitude * sec 1023 d_phase = TWOPI * amplitude * sec 1024 return [d_amplitude, d_frequency, d_phase] 1025 1026 @property 1027 def inverse(self): 1028 """One dimensional inverse of Tangent""" 1029 1030 return ArcTangent1D(amplitude=self.amplitude, frequency=self.frequency, phase=self.phase) 1031 1032 def bounding_box(self): 1033 """ 1034 Tuple defining the default ``bounding_box`` limits, 1035 ``(x_low, x_high)`` 1036 """ 1037 1038 bbox = [(-1/4 - self.phase) / self.frequency, (1/4 - self.phase) / self.frequency] 1039 1040 if self.frequency.unit is not None: 1041 bbox = bbox / self.frequency.unit 1042 1043 return bbox 1044 1045 1046class _InverseTrigonometric1D(_Trigonometric1D): 1047 """ 1048 Base class for one dimensional inverse trigonometric models 1049 """ 1050 1051 @property 1052 def input_units(self): 1053 if self.amplitude.unit is None: 1054 return None 1055 return {self.inputs[0]: self.amplitude.unit} 1056 1057 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 1058 return {'frequency': outputs_unit[self.outputs[0]] ** -1, 1059 'amplitude': inputs_unit[self.inputs[0]]} 1060 1061 1062class ArcSine1D(_InverseTrigonometric1D): 1063 """ 1064 One dimensional ArcSine model returning values between -pi/2 and pi/2 1065 only. 1066 1067 Parameters 1068 ---------- 1069 amplitude : float 1070 Oscillation amplitude for corresponding Sine 1071 frequency : float 1072 Oscillation frequency for corresponding Sine 1073 phase : float 1074 Oscillation phase for corresponding Sine 1075 1076 See Also 1077 -------- 1078 Sine1D, ArcCosine1D, ArcTangent1D 1079 1080 1081 Notes 1082 ----- 1083 Model formula: 1084 1085 .. math:: f(x) = ((arcsin(x / A) / 2pi) - p) / f 1086 1087 The arcsin function being used for this model will only accept inputs 1088 in [-A, A]; otherwise, a runtime warning will be thrown and the result 1089 will be NaN. To avoid this, the bounding_box has been properly set to 1090 accommodate this; therefore, it is recommended that this model always 1091 be evaluated with the ``with_bounding_box=True`` option. 1092 1093 Examples 1094 -------- 1095 .. plot:: 1096 :include-source: 1097 1098 import numpy as np 1099 import matplotlib.pyplot as plt 1100 1101 from astropy.modeling.models import ArcSine1D 1102 1103 plt.figure() 1104 s1 = ArcSine1D(amplitude=1, frequency=.25) 1105 r=np.arange(-1, 1, .01) 1106 1107 for amplitude in range(1,4): 1108 s1.amplitude = amplitude 1109 plt.plot(r, s1(r), color=str(0.25 * amplitude), lw=2) 1110 1111 plt.axis([-1, 1, -np.pi/2, np.pi/2]) 1112 plt.show() 1113 """ 1114 1115 @staticmethod 1116 def evaluate(x, amplitude, frequency, phase): 1117 """One dimensional ArcSine model function""" 1118 # Note: If frequency and x are quantities, they should normally have 1119 # inverse units, so that argument ends up being dimensionless. However, 1120 # np.sin of a dimensionless quantity will crash, so we remove the 1121 # quantity-ness from argument in this case (another option would be to 1122 # multiply by * u.rad but this would be slower overall). 1123 1124 argument = x / amplitude 1125 if isinstance(argument, Quantity): 1126 argument = argument.value 1127 arc_sine = np.arcsin(argument) / TWOPI 1128 1129 return (arc_sine - phase) / frequency 1130 1131 @staticmethod 1132 def fit_deriv(x, amplitude, frequency, phase): 1133 """One dimensional ArcSine model derivative""" 1134 1135 d_amplitude = - x / (TWOPI * frequency * amplitude**2 * np.sqrt(1 - (x / amplitude)**2)) 1136 d_frequency = (phase - (np.arcsin(x / amplitude) / TWOPI)) / frequency**2 1137 d_phase = - 1 / frequency * np.ones(x.shape) 1138 return [d_amplitude, d_frequency, d_phase] 1139 1140 def bounding_box(self): 1141 """ 1142 Tuple defining the default ``bounding_box`` limits, 1143 ``(x_low, x_high)`` 1144 """ 1145 1146 return -1 * self.amplitude, 1 * self.amplitude 1147 1148 @property 1149 def inverse(self): 1150 """One dimensional inverse of ArcSine""" 1151 1152 return Sine1D(amplitude=self.amplitude, frequency=self.frequency, phase=self.phase) 1153 1154 1155class ArcCosine1D(_InverseTrigonometric1D): 1156 """ 1157 One dimensional ArcCosine returning values between 0 and pi only. 1158 1159 1160 Parameters 1161 ---------- 1162 amplitude : float 1163 Oscillation amplitude for corresponding Cosine 1164 frequency : float 1165 Oscillation frequency for corresponding Cosine 1166 phase : float 1167 Oscillation phase for corresponding Cosine 1168 1169 See Also 1170 -------- 1171 Cosine1D, ArcSine1D, ArcTangent1D 1172 1173 1174 Notes 1175 ----- 1176 Model formula: 1177 1178 .. math:: f(x) = ((arccos(x / A) / 2pi) - p) / f 1179 1180 The arccos function being used for this model will only accept inputs 1181 in [-A, A]; otherwise, a runtime warning will be thrown and the result 1182 will be NaN. To avoid this, the bounding_box has been properly set to 1183 accommodate this; therefore, it is recommended that this model always 1184 be evaluated with the ``with_bounding_box=True`` option. 1185 1186 Examples 1187 -------- 1188 .. plot:: 1189 :include-source: 1190 1191 import numpy as np 1192 import matplotlib.pyplot as plt 1193 1194 from astropy.modeling.models import ArcCosine1D 1195 1196 plt.figure() 1197 s1 = ArcCosine1D(amplitude=1, frequency=.25) 1198 r=np.arange(-1, 1, .01) 1199 1200 for amplitude in range(1,4): 1201 s1.amplitude = amplitude 1202 plt.plot(r, s1(r), color=str(0.25 * amplitude), lw=2) 1203 1204 plt.axis([-1, 1, 0, np.pi]) 1205 plt.show() 1206 """ 1207 1208 @staticmethod 1209 def evaluate(x, amplitude, frequency, phase): 1210 """One dimensional ArcCosine model function""" 1211 # Note: If frequency and x are quantities, they should normally have 1212 # inverse units, so that argument ends up being dimensionless. However, 1213 # np.sin of a dimensionless quantity will crash, so we remove the 1214 # quantity-ness from argument in this case (another option would be to 1215 # multiply by * u.rad but this would be slower overall). 1216 1217 argument = x / amplitude 1218 if isinstance(argument, Quantity): 1219 argument = argument.value 1220 arc_cos = np.arccos(argument) / TWOPI 1221 1222 return (arc_cos - phase) / frequency 1223 1224 @staticmethod 1225 def fit_deriv(x, amplitude, frequency, phase): 1226 """One dimensional ArcCosine model derivative""" 1227 1228 d_amplitude = x / (TWOPI * frequency * amplitude**2 * np.sqrt(1 - (x / amplitude)**2)) 1229 d_frequency = (phase - (np.arccos(x / amplitude) / TWOPI)) / frequency**2 1230 d_phase = - 1 / frequency * np.ones(x.shape) 1231 return [d_amplitude, d_frequency, d_phase] 1232 1233 def bounding_box(self): 1234 """ 1235 Tuple defining the default ``bounding_box`` limits, 1236 ``(x_low, x_high)`` 1237 """ 1238 1239 return -1 * self.amplitude, 1 * self.amplitude 1240 1241 @property 1242 def inverse(self): 1243 """One dimensional inverse of ArcCosine""" 1244 1245 return Cosine1D(amplitude=self.amplitude, frequency=self.frequency, phase=self.phase) 1246 1247 1248class ArcTangent1D(_InverseTrigonometric1D): 1249 """ 1250 One dimensional ArcTangent model returning values between -pi/2 and 1251 pi/2 only. 1252 1253 Parameters 1254 ---------- 1255 amplitude : float 1256 Oscillation amplitude for corresponding Tangent 1257 frequency : float 1258 Oscillation frequency for corresponding Tangent 1259 phase : float 1260 Oscillation phase for corresponding Tangent 1261 1262 See Also 1263 -------- 1264 Tangent1D, ArcSine1D, ArcCosine1D 1265 1266 1267 Notes 1268 ----- 1269 Model formula: 1270 1271 .. math:: f(x) = ((arctan(x / A) / 2pi) - p) / f 1272 1273 Examples 1274 -------- 1275 .. plot:: 1276 :include-source: 1277 1278 import numpy as np 1279 import matplotlib.pyplot as plt 1280 1281 from astropy.modeling.models import ArcTangent1D 1282 1283 plt.figure() 1284 s1 = ArcTangent1D(amplitude=1, frequency=.25) 1285 r=np.arange(-10, 10, .01) 1286 1287 for amplitude in range(1,4): 1288 s1.amplitude = amplitude 1289 plt.plot(r, s1(r), color=str(0.25 * amplitude), lw=2) 1290 1291 plt.axis([-10, 10, -np.pi/2, np.pi/2]) 1292 plt.show() 1293 """ 1294 1295 @staticmethod 1296 def evaluate(x, amplitude, frequency, phase): 1297 """One dimensional ArcTangent model function""" 1298 # Note: If frequency and x are quantities, they should normally have 1299 # inverse units, so that argument ends up being dimensionless. However, 1300 # np.sin of a dimensionless quantity will crash, so we remove the 1301 # quantity-ness from argument in this case (another option would be to 1302 # multiply by * u.rad but this would be slower overall). 1303 1304 argument = x / amplitude 1305 if isinstance(argument, Quantity): 1306 argument = argument.value 1307 arc_cos = np.arctan(argument) / TWOPI 1308 1309 return (arc_cos - phase) / frequency 1310 1311 @staticmethod 1312 def fit_deriv(x, amplitude, frequency, phase): 1313 """One dimensional ArcTangent model derivative""" 1314 1315 d_amplitude = - x / (TWOPI * frequency * amplitude**2 * (1 + (x / amplitude)**2)) 1316 d_frequency = (phase - (np.arctan(x / amplitude) / TWOPI)) / frequency**2 1317 d_phase = - 1 / frequency * np.ones(x.shape) 1318 return [d_amplitude, d_frequency, d_phase] 1319 1320 @property 1321 def inverse(self): 1322 """One dimensional inverse of ArcTangent""" 1323 1324 return Tangent1D(amplitude=self.amplitude, frequency=self.frequency, phase=self.phase) 1325 1326 1327class Linear1D(Fittable1DModel): 1328 """ 1329 One dimensional Line model. 1330 1331 Parameters 1332 ---------- 1333 slope : float 1334 Slope of the straight line 1335 1336 intercept : float 1337 Intercept of the straight line 1338 1339 See Also 1340 -------- 1341 Const1D 1342 1343 Notes 1344 ----- 1345 Model formula: 1346 1347 .. math:: f(x) = a x + b 1348 """ 1349 slope = Parameter(default=1, description="Slope of the straight line") 1350 intercept = Parameter(default=0, description="Intercept of the straight line") 1351 linear = True 1352 1353 @staticmethod 1354 def evaluate(x, slope, intercept): 1355 """One dimensional Line model function""" 1356 1357 return slope * x + intercept 1358 1359 @staticmethod 1360 def fit_deriv(x, *params): 1361 """One dimensional Line model derivative with respect to parameters""" 1362 1363 d_slope = x 1364 d_intercept = np.ones_like(x) 1365 return [d_slope, d_intercept] 1366 1367 @property 1368 def inverse(self): 1369 new_slope = self.slope ** -1 1370 new_intercept = -self.intercept / self.slope 1371 return self.__class__(slope=new_slope, intercept=new_intercept) 1372 1373 @property 1374 def input_units(self): 1375 if self.intercept.unit is None and self.slope.unit is None: 1376 return None 1377 return {self.inputs[0]: self.intercept.unit / self.slope.unit} 1378 1379 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 1380 return {'intercept': outputs_unit[self.outputs[0]], 1381 'slope': outputs_unit[self.outputs[0]] / inputs_unit[self.inputs[0]]} 1382 1383 1384class Planar2D(Fittable2DModel): 1385 """ 1386 Two dimensional Plane model. 1387 1388 Parameters 1389 ---------- 1390 slope_x : float 1391 Slope of the plane in X 1392 1393 slope_y : float 1394 Slope of the plane in Y 1395 1396 intercept : float 1397 Z-intercept of the plane 1398 1399 Notes 1400 ----- 1401 Model formula: 1402 1403 .. math:: f(x, y) = a x + b y + c 1404 """ 1405 1406 slope_x = Parameter(default=1, description="Slope of the plane in X") 1407 slope_y = Parameter(default=1, description="Slope of the plane in Y") 1408 intercept = Parameter(default=0, description="Z-intercept of the plane") 1409 linear = True 1410 1411 @staticmethod 1412 def evaluate(x, y, slope_x, slope_y, intercept): 1413 """Two dimensional Plane model function""" 1414 1415 return slope_x * x + slope_y * y + intercept 1416 1417 @staticmethod 1418 def fit_deriv(x, y, *params): 1419 """Two dimensional Plane model derivative with respect to parameters""" 1420 1421 d_slope_x = x 1422 d_slope_y = y 1423 d_intercept = np.ones_like(x) 1424 return [d_slope_x, d_slope_y, d_intercept] 1425 1426 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 1427 return {'intercept': outputs_unit['z'], 1428 'slope_x': outputs_unit['z'] / inputs_unit['x'], 1429 'slope_y': outputs_unit['z'] / inputs_unit['y']} 1430 1431 1432class Lorentz1D(Fittable1DModel): 1433 """ 1434 One dimensional Lorentzian model. 1435 1436 Parameters 1437 ---------- 1438 amplitude : float or `~astropy.units.Quantity`. 1439 Peak value - for a normalized profile (integrating to 1), 1440 set amplitude = 2 / (np.pi * fwhm) 1441 x_0 : float or `~astropy.units.Quantity`. 1442 Position of the peak 1443 fwhm : float or `~astropy.units.Quantity`. 1444 Full width at half maximum (FWHM) 1445 1446 See Also 1447 -------- 1448 Gaussian1D, Box1D, RickerWavelet1D 1449 1450 Notes 1451 ----- 1452 Either all or none of input ``x``, position ``x_0`` and ``fwhm`` must be provided 1453 consistently with compatible units or as unitless numbers. 1454 1455 Model formula: 1456 1457 .. math:: 1458 1459 f(x) = \\frac{A \\gamma^{2}}{\\gamma^{2} + \\left(x - x_{0}\\right)^{2}} 1460 1461 where :math:`\\gamma` is half of given FWHM. 1462 1463 Examples 1464 -------- 1465 .. plot:: 1466 :include-source: 1467 1468 import numpy as np 1469 import matplotlib.pyplot as plt 1470 1471 from astropy.modeling.models import Lorentz1D 1472 1473 plt.figure() 1474 s1 = Lorentz1D() 1475 r = np.arange(-5, 5, .01) 1476 1477 for factor in range(1, 4): 1478 s1.amplitude = factor 1479 plt.plot(r, s1(r), color=str(0.25 * factor), lw=2) 1480 1481 plt.axis([-5, 5, -1, 4]) 1482 plt.show() 1483 """ 1484 1485 amplitude = Parameter(default=1, description="Peak value") 1486 x_0 = Parameter(default=0, description="Position of the peak") 1487 fwhm = Parameter(default=1, description="Full width at half maximum") 1488 1489 @staticmethod 1490 def evaluate(x, amplitude, x_0, fwhm): 1491 """One dimensional Lorentzian model function""" 1492 1493 return (amplitude * ((fwhm / 2.) ** 2) / ((x - x_0) ** 2 + 1494 (fwhm / 2.) ** 2)) 1495 1496 @staticmethod 1497 def fit_deriv(x, amplitude, x_0, fwhm): 1498 """One dimensional Lorentzian model derivative with respect to parameters""" 1499 1500 d_amplitude = fwhm ** 2 / (fwhm ** 2 + (x - x_0) ** 2) 1501 d_x_0 = (amplitude * d_amplitude * (2 * x - 2 * x_0) / 1502 (fwhm ** 2 + (x - x_0) ** 2)) 1503 d_fwhm = 2 * amplitude * d_amplitude / fwhm * (1 - d_amplitude) 1504 return [d_amplitude, d_x_0, d_fwhm] 1505 1506 def bounding_box(self, factor=25): 1507 """Tuple defining the default ``bounding_box`` limits, 1508 ``(x_low, x_high)``. 1509 1510 Parameters 1511 ---------- 1512 factor : float 1513 The multiple of FWHM used to define the limits. 1514 Default is chosen to include most (99%) of the 1515 area under the curve, while still showing the 1516 central feature of interest. 1517 1518 """ 1519 x0 = self.x_0 1520 dx = factor * self.fwhm 1521 1522 return (x0 - dx, x0 + dx) 1523 1524 @property 1525 def input_units(self): 1526 if self.x_0.unit is None: 1527 return None 1528 return {self.inputs[0]: self.x_0.unit} 1529 1530 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 1531 return {'x_0': inputs_unit[self.inputs[0]], 1532 'fwhm': inputs_unit[self.inputs[0]], 1533 'amplitude': outputs_unit[self.outputs[0]]} 1534 1535 1536class Voigt1D(Fittable1DModel): 1537 """ 1538 One dimensional model for the Voigt profile. 1539 1540 Parameters 1541 ---------- 1542 x_0 : float or `~astropy.units.Quantity` 1543 Position of the peak 1544 amplitude_L : float or `~astropy.units.Quantity`. 1545 The Lorentzian amplitude (peak of the associated Lorentz function) 1546 - for a normalized profile (integrating to 1), set 1547 amplitude_L = 2 / (np.pi * fwhm_L) 1548 fwhm_L : float or `~astropy.units.Quantity` 1549 The Lorentzian full width at half maximum 1550 fwhm_G : float or `~astropy.units.Quantity`. 1551 The Gaussian full width at half maximum 1552 method : str, optional 1553 Algorithm for computing the complex error function; one of 1554 'Humlicek2' (default, fast and generally more accurate than ``rtol=3.e-5``) or 1555 'Scipy', alternatively 'wofz' (requires ``scipy``, almost as fast and 1556 reference in accuracy). 1557 1558 See Also 1559 -------- 1560 Gaussian1D, Lorentz1D 1561 1562 Notes 1563 ----- 1564 Either all or none of input ``x``, position ``x_0`` and the ``fwhm_*`` must be provided 1565 consistently with compatible units or as unitless numbers. 1566 Voigt function is calculated as real part of the complex error function computed from either 1567 Humlicek's rational approximations (JQSRT 21:309, 1979; 27:437, 1982) following 1568 Schreier 2018 (MNRAS 479, 3068; and ``hum2zpf16m`` from his cpfX.py module); or 1569 `~scipy.special.wofz` (implementing 'Faddeeva.cc'). 1570 1571 Examples 1572 -------- 1573 .. plot:: 1574 :include-source: 1575 1576 import numpy as np 1577 from astropy.modeling.models import Voigt1D 1578 import matplotlib.pyplot as plt 1579 1580 plt.figure() 1581 x = np.arange(0, 10, 0.01) 1582 v1 = Voigt1D(x_0=5, amplitude_L=10, fwhm_L=0.5, fwhm_G=0.9) 1583 plt.plot(x, v1(x)) 1584 plt.show() 1585 """ 1586 1587 x_0 = Parameter(default=0, 1588 description="Position of the peak") 1589 amplitude_L = Parameter(default=1, # noqa: N815 1590 description="The Lorentzian amplitude") 1591 fwhm_L = Parameter(default=2/np.pi, # noqa: N815 1592 description="The Lorentzian full width at half maximum") 1593 fwhm_G = Parameter(default=np.log(2), # noqa: N815 1594 description="The Gaussian full width at half maximum") 1595 1596 sqrt_pi = np.sqrt(np.pi) 1597 sqrt_ln2 = np.sqrt(np.log(2)) 1598 sqrt_ln2pi = np.sqrt(np.log(2) * np.pi) 1599 _last_z = np.zeros(1, dtype=complex) 1600 _last_w = np.zeros(1, dtype=float) 1601 _faddeeva = None 1602 1603 def __init__(self, x_0=x_0.default, amplitude_L=amplitude_L.default, # noqa: N803 1604 fwhm_L=fwhm_L.default, fwhm_G=fwhm_G.default, method='humlicek2', # noqa: N803 1605 **kwargs): 1606 if str(method).lower() in ('wofz', 'scipy'): 1607 from scipy.special import wofz 1608 self._faddeeva = wofz 1609 elif str(method).lower() == 'humlicek2': 1610 self._faddeeva = self._hum2zpf16c 1611 else: 1612 raise ValueError(f'Not a valid method for Voigt1D Faddeeva function: {method}.') 1613 self.method = self._faddeeva.__name__ 1614 1615 super().__init__(x_0=x_0, amplitude_L=amplitude_L, fwhm_L=fwhm_L, fwhm_G=fwhm_G, **kwargs) 1616 1617 def _wrap_wofz(self, z): 1618 """Call complex error (Faddeeva) function w(z) implemented by algorithm `method`; 1619 cache results for consecutive calls from `evaluate`, `fit_deriv`.""" 1620 1621 if (z.shape == self._last_z.shape and 1622 np.allclose(z, self._last_z, rtol=1.e-14, atol=1.e-15)): 1623 return self._last_w 1624 1625 self._last_w = self._faddeeva(z) 1626 self._last_z = z 1627 return self._last_w 1628 1629 def evaluate(self, x, x_0, amplitude_L, fwhm_L, fwhm_G): # noqa: N803 1630 """One dimensional Voigt function scaled to Lorentz peak amplitude.""" 1631 1632 z = np.atleast_1d(2 * (x - x_0) + 1j * fwhm_L) * self.sqrt_ln2 / fwhm_G 1633 # The normalised Voigt profile is w.real * self.sqrt_ln2 / (self.sqrt_pi * fwhm_G) * 2 ; 1634 # for the legacy definition we multiply with np.pi * fwhm_L / 2 * amplitude_L 1635 return self._wrap_wofz(z).real * self.sqrt_ln2pi / fwhm_G * fwhm_L * amplitude_L 1636 1637 def fit_deriv(self, x, x_0, amplitude_L, fwhm_L, fwhm_G): # noqa: N803 1638 """Derivative of the one dimensional Voigt function with respect to parameters.""" 1639 1640 s = self.sqrt_ln2 / fwhm_G 1641 z = np.atleast_1d(2 * (x - x_0) + 1j * fwhm_L) * s 1642 # V * constant from McLean implementation (== their Voigt function) 1643 w = self._wrap_wofz(z) * s * fwhm_L * amplitude_L * self.sqrt_pi 1644 1645 # Schreier (2018) Eq. 6 == (dvdx + 1j * dvdy) / (sqrt(pi) * fwhm_L * amplitude_L) 1646 dwdz = -2 * z * w + 2j * s * fwhm_L * amplitude_L 1647 1648 return [-dwdz.real * 2 * s, 1649 w.real / amplitude_L, 1650 w.real / fwhm_L - dwdz.imag * s, 1651 (-w.real - s * (2 * (x - x_0) * dwdz.real - fwhm_L * dwdz.imag)) / fwhm_G] 1652 1653 @property 1654 def input_units(self): 1655 if self.x_0.unit is None: 1656 return None 1657 return {self.inputs[0]: self.x_0.unit} 1658 1659 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 1660 return {'x_0': inputs_unit[self.inputs[0]], 1661 'fwhm_L': inputs_unit[self.inputs[0]], 1662 'fwhm_G': inputs_unit[self.inputs[0]], 1663 'amplitude_L': outputs_unit[self.outputs[0]]} 1664 1665 @staticmethod 1666 def _hum2zpf16c(z, s=10.0): 1667 """Complex error function w(z) for z = x + iy combining Humlicek's rational approximations: 1668 1669 |x| + y > 10: Humlicek (JQSRT, 1982) rational approximation for region II; 1670 else: Humlicek (JQSRT, 1979) rational approximation with n=16 and delta=y0=1.35 1671 1672 Version using a mask and np.place; 1673 single complex argument version of Franz Schreier's cpfX.hum2zpf16m. 1674 Originally licensed under a 3-clause BSD style license - see 1675 https://atmos.eoc.dlr.de/tools/lbl4IR/cpfX.py 1676 """ 1677 1678 # Optimized (single fraction) Humlicek region I rational approximation for n=16, delta=1.35 1679 1680 AA = np.array([+46236.3358828121, -147726.58393079657j, # noqa: N806 1681 -206562.80451354137, 281369.1590631087j, 1682 +183092.74968253175, -184787.96830696272j, 1683 -66155.39578477248, 57778.05827983565j, 1684 +11682.770904216826, -9442.402767960672j, 1685 -1052.8438624933142, 814.0996198624186j, 1686 +45.94499030751872, -34.59751573708725j, 1687 -0.7616559377907136, 0.5641895835476449j]) # 1j/sqrt(pi) to the 12. digit 1688 1689 bb = np.array([+7918.06640624997, 0.0, 1690 -126689.0625, 0.0, 1691 +295607.8125, 0.0, 1692 -236486.25, 0.0, 1693 +84459.375, 0.0, 1694 -15015.0, 0.0, 1695 +1365.0, 0.0, 1696 -60.0, 0.0, 1697 +1.0]) 1698 1699 sqrt_piinv = 1.0 / np.sqrt(np.pi) 1700 1701 zz = z * z 1702 w = 1j * (z * (zz * sqrt_piinv - 1.410474)) / (0.75 + zz*(zz - 3.0)) 1703 1704 if np.any(z.imag < s): 1705 mask = abs(z.real) + z.imag < s # returns true for interior points 1706 # returns small complex array covering only the interior region 1707 Z = z[np.where(mask)] + 1.35j 1708 ZZ = Z * Z 1709 numer = (((((((((((((((AA[15]*Z + AA[14])*Z + AA[13])*Z + AA[12])*Z + AA[11])*Z + 1710 AA[10])*Z + AA[9])*Z + AA[8])*Z + AA[7])*Z + AA[6])*Z + 1711 AA[5])*Z + AA[4])*Z+AA[3])*Z + AA[2])*Z + AA[1])*Z + AA[0]) 1712 denom = (((((((ZZ + bb[14])*ZZ + bb[12])*ZZ + bb[10])*ZZ+bb[8])*ZZ + bb[6])*ZZ + 1713 bb[4])*ZZ + bb[2])*ZZ + bb[0] 1714 np.place(w, mask, numer / denom) 1715 1716 return w 1717 1718 1719class Const1D(Fittable1DModel): 1720 """ 1721 One dimensional Constant model. 1722 1723 Parameters 1724 ---------- 1725 amplitude : float 1726 Value of the constant function 1727 1728 See Also 1729 -------- 1730 Const2D 1731 1732 Notes 1733 ----- 1734 Model formula: 1735 1736 .. math:: f(x) = A 1737 1738 Examples 1739 -------- 1740 .. plot:: 1741 :include-source: 1742 1743 import numpy as np 1744 import matplotlib.pyplot as plt 1745 1746 from astropy.modeling.models import Const1D 1747 1748 plt.figure() 1749 s1 = Const1D() 1750 r = np.arange(-5, 5, .01) 1751 1752 for factor in range(1, 4): 1753 s1.amplitude = factor 1754 plt.plot(r, s1(r), color=str(0.25 * factor), lw=2) 1755 1756 plt.axis([-5, 5, -1, 4]) 1757 plt.show() 1758 """ 1759 1760 amplitude = Parameter(default=1, description="Value of the constant function") 1761 linear = True 1762 1763 @staticmethod 1764 def evaluate(x, amplitude): 1765 """One dimensional Constant model function""" 1766 1767 if amplitude.size == 1: 1768 # This is slightly faster than using ones_like and multiplying 1769 x = np.empty_like(x, subok=False) 1770 x.fill(amplitude.item()) 1771 else: 1772 # This case is less likely but could occur if the amplitude 1773 # parameter is given an array-like value 1774 x = amplitude * np.ones_like(x, subok=False) 1775 1776 if isinstance(amplitude, Quantity): 1777 return Quantity(x, unit=amplitude.unit, copy=False) 1778 return x 1779 1780 @staticmethod 1781 def fit_deriv(x, amplitude): 1782 """One dimensional Constant model derivative with respect to parameters""" 1783 1784 d_amplitude = np.ones_like(x) 1785 return [d_amplitude] 1786 1787 @property 1788 def input_units(self): 1789 return None 1790 1791 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 1792 return {'amplitude': outputs_unit[self.outputs[0]]} 1793 1794 1795class Const2D(Fittable2DModel): 1796 """ 1797 Two dimensional Constant model. 1798 1799 Parameters 1800 ---------- 1801 amplitude : float 1802 Value of the constant function 1803 1804 See Also 1805 -------- 1806 Const1D 1807 1808 Notes 1809 ----- 1810 Model formula: 1811 1812 .. math:: f(x, y) = A 1813 """ 1814 1815 amplitude = Parameter(default=1, description="Value of the constant function") 1816 linear = True 1817 1818 @staticmethod 1819 def evaluate(x, y, amplitude): 1820 """Two dimensional Constant model function""" 1821 1822 if amplitude.size == 1: 1823 # This is slightly faster than using ones_like and multiplying 1824 x = np.empty_like(x, subok=False) 1825 x.fill(amplitude.item()) 1826 else: 1827 # This case is less likely but could occur if the amplitude 1828 # parameter is given an array-like value 1829 x = amplitude * np.ones_like(x, subok=False) 1830 1831 if isinstance(amplitude, Quantity): 1832 return Quantity(x, unit=amplitude.unit, copy=False) 1833 return x 1834 1835 @property 1836 def input_units(self): 1837 return None 1838 1839 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 1840 return {'amplitude': outputs_unit[self.outputs[0]]} 1841 1842 1843class Ellipse2D(Fittable2DModel): 1844 """ 1845 A 2D Ellipse model. 1846 1847 Parameters 1848 ---------- 1849 amplitude : float 1850 Value of the ellipse. 1851 1852 x_0 : float 1853 x position of the center of the disk. 1854 1855 y_0 : float 1856 y position of the center of the disk. 1857 1858 a : float 1859 The length of the semimajor axis. 1860 1861 b : float 1862 The length of the semiminor axis. 1863 1864 theta : float 1865 The rotation angle in radians of the semimajor axis. The 1866 rotation angle increases counterclockwise from the positive x 1867 axis. 1868 1869 See Also 1870 -------- 1871 Disk2D, Box2D 1872 1873 Notes 1874 ----- 1875 Model formula: 1876 1877 .. math:: 1878 1879 f(x, y) = \\left \\{ 1880 \\begin{array}{ll} 1881 \\mathrm{amplitude} & : \\left[\\frac{(x - x_0) \\cos 1882 \\theta + (y - y_0) \\sin \\theta}{a}\\right]^2 + 1883 \\left[\\frac{-(x - x_0) \\sin \\theta + (y - y_0) 1884 \\cos \\theta}{b}\\right]^2 \\leq 1 \\\\ 1885 0 & : \\mathrm{otherwise} 1886 \\end{array} 1887 \\right. 1888 1889 Examples 1890 -------- 1891 .. plot:: 1892 :include-source: 1893 1894 import numpy as np 1895 from astropy.modeling.models import Ellipse2D 1896 from astropy.coordinates import Angle 1897 import matplotlib.pyplot as plt 1898 import matplotlib.patches as mpatches 1899 x0, y0 = 25, 25 1900 a, b = 20, 10 1901 theta = Angle(30, 'deg') 1902 e = Ellipse2D(amplitude=100., x_0=x0, y_0=y0, a=a, b=b, 1903 theta=theta.radian) 1904 y, x = np.mgrid[0:50, 0:50] 1905 fig, ax = plt.subplots(1, 1) 1906 ax.imshow(e(x, y), origin='lower', interpolation='none', cmap='Greys_r') 1907 e2 = mpatches.Ellipse((x0, y0), 2*a, 2*b, theta.degree, edgecolor='red', 1908 facecolor='none') 1909 ax.add_patch(e2) 1910 plt.show() 1911 """ 1912 1913 amplitude = Parameter(default=1, description="Value of the ellipse") 1914 x_0 = Parameter(default=0, description="X position of the center of the disk.") 1915 y_0 = Parameter(default=0, description="Y position of the center of the disk.") 1916 a = Parameter(default=1, description="The length of the semimajor axis") 1917 b = Parameter(default=1, description="The length of the semiminor axis") 1918 theta = Parameter(default=0, description="The rotation angle in radians of the semimajor axis (Positive - counterclockwise)") 1919 1920 @staticmethod 1921 def evaluate(x, y, amplitude, x_0, y_0, a, b, theta): 1922 """Two dimensional Ellipse model function.""" 1923 1924 xx = x - x_0 1925 yy = y - y_0 1926 cost = np.cos(theta) 1927 sint = np.sin(theta) 1928 numerator1 = (xx * cost) + (yy * sint) 1929 numerator2 = -(xx * sint) + (yy * cost) 1930 in_ellipse = (((numerator1 / a) ** 2 + (numerator2 / b) ** 2) <= 1.) 1931 result = np.select([in_ellipse], [amplitude]) 1932 1933 if isinstance(amplitude, Quantity): 1934 return Quantity(result, unit=amplitude.unit, copy=False) 1935 return result 1936 1937 @property 1938 def bounding_box(self): 1939 """ 1940 Tuple defining the default ``bounding_box`` limits. 1941 1942 ``((y_low, y_high), (x_low, x_high))`` 1943 """ 1944 1945 a = self.a 1946 b = self.b 1947 theta = self.theta.value 1948 dx, dy = ellipse_extent(a, b, theta) 1949 1950 return ((self.y_0 - dy, self.y_0 + dy), 1951 (self.x_0 - dx, self.x_0 + dx)) 1952 1953 @property 1954 def input_units(self): 1955 if self.x_0.unit is None: 1956 return None 1957 return {self.inputs[0]: self.x_0.unit, 1958 self.inputs[1]: self.y_0.unit} 1959 1960 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 1961 # Note that here we need to make sure that x and y are in the same 1962 # units otherwise this can lead to issues since rotation is not well 1963 # defined. 1964 if inputs_unit[self.inputs[0]] != inputs_unit[self.inputs[1]]: 1965 raise UnitsError("Units of 'x' and 'y' inputs should match") 1966 return {'x_0': inputs_unit[self.inputs[0]], 1967 'y_0': inputs_unit[self.inputs[0]], 1968 'a': inputs_unit[self.inputs[0]], 1969 'b': inputs_unit[self.inputs[0]], 1970 'theta': u.rad, 1971 'amplitude': outputs_unit[self.outputs[0]]} 1972 1973 1974class Disk2D(Fittable2DModel): 1975 """ 1976 Two dimensional radial symmetric Disk model. 1977 1978 Parameters 1979 ---------- 1980 amplitude : float 1981 Value of the disk function 1982 x_0 : float 1983 x position center of the disk 1984 y_0 : float 1985 y position center of the disk 1986 R_0 : float 1987 Radius of the disk 1988 1989 See Also 1990 -------- 1991 Box2D, TrapezoidDisk2D 1992 1993 Notes 1994 ----- 1995 Model formula: 1996 1997 .. math:: 1998 1999 f(r) = \\left \\{ 2000 \\begin{array}{ll} 2001 A & : r \\leq R_0 \\\\ 2002 0 & : r > R_0 2003 \\end{array} 2004 \\right. 2005 """ 2006 2007 amplitude = Parameter(default=1, description="Value of disk function") 2008 x_0 = Parameter(default=0, description="X position of center of the disk") 2009 y_0 = Parameter(default=0, description="Y position of center of the disk") 2010 R_0 = Parameter(default=1, description="Radius of the disk") 2011 2012 @staticmethod 2013 def evaluate(x, y, amplitude, x_0, y_0, R_0): 2014 """Two dimensional Disk model function""" 2015 2016 rr = (x - x_0) ** 2 + (y - y_0) ** 2 2017 result = np.select([rr <= R_0 ** 2], [amplitude]) 2018 2019 if isinstance(amplitude, Quantity): 2020 return Quantity(result, unit=amplitude.unit, copy=False) 2021 return result 2022 2023 @property 2024 def bounding_box(self): 2025 """ 2026 Tuple defining the default ``bounding_box`` limits. 2027 2028 ``((y_low, y_high), (x_low, x_high))`` 2029 """ 2030 2031 return ((self.y_0 - self.R_0, self.y_0 + self.R_0), 2032 (self.x_0 - self.R_0, self.x_0 + self.R_0)) 2033 2034 @property 2035 def input_units(self): 2036 if self.x_0.unit is None and self.y_0.unit is None: 2037 return None 2038 return {self.inputs[0]: self.x_0.unit, 2039 self.inputs[1]: self.y_0.unit} 2040 2041 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 2042 # Note that here we need to make sure that x and y are in the same 2043 # units otherwise this can lead to issues since rotation is not well 2044 # defined. 2045 if inputs_unit[self.inputs[0]] != inputs_unit[self.inputs[1]]: 2046 raise UnitsError("Units of 'x' and 'y' inputs should match") 2047 return {'x_0': inputs_unit[self.inputs[0]], 2048 'y_0': inputs_unit[self.inputs[0]], 2049 'R_0': inputs_unit[self.inputs[0]], 2050 'amplitude': outputs_unit[self.outputs[0]]} 2051 2052 2053class Ring2D(Fittable2DModel): 2054 """ 2055 Two dimensional radial symmetric Ring model. 2056 2057 Parameters 2058 ---------- 2059 amplitude : float 2060 Value of the disk function 2061 x_0 : float 2062 x position center of the disk 2063 y_0 : float 2064 y position center of the disk 2065 r_in : float 2066 Inner radius of the ring 2067 width : float 2068 Width of the ring. 2069 r_out : float 2070 Outer Radius of the ring. Can be specified instead of width. 2071 2072 See Also 2073 -------- 2074 Disk2D, TrapezoidDisk2D 2075 2076 Notes 2077 ----- 2078 Model formula: 2079 2080 .. math:: 2081 2082 f(r) = \\left \\{ 2083 \\begin{array}{ll} 2084 A & : r_{in} \\leq r \\leq r_{out} \\\\ 2085 0 & : \\text{else} 2086 \\end{array} 2087 \\right. 2088 2089 Where :math:`r_{out} = r_{in} + r_{width}`. 2090 """ 2091 2092 amplitude = Parameter(default=1, description="Value of the disk function") 2093 x_0 = Parameter(default=0, description="X position of center of disc") 2094 y_0 = Parameter(default=0, description="Y position of center of disc") 2095 r_in = Parameter(default=1, description="Inner radius of the ring") 2096 width = Parameter(default=1, description="Width of the ring") 2097 2098 def __init__(self, amplitude=amplitude.default, x_0=x_0.default, 2099 y_0=y_0.default, r_in=None, width=None, 2100 r_out=None, **kwargs): 2101 if (r_in is None) and (r_out is None) and (width is None): 2102 r_in = self.r_in.default 2103 width = self.width.default 2104 elif (r_in is not None) and (r_out is None) and (width is None): 2105 width = self.width.default 2106 elif (r_in is None) and (r_out is not None) and (width is None): 2107 r_in = self.r_in.default 2108 width = r_out - r_in 2109 elif (r_in is None) and (r_out is None) and (width is not None): 2110 r_in = self.r_in.default 2111 elif (r_in is not None) and (r_out is not None) and (width is None): 2112 width = r_out - r_in 2113 elif (r_in is None) and (r_out is not None) and (width is not None): 2114 r_in = r_out - width 2115 elif (r_in is not None) and (r_out is not None) and (width is not None): 2116 if np.any(width != (r_out - r_in)): 2117 raise InputParameterError("Width must be r_out - r_in") 2118 2119 if np.any(r_in < 0) or np.any(width < 0): 2120 raise InputParameterError(f"{r_in=} and {width=} must both be >=0") 2121 2122 super().__init__( 2123 amplitude=amplitude, x_0=x_0, y_0=y_0, r_in=r_in, width=width, 2124 **kwargs) 2125 2126 @staticmethod 2127 def evaluate(x, y, amplitude, x_0, y_0, r_in, width): 2128 """Two dimensional Ring model function.""" 2129 2130 rr = (x - x_0) ** 2 + (y - y_0) ** 2 2131 r_range = np.logical_and(rr >= r_in ** 2, rr <= (r_in + width) ** 2) 2132 result = np.select([r_range], [amplitude]) 2133 2134 if isinstance(amplitude, Quantity): 2135 return Quantity(result, unit=amplitude.unit, copy=False) 2136 return result 2137 2138 @property 2139 def bounding_box(self): 2140 """ 2141 Tuple defining the default ``bounding_box``. 2142 2143 ``((y_low, y_high), (x_low, x_high))`` 2144 """ 2145 2146 dr = self.r_in + self.width 2147 2148 return ((self.y_0 - dr, self.y_0 + dr), 2149 (self.x_0 - dr, self.x_0 + dr)) 2150 2151 @property 2152 def input_units(self): 2153 if self.x_0.unit is None: 2154 return None 2155 return {self.inputs[0]: self.x_0.unit, 2156 self.inputs[1]: self.y_0.unit} 2157 2158 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 2159 # Note that here we need to make sure that x and y are in the same 2160 # units otherwise this can lead to issues since rotation is not well 2161 # defined. 2162 if inputs_unit[self.inputs[0]] != inputs_unit[self.inputs[1]]: 2163 raise UnitsError("Units of 'x' and 'y' inputs should match") 2164 return {'x_0': inputs_unit[self.inputs[0]], 2165 'y_0': inputs_unit[self.inputs[0]], 2166 'r_in': inputs_unit[self.inputs[0]], 2167 'width': inputs_unit[self.inputs[0]], 2168 'amplitude': outputs_unit[self.outputs[0]]} 2169 2170 2171class Box1D(Fittable1DModel): 2172 """ 2173 One dimensional Box model. 2174 2175 Parameters 2176 ---------- 2177 amplitude : float 2178 Amplitude A 2179 x_0 : float 2180 Position of the center of the box function 2181 width : float 2182 Width of the box 2183 2184 See Also 2185 -------- 2186 Box2D, TrapezoidDisk2D 2187 2188 Notes 2189 ----- 2190 Model formula: 2191 2192 .. math:: 2193 2194 f(x) = \\left \\{ 2195 \\begin{array}{ll} 2196 A & : x_0 - w/2 \\leq x \\leq x_0 + w/2 \\\\ 2197 0 & : \\text{else} 2198 \\end{array} 2199 \\right. 2200 2201 Examples 2202 -------- 2203 .. plot:: 2204 :include-source: 2205 2206 import numpy as np 2207 import matplotlib.pyplot as plt 2208 2209 from astropy.modeling.models import Box1D 2210 2211 plt.figure() 2212 s1 = Box1D() 2213 r = np.arange(-5, 5, .01) 2214 2215 for factor in range(1, 4): 2216 s1.amplitude = factor 2217 s1.width = factor 2218 plt.plot(r, s1(r), color=str(0.25 * factor), lw=2) 2219 2220 plt.axis([-5, 5, -1, 4]) 2221 plt.show() 2222 """ 2223 2224 amplitude = Parameter(default=1, description="Amplitude A") 2225 x_0 = Parameter(default=0, description="Position of center of box function") 2226 width = Parameter(default=1, description="Width of the box") 2227 2228 @staticmethod 2229 def evaluate(x, amplitude, x_0, width): 2230 """One dimensional Box model function""" 2231 2232 inside = np.logical_and(x >= x_0 - width / 2., x <= x_0 + width / 2.) 2233 return np.select([inside], [amplitude], 0) 2234 2235 @property 2236 def bounding_box(self): 2237 """ 2238 Tuple defining the default ``bounding_box`` limits. 2239 2240 ``(x_low, x_high))`` 2241 """ 2242 2243 dx = self.width / 2 2244 2245 return (self.x_0 - dx, self.x_0 + dx) 2246 2247 @property 2248 def input_units(self): 2249 if self.x_0.unit is None: 2250 return None 2251 return {self.inputs[0]: self.x_0.unit} 2252 2253 @property 2254 def return_units(self): 2255 if self.amplitude.unit is None: 2256 return None 2257 return {self.outputs[0]: self.amplitude.unit} 2258 2259 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 2260 return {'x_0': inputs_unit[self.inputs[0]], 2261 'width': inputs_unit[self.inputs[0]], 2262 'amplitude': outputs_unit[self.outputs[0]]} 2263 2264 2265class Box2D(Fittable2DModel): 2266 """ 2267 Two dimensional Box model. 2268 2269 Parameters 2270 ---------- 2271 amplitude : float 2272 Amplitude 2273 x_0 : float 2274 x position of the center of the box function 2275 x_width : float 2276 Width in x direction of the box 2277 y_0 : float 2278 y position of the center of the box function 2279 y_width : float 2280 Width in y direction of the box 2281 2282 See Also 2283 -------- 2284 Box1D, Gaussian2D, Moffat2D 2285 2286 Notes 2287 ----- 2288 Model formula: 2289 2290 .. math:: 2291 2292 f(x, y) = \\left \\{ 2293 \\begin{array}{ll} 2294 A : & x_0 - w_x/2 \\leq x \\leq x_0 + w_x/2 \\text{ and} \\\\ 2295 & y_0 - w_y/2 \\leq y \\leq y_0 + w_y/2 \\\\ 2296 0 : & \\text{else} 2297 \\end{array} 2298 \\right. 2299 2300 """ 2301 2302 amplitude = Parameter(default=1, description="Amplitude") 2303 x_0 = Parameter(default=0, description="X position of the center of the box function") 2304 y_0 = Parameter(default=0, description="Y position of the center of the box function") 2305 x_width = Parameter(default=1, description="Width in x direction of the box") 2306 y_width = Parameter(default=1, description="Width in y direction of the box") 2307 2308 @staticmethod 2309 def evaluate(x, y, amplitude, x_0, y_0, x_width, y_width): 2310 """Two dimensional Box model function""" 2311 2312 x_range = np.logical_and(x >= x_0 - x_width / 2., 2313 x <= x_0 + x_width / 2.) 2314 y_range = np.logical_and(y >= y_0 - y_width / 2., 2315 y <= y_0 + y_width / 2.) 2316 2317 result = np.select([np.logical_and(x_range, y_range)], [amplitude], 0) 2318 2319 if isinstance(amplitude, Quantity): 2320 return Quantity(result, unit=amplitude.unit, copy=False) 2321 return result 2322 2323 @property 2324 def bounding_box(self): 2325 """ 2326 Tuple defining the default ``bounding_box``. 2327 2328 ``((y_low, y_high), (x_low, x_high))`` 2329 """ 2330 2331 dx = self.x_width / 2 2332 dy = self.y_width / 2 2333 2334 return ((self.y_0 - dy, self.y_0 + dy), 2335 (self.x_0 - dx, self.x_0 + dx)) 2336 2337 @property 2338 def input_units(self): 2339 if self.x_0.unit is None: 2340 return None 2341 return {self.inputs[0]: self.x_0.unit, 2342 self.inputs[1]: self.y_0.unit} 2343 2344 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 2345 return {'x_0': inputs_unit[self.inputs[0]], 2346 'y_0': inputs_unit[self.inputs[1]], 2347 'x_width': inputs_unit[self.inputs[0]], 2348 'y_width': inputs_unit[self.inputs[1]], 2349 'amplitude': outputs_unit[self.outputs[0]]} 2350 2351 2352class Trapezoid1D(Fittable1DModel): 2353 """ 2354 One dimensional Trapezoid model. 2355 2356 Parameters 2357 ---------- 2358 amplitude : float 2359 Amplitude of the trapezoid 2360 x_0 : float 2361 Center position of the trapezoid 2362 width : float 2363 Width of the constant part of the trapezoid. 2364 slope : float 2365 Slope of the tails of the trapezoid 2366 2367 See Also 2368 -------- 2369 Box1D, Gaussian1D, Moffat1D 2370 2371 Examples 2372 -------- 2373 .. plot:: 2374 :include-source: 2375 2376 import numpy as np 2377 import matplotlib.pyplot as plt 2378 2379 from astropy.modeling.models import Trapezoid1D 2380 2381 plt.figure() 2382 s1 = Trapezoid1D() 2383 r = np.arange(-5, 5, .01) 2384 2385 for factor in range(1, 4): 2386 s1.amplitude = factor 2387 s1.width = factor 2388 plt.plot(r, s1(r), color=str(0.25 * factor), lw=2) 2389 2390 plt.axis([-5, 5, -1, 4]) 2391 plt.show() 2392 """ 2393 2394 amplitude = Parameter(default=1, description="Amplitude of the trapezoid") 2395 x_0 = Parameter(default=0, description="Center position of the trapezoid") 2396 width = Parameter(default=1, description="Width of constant part of the trapezoid") 2397 slope = Parameter(default=1, description="Slope of the tails of trapezoid") 2398 2399 @staticmethod 2400 def evaluate(x, amplitude, x_0, width, slope): 2401 """One dimensional Trapezoid model function""" 2402 2403 # Compute the four points where the trapezoid changes slope 2404 # x1 <= x2 <= x3 <= x4 2405 x2 = x_0 - width / 2. 2406 x3 = x_0 + width / 2. 2407 x1 = x2 - amplitude / slope 2408 x4 = x3 + amplitude / slope 2409 2410 # Compute model values in pieces between the change points 2411 range_a = np.logical_and(x >= x1, x < x2) 2412 range_b = np.logical_and(x >= x2, x < x3) 2413 range_c = np.logical_and(x >= x3, x < x4) 2414 val_a = slope * (x - x1) 2415 val_b = amplitude 2416 val_c = slope * (x4 - x) 2417 result = np.select([range_a, range_b, range_c], [val_a, val_b, val_c]) 2418 2419 if isinstance(amplitude, Quantity): 2420 return Quantity(result, unit=amplitude.unit, copy=False) 2421 return result 2422 2423 @property 2424 def bounding_box(self): 2425 """ 2426 Tuple defining the default ``bounding_box`` limits. 2427 2428 ``(x_low, x_high))`` 2429 """ 2430 2431 dx = self.width / 2 + self.amplitude / self.slope 2432 2433 return (self.x_0 - dx, self.x_0 + dx) 2434 2435 @property 2436 def input_units(self): 2437 if self.x_0.unit is None: 2438 return None 2439 return {self.inputs[0]: self.x_0.unit} 2440 2441 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 2442 return {'x_0': inputs_unit[self.inputs[0]], 2443 'width': inputs_unit[self.inputs[0]], 2444 'slope': outputs_unit[self.outputs[0]] / inputs_unit[self.inputs[0]], 2445 'amplitude': outputs_unit[self.outputs[0]]} 2446 2447 2448class TrapezoidDisk2D(Fittable2DModel): 2449 """ 2450 Two dimensional circular Trapezoid model. 2451 2452 Parameters 2453 ---------- 2454 amplitude : float 2455 Amplitude of the trapezoid 2456 x_0 : float 2457 x position of the center of the trapezoid 2458 y_0 : float 2459 y position of the center of the trapezoid 2460 R_0 : float 2461 Radius of the constant part of the trapezoid. 2462 slope : float 2463 Slope of the tails of the trapezoid in x direction. 2464 2465 See Also 2466 -------- 2467 Disk2D, Box2D 2468 """ 2469 2470 amplitude = Parameter(default=1, description="Amplitude of the trapezoid") 2471 x_0 = Parameter(default=0, description="X position of the center of the trapezoid") 2472 y_0 = Parameter(default=0, description="Y position of the center of the trapezoid") 2473 R_0 = Parameter(default=1, description="Radius of constant part of trapezoid") 2474 slope = Parameter(default=1, description="Slope of tails of trapezoid in x direction") 2475 2476 @staticmethod 2477 def evaluate(x, y, amplitude, x_0, y_0, R_0, slope): 2478 """Two dimensional Trapezoid Disk model function""" 2479 2480 r = np.sqrt((x - x_0) ** 2 + (y - y_0) ** 2) 2481 range_1 = r <= R_0 2482 range_2 = np.logical_and(r > R_0, r <= R_0 + amplitude / slope) 2483 val_1 = amplitude 2484 val_2 = amplitude + slope * (R_0 - r) 2485 result = np.select([range_1, range_2], [val_1, val_2]) 2486 2487 if isinstance(amplitude, Quantity): 2488 return Quantity(result, unit=amplitude.unit, copy=False) 2489 return result 2490 2491 @property 2492 def bounding_box(self): 2493 """ 2494 Tuple defining the default ``bounding_box``. 2495 2496 ``((y_low, y_high), (x_low, x_high))`` 2497 """ 2498 2499 dr = self.R_0 + self.amplitude / self.slope 2500 2501 return ((self.y_0 - dr, self.y_0 + dr), 2502 (self.x_0 - dr, self.x_0 + dr)) 2503 2504 @property 2505 def input_units(self): 2506 if self.x_0.unit is None and self.y_0.unit is None: 2507 return None 2508 return {self.inputs[0]: self.x_0.unit, 2509 self.inputs[1]: self.y_0.unit} 2510 2511 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 2512 # Note that here we need to make sure that x and y are in the same 2513 # units otherwise this can lead to issues since rotation is not well 2514 # defined. 2515 if inputs_unit['x'] != inputs_unit['y']: 2516 raise UnitsError("Units of 'x' and 'y' inputs should match") 2517 return {'x_0': inputs_unit[self.inputs[0]], 2518 'y_0': inputs_unit[self.inputs[0]], 2519 'R_0': inputs_unit[self.inputs[0]], 2520 'slope': outputs_unit[self.outputs[0]] / inputs_unit[self.inputs[0]], 2521 'amplitude': outputs_unit[self.outputs[0]]} 2522 2523 2524class RickerWavelet1D(Fittable1DModel): 2525 """ 2526 One dimensional Ricker Wavelet model (sometimes known as a "Mexican Hat" 2527 model). 2528 2529 .. note:: 2530 2531 See https://github.com/astropy/astropy/pull/9445 for discussions 2532 related to renaming of this model. 2533 2534 Parameters 2535 ---------- 2536 amplitude : float 2537 Amplitude 2538 x_0 : float 2539 Position of the peak 2540 sigma : float 2541 Width of the Ricker wavelet 2542 2543 See Also 2544 -------- 2545 RickerWavelet2D, Box1D, Gaussian1D, Trapezoid1D 2546 2547 Notes 2548 ----- 2549 Model formula: 2550 2551 .. math:: 2552 2553 f(x) = {A \\left(1 - \\frac{\\left(x - x_{0}\\right)^{2}}{\\sigma^{2}}\\right) 2554 e^{- \\frac{\\left(x - x_{0}\\right)^{2}}{2 \\sigma^{2}}}} 2555 2556 Examples 2557 -------- 2558 .. plot:: 2559 :include-source: 2560 2561 import numpy as np 2562 import matplotlib.pyplot as plt 2563 2564 from astropy.modeling.models import RickerWavelet1D 2565 2566 plt.figure() 2567 s1 = RickerWavelet1D() 2568 r = np.arange(-5, 5, .01) 2569 2570 for factor in range(1, 4): 2571 s1.amplitude = factor 2572 s1.width = factor 2573 plt.plot(r, s1(r), color=str(0.25 * factor), lw=2) 2574 2575 plt.axis([-5, 5, -2, 4]) 2576 plt.show() 2577 """ 2578 2579 amplitude = Parameter(default=1, description="Amplitude (peak) value") 2580 x_0 = Parameter(default=0, description="Position of the peak") 2581 sigma = Parameter(default=1, description="Width of the Ricker wavelet") 2582 2583 @staticmethod 2584 def evaluate(x, amplitude, x_0, sigma): 2585 """One dimensional Ricker Wavelet model function""" 2586 2587 xx_ww = (x - x_0) ** 2 / (2 * sigma ** 2) 2588 return amplitude * (1 - 2 * xx_ww) * np.exp(-xx_ww) 2589 2590 def bounding_box(self, factor=10.0): 2591 """Tuple defining the default ``bounding_box`` limits, 2592 ``(x_low, x_high)``. 2593 2594 Parameters 2595 ---------- 2596 factor : float 2597 The multiple of sigma used to define the limits. 2598 2599 """ 2600 x0 = self.x_0 2601 dx = factor * self.sigma 2602 2603 return (x0 - dx, x0 + dx) 2604 2605 @property 2606 def input_units(self): 2607 if self.x_0.unit is None: 2608 return None 2609 return {self.inputs[0]: self.x_0.unit} 2610 2611 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 2612 return {'x_0': inputs_unit[self.inputs[0]], 2613 'sigma': inputs_unit[self.inputs[0]], 2614 'amplitude': outputs_unit[self.outputs[0]]} 2615 2616 2617class RickerWavelet2D(Fittable2DModel): 2618 """ 2619 Two dimensional Ricker Wavelet model (sometimes known as a "Mexican Hat" 2620 model). 2621 2622 .. note:: 2623 2624 See https://github.com/astropy/astropy/pull/9445 for discussions 2625 related to renaming of this model. 2626 2627 Parameters 2628 ---------- 2629 amplitude : float 2630 Amplitude 2631 x_0 : float 2632 x position of the peak 2633 y_0 : float 2634 y position of the peak 2635 sigma : float 2636 Width of the Ricker wavelet 2637 2638 See Also 2639 -------- 2640 RickerWavelet1D, Gaussian2D 2641 2642 Notes 2643 ----- 2644 Model formula: 2645 2646 .. math:: 2647 2648 f(x, y) = A \\left(1 - \\frac{\\left(x - x_{0}\\right)^{2} 2649 + \\left(y - y_{0}\\right)^{2}}{\\sigma^{2}}\\right) 2650 e^{\\frac{- \\left(x - x_{0}\\right)^{2} 2651 - \\left(y - y_{0}\\right)^{2}}{2 \\sigma^{2}}} 2652 """ 2653 2654 amplitude = Parameter(default=1, description="Amplitude (peak) value") 2655 x_0 = Parameter(default=0, description="X position of the peak") 2656 y_0 = Parameter(default=0, description="Y position of the peak") 2657 sigma = Parameter(default=1, description="Width of the Ricker wavelet") 2658 2659 @staticmethod 2660 def evaluate(x, y, amplitude, x_0, y_0, sigma): 2661 """Two dimensional Ricker Wavelet model function""" 2662 2663 rr_ww = ((x - x_0) ** 2 + (y - y_0) ** 2) / (2 * sigma ** 2) 2664 return amplitude * (1 - rr_ww) * np.exp(- rr_ww) 2665 2666 @property 2667 def input_units(self): 2668 if self.x_0.unit is None: 2669 return None 2670 return {self.inputs[0]: self.x_0.unit, 2671 self.inputs[1]: self.y_0.unit} 2672 2673 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 2674 # Note that here we need to make sure that x and y are in the same 2675 # units otherwise this can lead to issues since rotation is not well 2676 # defined. 2677 if inputs_unit[self.inputs[0]] != inputs_unit[self.inputs[1]]: 2678 raise UnitsError("Units of 'x' and 'y' inputs should match") 2679 return {'x_0': inputs_unit[self.inputs[0]], 2680 'y_0': inputs_unit[self.inputs[0]], 2681 'sigma': inputs_unit[self.inputs[0]], 2682 'amplitude': outputs_unit[self.outputs[0]]} 2683 2684 2685class AiryDisk2D(Fittable2DModel): 2686 """ 2687 Two dimensional Airy disk model. 2688 2689 Parameters 2690 ---------- 2691 amplitude : float 2692 Amplitude of the Airy function. 2693 x_0 : float 2694 x position of the maximum of the Airy function. 2695 y_0 : float 2696 y position of the maximum of the Airy function. 2697 radius : float 2698 The radius of the Airy disk (radius of the first zero). 2699 2700 See Also 2701 -------- 2702 Box2D, TrapezoidDisk2D, Gaussian2D 2703 2704 Notes 2705 ----- 2706 Model formula: 2707 2708 .. math:: f(r) = A \\left[\\frac{2 J_1(\\frac{\\pi r}{R/R_z})}{\\frac{\\pi r}{R/R_z}}\\right]^2 2709 2710 Where :math:`J_1` is the first order Bessel function of the first 2711 kind, :math:`r` is radial distance from the maximum of the Airy 2712 function (:math:`r = \\sqrt{(x - x_0)^2 + (y - y_0)^2}`), :math:`R` 2713 is the input ``radius`` parameter, and :math:`R_z = 2714 1.2196698912665045`). 2715 2716 For an optical system, the radius of the first zero represents the 2717 limiting angular resolution and is approximately 1.22 * lambda / D, 2718 where lambda is the wavelength of the light and D is the diameter of 2719 the aperture. 2720 2721 See [1]_ for more details about the Airy disk. 2722 2723 References 2724 ---------- 2725 .. [1] https://en.wikipedia.org/wiki/Airy_disk 2726 """ 2727 2728 amplitude = Parameter(default=1, description="Amplitude (peak value) of the Airy function") 2729 x_0 = Parameter(default=0, description="X position of the peak") 2730 y_0 = Parameter(default=0, description="Y position of the peak") 2731 radius = Parameter(default=1, 2732 description="The radius of the Airy disk (radius of first zero crossing)") 2733 _rz = None 2734 _j1 = None 2735 2736 @classmethod 2737 def evaluate(cls, x, y, amplitude, x_0, y_0, radius): 2738 """Two dimensional Airy model function""" 2739 2740 if cls._rz is None: 2741 from scipy.special import j1, jn_zeros 2742 cls._rz = jn_zeros(1, 1)[0] / np.pi 2743 cls._j1 = j1 2744 2745 r = np.sqrt((x - x_0) ** 2 + (y - y_0) ** 2) / (radius / cls._rz) 2746 2747 if isinstance(r, Quantity): 2748 # scipy function cannot handle Quantity, so turn into array. 2749 r = r.to_value(u.dimensionless_unscaled) 2750 2751 # Since r can be zero, we have to take care to treat that case 2752 # separately so as not to raise a numpy warning 2753 z = np.ones(r.shape) 2754 rt = np.pi * r[r > 0] 2755 z[r > 0] = (2.0 * cls._j1(rt) / rt) ** 2 2756 2757 if isinstance(amplitude, Quantity): 2758 # make z quantity too, otherwise in-place multiplication fails. 2759 z = Quantity(z, u.dimensionless_unscaled, copy=False) 2760 2761 z *= amplitude 2762 return z 2763 2764 @property 2765 def input_units(self): 2766 if self.x_0.unit is None: 2767 return None 2768 return {self.inputs[0]: self.x_0.unit, 2769 self.inputs[1]: self.y_0.unit} 2770 2771 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 2772 # Note that here we need to make sure that x and y are in the same 2773 # units otherwise this can lead to issues since rotation is not well 2774 # defined. 2775 if inputs_unit[self.inputs[0]] != inputs_unit[self.inputs[1]]: 2776 raise UnitsError("Units of 'x' and 'y' inputs should match") 2777 return {'x_0': inputs_unit[self.inputs[0]], 2778 'y_0': inputs_unit[self.inputs[0]], 2779 'radius': inputs_unit[self.inputs[0]], 2780 'amplitude': outputs_unit[self.outputs[0]]} 2781 2782 2783class Moffat1D(Fittable1DModel): 2784 """ 2785 One dimensional Moffat model. 2786 2787 Parameters 2788 ---------- 2789 amplitude : float 2790 Amplitude of the model. 2791 x_0 : float 2792 x position of the maximum of the Moffat model. 2793 gamma : float 2794 Core width of the Moffat model. 2795 alpha : float 2796 Power index of the Moffat model. 2797 2798 See Also 2799 -------- 2800 Gaussian1D, Box1D 2801 2802 Notes 2803 ----- 2804 Model formula: 2805 2806 .. math:: 2807 2808 f(x) = A \\left(1 + \\frac{\\left(x - x_{0}\\right)^{2}}{\\gamma^{2}}\\right)^{- \\alpha} 2809 2810 Examples 2811 -------- 2812 .. plot:: 2813 :include-source: 2814 2815 import numpy as np 2816 import matplotlib.pyplot as plt 2817 2818 from astropy.modeling.models import Moffat1D 2819 2820 plt.figure() 2821 s1 = Moffat1D() 2822 r = np.arange(-5, 5, .01) 2823 2824 for factor in range(1, 4): 2825 s1.amplitude = factor 2826 s1.width = factor 2827 plt.plot(r, s1(r), color=str(0.25 * factor), lw=2) 2828 2829 plt.axis([-5, 5, -1, 4]) 2830 plt.show() 2831 """ 2832 2833 amplitude = Parameter(default=1, description="Amplitude of the model") 2834 x_0 = Parameter(default=0, description="X position of maximum of Moffat model") 2835 gamma = Parameter(default=1, description="Core width of Moffat model") 2836 alpha = Parameter(default=1, description="Power index of the Moffat model") 2837 2838 @property 2839 def fwhm(self): 2840 """ 2841 Moffat full width at half maximum. 2842 Derivation of the formula is available in 2843 `this notebook by Yoonsoo Bach <https://nbviewer.jupyter.org/github/ysbach/AO_2017/blob/master/04_Ground_Based_Concept.ipynb#1.2.-Moffat>`_. 2844 """ 2845 return 2.0 * np.abs(self.gamma) * np.sqrt(2.0 ** (1.0 / self.alpha) - 1.0) 2846 2847 @staticmethod 2848 def evaluate(x, amplitude, x_0, gamma, alpha): 2849 """One dimensional Moffat model function""" 2850 2851 return amplitude * (1 + ((x - x_0) / gamma) ** 2) ** (-alpha) 2852 2853 @staticmethod 2854 def fit_deriv(x, amplitude, x_0, gamma, alpha): 2855 """One dimensional Moffat model derivative with respect to parameters""" 2856 2857 fac = (1 + (x - x_0) ** 2 / gamma ** 2) 2858 d_A = fac ** (-alpha) 2859 d_x_0 = (2 * amplitude * alpha * (x - x_0) * d_A / (fac * gamma ** 2)) 2860 d_gamma = (2 * amplitude * alpha * (x - x_0) ** 2 * d_A / 2861 (fac * gamma ** 3)) 2862 d_alpha = -amplitude * d_A * np.log(fac) 2863 return [d_A, d_x_0, d_gamma, d_alpha] 2864 2865 @property 2866 def input_units(self): 2867 if self.x_0.unit is None: 2868 return None 2869 return {self.inputs[0]: self.x_0.unit} 2870 2871 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 2872 return {'x_0': inputs_unit[self.inputs[0]], 2873 'gamma': inputs_unit[self.inputs[0]], 2874 'amplitude': outputs_unit[self.outputs[0]]} 2875 2876 2877class Moffat2D(Fittable2DModel): 2878 """ 2879 Two dimensional Moffat model. 2880 2881 Parameters 2882 ---------- 2883 amplitude : float 2884 Amplitude of the model. 2885 x_0 : float 2886 x position of the maximum of the Moffat model. 2887 y_0 : float 2888 y position of the maximum of the Moffat model. 2889 gamma : float 2890 Core width of the Moffat model. 2891 alpha : float 2892 Power index of the Moffat model. 2893 2894 See Also 2895 -------- 2896 Gaussian2D, Box2D 2897 2898 Notes 2899 ----- 2900 Model formula: 2901 2902 .. math:: 2903 2904 f(x, y) = A \\left(1 + \\frac{\\left(x - x_{0}\\right)^{2} + 2905 \\left(y - y_{0}\\right)^{2}}{\\gamma^{2}}\\right)^{- \\alpha} 2906 """ 2907 2908 amplitude = Parameter(default=1, description="Amplitude (peak value) of the model") 2909 x_0 = Parameter(default=0, description="X position of the maximum of the Moffat model") 2910 y_0 = Parameter(default=0, description="Y position of the maximum of the Moffat model") 2911 gamma = Parameter(default=1, description="Core width of the Moffat model") 2912 alpha = Parameter(default=1, description="Power index of the Moffat model") 2913 2914 @property 2915 def fwhm(self): 2916 """ 2917 Moffat full width at half maximum. 2918 Derivation of the formula is available in 2919 `this notebook by Yoonsoo Bach <https://nbviewer.jupyter.org/github/ysbach/AO_2017/blob/master/04_Ground_Based_Concept.ipynb#1.2.-Moffat>`_. 2920 """ 2921 return 2.0 * np.abs(self.gamma) * np.sqrt(2.0 ** (1.0 / self.alpha) - 1.0) 2922 2923 @staticmethod 2924 def evaluate(x, y, amplitude, x_0, y_0, gamma, alpha): 2925 """Two dimensional Moffat model function""" 2926 2927 rr_gg = ((x - x_0) ** 2 + (y - y_0) ** 2) / gamma ** 2 2928 return amplitude * (1 + rr_gg) ** (-alpha) 2929 2930 @staticmethod 2931 def fit_deriv(x, y, amplitude, x_0, y_0, gamma, alpha): 2932 """Two dimensional Moffat model derivative with respect to parameters""" 2933 2934 rr_gg = ((x - x_0) ** 2 + (y - y_0) ** 2) / gamma ** 2 2935 d_A = (1 + rr_gg) ** (-alpha) 2936 d_x_0 = (2 * amplitude * alpha * d_A * (x - x_0) / 2937 (gamma ** 2 * (1 + rr_gg))) 2938 d_y_0 = (2 * amplitude * alpha * d_A * (y - y_0) / 2939 (gamma ** 2 * (1 + rr_gg))) 2940 d_alpha = -amplitude * d_A * np.log(1 + rr_gg) 2941 d_gamma = (2 * amplitude * alpha * d_A * rr_gg / 2942 (gamma * (1 + rr_gg))) 2943 return [d_A, d_x_0, d_y_0, d_gamma, d_alpha] 2944 2945 @property 2946 def input_units(self): 2947 if self.x_0.unit is None: 2948 return None 2949 else: 2950 return {self.inputs[0]: self.x_0.unit, 2951 self.inputs[1]: self.y_0.unit} 2952 2953 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 2954 # Note that here we need to make sure that x and y are in the same 2955 # units otherwise this can lead to issues since rotation is not well 2956 # defined. 2957 if inputs_unit[self.inputs[0]] != inputs_unit[self.inputs[1]]: 2958 raise UnitsError("Units of 'x' and 'y' inputs should match") 2959 return {'x_0': inputs_unit[self.inputs[0]], 2960 'y_0': inputs_unit[self.inputs[0]], 2961 'gamma': inputs_unit[self.inputs[0]], 2962 'amplitude': outputs_unit[self.outputs[0]]} 2963 2964 2965class Sersic2D(Fittable2DModel): 2966 r""" 2967 Two dimensional Sersic surface brightness profile. 2968 2969 Parameters 2970 ---------- 2971 amplitude : float 2972 Surface brightness at r_eff. 2973 r_eff : float 2974 Effective (half-light) radius 2975 n : float 2976 Sersic Index. 2977 x_0 : float, optional 2978 x position of the center. 2979 y_0 : float, optional 2980 y position of the center. 2981 ellip : float, optional 2982 Ellipticity. 2983 theta : float, optional 2984 Rotation angle in radians, counterclockwise from 2985 the positive x-axis. 2986 2987 See Also 2988 -------- 2989 Gaussian2D, Moffat2D 2990 2991 Notes 2992 ----- 2993 Model formula: 2994 2995 .. math:: 2996 2997 I(x,y) = I(r) = I_e\exp\left\{-b_n\left[\left(\frac{r}{r_{e}}\right)^{(1/n)}-1\right]\right\} 2998 2999 The constant :math:`b_n` is defined such that :math:`r_e` contains half the total 3000 luminosity, and can be solved for numerically. 3001 3002 .. math:: 3003 3004 \Gamma(2n) = 2\gamma (2n,b_n) 3005 3006 Examples 3007 -------- 3008 .. plot:: 3009 :include-source: 3010 3011 import numpy as np 3012 from astropy.modeling.models import Sersic2D 3013 import matplotlib.pyplot as plt 3014 3015 x,y = np.meshgrid(np.arange(100), np.arange(100)) 3016 3017 mod = Sersic2D(amplitude = 1, r_eff = 25, n=4, x_0=50, y_0=50, 3018 ellip=.5, theta=-1) 3019 img = mod(x, y) 3020 log_img = np.log10(img) 3021 3022 3023 plt.figure() 3024 plt.imshow(log_img, origin='lower', interpolation='nearest', 3025 vmin=-1, vmax=2) 3026 plt.xlabel('x') 3027 plt.ylabel('y') 3028 cbar = plt.colorbar() 3029 cbar.set_label('Log Brightness', rotation=270, labelpad=25) 3030 cbar.set_ticks([-1, 0, 1, 2], update_ticks=True) 3031 plt.show() 3032 3033 References 3034 ---------- 3035 .. [1] http://ned.ipac.caltech.edu/level5/March05/Graham/Graham2.html 3036 """ 3037 3038 amplitude = Parameter(default=1, description="Surface brightness at r_eff") 3039 r_eff = Parameter(default=1, description="Effective (half-light) radius") 3040 n = Parameter(default=4, description="Sersic Index") 3041 x_0 = Parameter(default=0, description="X position of the center") 3042 y_0 = Parameter(default=0, description="Y position of the center") 3043 ellip = Parameter(default=0, description="Ellipticity") 3044 theta = Parameter(default=0, description="Rotation angle in radians (counterclockwise-positive)") 3045 _gammaincinv = None 3046 3047 @classmethod 3048 def evaluate(cls, x, y, amplitude, r_eff, n, x_0, y_0, ellip, theta): 3049 """Two dimensional Sersic profile function.""" 3050 3051 if cls._gammaincinv is None: 3052 from scipy.special import gammaincinv 3053 cls._gammaincinv = gammaincinv 3054 3055 bn = cls._gammaincinv(2. * n, 0.5) 3056 a, b = r_eff, (1 - ellip) * r_eff 3057 cos_theta, sin_theta = np.cos(theta), np.sin(theta) 3058 x_maj = (x - x_0) * cos_theta + (y - y_0) * sin_theta 3059 x_min = -(x - x_0) * sin_theta + (y - y_0) * cos_theta 3060 z = np.sqrt((x_maj / a) ** 2 + (x_min / b) ** 2) 3061 3062 return amplitude * np.exp(-bn * (z ** (1 / n) - 1)) 3063 3064 @property 3065 def input_units(self): 3066 if self.x_0.unit is None: 3067 return None 3068 return {self.inputs[0]: self.x_0.unit, 3069 self.inputs[1]: self.y_0.unit} 3070 3071 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 3072 # Note that here we need to make sure that x and y are in the same 3073 # units otherwise this can lead to issues since rotation is not well 3074 # defined. 3075 if inputs_unit[self.inputs[0]] != inputs_unit[self.inputs[1]]: 3076 raise UnitsError("Units of 'x' and 'y' inputs should match") 3077 return {'x_0': inputs_unit[self.inputs[0]], 3078 'y_0': inputs_unit[self.inputs[0]], 3079 'r_eff': inputs_unit[self.inputs[0]], 3080 'theta': u.rad, 3081 'amplitude': outputs_unit[self.outputs[0]]} 3082 3083 3084class KingProjectedAnalytic1D(Fittable1DModel): 3085 """ 3086 Projected (surface density) analytic King Model. 3087 3088 3089 Parameters 3090 ---------- 3091 amplitude : float 3092 Amplitude or scaling factor. 3093 r_core : float 3094 Core radius (f(r_c) ~ 0.5 f_0) 3095 r_tide : float 3096 Tidal radius. 3097 3098 3099 Notes 3100 ----- 3101 3102 This model approximates a King model with an analytic function. The derivation of this 3103 equation can be found in King '62 (equation 14). This is just an approximation of the 3104 full model and the parameters derived from this model should be taken with caution. 3105 It usually works for models with a concentration (c = log10(r_t/r_c) parameter < 2. 3106 3107 Model formula: 3108 3109 .. math:: 3110 3111 f(x) = A r_c^2 \\left(\\frac{1}{\\sqrt{(x^2 + r_c^2)}} - 3112 \\frac{1}{\\sqrt{(r_t^2 + r_c^2)}}\\right)^2 3113 3114 Examples 3115 -------- 3116 .. plot:: 3117 :include-source: 3118 3119 import numpy as np 3120 from astropy.modeling.models import KingProjectedAnalytic1D 3121 import matplotlib.pyplot as plt 3122 3123 plt.figure() 3124 rt_list = [1, 2, 5, 10, 20] 3125 for rt in rt_list: 3126 r = np.linspace(0.1, rt, 100) 3127 3128 mod = KingProjectedAnalytic1D(amplitude = 1, r_core = 1., r_tide = rt) 3129 sig = mod(r) 3130 3131 3132 plt.loglog(r, sig/sig[0], label='c ~ {:0.2f}'.format(mod.concentration)) 3133 3134 plt.xlabel("r") 3135 plt.ylabel(r"$\\sigma/\\sigma_0$") 3136 plt.legend() 3137 plt.show() 3138 3139 References 3140 ---------- 3141 .. [1] https://ui.adsabs.harvard.edu/abs/1962AJ.....67..471K 3142 """ 3143 3144 amplitude = Parameter(default=1, bounds=(FLOAT_EPSILON, None), description="Amplitude or scaling factor") 3145 r_core = Parameter(default=1, bounds=(FLOAT_EPSILON, None), description="Core Radius") 3146 r_tide = Parameter(default=2, bounds=(FLOAT_EPSILON, None), description="Tidal Radius") 3147 3148 @property 3149 def concentration(self): 3150 """Concentration parameter of the king model""" 3151 return np.log10(np.abs(self.r_tide/self.r_core)) 3152 3153 @staticmethod 3154 def evaluate(x, amplitude, r_core, r_tide): 3155 """ 3156 Analytic King model function. 3157 """ 3158 3159 result = amplitude * r_core ** 2 * (1/np.sqrt(x ** 2 + r_core ** 2) - 3160 1/np.sqrt(r_tide ** 2 + r_core ** 2)) ** 2 3161 3162 # Set invalid r values to 0 3163 bounds = (x >= r_tide) | (x < 0) 3164 result[bounds] = result[bounds] * 0. 3165 3166 return result 3167 3168 @staticmethod 3169 def fit_deriv(x, amplitude, r_core, r_tide): 3170 """ 3171 Analytic King model function derivatives. 3172 """ 3173 d_amplitude = r_core ** 2 * (1/np.sqrt(x ** 2 + r_core ** 2) - 3174 1/np.sqrt(r_tide ** 2 + r_core ** 2)) ** 2 3175 3176 d_r_core = 2 * amplitude * r_core ** 2 * (r_core/(r_core ** 2 + r_tide ** 2) ** (3/2) - 3177 r_core/(r_core ** 2 + x ** 2) ** (3/2)) * \ 3178 (1./np.sqrt(r_core ** 2 + x ** 2) - 1./np.sqrt(r_core ** 2 + r_tide ** 2)) + \ 3179 2 * amplitude * r_core * (1./np.sqrt(r_core ** 2 + x ** 2) - 3180 1./np.sqrt(r_core ** 2 + r_tide ** 2)) ** 2 3181 3182 d_r_tide = (2 * amplitude * r_core ** 2 * r_tide * 3183 (1./np.sqrt(r_core ** 2 + x ** 2) - 3184 1./np.sqrt(r_core ** 2 + r_tide ** 2)))/(r_core ** 2 + r_tide ** 2) ** (3/2) 3185 3186 # Set invalid r values to 0 3187 bounds = (x >= r_tide) | (x < 0) 3188 d_amplitude[bounds] = d_amplitude[bounds]*0 3189 d_r_core[bounds] = d_r_core[bounds]*0 3190 d_r_tide[bounds] = d_r_tide[bounds]*0 3191 3192 return [d_amplitude, d_r_core, d_r_tide] 3193 3194 @property 3195 def bounding_box(self): 3196 """ 3197 Tuple defining the default ``bounding_box`` limits. 3198 3199 The model is not defined for r > r_tide. 3200 3201 ``(r_low, r_high)`` 3202 """ 3203 3204 return (0 * self.r_tide, 1 * self.r_tide) 3205 3206 @property 3207 def input_units(self): 3208 if self.r_core.unit is None: 3209 return None 3210 return {self.inputs[0]: self.r_core.unit} 3211 3212 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 3213 return {'r_core': inputs_unit[self.inputs[0]], 3214 'r_tide': inputs_unit[self.inputs[0]], 3215 'amplitude': outputs_unit[self.outputs[0]]} 3216 3217 3218class Logarithmic1D(Fittable1DModel): 3219 """ 3220 One dimensional logarithmic model. 3221 3222 Parameters 3223 ---------- 3224 amplitude : float, optional 3225 tau : float, optional 3226 3227 See Also 3228 -------- 3229 Exponential1D, Gaussian1D 3230 """ 3231 3232 amplitude = Parameter(default=1) 3233 tau = Parameter(default=1) 3234 3235 @staticmethod 3236 def evaluate(x, amplitude, tau): 3237 return amplitude * np.log(x / tau) 3238 3239 @staticmethod 3240 def fit_deriv(x, amplitude, tau): 3241 d_amplitude = np.log(x / tau) 3242 d_tau = np.zeros(x.shape) - (amplitude / tau) 3243 return [d_amplitude, d_tau] 3244 3245 @property 3246 def inverse(self): 3247 new_amplitude = self.tau 3248 new_tau = self.amplitude 3249 return Exponential1D(amplitude=new_amplitude, tau=new_tau) 3250 3251 @tau.validator 3252 def tau(self, val): 3253 if np.all(val == 0): 3254 raise ValueError("0 is not an allowed value for tau") 3255 3256 @property 3257 def input_units(self): 3258 if self.tau.unit is None: 3259 return None 3260 return {self.inputs[0]: self.tau.unit} 3261 3262 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 3263 return {'tau': inputs_unit[self.inputs[0]], 3264 'amplitude': outputs_unit[self.outputs[0]]} 3265 3266 3267class Exponential1D(Fittable1DModel): 3268 """ 3269 One dimensional exponential model. 3270 3271 Parameters 3272 ---------- 3273 amplitude : float, optional 3274 tau : float, optional 3275 3276 See Also 3277 -------- 3278 Logarithmic1D, Gaussian1D 3279 """ 3280 amplitude = Parameter(default=1) 3281 tau = Parameter(default=1) 3282 3283 @staticmethod 3284 def evaluate(x, amplitude, tau): 3285 return amplitude * np.exp(x / tau) 3286 3287 @staticmethod 3288 def fit_deriv(x, amplitude, tau): 3289 ''' Derivative with respect to parameters''' 3290 d_amplitude = np.exp(x / tau) 3291 d_tau = -amplitude * (x / tau**2) * np.exp(x / tau) 3292 return [d_amplitude, d_tau] 3293 3294 @property 3295 def inverse(self): 3296 new_amplitude = self.tau 3297 new_tau = self.amplitude 3298 return Logarithmic1D(amplitude=new_amplitude, tau=new_tau) 3299 3300 @tau.validator 3301 def tau(self, val): 3302 ''' tau cannot be 0''' 3303 if np.all(val == 0): 3304 raise ValueError("0 is not an allowed value for tau") 3305 3306 @property 3307 def input_units(self): 3308 if self.tau.unit is None: 3309 return None 3310 return {self.inputs[0]: self.tau.unit} 3311 3312 def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): 3313 return {'tau': inputs_unit[self.inputs[0]], 3314 'amplitude': outputs_unit[self.outputs[0]]} 3315