Implement the get_important_pressure function

This function reduces the input data to only the starting and
equilibrium pressures of each "cycle".

diff --git a/get_important_pressure.m b/get_important_pressure.m
new file mode 100644 (file)
index 0000000..7e546f2
--- /dev/null
+++ b/get_important_pressure.m
@@ -0,0 +1,31 @@
+function [p_eq, p_after_open] = get_important_pressure(t, p,  graph = 1, threshold = 0.05)
+       log_inds = detect_open(t, p, threshold);
+       inds = [1 find(log_inds == 1) length(p)];
+       if (graph)
+               figure()
+               plot(t, p);
+               hold on;
+       end
+       for i = 2:length(inds)
+               t_single = t(inds(i - 1):(inds(i) - 1));
+               p_single = p(inds(i - 1):(inds(i) - 1));
+               [p_eq(i - 1) p_after_open(i - 1)] = process_single(t_single, p_single);
+               if (graph)
+                       if (i < length(inds))
+                               plot([t(inds(i) - 5) t(inds(i) + 4)], p_eq(i - 1) .* [1 1], "red");
+                       else
+                               plot([t(inds(i) - 5) t(end)], p_eq(i - 1) .* [1 1], "red");
+                       end
+                       if (i > 2)
+                               plot([t(inds(i - 1) - 4) t(inds(i - 1) + 5)], p_after_open(i - 1) .* [1 1], "green");
+                       end
+               end
+       end
+       p_after_open(1) = [];
+end
+
+function [p_eq p_after_open] = process_single(t, p)
+       p_after_open = max(p);
+       weights = ((t - t(1)) / (t(end) - t(1))) .^ 10;
+       p_eq = weighted_mean(p, weights);
+end

diff --git a/pseudo_exp.m b/pseudo_exp.m
index 64357629ba26bf9b4934eb2a7d6c23ae92a86b30..8d5a40a6fa2e1fe31d25285e2735d3a7b9fa7ae8 100644 (file)
--- a/pseudo_exp.m
+++ b/pseudo_exp.m
@@ -7,7 +7,7 @@ function [t, p, t_open] = pseudo_exp(p1, V1, V2, Vvz, n)
        for i = ceil(t_open)
                p2 = (p(i + 1) * (V2 - Vvz) + p1 * V1) / (V1 + V2 - Vvz);
                peq += 0.7 * (p2 - peq);
-               p(i + 1:end) = (p2 - peq) * exp((-(i + 1:length(p)) + i + 1) / 10) + peq;
+               p(i + 1:end) = (p2 - peq) * exp((-((i + 1):length(p)) + i + 1) / 10) + peq;
        end
        p += rand(size(t)) * p1 / 100;
 end

diff --git a/weighted_mean.m b/weighted_mean.m
new file mode 100644 (file)
index 0000000..646ec0e
--- /dev/null
+++ b/weighted_mean.m
@@ -0,0 +1,4 @@
+function res = weighted_mean(values, weights)
+       weighted_sum = sum(values .* weights);
+       res = weighted_sum / sum(weights);
+end