mxlib
c++ tools for analyzing astronomical data and other tasks by Jared R. Males. [git repo]
moffat.hpp
Go to the documentation of this file.
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 
33 namespace mx
34 {
35 namespace math
36 {
37 namespace func
38 {
39 
40 /** \addtogroup gen_math_moffats
41  * The Moffat Function\cite moffat_1969, a.k.a. the Moffat Profile, a.k.a. the Moffat Distribution, has the form
42  * \f[
43  I(x) = I_{pk}\left[ 1 + \frac{x^2}{\alpha^2}\right]^{-\beta}
44  * \f]
45  * With \f$\beta=1\f$ it is the
46  * Lorentzian or Cauchy distribution. See also https://en.wikipedia.org/wiki/Moffat_distribution and
47  * https://en.wikipedia.org/wiki/Cauchy_distribution.
48  *
49  * 1-D and 2-D symmetric forms are provided. Utilities are provided for normalizing and calculating the full-width at 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  * \test Scenario: compiling 1D Moffat function \ref tests_math_func_moffat1D "[test doc]"
69  *
70  * \ingroup gen_math_moffats
71  */
72 template<typename realT>
73 realT moffat( const realT x, ///< [in] is the x-position at which to evaluate the Moffat function
74  const realT I0, ///< [in] is the constant to add to the Moffat function
75  const realT Ipk, ///< [in] is the scaling factor (peak = A)
76  const realT x0, ///< [in] is the x-coordinate of the center
77  const realT alpha, ///< [in] is the width parameter of the Moffat function.
78  const realT beta ///< [in] is the shape parameter of the Moffat function
79  )
80 {
81  return I0 + Ipk * pow( static_cast<realT>(1) + pow(x-x0,2)/pow(alpha,2), -beta);
82 }
83 
84 extern template
85 float moffat<float>(const float x, const float I0, const float Ipk, const float x0, const float alpha, const float beta);
86 
87 extern template
88 double moffat<double>(const double x, const double I0, const double Ipk, const double x0, const double alpha, const double beta);
89 
90 extern template
91 long double moffat<long double>(const long double x, const long double I0, const long double Ipk, const long double x0, const long double alpha, const long double beta);
92 
93 #ifdef HASQUAD
94 extern template
95 __float128 moffat<__float128>(const __float128 x, const __float128 I0, const __float128 Ipk, const __float128 x0, const __float128 alpha, const __float128 beta);
96 #endif
97 
98 /// Find value at position (x,y) of the 2D arbitrarily-centered unnormalized symmetric Moffat function
99 /** The Moffat distribution is due to \cite moffat_1969. With \f$\beta=1\f$ it is the
100  * Lorentzian or Cauchy distribution. See also https://en.wikipedia.org/wiki/Moffat_distribution and
101  * https://en.wikipedia.org/wiki/Cauchy_distribution.
102  *
103  * Here we use the unnormalized general form, most useful for peak fitting.
104  *
105  * This function omputes:
106  * \f[
107  I(x) = I_0 + I_{pk}\left[ 1 + \frac{(x-x_0)^2 + (y-y_0)^2}{\alpha^2}\right]^{-\beta}
108  * \f]
109  *
110  * \returns the value of the 2D arbitrarily-centered symmetric Moffat function at (x,y)
111  *
112  * \tparam realT is type to use for arithmetic
113  *
114  * \test Scenario: compiling 2D Moffat function \ref tests_math_func_moffat2D "[test doc]"
115  *
116  * \ingroup gen_math_moffats
117  */
118 template<typename realT>
119 realT moffat2D( const realT x, ///< [in] the x-position at which to evaluate the Moffat function
120  const realT y, ///< [in] the y-positoin at which to evaluate the Moffat function
121  const realT I0, ///< [in] the constant to add to the Moffat function
122  const realT Ipk, ///< [in] the scaling factor (peak height is A-G0)
123  const realT x0, ///< [in] the x-coordinate of the center
124  const realT y0, ///< [in] the y-coordinate of the center
125  const realT alpha, ///< [in] the width parameter of the Moffat function.
126  const realT beta ///< [in] the shape parameter of the Moffat function.
127  )
128 {
129  return I0 + Ipk * pow( static_cast<realT>(1) + (pow(x-x0,2)+pow(y-y0,2))/pow(alpha,2), -beta);
130 }
131 
132 extern template
133 float moffat2D<float>(const float x, const float y, const float I0, const float Ipk, const float x0, const float y0, const float alpha, const float beta);
134 
135 extern template
136 double moffat2D<double>(const double x, const double y, const double I0, const double Ipk, const double x0, const double y0, const double alpha, const double beta);
137 
138 extern template
139 long double moffat2D<long double>(const long double x, const long double y, const long double I0, const long double Ipk, const long double x0, const long double y0, const long double alpha, const long double beta);
140 
141 #ifdef HASQUAD
142 extern template
143 __float128 moffat2D<__float128>(const __float128 x, const __float128 y, const __float128 I0, const __float128 Ipk, const __float128 x0, const __float128 y0, const __float128 alpha, const __float128 beta);
144 #endif
145 
146 /// Compute the full-width at half-maximum of a Moffat profile
147 /** This returns the value of
148  * \f[
149  * FWHM = 2 \alpha \sqrt{2^{1/\beta} - 1}
150  * \f]
151  *
152  * \returns the FWHM of the Moffat profile
153  *
154  * \tparam realT is the type to use for arithmetic
155  *
156  * \test Scenario: compiling Moffat FWHM \ref tests_math_func_moffatFWHM "[test doc]"
157  *
158  * \ingroup gen_math_moffats
159  */
160 template<typename realT>
161 realT moffatFWHM( realT alpha, ///< [in] the width parameter of the Moffat function.
162  realT beta ///< [in] the shape parameter of the Moffat function.
163  )
164 {
165  return 2*alpha*sqrt( pow(static_cast<realT>(2), static_cast<realT>(1)/beta) - 1);
166 }
167 
168 extern template
169 float moffatFWHM(float alpha, float beta);
170 
171 extern template
172 double moffatFWHM(double alpha, double beta);
173 
174 extern template
175 long double moffatFWHM(long double alpha, long double beta);
176 
177 #ifdef HASQUAD
178 extern template
179 __float128 moffatFWHM(__float128 alpha, __float128 beta);
180 #endif
181 
182 } //namespace func
183 } //namespace math
184 } //namespace mx
185 
186 #endif //math_func_moffat_hpp
187 
188 
189 
realT moffatFWHM(realT alpha, realT beta)
Compute the full-width at half-maximum of a Moffat profile.
Definition: moffat.hpp:161
realT moffat(const realT x, const realT I0, const realT Ipk, const realT x0, const realT alpha, const realT beta)
Find value at position (x) of the 1D arbitrarily-centered symmetric unnormalized Moffat function.
Definition: moffat.hpp:73
realT moffat2D(const realT x, const realT y, const realT I0, const realT Ipk, const realT x0, const realT y0, const realT alpha, const realT beta)
Find value at position (x,y) of the 2D arbitrarily-centered unnormalized symmetric Moffat function.
Definition: moffat.hpp:119
The mxlib c++ namespace.
Definition: mxError.hpp:107