GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelSpectralLogParabola.hpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelSpectralLogParabola.hpp - Log parabola spectral model class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2012-2018 by Michael Mayer *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GModelSpectralLogParabola.hpp
23  * @brief Log parabola spectral model class definition
24  * @author Michael Mayer
25  */
26 
27 #ifndef GMODELSPECTRALLOGPARABOLA_HPP
28 #define GMODELSPECTRALLOGPARABOLA_HPP
29 
30 
31 /* __ Includes ___________________________________________________________ */
32 #include <string>
33 #include <cmath>
34 #include "GModelSpectral.hpp"
35 #include "GFunction.hpp"
36 #include "GModelPar.hpp"
37 #include "GEnergy.hpp"
38 
39 /* __ Forward declarations _______________________________________________ */
40 class GRan;
41 class GTime;
42 class GXmlElement;
43 
44 
45 /***********************************************************************//**
46  * @class GModelSpectralLogParabola
47  *
48  * @brief LogParabola spectral model class
49  *
50  * This class implements a log parabola spectrum. The spectrum is defined by
51  *
52  * \f[
53  * S_{\rm E}(E | t) = {\tt m\_norm}
54  * \left( \frac{E}{\tt m\_pivot} \right)^{{\tt m\_index} +
55  * {\tt m\_curvature} \, \ln \frac{E}{\tt m\_pivot}}
56  * \f]
57  *
58  * where
59  * - \f${\tt m\_norm}\f$ is the normalization or prefactor,
60  * - \f${\tt m\_index}\f$ is the spectral index,
61  * - \f${\tt m\_curvature}\f$ is the spectral curvature, and
62  * - \f${\tt m\_pivot}\f$ is the pivot energy.
63  ***************************************************************************/
65 
66 public:
67  // Constructors and destructors
69  GModelSpectralLogParabola(const std::string& type,
70  const std::string& prefactor,
71  const std::string& index,
72  const std::string& pivot,
73  const std::string& curvature);
74  GModelSpectralLogParabola(const double& prefactor,
75  const double& index,
76  const GEnergy& pivot,
77  const double& curvature);
78  explicit GModelSpectralLogParabola(const GXmlElement& xml);
80  virtual ~GModelSpectralLogParabola(void);
81 
82  // Operators
84 
85  // Implemented pure virtual base class methods
86  virtual void clear(void);
87  virtual GModelSpectralLogParabola* clone(void) const;
88  virtual std::string classname(void) const;
89  virtual std::string type(void) const;
90  virtual double eval(const GEnergy& srcEng,
91  const GTime& srcTime = GTime(),
92  const bool& gradients = false) const;
93  virtual double flux(const GEnergy& emin,
94  const GEnergy& emax) const;
95  virtual double eflux(const GEnergy& emin,
96  const GEnergy& emax) const;
97  virtual GEnergy mc(const GEnergy& emin,
98  const GEnergy& emax,
99  const GTime& time,
100  GRan& ran) const;
101  virtual void read(const GXmlElement& xml);
102  virtual void write(GXmlElement& xml) const;
103  virtual std::string print(const GChatter& chatter = NORMAL) const;
104 
105  // Other methods
106  void type(const std::string& type);
107  double prefactor(void) const;
108  void prefactor(const double& prefactor);
109  double index(void) const;
110  void index(const double& index);
111  double curvature(void) const;
112  void curvature(const double& curvature);
113  GEnergy pivot(void) const;
114  void pivot(const GEnergy& pivot);
115 
116 protected:
117  // Protected methods
118  void init_members(void);
119  void copy_members(const GModelSpectralLogParabola& model);
120  void free_members(void);
121  void update_eval_cache(const GEnergy& energy) const;
122  void update_mc_cache(const GEnergy& emin, const GEnergy& emax,
123  const GTime& time) const;
124 
125  // Class to determine to the integral photon flux
126  class flux_kern : public GFunction {
127  public:
128  // Constructor
129  flux_kern(const double& norm,
130  const double& index,
131  const double& curvature,
132  const GEnergy& pivot) :
133  m_norm(norm),
134  m_index(index),
135  m_curvature(curvature),
136  m_pivot(pivot) {};
137 
138  // Method
139  double eval(const double& x) {
140  double xrel = x/m_pivot.MeV();
141  return m_norm*std::pow(xrel, m_index + m_curvature * std::log(xrel));
142  }
143 
144  // Members
145  protected:
146  double m_norm; //!< Normalization
147  double m_index; //!< Spectral index at pivot
148  double m_curvature; //!< Curvature
149  GEnergy m_pivot; //!< Pivot energy
150  };
151 
152  // Class to determine the integral energy flux, derviation of flux_kern
153  class eflux_kern : public flux_kern {
154  public:
155  // Constructor
156  eflux_kern(const double& norm,
157  const double& index,
158  const double& curvature,
159  const GEnergy& pivot) :
160  flux_kern(norm, index, curvature, pivot) {};
161 
162  // Method
163  double eval(const double& x) {
164  return x * flux_kern::eval(x);
165  }
166  };
167 
168 
169  // Protected members
170  std::string m_type; //!< Model type
171  GModelPar m_norm; //!< Normalization factor
172  GModelPar m_index; //!< Spectral index
173  GModelPar m_curvature; //!< Curvature
174  GModelPar m_pivot; //!< Pivot energy
175 
176  // Cached members used for pre-computations
177  mutable GEnergy m_last_energy; //!< Last energy value
178  mutable double m_last_index; //!< Last index parameter
179  mutable double m_last_curvature; //!< Last curvature parameters
180  mutable double m_last_pivot; //!< Last pivot parameter
181  mutable double m_last_e_norm; //!< Last E/Epivot value
182  mutable double m_last_log_e_norm; //!< Last ln(E/Epivot) value
183  mutable double m_last_exponent; //!< Last exponent
184  mutable double m_last_power; //!< Last power value
185 
186  // Cached members used for pre-computations
187  mutable double m_mc_emin; //!< Minimum energy
188  mutable double m_mc_emax; //!< Maximum energy
189  mutable double m_mc_exponent; //!< Exponent (index+1)
190  mutable double m_mc_pow_emin; //!< Power of minimum energy
191  mutable double m_mc_pow_ewidth; //!< Power of energy width
192  mutable double m_mc_norm; //!< Norm of powerlaw model at logparabola pivot energy
193 };
194 
195 
196 /***********************************************************************//**
197  * @brief Return class name
198  *
199  * @return String containing the class name ("GModelSpectralLogParabola").
200  ***************************************************************************/
201 inline
203 {
204  return ("GModelSpectralLogParabola");
205 }
206 
207 
208 /***********************************************************************//**
209  * @brief Return model type
210  *
211  * @return Model type.
212  *
213  * Returns the type of the log parabola spectral model.
214  ***************************************************************************/
215 inline
216 std::string GModelSpectralLogParabola::type(void) const
217 {
218  return (m_type);
219 }
220 
221 
222 /***********************************************************************//**
223  * @brief Set model type
224  *
225  * @param[in] type Model type.
226  *
227  * Set the type of the log parabola spectral model.
228  ***************************************************************************/
229 inline
230 void GModelSpectralLogParabola::type(const std::string& type)
231 {
232  m_type = type;
233  return;
234 }
235 
236 
237 /***********************************************************************//**
238  * @brief Return pre factor
239  *
240  * @return Pre factor (ph/cm2/s/MeV).
241  *
242  * Returns the pre factor.
243  ***************************************************************************/
244 inline
246 {
247  return (m_norm.value());
248 }
249 
250 
251 /***********************************************************************//**
252  * @brief Set pre factor
253  *
254  * @param[in] prefactor Pre factor (ph/cm2/s/MeV).
255  *
256  * Sets the pre factor.
257  ***************************************************************************/
258 inline
259 void GModelSpectralLogParabola::prefactor(const double& prefactor)
260 {
261  m_norm.value(prefactor);
262  return;
263 }
264 
265 
266 /***********************************************************************//**
267  * @brief Return spectral index
268  *
269  * @return Spectral index.
270  *
271  * Returns the spectral index.
272  ***************************************************************************/
273 inline
275 {
276  return (m_index.value());
277 }
278 
279 
280 /***********************************************************************//**
281  * @brief Set spectral index
282  *
283  * @param[in] index Spectral index.
284  *
285  * Sets the spectral index.
286  ***************************************************************************/
287 inline
288 void GModelSpectralLogParabola::index(const double& index)
289 {
290  m_index.value(index);
291  return;
292 }
293 
294 
295 /***********************************************************************//**
296  * @brief Return spectral curvature
297  *
298  * @return Spectral curvature.
299  *
300  * Returns the spectral curvature.
301  ***************************************************************************/
302 inline
304 {
305  return (m_curvature.value());
306 }
307 
308 
309 /***********************************************************************//**
310  * @brief Set spectral curvature
311  *
312  * @param[in] curvature Spectral curvature.
313  *
314  * Sets the spectral curvature.
315  ***************************************************************************/
316 inline
317 void GModelSpectralLogParabola::curvature(const double& curvature)
318 {
319  m_curvature.value(curvature);
320  return;
321 }
322 
323 
324 /***********************************************************************//**
325  * @brief Return pivot energy
326  *
327  * @return Pivot energy.
328  *
329  * Returns the pivot energy.
330  ***************************************************************************/
331 inline
333 {
334  GEnergy energy;
335  energy.MeV(m_pivot.value());
336  return energy;
337 }
338 
339 
340 /***********************************************************************//**
341  * @brief Set pivot energy
342  *
343  * @param[in] pivot Pivot energy.
344  *
345  * Sets the pivot energy.
346  ***************************************************************************/
347 inline
349 {
350  m_pivot.value(pivot.MeV());
351  return;
352 }
353 
354 #endif /* GMODELSPECTRALLOGPARABOLA_HPP */
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:821
Energy value class definition.
Abstract spectral model base class.
GModelSpectralLogParabola(void)
Void constructor.
double m_mc_exponent
Exponent (index+1)
GModelPar m_index
Spectral index.
virtual void clear(void)
Clear log parabola model.
void copy_members(const GModelSpectralLogParabola &model)
Copy class members.
virtual GModelSpectralLogParabola * clone(void) const
Clone log parabola model.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print LogParabola information.
XML element node class.
Definition: GXmlElement.hpp:47
void update_eval_cache(const GEnergy &energy) const
Update eval precomputation cache.
GEnergy pivot(void) const
Return pivot energy.
Random number generator class.
Definition: GRan.hpp:44
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
Time class.
Definition: GTime.hpp:54
LogParabola spectral model class.
double m_last_log_e_norm
Last ln(E/Epivot) value.
double m_last_curvature
Last curvature parameters.
virtual std::string classname(void) const
Return class name.
virtual std::string type(void) const
Return model type.
Model parameter class interface definition.
double prefactor(void) const
Return pre factor.
Model parameter class.
Definition: GModelPar.hpp:87
Single parameter function abstract base class definition.
double curvature(void) const
Return spectral curvature.
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1184
virtual double flux(const GEnergy &emin, const GEnergy &emax) const
Returns model photon flux between emin, emax
void init_members(void)
Initialise class members.
GModelPar m_norm
Normalization factor.
eflux_kern(const double &norm, const double &index, const double &curvature, const GEnergy &pivot)
virtual void read(const GXmlElement &xml)
Read model from XML element.
GChatter
Definition: GTypemaps.hpp:33
void free_members(void)
Delete class members.
double m_last_pivot
Last pivot parameter.
virtual GEnergy mc(const GEnergy &emin, const GEnergy &emax, const GTime &time, GRan &ran) const
Returns Monte Carlo energy between [emin, emax].
double m_mc_pow_emin
Power of minimum energy.
Abstract spectral model base class interface definition.
double index(void) const
Return spectral index.
virtual double eval(const GEnergy &srcEng, const GTime &srcTime=GTime(), const bool &gradients=false) const
Evaluate model.
double m_last_e_norm
Last E/Epivot value.
virtual void write(GXmlElement &xml) const
Write model into XML element.
virtual double eflux(const GEnergy &emin, const GEnergy &emax) const
Returns model energy flux between emin, emax
GEnergy m_last_energy
Last energy value.
double m_last_index
Last index parameter.
double value(void) const
Return parameter value.
Single parameter function abstract base class.
Definition: GFunction.hpp:44
flux_kern(const double &norm, const double &index, const double &curvature, const GEnergy &pivot)
GVector pow(const GVector &vector, const double &power)
Computes tanh of vector elements.
Definition: GVector.cpp:1332
double m_mc_pow_ewidth
Power of energy width.
virtual ~GModelSpectralLogParabola(void)
Destructor.
void update_mc_cache(const GEnergy &emin, const GEnergy &emax, const GTime &time) const
Update Monte Carlo pre computation cache.
virtual GModelSpectralLogParabola & operator=(const GModelSpectralLogParabola &model)
Assignment operator.
double m_last_power
Last power value.
double m_mc_norm
Norm of powerlaw model at logparabola pivot energy.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48