%% Draw the RG flow with minimal 2-loop RG, for d=6.005
function []=figure2c()
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.005; % number of spatial dimension
lambda1=0; lambda2=0; % coordinate deformation parameters
gI_guessS=1.4584; gII_guessS=1.4456; % guesses for the two-loop fixed point values for d=6.005
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_guessS+0.001,0,             gI_starU-0.01,  gI_starU+0.001,gI_starU-0.001, gI_starU+0.01,   gI_starU-0.01,      2]; % initial conditions for gI
gII_initials=[gII_guessS,     gII_starU+0.03,gII_starU+0.04, gII_starU+0.04,gII_starU+0.02, gII_starU+0.0716,gII_starU-0.0167833,3.45]; % initial conditions for gII
where_to_put_vectors=[5,2,0.5,1.8,1,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-3);
yy=zeros(1,N_flow-3);
vx=zeros(1,N_flow-3);
vy=zeros(1,N_flow-3);
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-2)
            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-2 || i_flow==N_flow-1
        plot(gI_trajectory,gII_trajectory,'k','LineWidth',4);
    end
    % Infinite limit cycle
    if i_flow==N_flow
        gI_LC=gI_trajectory;
        gII_LC=gII_trajectory;
        add_heights=ones(size(gI_LC));
        plot3(gI_LC,gII_LC,add_heights,'Color',[255,127,80]/255,'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', './figure2c.png');


end