%% Secondary buckling of pinned-pinned
%% Free standing fold
%------------------------------------------
clear all
close all
tot=100 % number of points along the branch
%% Initialize variables::
%----------------------------------------------------------------------
c_l=0.05; % constant
c_r=linspace(0.0001,2*c_l,tot);            % assign of "clearance" of wall-wall segment
alpha=zeros(tot,1);R_x=zeros(tot,1);       % inclination of force at left segment and the horizontal load
Len_L=zeros(tot,1);Len_R=zeros(tot,1);     % lengths of left and right segments 
Theta_L=zeros(tot,1);Theta_R=zeros(tot,1); % inclination of inflection points of left and right segments
for i=1:tot
lb = [0,0,0,0,0];  % set up lower bounds
rng default
%% Gradual assignment of initial guess for faster solution::
%--------------------------------------------------
if i==1
%% Known initial point from previous contact pattern
%--------------------------------------------------------------
x0 = [0.4,0.000,0.2,157,0.25]';
else
%% set up initial guess the solution of the previous solution::
%---------------------------------------------------------------
x0 = [x(1),x(2),x(3),x(4),x(5)]';
%% x(1): theta_o_l
%% x(2): theta_o_r
%% x(3): alpha_l
%% x(4): force
%% x(5): length_l
end
c=[c_r(i),c_l]; % predefined clearances of left and right segments
options = optimoptions(@lsqnonlin,'Algorithm','trust-region-reflective','FunctionTolerance',1e-18);
%options = optimoptions(@lsqnonlin,'Algorithm','levenberg-marquardt','FunctionTolerance',1e-18);
[x,res] = lsqnonlin(@(x)fbnd_secondary(x,c),x0,lb,[],options);
%% calculation of end-shortening delta::
%--------------------------------------------------------------
len_l=(0.5-x(5))/2; % length of left segment
fun_1=@(s) cos((x(1)-x(3))*cos(sqrt(x(4))*s)+x(3));
x_left(i) = integral(fun_1,0,x(5));
fun_2=@(s) cos(x(2)*cos(sqrt(x(4))*s));
x_right(i) = integral(fun_2,0,len_l);
delta(i)=abs(1-(2*x_left(i)+4*x_right(i)));
%% Horizonta force:
%--------------------------------------------------------------
R_x(i)=abs(x(4));
end
%% Plot of Horizontal force - endshortening::
%---------------------------------------------------
figure 
plot(delta(2:tot),R_x(2:tot)/(pi^2))
title ('Free standing fold')
xlabel('End-shortening, \delta')
ylabel ('Applied force, Q/\pi^{2}')