% (c) 2018 Aavishkar A. Patel

% This code solves the imaginary-time Dyson equations iteratively for the RM+U(1) gauge-field model
% on a discrete-time grid, by Fourier transforming back and forth between
% time and frequency representations of the two point functions.

% Gt is the time-domain (-\beta to \beta) fermion
% Green's function and Dt is the same for the boson. S is a shorthand for \Sigma,
% the self-energy, and w is a shorthand for \omega. 

% Values of couplings, number of grid points,
% and other parameters in the model. Also sets the linear update weights xf and xb.

t = 1.0; T = 1.0e-5; mu = 0.0; itermax = 600; points = 2^22; 
xf = 0.05; xb = 0.02; g = 1.0; MbyN = 0.5; mB = 1.0e-6;

% 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 (bosonic) Matsubara
% frequencies.

skipinit = false; % Set true if you want to reuse a previous solution as initial condition.

if(~skipinit)
    Gt = 0.5*sign(points+0.5-(1:2*points)); % Free fermion Green's function.
    Gw = -conj(0.5*tstep*fft(Gt));
    St = zeros(1,2*points);
    Sw = zeros(1,2*points);
    Pit = 2.0*t^2*MbyN*Gt.*fliplr(Gt); % Flipping the time domain array left to right implements \tau --> -\tau.
    Piw = conj(0.5*tstep*fft(Pit));
    Dw = zeros(1,2*points);
    Dw(2:2*points) = 1./((w1/tstep).^2/g^2+mB^2+Piw(2:2*points)-Piw(1)).*(1-mod(1:2*points-1,2));
    Dt = T*2*points*ifft(conj(Dw));
    Dthold = Dt;
end
  
% Begin iterations
    
    for n=1:itermax
        
        % Fermion Part
        St = t^2*Gt.*(1+Dt-Dt(1)); 
        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 for fermions as the matsubara frequencies are offset by a constant \pi T. 
        Gw(2:2*points) = (1./(1i*w1/tstep+mu-Sw(2:2*points))-tstep./(1i*w1)).*mod(1:2*points-1,2); % The mod restricts to fermion Matsubara frequencies.
        
        % Boson Part. Same strategy.
        Pit = 2.0*t^2*MbyN*Gt.*fliplr(Gt);
        Piw = conj(0.5*tstep*fft(Pit));
        Dw(2:2*points) = 1./((w1/tstep).^2/g^2+mB^2+Piw(2:2*points)-Piw(1)).*(1-mod(1:2*points-1,2));
        
        % Update Green's functions
        Gt = (1-xf)*Gt+xf*(0.5*sign(points+0.5-(1:2*points))-T*2*points*ifft(conj(Gw))); % Another -conj for the same reason, and linear weighted update is required for convergence.
        Dthold = (1-xb)*Dt+xb*T*2*points*ifft(conj(Dw));
        
        % Overshoot protection mechanism to avoid unstable regime and ensure convergence
        if(1-Dthold(1)<0)
            xb = 0.9*xb; % If boson overshoot detected, block the boson update and slow down its future updates, allowing the fermion to catch up.
        else
            Dt = Dthold; % If no overshoot detected, proceed as usual.
        end
        
        plot(-imag(Sw(2:2:100))); % Plot negative imaginary part of Matsubara fermion self-energy.
        pause(1/60); % Stop by using ctrl+c when converged. 
        disp([100*(1-Dt(1)) n]);
        
    end
    
Gw(2:2*points) = (1./(1i*w1/tstep+mu-Sw(2:2*points))).*mod(1:2*points-1,2);
Gw0 = Gw;
Gw0(2:2*points) = (1./(1i*w1/tstep+mu)).*mod(1:2*points-1,2);
Q = Gt(points); % Filling
FF = T*sum(log(Gw(2:2:2*points)./Gw0(2:2:2*points))) - T*log(1+exp(mu/T)) + mu*Q - T*sum(Sw.*Gw) + 0.5*t^2*T*sum(Gw.^2); % Fermion free energy.

Dw0 = Dw;
Dw0(2:2*points) = 1./((w1/tstep).^2/g^2+mB^2).*(1-mod(1:2*points-1,2));
FB = 0.5*T*sum(log(Dw0(3:2:2*points)./Dw(3:2:2*points)))/(2.0*MbyN) + T*log(1-exp(-mB*g/T))/(2.0*MbyN); % Boson free energy.

F = FF + FB; % Total free energy.