% Solves the nonlinear model for the L-H transition proposed in Diamond PRL 1994

iter=10000; % 
dt=0.01; % 

% Parameters 
global a b;
a = 1; 
b = 0;

% Initial conditions
X0 = [0.1,0.5]; % 

% Solution with ode45
[t,X] = ode45(@diffEq,[0,iter*dt],X0');
% X = RK4(X0,dt,iter);
% t = linspace(0,dt*iter,iter+1);

E = X(:,1);
U = X(:,2);

%% Plot evolution (f)

 figure
 plot(t,E,'r-','linewidth',2)
 hold on
 plot(t,U,'b-.','linewidth',2)
 legend('E','U')
 xlabel('\tau')
 title (['a = ',sprintf('%2.1f',a)])
 ylim([-0.1 1.4])
%  


% Save the figure in pdf format
% h = gcf;
% h.PaperPositionMode = 'auto';
% set(gca,'FontSize',16)
% fig_pos = h.PaperPosition;
% fig_pos(2)=fig_pos(2)-18.;
% h.PaperSize = [fig_pos(3) fig_pos(4)*1.1];
% set(gcf,'Units','Inches');
% pos = get(h,'Position');
% set(h,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)-0.5])
% cd ../../ExoFig/
% print(h,['lor_attr_dim_r',num2str(r)],'-dpdf','-r200','-bestfit')
% if a==0.5 && b==1 && iter==10000
%     print(h,'L_H_transition_L','-dpdf','-r200')
% elseif a==1.2 && b==1 && iter==10000
%     print(h,'L_H_transition_H','-dpdf','-r200')
% elseif a==3.0 && b==1 iter==10000
%     print(h,'L_H_transition_H_osc','-dpdf','-r200')
% end

% cd ../ExoNumericalFiles/L-H_transition/

%% Quiver plot 
 

% Setting up the grids
xinf = [0.1 0.1];
xsup = [2.1 2.1];
n = 16;
x1 = linspace(xinf(1),xsup(1),n);
x2 = linspace(xinf(2),xsup(2),n);
[X1,X2] = meshgrid(x1,x2);


% Computing the vector field of diffEq
dX1 = zeros(size(X1));
dX2 = zeros(size(X1));
t0=0;

%%

for i = 1:numel(X1)	% Loop on all elements of 2D matrix
    dX = diffEq(t0,[X1(i); X2(i)]);
    dX1(i) = dX(1);
    dX2(i) = dX(2);
end

% 
% % Normalisation (if desired)
% norm_dX1 = dX1./sqrt(dX1.^2+dX2.^2);
% norm_dX2 = dX2./sqrt(dX1.^2+dX2.^2);


% Quiver plot
figure;
quiver(X1,X2,0.7,'k');
quiver(X1,X2,dX1,dX2,2.0,'b');
title('$\dot{\vec{x}}=(\cos x_1+\sin x_2,\cos x_1)$','Interpreter','latex');
grid on;
xlabel('E');
ylabel('U');
title(['a = ',num2str(a), ' , b = ',num2str(b)])
axis equal; %tight;
xlim([xinf(1),xsup(1)]); ylim([xinf(2),xsup(2)]);
set(gca,'FontSize',14);
box on


%% Stream plot (g)
% figure
hold on; % To plot multiple trajectories
tf=20; % Final time
% gridsize = 21; % Here the grid of initial conditions can be of different density than the one for quiver
[X1ini,X2ini] = meshgrid(linspace(xinf(1),xsup(1),n),linspace(xinf(2),xsup(2),n)); % Grid of initial conditions
for i=1:numel(X1ini) % Loop on the 2D matrix of initial conditions
    [~,Xsol] = ode45(@diffEq,[t0,tf],[X1ini(i),X2ini(i)]); % Integration
    if mod(i,2)==0  % For good looks, can alternate colors
    	color='r';
    else
    	color='k';
    end
    plot(Xsol(:,1),Xsol(:,2),color) % Trajectory
    hold on
    plot(Xsol(end,1),Xsol(end,2),'bo','markersize',10) % Ending point
end
xlim([-0.1 xsup(1)]); ylim([0 xsup(2)]);
hold off;
xlabel('E');
ylabel('U');
% % 

% Save the figure in pdf format
% h = gcf;
% h.PaperPositionMode = 'auto';
% fig_pos = h.PaperPosition;
% fig_pos(2)=fig_pos(2)-18.;
% h.PaperSize = [fig_pos(3) fig_pos(4)*1.1];
% set(gcf,'Units','Inches');
% pos = get(h,'Position');
% set(h,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)-0.5])
% cd ../../ExoFig/
% % print(h,['lor_attr_dim_r',num2str(r)],'-dpdf','-r200','-bestfit')
% if a==1 && b==0  
%     print(h,'phase_space_trajectory_a1_b0','-dpdf','-r200')
% elseif a==3 && b==1
%     print(h,'phase_space_trajectory_a3_b1','-dpdf','-r200')
% end
% 
% cd ../ExoNumericalFiles/L-H_transition/

% %% Plot for all a
% a_vec = linspace(0,3,31);
% E_fin = zeros(1,length(a_vec));
% U_fin = zeros(1,length(a_vec));
% E_ana = zeros(1,length(a_vec));
% U_ana = zeros(1,length(a_vec));
% 
% for i = 1:length(a_vec)
%   [t,X_dum] = ode45(@(t,X) diffEq_param(t,X,a_vec(i),b),[0,iter*dt],X0');
%   E_fin(i)=X_dum(end,1);
%   U_fin(i)=X_dum(end,2);
% end
% 
% for i = 1:length(a_vec)
%     if a_vec(i)<=b
%         E_ana(i) = 1;
%         U_ana(i) = 0;
%     else
%         E_ana(i) = b/a_vec(i);
%         U_ana(i) = (a_vec(i)-b)/a_vec(i);
%     end
% end
% 
% figure
% plot(a_vec,E_fin,'k')
% hold on
% plot(a_vec,E_ana,'k--o')
% xlabel('a')
% ylabel('E')
% 
% figure
% plot(a_vec,U_fin,'b')
% hold on
% plot(a_vec,U_ana,'bo')
% xlabel('a')
% ylabel('U')

