%% Práctica sobre el MODELO DE LORENZ-HAKEN (1)
% Programa traducido de Mathematica a Matlab  y comentado 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];

% Estabilidad lineal (hasta bifurcación de Hopf rHB)
rHB=-s*(3+b+s)/(1+b-s);

%% COMPROBACIÓN DE:
%% a) Estabilidad lineal (condiciones iniciales próximas a la solución estacionaria)
% Evaluación para valores concretos
RHB=subs(rHB,{s,b},[3,1])
Solestac=subs(solestac,r,22)
%% Comentarios físicos:
% La solución estacionaria no es la única solución "estable" en la que el láser puede funcionar.
% Para r grande, mayor que la bifurcación de Hopf, tendremos soluciones
% caóticas en general, con posibles islas de estabilidad y soluciones
% oscilantes de frecuencia fija (Ver doblamiento de periodo el segundo
% archivo). También podremos abandonar la solución estacionaria si hacemos
% un encendido brusco (Hard-mode excitation), con condiciones iniciales
% lejanas a las de la solución estacionaria, de forma que ya no la "encuentra
% (se puede hacer la analogía mecánica de que intentas encestar una pelota
% en un pequeño pozo de potencial justo en la cúspide de una montaña). Por contra, si ponemos
% condiciones iniciales cercanas a la misma, ese pequeño ruido (cambio adiabático) no es
% suficiente para sacarlo de ella y cae a la línea de estabilidad.
%% 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; %no lo defino como vector "par" sino como variable global para ahorrar notación
Tfin=100;%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,[4.58258,4.58258,1]);
%%
% solnum es una matriz con tres columnas ordenadas (F,P,d). Las filas marcan
% distintos valores del tiempo. Se fijan los valores iniciales [F0,P0,d0]=[4.58258,4.58258,1] por ejemplo.
% Se utiliza el método de Runge-Kutta porque el de ODE45 no se comporta
% adecuadamente, es decir, no converge a la solución estable sino que es
% infinitamente oscilante para r=20, s=3, b=1 debido a error numérico.
%%
tmin=0;%mínimo del eje de abcisas del plot
tmax=Tfin;%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='Estabilidad lineal';
figure('Name',titul),plot(t(nmin:nmax), solnum(nmin:nmax,1).^2);
%el 1 representa la primera columna, es decir, la F
axis tight
title(titul)
xlabel('t')
ylabel('F^2')

%%
% Como se observa, el transitorio hasta llegar a la solución final
% estacionaria es muy largo. Si quisiésemos centrarnos en la estacionaria,
% deberíamos poner Tfin grande (500 aprox.) y establecer tmin=400. De esta
% manera podemos cortar el transitorio de los cálculos. Se definen en todos
% los plots estas variables para poder prescindir de los transitorios en
% caso de que no nos interesasen.

%% b) Hard-mode excitation -> biestabilidad generalizada
s=3; b=1; r=20;
Tfin=100;
h=0.01; t=0:h:Tfin;
[t,solnum] = rungekut('laser',t,[0.001,0.001,1]);%orden F0,P0,d0
close all
titul='Hard-mode excitation'
figure('Name',titul),plot(t, solnum(:,1).^2);
axis tight
title(titul)
xlabel('t')
ylabel('F^2')
%% Comentarios físicos:
% El Hard-mode excitation se alcanza dando unos valores iniciales alejados
% de la solución estacionaria. De esta manera, el sistema no es capaz de
% "encontrar" dicha solución estacionaria y cae al atractor caótico y
% permanece en ella de manera estable (ya no puede volver a la
% estacionaria). A diferencia de a), donde s,b y r son compartidos, el valor
% de F0 y P0 se cambia sustancialmente para sacar al sistema de la solución
% estacionaria. Para F0=P0=4.5826 (caso a), que son los valores de la solución
% estacionaria, el sistema es capaz de alcanzar la solución estacionaria
% pese al alto bombeo, ya cercano a rHB (bifurcación de Hopf), si se deja
% pasar suficientemente tiempo (t=300 ó 400). Por contra, en el caso b),
% cambiamos notablemente dichos valores iniciales, y provocamos el hard mode
% excitation, es decir, el salto a la solución caótica pese a estar por
% debajo de rHB, dado que la perturbación inicial es grande y se va amplificando cada oscilación en lugar de irse atenuando hacia la estable. Esta es la principal diferencia entre el cambio adiabático
% y el hard mode excitation. En otras palabras, si damos un "golpe" muy fuerte al láser
% (valores iniciales lejanos a la estacionaria) o subimos de golpe el bombeo
% manteniendo valores iniciales de otra solución, podemos salirnos de la
% solución estacionaria y caer a la caótica. De ahí la importancia del
% cambio adiabático (subir poco a poco, con condiciones iniciales siempre cercanas a la estacionaria) para mantenerse en la estacionaria
% hasta rHB. 
%%
% Cabe resaltar la diferencia entre estacionario y estable. La solución
% estacionaria indica que el valor de la variable es fijo en el tiempo, no
% cambia (derivada nula). La solución estacionaria es única y está
% calculada analíticamente (ver solestac). El resto de soluciones dinámicas
% del sistema son no estacionarias. Pero si el láser se mantiene en un
% subconjunto de soluciones (no estacionarias), oscilando por ejemplo alrededor de un valor medio, entonces decimos que la
% solución es estable. Es el caso de la forma autopulsada estable de la
% curva correspondiente al hard-mode excitation.
%% Representación de atractores
s=3; b=1; r=30;
Tfin=100;
h=0.01; t=0:h:Tfin;
[t,solnum] = rungekut('laser',t,[0.001,0.001,1]);%orden F0,P0,d0

close all
titul='Representación de atractores 2D';
figure('Name',titul),plot(solnum(:,1), -solnum(:,3));
axis tight
title(titul)
xlabel('F')
ylabel('-d')

titul='Representación de atractores 3D';
figure('Name',titul),plot3(solnum(:,1), solnum(:,2),-solnum(:,3));
axis tight
title(titul)
xlabel('F')
ylabel('P')
zlabel('-d')
%% Comentarios físicos:
% Aparecen dos "alas de mariposa" o cuencas de atracción, es decir, las
% soluciones dinámicas del sistema dan vueltas alrededor de una zona del
% espacio de soluciones. Esto es característico del atractor de Lorenz. El
% hecho de que aparezcan dos puntos de atracción se debe a que en los
% casos 2D y 3D se representa F y no F^2 como en el caso 1D. Si en el 1D
% representásemos F, se vería que saltaría a valores negativos y positivos
% aleatoriamente.
% Cabe señalar que la notación 1D, 2D y 3D se refiere a: 1D-> F(t);
% 2D->Representación paramétrica F(t), -d(t); 3D -> Paramétrica 3D F(t),
% P(t), -d(t).
%% Comprobación Sensibilidad a las Condiciones Iniciales
s=3; b=1; r=23;
Tfin=50;
h=0.01; t2=0:h:Tfin;
[t2,solnum2] = rungekut('laser',t2,[1,1,1]);%orden F0,P0,d0
% Evaluar los valores de F, P y D para el último t calculado:
F0=solnum2(end,1)
P0=solnum2(end,2)
d0=solnum2(end,3)

Tfin=20;
h=0.01; t3=0:h:Tfin; t4=0:h:Tfin;
[t3,solnum3] = rungekut('laser',t3,[F0,P0,d0]);
[t4,solnum4] = rungekut('laser',t4,[F0+0.01,P0+0.01,d0]);

close all
titul='Sensibilidad a las condiciones iniciales (s3)';
figure('Name',titul),plot(t3, solnum3(:,1).^2);
axis tight
title(titul)
xlabel('t')
ylabel('F^2')

titul='Sensibilidad a las condiciones iniciales (s4)';
figure('Name',titul),plot(t4, solnum4(:,1).^2);
axis tight
title(titul)
xlabel('t')
ylabel('F^2')

titul='Sensibilidad a las condiciones iniciales (3vs4)';
figure('Name',titul),plot(t3, solnum3(:,1).^2, t4, solnum4(:,1).^2);
axis tight
title(titul)
xlabel('t')
ylabel('F^2')
legend('solnum3','solnum4')

tmin=0;%mínimo del eje de abcisas del plot
tmax=2;%máximo del eje de abcisas del plot

nmin=find(t3>=tmin,1,'first');%forma de decirle al Matlab hasta donde plotear
nmax=find(t3<=tmax,1,'last');

titul='Sensibilidad a las condiciones iniciales (s3^2-s4^2)';
figure('Name',titul),plot(t3(nmin:nmax), solnum3(nmin:nmax,1).^2 - solnum4(nmin:nmax,1).^2);
axis tight
title(titul)
xlabel('t')
ylabel('F_3^2-F_4^2')

%%
% Cabe señalar que estos últimos plots (y los de aquí en adelante) son ligeramente distinto al de Mathematica. La
% razón es que utilizamos un valor de F0, P0 y d0 calculado mediante
% Runge-Kutta que es notablemente distinto al obtenido mediante
% NDSolve. Esta diferencia se debe a que estamos en una zona caótica
% (solución no es estacionaria), con lo que cualquier mínima diferencia
% entre uno y otro método numérico genera variaciones grandes en la curva y
% por tanto en el valor de F, P y d en Tfin. Si se representa solnum2 en
% Matlab y Mathematica, se aprecia que su aspecto es similar pero desplazado en el tiempo, y a
% partir de cierto valor de t, el orden de picos es difícil seguirlo. De hecho, esta observación 
% constituye una prueba más de la sensibilidad a las condiciones iniciales.

%%
close all
titul='Sensibilidad a las condiciones iniciales 3D (s3)';
figure('Name',titul),plot3(solnum3(:,1), solnum3(:,2),-solnum3(:,3));
axis tight
title(titul)
xlabel('F')
ylabel('P')
zlabel('-d')

titul='Sensibilidad a las condiciones iniciales 3D (s4)';
figure('Name',titul),plot3(solnum4(:,1), solnum4(:,2),-solnum4(:,3));
axis tight
title(titul)
xlabel('F')
ylabel('P')
zlabel('-d')

titul='Sensibilidad a las condiciones iniciales 3D (3vs4)';
figure('Name',titul),plot3(solnum3(:,1), solnum3(:,2),-solnum3(:,3), solnum4(:,1), solnum4(:,2),-solnum4(:,3));
axis tight
title(titul)
xlabel('F')
ylabel('P')
zlabel('-d')
legend('solnum3','solnum4')
%%
close all
%% Comentarios físicos:
% Se aprecia que el sistema tiene un comportamiento caótico debido a la
% gran sensibilidad en las condiciones iniciales. Pese a poner una
% perturbación muy pequeña entre dos soluciones, ambas acaban divergiendo
% si se deja pasar el tiempo necesario (más tiempo cuanto menor sea la
% perturbación). Esto es característico de los sistemas caóticos clásicos (caos determinista), donde la
% incertidumbre en las condiciones iniciales te impiden predecir el devenir
% del sistema. Cualquier pequeña variación o ruido va a provocar que el
% cabo de un tiempo largo dé una solución radicalmente diferente a la
% obtenida sin dicha variación.