@@ -131,23 +131,20 @@ namespace Splines {
131131
132132 integer const n{ npts - 1 };
133133
134- // integer ierr = 0;
135-
136134 // function definition is ok, go on.
137- real_type h1{ X[1 ] - X[0 ] };
138- real_type del1{ (Y[1 ]-Y[0 ])/h1 };
139- // real_type dsave = del1;
135+ real_type h1 { X[1 ] - X[0 ] };
136+ real_type del1 { (Y[1 ]-Y[0 ])/h1 };
140137
141138 // special case n=2 -- use linear interpolation.
142139 if ( n == 1 ) { Yp[0 ] = Yp[1 ] = del1; return ; }
143140
144- real_type h2{ X[2 ] - X[1 ] };
145- real_type del2{ (Y[2 ]-Y[1 ])/h2 };
141+ real_type h2 { X[2 ] - X[1 ] };
142+ real_type del2 { (Y[2 ]-Y[1 ])/h2 };
146143
147144 // Set Yp[0] via non-centered three-point formula, adjusted to be shape-preserving.
148- real_type hsum = h1 + h2;
149- real_type w1 = (h1 + hsum) /hsum;
150- real_type w2 = -h1/hsum;
145+ real_type hsum { h1 + h2 } ;
146+ real_type w1 { 1 +h1 /hsum } ;
147+ real_type w2 { -h1/hsum } ;
151148 Yp[0 ] = w1*del1 + w2*del2;
152149 real_type dmin, dmax;
153150 if ( signTest (Yp[0 ],del1) <= 0 ) {
@@ -173,20 +170,16 @@ namespace Splines {
173170 switch ( signTest (del1,del2) ) {
174171 case -1 :
175172 if ( Utils::is_zero (del2) ) break ;
176- // if ( signTest(dsave,del2) < 0 ) ++ierr;
177- // dsave = del2;
178173 break ;
179174 case 0 :
180- // ++ierr;
181- // dsave = del2;
182175 break ;
183176 case 1 : // use brodlie modification of butland formula.
184177 w1 = (1 +h1/hsum)/3 ;
185178 w2 = (1 +h2/hsum)/3 ;
186179 dmax = max_abs ( del1, del2 );
187180 dmin = min_abs ( del1, del2 );
188- real_type const drat1 = del1/dmax;
189- real_type const drat2 = del2/dmax;
181+ real_type const drat1{ del1/dmax } ;
182+ real_type const drat2{ del2/dmax } ;
190183 Yp[i] = dmin/(w1*drat1 + w2*drat2);
191184 break ;
192185 }
@@ -205,73 +198,23 @@ namespace Splines {
205198 // cout << "ierr = " << ierr << '\n';
206199 }
207200
208- #if 0
209- static // non usata per ora
210- void
211- Pchip_build_new(
212- real_type const X[],
213- real_type const Y[],
214- real_type Yp[],
215- integer npts
216- ) {
217-
218- first_derivative_build( X, Y, Yp, npts );
219-
220- integer n = npts > 0 ? npts - 1 : 0;
221-
222- // loop through interior points.
223- for ( integer i{1}; i < n ; ++i ) {
224-
225- real_type hL = X[i+0] - X[i-1];
226- real_type hR = X[i+1] - X[i+0];
227-
228- real_type SL = (Y[i+0]-Y[i-1])/hL;
229- real_type SR = (Y[i+1]-Y[i+0])/hR;
230-
231- real_type fp = Yp[i];
232-
233- real_type sigma = 0;
234- if ( SL*SR > 0 ) sigma = SR > 0 ? 1 : -1;
235- real_type absSL = SL < 0 ? -SL: SL;
236- real_type absSR = SR < 0 ? -SR: SR;
237- real_type Delta = 3*( absSL < absSR ? absSL : absSR );
238- if ( sigma > 0 ) {
239- if ( fp < 0 ) fp = 0;
240- if ( fp > Delta ) fp = Delta;
241- } else if ( sigma < 0 ) {
242- if ( fp > 0 ) fp = 0;
243- if ( fp < -Delta ) fp = -Delta;
244- } else {
245- fp = 0;
246- }
247- Yp[i] = fp;
248- }
249- }
250- #endif
251-
252201 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
253202
254203 void
255204 PchipSpline::build () {
256205 string msg{ fmt::format (" PchipSpline[{}]::build():" , m_name ) };
257- UTILS_ASSERT (
258- m_npts > 1 ,
259- " {} npts = {} not enought points\n " ,
260- msg, m_npts
261- );
206+ UTILS_ASSERT ( m_npts > 1 , " {} npts = {} not enought points\n " , msg, m_npts );
207+
262208 Utils::check_NaN ( m_X, msg+" X" , m_npts, __LINE__, __FILE__ );
263209 Utils::check_NaN ( m_Y, msg+" Y" , m_npts, __LINE__, __FILE__ );
210+
264211 integer ibegin { 0 };
265212 integer iend { 0 };
213+
266214 do {
267215 // cerca intervallo monotono strettamente crescente
268- while ( ++iend < m_npts && m_X[iend-1 ] < m_X[iend] ) {}
269- Pchip_build (
270- m_X+ibegin,
271- m_Y+ibegin,
272- m_Yp+ibegin,
273- iend-ibegin
274- );
216+ for ( ++iend; iend < m_npts && m_X[iend-1 ] < m_X[iend]; ++iend ) {}
217+ Pchip_build ( m_X+ibegin, m_Y+ibegin, m_Yp+ibegin, iend-ibegin );
275218 ibegin = iend;
276219 } while ( iend < m_npts );
277220
0 commit comments