github-starred/ultimatepp
2
0
Fork 0
mirror of https://github.com/ultimatepp/ultimatepp.git synced 2026-06-20 06:05:32 -06:00
Code Issues 33 Projects Packages Wiki Activity Actions

STEM4U: Simplified integral calculation for tine data sets

Browse source 
git-svn-id: svn://ultimatepp.org/upp/trunk@15332 f0d560ea-af0d-0410-9eb7-867de7ffcac7
This commit is contained in:
koldo 2020-10-30 13:04:13 +00:00
parent 86e4cd9fdc
commit d115c62e3a
1 changed files with 28 additions and 13 deletions
Download patch file Download diff file Expand all files Collapse all files

41
bazaar/STEM4U/Integral.h
View file

@ -8,24 +8,30 @@ enum IntegralType {TRAPEZOIDAL, SIMPSON_1_3, SIMPSON_3_8, HERMITE_3, HERMITE_5};
template <class Range, class T>
T Integral(const Range &y, const Range &x, IntegralType type = TRAPEZOIDAL) {
ASSERT(x.size() == y.size());
ASSERT(x.size() > 1);
if (y.size() <= 1)
return 0;
T ret = 0;
size_t n = x.size();
if(type == SIMPSON_1_3 && n == 2)
type = TRAPEZOIDAL;
else if (type == SIMPSON_3_8) {
if (n == 2)
type = TRAPEZOIDAL;
else if (n == 3)
type = SIMPSON_1_3;
}
if (type == TRAPEZOIDAL) {
for (int i = 1; i < n; ++i)
ret += Avg(y(i), y(i-1))*(x(i) - x(i-1));
} else if (type == SIMPSON_1_3) {
if (n < 3)
return Null;
int i;
for (i = 2; i < n; i += 2)
ret += (x[i] - x[i-2])/6.*(y[i-2] + 4*y[i-1] + y[i]);
if (i == n)
ret += Avg(y(n-1), y(n-2))*(x(n-1) - x(n-2));
} else if (type == SIMPSON_3_8) {
if (n < 4)
return Null;
int i;
for (i = 3; i < n; i += 3)
ret += (x[i] - x[i-3])/8.*(y[i-3] + 3*y[i-2] + 3*y[i-1] + y[i]);
@ -61,12 +67,25 @@ T Integral(Range &y, T dx, IntegralType type = TRAPEZOIDAL) {
return 0;
Eigen::Index n = y.size();
if(type == SIMPSON_1_3 && n == 2)
type = TRAPEZOIDAL;
else if (type == SIMPSON_3_8) {
if (n == 2)
type = TRAPEZOIDAL;
else if (n == 3)
type = SIMPSON_1_3;
} else if (type == HERMITE_3 && n == 2)
type = TRAPEZOIDAL;
else if (type == HERMITE_5) {
if (n == 2)
type = TRAPEZOIDAL;
else if (n == 3)
type = HERMITE_3;
}
if (type == TRAPEZOIDAL)
return (y.segment(1, n-2).sum() + (y(0) + y(n-1))/2)*dx;
else if (type == SIMPSON_1_3) {
if (n < 3)
return Null;
T ret0 = 0;
if ((n-1)%2) {
ret0 = Avg(y(n-1), y(n-2))*dx;
@ -74,8 +93,6 @@ T Integral(Range &y, T dx, IntegralType type = TRAPEZOIDAL) {
}
return ret0 + Calc1_3(y, dx, n);
} else if (type == SIMPSON_3_8) {
if (n < 4)
return Null;
T ret0 = 0;
int rem = (n-1)%3;
if (rem == 2) {
@ -94,13 +111,11 @@ T Integral(Range &y, T dx, IntegralType type = TRAPEZOIDAL) {
}
return ret0 + ret*dx*3./8;
} else if (type == HERMITE_3) {
if (n < 3)
return Null;
return (y.segment(1, n-2).sum() + (y(0) + y(n-1))/2)*dx
+ dx/24.*(3*y[0] - 4*y[1] + y[2] + y[n-3] - 4*y[n-2] + 3*y[n-1]);
} else if (type == HERMITE_5) {
if (n < 5)
return Null;
if (n == 4)
return ((y[0] + y[3])/2. + y[1] + y[2])*dx;
return (y.segment(1, n-2).sum() + (y(0) + y(n-1))/2)*dx
- dx/144.*(25*(y[0] - y[n-1]) - 48*(y[1] - y[n-2]) + 36*(y[2] - y[n-3])
- 16*(y[3] - y[n-4]) + 3*(y[4] - y[n-5]));