%% Práctica sobre el MODELO DE LORENZ-HAKEN (2b)
% Programa traducido de Mathematica a Matlab por Fernando Hueso González.
% Óptica cuántica (Seminarios de teoría semiclásica y simulaciones) - 
% 4º de Grado de Física, 31 de mayo de 2011.
global s b r
syms d s b r
format compact

% Solución estacionaria (analítica)
d = 1;
P = sqrt(r-1);
F = sqrt(r-1);
solestac = [d,P,F];

% Umbral de la bifurcación de Hopf
rHB=-s*(3+b+s)/(1+b-s);

%% Ecuaciones y resolución del modelo de Lorenz-Haken

s=3; b=1; r=20;
Tfin=1000; h=0.02; t=0:h:Tfin;
[t,solnum] = rungekut('laser',t,[0.001,0.001,1]);
tmin=900; tmax=1000; nmin=find(t>=tmin,1,'first'); nmax=find(t<=tmax,1,'last');

close all
titul='Hard-mode excitation 1D';
figure('Name',titul),plot(t(nmin:nmax), solnum(nmin:nmax,1).^2);
axis tight, title(titul), xlabel('t'), ylabel('F^2')

titul='Hard-mode excitation 2D';
figure('Name',titul),plot(solnum(nmin:nmax,1), -solnum(nmin:nmax,3));
axis tight, title(titul), xlabel('F'), ylabel('-d')

titul='Hard-mode excitation 3D';
figure('Name',titul),plot3(solnum(nmin:nmax,1), solnum(nmin:nmax,2),-solnum(nmin:nmax,3));
axis tight, title(titul), xlabel('F'), ylabel('P'), zlabel('-d')
%% Transformada de Fourier
tmin=500; tmax=Tfin; nmin=find(t>=tmin,1,'first'); nmax=find(t<=tmax,1,'last');
FT=fft(solnum(nmin:nmax,1).^2); L_FT=length(FT)

close all
titul='Transformada de Fourier';
figure('Name',titul),semilogy(abs(FT(1:1000)).^2);
axis tight, title(titul)
%% Mapas de máximos
% Verificar que has ejecutado previamente la orden solnum.
lista=solnum(nmin:nmax,1).^2; L=length(lista)

m=[];
for k=2:L-1
    if lista(k-1)<lista(k) && lista(k+1)<lista(k)
        parab=polyfit([-1;0;1], lista(k-1:k+1),2);
        x=-parab(2)/2/parab(1);
        y=polyval(parab,x);
        m=[m,y];
    end
end
L_m=length(m)
% implemento la interpolación parabólico en todos los cálculos, al haber
% comprobado que da mejor resultado.
close all
titul='Mapa de máximos';
figure('Name',titul),plot(m(1:end-1),m(2:end),'.')
xlabel('F^2max_i'), ylabel('F^2max_{i+1}'), title(titul);
%% Doblamiento de periodo
close all
s=3; b=1; r=99.0;
Tfin=100; h=0.01; t1=0:h:Tfin;
[t1,solnum1] = rungekut('laser',t1,[1,1,1]);

tmin=80; tmax=Tfin; nmin=find(t1>=tmin,1,'first'); nmax=find(t1<=tmax,1,'last');
titul='Doblamiento de periodo 1D';
figure('Name',titul),plot(t1(nmin:nmax), solnum1(nmin:nmax,1).^2);
axis tight, title(titul), xlabel('t'), ylabel('F^2')

titul='Doblamiento de periodo 2D';
figure('Name',titul),plot(solnum1(nmin:nmax,1), solnum1(nmin:nmax,2));
axis tight, title(titul), xlabel('F'), ylabel('P')
%%
tmin=80; tmax=Tfin;
nmin=find(t1>=tmin,1,'first'); nmax=find(t1<=tmax,1,'last');

lista=solnum1(nmin:nmax,1).^2; L=length(lista);
FT1=fft(solnum1(nmin:nmax,1).^2); L_FT1=length(FT1)

close all
titul='Transformada de Fourier';
figure('Name',titul),semilogy(abs(FT1(1:200)).^2);
axis tight, title(titul)
%Se ve la generación del 2º armónico al haberse duplicado el periodo

m=[];
for k=10:L-1
    if lista(k-1)<lista(k) && lista(k+1)<lista(k)
        parab=polyfit([-1;0;1], lista(k-1:k+1),2);
        x=-parab(2)/2/parab(1);
        y=polyval(parab,x);
        m=[m,y];
    end
end
L_m=length(m)

titul='Mapa de máximos';
figure('Name',titul),plot(m(1:end-1),m(2:end),'.')
title(titul); xlabel('F^2max_i'), ylabel('F^2max_{i+1}')
% Se observan puntos sueltos, lo cual es lógico al no estar en una zona
% caótica o autopulsada, sino que hay una oscilación de amplitud y
% frecuencias bien definidas.
%%
close all