%%   Date: 01/03/2021
% 
%    Author: Anand Rampadarath
%    Version 1
%    This model represents the Uterine smooth muscle cell as 
%    developed by Maggio et al 2012. It is a variant of the 
%    Hai-Murphy model including spatian dependence of the 
%    cross-bridges, similar to Mijailovich et al. 
%    We us the DM approximation to reduce the computational complexity.
% 
function [X0,T,R,force,stiffness,phosphorylation,m10,m20,nmp,nm]=ZahM(trange,N,M,xrange)
%N number of time points
%trange specifies the interval of time steps
%T time grid
%X0 space grid
%gamma=25 rho=1 s=airway generation lambda=force activation

   T=linspace(trange(1),trange(2),N);
   X0=linspace(xrange(1),xrange(2),M);
   dt=(trange(2)-trange(1))/(N); % dt must be <=0.1
  
   M10=0.309120293063990;
   M11=0.118325676501984;
   M12=0.066855519431672;
   M20=0.134379171029593;
   M21=-0.045132660622671;
   M22=0.101673063419951;
   C=0.6;
    
   R=[M10;M11;M12;M20;M21;M22;C];
  
%     g=@(r) mis1(r);
%     options = optimoptions('fsolve');
%     [R,~,exitflag]=fsolve(g,R,options);

 for i=1:N %loop over number of timesteps to update then run euler
  
        t=trange(1)+i*dt;
         X=X0;
      %V=velvar(t,f); 
%       [~,V]=lengthvar(t,f);
      
      [F1,~,~,~,~]=DM_funcs_1(t,R);
      [F2,~,~,~,~]=DM_funcs_1(t+(dt/2),R+(dt/2)*F1);
      [F3,~,~,~,~]=DM_funcs_1(t+(dt/2),R+(dt/2)*F2);
      [F4,q1,q2,p1,p2]=DM_funcs_1(t+dt,R+(dt*F3));

      Rnew = R + (dt/6)*(F1+(2*F2)+(2*F3)+F4);
      R=Rnew;
      
      % A rebuild of the distribution approximation for 
     n1=(R(1)/((sqrt(2*pi))*q1)).*exp(-((X-p1).^(2))/(2*(q1.^(2))));
     n2=(R(4)/((sqrt(2*pi))*q2)).*exp(-((X-p2).^(2))/(2*(q2.^(2))));
     nmp(i)=1-R(7)-R(1);
     nm(i)=1-R(4)-R(1)-nmp(i);
     m10(i)=R(1);
     m11(i)=R(2);
     m21(i)=R(5);
     m20(i)=R(4);
     force(i)=R(2)+R(5);
     stiffness(i)=R(1)+R(4);
     nMpave(i)=nmp(i);
     nAmpave(i)=R(1);
     phosphorylation(i)=(nMpave(i)+nAmpave(i));

 end
