% Cartoon movie of electric-field vectors induced by a moving charge
% This version moves the electron in a circular path, inducing circular
% polarization
clear; close all;
vid = VideoWriter('electricField3.mp4','MPEG-4');
vid.Quality = 100;
frmRate = 30;
vid.FrameRate = frmRate;
open(vid);
figure('units','pixels','position',[0 0 3840 2160],'ToolBar','none');
set(0,'defaultfigurecolor',[1 1 1]); % white background = [1 1 1]
set(gca,'linewidth',7);

% Read in image files for white x, y, RCP, and LCP labels 
[whitex,~,alphaMap] = imread('whitex.png');
[whitey,~,alphaMap2] = imread('whitey.png');
[RCP,~,alphaMap3] = imread('RCP.png');
[LCP,~,alphaMap4] = imread('LCP.png');


lp = [2 1 1];
myBlue = [0.5 0.55 0.88];
myBlue2 = [0.61 0.67 0.92];
pos1 = [-0.4 -0.4 1.91 1.9];

xMx = 32;
yMx = 32;

[xx,yy,zz] = sphere(70);
r = 0.5; % Radius of electron
ampl = 1.25; % Radius of circular path
Eampl = 4; % Maximum amplitude of E-field

% Vector of alpha values for fading in and out of the RCP and LCP labels 
alpRCPLCP = [zeros(1,150) 1/75:1/75:1 ones(1,150) 1:-1/75:1/75 zeros(1,150)];

[xa,ya,za] = cylinder([0.2 0]); % Used to generate the axis arrow heads 

lambda = 20; % Wavelength in terms of grid dimensions

x = -xMx:0.5:xMx;

subplot('position',pos1);
for frm = 0:599 % Time in units of 1/frmRate
    newplot
    hold on
    axis equal
    axis off
    set(gca, 'Projection','perspective');
    set(gcf,'Color','k'); 
    % Varying camera position 
    campos([69*cos(2*pi*frm/4800+pi/6) 67*sin(2*pi*frm/4800+pi/6) 29]);
    xlim([-xMx-2 xMx+2.1]);
    ylim([-yMx-2 yMx+2]);
    zlim([-7 7]); 
    % Arrow heads 
    axisArrowx = surf(xa,ya,za+xMx+0.25,'FaceColor','w','LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1);
    axisArrowy = surf(xa,ya,za+xMx,'FaceColor','w','LineStyle','none',...
        'FaceLighting','gouraud','DiffuseStrength',1);

    % Draw x- and y-axes
    plot3([-xMx xMx+0.25],[0 0],[0 0],...
        'color','w', 'LineWidth', 1);
    tickLimx = xMx - mod(xMx,5);
    for ii = -tickLimx:5:tickLimx
        plot3([ii ii],[-0.25 0.25],[0 0],...
            'color','w', 'LineWidth', 1);
    end
    rotate(axisArrowx,[0 1 0],90,[0,0,0]);
    surf([xMx+1.4 xMx+2.1; xMx+1.4 xMx+2.1],[0 0; 0 0],[-0.35 -0.35; 0.35 0.35],...
        whitex,'FaceColor','texturemap','FaceAlpha','texturemap',...
        'AlphaData',alphaMap,'EdgeColor','none','AlphaDataMapping',...
        'none','FaceLighting','none');

    plot3([0 0],[-yMx yMx],[0 0],...
        'color','w', 'LineWidth', 1);
    tickLimy = yMx - mod(yMx,5);
    for ii = -tickLimy:5:tickLimy
        plot3([-0.25 0.25],[ii ii],[0 0],...
            'color','w', 'LineWidth', 1);
    end
    rotate(axisArrowy,[1 0 0],-90,[0,0,0]);
    surf([0 0; 0 0],[yMx+1.2 yMx+1.9; yMx+1.2 yMx+1.9],[0.35 0.35; -0.41 -0.41],...
        whitey,'FaceColor','texturemap','FaceAlpha','texturemap',...
        'AlphaData',alphaMap2,'EdgeColor','none','AlphaDataMapping',...
        'none','FaceLighting','none');

    % Oscillating electron and acceleration arrow
    surf(xx*r,yy*r + ampl*sin(2*pi*frm/60),zz*r - ampl*cos(2*pi*frm/60),...
        'FaceColor',myBlue,'EdgeColor','none');
    elAccArrow = mArrow3([0 0 -ampl],...
        [0 0 3*ampl],'color','r','stemWidth',0.125,'tipWidth',0.25);
    rotate(elAccArrow,[1 0 0],360*frm/60,[0,0,0]);
    % R, the distance of electric-field position along x-axis from present
    % position of the electron
    R = (1 + x.^2).^0.5;
    % Angle of R relative to x-axis
    theta = atan2(1,x);
    accperp = cos(theta); % Acceleration of electron perpendicular to R 

    % Draw E-field arrows 
    for ii = 1:4*xMx+1
        if (accperp(ii) ~= 0)
            arrowStartx = x(ii);
            arrowStarty = 0;
            arrowStartz = 0;
            arrowEndx = x(ii);
            arrowEndy = 0;
            arrowEndz = -Eampl*accperp(ii);
            h = mArrow3([arrowStartx arrowStarty arrowStartz],...
                [arrowEndx arrowEndy arrowEndz],...
                'color',myBlue2,'stemWidth',0.07,'tipWidth',0.125);
            % Rotate around y to get perpendicular angle 
            rotate(h,[0 1 0],-(180/pi)*theta(ii),[x(ii),0,0]); 
            % Rotate around x to generate circular polarization 
            rotate(h,[1 0 0],360*(frm/60 - R(ii)/lambda),[x(ii),0,0]);
        end
    end
    alphaMap3a = im2double(alphaMap3)*alpRCPLCP(1,frm+1);
    RCPlabel = surf([xMx/2+4 xMx/2; xMx/2+4 xMx/2],[0 0; 0 0],[6.5 6.5; 5 5],...
        RCP,'FaceColor','texturemap','FaceAlpha','texturemap',...
        'AlphaData',alphaMap3a,'EdgeColor','none','AlphaDataMapping',...
        'none','FaceLighting','none');
    alphaMap4a = im2double(alphaMap4)*alpRCPLCP(1,frm+1);
    LCPlabel = surf([-xMx/2 -xMx/2-4; -xMx/2 -xMx/2-4],[0 0; 0 0],[6.5 6.5; 5 5],...
        LCP,'FaceColor','texturemap','FaceAlpha','texturemap',...
        'AlphaData',alphaMap4a,'EdgeColor','none','AlphaDataMapping',...
        'none','FaceLighting','none');

    light('Position',lp,'Style','infinite');

    frame = getframe(gcf);
    writeVideo(vid,frame);
    hold off
end
% Output the movie as an mpg file
close(vid);
