%
% << hopf >> 
%
function hopf
clear all
close all
global mu alpha omega seuil
%
disp('Bifurcation hopf. Patientez un peu SVP !')
% parametres 
tmax=1.5; nt=20;  nx0=11; nmu=11;
tmax=1.5; nt=20;  nx0=11; nmu=6;
%tmax=1; nt=10;  nx0=4; nmu=4;
%
alpha=1; omega=2;
tmin=0;
dt=(tmax-tmin)/(nt-1);
tspan=[tmin:dt:tmax];
%%%
% initialize x
x0max=2; x0min=-x0max;
x0maxe=x0max*1.05; x0mine=x0min*1.05;
dx0=(x0max-x0min)/(nx0-1);
x0v=x0min:dx0:x0max;
% seuil abs(x) < seuil pour calcul f(x)
seuil=x0max*1.5;
% initialize mu
mumax=2; mumin=-.5;
mumaxe=mumax*1.05; mumine=mumin*1.05;
dmu=(mumax-mumin)/(nmu-1);
muv=mumin:dmu:mumax;
%
% initialize mun
mun=mumine:.01:0;
% initialize mup
dmup=mumaxe/100; mup=0:dmup:mumaxe;
sqrtmup=sqrt(mup/alpha);
% initialize MUP for cylinder 
MUP=0:.1:1; scale=.4*mumaxe;
[X,Y,Z]=cylinder(MUP,20);
XX=X*sqrt(alpha*scale); YY=Y*sqrt(alpha*scale); ZZ=Z.^2*scale;
%
[XC,YC,ZC]=cylinder(1,20);
sc=sqrt(mumax/alpha);
ZZC=mumax+(mumaxe-mumax)*ZC; XXC=sc*XC; YYC=sc*YC;
% Figure d'attente
hold off
cbrouges(mumine,mumaxe,x0mine,x0maxe,mup,sqrtmup,mun,ZZ,XX,YY,ZZC,XXC,YYC)
for imu=1:nmu
	mu=muv(imu);
	for ix0=1:nx0
		plot3(muv(imu),x0v(ix0),x0v(ix0),'o','Color','k','LineWidth',2)
	end
end
disp('Changement de taille de l''image pendant l''attente !')
drawnow
%pause
% Resolution de l'ode
nla=0; nlamax=nt-1;
y0=x0max;
for imu=1:nmu
	for ix0=1:nx0
		mu=muv(imu);
		z0=[x0v(ix0) x0v(ix0)];
		[t,z]=ode23('ho',tspan,z0);
		for it=1:nt
			xa(imu,ix0,it)=z(it,1);
			ya(imu,ix0,it)=z(it,2);
		end
	end
end
% Boucle temporelle
hold off
for it=1:nt
	%% courbes rouges
	cbrouges(mumine,mumaxe,x0mine,x0maxe,mup,sqrtmup,mun,ZZ,XX,YY,ZZC,XXC,YYC)
	%
	if it==1
		for imu=1:nmu
			for ix0=1:nx0
				plot3(muv(imu), xa(imu,ix0,it),ya(imu,ix0,it),'o','Color','k','LineWidth',2)
			end
		end
	elseif  it==2
		for imu=1:nmu
			for ix0=1:nx0
				muvimu=muv(imu); mudeux=[muvimu muvimu];
				xadeux=[xa(imu,ix0,it) xa(imu,ix0,it-1)];
				yadeux=[ya(imu,ix0,it) ya(imu,ix0,it-1)];
				plot3(mudeux,xadeux,yadeux,'-b')
				plot3(muv(imu), xa(imu,ix0,it),ya(imu,ix0,it),'o','Color','k','LineWidth',2)
				plot3(muv(imu), xa(imu,ix0,it-1),ya(imu,ix0,it-1),'o','Color','b','LineWidth',1)
			end
		end
	else
		for imu=1:nmu
			for ix0=1:nx0
				muvimu=muv(imu); mutrois=[muvimu muvimu muvimu];
				xatrois=[xa(imu,ix0,it) xa(imu,ix0,it-1) xa(imu,ix0,it-2)];
				yatrois=[ya(imu,ix0,it) ya(imu,ix0,it-1) ya(imu,ix0,it-2)];
				plot3(mutrois,xatrois,yatrois,'-b')
				plot3(muv(imu), xa(imu,ix0,it),ya(imu,ix0,it),'o','Color','k','LineWidth',2)
				plot3(muv(imu), xa(imu,ix0,it-1),ya(imu,ix0,it-1),'o','Color','b','LineWidth',1)
				plot3(muv(imu), xa(imu,ix0,it-2),ya(imu,ix0,it-2),'o','Color','b','LineWidth',1)
			end
		end
	end
	pause(0)
	hold off
	if nla==0 
		%disp(' You can change the window size now !')
		%pause
		Mla=moviein(nlamax);
	end
	nla=nla+1;
	Mla(:,nla)=getframe; 
	%filew(gcf,it);
end
%%% N times, FPS frame per second
%pause
movie(Mla,5,5)
return
%
% << cbrouges >>
%
function cbrouges(mumine,mumaxe,x0mine,x0maxe,mup,sqrtmup,mun,ZZ,XX,YY,ZZC,XXC,YYC)
%% courbes rouges
%
plot(mup,sqrtmup,'r','LineWidth',2)
axis([mumine mumaxe x0mine x0maxe x0mine x0maxe])
grid
hold on
plot(mup,-sqrtmup,'r','LineWidth',2)
plot(mup,0*mup,'--r','LineWidth',2)
plot(mun,0*mun,'r','LineWidth',2)
surf(ZZ,XX,YY)
hold on
surf(ZZC,XXC,YYC)
xlabel('\mu','FontSize',20)
ylabel('x','FontSize',20)
zlabel('y','FontSize',20)
return
