[prev in list] [next in list] [prev in thread] [next in thread] 

List:       kde-commits
Subject:    [labplot/analysis_interpolation] src: added rational function interpolation
From:       Stefan Gerlach <stefan.gerlach () uni-konstanz ! de>
Date:       2016-04-30 18:29:59
Message-ID: E1awZeR-0006S1-OE () scm ! kde ! org
[Download RAW message or body]

Git commit 4cc356dfc00c1ca25e31622d4ba1c9f65cf1a168 by Stefan Gerlach.
Committed on 30/04/2016 at 18:29.
Pushed by sgerlach into branch 'analysis_interpolation'.

added rational function interpolation

M  +56   -5    src/backend/worksheet/plots/cartesian/XYInterpolationCurve.cpp
M  +1    -1    src/backend/worksheet/plots/cartesian/XYInterpolationCurve.h
M  +1    -0    src/backend/worksheet/plots/cartesian/XYInterpolationCurvePrivate.h
M  +1    -1    src/kdefrontend/dockwidgets/XYInterpolationCurveDock.cpp

http://commits.kde.org/labplot/4cc356dfc00c1ca25e31622d4ba1c9f65cf1a168

diff --git a/src/backend/worksheet/plots/cartesian/XYInterpolationCurve.cpp \
b/src/backend/worksheet/plots/cartesian/XYInterpolationCurve.cpp index \
                0c855b3..3567552 100644
--- a/src/backend/worksheet/plots/cartesian/XYInterpolationCurve.cpp
+++ b/src/backend/worksheet/plots/cartesian/XYInterpolationCurve.cpp
@@ -41,6 +41,7 @@
 #include "backend/lib/commandtemplates.h"
 
 #include <cmath>	// isnan
+#include <values.h>	// DBL_MIN
 #include <gsl_errno.h>
 #include <gsl/gsl_interp.h>
 #include <gsl/gsl_spline.h>
@@ -239,6 +240,53 @@ void XYInterpolationCurvePrivate::integ(double *x, double *y, \
unsigned n) {  y[n-1]=y[n-2]+vold;
 }
 
+// calculates rational interpolation of n points of xy-data at xn using \
Burlisch-Stoer method. result in v (error dv) +void \
XYInterpolationCurvePrivate::ratint(double *x, double *y, int n, double xn, double \
*v, double *dv) { +	int i,j,a=0,b=n-1;
+	while(b-a>1) {	// find interval using bisection
+		j=floor((a+b)/2.);
+		if(x[j]>xn)
+			b=j;
+		else
+			a=j;
+	}
+
+	int ns=a;// nearest index
+	if(fabs(xn-x[a])>fabs(xn-x[b]))
+		ns=b;
+
+	if(xn==x[ns]) {	// exact point
+		*v=y[ns];
+		*dv=0;
+		return;
+	}
+
+	double *c = (double*)malloc(n*sizeof(double));
+	double *d = (double*)malloc(n*sizeof(double));
+	for(i=0;i<n;i++)
+		c[i]=d[i]=y[i];
+
+	*v=y[ns--];
+
+	double t,dd;
+	for(int m=1;m<n;m++) {
+		for(i=0;i<n-m;i++) {
+			t=(x[i]-xn)*d[i]/(x[i+m]-xn);
+			dd=t-c[i+1];
+			if(dd==0.0) // pole
+				dd+=DBL_MIN;
+			dd=(c[i+1]-d[i])/dd;
+			d[i]=c[i+1]*dd;
+			c[i]=t*dd;
+		}
+
+		*dv = (2*(ns+1) < (n-m) ? c[ns+1] : d[ns--]);
+		*v += *dv;
+	}
+
+	free(c);free(d);
+}
+
 void XYInterpolationCurvePrivate::recalculate() {
 	QElapsedTimer timer;
 	timer.start();
@@ -427,10 +475,6 @@ void XYInterpolationCurvePrivate::recalculate() {
 			t = (x-xdata[a])/(xdata[b]-xdata[a]);
 			(*yVector)[i] = ydata[a]*pow(ydata[b]/ydata[a],t);
 			break;
-		case XYInterpolationCurve::Rational:
-			//TODO
-			//(*yVector)[i] = 
-			break;
 		case XYInterpolationCurve::PCH: {
 			t = (x-xdata[a])/(xdata[b]-xdata[a]);
 			double h1=2.*t*t*t-3.*t*t+1, h2=-2.*t*t*t+3.*t*t, h3=t*t*t-2*t*t+t, h4=t*t*t-t*t;
@@ -492,6 +536,13 @@ void XYInterpolationCurvePrivate::recalculate() {
 			(*yVector)[i] = ydata[a]*h1+ydata[b]*h2+(xdata[b]-xdata[a])*(m1*h3+m2*h4);
 		}
 			break;
+		case XYInterpolationCurve::Rational: {
+			double v,dv;
+			ratint(xdata, ydata, n, x, &v, &dv);
+			(*yVector)[i] = v;
+			//TODO: use dv
+		}
+			break;
 		}
 	}
 
@@ -499,8 +550,8 @@ void XYInterpolationCurvePrivate::recalculate() {
 	switch(type) {
 	case XYInterpolationCurve::Cosine:
 	case XYInterpolationCurve::Exponential:
-	case XYInterpolationCurve::Rational:
 	case XYInterpolationCurve::PCH:
+	case XYInterpolationCurve::Rational:
 		switch(evaluate) {
 		case XYInterpolationCurve::Function:
 			break;
diff --git a/src/backend/worksheet/plots/cartesian/XYInterpolationCurve.h \
b/src/backend/worksheet/plots/cartesian/XYInterpolationCurve.h index 6d77b76..cf3fb08 \
                100644
--- a/src/backend/worksheet/plots/cartesian/XYInterpolationCurve.h
+++ b/src/backend/worksheet/plots/cartesian/XYInterpolationCurve.h
@@ -37,7 +37,7 @@ class XYInterpolationCurve: public XYCurve {
 	Q_OBJECT
 
 	public:
-		enum InterpolationType \
{Linear,Polynomial,CSpline,CSplinePeriodic,Akima,AkimaPeriodic,Steffen,Cosine,Exponential,Rational,PCH};
 +		enum InterpolationType \
{Linear,Polynomial,CSpline,CSplinePeriodic,Akima,AkimaPeriodic,Steffen,Cosine,Exponential,PCH,Rational};
  enum CubicHermiteVariant {FiniteDifference,CatmullRom,Cardinal,KochanekBartels};
 		enum InterpolationEval {Function,Derivative,Derivative2,Integral};
 
diff --git a/src/backend/worksheet/plots/cartesian/XYInterpolationCurvePrivate.h \
b/src/backend/worksheet/plots/cartesian/XYInterpolationCurvePrivate.h index \
                d2f6f2d..4fba513 100644
--- a/src/backend/worksheet/plots/cartesian/XYInterpolationCurvePrivate.h
+++ b/src/backend/worksheet/plots/cartesian/XYInterpolationCurvePrivate.h
@@ -65,6 +65,7 @@ class XYInterpolationCurvePrivate: public XYCurvePrivate {
 		void deriv(double *x, double *y, unsigned n);
 		void deriv2(double *x, double *y, unsigned n);
 		void integ(double *x, double *y, unsigned n);
+		void ratint(double *x, double *y, int n, double xn, double *v, double *dv);
 };
 
 #endif
diff --git a/src/kdefrontend/dockwidgets/XYInterpolationCurveDock.cpp \
b/src/kdefrontend/dockwidgets/XYInterpolationCurveDock.cpp index 416a51f..90344d7 \
                100644
--- a/src/kdefrontend/dockwidgets/XYInterpolationCurveDock.cpp
+++ b/src/kdefrontend/dockwidgets/XYInterpolationCurveDock.cpp
@@ -93,8 +93,8 @@ void XYInterpolationCurveDock::setupGeneral() {
 #endif
 	uiGeneralTab.cbType->addItem(i18n("cosine"));
 	uiGeneralTab.cbType->addItem(i18n("exponential"));
-	uiGeneralTab.cbType->addItem(i18n("rational functions"));
 	uiGeneralTab.cbType->addItem(i18n("piecewise cubic Hermite (PCH)"));
+	uiGeneralTab.cbType->addItem(i18n("rational functions"));
 
 	uiGeneralTab.cbVariant->addItem(i18n("finite differences"));
 	uiGeneralTab.cbVariant->addItem(i18n("Catmull-Rom"));


[prev in list] [next in list] [prev in thread] [next in thread]