% Cartoon of a multielement spherical grating monochromator with the
% diffracting elements on individual Rowland circles 

clear; close all
vid = VideoWriter('multiRowland.mp4','MPEG-4');
vid.Quality = 100;
vid.FrameRate = 30;
open(vid);
figure('units','pixels','position',[0 0 2800 2000],'ToolBar','none');
set(0,'defaultfigurecolor',[1 1 1]); 

% red - yellow - green - blue - violet 
rgb = customcolormap([0 0.25 0.5 0.75 1], {'#ff0000','#ffff00','#00ff00','#0000ff','#7f00ff'});
fliprgb = flip(rgb);
colormap(fliprgb);

myBlue = [0.4 0.44 0.73];
lp = [-7 2 -2];
lp2 = [7 0 2];

rR = 0.5 % Rowland circle radius 

theta = 0:pi/180:2*pi;
circlex = 0.5*sin(theta); 
circlez = 0.5*cos(theta); 

% Set up spherical bowl for crystal analyzer surfaces
[X,Y,Z] = sphere(2000);
X = X/2; Y = Y/2; Z = Z/2; % Sphere radius 0.5
xnot = 0.496;
% Only surviving part of sphere for X > xnot - a small slice of the sphere 
X(X<xnot) = NaN ; Y(X<xnot) = NaN ; Z(X<xnot) = NaN;
cylRad = (rR^2-xnot^2)^0.5; % Radius of cylinder defining side wall of crystal analyzers
[a,b,c] = cylinder(cylRad,500);
xstht = 0.05; % Height of side wall of crystal analyzers 
detHt = 0.14;

t = 0:1;
r = t;
[xx,yy,zz] = cylinder(r,200);

pos1 = [0 0 1 1];
phiMin = 5; 
phiMax = 25; 

for phi = phiMin:0.125:phiMax % Rotate Rowland circle 
    
    colorIndex = floor(256*(phi-phiMin)/(phiMax-phiMin))+1 
    if (colorIndex > 256)
        colorIndex = 256; 
    end 
    
    Lsx = 2*rR*cosd(phi); % Sample - crystal analyzer distance
    rR0x = Lsx*cosd(phi) - rR; % x-value of centre of Rowland circle
    rR0z = -Lsx*sind(phi); % z-value of centre of Rowland circle
    del = rR - (rR^2 - cylRad^2)^0.5; % Height of lip cf middle of crystal analyzer surface
    Lax = Lsx*cosd(phi); % Distance in x-direction from middle of crystal analyzer surface to sample-detector axis
    cylOff = -rR0z + Lax + xstht - del;
    hold off
    subplot('position',pos1);
    newplot
    hold on
    axis off
    axis equal
    xlim([-0.4 1.1]);
    ylim([-0.7 0.7]);
    zlim([-0.9 0.5]);
    light('Position',lp,'Style','infinite');
    light('Position',lp2,'Style','infinite');
    view(-50,22);
    [V,F]=platonic_solid(4,0.025);
    ps1 = patch('Faces',F,'Vertices',V,'FaceColor',myBlue,'FaceAlpha',1,...
        'EdgeColor','none','FaceLighting','flat','DiffuseStrength',0.8,...
        'SpecularStrength',0.8);
    
    for xst = -3:3
        % Crystal analyzer surface 
        xstsurf = surf(X+rR0x,Y,Z+rR0z,'FaceColor',[0.4 0.4 0.4],'LineStyle','none',...
            'FaceLighting','flat','DiffuseStrength',0.4);
        rotate(xstsurf,[0 0 1],10*xst,[0,0,rR0z]); % Rotate around vertical z-axis
        % Side walls of crystal analyzers 
        xstside = surf(a,b,xstht*c - cylOff,'FaceColor',[0.8 0.8 0.8],'LineStyle','none',...
            'FaceLighting','flat','DiffuseStrength',0.4);
        rotate(xstside,[0 1 0],-90,[0,0,rR0z]);
        rotate(xstside,[0 0 1],10*xst,[0,0,rR0z]); 
        % Rowland circles 
        cc = plot3(circlex+rR0x,0*circlex,circlez+rR0z,'Color',[0.5 0.5 0.5],'LineStyle','-.','LineWidth',0.4);
        rotate(cc,[0 0 1],10*xst,[0,0,rR0z]); 
        % Fluorescent radiation cones 
        fluor1 = surf(cosd(phi)*cylRad*xx, cylRad*yy, zz*Lsx,'FaceColor',fliprgb(colorIndex,:),'FaceAlpha','0.08','LineStyle','none',...
            'FaceLighting','flat','DiffuseStrength',0.4); 
        rotate(fluor1,[0 1 0],90+phi,[0,0,0]); 
        rotate(fluor1,[0 0 1],10*xst,[0,0,0]); 
        fluor2 = surf(cosd(phi)*cylRad*xx, cylRad*yy, zz*Lsx+2*rR0z,'FaceColor',fliprgb(colorIndex,:),'FaceAlpha','0.08','LineStyle','none',...
            'FaceLighting','flat','DiffuseStrength',0.4); 
        rotate(fluor2,[0 1 0],90-phi,[0,0,2*rR0z]); 
        rotate(fluor2,[0 0 1],10*xst,[0,0,0]); 
    end
    % Detector
    detWall = surf(a/2,b/2,c*detHt - detHt + 2*rR0z,'FaceColor',[0.8 0.8 0.8],'LineStyle','none',...
            'FaceLighting','flat','DiffuseStrength',0.7);
    rotate(detWall,[0 1 0],90-phi,[0,0,2*rR0z]);
    plot3([0 0],[0 0],[-0.9 0.1],'k','LineStyle','-.','LineWidth',1.6);
    detEnd = fill3(cylRad*circlex,cylRad*circlez,0*circlez+2*rR0z-detHt,[0.5 0.5 0.5],'LineStyle','none',...
        'FaceLighting','flat','DiffuseStrength',0.5,'SpecularStrength',0.5); 
    rotate(detEnd,[0 1 0],90-phi,[0,0,2*rR0z]);
    
    frame = getframe(gcf);
    writeVideo(vid,frame);
    
    %hold off
end
% Output the movie as an mpg file
close(vid);