% Program to generate a cartoon of Bragg's law

clear; close all;

figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
vid = VideoWriter('braggsLaw1.mp4','MPEG-4');
vid.Quality = 100;
vid.FrameRate = 30;
open(vid);
set(0,'defaultfigurecolor',[1 1 1]);
set(gca,'linewidth',7);
set(gca, 'Projection','perspective');
lp1 = [-0.4 -10.7 1.4];
lp2 = [0.4 -0. 1.4];
myBlue = [0.3 0.34 0.61];
myGold = [0.7 0.6 0];
myPurple = [0.5 0 0.5];
axis equal

LL = 4;
d = 1;
lambda = 0.5;

xlim([-LL LL]);
ylim([-LL LL]);
zlim([-LL LL]);


[x,y,z] = sphere(40); % Create sphere coordinates
ar = 0.25; % radius of atoms

thetaStart = 20;
thetaEnd = 40;
thetaStep = 0.1;
nthsteps = (thetaEnd-thetaStart)/thetaStep;

fwd = [thetaStart:thetaStep:thetaEnd];

for theta = [fwd flip(fwd)] % from min angle to max and back again 
    hold off
    for jj = 0:3
        for ii=-3:3
            atom = surf(ii*d+x*ar,y*ar,z*ar-jj*d,'FaceAlpha',(4.7-jj)/5,'FaceColor',myBlue, ...
                'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1, ...
                'AmbientStrength',0.1);
            hold on
        end
        set(gca, 'Projection','perspective');
        set(gca,'View',[0,0]);
        axis equal
        xlim([-LL LL]);
        ylim([-LL LL]);
        zlim([-LL LL]);
        
        % incident x-rays at angle theta to atomic planes
        rayLength = 2*LL+jj*d*sind(theta);
        braggAng = asind(2*lambda/(2*d));
        offdiff = abs((braggAng-theta)/braggAng);
        if (offdiff <= 0.01)
            myCol = [1 0 0];
            myCol2 = myCol;
        else
            myCol = myGold;
            myCol2 = [0 0 0];
        end
        wavex = 0:LL/200:rayLength;
        wavey = 0.0*wavex-0.25;
        wavez = 0.2*sin(2*pi*(wavex/lambda))-jj;
        h = plot3(wavex-rayLength,wavey,wavez,'color',myCol, 'LineWidth', 3.5);
        rotate(h,[0 1 0],theta,[0,-0.25,-jj]);
        h.Color(4) = (4.7-jj)/5; % Semitransparent line
        
        % x-rays scattered by 2theta
        wavex = 0:LL/200:2*LL+jj*sind(theta);
        wavez = 0.2*sin(2*pi*(wavex+jj*d*sind(theta))/lambda);
        h = plot3(wavex,wavey,wavez-jj*d,'color',myCol, 'LineWidth', 3.5);
        rotate(h,[0 1 0],-theta,[0,-0.25,-d*jj]);
        h.Color(4) = (4-jj)/5; % Semitransparent line
        
        % Peak-indicating dot-dashed lines showing phase differences
        phline = plot3([2-3*lambda/4-jj*d*(sind(theta)-1) ...
            2-3*lambda/4-jj*d*(sind(theta)-1)],[-0.25 -0.25], ...
            [0.5-d*jj -0.5-d*jj],'color',myCol2, 'LineWidth',1.4,'LineStyle','-.');
        rotate(phline,[0 1 0],-theta,[0,-0.25,-d*jj]);
        
        % Plot Gaussian to represent Bragg peak signal
        gaussx = [-0.5:1/nthsteps:0.5];
        gaussy = 0*gaussx;
        gaussz = 1+exp(-gaussx.^2/0.01);
        gaIndex = (theta-thetaStart)/thetaStep+1;
        gauss = plot3(gaussx(1:gaIndex),gaussy(1:gaIndex),gaussz(1:gaIndex), ...
            'color',myCol2,'LineWidth',2);
        hText = text(0, 0, 0.83,'2\theta', 'FontSize',18,'HorizontalAlignment', 'center');
        
        
        % In and out arrows
        arrow1 = mArrow3([-LL -0.4 -1.5*d],[-LL+1 -0.4 -1.5*d],'color',myGold, ...
            'stemWidth',0.05,'tipWidth',0.1,'facealpha',1); %
        rotate(arrow1,[0 1 0],theta,[0,-0.4,-1.5*d]);
        arrow2 = mArrow3([LL-1 -0.4 -1.5*d],[LL -0.4 -1.5*d],'color',myGold, ...
            'stemWidth',0.05,'tipWidth',0.1,'facealpha',1); %
        rotate(arrow2,[0 1 0],-theta,[0,-0.4,-1.5*d]);
        
    end
    light('Position',lp1,'Style','infinite');
    light('Position',lp2,'Style','infinite');
    
    axis off
    frame = getframe(gcf);
    writeVideo(vid,frame);
    
end

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