%% Generate an undecorated version of Figure 1c

function []=Figure1c()
%% Global Parameters
n=0; % number of replicas
d=5; % number of spatial dimension

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

figure();
hold on

%% Draw the colormap of beta_gI^2+beta_gII^2
gIs=(-0.5):0.01:2; NI=size(gIs,2); % gI-values
gIIs=0:0.01:4; NII=size(gIIs,2);   % gII-values
[gIs,gIIs]=meshgrid(gIs,gIIs);

% Obtain beta_gI^2+beta_gII^2
beta2=zeros(NII,NI);
for iI=1:NI
    for iII=1:NII
        gI=gIs(iII,iI); gII=gIIs(iII,iI); gs=[gI;gII];
        [betags,~,~]=betas_gammas(d,gs,S,a1,a2);
        beta2(iII,iI)=sum(betags.^2);
    end
end

% Set color
Ncolor=1024;
color_map=1-hot(Ncolor);
color=ones(NII,NI,3);
logbeta2max=log(10^8);
logbeta2min=log(10^(-2));
for iI=1:NI
    for iII=1:NII
        logbeta2=log(beta2(iII,iI));
        i_color=round(Ncolor*((logbeta2-logbeta2min)/(logbeta2max-logbeta2min)));
        if i_color<=0
            i_color=1;
        end
        color(iII,iI,1:3)=color_map(i_color,1:3);
    end
end

% Plot
surf(gIs,gIIs,beta2,color,'EdgeColor','none','LineStyle','none')

%% Put dots at the fixed points
% Find the nontrivial fixed point
gI_guess=1.17; gII_guess=1.4; % guesses for the two-loop fixed point values for d=5
[gI_star,gII_star]=fixed_point_locater(n,d,gI_guess,gII_guess);

% Plot
add_height=2*(10^8); % so that it is visible
scatter3(0,0,add_height,100,[0.9,0,0],'filled') % Gaussian fixed point
scatter3(gI_star,gII_star,add_height,100,[0,0,0.9],'filled') % Nontrivial fixed point

%% Draw RG flow trajectories
dt=0.0001; % the infinitesimal ``time" step for trajectories, where ``time $t$" is related to energy scale $mu$ via mu=mu_0*e^(-t).
max_shift=0.001; % if the shift in coupling-space is larger than this, shrink the size of the step
Nt_max=100/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

% Set initial conditions for RG flow (last two delineate the basin of attraction)
gI_initials=[0.001,  0.001,   0.001,0.001,0.001,0.001,   -0.001, -0.001,-0.001,   0,    0.001*cos(0.526367*pi/2)];
gII_initials=[0.0006,0.001128,0.00117,0.00123,0.007,0.04, 0.000008, 0.003, 0.025,0.0001,0.001*sin(0.526367*pi/2)];
N_flow=size(gI_initials,2);

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);
    
    
    while(Nt<=Nt_max && length<length_max)
        % Evaluate the beta-functions
        [betags,~,~]=betas_gammas(d,gs,S,a1,a2);
        
        % Get the infinitesimal shift and its length
        delta_gs=-dt*betags;
        infinitesimal_length=sqrt(sum(delta_gs.^2));
        
        % Adjust the length of the displacement if necessary
        if infinitesimal_length>max_shift
            delta_gs=delta_gs*max_shift/infinitesimal_length;
            infinitesimal_length=max_shift;
        end
        
        % Advance the trajectory
        gs=gs+delta_gs;
        length=length+infinitesimal_length;
        Nt=Nt+1;
        
        % Record the trajectory
        gI_trajectory(Nt)=gs(1); gII_trajectory(Nt)=gs(2);
    end
    
    % Plot
    gI_trajectory=gI_trajectory(1:10:Nt); % make it a bit sparse
    gII_trajectory=gII_trajectory(1:10:Nt); % make it a bit sparse
    add_height=(10^8)*ones(size(gI_trajectory)); % so that it is visible
    if (i_flow==N_flow) ||  (i_flow==N_flow-1)
        % Make lines thicker for the edges of the basin of attraction
        plot3(gI_trajectory,gII_trajectory,add_height,'k','LineWidth',5);
    else  
        plot3(gI_trajectory,gII_trajectory,add_height,'k');
        plot3(-gI_trajectory,-gII_trajectory,add_height,'k'); % Mirror image
    end
end

xlabel('$g^{\rm I}$','Interpreter','latex')
ylabel('$g^{\rm II}$','Interpreter','latex')
xlim([-0.5,2])
ylim([0,4])
axis square
set(gcf, 'PaperPositionMode', 'auto');
print -depsc2 Figure1c.eps

end