% Program to plot different hybridized orbitals sp, sp2, sp3, sp + 2
% p-orbitals (as in ethyne), and sp2 + 1 p-orbital (as in ethene).
% The lobes are purely schematic and do not follow true quantum-mechanical
% wavefunctions, as these are anyway approximate for multi-electron atoms

clear; close all

disp('Choose one of the following orbitals: sp1, sp2, sp3,');
prompt = 'sp1pp (sp + 2 p-orbitals), sp2p (sp2 + 1 p-orbital) ';
str1 = input(prompt,'s');
str2 = '-orbital.mp4';
fileName = [str1,str2]; % Name of file depends on type of crystal sample
vid = VideoWriter(fileName,'MPEG-4');
vid.Quality = 100;
vid.FrameRate = 30;
open(vid);

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

orng = [1 0.5 0]; % Orange
lp = [-2 -1 1]; % Lighting angle
view(0,16); % Viewing angle
L = 0.7; % Long lobe
S = 0.5; % Short lobe
F = 1.4; % Fat lobe
T = 1.0; % Thin lobe
scale = 3; % Amount to reduce opposite lobe of hybridized orbitals
opac1 = 0.8; % Opacity of lobe
opac2 = 0.64; % Alternative opacity of lobe
sz = [5 5 3.5 5 5]; % Size of transparent lines to keep fixed xyz limits


% Set up parameters for drawing orbital lobes
% First index defines which lobe (1 to 4); the second index defines the
% five parameters associated with that lobe (phi, theta axis-rotation functions,
% length of lobe, width of lobe, opacity); the third index defines which of the five
% orbitals to plot (sp1, sp2, sp3, sp+pp, sp2+p)
%
% sp1
pars(1,:,1) = [90 0 L T opac2]; % Lobe-1 parameters, sp-orbital
pars(2,:,1) = [-90 0 L T opac2]; % Lobe-2 parameters, sp-orbital
pars(3,:,1) = [0 0 L T 0]; % Lobe-3 parameters, sp-orbital. Doesn't exist, so opacity = 0
pars(4,:,1) = [0 0 L T 0]; % Lobe-4 parameters, sp-orbital. Doesn't exist, so opacity = 0
% sp2
pars(1,:,2) = [90 0 L T opac1]; % Lobe-1 parameters, sp2-orbital
pars(2,:,2) = [90 120 L T opac1]; % Lobe-2 parameters, sp2-orbital
pars(3,:,2) = [90 240 L T opac1]; % Lobe-3 parameters, sp2-orbital
pars(4,:,2) = [0 0 L T 0]; % Lobe-4 parameters, sp2-orbital. Doesn't exist, so opacity = 0
% sp3
pars(1,:,3) = [0 0 S T opac1]; % Lobe-1 parameters, sp3-orbital
pars(2,:,3) = [109.5 0 S T opac1]; % Lobe-2 parameters, sp3-orbital
pars(3,:,3) = [109.5 120 S T opac1]; % Lobe-3 parameters, sp3-orbital
pars(4,:,3) = [109.5 240 S T opac1]; % Lobe-4 parameters, sp3-orbital
% sp1pp
pars(1,:,4) = [90 0 L T opac2]; % Lobe-1 parameters, sp-orbital
pars(2,:,4) = [-90 0 L T opac2]; % Lobe-2 parameters, sp-orbital
pars(3,:,4) = [0 0 S F opac1]; % Lobe-3 parameters, p-orbital
pars(4,:,4) = [90 90 S F opac1]; % Lobe-4 parameters, p-orbital
% sp2p
pars(1,:,5) = [90 0 L T opac1]; % Lobe-1 parameters, sp2-orbital
pars(2,:,5) = [90 120 L T opac1]; % Lobe-2 parameters, sp2-orbital
pars(3,:,5) = [90 240 L T opac1]; % Lobe-3 parameters, sp2-orbital
pars(4,:,5) = [0 0 S F opac1]; % Lobe-4 parameters, p-orbital

% Lobe equations
conecoord1 = 0:pi/320:pi/2; % coordinates of light cone schematic
conecoord2 = 0:pi/150:2*pi; % coordinates of light cone schematic
[conecoord1,conecoord2]=meshgrid(conecoord1,conecoord2);
xcoordcone1 = 4*(1 - sin(conecoord1)).*sin(conecoord1).*cos(conecoord2);
ycoordcone1 = 4*(1 - sin(conecoord1)).*sin(conecoord1).*sin(conecoord2);
zcoordcone1 = 7*(abs(cos(conecoord1))).^3;

if strcmp('sp1',str1)
    orb = 1;
elseif strcmp('sp2',str1)
    orb = 2;
elseif strcmp('sp3',str1)
    orb = 3;
elseif strcmp('sp1pp',str1)
    orb = 4;
elseif strcmp('sp2p',str1)
    orb = 5;
else
    disp('Unknown orbital');
end

%_________________________
% Plotting part of program
for i = 0:2:358 % Rotation around vertical axis
    % Lobes 1 - 4 in blue
    lobe1 = surf(xcoordcone1*pars(1,4,orb),ycoordcone1*pars(1,4,orb),...
        zcoordcone1*pars(1,3,orb),'FaceAlpha',pars(1,5,orb),'FaceColor',...
        'b','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(lobe1,[1 0 0],pars(1,1,orb),[0 0 0]);
    rotate(lobe1,[0 0 1],pars(1,2,orb)+i,[0 0 0]);
    hold on
    lobe2 = surf(xcoordcone1*pars(2,4,orb),ycoordcone1*pars(2,4,orb),...
        zcoordcone1*pars(2,3,orb),'FaceAlpha',pars(2,5,orb),'FaceColor',...
        'b','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(lobe2,[1 0 0],pars(2,1,orb),[0 0 0]);
    rotate(lobe2,[0 0 1],pars(2,2,orb)+i,[0 0 0]);

    lobe3 = surf(xcoordcone1*pars(3,4,orb),ycoordcone1*pars(3,4,orb),...
        zcoordcone1*pars(3,3,orb),'FaceAlpha',pars(3,5,orb),'FaceColor',...
        'b','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(lobe3,[1 0 0],pars(3,1,orb),[0 0 0]);
    rotate(lobe3,[0 0 1],pars(3,2,orb)+i,[0 0 0]);

    lobe4 = surf(xcoordcone1*pars(4,4,orb),ycoordcone1*pars(4,4,orb),...
        zcoordcone1*pars(4,3,orb),'FaceAlpha',pars(4,5,orb),'FaceColor',...
        'b','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(lobe4,[1 0 0],pars(4,1,orb),[0 0 0]);
    rotate(lobe4,[0 0 1],pars(4,2,orb)+i,[0 0 0]);

    %------------------------------------------------------------------------
    % Opposite lobes to lobes 1 - 4 in orange

    lobe5 = surf(xcoordcone1*pars(1,4,orb)/scale,ycoordcone1*pars(1,4,orb)/scale,...
        -zcoordcone1*pars(1,3,orb)/scale,'FaceAlpha',pars(1,5,orb),'FaceColor',...
        orng,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(lobe5,[1 0 0],pars(1,1,orb),[0 0 0]);
    rotate(lobe5,[0 0 1],pars(1,2,orb)+i,[0 0 0]);

    lobe6 = surf(xcoordcone1*pars(2,4,orb)/scale,ycoordcone1*pars(2,4,orb)/scale,...
        -zcoordcone1*pars(2,3,orb)/scale,'FaceAlpha',pars(2,5,orb),'FaceColor',...
        orng,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(lobe6,[1 0 0],pars(2,1,orb),[0 0 0]);
    rotate(lobe6,[0 0 1],pars(2,2,orb)+i,[0 0 0]);

    if strcmp('sp1pp',str1)
        scale2 = 1; % Unhybridized p-orbital
    else
        scale2 = scale; % Hybridized
    end
    lobe7 = surf(xcoordcone1*pars(3,4,orb)/scale2,ycoordcone1*pars(3,4,orb)/scale2,...
        -zcoordcone1*pars(3,3,orb)/scale2,'FaceAlpha',pars(3,5,orb),'FaceColor',...
        orng,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(lobe7,[1 0 0],pars(3,1,orb),[0 0 0]);
    rotate(lobe7,[0 0 1],pars(3,2,orb)+i,[0 0 0]);

    if strcmp('sp1pp',str1) || strcmp('sp2p',str1)
        scale2 = 1; % Unhybridized p-orbital
    else
        scale2 = scale; % Hybridized
    end
    lobe8 = surf(xcoordcone1*pars(4,4,orb)/scale2,ycoordcone1*pars(4,4,orb)/scale2,...
        -zcoordcone1*pars(4,3,orb)/scale2,'FaceAlpha',pars(4,5,orb),'FaceColor',...
        orng,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
    rotate(lobe8,[1 0 0],pars(4,1,orb),[0 0 0]);
    rotate(lobe8,[0 0 1],pars(4,2,orb)+i,[0 0 0]);

    % Dummy horizontal and vertical lines to stop image resizing
    % Four-element color definition [rgb alpha] only works with MATLAB 2021
    % onwards
    plot3([sz(orb) -sz(orb)],[0 0],[-sz(orb) -sz(orb)],...
        'color',[1 1 1 0], 'LineWidth', 0.001);
    plot3([0 0],[sz(orb) -sz(orb)],[-sz(orb) -sz(orb)],...
        'color',[1 1 1 0], 'LineWidth', 0.001);
    plot3([-sz(orb) -sz(orb)],[0 0],[sz(orb) -sz(orb)],...
        'color',[1 1 1 0], 'LineWidth', 0.001);

    axis off
    axis equal
    axis tight
    light('Position',lp,'Style','infinite');

    frame = getframe(gcf);
    writeVideo(vid,frame);

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