66 for(
unsigned int i = 0; i < point1.size(); i++ ) {
73 if( point2[i] < point1[i] ) {
83 if( std::abs(y2 - y1) < 1E-7 )
continue;
86 double beta1 = 1 + 2 * (y1 -
m_lower( i )) / (y2 - y1);
87 double beta2 = 1 + 2 * (
m_upper( i ) - y2) / (y2 - y1);
88 double expp =
m_nc + 1.;
90 double alpha1 = 2. - std::pow(beta1 , -expp);
91 double alpha2 = 2. - std::pow(beta2 , -expp);
96 if( u > 1. / alpha1 ) {
97 alpha1 = 1. / (2. - alpha1);
99 if( u > 1. / alpha2 ) {
100 alpha2 = 1. / (2. - alpha2);
102 betaQ1 = std::pow( alpha1, 1.0/expp );
103 betaQ2 = std::pow( alpha2, 1.0/expp );
106 point1[i] = 0.5 * ((y1 + y2) - betaQ1 * (y2 - y1));
107 point2[i] = 0.5 * ((y1 + y2) + betaQ2 * (y2 - y1));
113 point1[i] = std::max( point1[i],
m_lower( i ) );
114 point1[i] = std::min( point1[i],
m_upper( i ) );
115 point2[i] = std::max( point2[i],
m_lower( i ) );
116 point2[i] = std::min( point2[i],
m_upper( i ) );