VERB_code_2.3
polint.cpp
Go to the documentation of this file.
1 
6 //#include "../../Logging/Output.h"
7 #include "math.h"
8 #include <malloc/malloc.h>
9 #include <stdlib.h>
10 #include <cstdio>
11 
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 }
static const double m
mass of electron in grams
void polint(double *xa, double *ya, int n, double x, double *y, double *dy)
Definition: polint.cpp:15