% (c) 2017 Aavishkar A. Patel

% This code solves the imaginary time Dyson equations iteratively for the c-f-SYK model
% on a discrete-time grid, by Fourier transforming back and forth between
% time and frequency representations of the two point functions.

% Run, and hold down return until the plot converges.
% You can choose what to plot where the plot command appears.
% Gt is the time-domain (-\beta to \beta) SYK
% Green's function and Gct is the same for Gc. S is a shorthand for \Sigma,
% the self-energy, and w is a shorthand for \omega. 

% Values of couplings, bandwidth, number of grid points,
% and other parameters in the model. Also sets the linear update weight x.

J = 100.0; T = 1.0; mu = 0.0; itermax = 1000; points = 2^20; x = 0.05;
g = 80.0; Lambda = 3.0; muc = 0.0*Lambda; nu0 = 2.0*pi/Lambda; MbyN = 0.5;

% Appropriate time and frequency steps. Set up arrays for time and
% frequency domains. The time domain array is directly -\beta to \beta,
% left to right. The frequency domain array stores positive frequencies
% in the left half, going left to right in increasing magnitude,
% and negative frequencies on the right end, going right to left in increasing
% magnitude.
 
tstep = 1/(points*T);
wstep = pi/points;
w1 = zeros(1,2*points-1);
w1(1:points)=wstep*(1:points);
w1(points:2*points-1) = wstep*((points:2*points-1)-2*points);

% Initialize green's functions and self energies in both representations. 
% When the initialization is nonzero, it contains only fermionic Matsubara
% frequencies.

Gt = 0.5*sign(points+0.5-(1:2*points));
Gw = -conj(0.5*tstep*fft(Gt));
St = zeros(1,2*points);
Sw = zeros(1,2*points);
Gct = 0.5*sign(points+0.5-(1:2*points));
Gcw = -conj(0.5*tstep*fft(Gct));
Sct = zeros(1,2*points);
Scw = zeros(1,2*points);
  
% Begin iterations
    
    for n=1:itermax
        
        % SYK Part
        St = -J^2*(Gt.^2).*fliplr(Gt) - MbyN*g^2*Gt.*Gct.*(fliplr(Gct)); % Flipping the time domain array left to right implements \tau --> -\tau
        Sw = -conj(0.5*tstep*fft(St)); % -Conj is used because MATLAB's fft has the exponent convention +i\omega\tau, but we want -i\omega\tau for the forward transform. Also, since the time-domain is doubled (i.e. -\beta to \beta), the - sign is required as the matsubara frequencies are offset by a constant \pi T. 
        Gw(2:2*points) = 1./(-(exp(-1i*w1)-1)/tstep+mu-Sw(2:2*points)).*mod(1:2*points-1,2); % The i\omega term has to be "lattice-regulated" in order to get the correct values of G close to \tau = 0. The mod restricts to fermion Matsubara frequencies.
        Gt =  (1-x)*Gt-x*T*2*points*ifft(conj(Gw)); % Another -conj for the same reason, and linear weighted update is required for convergence.
        
        % Itinerant Part. Same strategy.
        Sct = -g^2*Gct.*Gt.*fliplr(Gt);
        Scw = -conj(0.5*tstep*fft(Sct));
        Gcw(2:2*points) = (nu0/(2.0*pi))*(log(Lambda+2*muc-2*(exp(-1i*w1)-1)/tstep-2*Scw(2:2*points))-log(2*muc-Lambda-2*(exp(-1i*w1)-1)/tstep-2*Scw(2:2*points))).*mod(1:2*points-1,2); % This is the finite-bandwidth equation from the paper.
        Gct = (1-x)*Gct-x*T*2*points*ifft(conj(Gcw));
        
        plot(real(Gct)); % Plot what you want here
        b = waitforbuttonpress; % You'll need to hold down a key to make the iteration go. Release when it converges. 
        
    end