--- /dev/null
+clear all, close all\r
+i=1;\r
+mezvak=0.07;\r
+cit=50;\r
+k=50;\r
+\r
+Vb=36.81;\r
+KK=616.26;%612.47\r
+Va=16.26;\r
+Vv=33.87;\r
+Vven=Va+Vb+KK;\r
+\r
+[fnamecelk,path]=uigetfile({'*.txt';'*.*'},'File Selector','MultiSelect','on');\r
+fnamecelk=cellstr(fnamecelk);\r
+\r
+\r
+for i=1:length(fnamecelk)\r
+ dataV{i}=importdata([path '/' fnamecelk{i}], '\t', 2);\r
+ mn=menu(['Jaký objem je souborem ' fnamecelk{i} ' kalibrován?'],'Horní(V1)','Spodní(V2)');\r
+ P{i,1}=dataV{i}.data(:,4-mn)';\r
+ P{i,4}=dataV{i}.data(:,mn+1)';\r
+%hledání skoků\r
+ vaku=mean(P{i,1}(P{i,1}<0.01));\r
+ P{i,4}(P{i,1}<mezvak)=[];%smazaní evakuace\r
+ P{i,1}(P{i,1}<mezvak)=[];\r
+ P{i,2}=-P{i,1}(2:end)+P{i,1}(1:end-1); %derivace odhaluje skoky\r
+ P{i,3}=P{i,2}./(P{i,1}(1:end-1)/max(P{i,1})); %nomralizace skoku\r
+ minind=find(P{i,3}<min(P{i,3})/cit);%hledání napouštění\r
+ %smazání vakua a napouštění\r
+ P{i,3}(P{i,3}>max(P{i,3})/cit)=1;\r
+ P{i,3}(P{i,3}<min(P{i,3})/cit)=-1;\r
+ extind=find(abs(P{i,3})==1);%hledání skoků\r
+% maxind=find(P{i,3}==1)\r
+% minind=find(P{i,3}==-1)\r
+ %skoky do 1 bodu\r
+ smaz=[];\r
+ for l=flip((2:length(extind)))\r
+ if (extind(l)-extind(l-1)<k)&&(P{i,3}(extind(l))>0)\r
+ smaz=[smaz extind(l)];%zapisuje "duální" hodnoty\r
+ end\r
+ end\r
+ for l=smaz\r
+ extind(extind==l)=[];\r
+ extind(extind>l)=extind(extind>l)-1;\r
+ P{i,1}(l)=[]; \r
+ P{i,4}(l)=[];\r
+ P{i,2}(l)=[];\r
+ P{i,3}(l)=[];\r
+ end\r
+ %ještě jednou z druhé strany\r
+ smaz=[];\r
+ for l=((1:length(extind)-1))\r
+ if (extind(l+1)-extind(l)<k)&&(P{i,3}(extind(l))>0)\r
+ smaz=[smaz extind(l)];%zapisuje "duální" hodnoty\r
+ end\r
+ end\r
+ for l=flip(smaz)\r
+ extind(extind==l)=[];\r
+ extind(extind>l)=extind(extind>l)-1;\r
+ P{i,1}(l)=[]; \r
+ P{i,2}(l)=[];\r
+ P{i,3}(l)=[];\r
+ P{i,4}(l)=[];\r
+ end\r
+ %a ještě jednou jen pro nárůst tlaku\r
+ smaz=[];\r
+ minind=extind(P{i,3}(extind)==-1);\r
+ for l=flip((1:length(minind)-1))\r
+ if minind(l+1)-minind(l)<k %100\r
+ smaz=[smaz minind(l)];%zapisuje "duální" hodnoty\r
+ end\r
+ end \r
+ \r
+ for l=smaz\r
+ extind(extind==l)=[];\r
+ extind(extind>l)=extind(extind>l)-1;\r
+ P{i,1}(l)=[]; \r
+ P{i,2}(l)=[];\r
+ P{i,3}(l)=[];\r
+ P{i,4}(l)=[];\r
+ end\r
+ P{i,1}(end-10:end)=[];\r
+ P{i,2}(end-10:end)=[];\r
+ P{i,3}(end-10:end)=[];\r
+ P{i,4}(end-10:end)=[];\r
+ \r
+ \r
+ figure(2)\r
+ plot(P{i,3})\r
+ figure(1)\r
+ while extind(end)>length(P{i,3})\r
+ extind(end)=[];\r
+ end\r
+ if extind(1)==1\r
+ P{i,1}(1:extind(1))=[];\r
+ P{i,2}(1:extind(1))=[];\r
+ P{i,3}(1:extind(1))=[];\r
+ P{i,4}(1:extind(1))=[];\r
+ extind=extind-extind(1);\r
+ extind(1)=[];\r
+ end\r
+% extind=find(abs(P{i,3})==1);\r
+ \r
+ %Ruční kontrola\r
+ kontr=1;\r
+ while kontr~=0\r
+ plot(P{i,1},'blue')\r
+ hold on\r
+ for l=extind\r
+ if P{i,3}(l)==1\r
+ plot([l l],[0 max(P{i,1})],'--green')\r
+ elseif P{i,3}(l)==-1\r
+ plot([l l],[0 max(P{i,1})],'--red')\r
+ else\r
+ error(['vector extind missaligned at position ' int2str(l)])\r
+ end\r
+ end\r
+ kontr=menu('Byla detekce bezchybná?','Ano','Ne, chci odebrat čáru','Ne, chci přebarvit čáru', 'Ne, chci přidat červenou čáru')-1;\r
+ if kontr==1\r
+ %msgbox('Klikni na čáry označující špatně, potvrď enterem.')\r
+ [xspa,~]=ginput;\r
+ for l=xspa'\r
+ [~,ind]=min(abs(l-extind));\r
+ extind(ind)=[];\r
+ end\r
+ elseif kontr==2\r
+ [xspa,~]=ginput;\r
+ for l=xspa'\r
+ [~,ind]=min(abs(l-extind));\r
+ P{i,3}(extind(ind))=-1*P{i,3}(extind(ind));\r
+ end\r
+ elseif kontr==3\r
+ [xspa,~]=ginput;\r
+ for l=xspa'\r
+ \r
+ end\r
+ end\r
+ hold off\r
+ end \r
+ close all\r
+ vecdel=extind-[0 extind(1:end-1)];\r
+ Y=mat2cell(P{i,1},1,[vecdel length(P{i,1})-extind(end)]);\r
+ \r
+ plot(P{i,1},'blue')\r
+ hold on\r
+ %třídění vektorů\r
+ extind=[0 extind length(P{i,1})];\r
+ for j=1:length(Y)\r
+ iter=0;\r
+ X=(extind(j)+1:extind(j+1));\r
+ Hran{i}(j,:)=[X(1) X(end)];\r
+ suctv=[1 -1];\r
+ while max(suctv)>2*mean(suctv)\r
+ iter=iter+1;\r
+ Koef{i}(j,:)=polyfit(X,Y{j},1);\r
+ Yodh=polyval(Koef{i}(j,:),X);\r
+ suctv=(Y{j}-Yodh).^2;\r
+ X(suctv>2*mean(suctv))=[];\r
+ Y{j}(suctv>2*mean(suctv))=[];\r
+ end\r
+ plot(Hran{i}(j,:),Koef{i}(j,1)*Hran{i}(j,:)+Koef{i}(j,2))\r
+ end\r
+ hold off\r
+ Phran{i}=[];\r
+ for j=1:length(Y)-1\r
+ if P{i,3}(Hran{i}(j,2))==1\r
+ p{i,1}(j)=Koef{i}(j,1)*Hran{i}(j,2)+Koef{i}(j,2);\r
+ p{i,2}(j)=Koef{i}(j+1,1)*Hran{i}(j+1,1)+Koef{i}(j+1,2);\r
+ Phran{i}=[Phran{i} [p{i,1}(j); p{i,2}(j)]];\r
+ end\r
+ end\r
+ \r
+ Vznak(i)=mn;\r
+ V{i}=((Vv+Vven)*Phran{i}(2,:)-Vv*Phran{i}(1,:)-Vven*vaku)./(Phran{i}(1,:)-Phran{i}(2,:));\r
+ extind(1)=[];\r
+ extind(end)=[];\r
+ extind(P{i,3}(extind)~=1)=[];\r
+ P{i,5}=abs(P{i,4}-P{i,1});\r
+ delP{i}=P{i,5}(extind);\r
+ \r
+end\r
+for i=1:length(V)\r
+ Vprum(i)=mean(V{i});\r
+ printf("Výsledek z měření %s:\n\tV%i = %f mm^3\n\t", fnamecelk{i}, Vznak(i), Vprum(i));\r
+end\r
+V1=mean(Vprum(Vznak==1))\r
+V2=mean(Vprum(Vznak==2))\r
+uA1=(sum((V1-[V{Vznak==1}]).^2)/(length([V{Vznak==1}])-1))^0.5\r
+uA2=(sum((V2-[V{Vznak==2}]).^2)/(length([V{Vznak==2}])-1))^0.5\r
--- /dev/null
+function [val, ind] = allmin(array)
+ dim = sum(size(array) > 1);
+ if dim == 0
+ val = array;
+ ind = 1;
+ elseif dim == 1
+ [val ,ind] = min(array);
+ else
+ ind = [];
+ [small_arr, single_ind] = min(array, [], dim);
+ [val, other_ind] = allmin(small_arr);
+ ind = single_ind(num2cell(other_ind){:});
+ ind = [other_ind ind];
+ end
+end
--- /dev/null
+function [x, hist] = alternating_min(objfun, x_0, tol = 0.01, max_iter = 10, max_subiter = 2)
+ input_shape = size(x_0);
+ arity = prod(input_shape);
+ mask = zeros(input_shape);
+ x = x_0;
+ iter = 0;
+ rel_change = tol + 1;
+ while (iter < max_iter && rel_change > tol)
+ rel_change = 0;
+ iter += 1;
+ for i = 1:arity
+ mask(i) = 1;
+ new_fun = @(single_par) objfun(mask * single_par + ~mask .* x);
+ old = x(i);
+ x(i) = steepest_grad(new_fun, x(i), tol, max_subiter);
+ mask(i) = 0;
+ if (abs((old - x(i))/old) > rel_change)
+ rel_change = abs((old - x(i))/old);
+ end
+ end
+ end
+end
--- /dev/null
+function x = altminimize(objfun, x_0, tol = 0.01, max_iter = 2)
+
+ s = size(x_0);
+ if s(1) ~= 1 && s(2) ~= 1
+ printf("x_0 must be given as an array of parameters!\n");
+ return ;
+ elseif (s(2) ~= 1)
+ x = x_0';
+ else
+ x = x_0;
+ end
+
+ Y = zeros(size(x));
+ change = (tol + 1) * ones(size(x));
+ while (max(abs(change)) > tol)
+ for i=1:length(x)
+ y0 = x(i);
+ x(i) = 0;
+ Y(i) = 1;
+ x(i) = minimize(@(y) objfun(x + y*Y), y0, tol);
+ change(i) = (y0 - x(i)) / min(y0, x(i))
+ Y(i) = 0;
+ end
+ end
+end
--- /dev/null
+function [t, x, p] = basic_model(pars)
+% Unpacking x_0
+ tstep = pars(1);
+ xstep = pars(2);
+ tstart = pars(3);
+ tend = pars(4);
+ D = pars(5);
+ P = pars(6);
+ l = pars(7);
+ d = pars(8);
+ V0 = pars(9);
+ V1 = pars(10);
+ V2 = pars(11);
+ pmax = pars(12);
+
+ if length(pars)<16
+ all = false;
+ else
+ all = pars(16);
+ end
+
+ if length(pars)<15
+ T = 298;
+ else
+ T = pars(15);
+ end
+
+ if length(pars)<14
+ pvac = 0;
+ else
+ pvac = pars(14);
+ end
+
+ if length(pars)<13
+ tu = 0.5;
+ else
+ tu = pars(13);
+ end
+
+% Defining helpful constants
+ R = 8.314;
+ A = pi*d^2/4;
+ K0 = P*R*T*A*tstep/(V1+V0)/2/xstep;
+ K1 = P*R*T*A*tstep/V1/2/xstep;
+ K2 = P*R*T*A*tstep/V2/2/xstep;
+ alpha = D*tstep/xstep^2;
+
+ m1 = tu/tstep+1;
+ m2 = floor((tend - tu)/tstep);
+ n = l/xstep+1;
+ m = m1+m2;
+
+% Initialization of resulting matrix
+ p = zeros(n,m);
+ p(1,1) = pmax;
+
+% Defining matrices of partial differential equations
+ A1 = zeros(n);
+ A1(2:n+1:n*n-2*n-1) = -alpha;
+ A1(n+2:n+1:n*n-n-1) = 1+2*alpha;
+ A1(2*n+2:n+1:n*n-1) = -alpha;
+ A1(1, 1:3) = [1+3*K0 -4*K0 K0];
+ A1(end, end-2:end) = [K2 -4*K2 1+3*K2];
+
+ A2 = zeros(n);
+ A2(2:n+1:n*n-2*n-1) = -alpha;
+ A2(n+2:n+1:n*n-n-1) = 1+2*alpha;
+ A2(2*n+2:n+1:n*n-1) = -alpha;
+ A2(1, 1:3) = [1+3*K1 -4*K1 K1];
+ A2(end, end-2:end) = [K2 -4*K2 1+3*K2];
+
+% Numeric solution of Fick's law
+ for i = 2:m1
+ p(:,i) = A1\p(:,i-1);
+ end
+ for i = m1+1:m
+ p(:,i) = A2\p(:,i-1);
+ end
+ t = (0:tstep:tend) + tstart;
+ x = 0:xstep:l;
+end
--- /dev/null
+clear all
+close all
+
+default_dir = "../Laborka/Zaloha/Matlab/";
+[fname,path] = uigetfile(default_dir, "Select 1 or more experiment files", "Multiselect", "on");
+fname = cellstr(fname);
+conditions.upper_vol = input("Enter the upper volume in mm^3: ")*1e-9; % 872e-9
+conditions.lower_vol = input("Enter the lower volume in mm^3: ")*1e-9; % 600e-9
+conditions.membrane_thickness = input("Enter the membrane thickness in μm: ")*1e-6; % 350e-6
+conditions.membrane_diameter = input("Enter the effective diameter of membrane in mm: ")*1e-3; % 32e-3
+conditions.membrane_mass = input("Enter the effective mass of the membrane in g: ")*1e-3; % 0.4368e-3
+for i = 1:length(fname)
+ simulate_cell([path fname{i}], conditions)
+end
--- /dev/null
+function draw_hist(hist, dims)
+ if length(dims) == 1
+ plot(hist(dims, :))
+ return
+ end
+ points = size(hist)(2);
+ C = [linspace(0, 1, points); ones(1, points); ones(1, points)]';
+ C = hsv2rgb(C);
+ if length(dims) == 2
+ scatter(hist(dims(1), :), hist(dims(2), :), 10, C);
+ elseif length(dims) == 3
+ scatter3(hist(dims(1), :), hist(dims(2), :), hist(dims(3), :), 10, C);
+ else
+ printf("The dimensionality is too high to plot!\n");
+ end
+end
--- /dev/null
+function [D, P] = estimateDP(t, p1, p2, t_start, l, d, V0, V1, V2, T)
+
+ R = 8.314;
+ A = pi * d^2 / 4;
+
+ t = t - t_start;
+ p1 = p1(t>=0);
+ p2 = p2(t>=0);
+ t = t(t>=0);
+
+ M = 20;
+ p1 = filtfilt(ones(M,1)/M,1,p1);
+ p2 = filtfilt(ones(M,1)/M,1,p2);
+ dp1 = p1(2:end) - p1(1:end - 1);
+ dp2 = p2(2:end) - p2(1:end - 1);
+ dt = t(2:end) - t(1:end - 1);
+ J1 = - V1/R/T/A * dp1/dt;
+ J2 = V2/R/T/A * dp2/dt;
+ P1 = J1 * l ./ (p1 - p2);
+ P2 = J2 * l ./ (p1 - p2);
+ P = mean([P1, P2]);
+
+ p_eq = (p1(end) + p2(end)) / 2;
+ n_abs = (max(p1)*V1 + min(p2)*V2 - (V1 + V2)*p_eq)/R/T;
+ Sol = n_abs / (A * l) / p_eq;
+ D = P / Sol;
+end
--- /dev/null
+function [t, x, p] = expl_sol_model(pars)
+% Unpacking x_0
+ tstep = pars(1);
+ xstep = pars(2);
+ tstart = pars(3);
+ tend = pars(4);
+ D = pars(5);
+ P = pars(6);
+ l = pars(7);
+ d = pars(8);
+ V0 = pars(9);
+ V1 = pars(10);
+ V2 = pars(11);
+ pmax = pars(12);
+
+ if length(pars)<16
+ all = false;
+ else
+ all = pars(16);
+ end
+
+ if length(pars)<15
+ T = 298;
+ else
+ T = pars(15);
+ end
+
+ if length(pars)<14
+ pvac = 0;
+ else
+ pvac = pars(14);
+ end
+
+ if length(pars)<13
+ tu = 0.5;
+ else
+ tu = pars(13);
+ end
+
+% Defining helpful constants
+ H = P/D;
+ pmax = pmax * H; # from here on p is used as a concentration
+ R = 8.314;
+ A = pi*d^2/4;
+ K0 = P*R*T*A*tstep/(V1+V0)/2/xstep;
+ K1 = P*R*T*A*tstep/V1/2/xstep;
+ K2 = P*R*T*A*tstep/V2/2/xstep;
+ alpha = D*tstep/xstep^2;
+
+ m1 = tu/tstep+1;
+ m2 = floor((tend - tu)/tstep);
+ n = l/xstep+1;
+ m = m1+m2;
+
+% Initialization of resulting matrix
+ p = zeros(n,m);
+ p(1,1) = pmax;
+
+% Defining matrices of partial differential equations
+ A1 = zeros(n);
+ A1(2:n+1:n*n-2*n-1) = -alpha;
+ A1(n+2:n+1:n*n-n-1) = 1+2*alpha;
+ A1(2*n+2:n+1:n*n-1) = -alpha;
+ A1(1, 1:3) = [1+3*K0 -4*K0 K0];
+ A1(end, end-2:end) = [K2 -4*K2 1+3*K2];
+
+ A2 = zeros(n);
+ A2(2:n+1:n*n-2*n-1) = -alpha;
+ A2(n+2:n+1:n*n-n-1) = 1+2*alpha;
+ A2(2*n+2:n+1:n*n-1) = -alpha;
+ A2(1, 1:3) = [1+3*K1 -4*K1 K1];
+ A2(end, end-2:end) = [K2 -4*K2 1+3*K2];
+
+% Numeric solution of Fick's law
+ for i = 2:m1
+ p(:,i) = A1\p(:,i-1);
+ end
+ for i = m1+1:m
+ p(:,i) = A2\p(:,i-1);
+ end
+ t = (0:tstep:tend) + tstart;
+ x = 0:xstep:l;
+
+p = p/H; # This converts concentration p to equilibrium pressure p
+
+end
--- /dev/null
+function g = grad(fun, args, step = 1e-6)
+ n = size(args);
+ g = zeros(n);
+
+ argdif = zeros(n);
+ for i = 1:max(n)
+ argdif(i) = args(i) * step;
+ f_plus = fun(args + argdif);
+ f_minus = fun(args - argdif);
+ g(i) = (f_plus - f_minus)/(2 * argdif(i));
+ argdif(i) = 0;
+ end
+end
+
--- /dev/null
+function x = halving_int(fun, a, b, prec = 1e-10)
+ fa = fun(a);
+ fb = fun(b);
+ x = (a + b) / 2;
+ fx = fun(x);
+ while (abs(fx) > prec)
+ if (fx * fa < 0)
+ b = x;
+ fb = fx;
+ else
+ a = x;
+ fa = fx;
+ end
+ x = (a + b) / 2
+ fx = fun(x)
+ end
+end
--- /dev/null
+function [Hess, grad, F] = hessgrad(fun, args, step = 0.00001)
+ n = length(args);
+ grad = zeros(n, 1);
+ Hess = zeros(n, n);
+
+ F = fun(args);
+
+ argdif = zeros(n, 1);
+ signledif = zeros(n, 2);
+ for i = 1:n
+ argdif(i) = args(i) * step;
+ singledif(i, 1) = fun(args + argdif);
+ singledif(i, 2) = fun(args - argdif);
+ grad(i) = (singledif(i, 1) - singledif(i, 2))/(2 * args(i) * step);
+ argdif(i) = 0;
+ endfor
+
+ doubledif = zeros(n, n);
+ for i=1:n
+ argdif(i) = args(i) * step;
+ for j=i+1:n
+ argdif(j) = args(j) * step;
+ doubledif(i, j) = fun(args + argdif);
+ doubledif(j, i) = fun(args - argdif);
+ argdif(j) = 0;
+ endfor
+ argdif(i) = 0;
+ endfor
+
+ for i=1:n
+ for j=i:n
+ if j==i
+ Hess(i, i) = (singledif(i, 1) - 2*F + singledif(i, 2))/(args(i)*step)^2;
+ else
+ Hess(i, j) = (doubledif(i, j) + doubledif(j, i) - singledif(i, 1) - singledif(i, 2) - singledif(j, 1) - singledif(j, 2) + 2*F)/(2*args(i)*args(j)*step^2);
+ Hess(j, i) = Hess(i, j);
+ endif
+ endfor
+ endfor
+end
--- /dev/null
+clear all
+close all
+
+pkg load optim
+%{
+classdef Membrane
+ properties
+ D double {mustBePositive}
+ P double {mustBePositive}
+ S double {mustBePositive}
+
+ r double {mustBePositive}
+ d double {mustBePositive}
+ A double {mustBePositive}
+ l double {mustBePositive}
+
+ material char
+ endproperties
+endclassdef
+
+classdef MembraneSimulation
+ properties
+ membrane Membrane
+
+ tstep double {mustBePositive}
+ xstep double {mustBePositive}
+
+ tend double {mustBePositive}
+ tc double {mustBeNonnegative} = 0.5
+
+ pmax double {mustBePositive}
+ pvac double {mustBeNonnegative} = 0
+
+ T double {mustBePositive} = 0
+
+ endproperties
+ methods
+ function [t, p] = simulate();
+ D = membrane.D
+ P = membrane.P
+ l = membrane.l
+ A = membrane.A
+ endfunction
+ endmethods
+endclassdef
+%}
+
+
+[fname,path]=uigetfile;
+data=importdata([path fname], '\t', 2);
+
+datasize=size(data.data);
+t=data.data(:,1)';
+p1=data.data(:,datasize(2))';
+p2=data.data(:,datasize(2)-1)';
+
+pvac = min([p1 p2]);
+pmax = max(p1);
+
+xm=find(p1>0.5*max(p1),1,'first');
+tstart = t(xm);
+
+l = 350e-6;
+d = 32e-3;
+A = d^2/4*pi;
+V0 = 616e-9;
+V1 = 872e-9;
+V2 = 684e-9;
+
+
+
+pars_cur = [1e-3, 3.5e-5, tstart, 100, 3e-10, 3e-13, 3.5e-4, 32e-3, V0, V1, V2, pmax, 0.3];
+
+function [tmod, p] = cellmod(pars)
+ % cellmod(tstep, xstep, tstart, tend, D, P, l, d, V0, V1, V2, pmax, tc=0.5, pvac=0, T = 298, all)
+
+ tstep = pars(1);
+ xstep = pars(2);
+ tstart = pars(3);
+ tend = pars(4);
+ D = pars(5);
+ P = pars(6);
+ l = pars(7);
+ d = pars(8);
+ V0 = pars(9);
+ V1 = pars(10);
+ V2 = pars(11);
+ pmax = pars(12);
+
+
+ if length(pars)<16
+ all = false;
+ else
+ all = pars(16);
+ end
+
+ if length(pars)<15
+ T = 298;
+ else
+ T = pars(15);
+ end
+
+ if length(pars)<14
+ pvac = 0;
+ else
+ pvac = pars(14);
+ end
+
+ if length(pars)<13
+ tc = 0.5;
+ else
+ tc = pars(13);
+ end
+
+ A = d^2/4*pi;
+
+ R = 8.314;
+
+ m = floor(l/xstep)+1;
+ xstep = l/m;
+
+ n1 = floor(tc/tstep)+1;
+ tstep1 = tc/n1;
+
+ if tend>tc
+ n2 = floor( (tend-tc)/tstep);
+ tstep2 = (tend-tc)/n2;
+ p = zeros(n1+n2-1, m);
+ else
+ tstep2 = tstep;
+ p = zeros(n1, m);
+ end
+
+ % criterium for numerical stability
+ crit = D*max([tstep1, tstep2])/xstep^2;
+ while crit>1 % Meaning error grows
+ if tstep1 > tstep2
+ tstep1 = tstep1/2;
+ n1 = n1*2;
+ else
+ tstep2 = tstep2/2;
+ if tend > tc
+ n2 = n2*2;
+ end
+ endif
+ crit = D*max([tstep1, tstep2])/xstep^2;
+ end
+
+ tstep1;
+ tstep2;
+
+ tmod(1) = 0;
+
+ p(1, :) = pvac;
+ p(1,1) = pmax;
+
+ for j=1:n1
+ p(j+1, 2:end-1) = p(j, 2:end-1) + ...
+ D*tstep1/xstep^2*( p(j, 1:end-2) - 2*p(j, 2:end-1) + p(j, 3:end));
+ p(j+1, 1) = p(j, 1) + tstep1/2/xstep*P*R*T*A/(V0+V1)*...
+ ( -3*p(j, 1) + 4*p(j, 2) - p(j, 3) );
+ p(j+1, m) = p(j, m) - tstep1/2/xstep*P*R*T*A/V2*...
+ ( p(j, m-2) - 4*p(j, m-1) + 3*p(j, m) );
+
+ tmod(j+1) = tmod(j) + tstep1;
+ end
+ if tend > tc
+ for j=n1+1:n1+n2-1
+ p(j+1, 2:end-1) = p(j, 2:end-1) + ...
+ D*tstep1/xstep^2*( p(j, 1:end-2) - 2*p(j, 2:end-1) + p(j, 3:end));
+ p(j+1, 1) = p(j, 1) + tstep2/2/xstep*P*R*T*A/V1*...
+ ( -3*p(j, 1) + 4*p(j, 2) - p(j, 3) );
+ p(j+1, m) = p(j, m) - tstep2/2/xstep*P*R*T*A/V2*...
+ ( p(j, m-2) - 4*p(j, m-1) + 3*p(j, m) );
+
+ tmod(j+1) = tmod(j) + tstep2;
+ end
+ end
+
+ tmod = [0 tstart tmod+tstart];
+ p = [pvac*ones(2,size(p,2)); p];
+
+ if all<0.5
+ p = [p(:,1) p(:,end)]; %Erasing all but the boundary (cell) pressures
+ endif
+end
+
+function Obj = objfun(model, pars, exp)
+ xexp = exp(:,1);
+ yexp = exp(:,2:end);
+
+ [xmod, ymod] = model(pars);
+
+ ymod_interpol = interp1(xmod, ymod, xexp);
+
+ valid = ~isnan(ymod_interpol) & (yexp>0);
+
+ Obj = 1/sum(sum(valid)) * sum(sum( (yexp(valid) - ymod_interpol(valid)).^2 ./ yexp(valid) ));
+end
+
+function grad = grad(field, pars, ind, step = 0.001)
+ grad = zeros(size(pars));
+ temp_pars = pars;
+ for i = ind
+ j = j+1;
+
+ temp_pars_pos = pars;
+ temp_pars_pos(i) = temp_pars_pos(i)*(1+step);
+
+ val_pos = field(temp_pars_pos);
+
+ temp_pars_neg = pars;
+ temp_pars_neg(i) = temp_pars_neg(i)*(1-step);
+
+ val_neg = field(temp_pars_neg);
+
+ grad(i) = (val_pos - val_neg)/(temp_pars_pos(i) - temp_pars_neg(i));
+ endfor
+end
+
+tolobj = 1e-4;
+
+curobjfun = @(pars) objfun(@cellmod, pars, [t; p1; p2]');
+
+
+pars_prev = pars_cur;
+pars_prev(5:6) = pars_prev(5:6)*1.01;
+gradF_prev = grad(curobjfun, pars_prev, [5 6]);
+F_prev = curobjfun(pars_prev);
+
+krit = tolobj+1;
+
+optstep = 1;
+
+while F_prev>tolobj
+
+
+ gradF_cur = grad(curobjfun, pars_cur, [5 6])
+
+ pars_cur = pars_cur - optstep * numhessian(curobjfun, pars_cur, [5 6]) \ gradF_cur
+ F_cur = curobjfun(pars_cur)
+%{
+ while F_cur > F_prev
+ optstep = optstep/2;
+ pars_cur = pars_prev - optstep * gradF;
+ F_cur = curobjfun(pars_cur);
+ endwhile
+
+ krit = abs(F_prev - F_cur)
+ F_cur
+ optstep
+ pars_cur
+ gradF_cur
+%}
+ F_prev = F_cur;
+ gradF_prev = gradF_cur;
+ pars_prev = pars_cur;
+
+end
+
+%{
+p1mod = p(:,1);
+p2mod = p(:,end);
+
+t = (0:0.01:tend);
+x = (0: 10e-6 : 350e-6);
+
+[X,T] = meshgrid(x,t);
+
+surf(X,T,p)
+%}
--- /dev/null
+classdef membrane
+ properties
+ name
+ radius
+
--- /dev/null
+function [E, t_calc, x, p] = objfun(model, pars, experimental)
+% experimental = [t_exp, p1_exp, p2_exp]
+ t_exp = experimental(1, :);
+ p1_exp = experimental(2, :);
+ p2_exp = experimental(3, :);
+
+ [t_calc, x, p] = model(pars);
+
+ first_ind = find(t_exp >= min(t_calc), 1, 'first');
+ last_ind = find (t_exp <= max(t_calc), 1, 'last');
+
+ t_exp = t_exp(first_ind:last_ind);
+ p1_exp = p1_exp(first_ind:last_ind);
+ p2_exp = p2_exp(first_ind:last_ind);
+
+ p1_calc = interp1(t_calc, p(1, :), t_exp);
+ p2_calc = interp1(t_calc, p(end, :), t_exp);
+
+ E = (sumsq((p1_calc - p1_exp)) + sumsq((p2_calc - p2_exp)))/(length(p1_calc) + length(p2_calc));
+end
--- /dev/null
+function simulate_cell(full_path, conditions)
+ %{
+ everything is in SI base units
+ conditions
+ mandatory:
+ upper_vol % V1
+ lower_vol % V2
+ membrane_mass % m
+ membrane_thickness % l
+ membrane_diameter % diameter of the part of membrane
+ % that is exposed inside the aparatus
+ voluntary:
+ temperature % temp
+ outer_vol % volume from which the gas expands inside
+ %}
+
+ %{
+ figure
+ surf(t, x, p);
+ shading flat
+ %}
+
+ pkg load signal
+
+ data=importdata(full_path, '\t', 2);
+ datasize=size(data.data);
+ t_exp=data.data(:,1)';
+ p1=data.data(:,datasize(2))'*1e5;
+ p2=data.data(:,datasize(2)-1)'*1e5;
+
+ a=0;
+ b=0;
+
+ pvac = min([p1 p2]);
+ pmax = max(p1);
+
+ xm=find(p1>0.5*max(p1),1,'first');
+ t_start = t_exp(xm);
+
+ R = 8.314;
+ % T = 295;
+ if (isfield(conditions, "temperature"))
+ T = conditions.temperature;
+ else
+ T = 295;
+ endif
+ if (isfield(conditions, "outer_vol"))
+ V0 = conditions.outer_vol;
+ else
+ V0 = inf;
+ endif
+ % V0 = 616e-6;
+ % V1 = 872e-9;
+ % V2 = 600e-9;
+ % l = 350e-6;
+ % d = 32e-3;
+ % m = 0.4368e-3;
+ V1 = conditions.upper_vol;
+ V2 = conditions.lower_vol;
+ l = conditions.membrane_thickness;
+ d = conditions.membrane_diameter;
+ m = conditions.membrane_mass;
+ tu = 0.5;
+
+ A = pi * d^2 / 4;
+
+ tc = 100;
+
+ %[D, P] = estimateDP(t_exp, p1, p2, t_start, l, d, V0, V1, V2, T);
+ D = 1.97093e-9;
+ P = 1.04363e-12;
+ H = P/D;
+
+ %printf("First estimate of membrane characteristics:\tD = %g\tP = %g\n\n", D, P);
+
+ tstep = tu/100;
+ xstep = l/50;
+
+ curobjfun = @(pars) objfun(@basic_model, [tstep, xstep, t_start, 300, pars(1), pars(2), l, d, V0, V1, V2, pmax], [t_exp; p1; p2]);
+
+ x = [D, P];
+ %x = steepest_grad(curobjfun, x, 0.01, 20);
+ x = alternating_min(curobjfun, x, 0.01, 10);
+ D = x(1);
+ P = x(2);
+ printf("Final estimate of membrane characteristics:\tD = %g\tP = %g\n\n", D, P);
+
+ [E, t, l, p] = curobjfun(x);
+
+ printf("The error of the final estimate is %g Pa^2.\n\n", E);
+
+ %mkdir(['Obrázky/PNG/NumFick/' Membrana '_' plyn '_p0=' tlak '_' dnes])
+
+ figure()
+ plot(t, p(1,:), 'r', t, p(end,:), 'r', t_exp, p1, 'b', t_exp, p2, 'b');
+ legend('p\_calc', '', 'p\_exp', '')
+ ylabel('p', 'Rotation', 0)
+ xlabel('t')
+
+ %print(['Obrázky/PNG/NumFick/' Membrana '_' plyn '_p0=' tlak '_' dnes '/tlak.png'],'-dpng','-r360')%,'-painters')
+
+ J1_mem = - P * (-3*p(1,:) + 4*p(2,:) - p(3,:)) / (2*xstep);
+ J2_mem = - P * (3*p(end,:) - 4*p(end-1,:) + p(end-2,:)) / (2*xstep);
+ figure()
+ plot(t, J1_mem, t, J2_mem)
+ hold on
+
+ before_closed = (t - t_start) < tu;
+ V_up = (V0 + V1) * before_closed + V1 * (~before_closed);
+
+ J1_der = -V_up(2:end-1).*(1/A/R/T * (p(1, 3:end) - p(1, 1:end-2)) / (2 * tstep));
+ J1_der(tu/tstep) = (J1_der(tu/tstep - 1) + J1_der(tu/tstep + 1)) / 2;
+ J1_der = [ (-V_up(1) * 1/A/R/T * (-3*p(1, 1) + 4*p(1, 2) - p(1, 3)) / (2 * tstep)) J1_der];
+ J1_der = [ J1_der (-V_up(end) * 1/A/R/T * (3*p(1, end) - 4*p(1, end - 1) + p(1, end - 2)) / (2 * tstep))];
+ J2_der = V2/A/R/T * (p(end, 3:end) - p(end, 1:end-2)) / (2 * tstep);
+ J2_der = [ (V2 * 1/A/R/T * (-3*p(end, 1) + 4*p(end, 2) - p(end, 3)) / (2 * tstep)) J2_der];
+ J2_der = [ J2_der (V2 * 1/A/R/T * (3*p(end, end) - 4*p(end, end - 1) + p(end, end - 2)) / (2 * tstep))];
+ plot(t, J1_der, t, J2_der)
+ legend('J1_x', 'J2_x', 'J1_t', 'J2_t')
+ ylabel('J', 'Rotation', 0)
+ xlabel('t')
+ hold off
+
+ %print(['Obrázky/PNG/NumFick/' Membrana '_' plyn '_p0=' tlak '_' dnes '/tok.png'],'-dpng','-r360')%,'-painters')
+
+ n = zeros(size(J1_mem));
+ n(1) = tstep * (J1_mem(1) - J2_mem(1));
+ for i=2:length(n)
+ n(i) = n(i-1) + tstep * (J1_mem(i) - J2_mem(i));
+ end
+ n = n * A / m;
+ n_der = zeros(size(J1_der));
+ n_der(1) = tstep * (J1_der(1) - J2_der(1));
+ for i=2:length(n_der)
+ n_der(i) = n_der(i-1) + tstep * (J1_der(i) - J2_der(i));
+ end
+ n_der = n_der * A / m;
+ n_bil = 1 / R / T * (max(p1) * V1 - p1 * V1 - p2 * V2) / m;
+ figure()
+ plot(t, n, t, n_der, t_exp, n_bil)
+ legend('n_x', 'n_t', 'n_{bil}')
+ ylabel('n/m', 'Rotation', 0)
+ xlabel('t')
+
+ %print(['Obrázky/PNG/NumFick/' Membrana '_' plyn '_p0=' tlak '_' dnes '/sorpce.png'],'-dpng','-r360')%,'-painters')
+
+ p_end = (p(1,end)*V1 + p(end,end)*V2)/(V1 + V2)/100000;
+ printf("Solubility of the gas in mol/kg/bar by method is:\n\tSpatial\t\tTemporal\tBilance\n\t%d\t%d\t%d\n", n(end)/p_end, n_der(end)/p_end, n_bil(end)/p_end)
+end
--- /dev/null
+function [x, hist] = steepest_grad(objfun, x_0, tol = 0.01, max_iter = 5)
+ s = size(x_0);
+ if s(1) ~= 1 && s(2) ~= 1
+ printf("x_0 must be given as an array of parameters!\n");
+ return ;
+ elseif (s(2) ~= 1)
+ x = x_0';
+ else
+ x = x_0;
+ end
+
+ F = objfun(x);
+ hist = [x; F];
+ reldel = tol + 1;
+ iter = 0;
+ g = - grad(objfun, x);
+ alpha = 0.9 * min(abs(x./g));
+
+ while (abs(reldel) > tol && iter < max_iter)
+ iter = iter + 1;
+ %printf("Iteration number:\t%i\n", iter)
+ % Armijo-Goldstein condition
+ tau = 0.5;
+ c = 0.1;
+ alpha = 0.9 * min(abs(x./g));
+ m = g*g';
+ t = - c * m;
+ F_new = objfun(x + alpha * g);
+ while (F - F_new < alpha * t)
+ %printf("Decreasing alpha.\n")
+ alpha = alpha * tau;
+ if (min(abs((alpha * g) ./ x)) < sqrt(eps))
+ printf("Alpha too small, reached a very sharp \"valey\".\n")
+ return ;
+ end
+ F_new = objfun(x + alpha * g);
+ end
+ reldel = 1 - F_new / F;
+ F = F_new;
+ hist = [hist [(x + alpha * g); F_new]];
+ if (any(isnan(x + alpha * g)))
+ printf("NaN encountered, ending minimization!\n");
+ break
+ end
+ x = x + alpha * g;
+ %{
+ printf("\nImprovement:\t\t%f %%\n", reldel * 100);
+ printf("Current parameters:");
+ for i = 1:length(x)-1
+ printf("\t%g, ", x(i));
+ end
+ printf("\t%g\n\n", x(end));
+ %}
+ g = - grad(objfun, x);
+ end
+end
projects / Uchytilap.git / commitdiff