From 9555b838910af32dea517c837f5c533ef2ac61ba Mon Sep 17 00:00:00 2001
From: =?utf8?q?Luk=C3=A1=C5=A1=20Ji=C5=99i=C5=A1t=C4=9B?=
 <jiriste@icpf.cas.cz>
Date: Tue, 14 May 2024 14:27:25 +0200
Subject: [PATCH] Implement the get_important_pressure function

This function reduces the input data to only the starting and
equilibrium pressures of each "cycle".
---
 get_important_pressure.m | 31 +++++++++++++++++++++++++++++++
 pseudo_exp.m             |  2 +-
 weighted_mean.m          |  4 ++++
 3 files changed, 36 insertions(+), 1 deletion(-)
 create mode 100644 get_important_pressure.m
 create mode 100644 weighted_mean.m

diff --git a/get_important_pressure.m b/get_important_pressure.m
new file mode 100644
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 6435762..8d5a40a 100644
--- 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
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
-- 
2.30.2