% Cartoon of the spin-orbit interaction

clear
close all
vid = VideoWriter('SpinOrbitCoupling.mp4','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]);
set(gca,'linewidth',7); 

% Plotting areas of two components of video 
pos1 = [0 -0.1 0.7 1.28]; % Cartoon electron orbiting nucleus 
pos2 = [0.7 0.25 0.25 0.5]; % Energy-level scheme 

lp = [7 5 8]; % Angle of lighting
myNickel = [0.44 0.46 0.44];
myBlue = [0.4 0.44 0.73]; 

sRotAlpRed(1,:) = 0:2:3*360-2; % Rotation 2-degree steps over three rotations 
alp1 = zeros(1,90); % No central magnet for first half rotation (90 frames)
alp2 = 0:0.05:1; % Central magnet fade in over 21 frames 
alp3 = ones(1,3*180-90-21); % Central magnet 100% visible
sRotAlpRed(2,:) = [alp1 alp2 alp3]; % Second row determines behaviour of transparency of central magnet 
alp1 = zeros(1,180); % No spin magnet for first rotation (180 frames)
alp2 = 0:0.05:1; % Spin magnet fade in over 21 frames 
alp3 = ones(1,3*180-180-21); % Spin magnet 100% visible
sRotAlpRed(3,:) = [alp1 alp2 alp3]; % Third row determines behaviour of transparency of spin magnet 
Red1 = zeros(1,360); % No magnets for first rotation then in second rotation the colour is flipped from red to black
Red2 = 0:0.05:1; % s-magnet colour changes from black to red over 21 frames
Red3 = ones(1,3*180-360-21); 
sRotAlpRed(4,:) = [Red1 Red2 Red3]; % Fourth row determines colour flip for spin magnet 

[a,b,c] = sphere(50); % Electron or nucleon

rA = 0:0.125:0.125; % For cone arrowhead
[X,Y,Z] = cylinder(rA,50); % Arrow head
arrowLength = 0.4; % Length of arrow head

r2 = 0.22; % Radius of electron 
fac = 1.4; % Field-line ellipse b/a ratio 
gr = (1 + sqrt(5))/2; % Golden ratio 
nuclCol = [0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1]; % Red or black nucleons 
magFalp = 0.5; % Transparency of magnetic field lines 

R=5; % outer radius of torus representinig electron orbit 

for elRot = sRotAlpRed(1,:) % Three rotations 
    index = (elRot+2)/2; 
    subplot('position',pos1); % Cartoon of electron orbiting nucleus 
    newplot
    nucIndex = 0;
    hold off
    % Create filled dodecahedron defined by spheres at the vertices for
    % nuclear core 
    r = 0.16;
    for i = (-1:2:1)*r*fac
        for j = (-1:2:1)*r*fac
            for k = (-1:2:1)*r*fac
                nucIndex = nucIndex+1;
                if (nuclCol(nucIndex)>=0.5)
                    borr = 'r';
                else
                    borr = 'k';
                end
                surf(r*a+i,r*b+j,r*c+k,...
                    'FaceColor',borr,'LineStyle','none','FaceAlpha',1-sRotAlpRed(2,index),...
                    'FaceLighting','flat','DiffuseStrength',1);
                if (nuclCol(nucIndex)>=0.5)
                    borr = 'r';
                else
                    borr = 'k';
                end
                surf(r*a+i/fac,r*b+j/fac,r*c+k/fac,...
                    'FaceColor','r','LineStyle','none','FaceAlpha',(1-sRotAlpRed(2,index)),...
                    'FaceLighting','flat','DiffuseStrength',1);
                hold on
            end
        end
    end
    for i = (-gr:2*gr:gr)*r*fac
        for j = (-1/gr:2/gr:1/gr)*r*fac
            nucIndex = nucIndex+1;
            if (nuclCol(nucIndex)>=0.5)
                borr = 'r';
            else
                borr = 'k';
            end
            surf(r*a,r*b+i,r*c+j,...
                'FaceColor',borr,'LineStyle','none','FaceAlpha',(1-sRotAlpRed(2,index)),...
                'FaceLighting','flat','DiffuseStrength',1);
            nucIndex = nucIndex+1;
            if (nuclCol(nucIndex)>=0.5)
                borr = 'r';
            else
                borr = 'k';
            end
            surf(r*a+j,r*b,r*c+i,...
                'FaceColor',borr,'LineStyle','none','FaceAlpha',(1-sRotAlpRed(2,index)),...
                'FaceLighting','flat','DiffuseStrength',1);
            nucIndex = nucIndex+1;
            if (nuclCol(nucIndex)>=0.5)
                borr = 'r';
            else
                borr = 'k';
            end
            surf(r*a+i,r*b+j,r*c,...
                'FaceColor',borr,'LineStyle','none','FaceAlpha',(1-sRotAlpRed(2,index)),...
                'FaceLighting','flat','DiffuseStrength',1);
            hold on
        end
    end
    
    % Make central magnet representing angular momentum l 
    magWidth = 0.7; 
    [V,F] = platonic_solid(2,magWidth); % Octahedron, size = 1
    V(:,3) = sqrt(3)*V(:,3)-magWidth;
    patch('Faces',F,'Vertices',V,'FaceColor','r','FaceAlpha',sRotAlpRed(2,index), ...
        'EdgeColor','none','FaceLighting','flat','DiffuseStrength',1,'AmbientStrength',1,'SpecularStrength',0);
    [V,F] = platonic_solid(2,magWidth); % Octahedron, size = 1
    V(:,3) = sqrt(3)*V(:,3)+magWidth;
    patch('Faces',F,'Vertices',V,'FaceColor','k','FaceAlpha',sRotAlpRed(2,index), ...
        'EdgeColor','none','FaceLighting','flat','DiffuseStrength',1,'AmbientStrength',1,'SpecularStrength',0);
    material shiny

    % Make electron magnet representing angular momentum s
    magWidth = 0.4; 
    [V,F] = platonic_solid(2,magWidth); % Octahedron, size = 1
    V(:,1) = V(:,1)-R;
    V(:,3) = sqrt(3)*V(:,3) + magWidth;
    colNow = [1-sRotAlpRed(4,index) 0 0]; 
    ps = patch('Faces',F,'Vertices',V,'FaceColor',colNow,'FaceAlpha',sRotAlpRed(3,index), ...
        'EdgeColor','none','FaceLighting','flat','DiffuseStrength',1,'AmbientStrength',1,'SpecularStrength',0);
    direction = [0 0 1];
    rotate(ps,direction,-elRot,[-R 0 0])
    rotate(ps,direction,elRot,[0 0 0])
    [V,F] = platonic_solid(2,magWidth); % Octahedron, size = 1
    V(:,1) = V(:,1)-R;
    V(:,3) = sqrt(3)*V(:,3)-magWidth;
    colNow = [sRotAlpRed(4,index) 0 0]; 
    ps = patch('Faces',F,'Vertices',V,'FaceColor',colNow,'FaceAlpha',sRotAlpRed(3,index), ...
        'EdgeColor','none','FaceLighting','flat','DiffuseStrength',1,'AmbientStrength',1,'SpecularStrength',0);
    direction = [0 0 1];
    rotate(ps,direction,-elRot,[-R 0 0])
    rotate(ps,direction,elRot,[0 0 0])
    material shiny

    
    % Create a very thin torus to describe the electron orbit
    r=0.032; % inner tube radius
    th=linspace(0,2*pi,90); % e.g. 90 partitions along perimeter of the tube
    phi=linspace(0,2*pi,180); % e.g. 180 partitions along azimuth of torus
    % we convert our vectors phi and th to [n x n] matrices with meshgrid command:
    [Phi,Th]=meshgrid(phi,th);
    % now we generate n x n matrices for x,y,z according to eqn of torus
    x=(R+r.*cos(Th)).*cos(Phi);
    y=(R+r.*cos(Th)).*sin(Phi);
    z=r.*sin(Th);
    surf(x,y,z,'FaceColor','k','LineStyle',...
        'none','FaceLighting','flat','DiffuseStrength',1); % plot surface
    
    % Place electron on its orbit, including a spin vector arrow 
    elec = surf(r2*a-R,r2*b,r2*c,...
        'FaceColor',myBlue,'LineStyle',...
        'none','FaceLighting','flat','DiffuseStrength',1);
    rotate(elec,[0 0 1],elRot,[0,0,0]); 
    spinArrow = mArrow3([-R 0 -0.64],[-R 0 0.865],'color',myBlue,'stemWidth',0.04,...
        'tipWidth',0.08); 
    rotate(spinArrow,[0 0 1],elRot,[0,0,0]); 

    % Create a set of very thin torii to describe the magnetic-field lines
    for ii = 0:360/7:360
        for jj = 1:4
            R2=3; % outer radius of torus
            r=0.025; % inner tube radius
            th=linspace(0,2*pi,90); % e.g. 90 partitions along perimeter of the tube
            phi=linspace(0,2*pi,180); % e.g. 180 partitions along azimuth of torus
            % we convert our vectors phi and th to [n x n] matrices with meshgrid command:
            [Phi,Th]=meshgrid(phi,th);
            % now we generate n x n matrices for x,y,z according to eqn of torus
            x=(R2*1.4*jj/2+r.*cos(Th)).*cos(Phi);
            y=(R2*jj/2+r.*cos(Th)).*sin(Phi);
            z=r.*sin(Th);
            % Magnetic flux lines 
            magFL = surf(x,y-R+0.73-jj/1.1-jj^2/25,z,'FaceColor',myBlue,'LineStyle',...
                'none','FaceAlpha',magFalp,'FaceLighting','flat','DiffuseStrength',1); % plot surface
            rotate(magFL,[0 1 0],90,[0,0,0]);
            rotate(magFL,[0 0 1],ii,[0,0,0]);
            % Arrow head showing direction of magnetic flux lines 
            arrowHead = surf(X,Y-R+0.73+(cosd(elRot)*R2/2-1/1.1-jj/25)*jj,-sind(elRot)*R2*1.4*jj/2+Z*arrowLength-arrowLength/2,...
                'FaceColor',myBlue,'LineStyle',...
                'none','FaceLighting','flat','FaceAlpha',magFalp,'DiffuseStrength',1);
            rotate(arrowHead,[1 0 0],-elRot,[0,-R+0.73+(cosd(elRot)*R2/2-1/1.1-jj/25)*jj,-sind(elRot)*R2*1.4*jj/2]);            
            rotate(arrowHead,[0 0 1],ii,[0,0,0]);            
            arrowHeadDum = surf(X,Y-R+0.73+(cosd(180)*R2/2-1/1.1-jj/25)*jj,-sind(180)*R2*1.4*jj/2+Z*arrowLength-arrowLength/2,...
                'FaceColor',myBlue,'LineStyle',...
                'none','FaceLighting','flat','FaceAlpha',0,'DiffuseStrength',1);
            rotate(arrowHeadDum,[1 0 0],-180,[0,-R+0.73+(cosd(180)*R2/2-1/1.1-jj/25)*jj,-sind(180)*R2*1.4*jj/2]);            
            rotate(arrowHeadDum,[0 0 1],ii,[0,0,0]);            
        end
    end
    axis equal
    axis off
    axis tight
    light('Position',lp,'Style','infinite');
    set(gca, 'Projection','perspective');
    set(gca,'View',[139,50]); 
    
    
%__________________________________________________________________________ 
% Spin-orbit split levels j = l - s and j = l + s 
    subplot('position',pos2);
    newplot 
    hold off
    line([0 0.4],[0 0],'Color',[0.5 0.5 0.5],'LineWidth',5); % Unperturbed energy level
    hold on
    text(0.2,0.05,'l','FontSize',20,'HorizontalAlignment', 'center', ...
            'VerticalAlignment', 'bottom','color','k');
        
    % j = l - s energy level 
    rect = findall(gcf,'Type', 'Rectangle'); 
    delete(rect); 
    if (elRot >= 360)
        jlms = plot([0.6 1.0], [-1 -1], 'color',[0.5 0.5 0.5],'Linewidth', 5); 
        connect1 = plot([0.4 0.6], [0 -1], 'color',[0.75 0.75 0.75],'Linewidth', 2.5); 
    end
    if (elRot >= 360) && (elRot < 400)
        jlms.Color(4) = 1 - (400-elRot)/40; % Semitransparent line
        connect1.Color(4) = 1 - (400-elRot)/40; % Semitransparent line
        textCol = [1 1 1]*(400-elRot)/40;
        text(0.8,-0.95,'j = l - s','FontSize',20,'HorizontalAlignment', 'center', ...
            'VerticalAlignment', 'bottom','color',textCol);
        rectangle('Position',[0.775,-0.55,0.05,0.2],'FaceColor',[0 0 0 1-(400-elRot)/40],'EdgeColor','none')
        rectangle('Position',[0.775,-0.75,0.05,0.2],'FaceColor',[1 0 0 1-(400-elRot)/40],'EdgeColor','none')
        rectangle('Position',[0.85,-0.55,0.025,0.1],'FaceColor',[1 0 0 1-(400-elRot)/40],'EdgeColor','none')
        rectangle('Position',[0.85,-0.65,0.025,0.1],'FaceColor',[0 0 0 1-(400-elRot)/40],'EdgeColor','none')
    elseif(elRot >= 400) 
        text(0.8,-0.95,'j = l - s','FontSize',20,'HorizontalAlignment', 'center', ...
            'VerticalAlignment', 'bottom','color','k');
        rectangle('Position',[0.775,-0.55,0.05,0.2],'FaceColor',[0 0 0],'EdgeColor','none')
        rectangle('Position',[0.775,-0.75,0.05,0.2],'FaceColor',[1 0 0],'EdgeColor','none')
        rectangle('Position',[0.85,-0.55,0.025,0.1],'FaceColor',[1 0 0],'EdgeColor','none')
        rectangle('Position',[0.85,-0.65,0.025,0.1],'FaceColor',[0 0 0],'EdgeColor','none')
   end   
    
    % j = l + s energy level 
    if (elRot >= 720)
        jlps = line([0.6 1],[1 1],'Color',[0.5 0.5 0.5],'LineWidth',5);
        connect2 = plot([0.4 0.6], [0 1], 'color',[0.75 0.75 0.75],'Linewidth', 2.5); 
    end
    if (elRot >= 720) && (elRot < 760)
        jlps.Color(4) = 1 - (760-elRot)/40; % Semitransparent line
        connect2.Color(4) = 1 - (760-elRot)/40; % Semitransparent line
        textCol = [1 1 1]*(760-elRot)/40;
        text(0.8,1.05,'j = l + s','FontSize',20,'HorizontalAlignment', 'center', ...
            'VerticalAlignment', 'bottom','color',textCol);
        rectangle('Position',[0.77,1.45,0.05,0.2],'FaceColor',[0 0 0 1-(760-elRot)/40],'EdgeColor','none')
        rectangle('Position',[0.77,1.25,0.05,0.2],'FaceColor',[1 0 0 1-(760-elRot)/40],'EdgeColor','none')
        rectangle('Position',[0.85,1.35,0.025,0.1],'FaceColor',[1 0 0 1-(760-elRot)/40],'EdgeColor','none')
        rectangle('Position',[0.85,1.45,0.025,0.1],'FaceColor',[0 0 0 1-(760-elRot)/40],'EdgeColor','none')
    elseif(elRot >= 760) 
        text(0.8,1.05,'j = l + s','FontSize',20,'HorizontalAlignment', 'center', ...
            'VerticalAlignment', 'bottom','color','k');
        rectangle('Position',[0.77,1.45,0.05,0.2],'FaceColor',[0 0 0],'EdgeColor','none')
        rectangle('Position',[0.77,1.25,0.05,0.2],'FaceColor',[1 0 0],'EdgeColor','none')
        rectangle('Position',[0.85,1.35,0.025,0.1],'FaceColor',[1 0 0],'EdgeColor','none')
        rectangle('Position',[0.85,1.45,0.025,0.1],'FaceColor',[0 0 0],'EdgeColor','none')
    end
    axis off
    xlim([0 1]);
    ylim([-2 2]);
    
    frame = getframe(gcf);
    writeVideo(vid,frame);
end
% Output the movie as an mpg file
close(vid);