GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelSpectralBrokenPlaw.hpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelSpectralBrokenPlaw.hpp - Broken power law spectrum class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2013-2018 by Anneli Schulz *
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 GModelSpectralBrokenPlaw.hpp
23  * @brief Broken power law spectrum class definition
24  * @author Anneli Schulz
25  */
26 
27 #ifndef GMODELSPECTRALBROKENPLAW_HPP
28 #define GMODELSPECTRALBROKENPLAW_HPP
29 
30 /* __ Includes ___________________________________________________________ */
31 #include <string>
32 #include "GModelSpectral.hpp"
33 #include "GModelPar.hpp"
34 #include "GEnergy.hpp"
35 
36 /* __ Forward declarations _______________________________________________ */
37 class GRan;
38 class GTime;
39 class GXmlElement;
40 
41 
42 /***********************************************************************//**
43  * @class GModelSpectralBrokenPlaw
44  *
45  * @brief Broken power law spectral model class
46  *
47  * This class implements a broken power law spectrum. The model is defined
48  * by
49  *
50  * \f[
51  * S_{\rm E}(E | t) = k_0 \times \left \{
52  * \begin{eqnarray}
53  * \left( \frac{E}{E_b} \right)^{\gamma_1} & {\rm if\,\,} E < E_b \\
54  * \left( \frac{E}{E_b} \right)^{\gamma_2} & {\rm otherwise}
55  * \end{eqnarray}
56  * \right .
57  * \f]
58  *
59  * where
60  * \f$k_0\f$ is the normalization or prefactor,
61  * \f$\gamma_1\f$ is the spectral index before the break,
62  * \f$\gamma_2\f$ is the spectral index after the break, and
63  * \f$E_b\f$ is the break energy.
64  ***************************************************************************/
66 
67 public:
68  // Constructors and destructors
70  GModelSpectralBrokenPlaw(const std::string& type,
71  const std::string& prefactor,
72  const std::string& index1,
73  const std::string& breakenergy,
74  const std::string& index2);
75  GModelSpectralBrokenPlaw(const double& prefactor,
76  const double& index1,
77  const GEnergy& breakenergy,
78  const double& index2);
79  explicit GModelSpectralBrokenPlaw(const GXmlElement& xml);
81  virtual ~GModelSpectralBrokenPlaw(void);
82 
83  // Operators
85 
86  // Implemented pure virtual methods
87  virtual void clear(void);
88  virtual GModelSpectralBrokenPlaw* clone(void) const;
89  virtual std::string classname(void) const;
90  virtual std::string type(void) const;
91  virtual double eval(const GEnergy& srcEng,
92  const GTime& srcTime = GTime(),
93  const bool& gradients = false) const;
94  virtual double flux(const GEnergy& emin,
95  const GEnergy& emax) const;
96  virtual double eflux(const GEnergy& emin,
97  const GEnergy& emax) const;
98  virtual GEnergy mc(const GEnergy& emin,
99  const GEnergy& emax,
100  const GTime& time,
101  GRan& ran) const;
102  virtual void read(const GXmlElement& xml);
103  virtual void write(GXmlElement& xml) const;
104  virtual std::string print(const GChatter& chatter = NORMAL) const;
105 
106  // Other methods
107  void type(const std::string& type);
108  double prefactor(void) const;
109  double index1(void) const;
110  double index2(void) const;
111  GEnergy breakenergy(void) const;
112  void prefactor(const double& prefactor);
113  void index1(const double& index1);
114  void index2(const double& index2);
115  void breakenergy(const GEnergy& breakenergy);
116 
117 protected:
118  // Protected methods
119  void init_members(void);
120  void copy_members(const GModelSpectralBrokenPlaw& model);
121  void free_members(void);
122  void update_eval_cache(const GEnergy& energy) const;
123  void update_mc_cache(const GEnergy& emin, const GEnergy& emax) const;
124 
125  // Protected members
126  std::string m_type; //!< Model type
127  GModelPar m_norm; //!< Normalization factor
128  GModelPar m_index1; //!< Spectral index1
129  GModelPar m_index2; //!< Spectral index2
130  GModelPar m_breakenergy; //!< Energy of spectral break
131 
132  // Cached members used for pre-computations
133  mutable GEnergy m_last_energy; //!< Last energy value
134  mutable double m_last_index1; //!< Last index1 parameter
135  mutable double m_last_index2; //!< Last index1 parameter
136  mutable double m_last_breakenergy; //!< Last breakenergy parameter
137  mutable double m_last_e_norm; //!< Last E/Ebreakenergy value
138  mutable double m_last_log_e_norm; //!< Last ln(E/Ebreakenergy) value
139  mutable double m_last_power; //!< Last power value
140  mutable double m_mc_emin; //!< Minimum energy
141  mutable double m_mc_emax; //!< Maximum energy
142  mutable double m_mc_exponent1; //!< Exponent (index1+1)
143  mutable double m_mc_exponent2; //!< Exponent (index2+1)
144  mutable double m_mc_pow_emin; //!< Power of minimum energy
145  mutable double m_mc_pow_ewidth; //!< Power of energy width
146  mutable std::vector<double> m_mc_cum; //!< Cumulative distribution
147  mutable std::vector<double> m_mc_min; //!< Lower boundary for MC
148  mutable std::vector<double> m_mc_max; //!< Upper boundary for MC
149  mutable std::vector<double> m_mc_exp; //!< Exponent for MC
150 };
151 
152 
153 /***********************************************************************//**
154  * @brief Return class name
155  *
156  * @return String containing the class name ("GModelSpectralBrokenPlaw").
157  ***************************************************************************/
158 inline
159 std::string GModelSpectralBrokenPlaw::classname(void) const
160 {
161  return ("GModelSpectralBrokenPlaw");
162 }
163 
164 
165 /***********************************************************************//**
166  * @brief Return model type
167  *
168  * @return "PowerLaw".
169  *
170  * Returns the type of the spectral broken power law model.
171  ***************************************************************************/
172 inline
173 std::string GModelSpectralBrokenPlaw::type(void) const
174 {
175  return (m_type);
176 }
177 
178 
179 /***********************************************************************//**
180  * @brief Set model type
181  *
182  * @param[in] type Model type.
183  *
184  * Set the type of the spectral broken power law model.
185  ***************************************************************************/
186 inline
187 void GModelSpectralBrokenPlaw::type(const std::string& type)
188 {
189  m_type = type;
190  return;
191 }
192 
193 
194 /***********************************************************************//**
195  * @brief Return pre factor
196  *
197  * @return Pre factor (ph/cm2/s/MeV).
198  *
199  * Returns the pre factor.
200  ***************************************************************************/
201 inline
203 {
204  return (m_norm.value());
205 }
206 
207 
208 /***********************************************************************//**
209  * @brief Set pre factor
210  *
211  * @param[in] prefactor Pre factor (ph/cm2/s/MeV).
212  *
213  * Sets the pre factor.
214  ***************************************************************************/
215 inline
216 void GModelSpectralBrokenPlaw::prefactor(const double& prefactor)
217 {
218  m_norm.value(prefactor);
219  return;
220 }
221 
222 
223 /***********************************************************************//**
224  * @brief Return power law index1
225  *
226  * @return Power law index1.
227  *
228  * Returns the power law index1.
229  ***************************************************************************/
230 inline
232 {
233  return (m_index1.value());
234 }
235 
236 
237 /***********************************************************************//**
238  * @brief Set power law index1
239  *
240  * @param[in] index1 Power law index1.
241  *
242  * Sets the power law index1.
243  ***************************************************************************/
244 inline
245 void GModelSpectralBrokenPlaw::index1(const double& index1)
246 {
247  m_index1.value(index1);
248  return;
249 }
250 
251 
252 /***********************************************************************//**
253  * @brief Return power law index2
254  *
255  * @return Power law index2.
256  *
257  * Returns the power law index2.
258  ***************************************************************************/
259 inline
261 {
262  return (m_index2.value());
263 }
264 
265 
266 /***********************************************************************//**
267  * @brief Set power law index2
268  *
269  * @param[in] index2 Power law index2.
270  *
271  * Sets the power law index2.
272  ***************************************************************************/
273 inline
274 void GModelSpectralBrokenPlaw::index2(const double& index2)
275 {
276  m_index2.value(index2);
277  return;
278 }
279 
280 
281 /***********************************************************************//**
282  * @brief Return breakenergy energy
283  *
284  * @return breakenergy energy.
285  *
286  * Returns the breakenergy energy.
287  ***************************************************************************/
288 inline
290 {
291  GEnergy energy;
292  energy.MeV(m_breakenergy.value());
293  return energy;
294 }
295 
296 
297 /***********************************************************************//**
298  * @brief Set breakenergy energy
299  *
300  * @param[in] breakenergy breakenergy energy.
301  *
302  * Sets the breakenergy energy.
303  ***************************************************************************/
304 inline
306 {
307  m_breakenergy.value(breakenergy.MeV());
308  return;
309 }
310 
311 #endif /* GMODELSPECTRALPLAW_HPP */
GModelPar m_index2
Spectral index2.
void copy_members(const GModelSpectralBrokenPlaw &model)
Copy class members.
double m_last_index2
Last index1 parameter.
std::vector< double > m_mc_min
Lower boundary for MC.
Energy value class definition.
double index2(void) const
Return power law index2.
void update_mc_cache(const GEnergy &emin, const GEnergy &emax) const
Update Monte Carlo pre computation cache.
Abstract spectral model base class.
virtual void read(const GXmlElement &xml)
Read model from XML element.
virtual std::string type(void) const
Return model type.
XML element node class.
Definition: GXmlElement.hpp:48
void free_members(void)
Delete class members.
Random number generator class.
Definition: GRan.hpp:44
double m_last_log_e_norm
Last ln(E/Ebreakenergy) value.
std::vector< double > m_mc_exp
Exponent for MC.
double prefactor(void) const
Return pre factor.
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
GModelPar m_index1
Spectral index1.
Time class.
Definition: GTime.hpp:55
double m_last_power
Last power value.
virtual double eval(const GEnergy &srcEng, const GTime &srcTime=GTime(), const bool &gradients=false) const
Evaluate function.
double m_last_index1
Last index1 parameter.
virtual GModelSpectralBrokenPlaw * clone(void) const
Clone broken power law model.
void update_eval_cache(const GEnergy &energy) const
Update eval precomputation cache.
virtual void clear(void)
Clear broken power law model.
double m_mc_exponent1
Exponent (index1+1)
std::vector< double > m_mc_cum
Cumulative distribution.
Model parameter class interface definition.
GModelPar m_norm
Normalization factor.
virtual void write(GXmlElement &xml) const
Write model into XML element.
Model parameter class.
Definition: GModelPar.hpp:87
GEnergy breakenergy(void) const
Return breakenergy energy.
GModelSpectralBrokenPlaw(void)
Void constructor.
std::vector< double > m_mc_max
Upper boundary for MC.
GModelPar m_breakenergy
Energy of spectral break.
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.
GChatter
Definition: GTypemaps.hpp:33
double index1(void) const
Return power law index1.
Broken power law spectral model class.
virtual double flux(const GEnergy &emin, const GEnergy &emax) const
Returns model photon flux between emin, emax
Abstract spectral model base class interface definition.
double m_last_breakenergy
Last breakenergy parameter.
virtual std::string classname(void) const
Return class name.
virtual ~GModelSpectralBrokenPlaw(void)
Destructor.
virtual GModelSpectralBrokenPlaw & operator=(const GModelSpectralBrokenPlaw &model)
Assignment operator.
double m_mc_pow_ewidth
Power of energy width.
double value(void) const
Return parameter value.
double m_mc_exponent2
Exponent (index2+1)
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_e_norm
Last E/Ebreakenergy value.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print brokenpowerlaw information.
void init_members(void)
Initialise class members.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48