BrokenPowerLaw.cxx

Go to the documentation of this file.
00001 
00012 #ifdef _MSC_VER
00013 #include "msdevstudio/MSconfig.h"
00014 #endif
00015 
00016 #include "BrokenPowerLaw.h"
00017 
00018 #include "FunctionHelper.h"
00019 
00020 #include <cmath>
00021 #include <cassert>
00022 
00023 #ifdef ITERATOR_MEMBER_DEFECT
00024 using namespace std;
00025 #else
00026 using std::max;
00027 using std::vector;
00028 #endif
00029 
00030 namespace hippodraw {
00031 
00032 BrokenPowerLaw::BrokenPowerLaw ( )
00033 {
00034   initialize ();
00035 }
00036 
00037 BrokenPowerLaw::BrokenPowerLaw ( double prefactor, double index1,
00038                                  double index2, double x_break)
00039 {
00040   initialize ();
00041 
00042   m_parms[0] = prefactor;
00043   m_parms[1] = index1;
00044   m_parms[2] = index2;
00045   m_parms[3] = x_break;
00046 }
00047 
00048 void BrokenPowerLaw::initialize ()
00049 {
00050   m_name = "BrokenPowerLaw";
00051   m_parm_names.push_back ( "Prefactor" );
00052   m_parm_names.push_back ( "Index1" );
00053   m_parm_names.push_back ( "Index2" );
00054   m_parm_names.push_back ( "Break" );
00055 
00056   resize ();
00057 }
00058 
00059 FunctionBase * BrokenPowerLaw::clone () const
00060 {
00061   return new BrokenPowerLaw ( *this );
00062 }
00063 
00064 double BrokenPowerLaw::operator () ( double x ) const
00065 {
00066    if (x < m_parms[3]) {
00067       return m_parms[0]*pow(x/m_parms[3], m_parms[1]);
00068    } else {
00069       return m_parms[0]*pow(x/m_parms[3], m_parms[2]);
00070    }
00071 }
00072 
00073 /* virtual */
00074 void 
00075 BrokenPowerLaw::
00076 initialParameters ( const FunctionHelper * helper )
00077 {
00078   double min_x = helper->minCoord ();
00079   double max_x = helper->maxCoord ();
00080   max_x = (min_x + max_x)/2.;
00081 
00082   double min_y, max_y;
00083   try {
00084      min_y = helper->valueAt(min_x);
00085      max_y = helper->valueAt(max_x);
00086      if (min_y != 0 && max_y != 0) {
00087         m_parms[1] = log( max_y/min_y ) / log( max_x/min_x );
00088         m_parms[2] = m_parms[1];
00089         m_parms[3] = helper->meanCoord();
00090         m_parms[0] = max_y/pow(max_x/m_parms[3], m_parms[1]);
00091         return;
00092      }
00093   } catch (...) {
00094 // do nothing
00095   }
00096 
00097 // default behavior
00098   min_y = max(helper->minValue(), 1.);
00099   max_y = helper->maxValue();
00100   m_parms[1] = log( max_y/min_y ) / log( max_x/min_x );
00101   m_parms[2] = m_parms[1];
00102   m_parms[3] = helper->meanCoord();
00103   m_parms[0] = max_y/pow(max_x/m_parms[3], m_parms[1]);
00104 }
00105 
00106 double BrokenPowerLaw::derivByParm ( int i, double x ) const
00107 {
00108   switch ( i ) {
00109   case 0 :
00110      return operator()(x)/m_parms[0];
00111      break;
00112 
00113   case 1 :
00114      if (x < m_parms[3]) {
00115         return operator()(x)*log(x/m_parms[3]);
00116      } else {
00117         return 0;
00118      }
00119      break;
00120 
00121   case 2 :
00122      if (x < m_parms[3]) {
00123         return 0;
00124      } else {
00125         return operator()(x)*log(x/m_parms[3]);
00126      }
00127      break;
00128 
00129   case 3 :
00130      if (x < m_parms[3]) {
00131         return -m_parms[1]*operator()(x)/m_parms[3];
00132      } else {
00133         return -m_parms[2]*operator()(x)/m_parms[3];
00134      }
00135      break;
00136 
00137   default:
00138     assert ( false );
00139     break;
00140   }
00141   return 0.0;
00142 }
00143 
00144 } // namespace hippodraw

Generated for HippoDraw Class Library by doxygen