clc
clear
D = 1e-3;
cdis = 2.0;
I_factor=1.25;
ops = optimset ('TolX',1e-9, 'TolFun', 1e-9, 'MaxFunEvals', 1e6, 'MaxIter',1e6);
ops2 = optimset ('TolX',1e-9, 'TolFun', 1e-9, 'MaxFunEvals', 1e6, 'MaxIter',1e6);
diamN = 0.2;
diamD = 5.0;
diamH = 1;
LH = 1.0;
l = 0.5;
LN = 2.0;
Lcrit_lambda = [];
Lcrit_lambda_simp=[];
Lbase = [0.05:0.01:0.5];
lambda_arr=[20:20:720];
I0_arr=[];
I0_step_arr=[];
hill_slope=300;
for lambda = lambda_arr
s = lambda/2/D;
a = coth(LN/lambda) + (diamN/diamH)^2 * coth(LH/lambda);
B = cosh(LN/lambda) * a - 1/sinh(LN/lambda);
Q = (D/lambda)*sinh(LN/lambda) * a;
Q = Q/B;
P = cosh(l/lambda)/sinh(LH/lambda);
P = P/B;
brk = coth(LH/lambda)+(diamH/diamN)^2*tanh(LN/lambda);
alf = cosh(l/lambda)/(sinh(LH/lambda)*cosh(LN/lambda)*brk);
bet_new = -(lambda/D)*(cosh(l/lambda)/(sinh(LH/lambda)^2 *brk));
bet_new = bet_new*(cosh(l/lambda)-sinh(LH/lambda)*cosh((LH-l)/lambda)*brk);
FA=(((alf * ((diamN/diamD)^2)) *P)/((1+((lambda/(2*D))*((diamN/diamD)^2))*Q)))+(2*D/lambda)*bet_new;
FB=(((2*alf * ((diamN/diamD)^2)) *P)/((1+(lambda/(2*D))*((diamN/diamD)^2)*Q)));
I0_crit_step = (2*D/lambda)*cdis/FA;
I0_crit_range=I0_crit_step*[0.9:0.2/500:1.1];
ln_range=length(I0_crit_range);
ch_2=zeros(ln_range,2);
for ii=1:ln_range
I0=I0_crit_range(ii);
iso_func = @(x) (-x + (lambda/(2*D))*I0*FA*hill(x,cdis,hill_slope));
ch_2(ii,1)=fsolve(iso_func,cdis/5,optimoptions('fsolve','Display','off'));
UP_step=(lambda/(2*D))*(FA*I0*1.01);
ch_2(ii,2)=fsolve(iso_func,UP_step,optimoptions('fsolve','Display','off'));
end
i_crit_ind=find(abs(ch_2(:,2)-ch_2(:,1))>max(ch_2(:,2))/1000,1);
I0_crit=I0_crit_range(i_crit_ind);
step_error=100*(I0_crit-I0_crit_step)/I0_crit;
I0 = I0_crit*I_factor;
I0_step=I0_crit_step*I_factor;
I0_arr=[I0_arr,I0];
I0_step_arr=[I0_step_arr,I0_step];
L_ax=[lambda/100:lambda/100:lambda];
peak = @(x) (-x + (lambda/(2*D))*I0*FA*hill(x,cdis,hill_slope));
cmin = fminbnd(peak,0,cdis*1.1, ops);
rhs=(2*D/lambda)*cmin/(FB*I0) - hill(cmin,cdis,hill_slope)*FA/FB;
fun = @(L) rhs - geoser(L,lambda);
MMM=lambda*log(1+I_factor*FB/FA);
Lcrit_lambda_simp=[Lcrit_lambda_simp,MMM];
LLL=fsolve(fun,MMM, optimoptions('fsolve','Display','off'));
Lcrit_lambda = [Lcrit_lambda,LLL];
end
figure(10)
F10=loglog(lambda_arr,Lcrit_lambda,'LineWidth',3);
hold on
F10a=loglog(lambda_arr,Lcrit_lambda_simp,'-.','LineWidth',2);
title("Phase Diagram")
xlabel("\lambda (\mu m)")
ylabel("Distance between spines (\mu m)")
lam_ind=6;
lambda=lambda_arr(lam_ind)
I0=I0_arr(lam_ind);
Lcrit=Lcrit_lambda(lam_ind);
I0_step=I0_step_arr(lam_ind);
Lcrit_step=Lcrit_lambda_simp(lam_ind);
s = lambda/2/D;
a = coth(LN/lambda) + (diamN/diamH)^2 * coth(LH/lambda);
B = cosh(LN/lambda) * a - 1/sinh(LN/lambda);
Q = (D/lambda)*sinh(LN/lambda) * a;
Q = Q/B;
P = cosh(l/lambda)/sinh(LH/lambda);
P = P/B;
brk = coth(LH/lambda)+(diamH/diamN)^2*tanh(LN/lambda);
alf = cosh(l/lambda)/(sinh(LH/lambda)*cosh(LN/lambda)*brk);
bet_new = -(lambda/D)*(cosh(l/lambda)/(sinh(LH/lambda)^2 *brk));
bet_new = bet_new*(cosh(l/lambda)-sinh(LH/lambda)*cosh((LH-l)/lambda)*brk);
FA=(((alf * ((diamN/diamD)^2)) *P)/((1+((lambda/(2*D))*((diamN/diamD)^2))*Q)))+(2*D/lambda)*bet_new;
FB=(((2*alf * ((diamN/diamD)^2)) *P)/((1+(lambda/(2*D))*((diamN/diamD)^2)*Q)));
L_arr=[Lcrit/2:Lcrit/400:Lcrit*2];
Ch_up=[];
Ch_down=[];
Ch_up_mod=[];
Ch_down_mod=[];
Ch_int_mod=[];
Lcrit_mod = 0;
for L=L_arr
Up=(lambda/(2*D))*I0*(FA+FB*geoser(L,lambda));
Down=(lambda/(2*D))*I0*(FB*geoser(L,lambda));
fun = @(x) (-x + (lambda/(2*D))*I0*(FA*hill(x,cdis,hill_slope) + FB*geoser(L,lambda)));
Ch_up_mod = [Ch_up_mod, fsolve(fun,Up,optimoptions('fsolve','Display','off'))];
dn_mod = fsolve(fun, Down,optimoptions('fsolve','Display','off'));
if (dn_mod < cdis) && (sign(fun(dn_mod+1e-6)) * sign(fun(dn_mod-1e-6)) < 0)
Ch_down_mod = [Ch_down_mod, dn_mod];
Ch_int_mod = [Ch_int_mod, fsolve(fun,cdis,optimoptions('fsolve','Display','off'))];
if Lcrit_mod == 0
Lcrit_mod = L - Lcrit/800;
end
else
Ch_down_mod = [Ch_down_mod, NaN];
Ch_int_mod = [Ch_int_mod, NaN];
end
Ch_up=[Ch_up,Up];
if L < Lcrit
Ch_down=[Ch_down,Up];
else
Ch_down=[Ch_down,Down];
end
end
figure(110)
plot(L_arr,Ch_up,'-b','LineWidth',3)
hold on
plot(L_arr,Ch_up_mod,':b','LineWidth',3)
plot(L_arr,Ch_down,'-r','LineWidth',3)
plot(L_arr,Ch_down_mod,':r','LineWidth',2)
plot(L_arr,Ch_int_mod,'--k','LineWidth',2)
title("\lambda=",lambda)
xlabel("L (\mu m)")
ylabel("C_h")
Ch_up_step=[];
Ch_down_step=[];
for L=L_arr
Up_step=(lambda/(2*D))*I0_step*(FA+FB*geoser(L,lambda));
if L>Lcrit_step
Down_step=(lambda/(2*D))*I0_step*(FB*geoser(L,lambda));
else
Down_step=Up_step;
end
Ch_up_step=[Ch_up_step,Up_step];
Ch_down_step=[Ch_down_step,Down_step];
end
figure(110)
plot(L_arr(1:20:length(L_arr)),Ch_up_step(1:20:length(L_arr)),'b+','LineWidth',3)
hold on
plot(L_arr(1:20:length(L_arr)),Ch_down_step(1:20:length(L_arr)),'rx','LineWidth',2)
Lstep_err=100*(Lcrit-Lcrit_step)/Lcrit;
disp(['Lcrit step % error= ',num2str(Lstep_err)])
function X = myFzero(fun, x0, tol)
X = []
for x = x0
guess = x
val = 2*tol
while abs(val) > tol
[root, val, info, out] = fsolve(fun, guess);
guess = guess + 1;
end
X = [X, root];
end
end
function s = hill(x,cdis,hs)
s = x.^hs./(x.^hs + cdis^hs);
end
function S= geoser(L, lam) % Infinite case start from 1
S=exp( -L /lam) ./(1 -exp(-L/lam));
end