% Cartoon movie of refraction of visible light by a prism

clear; close all;
vid = VideoWriter('prismRefraction.mp4','MPEG-4');
vid.Quality = 100;
vid.FrameRate = 30;
open(vid);
figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
set(0,'defaultfigurecolor',[0 0 0]); % black background = [0 0 0]. Normally white [1 1 1]
set(gcf,'color',[0 0 0]);
set(gca,'linewidth',7);
LL = 1;
ampl = 0.01;
nR = 1.5;

% Incident point of white beam on prism
incPtx = -LL/4;
incPty = -sqrt(3)/4;
% Angle of incident beam relative to horizontal
incAng = 18;
% End of straight line defined by incident beam behind vertical axis of prism
endIncx = 0;
endIncy = incPty+tand(incAng)*LL/4;

% Define 7 colors of dispersed rays
rayColors = [1 0 0; 1 0.5 0; 1 1 0; 0 1 0; 0 0.5 1; 0 0 1; 0.5 0 0.5];
for t = 0:2*pi/10:120*pi-pi/10
    newplot
    axis square
    axis equal
    axis tight
    ax = gca;
    ax.Visible = 'off';
    hold on
    xlim = ([-2*LL 2*LL]);
    ylim = ([-2*LL 0.5*LL]);
    
    alphaTemp = 0;
    hold on
    % Draw incident white beam
    for lw = 42:-4:2
        alphaTemp = alphaTemp+0.02;
        lineAlpha = exp(-lw^2/200);
        incRay = plot3([endIncx,endIncx-0.7*LL],[endIncy,endIncy-0.7*LL*tand(incAng)],...
            [-lw/1000,-lw/1000],'Color','w','LineWidth',lw);
        incRay.Color(4) = alphaTemp;
    end
    % Draw grey prism
    pgon = polyshape([0 0.88*LL/2 -0.88*LL/2],[0 -0.88*LL*sqrt(3)/2 -0.88*LL*sqrt(3)/2]);
    prism = plot(pgon);
    prism.FaceColor = [0.1 0.1 0.1];
    prism.FaceAlpha = 1;
    prism.EdgeColor = [0.2 0.2 0.2];
    prism.LineWidth = 2.5;
    
    for rayCol = 1:7
        hold on
        % Calculate where each dispersed ray passes out through the prism
        dispAng = 2*(7 - rayCol);
        m = tand(dispAng); % gradient of dispersed ray
        c = incPty - m*incPtx; % constant of y = mx + c describing ray
        % Now determine where ray passes out of prism
        % Given by where y = mx + c and y = -sqrt(3)x cross
        endPrisx = -c/(sqrt(3)+m);
        endPrisy = -sqrt(3)*endPrisx;
        % Now determine length of dispersed ray in prism and make an array 1000
        % steps long of this length
        dispRayLength = ((endPrisx-incPtx)^2 + (endPrisy-incPty)^2)^0.5;
        dispRay = (0:dispRayLength/1000:dispRayLength);
        lambda = 0.064*(14-rayCol)/14;
        lambdaRed = 0.064*6/7;
        k = 2*pi/lambda;
        dispWavey = ampl*sin(k*dispRay - lambdaRed*t/lambda);
        ray1 = plot(dispRay+incPtx,dispWavey+incPty,'LineWidth',2.5,'Color',rayColors(rayCol,:));
        rotate(ray1,[0 0 1],dispAng,[incPtx,incPty,0]);
        
        % Determine phase of wave as it exits prism
        phiPrism = 2*pi*mod(dispRayLength,lambda)/lambda
        % shift in start of ray emerging from prism
        delStart = -ampl*sin(phiPrism)*sind(2.6*dispAng-incAng);
        outRayLength = 0.7*LL-(endPrisx+delStart);
        outRay = (0:outRayLength/1000:outRayLength);
        
        % Draw dispersed waves after emerging from prism
        lambda = nR*0.064*(14-rayCol)/14;
        k = 2*pi/lambda;
        outWavey = ampl*sin(phiPrism+k*dispRay - nR*lambdaRed*t/lambda);
        ray2 = plot(outRay+(endPrisx+delStart),outWavey+endPrisy,...
            'LineWidth',2.5,'Color',rayColors(rayCol,:));
        rotate(ray2,[0 0 1],1.6*dispAng-incAng,[endPrisx+delStart,endPrisy,0]);
    end
    % Store the frame
    frame = getframe(gcf);
    writeVideo(vid,frame);
    hold off
end
% Output the movie as an mpg file
close(vid);


