Improve "fill" pressure estimate

diff --git a/get_important_pressure.m b/get_important_pressure.m
index 2093515a089b585cdaad45067e16bf42fb9e4bf9..a75a82a09d0ea94a84e98cb0f6dae8a31b0d4d52 100644 (file)
--- a/get_important_pressure.m
+++ b/get_important_pressure.m
@@ -13,7 +13,7 @@ function [p_eq, p_after_fill] = get_important_pressure(t, p,  graph = 1, thresho
                p_single = p(open_inds(i - 1):(fill_inds(i) - 1));
                p_eq(i - 1)  = get_equilibrium_pressure(t_single, p_single);
                if (i > 2)
-                       p_after_fill(i - 1) = get_equilibrium_pressure(t(fill_inds(i - 1):open_inds(i - 1) - 1), p(fill_inds(i - 1):open_inds(i - 1) - 1));
+                       p_after_fill(i - 1) = get_fill_pressure(p(fill_inds(i - 1):open_inds(i - 1) - 1));
                end
                if (graph)
                        if (i < length(fill_inds))
@@ -33,3 +33,19 @@ function p_eq = get_equilibrium_pressure(t, p)
        weights = ((t - t(1)) / (t(end) - t(1))) .^ 10;
        p_eq = weighted_mean(p, weights);
 end
+
+function p_fill = get_fill_pressure(p)
+       p_fill = mean(p);
+       while (numel(p) > 2)
+               p_mean = mean(p);
+               res_sq = (p - p_mean) .^ 2;
+               p_var = mean(res_sq);
+               remove_inds = res_sq > 3 * p_var;
+               if (all(~remove_inds))
+                       break ;
+               end
+               [~, remove_inds] = max(res_sq);
+               p(remove_inds) = [];
+       end
+       p_fill = p_mean;
+end