% Cartoon movie of electric-field vectors induced by a moving charge 
% This version moves the electron perpendicularly out of the plane of the field grid 

clear; close all;
vid = VideoWriter('electricField1.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);
lp = [2 1 1];
myBlue = [0.5 0.55 0.88];
myBlue2 = [0.61 0.67 0.92];
pos1 = [-0.1 -0.2 1.2 1.5];

[xx,yy,zz] = sphere(70); 
r = 0.5; % Radius of electron
ampl = 1.25;
lambda = 16; % Wavelength in terms of grid dimensions
% Want c to assume a value so that waves move one wavelength/second
c = lambda/frmRate; % c = f lambda; f = 1s = frmRate frames

xMx = 28;
yMx = 28;
x = linspace(-xMx,xMx,2*xMx+1);
y = linspace(-yMx,yMx,2*yMx+1);
[X,Y] = meshgrid(x,y);
Z = 0.0*X;
R = (X.^2 + Y.^2).^0.5; % Matrix of distances of each element in the meshgrid from (0,0)
phi = atan2(Y,X); % In-plane angle of base of E-field vectors in x-y plane

subplot('position',pos1);
for frm = 0:29 % Time in units of 1/frmRate
    newplot
    hold on
    plot3(X,Y,Z,'.y')
    axis equal
    axis off
    set(gca, 'Projection','perspective');
    set(gca,'View',[44,40]);
    set(gcf,'Color','k');
    campos([43 41 25]);
    xlim([-xMx-0.1 xMx+0.1]);
    ylim([-yMx-0.1 yMx+0.1]);
    zlim([-4*ampl 4*ampl]);

    accz = sin(2*pi*c*frm/lambda); % Acceleration of electron in z-direction
    posz = -sin(2*pi*c*frm/lambda); % Position of electron in z-direction
    theta = atan2(posz,R); % Out-of-plane angle of E-field vectors away from vertical
    accperp = -ampl*sin(2*pi*c*(frm - abs(R)/c)/lambda).*R./(posz^2 + R.^2).^0.5;
    arrowsEnds = accperp.*cos(theta); % z-position of arrow end
    Efield = accperp./R;
    delxy = accperp.*sin(theta); % In-plane shift of arrow end w.r.t. arrow start
    delx = delxy.*cos(phi);
    dely = delxy.*sin(phi);
    for ii = 1:2*xMx+1
        for jj = 1:2*yMx+1
            if (X(jj,ii) ~= 0) || (Y(jj,ii) ~= 0)
                arrowStartx = X(jj,ii);
                arrowStarty = Y(jj,ii);
                arrowStartz = Z(jj,ii);
                arrowEndx = arrowStartx + delx(jj,ii);
                arrowEndy = arrowStarty + dely(jj,ii);
                arrowEndz = arrowsEnds(jj,ii);
                % Determine whether arrow has nonzero length because
                % mArrow3 totally craps out if it is zero. Took me forever
                % to debug!! 
                delxa = arrowEndx - arrowStartx; 
                delya = arrowEndy - arrowStarty; 
                delza = arrowEndz - arrowStartz; 
                arrowLen = (delxa^2 + delya^2 + delza^2)^0.5; 
                if (arrowLen > 0)
                h = mArrow3([arrowStartx arrowStarty arrowStartz],...
                    [arrowEndx arrowEndy arrowEndz],...
                    'color',myBlue2,'stemWidth',0.07,'tipWidth',0.1,...
                    'facealpha',(abs(Efield(jj,ii))/ampl)^0.35);
                end 
            end
        end
    end 
    % Oscillating electron 
    surf(xx*r,yy*r,zz*r+posz,'FaceColor',myBlue,'EdgeColor','none');
    light('Position',lp,'Style','infinite');

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

