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){
19 double max_f,x_lower,y_lower,x_upper,y_upper;
20 int exp_lower,exp_upper;
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];
43 if (extrapolation_type == 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]);
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]);
59 }
else if (extrapolation_type == 1) {
63 x_lower=x0[0]-(x0[0]-xi[0]);
68 if(x0[nx-1]<xi[nxi-1]) {
70 x_upper=x0[nx-1]+(xi[nxi-1]-x0[nx-1]);
75 }
else if (extrapolation_type == 2) {
85 cout <<
"extrapolation" << endl;
90 if(x0[nx-1]<xi[nxi-1]) {
91 cout <<
"extrapolation" << endl;
99 if(exp_lower==1 && exp_upper==1) {
102 for(
int i=0;i<nx;i++) {
109 }
else if(exp_lower==1 && exp_upper==0) {
112 for(
int i=0;i<nx;i++) {
117 }
else if(exp_lower==0 && exp_upper==1){
118 for(
int i=0;i<nx;i++) {
126 for(
int i=0;i<nx;i++) {
135 for (
int i=0;i<nx-1;i++) {
138 diff[i+2]=dy[i]/dx[i];
140 cout<<
"input x-array must be in strictly ascending order"<<endl;
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];
157 for (
int i=0;i<nx+2;i++) {
158 ds[i]=fabs(diff[i+1]-diff[i]);
161 for (
int i=0;i<nx;i++) {
172 for (
int i=0;i<nx;i++) {
181 for (
int i=0;i<nx;i++) {
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]];
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);
209 for(
int j=0;j<nxi;j++){
210 if(xi[j]>=x[i] && xi[j]<x[i+1]){
215 for(
int k=j-1;k<nxi;k++) {
240 for (
int i=0;i<nxi;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]] )
284 double akima2(
double *x,
double *y,
int nx,
double *xx,
double *yy,
int nnx,
double ub,
double lb) {
288 double* d =
new double[nx];
289 double* w =
new double[nx];
290 double* dx =
new double[nx];
291 double* dy =
new double[nx];
293 double* p0 =
new double[nx];
294 double* p1 =
new double[nx];
295 double* p2 =
new double[nx];
296 double* p3 =
new double[nx];
298 double* yder =
new double[nx];
300 double dm1, dm2, dnx, wm1, wnx;
302 for (i = 0; i <= nx-2; i++) {
303 dx[i] = x[i+1] - x[i];
304 dy[i] = y[i+1] - y[i];
307 for (i = 0; i <= nx-2; i++) {
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];
316 for (i = 1; i <= nx-2; i++) {
317 w[i] = fabs(d[i] - d[i-1]);
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]);
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) {
335 }
else if (fabs(d[i+1] - d[i]) <
MERR) {
338 yder[i] = (w[i+1]*d[i-1] + w[i-1]*d[i])/(w[i+1] + w[i-1]);
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) {
346 }
else if (fabs(d[1] - d[0]) <
MERR) {
349 yder[0] = (w[1]*dm1 + wm1*d[0])/(w[1] + wm1);
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) {
356 }
else if (fabs(d[2] - d[1]) <
MERR) {
359 yder[1] = (w[2]*d[0] + w[0]*d[1])/(w[2] + w[0]);
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];
369 yder[nx-1] = (wnx*d[nx-2] + w[nx-2]*d[nx-1])/(wnx + w[nx-2]);
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];
379 yder[nx-2] = (w[nx-1]*d[nx-3] + w[nx-3]*d[nx-2])/(w[nx-1] + w[nx-3]);
383 for (i = 0; i <= nx-2; 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];
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) {
398 }
else if (fabs(x[i+1] - xx[ii]) <
MERR) {
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);