% Program to plot the dispersion curve and density of states of a 1-D array
% of atoms according to the Kronig-Penney model, a simplified model in
% solid-state physics for electrons in a one-dimensional pulselike
% periodic potential. See e.g.,
% https://en.wikipedia.org/wiki/Particle_in_a_one-dimensional_lattice for
% the physics background.

close all
clear 

figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
set(0,'defaultfigurecolor',[1 1 1]);

% Locations and sizes of subplots
pos1 = [0.05 0.25 0.28 0.4];
pos2 = [0.38 0.25 0.28 0.4];
pos3 = [0.71 0.25 0.28 0.4];

nFrames = 300; % Number of frames in the movie 
potLength = 8; % Length of potential plot in x-direction in Angstroms 

% Constants
pi = 3.14159265;
hbar = 1.05457E-34;
ec = 1.6022E-19;

% Functions for hyperbolic cosine and sine
cosh = @(x) (exp(x) + exp(-x)) / 2;
sinh = @(x) (exp(x) - exp(-x)) / 2;

prompt = ['Which parameter do you want to scan? \n' ...
    'V = potential difference between maximum and minimum in periodic perturbation \n' ...
    'P = Periodicity \n' ...
    'D = duty cycle (peak width/periodicity): '];
varPar = input(prompt, 's');

prompt = 'Lower energy limit for plot [eV]? Default = -1 eV: ';
E1 = input(prompt);
if isempty(E1)
    E1 = -1;
end

prompt = 'Upper energy limit for plot [eV]? Default = 50 eV: ';
E2 = input(prompt);
if isempty(E2)
    E2 = 50;
end

prompt = 'Lower potential of periodic square wave [V]? Default = 0 V: ';
V1 = input(prompt);
if isempty(V1)
    V1 = 0;
end

% Variable value of pulse potential height
if (varPar == 'v' || varPar == 'V') % Potential is the variable
    prompt = 'Starting upper potential of periodic square wave [V]? : ';
    V21 = input(prompt);
    prompt = 'Ending upper potential of periodic square wave [V]? : ';
    V22 = input(prompt);
    Vmax = V22; 
    loopVals = linspace(V21,V22,nFrames);
    vid = VideoWriter('kronigPenney_potential.mp4','MPEG-4');
else % Not the variable
    prompt = 'Upper potential of periodic square wave [V]? Default = 10 V: ';
    V2 = input(prompt);
    Vmax = V2; 
    if isempty(V2)
        V2 = 10;
        Vmax = V2; 
    end
end

% Variable value of pulse potential periodicity
if (varPar == 'p' || varPar == 'P') % Periodicity is the variable
    prompt = 'Starting periodicity [Å]? : ';
    P1 = input(prompt)*1e-10;
    prompt = 'Ending periodicity [Å]? : ';
    P2 = input(prompt)*1e-10;
    loopVals = linspace(P1,P2,nFrames);
    vid = VideoWriter('kronigPenney_periodicity.mp4','MPEG-4');
else % Not the variable
    prompt = 'Periodicity of 1-D potential [Å]? Default = 1.6 Å: ';
    a = input(prompt);
    if isempty(a)
        a = 1.6e-10;
    else
        a = a*1e-10; % Convert to m
    end
end

% Variable value of pulse potential duty cycle
if (varPar == 'd' || varPar == 'D') % Duty cycle is the variable
    prompt = 'Starting pulse width [Å]? : ';
    c1 = input(prompt)*1e-10;
    prompt = 'Ending pulse width [Å]? : ';
    c2 = input(prompt)*1e-10;
    b1 = a - c1; % Starting width of lower potential in m
    b2 = a - c2; % Ending width of lower potential in m
    loopVals = linspace(b1,b2,nFrames);
    vid = VideoWriter('kronigPenney_dutyCycle.mp4','MPEG-4');
else % Not the variable
    prompt = 'Width of peak voltage of 1-D potential [Å]? Default = 0.6 Å: ';
    c = input(prompt);
    if isempty(c)
        c = 0.6e-10;
    else
        c = c*1e-10; % Convert to m
    end
end 

m = 9.11e-31; % Mass of electron in kg
N = 5000; % Number of energy points to plot

vid.Quality = 100;
vid.FrameRate = 30;
open(vid); 

% Calculate dispersion and DoS curves 
for ii = loopVals
    % Initialize arrays
    dispersion = zeros(N, 2);
    density_of_states = zeros(N, 2);
    if (varPar == 'd' || varPar == 'D') % Duty cycle is the variable 
        b = ii; % Change b-value for each loop 
        V2 = Vmax; 
    end 
    if (varPar == 'v' || varPar == 'V') % Peak potential is the variable 
        V2 = ii; % Change upper-voltage-value for each loop 
        b = a - c; 
    end 
    if (varPar == 'p' || varPar == 'P') % Periodicity is the variable 
        a = ii; % Change periodicity for each loop 
        b = ii - c; % Width of lower potential in m
        V2 = Vmax; 
    end 
    potLoop = ceil(potLength*1e-10/a);
    potential = zeros(4*potLoop,2);
    for jj = 0:potLoop-1
        potential(1+jj*4,1) = 0 + a*jj*1e9;
        potential(1+jj*4,2) = V1;
        potential(2+jj*4,1) = b*1e9 + a*jj*1e9;
        potential(2+jj*4,2) = V1;
        potential(3+jj*4,1) = b*1e9 + a*jj*1e9;
        potential(3+jj*4,2) = V2;
        potential(4+jj*4,1) = a*1e9 + a*jj*1e9;
        potential(4+jj*4,2) = V2;
    end

    % Convert potentials to Joules 
    V1 = V1 * ec;
    V2 = V2 * ec;

    maxDOE = 0.0; 

    % Loop through energy levels
    for i = 1:N % i-loop 
        E = E1 + (E2 - E1) * (i - 1) / (N - 1);
        EJ = E * ec; % energy in Joules
        k1 = sqrt(2 * m * abs(EJ - V1)) / hbar;
        k2 = sqrt(2 * m * abs(EJ - V2)) / hbar;
        dk1dE = m / (hbar^2 * k1);
        dk2dE = m / (hbar^2 * k2);

        if EJ < V1
            alpha = 2 * cosh(k2 * (a - b)) * cosh(k1 * b) + (k2 / k1 + k1 / k2) * sinh(k2 * (a - b)) * sinh(k1 * b);
            dispersion(i, :) = [0, E];
            density_of_states(i, :) = [0, E];
            if (max(density_of_states(i,2)) > maxDOE)
                maxDOE = max(density_of_states(i,2));
            end 
        else
            if EJ < V2
                dk2dE = -m / (hbar^2 * k2);
                alpha = 2 * cosh(k2 * (a - b)) * cos(k1 * b) + (k2 / k1 - k1 / k2) * sinh(k2 * (a - b)) * sin(k1 * b);
                dalphadk1 = -2 * b * cosh(k2 * (a - b)) * sin(k1 * b) + (-k2 / (k1^2) - 1 / k2) * sinh(k2 * (a - b)) * sin(k1 * b) + b * (k2 / k1 - k1 / k2) * sinh(k2 * (a - b)) * cos(k1 * b);
                dalphadk2 = 2 * (a - b) * sinh(k2 * (a - b)) * cos(k1 * b) + (1 / k1 + k1 / (k2^2)) * sinh(k2 * (a - b)) * sin(k1 * b) + (a - b) * (k2 / k1 - k1 / k2) * cosh(k2 * (a - b)) * sin(k1 * b);
            else
                alpha = 2 * cos(k2 * (a - b)) * cos(k1 * b) - (k2 / k1 + k1 / k2) * sin(k2 * (a - b)) * sin(k1 * b);
                dalphadk1 = -2 * b * cos(k2 * (a - b)) * sin(k1 * b) - (-k2 / (k1^2) + 1 / k2) * sin(k2 * (a - b)) * sin(k1 * b) - b * (k2 / k1 + k1 / k2) * sin(k2 * (a - b)) * cos(k1 * b);
                dalphadk2 = -2 * (a - b) * sin(k2 * (a - b)) * cos(k1 * b) - (1 / k1 - k1 / (k2^2)) * sin(k2 * (a - b)) * sin(k1 * b) - (a - b) * (k2 / k1 + k1 / k2) * cos(k2 * (a - b)) * sin(k1 * b);
            end

            dkda = 1 / (a * sqrt(4 - alpha^2));

            if alpha > 0
                if alpha^2 < 4
                    dispersion(i, :) = [atan(sqrt(4 - alpha^2) / alpha), E];
                    density_of_states(i, :) = [abs(ec * 1E-9 * (2 / pi) * dkda * (dalphadk1 * dk1dE + dalphadk2 * dk2dE)), E];
                else
                    density_of_states(i, :) = [0, E];
                    dispersion(i, :) = [0, E];
                end
            else
                if alpha^2 < 4
                    dispersion(i, :) = [pi + atan(sqrt(4 - alpha^2) / alpha), E];
                    density_of_states(i, :) = [abs(ec * 1E-9 * (2 / pi) * dkda * (dalphadk1 * dk1dE + dalphadk2 * dk2dE)), E];
                else
                    density_of_states(i, :) = [0, E];
                    dispersion(i, :) = [pi, E];
                end
            end
            if (max(density_of_states(i,2)) > maxDOE)
                maxDOE = max(density_of_states(i,2));
            end 
        end 
    end % end of i-loop 

    % Plotting
    hold off
    subplot('position',pos1);
    newplot
    plot(potential(:, 1), potential(:, 2), 'r','LineWidth', 2.0);
    xlabel('Position [nm]');
    ylabel('Potential [eV]');
    xlim([0 potLength/10]);
    ylim([V1,Vmax*1.05]);
    set(gca,'fontsize',18);
    set(gca,'TickLength',[0.02, 1]);
    set(gca,'linewidth',2);
    box('off')
    if (varPar == 'd' || varPar == 'D') % Duty cycle is the variable 
        str1 = 'Potential = ';
        str2 = num2str(V2/ec,4); 
        str3 = ' eV';
        strTot = [str1 str2 str3]; 
        text(0.0, 1.1, strTot, 'FontSize',18,'units','normalized');
        str1 = 'Periodicity = ';
        str2 = num2str(a/1e-10,4); 
        str3 = ' Å';
        strTot = [str1 str2 str3]; 
        text(0.35, 1.1, strTot, 'FontSize',18,'units','normalized');
        str1 = 'Duty cycle = ';
        str2 = num2str((a-ii)/a,4); 
        strTot = [str1 str2]; 
        text(0.7, 1.1, strTot, 'FontSize',18,'units','normalized');
    end 
    if (varPar == 'v' || varPar == 'V') % Peak potential is the variable 
        str1 = 'Potential = ';
        str2 = num2str(ii,4); 
        str3 = ' eV';
        strTot = [str1 str2 str3]; 
        text(0.0, 1.1, strTot, 'FontSize',18,'units','normalized');
        str1 = 'Periodicity = ';
        str2 = num2str(a/1e-10,4); 
        str3 = ' Å';
        strTot = [str1 str2 str3]; 
        text(0.37, 1.1, strTot, 'FontSize',18,'units','normalized');
        str1 = 'Duty cycle = ';
        str2 = num2str((a-b)/a,4); 
        strTot = [str1 str2]; 
        text(0.7, 1.1, strTot, 'FontSize',18,'units','normalized');
    end 
    if (varPar == 'p' || varPar == 'P') % Periodicity is the variable 
        str1 = 'Potential = ';
        str2 = num2str(V2/ec,4); 
        str3 = ' eV';
        strTot = [str1 str2 str3]; 
        text(0.0, 1.1, strTot, 'FontSize',18,'units','normalized');
        str1 = 'Periodicity = ';
        str2 = num2str(ii/1e-10,4); 
        str3 = ' Å';
        strTot = [str1 str2 str3]; 
        text(0.32, 1.1, strTot, 'FontSize',18,'units','normalized');
        str1 = 'Duty cycle = ';
        str2 = num2str((a-b)/a,4); 
        strTot = [str1 str2]; 
        text(0.7, 1.1, strTot, 'FontSize',18,'units','normalized');
    end 


    hold off
    subplot('position',pos2);
    newplot
    plot(1.001*dispersion(:, 1)-0.001, dispersion(:, 2), 'g','LineWidth', 2.0);
    xlabel('ka');
    ylabel('E [eV]');
    xticks([0 pi/4 pi/2 3*pi/4 pi])
    xticklabels({'0','\pi/4','\pi/2','3\pi/4','\pi'})
    xlim([0,pi]);
    ylim([E1,E2]);
    set(gca,'fontsize',18);
    set(gca,'TickLength',[0.02, 1]);
    set(gca,'linewidth',2);

    hold off
    subplot('position',pos3);
    newplot
    lastPt = [0 50];
    density_of_states = [density_of_states; lastPt];
    plot(density_of_states(:, 1), density_of_states(:, 2), 'b','LineWidth', 2.0);
    hold on 
    fill(density_of_states(:, 1), density_of_states(:, 2), 'b','FaceAlpha',0.25);
    xlabel('Density of States [eV^{-1}nm^{-1}]');
    ylabel('E [eV]'); 
    xlim([0,2.5]);
    ylim([E1,E2]);
    set(gca,'fontsize',18);
    set(gca,'TickLength',[0.02, 1]);
    set(gca,'linewidth',2);


    frame = getframe(gcf);
    writeVideo(vid,frame);
end % loopVals loop
close(vid) 