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

clear;
format long
    
% PARAMETERS OF THE PROBLEM

mq = (1)^(1/2); %mass of the created particles

E0 = 10; % initial electric field
w0 = 1; %only for the sauter pulse

d = 0.1; %mode distance

t0 = 0; %if you want the Sauter pulse start at t0=-10 or so
tmax= 10;

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

 %profile (this is to save the results correctly)
  prof='';
  %prof='AC_';
  %%prof='S_';
  
  %Note: to change from one profile into another, you have to change:
  % the profile line (23)
  % line (60), initial condition for the electric field
  % line (110), definition of Eq
  %in the function fermions_2, the definition of j0

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

% 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; 
    W=(K^2+mq^2).^(1/2);
    prefactor1=((W-K)/(2*W))^(1/2);
    prefactor2=-((W+K)/(2*W))^(1/2);
    
   %c.i.
   Za(n+1)=prefactor1; %Re(h1)
   Za(n+2)=0; %Im(h1)
   Za(n+3)=prefactor2; %Re(h2)
   Za(n+4)=0; %Im(h2)
    
    end
Za(N-1)=0; %Potential vector
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); %we fix the error tolerance
RelTol=['e10e10'];

[t,y] = ode45(@(t,y)fermions_2(t,y,mq,d,Kmin,Kmax,N,E0,w0), tspan, Za,opts); 
% here, we solve the system of equations "fermions_2"
% imput: mq,d,Kmin,Kmax,N,E0,w0,tspan,Za, opts
%output; t, y (y is the vector of solutions, the first N-2 components
%correspond to the solutions for the mode equations. the component N-1 is 
% the potential vector, and the component N is - the electric field)

%-----------------------------------------%
%SOME DEFINITIONS
Sizey=size(y);

%here, we add two more components to the vector y.
y(:,N+1)=zeros(Sizey(1),1); %electric current
y(:,N+2)=zeros(Sizey(1),1); %energy density
np=zeros(Sizey(1),1);  %particle number 


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).*(y(:,n+3).^2 + y(:,n+4).^2 - (y(:,n+1).^2 + y(:,n+2).^2) - K./W );%pepe
     y(:,N+2) = y(:,N+2) + d/(2*pi).*((y(:,n+1).^2 + y(:,n+2).^2 - y(:,n+3).^2 - y(:,n+4).^2).*(K - y(:,N-1)) + 2*mq*(y(:,n+1).*y(:,n+3) + y(:,n+2).*y(:,n+4)) + W - K.*y(:,N-1)./W );%PEPECorregida

     %for the particle number
     w=((K - y(:,N-1)).^2 + mq^2).^(1/2);
     
     prefactor1=((w - (K - y(:,N-1)))./(2*w)).^(1/2);
     prefactor2=-((w + (K - y(:,N-1)))./(2*w)).^(1/2);
     Reg1=prefactor1;
     Reg2=prefactor2;
     
     np(:,1) = np(:,1) + d/(2*pi)*((y(:,n+4).^2+y(:,n+3).^2).*Reg1.^2 +(y(:,n+1).^2+y(:,n+2).^2).*Reg2.^2 - 2*(y(:,n+4).*y(:,n+2)+y(:,n+1).*y(:,n+3)).*Reg1.*Reg2 );%pepe

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
    
 
%% Other definitions
%%% spectrum %%%  


kk=zeros(1,(N-2)/4); %momentum 
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);
  
  
    %Definition of particle in t->+inf
        prefactor1= ((w - (K - y(i,N-1)))./(2*w)).^(1/2);
        prefactor2= - ((w + (K - y(i,N-1)))./(2*w)).^(1/2);
        Reg1=prefactor1;
        Reg2=prefactor2;
          
    
    kk(1,n/4+1) = K;    
    beta2(i,n/4+1) =((y(i,n+4).^2+y(i,n+3).^2)*Reg1.^2 +(y(i,n+1).^2+y(i,n+2).^2)*Reg2.^2 - 2*(y(i,n+4).*y(i,n+2)+y(i,n+1).*y(i,n+3)).*Reg1.*Reg2 );%pepe
    
    %spectrum of the electric current
    f(i,n/4+1) = -K./(2*pi).*((y(i,n+3).^2+y(i,n+4).^2) - (y(i,n+1).^2+y(i,n+2).^2) - K./W + mq^2*y(i,N-1)./(W.^3));

  end
 
end

%%Here we plot the spectrum
b=500;%with this you decide how many t_i you want to represent
  figure(6);
  
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(7);
  
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


% %% 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];
% 
%      filename=[prof,'E0=',num2str(E0),'_K=',num2str(Kmax),'_d=',num2str(d),'_mq2=',num2str(mq^2),'_',RelTol];
%      mkdir('fsolutions',filename);%crea carpeta con el nombre que le hemos dado arriba
%      filename=['fsolutions/',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'];
% 
% 
%      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);
%      
%       %% 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('fsolutions',filename);%it creates a folder with the name that we write in the line above
%      filename=['fsolutions/',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);
%    
      