VERB_code_2.2
2
Main Page
Related Pages
Namespaces
Classes
Files
File List
File Members
All
Classes
Namespaces
Files
Functions
Variables
Typedefs
Macros
Pages
polint.cpp
Go to the documentation of this file.
1
/**
2
* Some other interpolation from numerical recepies.
3
* \file polint.cpp
4
*/
5
6
//#include "../../Logging/Output.h"
7
#include "math.h"
8
#include <malloc.h>
9
#include <stdlib.h>
10
#include <cstdio>
11
12
/**
13
* Some other interpolation from numerical recepies.
14
*/
15
void
polint
(
double
*xa,
double
*ya,
int
n,
double
x,
double
*y,
double
*dy )
16
{
17
double
*c = NULL;
18
double
*d = NULL;
19
double
den;
20
double
dif;
21
double
dift;
22
double
ho;
23
double
hp;
24
double
w;
25
26
int
i;
27
int
m
;
28
int
ns;
29
30
if
( (c = (
double
*)malloc( n *
sizeof
(
double
) )) == NULL || (d = (
double
*)malloc( n *
sizeof
(
double
) )) == NULL ) {
31
//fprintf( stderr, "polint error: allocating workspace\n" );
32
//Output::echo(0, "polint error: allocating workspace\n" );
33
//fprintf( stderr, "polint error: setting y = 0 and dy = 1e9\n" );
34
//Output::echo(0, "polint error: setting y = 0 and dy = 1e9\n" );
35
*y = 0.0;
36
*dy = 1.e9;
37
38
if
( c != NULL )
39
free( c );
40
if
( d != NULL )
41
free( d );
42
return
;
43
}
44
45
ns = 0;
46
dif = fabs(x-xa[0]);
47
for
( i = 0; i < n; ++i ) {
48
dift = fabs( x-xa[i] );
49
if
( dift < dif ) {
50
ns = i;
51
dif = dift;
52
}
53
c[i] = ya[i];
54
d[i] = ya[i];
55
}
56
*y = ya[ns];
57
ns = ns-1;
58
//if (ns < 0) ns = 0; ////////////////
59
for
( m = 0; m < n-1; ++
m
) {
60
for
( i = 0; i < n-m-1; ++i ) {
61
ho = xa[i]-x;
62
hp = xa[i+m+1]-x;
63
w = c[i+1]-d[i];
64
den = ho-hp;
65
if
( fabs(den) < 1.e-25 ) {
66
fprintf( stderr,
"polint error: den = 0\n"
);
67
//Output::echo(0, "polint error: den = 0\n" );
68
fprintf( stderr,
"polint error: setting y = 0 and dy = 1e9\n"
);
69
//Output::echo(0, "polint error: setting y = 0 and dy = 1e9\n" );
70
exit(0);
71
*y = 0.0;
72
*dy = 1.e9;
73
74
if
( c != NULL )
75
free( c );
76
if
( d != NULL )
77
free( d );
78
return
;
79
}
80
den = w/den;
81
d[i] = hp*den;
82
c[i] = ho*den;
83
}
84
if
( 2*(ns+1) < n-m-1 ) {
85
*dy = c[ns+1];
86
}
else
{
87
*dy = d[ns];
88
ns = ns-1;
89
//if (ns < 0) ns = 0; //////////
90
}
91
*y = (*y)+(*dy);
92
}
93
94
95
if
( c != NULL )
96
free( c );
97
if
( d != NULL )
98
free( d );
99
return
;
100
}
Matrix
Interpolation
polint.cpp
Generated on Sat Nov 16 2013 09:31:46 for VERB_code_2.2 by
1.8.4