73 static const double sqrtTwo = sqrt(2.0);
80 std::vector<double> punctures(_punctures.size());
81 for (
size_t i=0;i<_punctures.size();i++) punctures[i]=_punctures[i].getValue();
90 for (
size_t i=0;i<punctures.size()/2;i++) {
91 std::sort(punctures.begin()+2*i, punctures.begin()+2*i+2);
92 double min1=punctures[2*i];
93 double max1=punctures[2*i+1];
94 for (
size_t j=i+1;j<punctures.size()/2;j++) {
95 std::sort(punctures.begin()+2*j, punctures.begin()+2*j+2);
96 double min2=punctures[2*j];
97 double max2=punctures[2*j+1];
98 if ((min2>min1 && max1>min2) || (min1>min2 && max2<min1)) {
99 punctures[2*i] =std::min(min1,min2);
100 punctures[2*i+1]=std::max(max1,max2);
101 std::vector<double>::iterator t0=punctures.begin()+2*j, t1=t0+2;
102 punctures.erase(t0,t1);
111 double expG=0,norm=0;
112 for (
size_t i=0;i<punctures.size()/2;i++) {
114 double a = punctures[2*i];
115 double b = punctures[2*i+1];
117 double alpha = (
a/xsigma + xsigma/tau)/sqrtTwo;
118 double beta = (
b/xsigma + xsigma/tau)/sqrtTwo;
119 double delta = 1/sqrtTwo/xsigma;
122 norm += 2*tau*exp(1/(4*delta*delta*tau*tau))*(exp(-alpha/(delta*tau)) - exp(-beta/(delta*tau)));
124 expG += ((erfc(alpha-delta*x) - erfc(beta-delta*x))*exp(-x/tau) );
130 if (norm==0)
return norm;