%% Draw the RG flow with minimal 2-loop RG, for d=6.05
function []=figure2d()
dt=0.0001; % 0.0001 timestep in production (overkill)
n=0; % number of replicas
direction_sign=(-1); % forward=IR=-1, backward=UV=+1
d=6.05; % number of spatial dimension
lambda1=0; lambda2=0; % coordinate deformation parameters
gI_guessS=1.4709; gII_guessS=1.4529; % guesses for the two-loop fixed point values for d=6.05
gI_starU=0; gII_starU=sqrt((288/1931)*(-3/2+sqrt((9/4)+(d-6)*(1931/288)))); % exact values for the fixed point ``U" (in Moore-Read notation)
gI_guessZ=0.399*sqrt(d-6); gII_guessZ=0.028*sqrt(d-6); % guess for the fixed point ``Z" (in Moore-Read notation)
gI_initials= [gI_starU-0.01,  gI_starU+0.001,gI_guessZ+0.01, gI_guessS-0.01,  gI_guessS-0.01, 0,             gI_starU-0.01,           gI_guessS-0.01,     gI_guessS-0.01     ]; % initial conditions for gI
gII_initials=[gII_starU+0.04, gII_starU+0.02,gII_guessZ,     gII_guessS-0.003,gII_guessS-0.03,gII_starU+0.03,gII_starU-0.018248498465,gII_guessS-0.004544,gII_guessS-0.00015]; % initial conditions for gII
where_to_put_vectors=[0.5,2,1,3,3.5,100,100,100,100];
Nt_max=500/dt; % stop the trajectory when number of step reaches this number
length_max=100; % stop the trajectory when the length of the trajectory reaches this number


%% Plot
close all
figure()
hold on

% Get the combinatorial factors
[S,a1,a2]=combinatorial_factors(n);

% Plot Gaussian fixed point
scatter3(0,0,1,100,[1,0,0],'filled')

% Plot the unstable fixed point, ``U"
scatter3(gI_starU,gII_starU,1,100,[0,0.8,0],'filled')

% Find and plot the nontrivial fixed point
[gI_starS,gII_starS,~]=fixed_point_locater_with_Hessian(n,d,gI_guessS,gII_guessS,lambda1,lambda2);
scatter3(gI_starS,gII_starS,1,100,[0,0,1],'filled')

% Find and plot the two-direction-unstable fixed point, ``Z," that is visible already at one-loop
[gI_starZ,gII_starZ,~]=fixed_point_locater_with_Hessian(n,d,gI_guessZ,gII_guessZ,lambda1,lambda2);
scatter3(gI_starZ,gII_starZ,1,100,[0,0.8,0],'filled')

% RG flow
N_flow=size(gI_initials,2);
xx=zeros(1,N_flow-4);
yy=zeros(1,N_flow-4);
vx=zeros(1,N_flow-4);
vy=zeros(1,N_flow-4);
i_vec=0;
for i_flow=1:N_flow
    % Set the initial condition
    gs=[gI_initials(i_flow);gII_initials(i_flow)];
    length=0;
    Nt=1;
    gI_trajectory=gs(1); gII_trajectory=gs(2);
    length_vercor=where_to_put_vectors(i_flow);
    put_vector_flag=0;
    
    while(Nt<=Nt_max && length<length_max)
        gs_old=gs;
        
        %% Forth-order Runge-Kutta
        gs=gs_old;
        [betags,~,~]=betas_gammas(d,gs,S,a1,a2,lambda1,lambda2);
        k1=direction_sign*betags;
        gs=gs_old+dt*(k1/2);
        [betags,~,~]=betas_gammas(d,gs,S,a1,a2,lambda1,lambda2);
        k2=direction_sign*betags;
        gs=gs_old+dt*(k2/2);
        [betags,~,~]=betas_gammas(d,gs,S,a1,a2,lambda1,lambda2);
        k3=direction_sign*betags;
        gs=gs_old+dt*k3;
        [betags,~,~]=betas_gammas(d,gs,S,a1,a2,lambda1,lambda2);
        k4=direction_sign*betags;
        delta_gs=dt*(k1+2*k2+2*k3+k4)/6;
        infinitesimal_length=sqrt(sum(delta_gs.^2));      
        
        %% Advance the trajectory
        gs=gs_old+delta_gs;
        length=length+infinitesimal_length;
        Nt=Nt+1;

        % Record the trajectory
        gI_trajectory(Nt)=gs(1); gII_trajectory(Nt)=gs(2);
        
        % Put flow-vector on the curve, once
        if (length>length_vercor) && (put_vector_flag==0) && (i_flow<N_flow-3)
            i_vec=i_vec+1;
            put_vector_flag=1;
            xx(i_vec)=gs(1);
            yy(i_vec)=gs(2);
            vx(i_vec)=delta_gs(1)/infinitesimal_length;
            vy(i_vec)=delta_gs(2)/infinitesimal_length;
        end
        
        
    end
    plot(gI_trajectory,gII_trajectory,'k','LineWidth',2);
    % Basin of attraction
    if i_flow>=N_flow-3
        plot(gI_trajectory,gII_trajectory,'k','LineWidth',4);
    end
    xlim([-0.5,2])
    ylim([0,4])
    drawnow;
end

% Put a directional vector on each flow
quiver(xx,yy,vx,vy,0.2,'k','LineWidth',3)

xlim([-0.5,2])
ylim([0,4])
xlabel('$g^{\mathrm{I}}$','Interpreter','latex');
ylabel('$g^{\mathrm{II}}$','Interpreter','latex')
print(gcf, '-dpng', './figure2d.png');


end