1 #include "splines/oskar_dierckx_fprati.h"
2 
3 #ifdef __cplusplus
4 extern "C" {
5 #endif
6 
oskar_dierckx_fprati(double * p1,double * f1,double p2,double f2,double * p3,double * f3)7 double oskar_dierckx_fprati(double *p1, double *f1, double p2, double f2,
8         double *p3, double *f3)
9 {
10     /* Local variables */
11     double p = 0.0, h1 = 0.0, h2 = 0.0, h3 = 0.0;
12 
13     if (*p3 > 0.0)
14     {
15         /* value of p in case p3 ^= infinity. */
16         h1 = *f1 * (f2 - *f3);
17         h2 = f2 * (*f3 - *f1);
18         h3 = *f3 * (*f1 - f2);
19         p = -(*p1 * p2 * h3 + p2 * *p3 * h1 + *p3 * *p1 * h2) /
20                 (*p1 * h1 + p2 * h2 + *p3 * h3);
21     }
22     else
23     {
24         /* value of p in case p3 = infinity. */
25         p = (*p1 * (*f1 - *f3) * f2 - p2 * (f2 - *f3) * *f1) /
26                 ((*f1 - f2) * *f3);
27     }
28 
29     /* adjust the value of p1,f1,p3 and f3 such that f1 > 0 and f3 < 0. */
30     if (f2 < 0.0)
31     {
32         *p3 = p2;
33         *f3 = f2;
34     }
35     else
36     {
37         *p1 = p2;
38         *f1 = f2;
39     }
40     return p;
41 }
42 
43 #ifdef __cplusplus
44 }
45 #endif
46