VERB_code_2.2  2
 All Classes Namespaces Files Functions Variables Typedefs Macros Pages
akima.cpp
Go to the documentation of this file.
1 /*
2  * Akima.cpp
3  *
4  * Created on: Sep 3, 2010
5  * Author: subbotin
6  */
7 
8 #include "akima.h"
9 #include <iostream>
10 #include <stdlib.h>
11 
12 using namespace std;
13 
14 #define MERR 1e-19
15 
16 
17 bool akima(double *x0, double *y0, int nx, double *xi,double *yi, int nxi, int extrapolation_type, double ub, double lb, double ub_x, double lb_x){
18 
19  double max_f,x_lower,y_lower,x_upper,y_upper;
20  int exp_lower,exp_upper;
21 
22  double* x=new double[nx+2];
23  double* y=new double[nx+2];
24  double* dx=new double[nx+2];
25  double* dy=new double[nx+2];
26  double* diff=new double[nx+6];
27  double* ds=new double[nx+5];
28  double* f=new double[nx+2];
29  double* f1=new double[nx+2];
30  double* f2=new double[nx+2];
31  double* b=new double[nx+2];
32  double* c=new double[nx+2];
33  double* d=new double[nx+2];
34  double* wj=new double[nxi];
35  int* max_id=new int[nx+2];
36  int* index=new int[nxi];
37  int* bb=new int[nxi];
38 
39  /////////////////////////////////////////////////
40  // cout<<input x and y arrays must be of same length<<endl;
41  /////////////////////////////////////////////////
42 
43  if (extrapolation_type == 0) {
44  ////////// linear extrapolation //////////////
45  if(x0[0]>xi[0]) {
46  y_lower=y0[0]-(y0[1]-y0[0])/(x0[1]-x0[0])*(x0[0]-xi[0]);
47  x_lower=x0[0]-(x0[0]-xi[0]);
48  exp_lower=1;
49  } else {
50  exp_lower=0;
51  }
52  if(x0[nx-1]<xi[nxi-1]) {
53  y_upper=y0[nx-1]+(y0[nx-1]-y0[nx-2])/(x0[nx-1]-x0[nx-2])*(xi[nxi-1]-x0[nx-1]);
54  x_upper=x0[nx-1]+(xi[nxi-1]-x0[nx-1]);
55  exp_upper=1;
56  } else {
57  exp_upper=0;
58  }
59  } else if (extrapolation_type == 1) {
60  ////////// fixed boundaries //////////////
61  if(x0[0]>xi[0]) {
62  y_lower=lb;
63  x_lower=x0[0]-(x0[0]-xi[0]);
64  exp_lower=1;
65  } else {
66  exp_lower=0;
67  }
68  if(x0[nx-1]<xi[nxi-1]) {
69  y_upper=ub;
70  x_upper=x0[nx-1]+(xi[nxi-1]-x0[nx-1]);
71  exp_upper=1;
72  } else {
73  exp_upper=0;
74  }
75  } else if (extrapolation_type == 2) {
76  ////////// fixed boundaries at fixed x //////////////
77  y_lower=lb;
78  x_lower=lb_x;
79  exp_lower=1;
80  y_upper=ub;
81  x_upper=ub_x;
82  } else {
83  ////////// no extrapolation //////////////
84  if(x0[0]>xi[0]) {
85  cout << "extrapolation" << endl;
86  exit(0);
87  } else {
88  exp_lower=0;
89  }
90  if(x0[nx-1]<xi[nxi-1]) {
91  cout << "extrapolation" << endl;
92  exit(0);
93  } else {
94  exp_upper=0;
95  }
96  }
97 
98  ///// expand an input array x0 by adding two extrapolated boundary conditions /////////
99  if(exp_lower==1 && exp_upper==1) {
100  x[0]=x_lower;
101  y[0]=y_lower;
102  for(int i=0;i<nx;i++) {
103  x[i+1]=x0[i];
104  y[i+1]=y0[i];
105  }
106  x[nx+1]=x_upper;
107  y[nx+1]=y_upper;
108  nx=nx+2;
109  } else if(exp_lower==1 && exp_upper==0) {
110  x[0]=x_lower;
111  y[0]=y_lower;
112  for(int i=0;i<nx;i++) {
113  x[i+1]=x0[i];
114  y[i+1]=y0[i];
115  }
116  nx=nx+1;
117  } else if(exp_lower==0 && exp_upper==1){
118  for(int i=0;i<nx;i++) {
119  x[i]=x0[i];
120  y[i]=y0[i];
121  }
122  x[nx]=x_upper;
123  y[nx]=y_upper;
124  nx=nx+1;
125  } else {
126  for(int i=0;i<nx;i++) {
127  x[i]=x0[i];
128  y[i]=y0[i];
129  }
130  nx=nx;
131  }
132  ///////////////////////////////////////
133 
134  ///////////////// slopes //////////////////////
135  for (int i=0;i<nx-1;i++) { // 0 to nx-2
136  dx[i]=x[i+1]-x[i];
137  dy[i]=y[i+1]-y[i];
138  diff[i+2]=dy[i]/dx[i];
139  if(dx[i]<0.0){
140  cout<<"input x-array must be in strictly ascending order"<<endl;
141  //return false;
142  exit(0);
143 
144  }
145  }
146 
147  //diff[i+2]=(y[i+1]-y[i])/(x[i+1]-x[i]);
148 
149 
150  diff[1]=2.0*diff[2]-diff[3];
151  diff[0]=2.0*diff[1]-diff[2];
152  diff[nx+1]=2.0*diff[nx]-diff[nx-1];
153  diff[nx+2]=2.0*diff[nx+1]-diff[nx];
154  /////////////////////////////////////////////
155 
156 
157  for (int i=0;i<nx+2;i++) {
158  ds[i]=fabs(diff[i+1]-diff[i]);
159  }
160 
161  for (int i=0;i<nx;i++) {
162  f1[i]=ds[i+2];
163  f2[i]=ds[i];
164  f[i]=f1[i]+f2[i];
165  ///////////////////////////////////
166  b[i]=diff[i+1];
167  ///////////////////////////////////
168  }
169 
170  // find the maximum value of f
171  max_f=f[0];
172  for (int i=0;i<nx;i++) {
173  if(f[i]>max_f){
174  max_f=f[i];
175  }
176  }
177 
178 
179  int j=0;
180  int max_count=0;
181  for (int i=0;i<nx;i++) {
182  if(f[i]>max_f*1e-8){
183  max_id[j]=i;
184  j=j+1;
185  max_count=j;
186  }
187  }
188 
189  /*
190  for (int i=0;i<nx;i++) {
191  b[i]=diff[i+1];
192  }
193  */
194 
195  if(max_count!=0){
196  for (int i=0;i<max_count;i++) {
197  b[max_id[i]]=( f1[max_id[i]]*diff[max_id[i]+1]+f2[max_id[i]]*diff[max_id[i]+2] ) / f[max_id[i]];
198  }
199  }
200 
201 
202  for (int i=0;i<nx-1;i++) {
203  c[i]=(3.0*diff[i+2]-2.0*b[i]-b[i+1])/dx[i];
204  d[i]=(b[i]+b[i+1]-2.0*diff[i+2])/pow(dx[i],2);
205  }
206 
207  ///////////////// data binning ////////////////////
208  int i=0;
209  for(int j=0;j<nxi;j++){
210  if(xi[j]>=x[i] && xi[j]<x[i+1]){
211  index[j]=i;
212  //cout<<i<<" "<<x[i]<<" "<<x[i+1]<<" "<<xi[j]<<" "<<endl;
213  } else {
214  if(i==nx-1) {
215  for(int k=j-1;k<nxi;k++) {
216  index[k]=i;
217  // cout<<k<<" "<<xi[k]<<" "<<index[k]<<endl;
218  }
219  break;
220  } else {
221  i++;
222  j--;
223  }
224  }
225  }
226 
227 
228  /*
229  for (int i=0;i<nxi;i++) {
230  if(index[i]>nx-2) {
231  index[i]=nx-2;
232  }
233  //cout<<nx-1<<" index "<<index[i]<<endl;
234  }
235  */
236 
237  ///////////////////////////
238 
239 
240  for (int i=0;i<nxi;i++) {
241  //
242  if(index[i]>nx-2) {
243  index[i]=nx-2;
244  }
245  //
246  bb[i]=index[i];
247  wj[i]=xi[i]-x[bb[i]];
248  yi[i]=( (wj[i]*d[bb[i]]+c[bb[i]])*wj[i]+b[bb[i]] )
249  *wj[i]+y[bb[i]];
250  }
251 
252  delete []x;
253  delete []y;
254  delete []dx;
255  delete []dy;
256  delete []diff;
257  delete []ds;
258  delete []f;
259  delete []f1;
260  delete []f2;
261  delete []b;
262  delete []c;
263  delete []d;
264  delete []index;
265  delete []bb;
266  delete []wj;
267 
268  return true;
269 }
270 
271 
272 
273 
274 
275 
276 
277 /**
278  * Akima
279  * /param s - function
280  * /param x - grid
281  * /param k - number of points
282  */
283 
284 double akima2(double *x, double *y, int nx, double *xx, double *yy, int nnx, double ub, double lb) {
285  int i, ii;
286 
287  // derivatives
288  double* d = new double[nx];
289  double* w = new double[nx];
290  double* dx = new double[nx];
291  double* dy = new double[nx];
292 
293  double* p0 = new double[nx];
294  double* p1 = new double[nx];
295  double* p2 = new double[nx];
296  double* p3 = new double[nx];
297 
298  double* yder = new double[nx];
299 
300  double dm1, dm2, dnx, wm1, wnx;
301 
302  for (i = 0; i <= nx-2; i++) {
303  dx[i] = x[i+1] - x[i];
304  dy[i] = y[i+1] - y[i];
305  }
306 
307  for (i = 0; i <= nx-2; i++) {
308  d[i] = dy[i]/dx[i]; // d defined for i = 0 .. k-2
309  }
310 
311  dm1 = 2 * d[0] - d[1];
312  dm2 = 2 * dm1 - d[0];
313  d[nx-1] = 2 * d[nx-2] - d[nx-3];
314  dnx = 2 * d[nx-1] - d[nx-2];
315 
316  for (i = 1; i <= nx-2; i++) {
317  w[i] = fabs(d[i] - d[i-1]); // w defined for i = 1 .. k-2 !!!
318  }
319  w[0] = fabs(d[0] - dm1);
320  wm1 = fabs(dm1 - dm2);
321  w[nx-1] = fabs(d[nx-1] - d[nx-2]);
322  wnx =fabs(dnx - d[nx-1]);
323 
324  //t = (|m4 - m3|*m2 + |m2 - m1|*m3)/(|m4 - m3| + |m2 - m1|)
325  //t = (|d[4] - d[3]|*d[2] + |d[2] - d[1]|*d[3])/(|d[4] - d[3]| + |d[2] - d[1]|)
326  //t3 = (|d[i+1] - d[i]|*d[i-1] + |d[i-1] - d[i-2]|*d[i])/(|d[i+1] - d[i]| + |d[i-1] - d[i-2]|)
327  //t3 = (|w[i+1]|*d[i-1] + |w[i-1]|*d[i])/(|w[i+1]| + |w[i-1]|)
328  //t3 = (w[i+1]*d[i-1] + w[i-1]*d[i])/(w[i+1] + w[i-1]);
329 
330  for (i = 2; i <= nx-3; i++) {
331  if (fabs(d[i-2] - d[i-1]) < MERR && fabs(d[i] - d[i+1]) < MERR && fabs(d[i-1] - d[i]) > MERR) {
332  yder[i] = 0.5 * (d[i-1] + d[i]);
333  } else if (fabs(d[i-2] - d[i-1]) < MERR) {
334  yder[i] = d[i-1];
335  } else if (fabs(d[i+1] - d[i]) < MERR) {
336  yder[i] = d[i];
337  } else {
338  yder[i] = (w[i+1]*d[i-1] + w[i-1]*d[i])/(w[i+1] + w[i-1]);
339  }
340  }
341 
342  if (fabs(dm2 - dm1) < MERR && fabs(d[0] - d[1]) < MERR && fabs(dm1 - d[0]) > MERR) {
343  yder[0] = 0.5 * (dm1 + d[0]);
344  } else if (fabs(dm2 - dm1) < MERR) {
345  yder[0] = dm1;
346  } else if (fabs(d[1] - d[0]) < MERR) {
347  yder[0] = d[0];
348  } else {
349  yder[0] = (w[1]*dm1 + wm1*d[0])/(w[1] + wm1);
350  }
351 
352  if (fabs(dm1 - d[0]) < MERR && fabs(d[1] - d[2]) < MERR && fabs(d[0] - d[1]) > MERR) {
353  yder[1] = 0.5 * (d[0] + d[1]);
354  } else if (fabs(dm1 - d[0]) < MERR) {
355  yder[1] = d[0];
356  } else if (fabs(d[2] - d[1]) < MERR) {
357  yder[1] = d[1];
358  } else {
359  yder[1] = (w[2]*d[0] + w[0]*d[1])/(w[2] + w[0]);
360  }
361 
362  if (fabs(d[nx-3] - d[nx-2]) < MERR && fabs(d[nx-1] - dnx) < MERR && fabs(d[nx-2] - d[nx-1]) > MERR) {
363  yder[nx-1] = 0.5 * (d[nx-2] + d[nx-1]);
364  } else if (fabs(d[nx-3] - d[nx-2]) < MERR) {
365  yder[nx-1] = d[nx-2];
366  } else if (fabs(dnx - d[nx-1]) < MERR) {
367  yder[nx-1] = d[nx-1];
368  } else {
369  yder[nx-1] = (wnx*d[nx-2] + w[nx-2]*d[nx-1])/(wnx + w[nx-2]);
370  }
371 
372  if (fabs(d[nx-4] - d[nx-3]) < MERR && fabs(d[nx-2] - d[nx-1]) < MERR && fabs(d[nx-3] - d[nx-2]) > MERR) {
373  yder[nx-2] = 0.5 * (d[nx-3] + d[nx-2]);
374  } else if (fabs(d[nx-4] - d[nx-3]) < MERR) {
375  yder[nx-2] = d[nx-3];
376  } else if (fabs(d[nx-1] - d[nx-2]) < MERR) {
377  yder[nx-2] = d[nx-2];
378  } else {
379  yder[nx-2] = (w[nx-1]*d[nx-3] + w[nx-3]*d[nx-2])/(w[nx-1] + w[nx-3]);
380  }
381 
382  //for (i = 2; i <= nx-3; i++) {
383  for (i = 0; i <= nx-2; i++) {
384  p0[i] = y[i];
385  p1[i] = yder[i];
386  p2[i] = (3*d[i] - 2*yder[i] - yder[i+1]) / dx[i];
387  p3[i] = (yder[i] + yder[i+1] - 2*d[i]) / dx[i] / dx[i];
388  //p2[i] = (3*d[i] - 2*w[i] - w[i+1]) / dx[i];
389  //p3[i] = (w[i] + w[i+1] - 2*d[i]) / dx[i] / dx[i];
390  }
391 
392  i = 0;
393  // assume grids values are 'sorted'
394  for (ii = 0; ii <= nnx-1; ii++) {
395  while (i <= nx-2 && x[i+1] <= xx[ii]) i++;
396  if (fabs(x[i] - xx[ii]) < MERR) {
397  yy[ii] = y[i];
398  } else if (fabs(x[i+1] - xx[ii]) < MERR) {
399  yy[ii] = y[i+1];
400  } else if (xx[ii] > x[i] && xx[ii] < x[i+1]) {
401  yy[ii] = p0[i] + p1[i]*(xx[ii] - x[i]) + p2[i]*pow((xx[ii] - x[i]), 2) + p3[i]*pow((xx[ii] - x[i]), 3);
402  } else {
403  // something wrong
404  exit(0);
405  }
406  }
407 
408 
409  return -1;
410 }
411 
412