% (c) 2018 Aavishkar A. Patel

% This code solves the zero temperature, real-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 spectral functions.

% Values of couplings, number of grid points,
% and other parameters in the model. Also sets the update parameter xf.
% The exponent x is inputted from the analytic solution in order to provide
% an initial condition that is close to the actual solution, otherwise the
% numerical procedure does not converge.

t = 1e0; itermax = 1e4; points = 2^18; xf = 1e-3;
g = 1e0; MbyN = 0.5; pad = 1/2; x = 0.23; mu = 0.0; mB = 0.0;

% Appropriate time and frequency steps. Set up arrays for time and
% frequency domains. The time domain array is 0 to points/pad*tstep, left to right,
% with the remainder containing negative time and zero padding. 
% 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. The effective dimensionful frequency step is wstep/tstep.
 
tstep = 1/2^4;
wstep = pi/(2.0*points);

w1 = zeros(1,4*points);
w1(1:2*points)=wstep*(0:2*points-1);
w1(2*points+1:4*points) = wstep*((2*points:4*points-1)-4*points);

% T = 0 Bose and Fermi functions.

nB = zeros(1,4*points);
nB(1)=-0.5;
nB(2:2*points)=0;
nB(2*points+1:4*points)=-1;
nF = zeros(1,4*points);
nF(1)=0.5;
nF(2:2*points)=0;
nF(2*points+1:4*points)=1;

% Initialize some guesses for the spectral functions. The fermion spectral
% function must be close to the actual answer, or else the program does not
% converge.

A = 1./(0.001+0.001*abs(w1/tstep).^x+abs(w1/tstep).^2);
A = A/((wstep/tstep)*(1/(2.0*pi))*sum(A));
n2 = (wstep/tstep)*(1/(2.0*pi))*fft(A.*nF);
n2 = n2(1:points/pad);
B = zeros(1,4*points);

% Nonlinear update. The singular low frequency regions must be updated
% slower than the high frequency regions, or else the program does not
% converge.

xupd = xf+(1-xf)*abs(tanh(w1/tstep)).^(1-2*x);
        
    for n=1:itermax
        
        n2 = (wstep/tstep)*(1/(2.0*pi))*fft(A.*nF);
        n2s = (wstep/tstep)*(1/(2.0*pi))*fft(A.*(1-nF));
        n2 = n2(1:points/pad);
        n2s = n2s(1:points/pad);
        PiR = 4.0*t^2*MbyN*tstep*4.0*points*ifft(imag(conj(n2).*n2s),4*points);
        PiR = PiR-real(PiR(1));
              
        B(2:4*points) = 2*imag(1./((w1(2:4*points)/tstep).^2/g^2+mB^2+PiR(2:4*points)));
       
        n3 = (wstep/tstep)*(1/(2.0*pi))*fft(B.*nB);
        n1 = (wstep/tstep)*(1/(2.0*pi))*sum(B.*nB);
        n3 = n3(1:points/pad);

        ISwr = imag(-1i*t^2*tstep*4.0*points*ifft(n2s.*conj(n3)+n2.*n3,4*points))-0.5*(1-n1)*t^2*A;
        
        A = (1-xupd).*A-2*xupd.*imag(1./(w1/tstep+mu-1i*hilbert(ISwr)));
        
        plot([A(4*points-2^12+2:4*points) A(1:2^12)]) % Plot fermion spectral function.
        pause(1/60); % Stop by using ctrl+c when converged.
        disp([1-n1 n]);
        
    end