BrokenPowerLaw.cxx
Go to the documentation of this file.
1 
12 #ifdef _MSC_VER
13 #include "msdevstudio/MSconfig.h"
14 #endif
15 
16 #include "BrokenPowerLaw.h"
17 
18 #include "FunctionHelper.h"
19 
20 #include <cmath>
21 #include <cassert>
22 
23 #ifdef ITERATOR_MEMBER_DEFECT
24 using namespace std;
25 #else
26 using std::max;
27 using std::vector;
28 #endif
29 
30 namespace hippodraw {
31 
32 BrokenPowerLaw::BrokenPowerLaw ( )
33 {
34  initialize ();
35 }
36 
37 BrokenPowerLaw::BrokenPowerLaw ( double prefactor, double index1,
38  double index2, double x_break)
39 {
40  initialize ();
41 
42  m_parms[0] = prefactor;
43  m_parms[1] = index1;
44  m_parms[2] = index2;
45  m_parms[3] = x_break;
46 }
47 
48 void BrokenPowerLaw::initialize ()
49 {
50  m_name = "BrokenPowerLaw";
51  m_parm_names.push_back ( "Prefactor" );
52  m_parm_names.push_back ( "Index1" );
53  m_parm_names.push_back ( "Index2" );
54  m_parm_names.push_back ( "Break" );
55 
56  resize ();
57 }
58 
60 {
61  return new BrokenPowerLaw ( *this );
62 }
63 
64 double BrokenPowerLaw::operator () ( double x ) const
65 {
66  if (x < m_parms[3]) {
67  return m_parms[0]*pow(x/m_parms[3], m_parms[1]);
68  } else {
69  return m_parms[0]*pow(x/m_parms[3], m_parms[2]);
70  }
71 }
72 
73 /* virtual */
74 void
75 BrokenPowerLaw::
76 initialParameters ( const FunctionHelper * helper )
77 {
78  double min_x = helper->minCoord ();
79  double max_x = helper->maxCoord ();
80  max_x = (min_x + max_x)/2.;
81 
82  double min_y, max_y;
83  try {
84  min_y = helper->valueAt(min_x);
85  max_y = helper->valueAt(max_x);
86  if (min_y != 0 && max_y != 0) {
87  m_parms[1] = log( max_y/min_y ) / log( max_x/min_x );
88  m_parms[2] = m_parms[1];
89  m_parms[3] = helper->meanCoord();
90  m_parms[0] = max_y/pow(max_x/m_parms[3], m_parms[1]);
91  return;
92  }
93  } catch (...) {
94 // do nothing
95  }
96 
97 // default behavior
98  min_y = max(helper->minValue(), 1.);
99  max_y = helper->maxValue();
100  m_parms[1] = log( max_y/min_y ) / log( max_x/min_x );
101  m_parms[2] = m_parms[1];
102  m_parms[3] = helper->meanCoord();
103  m_parms[0] = max_y/pow(max_x/m_parms[3], m_parms[1]);
104 }
105 
106 double BrokenPowerLaw::derivByParm ( int i, double x ) const
107 {
108  switch ( i ) {
109  case 0 :
110  return operator()(x)/m_parms[0];
111  break;
112 
113  case 1 :
114  if (x < m_parms[3]) {
115  return operator()(x)*log(x/m_parms[3]);
116  } else {
117  return 0;
118  }
119  break;
120 
121  case 2 :
122  if (x < m_parms[3]) {
123  return 0;
124  } else {
125  return operator()(x)*log(x/m_parms[3]);
126  }
127  break;
128 
129  case 3 :
130  if (x < m_parms[3]) {
131  return -m_parms[1]*operator()(x)/m_parms[3];
132  } else {
133  return -m_parms[2]*operator()(x)/m_parms[3];
134  }
135  break;
136 
137  default:
138  assert ( false );
139  break;
140  }
141  return 0.0;
142 }
143 
144 } // namespace hippodraw

Generated for HippoDraw Class Library by doxygen