% Thermo-electric current

% Physical constants
hbar = 1.054e-34;
q = 1.602e-19;

% Parameters (eV)
kBT_1 = 0.025;
kBT_2 = 0.026;
mu = 0;
gamma_1 = 0.005;
gamma_2 = gamma_1;
gamma_sum = gamma_1 + gamma_2;

% Channel energy levels, varying between -0.25eV and 0.25eV
epsilon = linspace(-0.25, 0.25, 101);
depsilon = epsilon(2) - epsilon(1);

% Energy grid
E = linspace(-1, 1, 501);
dE = E(2) - E(1);

% Contact fermi functions
f_1 = 1 ./ (1 + exp((E - mu)./kBT_1));
f_2 = 1 ./ (1 + exp((E - mu)./kBT_2));

% Iterate through channel energy levels
for n = 1:length(epsilon)
    % Compute energy level density functions - integral normalized to unity
    D = (gamma_sum./(2*pi))./((E-epsilon(n)).^2+((gamma_sum./2).^2));
    %D = (gamma./(2*pi))./((E-epsilon(n)).ˆ2+((gamma_sum./2).ˆ2));
    D = D./(dE*sum(D));
    
    % Compute number of channel electrons
    
    N(n) = dE*sum( ((gamma_1./gamma_sum).*f_1 + (gamma_2./gamma_sum).*f_2).*D );

    % Compute the current in Amps; factor of q to resolve units

    I(n) = q*(q/hbar)*dE*sum((f_1 - f_2).*D.*gamma_1.*gamma_2./gamma_sum);

    %plot f_1 - f_2 and D/2

    if (abs(epsilon(n) + 0.05) <= depsilon / 2) & (epsilon(n) <= 0)
        figure(3);
        h = plot(f_1-f_2, E, 'x', D/2500, E, 'k-');
        set(gca, 'Fontsize', [18]);
        axis([-0.01 0.02 -1 1]);
        xlabel('f1(E) - f2(E), D(E)/2500');
        ylabel('ENERGY  [eV]');
        legend('f1-f2', 'D(E)/2500');
        title('CHANNEL LEVEL = -0.05 eV');
    elseif (abs(epsilon(n)) <= depsilon / 2) & (epsilon(n) <= 0)
        figure(4);
        h = plot(f_1-f_2, E, 'x', D/2500, E, 'k-');
        set(gca, 'Fontsize', [18]);
        axis([-0.01 0.02 -1 1]);
        xlabel('f1(E) - f2(E), D(E)/2500');
        ylabel('ENERGY  [eV]');
        legend('f1-f2', 'D(E)/2500');
        title('CHANNEL LEVEL = 0 eV');
    end

end


% Final plots
figure(1);