- Author:
- aram148 <a.rampadarath@auckland.ac.nz>
- Date:
- 2021-11-04 16:02:42+13:00
- Desc:
- Updated USMC-Bursztyn model
- Permanent Source URI:
- http://models.cellml.org/workspace/6b0/rawfile/4c1ab73f48d7150c2094cf8017425482a1594ccf/USMC/DM+Mijaailovich_huxley_copy/Zah.m
function [R,T,force,stiffness,em10,n1,P1,Q1,tt]=Zah(trange,N,M,xrange,f)
%[R,T,force,stiffness,em11,em10,em12,n1,P1,Q1,tt,vee]
%N number of time points
%trange specifies the interval of time steps
%T time grid
%X0 space grid
T=linspace(trange(1),trange(2),N);
X0=linspace(xrange(1),xrange(2),M);
%dt=timestep
dt=(trange(2)-trange(1))/(N);
%Initial conditions
M10=1;
M11=0.5;
M12=0.4;
R=[M10;M11;M12];
g=@(r) mis(r,N,M,f);%to pass into fsolve to find equilibria solutions
options = optimoptions('fsolve');
[R,~]=fsolve(g,R,options);
for j=1:N
m10=R(1);
m11=R(2);
m12=R(3);
% t=0+j*dt;
e=dt;
tt(j)=0+j*e;
X=X0;%-lengthvar(t,f);%update X for each L as L updates for each T
% L=lengthvar(t,f);
% L=(-25*sin(2*pi*f*t));
% L=0;
[alpha,V] =lengthvar(tt(j),f);
vee(j)=V;
fp1=0.06;
g1=0.10;
g2=0.0040;
g3=0.01;
%Mean for cdf
p1=m11/m10;
P1(j)=m11/m10;
%p1=0.5;
%p2=M21(i)/M20(i);
%Standard deviation for cdf
q1=sqrt(((m12/m10)-((m11/m10).^(2))));
Q1(j)=sqrt(((m12/m10)-((m11/m10).^(2))));
r0=-p1./q1;
% r0=0;
r1=(1-p1)./q1;
%r1=1;
rinf=(Inf-p1)./q1;
%rinf=Inf;
phi0=cdf('Normal',r0,p1,q1);
% phi0=erf(-p1./q1);
phi1=cdf('Normal',r1,p1,q1);
%phi1=erf((1-p1)./q1);
phinf=cdf('Normal',rinf,p1,q1);
%phinf=erf((Inf-p1)./q1);
I0=-(exp(-((-p1./q1).^(2))/2))/(sqrt(2*pi));
I1=-(exp(-((((1-p1)./q1)).^(2))/2))/(sqrt(2*pi));
I2=-(exp(-(((Inf-p1)./q1)).^(2))/2)/(sqrt(2*pi));
%Functions for the rhs of the first PDE M1_lambda
J0=phi0;
J01=phi1;
J0inf=phinf;
J10=((p1.*phi0)+(q1.*I0));
J11=(p1.*phi1)+(q1.*I1);
J12=(p1.*phinf)+(q1.*I2);
J20=((p1.^(2)).*phi0)+((2*p1.*q1).*I0)+((q1.^(2)).*(phi0+(r0*I0)));
J21=((p1.^(2)).*phi1)+((2*p1.*q1).*I1)+((q1.^(2)).*(phi1+(r1*I1)));
J22=((p1.^(2)))+((2*p1.*q1).*I2)+((q1.^(2)));
J30=(p1.^(3).*phi0)+((3.*p1.^(2).*q1).*I0)+((3.*p1.*q1.^(2)).*((phi0)+(r0*I0)))+((q1.^(3).*(2+(r0.^(2))).*I0));
J31=(p1.^(3).*phi1)+((3.*p1.^(2).*q1).*I1)+((3.*p1.*q1.^(2)).*((phi1)+(r1*I1)))+(q1.^(3).*(2+(r1.^(2))).*I1);
J32=(p1.^(3))+((3.*p1.^(2).*q1).*I2)+((3.*p1.*q1.^(2)));
%Functions defined for the RHS of the second PDE M2_lambda
%Components for the matrix F that will represent each moment,
%M1_lambda and M2_lambda
%
b0=fp1/2;
b1=fp1/3;
b2=fp1/4;
F0=((g2*J0)+((fp1+g1)*(J11-J10))+(g1*(p1-J11))+(g3*(p1-J11-1+J01)))*m10;
F1=((g2*J10)+((fp1+g1)*(J21-J20))+(g1*(p1.^(2)+q1.^(2)-J21))+(g3*((p1.^(2)+q1.^(2))-J21-p1+J11)))*m10;
F2=((g2*J20)+((fp1+g1)*(J31-J30))+(g1*(p1.^(3)+3*p1*q1.^(2)-J31))+(g3*((p1.^(3)+3*p1*q1.^(2))-J31-(p1.^(2)+q1.^(2))+J21)))*m10;
% Putting together the matrix F for the RHS of the system of
% equations
% F=[b0-F0*m10;b1-F1*m10;b2-F2*m10];
F=[b0-F0;b1-F1-V*m10;b2-F2-2*V*m11];
Rnew=R+dt*F;
R=Rnew;
%m10=r(1);
% m=1;
n1=(R(1)/((sqrt(2*pi))*q1)).*exp(-((X-p1).^(2))./(2*(q1.^(2))));
force(j)=trapz(X,X.*(n1));
stiffness(j)=trapz(X,n1);
em11(j)= R(2);
em10(j)=R(1);
em12(j)=R(3);
% % m=ones(1,N/5000);
% t=zeros(1,N/5000);
% % em11=zeros(1,N/5000);
% em10=zeros(1,N/5000);
% % em12=zeros(1,N/5000);
%
% if (mod(j,5000)==0)% && j>=16500 && j<=16700)
%
% em10(m)=R(2);
%
% % n=m+1;
% % % t(m)=tt;
% % em11(1,n)= R(2);
% % em10(n)=R(1);
% % em12(n)=R(3);
% % n=n+1;
% end
% mnew=m+1;
% m=mnew;
% %
% % figure(1)
% % plot(X0,n1)
% % drawnow
% % hold on
% % q1
% %
% % pause
% % % [
% % figure(2)
% % plot(T,R(1))
% %
% % drawnow
% % hold on
% % pause(2)
end
