% Cartoon of HTSU10 undulator. This cartoon shows how the flux-density
% curves (5 - 60 keV) change as the magnetic-field amplitude between the 
% magnet arrays on axis increases. This is accompanied by an increase in 
% integrated intensity (how opaque the x-ray beam is); an increase in 
% electron-beam oscillation amplitude; an increase in solenoid current 
% (shown by the ammeter); and a reddening of the HTS magnet half-moon 
% blocks. The spectra are plotted with a logarithmic ordinate axis. Note 
% that in practice, one never applies the solenoid field during operation 
% of the storage ring. 
% 
% The parameters used to generate the undulator spectra in SPECTRA are: 
% Storage-ring energy = 2.7 GeV 
% Current = 400 mA 
% Circumference = 288 m 
% Bunches = 480 
% Natural emittance = 150 pm.rad 
% Coupling constant = 0.002 
% Energy spread = 0.0011 
% beta_x,y = 2.5, 1.28 
% All other accelerator constants = 0 
% Linear undulator 
% B = 0.01 - 2 Tesla in 0.01-T steps (this appears not to be a variable 
% parameter, hence I took 200 linear steps in K from 0.01 to 1.86746, as 
% this maximum value corresponds to B = 2 T). 
% lambda_u = 10 mm 
% Device length = 1 m
% End Correction Magnet: activated 
% Far field and ideal conditions: energy dependence of partial flux through
% a rectangular slit 1 x 1 mm2 at 30 m, slit position @ 0,0. 
% 5000 - 60000 eV, 1-eV steps 
% 
% Many thanks go to Marco Calvi, PSI for assistance in generating the
% spectra. 


clear; clc

vid = VideoWriter('HTSU10a.mp4','MPEG-4');
vid.Quality = 100;
vid.FrameRate = 30;
open(vid);
figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
set(0,'defaultfigurecolor',[1 1 1]);
set(gca,'linewidth',7);
lp = [2 -1 0.4]; % Angle of lighting
lp2 = [0 1 0]; % Second angle of lighting
grey = [0.64 0.64 0.64];
myBlue = [0.4 0.44 0.73];
myGold = [1.0 0.83 0];
myCream = [1.0 0.99 0.95];
% Custom rainbow color map for undulator spectrum
rgb = customcolormap([0 0.14 0.28 0.43 0.57 0.71 0.86 1],...
    {'#800000','#ff0000', '#ffa500', '#ffff00', '#00ff00',...
    '#00ffff', '#0000ff', '#7f00ff'},55001);
rgb = flip(rgb);
% Custom color map for increasing activation of magnets by solenoid B-field
greyRed = customcolormap([0 1], [0.64 0.64 0.64; 0.7 0 0],200);
greyRed = flip(greyRed);

% Maximize the area plotted
pos1 = [-0.1 -0.55 1.25 1.9];

str1 = 'HTSU10/P_'; % Folder HTSU10, files beginning with 'P_'
str3 = '.txt';

% Find max and min intensity of all ID spectral data
% Max
filename = [str1 '200' str3];
IDspect = importdata(filename,' ',2); % Ignore first two lines
IDspectrum = IDspect.data;
maxY = max(IDspectrum(:,2));

subplot('position',pos1);
for jj = 1:200
    newplot
    % -------------------------------------------------------------------------
    % Create coordinates for solenoid and draw it
    r1 = 1; % Inner radius of winding
    r2 = 1.02; % Outer radius of winding
    rWire = (r2 - r1)/2; % Radius of winding wire
    L = 5; % Length of solenoid
    loop = 50.5; % Number of windings
    nNodes = 10; % Resolution of drawing the wire surface
    nCS = 90*loop; % Number of cross-sections of wire going to make up coiled surface

    totAng = loop*360; % Total angle subtended by winding the wire
    theta = linspace(0,2*pi,nNodes);
    r = 0.5*(r2-r1)*ones(1,nNodes);
    z = r.*cos(theta);
    y = r.*sin(theta) + r1 + rWire;
    x = zeros(1,nNodes);

    dz = linspace(0,L,nCS) - L/2;
    phi = linspace(0,totAng,nCS);
    X = [];
    Y = [];
    Z = [];
    for i = 1:nCS
        temp_CS = rotz(phi(i)) * [x;y;z];
        X = [X;temp_CS(1,:)];
        Y = [Y;temp_CS(2,:)];
        Z = [Z;temp_CS(3,:) + dz(i)];
    end
    solenoid = surf(X,Y,Z,'lineStyle','none','faceColor',grey,...
        'FaceAlpha',0.7,'FaceLighting','gouraud');
    rotate(solenoid,[0 1 0],90,[0,0,0]);

    hold on

    % Straight end wire of windings
    wireL = 1.6;
    [xW,yW,zW] = cylinder(rWire,nNodes);
    xW = xW+L/2;
    yW = yW-(r2+r1)/2;
    zW = -wireL*zW;
    endWire = surf(xW,yW,zW,'lineStyle','none','faceColor',grey,...
        'FaceAlpha',0.7,'FaceLighting','gouraud');

    % Ammeter
    amHt = 0.2;
    [xA,yA,zA] = cylinder(amHt,180); % Curved surface
    xA = xA + L/2;
    yA = yA - (r2+r1)/2;
    zA = amHt*zA - wireL/2 - amHt/2-0.01;
    ammeter = surf(xA,yA,zA,'lineStyle','none','faceColor',grey,...
        'FaceLighting','flat','DiffuseStrength',0.7);
    rotate(ammeter,[1 0 0],90,[L/2,-(r2+r1)/2,-wireL/2])
    % Ammeter face
    th = 0:360;
    xF = 0.999*amHt*cosd(th)+L/2;
    zF = 0.999*amHt*sind(th)-wireL/2;
    yF = 0.0*xF-r1-amHt/2+0.02;
    fill3(xF,yF,zF,myCream,'lineStyle','none',...
        'FaceLighting','none','DiffuseStrength',1)
    for ii = -120:20:120
        plot3([0.19*sind(ii) 0.16*sind(ii)]+L/2,...
            [-r1-amHt/2+0.019 -r1-amHt/2+0.019],...
            [0.19*cosd(ii) 0.16*cosd(ii)]-wireL/2,...
            'Color','k','lineWidth',2.5)
    end
    % Write 'A'
    plot3([L/2 L/2-0.035],[-r1-amHt/2+0.019 -r1-amHt/2+0.019],...
        [-wireL/2-0.05 -wireL/2-0.1],'Color','k','lineWidth',2)
    plot3([L/2 L/2+0.035],[-r1-amHt/2+0.019 -r1-amHt/2+0.019],...
        [-wireL/2-0.05 -wireL/2-0.1],'Color','k','lineWidth',2)
    plot3([L/2-0.0175 L/2+0.0175],[-r1-amHt/2+0.019 -r1-amHt/2+0.019],...
        [-wireL/2-0.075 -wireL/2-0.075],'Color','k','lineWidth',2)
    % Ammeter pointer
    dialArrow = mArrow3([L/2,-r1-amHt/2,-wireL/2-0.1],...
        [L/2,-r1-amHt/2,-wireL/2+0.14],...
        'tipWidth',0.0125,'color','b','stemWidth',0.007);
    rotate(dialArrow,[0 1 0],-120+(jj*240/200),[L/2,-r1-amHt/2,-wireL/2])
    [a,b,c] = sphere(20);
    surf(0.02*a+L/2,0.01*b-r1-amHt/2,0.02*c-wireL/2,'lineStyle',...
        'none','faceColor','k','FaceLighting','gouraud','DiffuseStrength',1)

    % -------------------------------------------------------------------------
    % Create HTS half moons (not quite half moons)
    theta = 16:164; % Not quite 180 degrees for main semicircle
    phi = -90:16; % Range of rounding angles
    moonRad = 0.4; % Radius of HTS insert
    cornRad = 0.04; % Radius of rounded corners
    moonWidth = 0.125; % Width of HTS insert
    dum = 0.0*theta;
    % Flat face
    cornOriginy = -(moonRad-cornRad)*cosd(16); % y-position of origin of rounded corner
    cornOriginz = (moonRad-cornRad)*sind(16); % z-position of origin of rounded corner
    yCorn = -cornRad*cosd(phi)+cornOriginy;
    xCorn = 0.0*yCorn;
    zCorn = cornRad*sind(phi)+cornOriginz;
    yMoon = [yCorn -moonRad*cosd(theta) -flip(yCorn)];
    xMoon = [xCorn 0.0*theta xCorn];
    zMoon = [zCorn moonRad*sind(theta) flip(zCorn)];
    % Curved face
    [x2,y2,z2] = cylinder(moonRad,360);
    for kk = 1:2
        for ll = 1:361
            if (x2(kk,ll)<cornOriginz)
                x2(kk,ll) = NaN;
                y2(kk,ll) = NaN;
                z2(kk,ll) = NaN;
            end
        end
    end
    [x3,y3,z3] = cylinder(cornRad,90);

    for ii = -floor(L/2):2*moonWidth:floor(L/2)
        % Top magnets
        fill3(xMoon+ii-moonWidth,yMoon,zMoon,greyRed(jj,:),'lineStyle',...
            'none') % Back face top magnets
        fill3(xMoon+ii,yMoon,zMoon,greyRed(jj,:),'lineStyle',...
            'none') % Front face top magnets
        moonMain = surf(x2,y2,z2*moonWidth+ii,'FaceAlpha',1,'FaceColor',...
            greyRed(jj,:),'LineStyle',...
            'none','FaceLighting','gouraud','DiffuseStrength',1);
        rotate(moonMain,[0 1 0],-90,[0 0 0])
        moonCorn1 = surf(x3+cornOriginz,y3+cornOriginy,z3*moonWidth+ii,...
            'FaceAlpha',1,'FaceColor',greyRed(jj,:),'LineStyle',...
            'none','FaceLighting','gouraud','DiffuseStrength',1);
        rotate(moonCorn1,[0 1 0],-90,[0 0 0])
        moonCorn2 = surf(x3+cornOriginz,y3-cornOriginy,z3*moonWidth+ii,...
            'FaceAlpha',1,'FaceColor',greyRed(jj,:),'LineStyle',...
            'none','FaceLighting','gouraud','DiffuseStrength',1);
        rotate(moonCorn2,[0 1 0],-90,[0 0 0])

        % Bottom magnets
        fill3(xMoon+ii,yMoon,-zMoon,greyRed(jj,:),'lineStyle',...
            'none') % Back face top magnets
        fill3(xMoon+ii+moonWidth,yMoon,-zMoon,greyRed(jj,:),'lineStyle',...
            'none') % Front face top magnets
        moonMain = surf(-x2,y2,z2*moonWidth+ii-moonWidth,'FaceAlpha',1,...
            'FaceColor',greyRed(jj,:),'LineStyle',...
            'none','FaceLighting','gouraud','DiffuseStrength',1);
        rotate(moonMain,[0 1 0],-90,[0 0 0])
        moonCorn1 = surf(-(x3+cornOriginz),y3+cornOriginy,...
            z3*moonWidth+ii-moonWidth,'FaceColor',greyRed(jj,:),'LineStyle',...
            'none','FaceLighting','gouraud','DiffuseStrength',1);
        rotate(moonCorn1,[0 1 0],-90,[0 0 0])
        moonCorn2 = surf(-(x3+cornOriginz),y3-cornOriginy,...
            z3*moonWidth+ii-moonWidth,'FaceColor',greyRed(jj,:),'LineStyle',...
            'none','FaceLighting','gouraud','DiffuseStrength',1);
        rotate(moonCorn2,[0 1 0],-90,[0 0 0])
        fill3(ii+moonWidth*[0 1 1 0],cornOriginy*[-1 -1 1 1],...
            -(cornOriginz-cornRad)*[1 1 1 1],greyRed(jj,:),'LineStyle',...
            'none','FaceLighting','gouraud','DiffuseStrength',1)
    end

    % Electron beam
    ebeamXrange = -2*L/3:4*L/3000:2*L/3;
    ebeamx = ones(1,size(ebeamXrange,2));
    ebeamAlp = jj*0.025/200;
    % Suppress e-beam oscillations outside of magnet array
    squareFn = 1./(1+exp(-25*(ebeamXrange+floor(L/2))))*...
        1./(1+exp(25*(ebeamXrange-floor(L/2))));
    % e-beam itself
    ebeamy = ebeamAlp*cos(2*pi*ebeamXrange/(2*moonWidth)).*squareFn;
    [ebX,ebY,ebZ] = cylinder(ebeamx);
    ebeam = surf(0.005*ebX,0.01*ebY+ebeamy',4*L*ebZ/3-2*L/3,'lineStyle',...
        'none','faceColor',myBlue,'FaceLighting','gouraud','DiffuseStrength',1);
    rotate(ebeam,[0 1 0],90,[0 0 0])

    % -------------------------------------------------------------------------
    % Create light cone
    conecoord1 = 0:0.01:pi;     % coordinates of light cone schematic
    conecoord2 = 0:0.02:2*pi;   % coordinates of light cone schematic
    [conecoord1,conecoord2]=meshgrid(conecoord1,conecoord2);

    zcoordcone1 = 0.2*(1 - sin(conecoord1)).*sin(conecoord1).*cos(conecoord2);
    ycoordcone1 = 0.2*(1 - sin(conecoord1)).*sin(conecoord1).*sin(conecoord2);
    xcoordcone1 = 8*abs(cos(conecoord1)).^2;
    coneAlpha = jj/500; % Light cone brightness
    lightCone = surf(xcoordcone1,ycoordcone1,zcoordcone1,...
        'FaceAlpha',coneAlpha,'FaceColor',myGold,'LineStyle','none','FaceLighting',...
        'gouraud','DiffuseStrength',1);


    % -------------------------------------------------------------------------
    % Plot undulator spectrum
    str2 = num2str(jj);
    filename = [str1 str2 str3];
    IDspect = importdata(filename,' ',2); % Ignore first two lines
    IDspectrum = IDspect.data;

    xdata = IDspectrum(:,1);
    ydata = 0.0*xdata(:,1);
    zdata = IDspectrum(:,2); % Bring zData minimum to zero
    zdata = log(zdata)/log(maxY);
    % Fit polynomial expression to baseline of the high-energy tail of the 
    % spectrum. This is done to avoid annoying jumps in the vertical 
    % position of thespectrum when a new peak drifts into the end of that 
    % spectrum. The final y-value of the polynomial fit is used as minY
    ii = 1;
    if (zdata(end) ~= min(zdata))
        % Discover local minima in second half of spectrum and make an array of them
        for kk = 40000:54995
            if (zdata(kk-6) > zdata(kk)) && (zdata(kk+6) > zdata(kk))
                minArray(ii,1) = xdata(kk);
                minArray(ii,2) = zdata(kk);
                ii = ii+1;
            end
        end
        % Polynomial fit of high-energy part of spectrum 
        p1 = polyfit(minArray(:,1),minArray(:,2),1); 
        fit1 = polyval(p1,xdata(32000:end));
        minY = fit1(end);
    else
        minY = min(zdata);
    end
    xdata = 4*IDspectrum(:,1)/max(IDspectrum(:,1))+0.8*L;
    z2data = 2*(zdata-minY)+0.05;
    scatter3(xdata,ydata,z2data,7,rgb,'filled');
    clear minArray 
    % -------------------------------------------------------------------------
    % Plotting rendering gumpf
    axis tight
    axis equal
    axis off
    light('Position',lp,'Style','infinite');
    light('Position',lp2,'Style','infinite');
    set(gca, 'Projection','perspective');
    set(gca,'View',[64,10]);
    frame = getframe(gcf);
    writeVideo(vid,frame);
    hold off
end % jj-loop for the frames
% Output the movie as an mpg file
close(vid)
