Line data Source code
1 : /** \file fitMoffat.hpp
2 : * \author Jared R. Males
3 : * \brief Tools for fitting Moffat functions to data.
4 : * \ingroup fitting_files
5 : *
6 : */
7 :
8 : //***********************************************************************//
9 : // Copyright 2015, 2016, 2017 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 fitMoffat_hpp
28 : #define fitMoffat_hpp
29 :
30 : #include "levmarInterface.hpp"
31 : #include "../func/moffat.hpp"
32 :
33 : #include "../../improc/eigenImage.hpp"
34 :
35 : namespace mx
36 : {
37 : namespace math
38 : {
39 : namespace fit
40 : {
41 :
42 : /// Wrapper for a native array to pass to \ref levmarInterface, with Moffat details.
43 : /** \ingroup moffat_peak_fit
44 : */
45 : template <typename realT>
46 : struct array2FitMoffat1D
47 : {
48 : realT *data{ nullptr }; ///< Pointer to the array of y values
49 : realT *coords{ nullptr }; ///< Pointer to the array of x values (optional)
50 : size_t nx{ 0 }; ///< X dimension of the array
51 :
52 : realT m_I0{ 0 };
53 : realT m_I{ 0 };
54 : realT m_x0{ 0 };
55 : realT m_alpha{ 0 };
56 : realT m_beta{ 0 };
57 :
58 : int m_I0_idx{ 0 };
59 : int m_I_idx{ 1 };
60 : int m_x0_idx{ 2 };
61 : int m_alpha_idx{ 4 };
62 : int m_beta_idx{ 5 };
63 :
64 : int m_nparams = 5;
65 :
66 : /// Set whether each parameter is fixed.
67 : /** Sets the parameter indices appropriately.
68 : */
69 0 : void setFixed( bool I0, bool I, bool x0, bool alpha, bool beta )
70 : {
71 0 : int idx = 0;
72 :
73 0 : if( I0 )
74 0 : m_I0_idx = -1;
75 : else
76 0 : m_I0_idx = idx++;
77 :
78 0 : if( I )
79 0 : m_I_idx = -1;
80 : else
81 0 : m_I_idx = idx++;
82 :
83 0 : if( x0 )
84 0 : m_x0_idx = -1;
85 : else
86 0 : m_x0_idx = idx++;
87 :
88 0 : if( alpha )
89 0 : m_alpha_idx = -1;
90 : else
91 0 : m_alpha_idx = idx++;
92 :
93 0 : if( beta )
94 0 : m_beta_idx = -1;
95 : else
96 0 : m_beta_idx = idx++;
97 :
98 0 : m_nparams = idx;
99 0 : }
100 :
101 0 : realT I0( realT *p )
102 : {
103 0 : if( m_I0_idx < 0 )
104 : {
105 0 : return m_I0;
106 : }
107 : else
108 : {
109 0 : return p[m_I0_idx];
110 : }
111 : }
112 :
113 0 : void I0( realT *p, realT nI0 )
114 : {
115 0 : if( m_I0_idx < 0 )
116 : {
117 0 : m_I0 = nI0;
118 : }
119 : else
120 : {
121 0 : p[m_I0_idx] = nI0;
122 : }
123 0 : }
124 :
125 0 : realT I( realT *p )
126 : {
127 0 : if( m_I_idx < 0 )
128 : {
129 0 : return m_I;
130 : }
131 : else
132 : {
133 0 : return p[m_I_idx];
134 : }
135 : }
136 :
137 0 : void I( realT *p, realT nI )
138 : {
139 0 : if( m_I_idx < 0 )
140 : {
141 0 : m_I = nI;
142 : }
143 : else
144 : {
145 0 : p[m_I_idx] = nI;
146 : }
147 0 : }
148 :
149 0 : realT x0( realT *p )
150 : {
151 0 : if( m_x0_idx < 0 )
152 : {
153 0 : return m_x0;
154 : }
155 : else
156 : {
157 0 : return p[m_x0_idx];
158 : }
159 : }
160 :
161 0 : void x0( realT *p, realT nx0 )
162 : {
163 0 : if( m_x0_idx < 0 )
164 : {
165 0 : m_x0 = nx0;
166 : }
167 : else
168 : {
169 0 : p[m_x0_idx] = nx0;
170 : }
171 0 : }
172 :
173 0 : realT alpha( realT *p )
174 : {
175 0 : if( m_alpha_idx < 0 )
176 : {
177 0 : return m_alpha;
178 : }
179 : else
180 : {
181 0 : return p[m_alpha_idx];
182 : }
183 : }
184 :
185 0 : void alpha( realT *p, realT nalpha )
186 : {
187 0 : if( m_alpha_idx < 0 )
188 : {
189 0 : m_alpha = nalpha;
190 : }
191 : else
192 : {
193 0 : p[m_alpha_idx] = nalpha;
194 : }
195 0 : }
196 :
197 0 : realT beta( realT *p )
198 : {
199 0 : if( m_beta_idx < 0 )
200 : {
201 0 : return m_beta;
202 : }
203 : else
204 : {
205 0 : return p[m_beta_idx];
206 : }
207 : }
208 :
209 0 : void beta( realT *p, realT nbeta )
210 : {
211 0 : if( m_beta_idx < 0 )
212 : {
213 0 : m_beta = nbeta;
214 : }
215 : else
216 : {
217 0 : p[m_beta_idx] = nbeta;
218 : }
219 0 : }
220 :
221 0 : int nparams()
222 : {
223 0 : return m_nparams;
224 : }
225 : };
226 :
227 : /** \defgroup moffat_peak_fit Moffat Functions
228 : * \brief Fitting Moffat functions to data.
229 : *
230 : * The Moffat Function\cite moffat_1969, a.k.a. the Moffat Profile, a.k.a. the Moffat Distribution, has the form
231 : * \f[
232 : I(x) = I_{pk}\left[ 1 + \frac{x^2}{\alpha^2}\right]^{-\beta}
233 : * \f]
234 : * With \f$\beta=1\f$ it is the
235 : * Lorentzian or Cauchy distribution. See also https://en.wikipedia.org/wiki/Moffat_distribution and
236 : * https://en.wikipedia.org/wiki/Cauchy_distribution.
237 : *
238 : * For utilities for working with Moffat functions see \ref gen_math_moffats.
239 : *
240 : * \ingroup peak_fit
241 : */
242 :
243 : ///\ref levmarInterface fitter structure for the 1D Moffat.
244 : /** \ingroup moffat_peak_fit
245 : *
246 : */
247 : template <typename _realT>
248 : struct moffat1D_fitter
249 : {
250 : typedef _realT realT;
251 :
252 : static const int nparams = 6;
253 :
254 0 : static void func( realT *p, realT *hx, int m, int n, void *adata )
255 : {
256 0 : array2FitMoffat1D<realT> *arr = (array2FitMoffat1D<realT> *)adata;
257 :
258 : size_t idx_mat, idx_dat;
259 :
260 0 : idx_dat = 0;
261 :
262 0 : realT I0 = arr->I0( p );
263 0 : realT I = arr->I( p );
264 0 : realT x0 = arr->x0( p );
265 0 : realT alpha = arr->alpha( p );
266 0 : realT beta = arr->beta( p );
267 :
268 0 : if( arr->coords )
269 : {
270 0 : for( int i = 0; i < arr->nx; ++i )
271 : {
272 0 : hx[i] = func::moffat<realT>( arr->coords[i], I0, I, x0, alpha, beta ) - arr->data[i];
273 : }
274 : }
275 : else
276 : {
277 0 : for( int i = 0; i < arr->nx; ++i )
278 : {
279 0 : hx[i] = func::moffat<realT>( i, I0, I, x0, alpha, beta ) - arr->data[i];
280 : }
281 : }
282 0 : }
283 : };
284 :
285 : /// Class to manage fitting a 1D Moffat to data via the \ref levmarInterface
286 : /** Fits \ref gen_math_moffats to a 1-dimensional array of data.
287 : *
288 : * This class allows for treating any of the parameters as fixed.
289 : *
290 : * \tparam fitterT a type meeting the requirements specified in \ref levmarInterface.
291 : *
292 : * \ingroup moffat_peak_fit
293 : *
294 : */
295 : template <typename _realT>
296 : class fitMoffat1D : public levmarInterface<moffat1D_fitter<_realT>>
297 : {
298 :
299 : public:
300 : typedef _realT realT;
301 : typedef moffat1D_fitter<realT> fitterT;
302 :
303 : protected:
304 : array2FitMoffat1D<realT> arr;
305 :
306 0 : void initialize()
307 : {
308 0 : this->allocate_params( arr.nparams() );
309 0 : this->adata = &arr;
310 0 : }
311 :
312 : public:
313 0 : fitMoffat1D()
314 0 : {
315 0 : initialize();
316 0 : }
317 :
318 0 : ~fitMoffat1D()
319 : {
320 0 : }
321 :
322 : /// Set whether each parameter is fixed.
323 : /** Sets the parameter indices appropriately.
324 : */
325 0 : void setFixed( bool I0, ///< [in] if true, then I0 will be not be part of the fit
326 : bool I, ///< [in] if true, then I will be not be part of the fit
327 : bool x0, ///< [in] if true, then x0 will be not be part of the fit
328 : bool alpha, ///< [in] if true, then alpha will be not be part of the fit
329 : bool beta ///< [in] if true, then beta will be not be part of the fit
330 : )
331 : {
332 0 : arr.setFixed( I0, I, x0, alpha, beta );
333 0 : this->allocate_params( arr.nparams() );
334 0 : }
335 :
336 : /// Set the initial guess for a symmetric Moffat.
337 0 : void setGuess( realT I0, ///< [in] the constant background level
338 : realT I, ///< [in] the peak scaling
339 : realT x0, ///< [in] the center x-coordinate
340 : realT alpha, ///< [in] the width parameter
341 : realT beta ///< [in] the shape parameter
342 : )
343 : {
344 0 : arr.I0( this->p, I0 );
345 0 : arr.I( this->p, I );
346 0 : arr.x0( this->p, x0 );
347 0 : arr.alpha( this->p, alpha );
348 0 : arr.beta( this->p, beta );
349 0 : }
350 :
351 : /// Set the data aray.
352 0 : void setArray( realT *data, ///< [in] pointer to an nx sized array of the y-values to be fit
353 : int nx ///< [in] the number of pixels in the x direction of the data array
354 : )
355 : {
356 0 : arr.data = data;
357 0 : arr.nx = nx;
358 :
359 0 : this->n = nx;
360 0 : }
361 :
362 : /// Set the data aray and the coordinates.
363 0 : void setArray( realT *data, ///< [in] pointer to an nx sized array of the y-values to be fit
364 : realT *coords, ///< [in] point to an nx sized array of the x-values of the data
365 : int nx ///< [in] the number of pixels in the x direction of the data array
366 : )
367 : {
368 0 : arr.data = data;
369 0 : arr.coords = coords;
370 0 : arr.nx = nx;
371 :
372 0 : this->n = nx;
373 0 : }
374 :
375 : /// Get the current value of I0, the constant.
376 : /**
377 : * \returns the current value of I0
378 : */
379 0 : realT I0()
380 : {
381 0 : return arr.I0( this->p );
382 : }
383 :
384 : /// Get the current value of I, the peak scaling.
385 : /**
386 : * \returns the current value of I
387 : */
388 0 : realT I()
389 : {
390 0 : return arr.I( this->p );
391 : }
392 :
393 : /// Get the center x-
394 : /**
395 : * \returns the current value of x0
396 : */
397 0 : realT x0()
398 : {
399 0 : return arr.x0( this->p );
400 : }
401 :
402 : /// Return the width parameter
403 : /**
404 : * \returns the current value of alpha
405 : */
406 0 : realT alpha()
407 : {
408 0 : return arr.alpha( this->p );
409 : }
410 :
411 : /// Return the shape parameter
412 : /**
413 : * \returns the current value of beta
414 : */
415 0 : realT beta()
416 : {
417 0 : return arr.beta( this->p );
418 : }
419 :
420 : /// Return the full-width at half maximum
421 : /**
422 : * \returns the FWHM calculated from alpha and beta
423 : */
424 0 : realT fwhm()
425 : {
426 0 : return func::moffatFWHM( alpha(), beta() );
427 : }
428 : };
429 :
430 : // forward
431 : template <typename realT>
432 : struct array2FitMoffat;
433 :
434 : // forward
435 : template <typename _realT>
436 : struct moffat2D_sym_fitter;
437 :
438 : /// Class to manage fitting a 2D Moffat to data via the \ref levmarInterface
439 : /** Fits \ref gen_math_moffats to a 2-dimensional array of data.
440 : *
441 : * This class allows for treating any of the parameters as fixed.
442 : *
443 : * \tparam fitterT a type meeting the requirements specified in \ref levmarInterface.
444 : *
445 : * \ingroup moffat_peak_fit
446 : *
447 : */
448 : template <typename fitterT>
449 : class fitMoffat2D : public levmarInterface<fitterT>
450 : {
451 :
452 : public:
453 : typedef typename fitterT::realT realT;
454 :
455 : protected:
456 : array2FitMoffat<realT> arr;
457 :
458 0 : void initialize()
459 : {
460 0 : this->allocate_params( arr.nparams() );
461 0 : this->adata = &arr;
462 0 : }
463 :
464 : public:
465 0 : fitMoffat2D()
466 0 : {
467 0 : initialize();
468 0 : }
469 :
470 0 : ~fitMoffat2D()
471 : {
472 0 : }
473 :
474 : /// Set whether each parameter is fixed.
475 : /** Sets the parameter indices appropriately.
476 : */
477 0 : void setFixed( bool I0, ///< [in] if true, then I0 will be not be part of the fit
478 : bool I, ///< [in] if true, then I will be not be part of the fit
479 : bool x0, ///< [in] if true, then x0 will be not be part of the fit
480 : bool y0, ///< [in] if true, then y0 will be not be part of the fit
481 : bool alpha, ///< [in] if true, then alpha will be not be part of the fit
482 : bool beta ///< [in] if true, then beta will be not be part of the fit
483 : )
484 : {
485 0 : arr.setFixed( I0, I, x0, y0, alpha, beta );
486 0 : this->allocate_params( arr.nparams() );
487 0 : }
488 :
489 : /// Set the initial guess for a symmetric Moffat.
490 0 : void setGuess( realT I0, ///< [in] the constant background level
491 : realT I, ///< [in] the peak scaling
492 : realT x0, ///< [in] the center x-coordinate
493 : realT y0, ///< [in] the center y-coordinate
494 : realT alpha, ///< [in] the width parameter
495 : realT beta ///< [in] the shape parameter
496 : )
497 : {
498 0 : arr.I0( this->p, I0 );
499 0 : arr.I( this->p, I );
500 0 : arr.x0( this->p, x0 );
501 0 : arr.y0( this->p, y0 );
502 0 : arr.alpha( this->p, alpha );
503 0 : arr.beta( this->p, beta );
504 0 : }
505 :
506 : /// Set the data aray.
507 0 : void setArray( realT *data, ///< [in] pointer to an nx X ny array of data to be fit
508 : int nx, ///< [in] the number of pixels in the x direction of the data array
509 : int ny ///< [in] the number of pixels in the y direction of the data array
510 : )
511 : {
512 0 : arr.data = data;
513 0 : arr.nx = nx;
514 0 : arr.ny = ny;
515 :
516 0 : this->n = nx * ny;
517 0 : }
518 :
519 : /// Get the current value of I0, the constant.
520 : /**
521 : * \returns the current value of I0
522 : */
523 0 : realT I0()
524 : {
525 0 : return arr.I0( this->p );
526 : }
527 :
528 : /// Get the current value of I, the peak scaling.
529 : /**
530 : * \returns the current value of I
531 : */
532 0 : realT I()
533 : {
534 0 : return arr.I( this->p );
535 : }
536 :
537 : /// Get the center x-
538 : /**
539 : * \returns the current value of x0
540 : */
541 0 : realT x0()
542 : {
543 0 : return arr.x0( this->p );
544 : }
545 :
546 : /// Get the center y-coordinate
547 : /**
548 : * \returns the current value of y0
549 : */
550 0 : realT y0()
551 : {
552 0 : return arr.y0( this->p );
553 : }
554 :
555 : /// Return the width parameter
556 : /**
557 : * \returns the current value of alpha
558 : */
559 0 : realT alpha()
560 : {
561 0 : return arr.alpha( this->p );
562 : }
563 :
564 : /// Return the shape parameter
565 : /**
566 : * \returns the current value of beta
567 : */
568 0 : realT beta()
569 : {
570 0 : return arr.beta( this->p );
571 : }
572 :
573 : /// Return the full-width at half maximum
574 : /**
575 : * \returns the FWHM calculated from alpha and beta
576 : */
577 0 : realT fwhm()
578 : {
579 0 : return func::moffatFWHM( alpha(), beta() );
580 : }
581 : };
582 :
583 : /// Wrapper for a native array to pass to \ref levmarInterface, with Moffat details.
584 : /** \ingroup moffat_peak_fit
585 : */
586 : template <typename realT>
587 : struct array2FitMoffat
588 : {
589 : realT *data{ nullptr }; ///< Pointer to the array
590 : size_t nx{ 0 }; ///< X dimension of the array
591 : size_t ny{ 0 }; ///< Y dimension of the array
592 :
593 : realT m_I0{ 0 };
594 : realT m_I{ 0 };
595 : realT m_x0{ 0 };
596 : realT m_y0{ 0 };
597 : realT m_alpha{ 0 };
598 : realT m_beta{ 0 };
599 :
600 : int m_I0_idx{ 0 };
601 : int m_I_idx{ 1 };
602 : int m_x0_idx{ 2 };
603 : int m_y0_idx{ 3 };
604 : int m_alpha_idx{ 4 };
605 : int m_beta_idx{ 5 };
606 :
607 : int m_nparams = 6;
608 :
609 : /// Set whether each parameter is fixed.
610 : /** Sets the parameter indices appropriately.
611 : */
612 0 : void setFixed( bool I0, bool I, bool x0, bool y0, bool alpha, bool beta )
613 : {
614 0 : int idx = 0;
615 :
616 0 : if( I0 )
617 0 : m_I0_idx = -1;
618 : else
619 0 : m_I0_idx = idx++;
620 :
621 0 : if( I )
622 0 : m_I_idx = -1;
623 : else
624 0 : m_I_idx = idx++;
625 :
626 0 : if( x0 )
627 0 : m_x0_idx = -1;
628 : else
629 0 : m_x0_idx = idx++;
630 :
631 0 : if( y0 )
632 0 : m_y0_idx = -1;
633 : else
634 0 : m_y0_idx = idx++;
635 :
636 0 : if( alpha )
637 0 : m_alpha_idx = -1;
638 : else
639 0 : m_alpha_idx = idx++;
640 :
641 0 : if( beta )
642 0 : m_beta_idx = -1;
643 : else
644 0 : m_beta_idx = idx++;
645 :
646 0 : m_nparams = idx;
647 0 : }
648 :
649 0 : realT I0( realT *p )
650 : {
651 0 : if( m_I0_idx < 0 )
652 : {
653 0 : return m_I0;
654 : }
655 : else
656 : {
657 0 : return p[m_I0_idx];
658 : }
659 : }
660 :
661 0 : void I0( realT *p, realT nI0 )
662 : {
663 0 : if( m_I0_idx < 0 )
664 : {
665 0 : m_I0 = nI0;
666 : }
667 : else
668 : {
669 0 : p[m_I0_idx] = nI0;
670 : }
671 0 : }
672 :
673 0 : realT I( realT *p )
674 : {
675 0 : if( m_I_idx < 0 )
676 : {
677 0 : return m_I;
678 : }
679 : else
680 : {
681 0 : return p[m_I_idx];
682 : }
683 : }
684 :
685 0 : void I( realT *p, realT nI )
686 : {
687 0 : if( m_I_idx < 0 )
688 : {
689 0 : m_I = nI;
690 : }
691 : else
692 : {
693 0 : p[m_I_idx] = nI;
694 : }
695 0 : }
696 :
697 0 : realT x0( realT *p )
698 : {
699 0 : if( m_x0_idx < 0 )
700 : {
701 0 : return m_x0;
702 : }
703 : else
704 : {
705 0 : return p[m_x0_idx];
706 : }
707 : }
708 :
709 0 : void x0( realT *p, realT nx0 )
710 : {
711 0 : if( m_x0_idx < 0 )
712 : {
713 0 : m_x0 = nx0;
714 : }
715 : else
716 : {
717 0 : p[m_x0_idx] = nx0;
718 : }
719 0 : }
720 :
721 0 : realT y0( realT *p )
722 : {
723 0 : if( m_y0_idx < 0 )
724 : {
725 0 : return m_y0;
726 : }
727 : else
728 : {
729 0 : return p[m_y0_idx];
730 : }
731 : }
732 :
733 0 : void y0( realT *p, realT ny0 )
734 : {
735 0 : if( m_y0_idx < 0 )
736 : {
737 0 : m_y0 = ny0;
738 : }
739 : else
740 : {
741 0 : p[m_y0_idx] = ny0;
742 : }
743 0 : }
744 :
745 0 : realT alpha( realT *p )
746 : {
747 0 : if( m_alpha_idx < 0 )
748 : {
749 0 : return m_alpha;
750 : }
751 : else
752 : {
753 0 : return p[m_alpha_idx];
754 : }
755 : }
756 :
757 0 : void alpha( realT *p, realT nalpha )
758 : {
759 0 : if( m_alpha_idx < 0 )
760 : {
761 0 : m_alpha = nalpha;
762 : }
763 : else
764 : {
765 0 : p[m_alpha_idx] = nalpha;
766 : }
767 0 : }
768 :
769 0 : realT beta( realT *p )
770 : {
771 0 : if( m_beta_idx < 0 )
772 : {
773 0 : return m_beta;
774 : }
775 : else
776 : {
777 0 : return p[m_beta_idx];
778 : }
779 : }
780 :
781 0 : void beta( realT *p, realT nbeta )
782 : {
783 0 : if( m_beta_idx < 0 )
784 : {
785 0 : m_beta = nbeta;
786 : }
787 : else
788 : {
789 0 : p[m_beta_idx] = nbeta;
790 : }
791 0 : }
792 :
793 0 : int nparams()
794 : {
795 0 : return m_nparams;
796 : }
797 : };
798 :
799 : ///\ref levmarInterface fitter structure for the symmetric Moffat.
800 : /** \ingroup moffat_peak_fit
801 : *
802 : */
803 : template <typename _realT>
804 : struct moffat2D_sym_fitter
805 : {
806 : typedef _realT realT;
807 :
808 : static const int nparams = 6;
809 :
810 0 : static void func( realT *p, realT *hx, int m, int n, void *adata )
811 : {
812 0 : array2FitMoffat<realT> *arr = (array2FitMoffat<realT> *)adata;
813 :
814 : size_t idx_mat, idx_dat;
815 :
816 0 : idx_dat = 0;
817 :
818 0 : realT I0 = arr->I0( p );
819 0 : realT I = arr->I( p );
820 0 : realT x0 = arr->x0( p );
821 0 : realT y0 = arr->y0( p );
822 0 : realT alpha = arr->alpha( p );
823 0 : realT beta = arr->beta( p );
824 :
825 0 : for( int i = 0; i < arr->nx; ++i )
826 : {
827 0 : for( int j = 0; j < arr->ny; ++j )
828 : {
829 0 : idx_mat = i + j * arr->nx;
830 :
831 0 : hx[idx_dat] = func::moffat2D<realT>( i, j, I0, I, x0, y0, alpha, beta ) - arr->data[idx_mat];
832 :
833 0 : ++idx_dat;
834 : }
835 : }
836 0 : }
837 : };
838 :
839 : /// Alias for the fitMoffat2D type fitting the symmetric Moffat profile.
840 : /** \ingroup moffat_peak_fit
841 : */
842 : template <typename realT>
843 : using fitMoffat2Dsym = mx::math::fit::fitMoffat2D<mx::math::fit::moffat2D_sym_fitter<realT>>;
844 :
845 : } // namespace fit
846 : } // namespace math
847 :
848 : } // namespace mx
849 :
850 : #endif // fitMoffat_hpp
|
