clear all

global FD
figno=18;

TL=[pi/8,pi/4];
TL=linspace(0.001,pi/2+pi/4,801);
TL=pi/4;

for ia =1:length(TL)

thetaL=TL(ia);

dtheta0=sqrt(2*(1-cos(thetaL)));

xt=[0;dtheta0];
Nt=1001;
tt=linspace(0,4*pi,Nt);
dt=tt(2)-tt(1);

x(1)=xt(1);y(1)=xt(2);
iter=0;
for it =2:Nt

    k1=dt*dpendulum(xt);
    k2=dt*dpendulum(xt+k1/2);
    k3=dt*dpendulum(xt+k2/2);
    k4=dt*dpendulum(xt+k3);

    xt=xt+(1/6)*(k1+2*k2+2*k3+k4);

    x(it)=xt(1);y(it)=xt(2);
    
    if y(it)*y(it-1)<0
        iter=iter+1;
        if iter<4
        xL(ia,iter)=x(it);
        tL(ia,iter)=tt(it);
        THETA(ia,iter)=thetaL;
        end
    end

end
fc(ia,:)=tL(ia,:)/(pi/2)



% figure(67);hold on
% plot(x,y,'r','linewidth',2)
% xlabel('$\theta$','interpreter','latex')
% ylabel('$\dot \theta$','interpreter','latex')
% set(gca,'Fontsize',24)
% box on
% grid on


figure(67);
 plot(fc,THETA,'.k',fc,-THETA,'.k','linewidth',2)
 xlabel('$\sqrt{F_c/D} L/(\pi/2)$','interpreter','latex')
 ylabel('$\theta_L$','interpreter','latex')
 set(gca,'Fontsize',24)
box on
 grid on


 %determination of the shapethe shape
 Ns=401;
for ip=1:3


ss=linspace(0,tL(1,ip),Ns)
ds=ss(2)-ss(1);
FD=1/tL(ia,ip);
xt=[0;sqrt(2*(1-cos(thetaL)));0;0];

 for ib =2:Ns

    k1=ds*feuler_shape(xt);
    k2=ds*feuler_shape(xt+k1/2);
    k3=ds*feuler_shape(xt+k2/2);
    k4=ds*feuler_shape(xt+k3);

    xt=xt+(1/6)*(k1+2*k2+2*k3+k4);

    xx(ib)=xt(3);yy(ib)=xt(4);

 end
ip
figure(figno)
 plot(xx,yy,'r',-xx,yy,'b','linewidth',2)
 xlabel('$x$','interpreter','latex')
 ylabel('$y$','interpreter','latex')
 set(gca,'Fontsize',24)
box on
 grid on
 hold on
 pause



end
end

figure(figno);hold off











