%%%% 1+1 SCALAR BACKREACTION EQUATIONS %%%%
%% Silvia Pla Garcia %%

clear;
format long;
    
% PARAMETERS OF THE PROBLEM
      
  mq = (10)^(1/2); %mass of the created particles
  E0=50; % initial electric field
  w0=1;   %only for the sauter pulse
  d = 0.1;  %mode distance
  
  %note: choose E0 bigger or equal than mq^2
  
  t0 = 0;
  tmax=15;
  
  %profile 
  %(useful to save data. '' stands for the delta profile, 
  %'AC_' for the asymptotically constant profile and 'S_' for the sauter pulse)
  
  prof='';
  %prof='AC_';
  %prof='S_';

% relevant momentum
Kmin = -150;  % lower cutoff
Kmax = 100;  % upper cutoff

%-----------------------------------------%

% time interval
tspan=t0:0.01:tmax;


% Number of equations
N = 4*(Kmax - Kmin)/d + 6;


%initial conditions

    Za=zeros(1,N); %we create the vector of initial conditions
    for n=0:4:N-3 

     K=(n*d)/4+Kmin; 

    Za(n+1)=0; %Re(dh)
    Za(n+2)=-(K^2+mq^2)^(1/4); %Im(dh)
    Za(n+3)=1/(K^2 + mq^2)^(1/4); %Re(h)
    Za(n+4)=0; %Im(h)
    end

    Za(N-1)=0; %Potential vector A(0)
    Za(N)=-E0;  %initial electric field (no classical source)
    %Za(N)=0; %initial electric field (when we add a classical source)
    
%-----------------------------------------%
%ODE SOLVER. HERE, WE SOLVE THE EQUATIONS

    opts = odeset('RelTol',1e-10,'AbsTol',1e-10); %change this to improve the accuracy 
    RelTol=['e10e10'];

[t,y] = ode45(@(t,y)scalars_2(t,y,mq,d,Kmin,Kmax,N,E0,w0), tspan, Za,opts);

%-----------------------------------------%
%%Some definitions
    sizey=size(y);
    y(:,N+1)=zeros(sizey(1),1); %electric current
    y(:,N+2)=zeros(sizey(1),1); %energy density
    y(:,N+3)=zeros(sizey(1),1); %renormalized two point function
    y(:,N+4)=zeros(sizey(1),1); %unrenormalized two point function
    np=zeros(sizey(1),1);

    for n=0:4:N-3

        K=(n*d)/4 + Kmin;
        W=(K^2 + mq^2).^(1/2);

       y(:,N+1) = y(:,N+1) + d./(2*pi).*((K - y(:,N-1)).*((y(:,n+3)).^2 + (y(:,n+4)).^2) - K./(W));
       y(:,N+2) = y(:,N+2) + d./(4*pi).*((y(:,n+1)).^2 + (y(:,n+2)).^2 + ((K-y(:,N-1)).^2 + mq^2).*(y(:,n+3).^2 + y(:,n+4).^2) - 2.*W + 2*K*y(:,N-1)./(W));%Pepe
       y(:,N+3) = y(:,N+3) + d./(2*pi).*(y(:,n+3).^2 + y(:,n+4).^2 - 1./W);
       y(:,N+4) = y(:,N+4) + d./(2*pi).*(y(:,n+3).^2 + y(:,n+4).^2);
       
       %particle number
       w=((K - y(:,N-1)).^2 + mq^2).^(1/2);
       dw=-(K-y(:,N-1)).*y(:,N)./w;
  
    
    %Particle definition at t->+inf
        Reg=1./(w).^(1/2);
        Redg=0;%-dw./(2*w.^(3/2));
        Imdg=-(w).^(1/2);

     np(:,1) = np(:,1) + d/(2*pi)*1/4*(Reg.^2.*(y(:,n+2).^2+y(:,n+1).^2)+(y(:,n+3).^2+y(:,n+4).^2).*(Redg.^2+Imdg.^2)-2.*Reg.*(y(:,n+2).*y(:,n+4).*Redg - Imdg.*y(:,n+4).*y(:,n+1) + Imdg.*y(:,n+2).*y(:,n+3) + Redg.*y(:,n+1).*y(:,n+3)));
     
    end
    
    y(:,N+1) = y(:,N+1) + y(:,N-1)/(pi);
    y(:,N+2) = y(:,N+2)- y(:,N-1).^2/(2*pi);


    %we define the quantum electric field 
    
    Eq = -y(:,N); % NO classical source
    %Eq = -y(:,N) + E0./(1+t)-E0; %assymptotically constant profile  
    %Eq = -y(:,N) - E0*sech(w0*t).^2; %sauter pulse
    
    %% PLOTS

figure(1);
    
    subplot(2,1,1)
    plot(t,-y(:,N),'LineWidth',1) %electric field
    hold on
    xlabel('$$qt$$','Interpreter','latex','fontsize',16);
    ylabel('$$E/q$$','Interpreter','latex','fontsize',16);
    hold on
    

    subplot(2,1,2)
    plot(t,y(:,N+1),'LineWidth',1) %electric current
    hold on
    xlabel('$$ qt$$','Interpreter','latex','fontsize',16);
    ylabel('$$\langle J_Q\rangle/q^2 $$','Interpreter','latex','fontsize',16);
    hold on
    
figure(2)

    plot(t,2*np(:,1),'LineWidth',1)
    hold on
    xlabel('$$ qt$$','Interpreter','latex','fontsize',16);
    ylabel('$$\langle  N \rangle $$','Interpreter','latex','fontsize',16);
    hold on
    
figure(3)

    plot(t,1/2*y(:,N).^2,'LineWidth',1)
    hold on
    plot(t,y(:,N+2),'LineWidth',1)
    hold on
    plot(t,1/2*y(:,N).^2+y(:,N+2),'LineWidth',1)
    hold on
    xlabel('$$ qt$$','Interpreter','latex','fontsize',16);
    ylabel('$$ \rho /q^2 $$','Interpreter','latex','fontsize',16);
    hold on
    
  

%% spectrum 

kk=zeros(1,(N-2)/4); 
beta2=zeros(sizey(1),(N-2)/4); %n_k

for i=1:1:length(t);
    
  for n=0:4:N-3
        
     K=(n*d)/4 + Kmin; 
     W=(K^2 + mq^2).^(1/2);
     
     w=((K - y(i,N-1)).^2 + mq^2).^(1/2);
    dw=-(K-y(i,N-1)).*y(i,N)./w;
  
    
    %Particle definition at t->+inf
        Reg=1./(w).^(1/2);
        Redg=0;%-dw./(2*w.^(3/2));
        Imdg=-(w).^(1/2);
             
    
    kk(1,n/4+1) = K;   
    beta2(i,n/4+1)= 1/4*(Reg.^2.*(y(i,n+2).^2+y(i,n+1).^2)+(y(i,n+3).^2+y(i,n+4).^2).*(Redg.^2+Imdg.^2)-2.*Reg.*(y(i,n+2).*y(i,n+4).*Redg - Imdg.*y(i,n+4).*y(i,n+1) + Imdg.*y(i,n+2).*y(i,n+3) + Redg.*y(i,n+1).*y(i,n+3)));

    %spectrum of the electric current
    f(i,n/4+1) = K./(2*pi).*((K - y(i,N-1)).*((y(i,n+3)).^2 + (y(i,n+4)).^2) - K./(W) + mq^2*y(i,N-1)./(W.^3));

  end
 
end

%% Here we plot the spectrum
% the representation of the spectrum of some quantities allow us to
% determine if we have choosen the correct cut off for the mode integrals 

b=500;%with this you decide how many t_i you want to represent
  figure(4);
  
for i=1:b:sizey(1)
  plot(kk,beta2(i,:))
  hold on
  xlabel('$$ k/q$$','Interpreter','latex','fontsize',16);
  ylabel('$$|\beta|^2_k$$','Interpreter','latex','fontsize',16);
  title({'$$\langle n_k \rangle $$'},'Interpreter','latex','fontsize',16);
  
  hold on

end

figure(5)
for i=1:b:sizey(1)

  plot(kk,f(i,:))
  hold on
  xlabel('$$ k/q$$','Interpreter','latex','fontsize',16);
  ylabel('$$\langle J_k\rangle/q^2$$','Interpreter','latex','fontsize',16);
  hold on

end

%uncomment the text below to save the data...

% %% Save data in a .mat file (doesn't work in the WFU cluster)
% 
%     filename=[prof,'E0=',num2str(E0),'_K=',num2str(Kmax),'_d=',num2str(d),'_mq2=',num2str(mq^2),'_',RelTol,'.mat'];
%     save(filename);
%     
% %% Save relevant data in a tex file
% %here, I do not save the modes
%  
%      SolE=[t(:,1),-y(:,N)];
%      SolEq=[t(:,1),Eq];
%      SolA =[t(:,1),y(:,N-1)];
%      Solj =[t(:,1),y(:,N+1)];
%      SolEn =[t(:,1),y(:,N+2)];
%      Solnp =[t(:,1),2*np];
%      SolphiphiR =[t(:,1),y(:,N+3)];
%      Solphiphi =[t(:,1),y(:,N+4)];
% 
%      filename=[prof,'E0=',num2str(E0),'_K=',num2str(Kmax),'_d=',num2str(d),'_mq2=',num2str(mq^2),'_',RelTol];
%      mkdir('solutions',filename);%crea carpeta con el nombre que le hemos dado arriba
%      filename=['solutions/',prof,'E0=',num2str(E0),'_K=',num2str(Kmax),'_d=',num2str(d),'_mq2=',num2str(mq^2),'_',RelTol];
% 
% 
%      filenameE=[filename,'/E.txt'];
%      filenameEq=[filename,'/Eq.txt'];
%      filenameA=[filename,'/A.txt'];
%      filenamej=[filename,'/j.txt'];
%      filenameEn=[filename,'/En.txt'];
%      filenamenp=[filename,'/np.txt'];
%      filenamephiphiR=[filename,'/phiphiR.txt'];
%      filenamephiphi=[filename,'/phiphi.txt'];
% 
% 
%      fileID = fopen(filenameE,'w');
%      fprintf(fileID,'%.14f %.14f\n',transpose(SolE));
%      fclose(fileID);
% 
%      fileID = fopen(filenameEq,'w');
%      fprintf(fileID,'%.14f %.14f\n',transpose(SolEq));
%      fclose(fileID);
% 
%      fileID = fopen(filenameA,'w');
%      fprintf(fileID,'%.14f %.14f\n',transpose(SolA));
%      fclose(fileID);
% 
%      fileID = fopen(filenamej,'w');
%      fprintf(fileID,'%.14f %.14f\n',transpose(Solj));
%      fclose(fileID);
%      
%       fileID = fopen(filenameEn,'w');
%      fprintf(fileID,'%.14f %.14f\n',transpose(SolEn));
%      fclose(fileID);
%      
%       fileID = fopen(filenamenp,'w');
%      fprintf(fileID,'%.14f %.14f\n',transpose(Solnp));
%      fclose(fileID);
%      
%       fileID = fopen(filenamephiphiR,'w');
%      fprintf(fileID,'%.14f %.14f\n',transpose(SolphiphiR));
%      fclose(fileID);
%      
%      fileID = fopen(filenamephiphi,'w');
%      fprintf(fileID,'%.14f %.14f\n',transpose(Solphiphi));
%      fclose(fileID);
%       
%      
%      %% Save modes in a tex file
%      
%      Sol=y;
%      filename=[prof,'E0=',num2str(E0),'_K=',num2str(Kmax),'_d=',num2str(d),'_mq2=',num2str(mq^2),'_',RelTol];
%      mkdir('solutions',filename);%crea carpeta con el nombre que le hemos dado arriba
%      filename=['solutions/',prof,'E0=',num2str(E0),'_K=',num2str(Kmax),'_d=',num2str(d),'_mq2=',num2str(mq^2),'_',RelTol];
% 
%      
%      filenamey=[filename,'/ymod.txt'];
%      filenamet=[filename,'/t.txt'];
%         
%      fileID = fopen(filenamet,'w');
%      fprintf(fileID,'%.14f\n',transpose(t));
%      fclose(fileID);
%      
%      
%        s0="%.14f ";
%     s1=s0;
%     
%     for i=1:1:N-2
%         s1=strcat(s1,s0);
%     end
%         s1=strcat(s1,'%.14f\n');
%       
%       fileID = fopen(filenamey,'w');
%       fprintf(fileID,s1,transpose(Sol));
%      fclose(fileID);
%    