1 /* Given a set of data points, generate a spline curve that passes thru all points.
2 *
3 * Author - Marlow etc. Modified by - P. Ward Date - 3.10.1973
4 *
5 * Extracted from the ROOT package by O.Couet.
6 *
7 * Re-structured to remove most gotos by Dino Ferrero Merlino
8 * (the numerical part has not been touched)
9 *
10 * Changed to be ANSI-C compliant by Oliver Koch.
11 *
12 * Integrated into ploticus by Steve Grubb 10/3/03
13 *
14 * Usage:
15 *
16 * Helps to compute a smooth representation of a curve.
17 * The history of this code is quite long...
18 *
19 * This routine draws a smooth tangentially continuous curve through
20 * the sequence of data points P(I) I=1,N where P(I)=(X(I),Y(I))
21 * the curve is approximated by a polygonal arc of short vectors .
22 * the data points can represent open curves, P(1) != P(N) or closed
23 * curves P(2) == P(N) . If a tangential discontinuity at P(I) is
24 * required , then set P(I)=P(I+1) . loops are also allowed .
25 *
26 * Reference Marlow and Powell,Harwell report No.R.7092.1972
27 * MCCONALOGUE,Computer Journal VOL.13,NO4,NOV1970PP392 6
28 *
29 * The routine was then adapted as the IGRAP1 routine of HIGZ.
30 * The FORTRAN code was translated with f2c and adapted to ROOT.
31 *
32 *
33 * TODO:
34 *
35 */
36
37 #include <stdio.h>
38 #include <stdlib.h>
39 #include <math.h>
40
41 #define false 0;
42 #define true 1;
43
44 static int npoints, k, kp, NpointsMax, banksize, npt;
45 static double xratio, yratio, sxmin, symin;
46 static int closed, flgis, iw;
47 static double P1, P2, P3, P4, P5, P6, W1, W2, W3, A, B, C, R, S, T, Z;
48 static double CO, SO, CT, ST, CTU, STU, XNT, DX1, DY1, DX2, DY2, XO, YO, DX, DY, XT, YT;
49 static double XA, XB, YA, YB, U1, U2, U3, TJ, SB, STH;
50 static double wsign, tsquare, tcube, *qlx, *qly, *x, *y;
51 static const int maxiterations = 20;
52 static const double delta = 0.00055;
53
54 static int QpSmootherInit( int nPoints, double in[][2] );
55 static int QpSmootherDone();
56 static int continousTangentAtEndPoints() ;
57 static int scaleRatio() ;
58 static int skipConsecutiveEqualPoints(int start) ;
59 static void computeDirectionCosine (int kMinus1, int kP, int kPlus1) ;
60 static int resizeSmoothArrays() ;
61 static void computeStraightLine(int kP) ;
62 static void newPoint() ;
63 static void dirCosL150() ;
64 static void dirCosL130() ;
65 static void computeCubicCoeffs() ;
66 static int smooth() ;
67 static void zero(int *kCode, double AZ, double BZ, double E2, double *X, double *Y, int maxiterations);
68 static void checkMaxIterations(int *kCode, int *IT, int J1, double *X, double X2, int maxiterations);
69 /* static int getNPoint() ; */
70 static double* getX ();
71 static double* getY ();
72 static double min( double f, double g );
73 static double max( double f, double g );
74
75
76 /* --------------------------------------------------- */
77 /* SMOOTHFIT - ploticus entry point */
78 int
PL_smoothfit(in,ninp,out,noutp)79 PL_smoothfit( in, ninp, out, noutp )
80 double in[][2]; /* input vector */
81 int ninp; /* # input points */
82 double out[][2]; /* output vector */
83 int *noutp; /* # output points */
84 {
85 int i, stat;
86 double *resultx, *resulty;
87
88 stat = QpSmootherInit( ninp, in );
89 if( stat ) return( stat );
90
91 *noutp = smooth();
92
93 /* *noutp = getNPoint(); */
94
95 resultx = getX();
96 resulty = getY();
97
98 for( i = 0; i < *noutp; i++ ) {
99 out[i][0] = resultx[i];
100 out[i][1] = resulty[i];
101 }
102
103 stat = QpSmootherDone();
104
105 return( 0 );
106 }
107
108
109 /* ------------------------------- */
110
111 /* void QpSmootherInit(int nPoints, double *_x, double *_y) */
112
QpSmootherInit(nPoints,in)113 static int QpSmootherInit( nPoints, in )
114 int nPoints;
115 double in[][2];
116 {
117 int i;
118
119 C = T = CO = SO = CT = ST = CTU = STU = DX1 = DY1 = DX2 = DY2 = 0;
120 XT = YT = XA = XB = YA = YB = U1 = U2 = U3 = TJ = SB = 0;
121 /* make a local copy of the original polylines (smoother will scale them) */
122 npoints = nPoints;
123 x = (double *) malloc(sizeof(double) * npoints);
124 y = (double *) malloc(sizeof(double) * npoints);
125 if (x == NULL || y == NULL) return( 5720 ); /* Not enough memory for qpsmoother */
126 else {
127 for (i=0;i<npoints;i++) {
128 x[i] = in[i][0];
129 y[i] = in[i][1];
130 }
131
132 NpointsMax = npoints*10;
133 banksize = NpointsMax;
134
135 qlx = (double *) malloc(sizeof(double) * NpointsMax);
136 qly = (double *) malloc(sizeof(double) * NpointsMax);
137 if (qlx == NULL || qly == NULL) return( 5720 ); /* Not enough memory for qpsmoother */
138 else {
139 qlx[0] = x[0];
140 qly[0] = y[0];
141 /* Scale data and check whether endpoints collide */
142 closed = scaleRatio();
143 }
144 }
145 return( 0 );
146 }
147
148 /* ------------------------------- */
QpSmootherDone()149 static int QpSmootherDone() {
150 free(x);
151 free(y);
152 free(qlx);
153 free(qly);
154 return( 0 );
155 }
156
157 /* ------------------------------- */
continousTangentAtEndPoints()158 static int continousTangentAtEndPoints() {
159 int tangent = false;
160
161 if (x[0] != x[npoints-1] || y[0] != y[npoints-1])
162 tangent = true;
163 if (x[npoints-2] == x[npoints-1] && y[npoints-2] == y[npoints-1])
164 tangent = true;
165 if (x[0] == x[1] && y[0] == y[1])
166 tangent = true;
167 return tangent;
168 }
169
170 /* ------------------------------- */
171 /** Scale data to the range 0-ratio_signs in X, 0-1 in Y
172 where ratio_signs is the ratio between the number of changes
173 of sign in Y divided by the number of changes of sign in X.
174 Returns whether the endpoints are overlapping or not.
175 */
176
scaleRatio()177 static int scaleRatio() {
178 double ratio_signs;
179 double sxmax, symax;
180 int i;
181 double six;
182 double siy;
183 double dx1n;
184 double dy1n;
185
186 closed = false;
187 sxmin = x[0];
188 sxmax = x[0];
189 symin = y[0];
190 symax = y[0];
191 six = 1;
192 siy = 1;
193 for (i=1;i<npoints;i++) {
194 if (i > 1) {
195 if ((x[i]-x[i-1])*(x[i-1]-x[i-2]) < 0) six++;
196 if ((y[i]-y[i-1])*(y[i-1]-y[i-2]) < 0) siy++;
197 }
198 if (x[i] < sxmin) sxmin = x[i];
199 if (x[i] > sxmax) sxmax = x[i];
200 if (y[i] < symin) symin = y[i];
201 if (y[i] > symax) symax = y[i];
202 }
203 dx1n = fabs(x[npoints-1]-x[0]);
204 dy1n = fabs(y[npoints-1]-y[0]);
205 if (dx1n < 0.01*(sxmax-sxmin) && dy1n < 0.01*(symax-symin))
206 closed = true;
207 if (sxmin == sxmax)
208 xratio = 1;
209 else {
210 if (six > 1)
211 ratio_signs = siy/six;
212 else
213 ratio_signs = 20;
214 xratio = ratio_signs/(sxmax-sxmin);
215 }
216 if (symin == symax)
217 yratio = 1;
218 else
219 yratio = 1/(symax-symin);
220
221 for (i=0;i<npoints;i++) {
222 x[i] = (x[i]-sxmin)*xratio;
223 y[i] = (y[i]-symin)*yratio;
224 }
225
226 return closed;
227 }
228
229 /* ------------------------------- */
skipConsecutiveEqualPoints(int start)230 static int skipConsecutiveEqualPoints(int start) {
231 k = start;
232 while((k < npoints) && (x[k] == x[k-1] && y[k] == y[k-1])) {
233 k++;
234 }
235 return k;
236 }
237
238 /* ------------------------------- */
239 /* Calculate the direction cosines at P(k) obtained from P(kMinus1),P(k),P(kPlus1).
240 * If k=1 then N-1 is used for kMinus1. If k=N then 2 is used for kPlus1
241 */
242
computeDirectionCosine(int kMinus1,int kP,int kPlus1)243 static void computeDirectionCosine (int kMinus1, int kP, int kPlus1) {
244 double DK1,DK2;
245 DX1 = x[kP] - x[kMinus1];
246 DY1 = y[kP] - y[kMinus1];
247 DX2 = x[kPlus1] - x[kP];
248 DY2 = y[kPlus1] - y[kP];
249 DK1 = DX1*DX1 + DY1*DY1;
250 DK2 = DX2*DX2 + DY2*DY2;
251 CTU = DX1*DK2 + DX2*DK1;
252 STU = DY1*DK2 + DY2*DK1;
253 XNT = CTU*CTU + STU*STU;
254 /* If both ctu and stu are zero,then default.This can
255 * occur when P(k)=P(k+1). I.E. A loop.
256 */
257 if (XNT < 1.E-25) {
258 CTU = DY1;
259 STU =-DX1;
260 XNT = DK1;
261 }
262 /* Normalise direction cosines. */
263 CT = CTU/sqrt(XNT);
264 ST = STU/sqrt(XNT);
265 }
266
267
268 /* ------------------------------- */
resizeSmoothArrays()269 static int resizeSmoothArrays() {
270 int i;
271 int newsize = banksize + NpointsMax;
272 double *qt1=qlx;
273 double *qt2=qly;
274 qlx = (double *) malloc(sizeof(double) * newsize);
275 qly = (double *) malloc(sizeof(double) * newsize);
276 for (i=0;i<banksize;i++) {
277 qlx[i] = qt1[i];
278 qly[i] = qt2[i];
279 }
280 free(qt1);
281 free(qt2);
282 banksize = newsize;
283 return newsize;
284 }
285
286 /* ------------------------------- */
computeStraightLine(int kP)287 static void computeStraightLine(int kP) {
288 /* Compute a straight line between P(k-1) and P(k) in Fortran
289 * (XT,YT) will contain the coordinate of the endpoint
290 */
291 XT = x[kP] ;
292 YT = y[kP] ;
293 }
294
295 /* ------------------------------- */
newPoint()296 static void newPoint() {
297 /* Convert coordinates to original system */
298 qlx[npt] = sxmin + XT/xratio;
299 qly[npt] = symin + YT/yratio;
300 npt++;
301 /* If the arrays are filled up, enlarge them. */
302 if (npt >= banksize)
303 banksize = resizeSmoothArrays();
304 }
305
306 /* ------------------------------- */
307 /* Direction cosines at P(k) obtained from P(k-2),P(k-1),P(k). */
dirCosL150()308 static void dirCosL150() {
309 W3 = 2*(DX1*DY2-DX2*DY1);
310 CT = CTU-W3*DY2;
311 ST = STU+W3*DX2;
312 XNT = 1/sqrt(CT*CT+ST*ST);
313 CT = CT*XNT;
314 ST = ST*XNT;
315 flgis = 0;
316 }
317
318 /* ------------------------------- */
319 /* Direction cosines at P(k-1) obtained from P(k-1),P(k),P(k+1). */
dirCosL130()320 static void dirCosL130() {
321 W3 = 2*(DX1*DY2-DX2*DY1);
322 CO = CTU+W3*DY1;
323 SO = STU-W3*DX1;
324 XNT = 1/sqrt(CO*CO+SO*SO);
325 CO = CO*XNT;
326 SO = SO*XNT;
327 flgis = 1;
328 }
329
330
331 /* ------------------------------- */
332
333 /* For the arc between P(k-1) and P(k) with direction cosines CO,
334 * SO and CT,ST respectively, calculate the coefficients of the
335 * parametric cubic represented by X(T) and Y(T) where
336 * X(T)=XA*T**3 + XB*T**2 + CO*T + XO
337 * Y(T)=YA*T**3 + YB*T**2 + SO*T + YO
338 */
339
computeCubicCoeffs()340 static void computeCubicCoeffs() {
341 double CC ;
342 double ERR ;
343 int zeroLoop;
344
345 XO = x[k-2];
346 YO = y[k-2];
347 DX = x[k-1] - XO;
348 DY = y[k-1] - YO;
349
350 /* Initialise the values of X(TI),Y(TI) in XT and YT respectively. */
351
352 XT = XO;
353 YT = YO;
354 C = DX*DX+DY*DY;
355 A = CO+CT;
356 B = SO+ST;
357 R = DX*A+DY*B;
358 T = C*6/(sqrt(R*R+2*(7-CO*CT-SO*ST)*C)+R);
359 tsquare = T*T;
360 tcube = T*tsquare;
361 XA = (A*T-2*DX)/tcube;
362 XB = (3*DX-(CO+A)*T)/tsquare;
363 YA = (B*T-2*DY)/tcube;
364 YB = (3*DY-(SO+B)*T)/tsquare;
365
366 /* If the curve is close to a straight line then use a straight
367 * line between (XO,YO) and (XT,YT).
368 */
369
370 if (.75*max(fabs(DX*SO-DY*CO),fabs(DX*ST-DY*CT)) <= delta) {
371 computeStraightLine(k-1);
372 newPoint();
373 } else {
374
375 /* Calculate a set of values 0 == T(0).LTCT(1) < ... < T(M)=TC
376 * such that polygonal arc joining X(T(J)),Y(T(J)) (J=0,1,..M)
377 * is within the required accuracy of the curve
378 */
379
380 TJ = 0;
381 U1 = YA*XB-YB*XA;
382 U2 = YB*CO-XB*SO;
383 U3 = SO*XA-YA*CO;
384
385 /* Given T(J), calculate T(J+1). The values of X(T(J)), */
386 /* Y(T(J)) T(J) are contained in XT,YT and TJ respectively. */
387
388 do {
389 S = T - TJ;
390 iw = -2;
391
392 /* Define iw here later. */
393
394 P1 = (2*U1)*TJ-U3;
395 P2 = (U1*TJ-U3)*3*TJ+U2;
396 P3 = 3*TJ*YA+YB;
397 P4 = (P3+YB)*TJ+SO;
398 P5 = 3*TJ*XA+XB;
399 P6 = (P5+XB)*TJ+CO;
400 CC = 0.8209285;
401 ERR = 0.1209835;
402 iw -= 2;
403
404
405 L200:
406 /* Test D(TJ,THETA). A is set to (Y(TJ+S)-Y(TJ))/S.B is
407 * set to (X(TJ+S)-X(TJ))/S.
408 */
409 A = (S*YA+P3)*S+P4;
410 B = (S*XA+P5)*S+P6;
411
412 /* Set Z to PSI(D/delta)-CC. */
413
414 W1 = -S*(S*U1+P1);
415 W2 = S*S*U1-P2;
416 W3 = 1.5*W1+W2;
417
418 /* Set the estimate of (THETA-TJ)/S.Then set the numerator
419 * of the expression (EQUATION 4.4)/S. Then set the square
420 * of D(TJ,TJ+S)/delta. Then replace Z by PSI(D/delta)-CC.
421 */
422
423 if (W3 > 0)
424 wsign = fabs(W1);
425 else
426 wsign = -fabs(W1);
427 STH = 0.5+wsign/(3.4*fabs(W1)+5.2*fabs(W3));
428 Z = S*STH*(S-S*STH)*(W1*STH+W1+W2);
429 Z = Z*Z/((A*A+B*B)*(delta*delta));
430 Z = (Z+2.642937)*Z/((.3715652*Z+3.063444)*Z+.2441889)-CC;
431
432 /* Branch if Z has been calculated */
433 zeroLoop=true;
434 if (iw <= 0) {
435 if (Z > ERR) {
436 kp = 0;
437 C = Z;
438 SB = S;
439 } else {
440 L220:
441 if (iw+2 == 0) {
442 /* goto L190; */
443 iw -= 2;
444 goto L200;
445 }
446
447 if (iw+2 > 0) {
448 /* Update TJ,XT and YT. */
449 XT = XT + S*B;
450 YT = YT + S*A;
451 TJ = S + TJ;
452 } else {
453 /* Last part of arc. */
454 XT = x[k-1];
455 YT = y[k-1];
456 S = 0;
457 }
458 newPoint();
459 zeroLoop=false; /* passato qui. S>0 */
460 }
461 }
462
463 /* Z(S). find a value of S where 0 <= S <= SB such that */
464 /* fabs(Z(S)) < ERR */
465
466 while(zeroLoop) {
467 zero(&kp,0,SB,ERR,&S,&Z,maxiterations);
468 if (kp == 2) {
469 /* goto L210; */
470 iw -= 2;
471 goto L220;
472 }
473 if (kp > 2) {
474 fprintf( stderr, "Attempt to plot outside plot limits\n" );
475 XT = x[k-1];
476 YT = y[k-1];
477 newPoint();
478 /* goto L310; */
479 break;
480 }
481 if (iw > 0) goto L200;
482
483 /* Set Z=Z(S) for S=0. */
484
485 if (iw < 0) {
486 Z = -CC;
487 iw = 0;
488 } else {
489 /* Set Z=Z(S) for S=SB. */
490 Z = C;
491 iw = 1;
492 }
493 } /* End while (zeroLoop) */
494 } while ( S > 0 );
495 }
496 }
497
498 /* ------------------------------- */
smooth()499 static int smooth() {
500 npt = 1;
501 k = 1;
502
503 /* flgis is true in the main loop, but is false if there is
504 * a deviation from the main loop (cusp or endpoint? DinoFM)
505 */
506
507 /* The curve is open if :
508 * the two endpoints are not overlapping (closed == false)
509 * AND there is a continuous tangent at the endpoints */
510
511 if (!closed && continousTangentAtEndPoints()) {
512 closed = true;
513 flgis = 0;
514 k++;
515 } else {
516 closed = false;
517 flgis = 1;
518 /* Calculate direction cosines at P(0) using P(N-1),P(0),P(1). */
519 computeDirectionCosine(npoints-2,0,1);
520 /* Carry forward the direction cosines from the previous arc. */
521 CO = CT;
522 SO = ST;
523 k++;
524 }
525
526 while (k <= npoints) {
527 k=skipConsecutiveEqualPoints(k);
528 /* More arcs to compute */
529 if (k < npoints) {
530 /* Test whether P(k) is a cusp. */
531 if (x[k-1] == x[k] && y[k-1] == y[k]) {
532 if (flgis) {
533 dirCosL150();
534 computeCubicCoeffs();
535 } else {
536 /* Make a straight line from P(k-1) to P(k). */
537 computeStraightLine(k-1);
538 newPoint();
539 }
540 } else {
541 /* k is always >= 2 ! */
542 computeDirectionCosine(k-2,k-1,k);
543 if (!flgis)
544 dirCosL130();
545 computeCubicCoeffs();
546 }
547 } else {
548 /* k == npoints. Last arc ! */
549 if (!closed) {
550 /* k is always >= 2 ! */
551 computeDirectionCosine(k-2,k-1,1);
552 if (!flgis)
553 dirCosL130();
554 computeCubicCoeffs();
555 } else {
556 if (flgis) {
557 dirCosL150();
558 computeCubicCoeffs();
559 } else {
560 /* Make a straight line from P(k-1) to P(k). */
561 computeStraightLine(k-1);
562 newPoint();
563 }
564 }
565 }
566 /* Branch if the next section of the curve begins at a cusp. */
567 if (flgis) {
568 /* Carry forward the direction cosines from the previous arc. */
569 CO = CT;
570 SO = ST;
571 }
572 k++;
573 }
574 return npt;
575 }
576
577 /* ------------------------------- */
zero(int * kCode,double AZ,double BZ,double E2,double * X,double * Y,int maxiterations)578 static void zero(int *kCode, double AZ, double BZ, double E2, double *X, double *Y, int maxiterations)
579 {
580 static double AA, BB, YAA, ytest, Y1, X1, H;
581 static int J1, IT, J3, J2;
582 double YBB, X2;
583 int exitLoop;
584 YBB = 0;
585 exitLoop=false;
586
587 /* Calculate Y(X) at X=AZ. */
588 if ((*kCode) <= 0) {
589 AA = AZ;
590 BB = BZ;
591 (*X) = A;
592 J1 = 1;
593 IT = 1;
594 (*kCode) = J1;
595 return;
596 }
597
598 /* Test whether Y(X) is sufficiently small. */
599
600 if (fabs((*Y)) <= E2) { (*kCode) = 2; return; }
601
602 /* Calculate Y(X) at X=BZ. */
603
604 if (J1 == 1) {
605 YAA = (*Y);
606 (*X) = BB;
607 J1 = 2;
608 return;
609 }
610
611 /* Test whether the signs of Y(AZ) and Y(BZ) are different.
612 * if not, begin the binary subdivision.
613 */
614
615 if ((J1 == 2) && (YAA*(*Y) >= 0)) {
616 X1 = AA;
617 Y1 = YAA;
618 J1 = 3;
619 H = BB - AA;
620 J2 = 1;
621 X2 = AA + 0.5*H;
622 J3 = 1;
623 checkMaxIterations(&k,&IT,J1,X,X2,maxiterations);
624 return;
625 }
626
627 /* Test whether a bracket has been found .
628 * If not,continue the search
629 */
630
631 if ((J1 != 2) || (YAA*(*Y) >= 0)) {
632 if (J1 > 3) goto L170;
633 if (YAA*(*Y) >= 0) {
634 if (J3 >= J2) {
635 H = 0.5*H; J2 = 2*J2;
636 AA = X1; YAA = Y1; X2 = AA + 0.5*H; J3 = 1;
637 }
638 else {
639 AA = (*X); YAA = (*Y); X2 = (*X) + H; J3++;
640 }
641 checkMaxIterations(&k,&IT,J1,X,X2,maxiterations);
642 return;
643 }
644 }
645
646 /* The first bracket has been found.calculate the next X by the
647 * secant method based on the bracket.
648 */
649
650 BB = (*X);
651 YBB = (*Y);
652 J1 = 4;
653
654 if (fabs(YAA) > fabs(YBB)) {
655 X1 = AA; Y1 = YAA; (*X) = BB; (*Y) = YBB;
656 } else {
657 X1 = BB; Y1 = YBB; (*X) = AA; (*Y) = YAA;
658 }
659
660 /* Use the secant method based on the function values Y1 and Y.
661 * check that X2 is inside the interval (A,B).
662 */
663
664
665 do {
666 X2 = (*X)-(*Y)*((*X)-X1)/((*Y)-Y1);
667 X1 = (*X);
668 Y1 = (*Y);
669 ytest = 0.5*min(fabs(YAA),fabs(YBB));
670 if ((X2-AA)*(X2-BB) < 0) {
671 checkMaxIterations(kCode,&IT,J1,X,X2,maxiterations);
672 return;
673 }
674
675 /* Calculate the next value of X by bisection . Check whether
676 * the maximum accuracy has been achieved.
677 */
678
679 X2 = 0.5*(AA+BB);
680 ytest = 0;
681 if ((X2-AA)*(X2-BB) >= 0) {
682 (*kCode) = 2; return;
683 }
684 checkMaxIterations(kCode,&IT,J1,X,X2,maxiterations);
685 return;
686
687 /* Revise the bracket (A,B). */
688
689 L170:
690 if (J1 != 4) return;
691 if (YAA*(*Y) < 0) {
692 BB = (*X); YBB = (*Y);
693 } else {
694 AA = (*X); YAA = (*Y);
695 }
696
697 /* Use ytest to decide the method for the next value of X. */
698
699 if (ytest <= 0) {
700 if (fabs(YAA) > fabs(YBB)) {
701 X1 = AA; Y1 = YAA; (*X) = BB; (*Y) = YBB;
702 } else {
703 X1 = BB; Y1 = YBB; (*X) = AA; (*Y) = YAA;
704 }
705 } else {
706 if (fabs((*Y))-ytest > 0)
707 exitLoop=true;
708 }
709 } while (!exitLoop);
710
711 X2 = 0.5*(AA+BB);
712 ytest = 0;
713 if ((X2-AA)*(X2-BB) >= 0) {
714 k = 2;
715 return;
716 }
717 checkMaxIterations(&k,&IT,J1,X,X2,maxiterations);
718 }
719
720 /* ------------------------------- */
721 /* Check whether (maxiterations) function values have been calculated. */
722
checkMaxIterations(int * kCode,int * IT,int J1,double * X,double X2,int maxiterations)723 static void checkMaxIterations(int *kCode, int *IT, int J1, double *X, double X2, int maxiterations)
724 {
725 (*IT)++;
726 if ((*IT) >= maxiterations) {
727 (*kCode) = J1;
728 } else {
729 (*X) = X2;
730 }
731 }
732
733 /* ------------------------------- */
734 /*
735 * static int getNPoint() {
736 * return(npoints);
737 * }
738 */
739
740 /* ------------------------------- */
getX()741 static double* getX () {
742 return(qlx);
743 }
744
745 /* ------------------------------- */
getY()746 static double* getY () {
747 return(qly);
748 }
749
750 /* ------------------------------- */
min(f,g)751 static double min( f, g )
752 double f, g;
753 {
754 if( f < g ) return( f );
755 else return( g );
756 }
757
758 /* ------------------------------- */
max(f,g)759 static double max( f, g )
760 double f, g;
761 {
762 if( f > g ) return( f );
763 else return( g );
764 }
765