% Cartoon of CXDI in forward scattering direction. Scattered object is an
% ellipsoid with semi-axes of 11, 5, and 5 scattering centers and a random
% electron-density distribution - kind of like a paramecium? 

clear
close all
figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
set(0,'defaultfigurecolor',[1 1 1]);

vid = VideoWriter('CXDI.mp4','MPEG-4');
vid.Quality = 100;
vid.FrameRate = 30;
open(vid);

myGold = [0.8 0.64 0];
lp = [2 3 1];
lp2 = [1 -2 -1];

% Set up coordinates of front cap of Ewald sphere - I don't want to plot
% the entire sphere, as it is too large and its semiopaque nature would 
% compromise the interesting part of the cartoon 
[Xe,Ye,Ze] = sphere(2000);
ESrad = 500; % Ewald sphere radius
Xe = ESrad*Xe - ESrad; % Shift sphere so front of sphere at origin (0,0,0) 
Ye = ESrad*Ye; 
Ze = ESrad*Ze; 
xnot = -80; % Value of x, below which the Ewald sphere is removed 
% Only surviving part of sphere for Xe > xnot - a small slice of the sphere
Xe(Xe<xnot) = NaN ; Ye(Xe<xnot) = NaN ; Ze(Xe<xnot) = NaN;

% Collection of random numbers to be sculpted into an ellipsoidal shape 
X = rand(23,11,11); 
XX = X; % Use XX as replica for FT purposes (FT cannnot handle NaNs)

% Shave off elements outside ellipsoidal surface. For the FT (XX), make 
% these zero; for plotting (X), make these NaN 
for ii = 1:23
    for jj = 1:11
        for kk = 1:11
            % Radius of ellipsoid
            R = (((ii-12)/12)^2 + ((jj-6)/6)^2 + ((kk-6)/6)^2)^0.5;
            if (R > 1)
                X(ii,jj,kk) = NaN;
                XX(ii,jj,kk) = 0;
            end
            hold on
        end
    end
end 
% Reshape matrices of X for scatter3 plotting purposes 
[xR, yR, zR] = meshgrid(1:size(X,2), 1:size(X,1), 1:size(X,3));
xxR = reshape(xR,[],1)-6;
yyR = reshape(yR,[],1)-12;
zzR = reshape(zR,[],1)-6;
X = reshape(X,[11*11*23,1]);

% ________________________________________________________________________

% Generate FT of XX padded to pad^3 elements 
pad = 100;
Y = abs(fftn(XX,[pad,pad,pad]));
% Shift FT so maximum is in the middle 
Y = fftshift(Y); 
maxY = max(Y(:));

% Manipulate FT data for plotting purposes 
for ii = 1:pad
    for jj = 1:pad
        for kk = 1:pad
            R = ((ii-(pad+1)/2)^2 + (jj-(pad+1)/2)^2 + (kk-(pad+1)/2)^2)^0.5; 
            % Only plot out FT to a radius equal to pad/2 and for signal
            % stronger than 2% of the FT maximum 
            if (Y(ii,jj,kk) < 0.02*maxY) || (R > pad/2)
                Y(ii,jj,kk) = NaN;
            end
        end
    end
end 
% Reshape manipulated FT data for plotting with scatter3 
[x, y, z] = meshgrid(1:size(Y,1), 1:size(Y,2), 1:size(Y,3));
xx = reshape(x,[],1)-pad/2;
yy = reshape(y,[],1)-pad/2;
zz = reshape(z,[],1)-pad/2;
Y = reshape(Y,[pad^3,1]); 
cmin = min(log(Y));
cmax = max(log(Y));

% ________________________________________________________________________

% Loop for rotation of object and its FT by 360 degrees 
for theta = 0:359
    subplot('position',[-0.2 -0.2 1.4 1.55]);
    hold off
    newplot
    hold off 
    % Coordinates at theta of rotating paramecium ellipsoid 
    xxRnow = 3*(xxR*cosd(theta) - yyR*sind(theta))-500;
    yyRnow = 3*(xxR*sind(theta) + yyR*cosd(theta));
    hold on 
    
    % Plot ellipsoid paramecium 
    scatter3(xxRnow,yyRnow,3*zzR,50,X*cmax,'filled','MarkerFaceAlpha',... 
        0.25,'MarkerEdgeAlpha',0)    
    colormap jet
    
    % Incident parallel radiation on object 
    mArrow3([-640 0 0],[-500 0 0],...
        'color',myGold,'stemWidth',2.5,...
        'tipWidth',7,'facealpha',0.25);
    
    % Plot Ewald sphere
    ewaldSurf = surf(Xe,Ye,Ze,'FaceColor',myGold,'LineStyle','none',...
        'FaceLighting','flat','DiffuseStrength',1,'SpecularStrength',0.55);
    % Vary alpha of Ewald sphere in x-direction with offset of 50 so it
    % remains semitransparent for all elements of surf = ewaldSphere 
    alpha(ewaldSurf,ewaldSurf.XData-50);    
    
    % Scattering pattern - transform x- and y-coordinates according to
    % rotation angle theta then plot out with scatter3 
    xxnow = xx*cosd(theta) - yy*sind(theta);
    yynow = xx*sind(theta) + yy*cosd(theta);
    
    scatter3(xxnow,yynow,zz,Y.^0.5/5,log(Y(:)),'filled',... 
        'MarkerFaceAlpha',0.2,'MarkerEdgeAlpha',0)
    colormap turbo
    caxis([cmin cmax])
    
    % Plot part of scattering pattern that sits on the Ewald sphere 
    % highlighted in red 
    clear xxEw yyEw zzEw YEw 
    mm = 1;
    for ii = 1:size(xx)
        if (abs(xxnow(ii)) < 0.64)
            xxEw(mm) = xxnow(ii);
            yyEw(mm) = yynow(ii);
            zzEw(mm) = zz(ii);
            YEw(mm) = Y(ii);
            mm = mm+1;
        end
    end
    scatter3(xxEw,yyEw,zzEw,YEw/10,'r','filled','MarkerFaceAlpha',1,... 
        'MarkerEdgeAlpha',0)
    
    % Cone of forward-scattered x-rays
    t = 0:1; 
    % Generate cone coordinates with base radius = pad/2 
    [xC,yC,zC] = cylinder(t,pad/2); 
    % Transparency of cone representing scattered radiation proportional to
    % the sum of the scattered radiation that sits on the Ewald sphere, YEw
    coneAlp = sum(YEw,"omitnan")/2e6; 
    coneRad = surf(50*xC-510,50*yC,505*zC,'FaceAlpha',coneAlp,... 
        'FaceColor',myGold,'LineStyle','none','FaceLighting','flat',... 
        'DiffuseStrength',1);
    rotate(coneRad,[0 1 0],90,[-510,0,0]);
    
    % Finish off defining plotting appearance 
    axis off
    axis equal
    xlim([-640 70]);
    ylim([-280 280]);
    zlim([-280 280]);
    set(gca, 'Projection','perspective'); 
    view(35,22)
    light('Position',lp,'Style','infinite');
    light('Position',lp2,'Style','infinite'); 
    % Capture frame and write to video 
    frame = getframe(gcf);
    writeVideo(vid,frame);
end 
close(vid);
