% Program to plot the interference between two point oscillators as a
% function of their separation. The time-averaged far-field intensity is
% plotted on the right.

clear all; close all;

vid = VideoWriter('interference.mp4','MPEG-4');
vid.Quality = 70;
vid.FrameRate = 30;
vid.VideoCompressionMethod
open(vid);
figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
set(0,'defaultfigurecolor',[1 1 1]);
set(gca,'linewidth',10);
myBlue = [0.4 0.44 0.73];

axis square
axis equal
axis off
hold on

lp = [0.2 0.5 0.8]; % Light projection angle
ec = [0.35 0.39 0.64]; % Phil blue
fc = [1.0 0.83 0]; % Gold
alphaVal = 0.82;

set(gca,'View',[0,90]);
%set(gca, 'Projection','perspective');
%light('Position',lp,'Style','infinite');

per = 2*pi; % periodicity

% Limits of real-space volume in Angstroms
xlim([-4*per 35*per]);
ylim([-11*per 11*per]);
zlim([-4.1 4.1]);

% surface mesh grid for propagating circular wave
leftLim = 4;
rightLim = 35;
topBottomLim = 11;
x = -leftLim*per:per/25:rightLim*per;
y = -topBottomLim*per:per/25:topBottomLim*per;
[X,Y] = meshgrid(x,y);
yrange = -topBottomLim:0.01:topBottomLim;
theta = atan(yrange/rightLim);

% electron
[u,v,w] = sphere; % Create sphere coordinates
er = 1.25; % radius of electron
pos1 = [0.0 0.0 0.8 1];
pos2 = [0.82 0.113 0.17 0.774];

% single electron oscillating up and down in centre
for i = 0:2*pi/30:10*pi-2*pi/60
    subplot('Position',pos1); % Main plot
    newplot
    clear surf1
    set(gca,'View',[0,90]);
    set(gca, 'Projection','perspective');
    %light('Position',lp,'Style','infinite');
    axis square
    axis equal
    axis off
    
    xlim([-4*per 35*per]);
    ylim([-11*per 11*per]);
    zlim([-4.1 4.1]);
    hold on
    light('Position',lp,'Style','infinite');
    
    F = (0.5*cos((X.^2+(Y).^2).^0.5-i)+0.5*cos((X.^2+(Y).^2).^0.5-i)).^2;
    elec1 = surf(u*er,v*er,w*er+2.5*cos(i),'FaceAlpha',1.0, ...
        'FaceColor',ec,'LineStyle','none');
    surf1 = surf(X,Y,F,'FaceColor',fc,'EdgeColor','none','FaceAlpha',alphaVal);
    drawnow;
    hold off
    subplot('Position',pos2); % intensity plot
    ith = (cos(pi*0*sin(theta)).*sinc(pi*sin(theta))).^2;
    plot(ith,yrange,'Color',myBlue,'Linewidth',4);
    axis tight
    axis off
    frame = getframe(gcf);
    writeVideo(vid,frame);
end

% two electrons moving apart from one another
for i = 0:2*pi/30:16*pi-2*pi/60
    subplot('Position',pos1); % Main plot
    newplot
    clear surf1
    set(gca,'View',[0,90]);
    set(gca, 'Projection','perspective');
    %light('Position',lp,'Style','infinite');
    axis square
    axis equal
    axis off
    
    xlim([-4*per 35*per]);
    ylim([-11*per 11*per]);
    zlim([-4.1 4.1]);
    hold on
    light('Position',lp,'Style','infinite');
    
    F = (0.5*cos((X.^2+(Y-i/2).^2).^0.5-i)+0.5*cos((X.^2+(Y+i/2).^2).^0.5-i)).^2;
    elec1 = surf(u*er,v*er+i/2,w*er+2.5*cos(i),'FaceAlpha',1.0, ...
        'FaceColor',ec,'LineStyle','none');
    elec2 = surf(u*er,v*er-i/2,w*er+2.5*cos(i),'FaceAlpha',1.0, ...
        'FaceColor',ec,'LineStyle','none');
    surf1 = surf(X,Y,F,'FaceColor',fc,'EdgeColor','none','FaceAlpha',alphaVal);
    drawnow;
    hold off
    subplot('Position',pos2); % intensity plot
    ith = (cos(pi*(i/(2*pi))*sin(theta)).*sinc(pi*sin(theta))).^2;
    plot(ith,yrange,'Color',myBlue,'Linewidth',4);
    axis tight
    axis off
    frame = getframe(gcf);
    writeVideo(vid,frame);
end

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

