private import math; typedef real[] splinetype(real[], real[]); restricted real[] Spline(real[] x, real[] y); restricted splinetype[] Spline; string morepoints="interpolation requires at least 2 points"; string differentlengths="arrays have different lengths"; void checklengths(int x, int y, string text=differentlengths) { if(x != y) abort(text+": "+string(x)+" != "+string(y)); } void checkincreasing(real[] x) { if(!increasing(x,true)) abort("strictly increasing array expected"); } // Linear interpolation real[] linear(real[] x, real[] y) { int n=x.length; checklengths(n,y.length); real[] d=new real[n]; for(int i=0; i < n-1; ++i) d[i]=(y[i+1]-y[i])/(x[i+1]-x[i]); d[n-1]=d[n-2]; return d; } // Standard cubic spline interpolation with not-a-knot condition: // s'''(x_2^-)=s'''(x_2^+) et s'''(x_(n_2)^-)=s'''(x_(n-2)^+) // if n=2, linear interpolation is returned // if n=3, an interpolation polynomial of degree <= 2 is returned: // p(x_1)=y_1, p(x_2)=y_2, p(x_3)=y_3 real[] notaknot(real[] x, real[] y) { int n=x.length; checklengths(n,y.length); checkincreasing(x); real[] d; if(n > 3) { real[] a=new real[n]; real[] b=new real[n]; real[] c=new real[n]; real[] g=new real[n]; b[0]=x[2]-x[1]; c[0]=x[2]-x[0]; a[0]=0; g[0]=((x[1]-x[0])^2*(y[2]-y[1])/b[0]+b[0]*(2*b[0]+3*(x[1]-x[0]))* (y[1]-y[0])/(x[1]-x[0]))/c[0]; for(int i=1; i < n-1; ++i) { a[i]=x[i+1]-x[i]; c[i]=x[i]-x[i-1]; b[i]=2*(a[i]+c[i]); g[i]=3*(c[i]*(y[i+1]-y[i])/a[i]+a[i]*(y[i]-y[i-1])/c[i]); } c[n-1]=0; b[n-1]=x[n-2]-x[n-3]; a[n-1]=x[n-1]-x[n-3]; g[n-1]=((x[n-1]-x[n-2])^2*(y[n-2]-y[n-3])/b[n-1]+ b[n-1]*(2*b[n-1]+3(x[n-1]-x[n-2]))* (y[n-1]-y[n-2])/(x[n-1]-x[n-2]))/a[n-1]; d=tridiagonal(a,b,c,g); } else if(n == 2) { real val=(y[1]-y[0])/(x[1]-x[0]); d=new real[] {val,val}; } else if(n == 3) { real a=(y[1]-y[0])/(x[1]-x[0]); real b=(y[2]-y[1])/(x[2]-x[1]); real c=(b-a)/(x[2]-x[0]); d=new real[] {a+c*(x[0]-x[1]),a+c*(x[1]-x[0]),a+c*(2*x[2]-x[0]-x[1])}; } else abort(morepoints); return d; } // Standard cubic spline interpolation with periodic condition // s'(a)=s'(b), s''(a)=s''(b), assuming that f(a)=f(b) // if n=2, linear interpolation is returned real[] periodic(real[] x, real[] y) { int n=x.length; checklengths(n,y.length); checkincreasing(x); if(abs(y[n-1]-y[0]) > sqrtEpsilon*norm(y)) abort("function values are not periodic"); real[] d; if(n > 2) { real[] a=new real[n-1]; real[] b=new real[n-1]; real[] c=new real[n-1]; real[] g=new real[n-1]; c[0]=x[n-1]-x[n-2]; a[0]=x[1]-x[0]; b[0]=2*(a[0]+c[0]); g[0]=3*c[0]*(y[1]-y[0])/a[0]+3*a[0]*(y[n-1]-y[n-2])/c[0]; for(int i=1; i < n-1; ++i) { a[i]=x[i+1]-x[i]; c[i]=x[i]-x[i-1]; b[i]=2*(a[i]+c[i]); g[i]=3*(c[i]*(y[i+1]-y[i])/a[i]+a[i]*(y[i]-y[i-1])/c[i]); } d=tridiagonal(a,b,c,g); d.push(d[0]); } else if(n == 2) { d=new real[] {0,0}; } else abort(morepoints); return d; } // Standard cubic spline interpolation with the natural condition // s''(a)=s''(b)=0. // if n=2, linear interpolation is returned // Don't use the natural type unless the underlying function // has zero second end points derivatives. real[] natural(real[] x, real[] y) { int n=x.length; checklengths(n,y.length); checkincreasing(x); real[] d; if(n > 2) { real[] a=new real[n]; real[] b=new real[n]; real[] c=new real[n]; real[] g=new real[n]; b[0]=2*(x[1]-x[0]); c[0]=x[1]-x[0]; a[0]=0; g[0]=3*(y[1]-y[0]); for(int i=1; i < n-1; ++i) { a[i]=x[i+1]-x[i]; c[i]=x[i]-x[i-1]; b[i]=2*(a[i]+c[i]); g[i]=3*(c[i]*(y[i+1]-y[i])/a[i]+a[i]*(y[i]-y[i-1])/c[i]); } c[n-1]=0; a[n-1]=x[n-1]-x[n-2]; b[n-1]=2*a[n-1]; g[n-1]=3*(y[n-1]-y[n-2]); d=tridiagonal(a,b,c,g); } else if(n == 2) { real val=(y[1]-y[0])/(x[1]-x[0]); d=new real[] {val,val}; } else abort(morepoints); return d; } // Standard cubic spline interpolation with clamped conditions f'(a), f'(b) splinetype clamped(real slopea, real slopeb) { return new real[] (real[] x, real[] y) { int n=x.length; checklengths(n,y.length); checkincreasing(x); real[] d; if(n > 2) { real[] a=new real[n]; real[] b=new real[n]; real[] c=new real[n]; real[] g=new real[n]; b[0]=x[1]-x[0]; g[0]=b[0]*slopea; c[0]=0; a[0]=0; for(int i=1; i < n-1; ++i) { a[i]=x[i+1]-x[i]; c[i]=x[i]-x[i-1]; b[i]=2*(a[i]+c[i]); g[i]=3*(c[i]*(y[i+1]-y[i])/a[i]+a[i]*(y[i]-y[i-1])/c[i]); } c[n-1]=0; a[n-1]=0; b[n-1]=x[n-1]-x[n-2]; g[n-1]=b[n-1]*slopeb; d=tridiagonal(a,b,c,g); } else if(n == 2) { d=new real[] {slopea,slopeb}; } else abort(morepoints); return d; }; } // Piecewise Cubic Hermite Interpolating Polynomial (PCHIP) // Modified MATLAB code // [1] Fritsch, F. N. and R. E. Carlson, // "Monotone Piecewise Cubic Interpolation," // SIAM J. Numerical Analysis, Vol. 17, 1980, pp.238-246. // [2] Kahaner, David, Cleve Moler, Stephen Nash, // Numerical Methods and Software, Prentice Hall, 1988. real[] monotonic(real[] x, real[] y) { int n=x.length; checklengths(n,y.length); checkincreasing(x); real[] d=new real[n]; if(n > 2) { real[] h=new real[n-1]; real[] del=new real[n-1]; for(int i=0; i < n-1; ++i) { h[i]=x[i+1]-x[i]; del[i]=(y[i+1]-y[i])/h[i]; } int j=0; int k[]=new int[]; for(int i=0; i < n-2; ++i) if((sgn(del[i])*sgn(del[i+1])) > 0) {k[j]=i; j=j+1;} real[] hs=new real[j]; for(int i=0; i < j; ++i) hs[i]=h[k[i]]+h[k[i]+1]; real w1[]=new real[j]; real w2[]=new real[j]; real dmax[]=new real[j]; real dmin[]=new real[j]; for(int i=0; i < j; ++i) { w1[i]=(h[k[i]]+hs[i])/(3*hs[i]); w2[i]=(h[k[i]+1]+hs[i])/(3*hs[i]); dmax[i]=max(abs(del[k[i]]),abs(del[k[i]+1])); dmin[i]=min(abs(del[k[i]]),abs(del[k[i]+1])); } for(int i=0; i < n; ++i) d[i]=0; for(int i=0; i < j; ++i) d[k[i]+1]=dmin[i]/(w1[i]*(del[k[i]]/dmax[i])+w2[i]*(del[k[i]+1]/dmax[i])); d[0]=((2*h[0]+h[1])*del[0]-h[0]*del[1])/(h[0]+h[1]); if(sgn(d[0]) != sgn(del[0])) {d[0]=0;} else if((sgn(del[0]) != sgn(del[1])) && (abs(d[0]) > abs(3*del[0]))) d[0]=3*del[0]; d[n-1]=((2*h[n-2]+h[n-3])*del[n-2]-h[n-2]*del[n-2])/(h[n-2]+h[n-3]); if(sgn(d[n-1]) != sgn(del[n-2])) {d[n-1]=0;} else if((sgn(del[n-2]) != sgn(del[n-3])) && (abs(d[n-1]) > abs(3*del[n-2]))) d[n-1]=3*del[n-2]; } else if(n == 2) { d[0]=d[1]=(y[1]-y[0])/(x[1]-x[0]); } else abort(morepoints); return d; } // Return standard cubic spline interpolation as a guide guide hermite(real[] x, real[] y, splinetype splinetype=null) { int n=x.length; if(n == 0) return nullpath; guide g=(x[0],y[0]); if(n == 1) return g; if(n == 2) return g--(x[1],y[1]); if(splinetype == null) splinetype=(x[0] == x[x.length-1] && y[0] == y[y.length-1]) ? periodic : notaknot; real[] dy=splinetype(x,y); for(int i=1; i < n; ++i) { pair z=(x[i],y[i]); real dx=x[i]-x[i-1]; g=g..controls((x[i-1],y[i-1])+dx*(1,dy[i-1])/3) and (z-dx*(1,dy[i])/3)..z; } return g; }