% Program to generate a cartoon of two 1s atomic orbitals converging to
% produce bonding and antibonding orbitals

clear; close all;

vid = VideoWriter('bondingAntibonding.mp4','MPEG-4');
vid.Quality = 100;
vid.FrameRate = 60;
open(vid);
figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
set(0,'defaultfigurecolor',[1 1 1]);
set(gca,'linewidth',10); 

myBlue = [0.25 0.28 0.5];
grey = [0.5 0.5 0.5]; 

resn = 0.01;
plotHalfWidth = 3; % Relative halfwidth from atom centre that wavefunction is plotted
plottedDataPts = 1 + plotHalfWidth/resn;
rao = 8; % Right atom origin at start
lao = -rao; % Left atom origin at start
vmoveMax = 7.1; 
L = (0:resn:2*rao+plotHalfWidth);
bondLength = 2.5; 
nuclearRadius = 0.2; 
I0 = 3.5; % Height of isolated wavefunction

theta = (0:pi/45:2*pi);
x = nuclearRadius*sin(theta);
y = nuclearRadius*cos(theta);

% Creation of bonding orbital 
% ___________________________

for move = 0:0.025:rao-bondLength/2 % For loop converging two atomic orbitals 
    hold off
    % Decay of atom wavefunction to right
    Int1 = exp(-(L));
    % Decay of atom wavefunction to left
    Int2 = exp(L-plotHalfWidth);
    % Component due to right atom wavefunction bleeding in to right part of
    % left atom wavefunction
    Int3 = exp(L-2*(rao-move));
    % Component due to right atom wavefunction bleeding in to left part of
    % left atom wavefunction
    Int4 = exp(L-2*(rao-move)-plotHalfWidth);
    
    IntLR1 = I0*(Int1 + Int3); % Righthand side of left wavefunction 
    IntLL1 = I0*(Int2 + Int4); % Lefthand side of left wavefunction     
    decayRange1 = L(1:plottedDataPts)+lao+move; 
    m1 = (decayRange1 <= 0);
    plot(decayRange1(m1),IntLR1(m1),'lineWidth',4,'Color',myBlue);
    hold on
    line([lao-plotHalfWidth lao+plotHalfWidth],[0 0],'Color',grey,'linewidth',8) 
    line([rao-plotHalfWidth rao+plotHalfWidth],[0 0],'Color',grey,'linewidth',8)  
    plot(L(1:plottedDataPts)+lao+move-plotHalfWidth,IntLL1(1:plottedDataPts),...
        'lineWidth',4,'Color',myBlue);

    % Repeat mirror image wavefunction on right hand side 
    plot(-decayRange1(m1),IntLR1(m1),'lineWidth',4,'Color',myBlue);
    plot(-L(1:plottedDataPts)-lao-move+plotHalfWidth,IntLL1(1:plottedDataPts),...
        'lineWidth',4,'Color',myBlue);

    presentAtomOrigin = L(1)+lao+move;
    fill(x+presentAtomOrigin,y,myBlue,'linestyle','none')
    fill(x-presentAtomOrigin,y,myBlue,'linestyle','none')
    axis equal 
    axis off 
    xlim([lao-plotHalfWidth rao+plotHalfWidth]);
    ylim([-7.3 11]);
    frame = getframe(gcf);
    writeVideo(vid,frame);
end
pauseMov(vid,60) 
for vmove = 0:vmoveMax/120:vmoveMax % Move resulting molecular orbital down by vmoveMax 
    hold off 
    plot(decayRange1(m1),IntLR1(m1)-vmove,'lineWidth',4,'Color',myBlue); 
    hold on
    line([lao-plotHalfWidth lao+plotHalfWidth],[0 0],'Color',grey,'linewidth',8) 
    line([rao-plotHalfWidth rao+plotHalfWidth],[0 0],'Color',grey,'linewidth',8)  
    plot(L(1:plottedDataPts)+lao+move-plotHalfWidth,IntLL1(1:plottedDataPts)-vmove,...
        'lineWidth',4,'Color',myBlue);
    if (vmove == vmoveMax) 
        plotx1 = decayRange1(m1); % x-values associated with IntLR1 in final molecular orbital 
        plotx2 = L(1:plottedDataPts)+lao+move-plotHalfWidth; % x-values associated with IntLL1 in final molecular orbital
        ploty1 = IntLR1(m1); 
    end 

    % Repeat mirror image wavefunction on right hand side 
    plot(-decayRange1(m1),IntLR1(m1)-vmove,'lineWidth',4,'Color',myBlue);
    plot(-L(1:plottedDataPts)-lao-move+plotHalfWidth,IntLL1(1:plottedDataPts)-vmove,...
        'lineWidth',4,'Color',myBlue);
    fill(x+presentAtomOrigin,y-vmove,myBlue,'linestyle','none')
    fill(x-presentAtomOrigin,y-vmove,myBlue,'linestyle','none')
    axis equal 
    axis off 
    xlim([lao-plotHalfWidth rao+plotHalfWidth]);
    ylim([-7.3 11]);
    frame = getframe(gcf);
    writeVideo(vid,frame);
end 
for fadeIn = 0:0.025:1 % Fade in molecular orbital energy level 
    hold off 
    plot(decayRange1(m1),IntLR1(m1)-vmove,'lineWidth',4,'Color',myBlue); 
    hold on
    line([lao-plotHalfWidth lao+plotHalfWidth],[0 0],'Color',grey,'linewidth',8) 
    line([rao-plotHalfWidth rao+plotHalfWidth],[0 0],'Color',grey,'linewidth',8)  
    plot(L(1:plottedDataPts)+lao+move-plotHalfWidth,IntLL1(1:plottedDataPts)-vmove,...
        'lineWidth',4,'Color',myBlue);

    % Repeat mirror image wavefunction on right hand side 
    plot(-decayRange1(m1),IntLR1(m1)-vmove,'lineWidth',4,'Color',myBlue);
    plot(-L(1:plottedDataPts)-lao-move+plotHalfWidth,IntLL1(1:plottedDataPts)-vmove,...
        'lineWidth',4,'Color',myBlue);
    bond = line([-plotHalfWidth-bondLength/2 plotHalfWidth+bondLength/2],...
        [-vmove -vmove],'Color',grey,'linewidth',8); 
    bond.Color(4) = fadeIn;
    fill(x+presentAtomOrigin,y-vmove,myBlue,'linestyle','none')
    fill(x-presentAtomOrigin,y-vmove,myBlue,'linestyle','none')
    axis equal 
    axis off 
    xlim([lao-plotHalfWidth rao+plotHalfWidth]);
    ylim([-7.3 11]);
    frame = getframe(gcf);
    writeVideo(vid,frame);
end 
pauseMov(vid,60)

% Creation of antibonding orbital 
% _______________________________


for fadeIn = 0:1/90:1 % For loop fading in out-of-phase orbitals 
    hold off
    % Decay of atom wavefunction to right
    Int1 = exp(-(L));
    % Decay of atom wavefunction to left
    Int2 = exp(L-plotHalfWidth);
    % Component due to right atom wavefunction bleeding in to right part of
    % left atom wavefunction
    Int3 = -exp(L-2*(rao));
    % Component due to right atom wavefunction bleeding in to left part of
    % left atom wavefunction
    Int4 = -exp(L-2*(rao)-plotHalfWidth);
    
    IntLR2 = I0*(Int1 + Int3); % Righthand side of left wavefunction 
    IntLL2 = I0*(Int2 + Int4); % Lefthand side of left wavefunction     
    decayRange1 = L(1:plottedDataPts)+lao; 
    m1 = (decayRange1 <= 0);
    plot(decayRange1(m1),IntLR2(m1),'lineWidth',4,'Color',[0.25 0.28 0.5 fadeIn]);
    hold on
    line([lao-plotHalfWidth lao+plotHalfWidth],[0 0],'Color',grey,'linewidth',8) 
    line([rao-plotHalfWidth rao+plotHalfWidth],[0 0],'Color',grey,'linewidth',8)  
    plot(L(1:plottedDataPts)+lao-plotHalfWidth,IntLL2(1:plottedDataPts),...
        'lineWidth',4,'Color',[0.25 0.28 0.5 fadeIn]);

    % Repeat mirror image wavefunction on right hand side 
    plot(-decayRange1(m1),-IntLR2(m1),'lineWidth',4,'Color',[1 0 0 fadeIn]);
    plot(-L(1:plottedDataPts)-lao+plotHalfWidth,-IntLL2(1:plottedDataPts),...
        'lineWidth',4,'Color',[1 0 0 fadeIn]);

    presentAtomOrigin = L(1)+lao;
    core1 = fill(x+presentAtomOrigin,y,myBlue,'linestyle','none'); 
    set(core1,'facealpha',fadeIn);
    core2 = fill(x-presentAtomOrigin,y,myBlue,'linestyle','none'); 
    set(core2,'facealpha',fadeIn);
    
    % Replot lower bonding orbital 
    plot(plotx1,ploty1-vmoveMax,'lineWidth',4,'Color',myBlue); 
    plot(plotx2,IntLL1(1:plottedDataPts)-vmoveMax,...
        'lineWidth',4,'Color',myBlue);

    % Repeat mirror image wavefunction on right hand side 
    plot(-plotx1,ploty1-vmoveMax,'lineWidth',4,'Color',myBlue);
    plot(-plotx2,IntLL1(1:plottedDataPts)-vmoveMax,...
        'lineWidth',4,'Color',myBlue);
    bond = line([-plotHalfWidth-bondLength/2 plotHalfWidth+bondLength/2],...
        [-vmoveMax -vmoveMax],'Color',grey,'linewidth',8); 
    fill(x+bondLength/2,y-vmoveMax,myBlue,'linestyle','none')
    fill(x-bondLength/2,y-vmoveMax,myBlue,'linestyle','none')

    axis equal 
    axis off 
    xlim([lao-plotHalfWidth rao+plotHalfWidth]);
    ylim([-7.3 11]);
    frame = getframe(gcf);
    writeVideo(vid,frame);
end
for move = 0:0.025:rao-bondLength/2 % For loop converging two atomic orbitals 
    hold off
    % Decay of atom wavefunction to right
    Int1 = exp(-(L));
    % Decay of atom wavefunction to left
    Int2 = exp(L-plotHalfWidth);
    % Component due to right atom wavefunction bleeding in to right part of
    % left atom wavefunction
    Int3 = -exp(L-2*(rao-move));
    % Component due to right atom wavefunction bleeding in to left part of
    % left atom wavefunction
    Int4 = -exp(L-2*(rao-move)-plotHalfWidth);
    
    IntLR2 = I0*(Int1 + Int3); % Righthand side of left wavefunction 
    IntLL2 = I0*(Int2 + Int4); % Lefthand side of left wavefunction     
    decayRange1 = L(1:plottedDataPts)+lao+move; 
    m1 = (decayRange1 <= 0);
    plot(decayRange1(m1),IntLR2(m1),'lineWidth',4,'Color',myBlue);
    hold on
    line([lao-plotHalfWidth lao+plotHalfWidth],[0 0],'Color',grey,'linewidth',8) 
    line([rao-plotHalfWidth rao+plotHalfWidth],[0 0],'Color',grey,'linewidth',8)  
    plot(L(1:plottedDataPts)+lao+move-plotHalfWidth,IntLL2(1:plottedDataPts),...
        'lineWidth',4,'Color',myBlue);

    % Repeat mirror image wavefunction on right hand side 
    plot(-decayRange1(m1),-IntLR2(m1),'lineWidth',4,'Color','r');
    plot(-L(1:plottedDataPts)-lao-move+plotHalfWidth,-IntLL2(1:plottedDataPts),...
        'lineWidth',4,'Color','r');

    presentAtomOrigin = L(1)+lao+move;
    fill(x+presentAtomOrigin,y,myBlue,'linestyle','none')
    fill(x-presentAtomOrigin,y,myBlue,'linestyle','none') 
    
    % Replot lower bonding orbital 
    plot(plotx1,ploty1-vmoveMax,'lineWidth',4,'Color',myBlue); 
    plot(plotx2,IntLL1(1:plottedDataPts)-vmoveMax,...
        'lineWidth',4,'Color',myBlue);

    % Repeat mirror image wavefunction on right hand side 
    plot(-plotx1,ploty1-vmoveMax,'lineWidth',4,'Color',myBlue);
    plot(-plotx2,IntLL1(1:plottedDataPts)-vmoveMax,...
        'lineWidth',4,'Color',myBlue);
    bond = line([-plotHalfWidth-bondLength/2 plotHalfWidth+bondLength/2],...
        [-vmoveMax -vmoveMax],'Color',grey,'linewidth',8); 
    fill(x+bondLength/2,y-vmoveMax,myBlue,'linestyle','none')
    fill(x-bondLength/2,y-vmoveMax,myBlue,'linestyle','none')

    axis equal 
    axis off 
    xlim([lao-plotHalfWidth rao+plotHalfWidth]);
    ylim([-7.3 11]);
    frame = getframe(gcf);
    writeVideo(vid,frame);
    if (move==0) 
        pauseMov(vid,240) 
    end     
end
pauseMov(vid,60) 
for vmove = 0:vmoveMax/120:vmoveMax % Move resulting molecular orbital up by vmoveMax 
    hold off 
    plot(decayRange1(m1),IntLR2(m1)+vmove,'lineWidth',4,'Color',myBlue); 
    hold on
    line([lao-plotHalfWidth lao+plotHalfWidth],[0 0],'Color',grey,'linewidth',8) 
    line([rao-plotHalfWidth rao+plotHalfWidth],[0 0],'Color',grey,'linewidth',8)  
    plot(L(1:plottedDataPts)+lao+move-plotHalfWidth,IntLL2(1:plottedDataPts)+vmove,...
        'lineWidth',4,'Color',myBlue);

    % Repeat mirror image wavefunction on right hand side 
    plot(-decayRange1(m1),-IntLR2(m1)+vmove,'lineWidth',4,'Color','r');
    plot(-L(1:plottedDataPts)-lao-move+plotHalfWidth,-IntLL2(1:plottedDataPts)+vmove,...
        'lineWidth',4,'Color','r');
    fill(x+presentAtomOrigin,y+vmove,myBlue,'linestyle','none')
    fill(x-presentAtomOrigin,y+vmove,myBlue,'linestyle','none') 
    
    % Replot lower bonding orbital 
    plot(plotx1,ploty1-vmoveMax,'lineWidth',4,'Color',myBlue); 
    plot(plotx2,IntLL1(1:plottedDataPts)-vmoveMax,...
        'lineWidth',4,'Color',myBlue);

    % Repeat mirror image wavefunction on right hand side 
    plot(-plotx1,ploty1-vmoveMax,'lineWidth',4,'Color',myBlue);
    plot(-plotx2,IntLL1(1:plottedDataPts)-vmoveMax,...
        'lineWidth',4,'Color',myBlue);
    bond = line([-plotHalfWidth-bondLength/2 plotHalfWidth+bondLength/2],...
        [-vmoveMax -vmoveMax],'Color',grey,'linewidth',8); 
    fill(x+bondLength/2,y-vmoveMax,myBlue,'linestyle','none')
    fill(x-bondLength/2,y-vmoveMax,myBlue,'linestyle','none')

    axis equal 
    axis off 
    xlim([lao-plotHalfWidth rao+plotHalfWidth]);
    ylim([-7.3 11]);
    frame = getframe(gcf);
    writeVideo(vid,frame);
end 
for fadeIn = 0:0.025:1 % Fade in molecular orbital energy level 
    hold off 
    plot(decayRange1(m1),IntLR2(m1)+vmove,'lineWidth',4,'Color',myBlue); 
    hold on
    line([lao-plotHalfWidth lao+plotHalfWidth],[0 0],'Color',grey,'linewidth',8) 
    line([rao-plotHalfWidth rao+plotHalfWidth],[0 0],'Color',grey,'linewidth',8)  
    plot(L(1:plottedDataPts)+lao+move-plotHalfWidth,IntLL2(1:plottedDataPts)+vmove,...
        'lineWidth',4,'Color',myBlue);

    % Repeat mirror image wavefunction on right hand side 
    plot(-decayRange1(m1),-IntLR2(m1)+vmove,'lineWidth',4,'Color','r');
    plot(-L(1:plottedDataPts)-lao-move+plotHalfWidth,-IntLL2(1:plottedDataPts)+vmove,...
        'lineWidth',4,'Color','r');
    bond = line([-plotHalfWidth-bondLength/2 plotHalfWidth+bondLength/2],...
        [vmove vmove],'Color',grey,'linewidth',8); 
    bond.Color(4) = fadeIn;
    fill(x+presentAtomOrigin,y+vmove,myBlue,'linestyle','none')
    fill(x-presentAtomOrigin,y+vmove,myBlue,'linestyle','none') 
    
    % Replot lower bonding orbital 
    plot(plotx1,ploty1-vmoveMax,'lineWidth',4,'Color',myBlue); 
    plot(plotx2,IntLL1(1:plottedDataPts)-vmoveMax,...
        'lineWidth',4,'Color',myBlue);

    % Repeat mirror image wavefunction on right hand side 
    plot(-plotx1,ploty1-vmoveMax,'lineWidth',4,'Color',myBlue);
    plot(-plotx2,IntLL1(1:plottedDataPts)-vmoveMax,...
        'lineWidth',4,'Color',myBlue);
    bond = line([-plotHalfWidth-bondLength/2 plotHalfWidth+bondLength/2],...
        [-vmoveMax -vmoveMax],'Color',grey,'linewidth',8); 
    fill(x+bondLength/2,y-vmoveMax,myBlue,'linestyle','none')
    fill(x-bondLength/2,y-vmoveMax,myBlue,'linestyle','none')

    axis equal 
    axis off 
    xlim([lao-plotHalfWidth rao+plotHalfWidth]);
    ylim([-7.3 11]);
    frame = getframe(gcf);
    writeVideo(vid,frame);
end 

%Output the movie as an mpg file
close(vid);

