Line data Source code
1 : /** \file moffat.hpp
2 : * \author Jared R. Males
3 : * \brief Declarations for utilities related to the Moffat function.
4 : * \ingroup gen_math_files
5 : *
6 : */
7 :
8 : //***********************************************************************//
9 : // Copyright 2020 Jared R. Males (jaredmales@gmail.com)
10 : //
11 : // This file is part of mxlib.
12 : //
13 : // mxlib is free software: you can redistribute it and/or modify
14 : // it under the terms of the GNU General Public License as published by
15 : // the Free Software Foundation, either version 3 of the License, or
16 : // (at your option) any later version.
17 : //
18 : // mxlib is distributed in the hope that it will be useful,
19 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
20 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 : // GNU General Public License for more details.
22 : //
23 : // You should have received a copy of the GNU General Public License
24 : // along with mxlib. If not, see <http://www.gnu.org/licenses/>.
25 : //***********************************************************************//
26 :
27 : #ifndef math_func_moffat_hpp
28 : #define math_func_moffat_hpp
29 :
30 : #include <cmath>
31 :
32 : namespace mx
33 : {
34 : namespace math
35 : {
36 : namespace func
37 : {
38 :
39 : /** \addtogroup gen_math_moffats
40 : * The Moffat Function\cite moffat_1969, a.k.a. the Moffat Profile, a.k.a. the Moffat Distribution, has the form
41 : * \f[
42 : I(x) = I_{pk}\left[ 1 + \frac{x^2}{\alpha^2}\right]^{-\beta}
43 : * \f]
44 : * With \f$\beta=1\f$ it is the
45 : * Lorentzian or Cauchy distribution. See also https://en.wikipedia.org/wiki/Moffat_distribution and
46 : * https://en.wikipedia.org/wiki/Cauchy_distribution.
47 : *
48 : * 1-D and 2-D symmetric forms are provided. Utilities are provided for normalizing and calculating the full-width at
49 : half-maximum.
50 : */
51 :
52 : /// Find value at position (x) of the 1D arbitrarily-centered symmetric unnormalized Moffat function
53 : /** The Moffat distribution is due to \cite moffat_1969. With \f$\beta=1\f$ it is the
54 : * Lorentzian or Cauchy distribution. See also https://en.wikipedia.org/wiki/Moffat_distribution and
55 : * https://en.wikipedia.org/wiki/Cauchy_distribution.
56 : *
57 : * Here we use the unnormalized general form, most useful for peak fitting.
58 : *
59 : * This function computes:
60 : * \f[
61 : I(x) = I_0 + I_{pk}\left[ 1 + \frac{(x-x_0)^2}{\alpha^2}\right]^{-\beta}
62 : * \f]
63 : *
64 : * \returns the value of the 1D arbitrarily-centered unnormalized Moffat at (x)
65 : *
66 : * \tparam realT is type to use for arithmetic
67 : *
68 : *
69 : * \ingroup gen_math_moffats
70 : */
71 : template <typename realT>
72 1 : realT moffat( const realT x, ///< [in] is the x-position at which to evaluate the Moffat function
73 : const realT I0, ///< [in] is the constant to add to the Moffat function
74 : const realT Ipk, ///< [in] is the scaling factor (peak = A)
75 : const realT x0, ///< [in] is the x-coordinate of the center
76 : const realT alpha, ///< [in] is the width parameter of the Moffat function.
77 : const realT beta ///< [in] is the shape parameter of the Moffat function
78 : )
79 : {
80 1 : return I0 + Ipk * pow( static_cast<realT>( 1 ) + pow( x - x0, 2 ) / pow( alpha, 2 ), -beta );
81 : }
82 :
83 : extern template float
84 : moffat<float>( const float x, const float I0, const float Ipk, const float x0, const float alpha, const float beta );
85 :
86 : extern template double moffat<double>(
87 : const double x, const double I0, const double Ipk, const double x0, const double alpha, const double beta );
88 :
89 : extern template long double moffat<long double>( const long double x,
90 : const long double I0,
91 : const long double Ipk,
92 : const long double x0,
93 : const long double alpha,
94 : const long double beta );
95 :
96 : #ifdef HASQUAD
97 : extern template __float128 moffat<__float128>( const __float128 x,
98 : const __float128 I0,
99 : const __float128 Ipk,
100 : const __float128 x0,
101 : const __float128 alpha,
102 : const __float128 beta );
103 : #endif
104 :
105 : /// Find value at position (x,y) of the 2D arbitrarily-centered unnormalized symmetric Moffat function
106 : /** The Moffat distribution is due to \cite moffat_1969. With \f$\beta=1\f$ it is the
107 : * Lorentzian or Cauchy distribution. See also https://en.wikipedia.org/wiki/Moffat_distribution and
108 : * https://en.wikipedia.org/wiki/Cauchy_distribution.
109 : *
110 : * Here we use the unnormalized general form, most useful for peak fitting.
111 : *
112 : * This function omputes:
113 : * \f[
114 : I(x) = I_0 + I_{pk}\left[ 1 + \frac{(x-x_0)^2 + (y-y_0)^2}{\alpha^2}\right]^{-\beta}
115 : * \f]
116 : *
117 : * \returns the value of the 2D arbitrarily-centered symmetric Moffat function at (x,y)
118 : *
119 : * \tparam realT is type to use for arithmetic
120 : *
121 : *
122 : * \ingroup gen_math_moffats
123 : */
124 : template <typename realT>
125 1 : realT moffat2D( const realT x, ///< [in] the x-position at which to evaluate the Moffat function
126 : const realT y, ///< [in] the y-positoin at which to evaluate the Moffat function
127 : const realT I0, ///< [in] the constant to add to the Moffat function
128 : const realT Ipk, ///< [in] the scaling factor (peak height is A-G0)
129 : const realT x0, ///< [in] the x-coordinate of the center
130 : const realT y0, ///< [in] the y-coordinate of the center
131 : const realT alpha, ///< [in] the width parameter of the Moffat function.
132 : const realT beta ///< [in] the shape parameter of the Moffat function.
133 : )
134 : {
135 1 : return I0 + Ipk * pow( static_cast<realT>( 1 ) + ( pow( x - x0, 2 ) + pow( y - y0, 2 ) ) / pow( alpha, 2 ), -beta );
136 : }
137 :
138 : extern template float moffat2D<float>( const float x,
139 : const float y,
140 : const float I0,
141 : const float Ipk,
142 : const float x0,
143 : const float y0,
144 : const float alpha,
145 : const float beta );
146 :
147 : extern template double moffat2D<double>( const double x,
148 : const double y,
149 : const double I0,
150 : const double Ipk,
151 : const double x0,
152 : const double y0,
153 : const double alpha,
154 : const double beta );
155 :
156 : extern template long double moffat2D<long double>( const long double x,
157 : const long double y,
158 : const long double I0,
159 : const long double Ipk,
160 : const long double x0,
161 : const long double y0,
162 : const long double alpha,
163 : const long double beta );
164 :
165 : #ifdef HASQUAD
166 : extern template __float128 moffat2D<__float128>( const __float128 x,
167 : const __float128 y,
168 : const __float128 I0,
169 : const __float128 Ipk,
170 : const __float128 x0,
171 : const __float128 y0,
172 : const __float128 alpha,
173 : const __float128 beta );
174 : #endif
175 :
176 : /// Compute the full-width at half-maximum of a Moffat profile
177 : /** This returns the value of
178 : * \f[
179 : * FWHM = 2 \alpha \sqrt{2^{1/\beta} - 1}
180 : * \f]
181 : *
182 : * \returns the FWHM of the Moffat profile
183 : *
184 : * \tparam realT is the type to use for arithmetic
185 : *
186 : *
187 : * \ingroup gen_math_moffats
188 : */
189 : template <typename realT>
190 1 : realT moffatFWHM( realT alpha, ///< [in] the width parameter of the Moffat function.
191 : realT beta ///< [in] the shape parameter of the Moffat function.
192 : )
193 : {
194 1 : return 2 * alpha * sqrt( pow( static_cast<realT>( 2 ), static_cast<realT>( 1 ) / beta ) - 1 );
195 : }
196 :
197 : extern template float moffatFWHM( float alpha, float beta );
198 :
199 : extern template double moffatFWHM( double alpha, double beta );
200 :
201 : extern template long double moffatFWHM( long double alpha, long double beta );
202 :
203 : #ifdef HASQUAD
204 : extern template __float128 moffatFWHM( __float128 alpha, __float128 beta );
205 : #endif
206 :
207 : } // namespace func
208 : } // namespace math
209 : } // namespace mx
210 :
211 : #endif // math_func_moffat_hpp
|
