%% Generate an undecorated version of Figure 1a

function []=Figure1a()
%% Global Parameters
n=0; % number of replicas
d=5; % number of spatial dimension (set to 5 so that epsilon is normalized to 1)

% Obtain combinatorial factors
[S,a1,~]=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_oneloop(d,gs,S,a1);
        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 a dot at the Gaussian fixed points
% 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

%% 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
gI_initials=[  0.001,  0.001, 0.001,   0.001, 0.001, -0.001,   -0.001,-0.001,0];
gII_initials=[ 0.0006, 0.0011, 0.0015, 0.003, 0.022,  0.000008, 0.003, 0.025,0.0001];
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_oneloop(d,gs,S,a1);
        
        % 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
    plot3(gI_trajectory,gII_trajectory,add_height,'k');
    plot3(-gI_trajectory,-gII_trajectory,add_height,'k'); % Mirror image
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 Figure1a.eps

end