diff options
| -rw-r--r-- | hadamard_pt.hpp | 29 |
1 files changed, 24 insertions, 5 deletions
diff --git a/hadamard_pt.hpp b/hadamard_pt.hpp index a57c670..90fa661 100644 --- a/hadamard_pt.hpp +++ b/hadamard_pt.hpp @@ -85,11 +85,13 @@ public: } std::vector<double> f(Ms.size()); + std::vector<double> δf(Ms.size()); std::vector<unsigned> f_keep; double f_last = 0; for (unsigned i = 0; i < Ms.size(); i++) { f[i] = nu[i] / (double)(nu[i] + nd[i]); + δf[i] = f[i] * sqrt(1.0/nu[i] + 1.0/(nu[i]+nd[i])); if (f[i] >= f_last) { f_keep.push_back(i); f_last = f[i]; @@ -101,9 +103,21 @@ public: f[j] = f[f_keep[i]] + (f[f_keep[i + 1]] - f[f_keep[i]]) / (Ms[f_keep[i + 1]].β - Ms[f_keep[i]].β) * (Ms[j].β - Ms[f_keep[i]].β); + δf[j] = sqrt(pow(δf[f_keep[i]], 2) + pow(sqrt(pow(δf[f_keep[i + 1]], 2) + pow(δf[f_keep[i]], 2)) / + (Ms[f_keep[i + 1]].β - Ms[f_keep[i]].β) * + (Ms[j].β - Ms[f_keep[i]].β), 2)); } } + // measure the difference between f and a straight line + double Δf² = 0; + double δΔf⁴ = 0; + + for (unsigned j = 1; j < f.size(); j++) { + Δf² += pow(f[j] - (j + 1.0) / (Ms.size() + 1.0), 2); + δΔf⁴ += pow(2 * δf[j] * f[j], 2); + } + std::vector<double> η(Ms.size() - 1); for (unsigned i = 0; i < Ms.size() - 1; i++) { @@ -132,18 +146,23 @@ public: double err = 0; for (unsigned i = 0; i < T1.size(); i++) { err += fabs(T1[i] - 1 / Ms[Ms.size() - i - 2].β); - Ms[Ms.size() - i - 2].β = 1 / T1[i]; + Ms[Ms.size() - i - 2].β = 1.0 / T1[i]; } - double relErr = err / T1.size() * Ms.size() / (1 / Ms.front().β - 1 / Ms.back().β); + double Δf = sqrt(Δf²); + double δΔf = sqrt(δΔf⁴) / 2 / Δf; + double relErr = Δf / Ms.size(); + // double relErr = err / T1.size() * Ms.size() / (1 / Ms.front().β - 1 / Ms.back().β); if (relErr < ε2) { return f; } else { - std::cout << "Difference from linear in transit flow: " << relErr << ", continuing tuning.\n"; + std::cout << "RMS difference from ideal transit flow is " << relErr << " ± " << δΔf / Ms.size() << ", continuing tuning.\n"; + if (5 * δΔf > Δf || δΔf != δΔf) { + n *= 2; + std::cout << "Error in RMS difference too close to difference, increasing tuning to " << n << ".\n"; + } } - - n += n0; } } |