%% Práctica sobre el MODELO DE LORENZ-HAKEN (2)
% 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
% Ecuaciones diferenciales del láser definidas en el archivo laser.m
%F'=s*(P-F)
%P'=-P+d*F
%d'=b*(r-d-F*P)

s=3; b=1; r=20;%variables globales
Tfin=1000;%tiempo hasta el que calculas

h=0.01; t=0:h:Tfin; %Define el paso h que utiliza Runge-Kutta
[t,solnum] = rungekut('laser',t,[0.001,0.001,1]);

tmin=900;%mínimo del eje de abcisas del plot
tmax=1000;%máximo del eje de abcisas del plot

nmin=find(t>=tmin,1,'first');%forma de decirle al Matlab hasta donde plotear
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')
%%
% En esta ocasión, se fija un tmin>0 para saltarnos el transitorio hasta
% que se cae en el atractor caótico. De hecho, es fácil comprobar que si
% representamos de 0 a 100 en la escala de tiempos, se ven oscilaciones que
% van aumentando progresivamente (figura 1D) hasta llegar al régimen puramente caótico.
% Por tanto, el transitorio sería el intervalo entre 0 y 60 aproximadamente
% hasta ver el carácter autopulsado.
% Se observa un carácter similar autopulsado al del plot de Mathematica,
% pero con valores distintos dado que estamos en la solución caótica y hay
% dependencia grande con el método numérico utilizado. De hecho, si se
% cambia el paso h en Runge Kutta, el aspecto varía también
% considerablemente (análogo a la sensibilidad de las condiciones iniciales).
% En resumen, se observa claramente la estabilidad (la solución NO es estacionaria con un valor fijo, pero sí
% es estable u oscilante acotadamente) de la solución en torno al atractor
% caótico (atractor de Lorenz).
% Aparte, cabe apuntar que se ha elegido un paso de Runge-Kutta de 0.01
% para que el cálculo no se haga muy largo (al menos en mi ordenador, que
% es poco potente). Se puede aumentar la precisión bajando h, pero se
% necesitaría más tiempo de cálculo. Por ejemplo, si quieres llegar a
% t=2000, conviene aumentar h hasta 0.04, sin perder excesiva resolución o precisión en
% la curva.
%% Comentarios físicos:
% Cabe señalar que el atractor no "llena" todo el cubo del espacio de
% soluciones, sino que la región es aproximadamente plana, topología
% característica del atractor de Lorenz. De manera similar a los fractales,
% se podría definir una dimensión fraccionaria que lo describa.
% Aparte, las soluciones forman una especie de disco en cuyo centro debe
% estar la solución estacionaria, pero que está topológicamente aislada,
% por lo que no es alcanzable si el sistema ha caído en el atractor de Lorenz.
%% 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)
%%
% Se observa que el plot tiene los mismos picos que en Mathematica, como
% debería ser. Algunos valores se van ligeramente por la diferencia del
% método numérico, pero son valores sueltos y poco frecuentes respecto al comportamiento
% general.
% Cabe señalar que se escogen sólo las primeras frecuencias (de 1 a 1000),
% que es donde está lo interesante en la transformada. El hecho de que el
% espectro esté "lleno" de muchos picos es indicativo del comportamiento caótico del sistema.
%% 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)
        m=[m,lista(k)];
    end
end
L_m=length(m)

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)
%%
% Se observa que el plot es muy similar al obtenido en Mathematica. La
% posición de los picos no es exactamente la misma debido a las diferencias
% entre los distintos metódos numéricos (NDSolve frente a Runge-Kutta).
%% Comentarios físicos:
% La forma de este mapa de máximos, con un pico pronunciado, es la
% signatura de que estamos estudiando el modelo de Lorenz. En la gráfica de
% F^2(t) se observa que los consecutivos máximos de oscilación parecen
% tener una cierta relación entre ellos (véase el comportamiento autopulsado). El
% diagrama pone de manifiesto esta relación entre máximos consecutivos, con
% una curva que es característica del atractor de Lorenz.
%% Doblamiento de periodo
close all
s=3; b=1; r=140.0;
Tfin=100;%tiempo hasta el que calculas

h=0.01; t1=0:h:Tfin; %Define el paso h que utiliza Runge-Kutta
[t1,solnum1] = rungekut('laser',t1,[1,1,1]);

tmin=80;%mínimo del eje de abcisas del plot
tmax=Tfin;%máximo del eje de abcisas del plot

nmin=find(t1>=tmin,1,'first');%forma de decirle al Matlab hasta donde plotear
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')
%%
% Se aprecia que la forma es muy similar al plot de Mathematica, al estar
% en una solución estable con periodo único (solución oscilante con misma amplitud y frecuencia).
% De nuevo, elegimos un valor de plot con tmin>0 para recortar el
% transitorio del cálculo.
%%
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)

m=[];
for k=60:L-1
    if lista(k-1)<lista(k) && lista(k+1)<lista(k)
        m=[m,lista(k)];
    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}');

titul='Mapa de máximos (reescalado)';
figure('Name',titul),plot(m(1:end-1),m(2:end),'.')
title(titul);
xlabel('F^2max_i'), ylabel('F^2max_{i+1}')
axis([0,500,0,500]);
%%
% Se comprueba que hay coincidencia entre estos plots y el de Mathematica.
% A continuación, se hace un nuevo cálculo de los máximos con ajuste
% parabólico, para mejorar la dispersión que se obtiene en el mapa de
% máximos no escalado.
% Aparte, se verifica que según se va duplicando el periodo, van
% apareciendo en la transformada de Fourier pico más pequeños, lo que
% indica que la función original es periódica (pero no armónica). A partir
% de multiplicar por 8 ó 16 el periodo inicial, es difícil distinguir cada
% vuelta en el plot 2D.
% La cola que aparece al final de la transformada de Fourier es ruido
% debido ya al método numérico, por lo que es conveniente cortarlo.
%% Comentarios físicos:
% Bajando de r=140 hacia r=100, etc. iremos duplicando el periodo hasta
% llegar al comportamiento caótico del atractor (r_1=140, r_2=100,
% r_4=96.8, r_8=96.1, r_16=96.04). Este proceso se denomina
% ruta al caos por duplicación de periodo o de Feigenbaum. En la transformada de Fourier
% van apareciendo armónicos a distancias equiespaciadas (señal de movimiento uniforme), como múltiplos del
% fundamental. Esto es una prueba de que el periodo se duplica,
% cuadruplica, etc. También se puede razonar a partir de la intensidad del
% pico. Cuando bajas de r=91, pasas a la zona caótica, el periodo es
% infinito y se juntan todos los picos dando una transformada ruidosa. Si
% fuese igualmente probable (suceso azaroso), saldría una zona casi plana.
% Pero ése no es nuestro caso, sino que tenemos un comportamiento caótico,
% que muestra un cierto pico fundamental, correspondiente a lo que tarda en
% cambiar entre uno y otro atractor (las dos alas), y además modulado por las infinitas
% frecuencias.
% En r=78 se observa un nuevo atractor, pero con distinta topología. Hay 1
% general, y otro "interior" en cada solución. 
%% Cálculo con interpolación parabólica
lista=solnum1(nmin:nmax,1).^2;
L=length(lista);

m=[];
for k=60: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
close all
titul='Mapa de máximos con interpolación (reescalado)';
figure('Name',titul),plot(m(1:end-1),m(2:end),'.')
title(titul);
xlabel('F^2max_i'), ylabel('F^2max_{i+1}');
axis([395,397.2,396,397.2]);
%%
% Se ve, escogiendo la misma escala que en el mapa de máximos no escalado
% anterior, que la dispersión de los puntos es mucho menor, al haber utilizado un método
% más sofisticado de cálculo de los máximos (interpolación parabólica). Por
% tanto, el hecho de que haya dispersión de los puntos no se debe a que
% haya distintas frecuencias (desdoblamiento de periodo) sino al propio
% error numérico acumulado en el procedimiento de cálculo de máximos, que
% puede  ser más o menos preciso.

%%
% De cara a una mejora del programa, sería apropiado definir ficheros
% específicos con las funciones de búsqueda de máximos, realización de
% transformada de Fourier, etc, para llamarlos desde este programa general,
% y así evitar tener que estar copiar y pegando líneas de código. Además,
% si quisiésemos que el mapa de máximos siempre se hiciese con
% interpolación parabólica, bastaría con cambiar dicho fichero aparte y no
% cada orden si se pega específicamente en cada apartado del programa.
%%
close all