From 42863b8e03208e6fbba1c2e5d732dad130eafe35 Mon Sep 17 00:00:00 2001
From: "Joachim Wuttke (l)" <j.wuttke@fz-juelich.de>
Date: Tue, 22 Sep 2015 09:32:33 +0200
Subject: [PATCH] correct intpol for extrapolation by 0
---
pub/lib/node.cpp | 4 ++--
pub/lib/slice.cpp | 17 ++++++++++++-----
2 files changed, 14 insertions(+), 7 deletions(-)
diff --git a/pub/lib/node.cpp b/pub/lib/node.cpp
index 8ad5ed39..d54f5d35 100644
--- a/pub/lib/node.cpp
+++ b/pub/lib/node.cpp
@@ -1372,8 +1372,8 @@ RObjVecDbl CNodeDirac::convolve( const CContext& ctx, double theshift, const PSp
PObjVecDbl ret;
PObjVecEnu dret;
if ( ctx.want_error ) {
- dret = PObjVecEnu( new CObjVecEnu(n) );
- ret = dret;
+ dret = PObjVecEnu( new CObjVecEnu(n) );
+ ret = dret;
} else {
ret = PObjVecDbl( new CObjVecDbl(n) );
dret = NULL;
diff --git a/pub/lib/slice.cpp b/pub/lib/slice.cpp
index fda72fad..6db9a301 100644
--- a/pub/lib/slice.cpp
+++ b/pub/lib/slice.cpp
@@ -188,6 +188,11 @@ void CSpec::intpol( const vector<double>& xx, vector<double>* yy, vector<double>
(*dd)[ii] = 0;
}
int i = 0, ii = 0;
+ while( xx[ii]<x[0] ) {
+ ++ii;
+ if( ii>=nn )
+ return;
+ }
while ( i<n-1 && xx[ii]<x[n-1] ) {
while ( xx[ii]>x[i+1] && i<n-1 )
++i;
@@ -196,12 +201,14 @@ void CSpec::intpol( const vector<double>& xx, vector<double>* yy, vector<double>
if( ii>=nn )
return;
(*yy)[ii] = ( (x[i+1]-xx[ii])*y[i] + (xx[ii]-x[i])*y[i+1] )/(x[i+1]-x[i]);
- if( !isfinite((*yy)[ii]) ) {
- std::cout << "TROUBLE " << i << " " << ii << "\n";
+ //std::cout << "DEBUG INTPOL i=" << i << " ii=" << ii << " x=" << x[i] << ".." << x[i+1] <<
+ // " xx=" << xx[ii] << " y=" << y[i] << ".." << y[i+1] << " -> " << (*yy)[ii] << "\n";
+ if( dd ) {
+ double var = ( (x[i+1]-xx[ii])*SQR(dy[i]) + (xx[ii]-x[i])*SQR(dy[i+1]) ) / (x[i+1]-x[i]);
+ if ( var< 0 )
+ throw "interpolation failed: variance<0";
+ (*dd)[ii] = sqrt( var );
}
- if( dd )
- (*dd)[ii] = sqrt( (x[i+1]-xx[ii])*SQR(dy[i]) + (xx[ii]-x[i])*SQR(dy[i+1]) ) /
- (x[i+1]-x[i]);
++ii;
}
}
--
GitLab