clear all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Program to calculate the stability diagram in orthodox theory
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Vg_lim = [-3.41 3.41]*10^-3; %eV
Vb_lim = [-4 4]*10^-3; %eV

U = 2.34*10^-3; %eV

gamma = [1 1]; % 1/hbar

T = (10)*8.6*10^-5; %(K) eV

finesse = 200;

Vg_vec = linspace(Vg_lim(1),Vg_lim(2),finesse);
Vb_vec = linspace(Vb_lim(1),Vb_lim(2),finesse);

for k =1:length(Vg_vec);
    for j = 1:length(Vb_vec);
        Vb = Vb_vec(j); % V
        
        eps = 10^-6;
        
        VL=+Vb/2; % V       
        VR=-Vb/2; % V
        
        Vg=Vg_vec(k);% V
  
        
        
        n0 = -Vg/U;
       
        GammaL_in = 0;
        GammaL_out = 0;
        
        P_pos = 0;
        P_temp = 1;
        
        i = 1;
        while P_temp>eps
 
            DeltaE_L = U*(floor(n0)+i-1) + Vg - VL;
            DeltaE_R = U*(floor(n0)+i-1) + Vg - VR;
            GammaL_in(i)  = 2*gamma(1)*(+DeltaE_L)/(exp(DeltaE_L/T)-1);
            GammaR_in     = 2*gamma(2)*(+DeltaE_R)/(exp(DeltaE_R/T)-1);
            GammaL_out(i) = 2*gamma(1)*(-DeltaE_L)/(exp(-DeltaE_L/T)-1);
            GammaR_out    = 2*gamma(2)*(-DeltaE_R)/(exp(-DeltaE_R/T)-1);
            
            Gamma_in = GammaL_in(i)+ GammaR_in; 
            Gamma_out = GammaL_out(i) + GammaR_out;
            
            P_temp = P_temp*Gamma_in/Gamma_out;
            
            P_pos(i) = P_temp;
            i = i+1;
        end
        
        I_pos = [1,P_pos(1:end-1)]*GammaL_in' - P_pos*GammaL_out';
        
        GammaL_in = 0;
        GammaL_out = 0;
                
        P_neg = 0;
        P_temp = 1;
        i = 1;
        while P_temp>eps
            
            DeltaE_L = U*(floor(n0)-i) + Vg - VL;          
            DeltaE_R = U*(floor(n0)-i) + Vg - VR;
            
            GammaL_in(i) = 2*gamma(1)*(+DeltaE_L)/(exp(DeltaE_L/T)-1); %2eV/h
            GammaR_in    = 2*gamma(2)*(+DeltaE_R)/(exp(DeltaE_R/T)-1);
            GammaL_out(i)= 2*gamma(1)*(-DeltaE_L)/(exp(-DeltaE_L/T)-1);
            GammaR_out   = 2*gamma(2)*(-DeltaE_R)/(exp(-DeltaE_R/T)-1);
            
            Gamma_in = GammaL_in(i)+GammaR_in;
            Gamma_out = GammaL_out(i) + GammaR_out;
            
            P_temp = P_temp*Gamma_out/Gamma_in;
            P_neg(i) = P_temp;
            i = i+1;
        end

        I_neg = P_neg*GammaL_in' - [1,P_neg(1:end-1)]*GammaL_out';
 
        P = [fliplr(P_neg),1,P_pos];

        Unitconv = 1.602*10^-19/(0.6582*10^-15);
        IL(j,k) = (I_pos+I_neg)/sum(P)*Unitconv; %A
        P = P/sum(P);
    end
end

%%
%calculation of the differential conductance
G = diff(IL,1,1)/diff(Vb_vec(1:2));
%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%
% Graphical section
%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Plot of the current as a function of the bias at zero gate voltage
figure(1)
%subplot(2,2,1)
plot(Vb_vec*10^3,IL(:,20)*10^9)
xlabel('V_L-V_R (meV)','Fontsize',18)
ylabel('I (nA)','Fontsize',18)
title('Current (V_g=0 V)','Fontsize',18)

%Plot of the differential conductance as a function of the gate volage
%at 0 bias.
figure(2)
%subplot(2,2,2)
plot(Vg_vec*10^3,G(floor(length(Vb_vec)/2),:)*10^6)
xlabel('V_g (meV)','Fontsize',18)
ylabel('G (\muS)','Fontsize',18)
%ylim([0 2*max(max(G))*10^6])
title('Differential conductance (V=0 V)','Fontsize',18)
hold on


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Plot of the stationary current as a function of the bias voltage and
%of the gate voltage (the so called stability diagram)
figure(3)
%subplot(2,2,3)
surf(Vg_vec*10^3,Vb_vec*10^3,IL)
%note the conversion to meV. For some funny reason
%Matlab wants the input of the surf comand like (X,Y,f(Y,X))!!!!
xlabel('V_g (meV)','Fontsize',18)
ylabel('V_L-V_R (meV)','Fontsize',18)
title('Current','Fontsize',18)
box on
axis tight
shading flat
view([0 90])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Plot of the differential conductance as a function of the bias voltage
%and of the gate voltage (the so called stability diagram)
figure(4)
%subplot(2,2,4)
surf(Vg_vec*10^3,Vb_vec(1:end-1)*10^3,G)
%note the conversion to meV. For some funny reason
%Matlab wants the input of the surf comand like (X,Y,f(Y,X))!!!!
xlabel('V_g (meV)','Fontsize',18)
ylabel('V_L-V_R (meV)','Fontsize',18)
title('Differential conductance','Fontsize',18)
box on
axis tight
shading flat
view([0 90])
colorbar
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%