00001
00013 #ifdef _MSC_VER
00014 #include "msdevstudio/MSconfig.h"
00015 #endif
00016
00017 #include "Erfc.h"
00018 #include "FunctionHelper.h"
00019
00020 #include <cmath>
00021
00022 #include <cassert>
00023
00024 #ifdef ITERATOR_MEMBER_DEFECT
00025 using namespace std;
00026 #else
00027 using std::exp;
00028 using std::vector;
00029 #endif
00030
00031 namespace hippodraw {
00032
00033 Erfc::Erfc ( )
00034 {
00035 initialize ();
00036 }
00037
00038 Erfc::Erfc(double m, double s, double n)
00039 {
00040 initialize();
00041
00042 m_parms[MEAN] = m;
00043 m_parms[SIGMA] = s;
00044 m_parms[NORM] = n;
00045 }
00046
00047 void Erfc::initialize ()
00048 {
00049 m_name = "Erfc";
00050
00051 m_parm_names.push_back ( "Mean" );
00052 m_parm_names.push_back ( "Sigma" );
00053 m_parm_names.push_back ( "Normalization" );
00054
00055 resize ();
00056 }
00057
00058 FunctionBase * Erfc::clone () const
00059 {
00060 return new Erfc ( *this );
00061 }
00062
00065 double Erfc::operator () (double value) const
00066 {
00067 double result = 1;
00068
00069 double x = calcRed(value);
00070
00071 result = calcErfc(x);
00072
00073 result *= m_parms[NORM];
00074
00075 return result;
00076 }
00077
00078
00079
00080 void Erfc::initialParameters ( const FunctionHelper * helper )
00081 {
00082 m_parms[MEAN] = helper->meanCoord();
00083 m_parms[SIGMA] = 1;
00084 m_parms[NORM] = helper->maxValue () / 2;
00085 }
00086
00087
00088 double Erfc::derivByParm ( int i, double value ) const
00089 {
00090 double result = 0;
00091
00092 switch ( i ) {
00093 case MEAN :
00094 result = -1.0/m_parms[SIGMA] * m_parms[NORM] * derivByRed ( calcRed(value) );
00095 return result;
00096 break;
00097
00098 case SIGMA :
00099 result = - calcRed(value)/m_parms[SIGMA] * m_parms[NORM] * derivByRed ( calcRed(value) );
00100 return result;
00101 break;
00102
00103 case NORM :
00104 result = operator () ( value ) / m_parms[NORM];
00105 return result;
00106 break;
00107
00108 default :
00109 assert ( false );
00110 break;
00111 }
00112
00113 return 0.0;
00114 }
00115
00116
00117
00118 double Erfc::derivByRed( const double x ) const
00119 {
00120
00121 static const double m_2_sqrtpi = 1.12837916709551257390;
00122
00123 double result = - m_2_sqrtpi;
00124
00125 result *= exp(-x*x);
00126 return result;
00127 }
00128
00129
00130
00131 double Erfc::calcErfc(double x) const
00132 {
00133
00134
00135
00136
00137
00138
00139 const double
00140 a1 = -1.26551223, a2 = 1.00002368,
00141 a3 = 0.37409196, a4 = 0.09678418,
00142 a5 = -0.18628806, a6 = 0.27886807,
00143 a7 = -1.13520398, a8 = 1.48851587,
00144 a9 = -0.82215223, a10 = 0.17087277;
00145
00146 double v = 1;
00147
00148 double z = x>=0 ? x : -x;
00149
00150 if (z <= 0) return v;
00151
00152 double t = 1/(1+0.5*z);
00153
00154 v = t*exp((-z*z) +a1+t*(a2+t*(a3+t*(a4+t*(a5+t*(a6+t*(a7+t*(a8+t*(a9+t*a10)))))))));
00155
00156 if (x < 0) v = 2-v;
00157
00158 return v;
00159 }
00160
00161 }