Cartoon of a generic beamline
Matlab codes
% Cartoon movie of a generic beamline
% Contains in order of moving downstream with the photons...
% ID source
% Blade BPM
% Beam-defining aperture
% High-pass filter window
% Beam-defining slits
% Collimating mirror
% Double-crystal monochromator
% Second focussing mirror
% Pinhole
% FZP
% I0 monitor
% Crystal sample
% Detector
clear; close all;
% Source to fixed elements of BL
s2BDA = 1.1; % Horizontal distance source to front face of BDA
s2HPF = 1.55; % Horizontal distance source to high-pass filter
s2BDS = 1.84; % Horizontal distance source to beam-defining slits
s2M1 = 2.35; % Horizontal distance source to centre of first mirror
s2X1 = 3.0; % Horizontal distance source to first crystal
s2M2 = 4.3; % Horizontal distance source to second mirror
s2FZP = 5.25; % Horizontal distance source to FZP
s2pin = 4.95; % Horizontal distance source to pinhole
s2I0 = 5.35; % Horizontal distance source to I0 monitor
s2sample = 5.65; % Horizontal distance source to sample
s2detector = 6.0; % Horizontal distance source to detector
% Other fixed geometrical parameters
zOffset = 0.44; % Vertical offset of source radiation and radiation after second mirror
m2ht = zOffset; % Redundant, but nice to keep naming convention consistent
mirTilt = 3; % Tilt of mirrors
x1ht = (s2X1-s2M1)*tand(2*mirTilt); % Height of first crystal
bb = m2ht - x1ht; % Difference in height between 2nd mirror and first crystal
aa = s2M2 - s2X1; % Difference in horizontal distance between 2nd mirror and 1st crystal
% Create rainbow rgb colour scheme
% dark red = (0.5 0 0)
% red = (1 0 0)
% orange = 1 0.5 0
% yellow = 1 1 0
% green = 0 1 0
% cyan = 0 0.5 1
% blue = 0 0 1
% purple = 0.5 0 1
% violet = 0.5 0 0.5
z45 = zeros(1,45); % 45 zeros
o45 = z45+1; % 45 ones
op545 = z45+0.5; % 45 lots of 0.5
o21 = 0:1/44:1; % 0 to 1 in 44 steps
o2op5 = 0:0.5/44:0.5; % 0 to 0.5 in 44 steps
op521 = o2op5+0.5; % 0.5 to 1 in 44 steps
rcomp = [op521 o45 o45 flip(o21) z45 z45 o2op5 op545];
gcomp = [z45 o2op5 op521 o45 flip(op521) flip(o2op5) z45 z45];
bcomp = [z45 z45 z45 z45 o21 o45 o45 flip(op521)];
% Set up colour scheme for ID spectrum
Emin = 3500;
Emax = 18500;
% Positions of odd-harmonic maxima (3, 5, 7, 9, 11)
oddMax = [5000 8330 11670 15000 18340];
Esteps = 1+(oddMax(5) - oddMax(1))/10; % # steps across entire used spectrum
subred = 0.1:0.4/149:0.5; % rcomp for 3500 to 5000 eV: energy under minimum used
supvio = 0.5:-0.25/15:0.25; % rcomp and bcomp for 18340 to 18500 eV: energy above maximum used
z16 = zeros(1,16); % 16 zeros
z150 = zeros(1,150); % 150 zeros
colLength = 1:Esteps;
colIndx=1:size(rcomp,2); % size of colour index used between 5000 and 18340 eV
temp = griddedInterpolant(colIndx',rcomp');
IDcolx = linspace(1,size(rcomp,2),Esteps);
% Take interpolated used values and embed them in concatenation with values
% below and above used energy range
IDcol(:,1) = [subred temp(IDcolx) supvio];
temp = griddedInterpolant(colIndx',gcomp');
IDcolx = linspace(1,size(gcomp,2),Esteps);
IDcol(:,2) = [z150 temp(IDcolx) z16];
temp = griddedInterpolant(colIndx',bcomp');
IDcolx = linspace(1,size(bcomp,2),Esteps);
IDcol(:,3) = [z150 temp(IDcolx) supvio];
vid = VideoWriter('beamline.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);
myBlue = [0.4 0.44 0.73];
myCopper = [0.73 0.45 0.20];
myGold = [0.7969 0.6406 0.2383];
midGrey = [0.5 0.5 0.5];
lightGrey = [0.8 0.8 0.8];
darkGrey = [0.1 0.1 0.1];
lp = [-1 0 1];
lp2 = [1 0.1 0.5];
M = dlmread('080.dat',' ', 2, 0); % Read in ID spectrum file, ignore first two lines
% Select E and flux between 3.5 and 18.5 keV: a bit below and above used
% harmonic values of oddMax
IDspect = M(M(:,1) >= 3500 & M(:,1) <= 18500,[1 3]);
ydata = 0.0*IDspect(:,1);
zdata = log10(IDspect(:,2))/16 - 0.3;
% What about xdata?? Chill dude! Coming up in the ii loop...
cycle = 0; % # of cycles of undulator gap
for ii = 1:360
hold off;
newplot;
light('Position',lp,'Style','infinite');
light('Position',lp2,'Style','infinite');
set(gca, 'Projection','perspective');
set(gca,'view',[-69,10]);
axis equal
axis off
axis tight
set(gca, 'Projection','perspective');
% -------------------------------------------------------------
% Undulator magnet array
magwd = 0.07; % magnet width x
magle = 0.16; % magnet length y
maght = 0.05; % magnet height z
for jj = -0.4:0.2:0.4
% Increase in gap from frame to frame. the mod function resets the
% gap size to the minimum every 90 frames. This also determines the
% expansion and reset of the undulator spectrum plotted above the
% ID.
gInc = (0.03/90)*mod(ii-1,90);
plotcube([magwd magle maght],[jj -magle/2 0.03+gInc],1,midGrey,0);
hold on
plotcube([magwd magle maght],[jj+0.1 -magle/2 0.03+gInc],1,lightGrey,0);
plotcube([magwd magle maght],[jj -magle/2 -0.03-maght-gInc],1,lightGrey,0);
plotcube([magwd magle maght],[jj+0.1 -magle/2 -0.03-maght-gInc],1,midGrey,0);
end
% Dummy vertical line to stop image resizing
plot3([-0.55 -0.55],[magle+0.01 magle+0.01],[-0.111 0.7],...
'color','w', 'LineWidth', 0.1);
% Plot expanding ID spectrum
expandIndx = mod(ii-1,90);
if (expandIndx == 0)
cycle = cycle+1;
end
% Expansion factor for each step in ii and within an expansion cycle
% of 90 steps
Efact = (oddMax(cycle+1)/oddMax(cycle))^(expandIndx/90);
IDspectExp(:,1) = IDspect(:,1)*Efact; % Energy values as undulator gap opens up
IDspectExp(:,2) = IDspect(:,2);
% As spectrum expands, select only that part between 3500 and 18500 eV
IDspectPlot = IDspectExp(IDspectExp(:,1) >= 3500 & IDspectExp(:,1) <= 18500,[1 2]);
xdata = -0.7+6.4e-5*(IDspectPlot(:,1));
ydata = 0.0*IDspectPlot(:,1);
zdata = log10(IDspectPlot(:,2))/16 - 0.34;
rgbcolInd = 1+floor((IDspectPlot(:,1)-3500)/10);
scatter3(xdata,ydata,zdata,4,IDcol(rgbcolInd,:),'filled');
% Draw arrow
energyNow = oddMax(cycle)*Efact;
arrowcolInd = 1+floor((energyNow-3500)/10);
arrowX = -0.7+6.4e-5*energyNow;
mArrow3([arrowX 0 0.69-0.07*ii/360],[arrowX 0 0.60-0.07*ii/360],...
'color',IDcol(arrowcolInd,:),'stemWidth',0.0035,...
'tipWidth',0.01,'facealpha',0.7);
% -------------------------------------------------------------
% Create a cone of radiation emanating from an undulator
t = 0:3;
r = t/100;
[x,y,z] = cylinder(r,50);
for jj = 0:0.1:1
sourceRad = surf(x/2,(2-jj)*y,s2BDA*z+jj/5,'FaceAlpha',0.08,'FaceColor',...
myGold,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(sourceRad,[0 1 0],90,[0,0,0]);
end
% -------------------------------------------------------------
% Blade BPM
plotcube([0.01 0.05 0.08],[0.8 -0.025 0.04],1,darkGrey,0);
plotcube([0.01 0.05 0.08],[0.8 -0.025 -0.12],1,darkGrey,0);
plotcube([0.01 0.08 0.05],[0.8 -0.12 -0.025],1,darkGrey,0);
plotcube([0.01 0.08 0.05],[0.8 0.04 -0.025],1,darkGrey,0);
% -------------------------------------------------------------
% Beam-defining aperture
rBDA = 0.08;
xBDA = 0.05;
yBDA = 0.02;
BDAlength = 0.2;
rPipe = 0.01;
lPipe = 0.08;
th = 0:pi/60:2*pi;
x = rBDA*sin(th);
y = rBDA*cos(th);
% x-coords for end face with rectangular hole via concatenation
faceBDAx = [x 0 xBDA/2 xBDA/2 -xBDA/2 -xBDA/2 0 0];
% y-coords for end face with rectangular hole via concatenation
faceBDAy = [y yBDA/2 yBDA/2 -yBDA/2 -yBDA/2 yBDA/2 yBDA/2 rBDA];
faceBDAz = 0.0*faceBDAx;
faceBDA1 = fill3(faceBDAz+s2BDA,faceBDAx,faceBDAy,myCopper,'LineStyle','none',...
'FaceLighting','gouraud','DiffuseStrength',1);
faceBDA2 = fill3(faceBDAz+s2BDA+BDAlength,faceBDAx,faceBDAy,myCopper,'LineStyle','none',...
'FaceLighting','gouraud','DiffuseStrength',1);
[x,y,z] = cylinder(rBDA,50);
BDAcyl = surf(x,y,s2BDA+z*BDAlength,'FaceAlpha',1,'FaceColor',myCopper,'LineStyle',...
'none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(BDAcyl,[0 1 0],90,[0,0,0]);
[x,y,z] = cylinder(rPipe,50);
pipe1 = surf(x+s2BDA+BDAlength/2-BDAlength/4,y,rBDA+z*lPipe,'FaceAlpha',1,...
'FaceColor',myCopper,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(pipe1,[1 0 0],180,[s2BDA+BDAlength/2-BDAlength/4,0,0]);
pipe2 = surf(x+s2BDA+BDAlength/2+BDAlength/4,y,rBDA+z*lPipe,'FaceAlpha',1,...
'FaceColor',myCopper,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(pipe2,[1 0 0],180,[s2BDA+BDAlength/2+BDAlength/4,0,0]);
% -------------------------------------------------------------
% High-pass filter window and continued beam
HPFwi = 0.14; % Breadth of high-pass filter
HPFhe = 0.10; % Height of high-pass filter
HPFfr = 0.02; % Width of high-pass filter frame
HPFth = 0.01; % Thickness of high-pass filter frame
plotcube([HPFth HPFwi HPFfr],[s2HPF -HPFwi/2 HPFhe/2-HPFfr],1,midGrey,0); %top edge
plotcube([HPFth HPFwi HPFfr],[s2HPF -HPFwi/2 -HPFhe/2],1,midGrey,0); % bottom edge
plotcube([HPFth HPFfr HPFhe],[s2HPF -HPFwi/2 -HPFhe/2],1,midGrey,0); %left side
plotcube([HPFth HPFfr HPFhe],[s2HPF HPFwi/2-HPFfr -HPFhe/2],1,midGrey,0); %left side
plotcube([HPFth/10 HPFwi HPFhe],[s2HPF+HPFth/2 -HPFwi/2 -HPFhe/2],0.2,[0 0 1],0); %filter
% Create continued cone of radiation from an undulator
t = s2BDA+BDAlength:0.01:s2BDS;
coneLength = s2BDS - s2BDA - BDAlength;
r = t/100;
[x,y,z] = cylinder(r,50);
contRad = surf(x/2,2*y,s2BDA+BDAlength+(coneLength)*z,'FaceAlpha',0.4,...
'FaceColor',myGold,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(contRad,[0 1 0],90,[0,0,0]);
% -------------------------------------------------------------
% Beam-defining slits
BDSth = 0.01;
BDSwi = 0.05;
BDShe = 0.05;
plotcube([BDSth BDSwi BDShe],[s2BDS+BDSth/4 -1.4*BDSwi -0.5*BDShe],1,midGrey,0);
plotcube([BDSth BDSwi BDShe],[s2BDS+BDSth/4 0.4*BDSwi -0.5*BDShe],1,midGrey,0);
plotcube([BDSth BDSwi BDShe],[s2BDS-BDSth -0.5*BDSwi -1.1*BDShe],1,midGrey,0);
plotcube([BDSth BDSwi BDShe],[s2BDS-BDSth -0.5*BDSwi 0.1*BDShe],1,midGrey,0);
% Create continued cone of radiation from an undulator
t = s2BDS:0.01:s2M1*1.07;
coneLength = 1.07*s2M1 - s2BDS;
r = t/100;
[x,y,z] = cylinder(r,50);
contRad2 = surf(x/2,1.25*y,s2BDS+(coneLength)*z,'FaceAlpha',0.4,...
'FaceColor',myGold,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(contRad2,[0 1 0],90,[0,0,0]);
% -------------------------------------------------------------
% Collimating mirror
aminor = 1; % Torus minor radius
Rmajor = 70.; % Torus major radius
% Curved mirror surface
theta = linspace(-pi/36, pi/36, 256) ; % Poloidal angle
phi = linspace(-pi/720, pi/720, 256) ; % Toroidal angle
[p, t] = meshgrid(phi, theta);
x = (Rmajor + aminor.*cos(t)) .* cos(p);
z = (Rmajor + aminor.*cos(t)) .* sin(p);
y = aminor.*sin(t);
mirrorSurf = surf(10*(x-max(max((x))))+s2M1, y, z,'FaceAlpha',1,'FaceColor',...
lightGrey,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',0.05);
rotate(mirrorSurf,[0 1 0],90-mirTilt,[s2M1,0,0]);
% Side walls
curvedx = 10*((Rmajor+aminor*cos(pi/36))*cos(phi)-max(max((x))))+s2M1;
curvedx0 = 10*((Rmajor+aminor*cos(0/36))*cos(phi)-max(max((x))))+s2M1;
curvedSidex = [curvedx max(curvedx0)+0.1 max(curvedx0)+0.1 curvedx(1)];
mirL = max((Rmajor + aminor*cos(pi/36))*sin(phi));
curvedz = (Rmajor + aminor*cos(pi/36))*sin(phi);
curvedSidez = [curvedz max(curvedz) min(curvedz) min(curvedz)];
curvedSidey = 0.0*curvedSidex+aminor*sin(pi/36);
curvedSide1 = fill3(curvedSidex,curvedSidey,curvedSidez,midGrey,'LineStyle','none',...
'FaceLighting','gouraud','DiffuseStrength',1);
curvedSide2 = fill3(curvedSidex,curvedSidey-2*aminor*sin(pi/36),curvedSidez,...
midGrey,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(curvedSide1,[0 1 0],90-mirTilt,[s2M1,aminor*sin(pi/36),0]);
rotate(curvedSide2,[0 1 0],90-mirTilt,[s2M1,-aminor*sin(pi/36),0]);
% End faces
endy = aminor*sin(theta);
endFacey = [endy max(endy) min(endy) min(endy)];
endx = 10*((Rmajor+aminor*cos(theta))*cos(pi/720)-max(max((x))))+s2M1;
endFacex = [endx max(curvedx0)+0.1 max(curvedx0)+0.1 endx(1)];
endFacez = 0.0*endFacex+max(curvedz);
endFace1 = fill3(endFacex,endFacey,endFacez,midGrey,'LineStyle','none',...
'FaceLighting','gouraud','DiffuseStrength',1);
rotate(endFace1,[0 1 0],90-mirTilt,[s2M1,0,0]);
endFace2 = fill3(endFacex,endFacey,endFacez-2*max(curvedz),midGrey,'LineStyle','none',...
'FaceLighting','gouraud','DiffuseStrength',1);
rotate(endFace2,[0 1 0],90-mirTilt,[s2M1,0,0]);
bottomFacex = [s2M1+0.1 s2M1+0.1 s2M1+0.1 s2M1+0.1];
bottomFacey = [max(endy) max(endy) min(endy) min(endy)];
bottomFacez = [-mirL +mirL +mirL -mirL];
bottomFace = fill3(bottomFacex,bottomFacey,bottomFacez,midGrey,'LineStyle','none',...
'FaceLighting','gouraud','DiffuseStrength',1);
rotate(bottomFace,[0 1 0],90-mirTilt,[s2M1,0,0]);
% -------------------------------------------------------------
% DCM
% Sizes of crystal sides
Xx = 0.125;
Xy = 0.125;
Xz = 0.04;
% Create continued cone of parallel radiation from an undulator
parBeamLength = s2X1-0.93*s2M1;
r = 2.5/100;
[x,y,z] = cylinder(r,50);
contRad3 = surf(x/2+s2M1,1.25*y,1.04*parBeamLength*z-0.07*s2M1,...
'FaceAlpha',0.4,'FaceColor',myGold,'LineStyle',...
'none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(contRad3,[0 1 0],90-2*mirTilt,[s2M1,0,0]);
% Create 1st crystal
thBragg = 30-ii/18;
% Calculate position of second crystal
cc = tand(2*mirTilt);
dd = tand(2*(mirTilt+thBragg));
xx = (cc*aa-bb)/(cc-dd);
s2X2 = s2X1 + xx;
yy = dd*xx;
x2ht = x1ht + yy;
% top and bottom first DCM crystal faces
xFacex = [-0.5 0.5 0.5 -0.5];
xFacey = [-0.5 -0.5 0.5 0.5];
xFacez = [0 0 0 0];
x1Face1 = fill3(Xx*xFacex+s2X1,Xy*xFacey,Xz*xFacez+x1ht,myBlue,'LineStyle',...
'none','FaceLighting','gouraud','DiffuseStrength',0.1);
rotate(x1Face1,[0 1 0],-thBragg-2*mirTilt,[s2X1,0,x1ht]);
x1Face2 = fill3(Xx*xFacex+s2X1,Xy*xFacey,Xz*xFacez+x1ht-Xz,myBlue,'LineStyle',...
'none','FaceLighting','gouraud','DiffuseStrength',0.1);
rotate(x1Face2,[0 1 0],-thBragg-2*mirTilt,[s2X1,0,x1ht]);
% side faces
xSidex = [-0.5 0.5 0.5 -0.5];
xSidey = [0.0 0.0 0.0 0.0];
xSidez = [0.5 0.5 -0.5 -0.5];
x1Side1 = fill3(Xx*xSidex+s2X1,Xy*xSidey+Xy/2,Xz*xSidez-Xz/2+x1ht,myBlue,...
'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(x1Side1,[0 1 0],-thBragg-2*mirTilt,[s2X1,0,x1ht]);
x1Side2 = fill3(Xx*xSidex+s2X1,Xy*xSidey-Xy/2,Xz*xSidez-Xz/2+x1ht,myBlue,...
'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(x1Side2,[0 1 0],-thBragg-2*mirTilt,[s2X1,0,x1ht]);
% end faces
xEndx = [0 0 0 0];
xEndy = [-0.5 0.5 0.5 -0.5];
xEndz = [-0.5 -0.5 0.5 0.5];
x1End1 = fill3(Xx*xEndx+s2X1+Xx/2,Xy*xEndy,Xz*xEndz-Xz/2+x1ht,myBlue,...
'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(x1End1,[0 1 0],-thBragg-2*mirTilt,[s2X1,0,x1ht]);
x1End2 = fill3(Xx*xEndx+s2X1-Xx/2,Xy*xEndy,Xz*xEndz-Xz/2+x1ht,myBlue,...
'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(x1End2,[0 1 0],-thBragg-2*mirTilt,[s2X1,0,x1ht]);
% Create 1st rotation stage
gonioHt = 0.08;
r = 0.125;
tt = 0:pi/45:2*pi;
iir = ii*pi/180;
xgon = r*cos(tt-iir/18)+s2X1;
zgon = r*sin(tt-iir/18)+x1ht;
ygon = 0.0*xgon+Xy/2;
xinn = 0.88*r*cos(tt-iir/18)+s2X1;
zinn = 0.88*r*sin(tt-iir/18)+x1ht;
yinn = 0.0*xgon+Xy/2-0.0001;
fill3(xgon,ygon,zgon,lightGrey,'LineStyle','none',...
'FaceLighting','gouraud','DiffuseStrength',1);
for tt = 1:2:89
plot3([xgon(tt) xinn(tt)],[ygon(tt) yinn(tt)],[zgon(tt) zinn(tt)],...
'color','w', 'LineWidth', 1);
end
[x,y,z] = cylinder(r,50);
gonioSides = surf(x+s2X1,y+Xy/2,gonioHt*z+x1ht,'FaceAlpha',1,'FaceColor',...
midGrey,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(gonioSides,[1 0 0],-90,[s2X1,Xy/2,x1ht]);
% Create continued cone of diffracted parallel radiation from an undulator
XVoffset = yy;
parBeamLength = XVoffset/sind(2*thBragg+2*mirTilt); % Length of beam between crystals
xOffset = parBeamLength*cosd(2*thBragg+2*mirTilt); % Horizontal offset of DCM crystals
r = 2.5/100;
[x,y,z] = cylinder(r,50);
contRad4 = surf(x/2+s2X1,1.25*y,1.1*parBeamLength*z-0.05*parBeamLength...
+x1ht,'FaceAlpha',0.4,'FaceColor',IDcol(arrowcolInd,:),'LineStyle',...
'none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(contRad4,[0 1 0],90-2*thBragg-2*mirTilt,[s2X1,0,x1ht]);
% Create 2nd crystal
% top and bottom faces
x2Face1 = fill3(Xx*xFacex+s2X2,Xy*xFacey,Xz*xFacez+XVoffset+x1ht,myBlue,...
'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',0.1);
rotate(x2Face1,[0 1 0],-thBragg-2*mirTilt,[s2X2,0,XVoffset+x1ht]);
x2Face2 = fill3(Xx*xFacex+s2X2,Xy*xFacey,Xz*xFacez+XVoffset+x1ht+Xz,myBlue,...
'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',0.1);
rotate(x2Face2,[0 1 0],-thBragg-2*mirTilt,[s2X2,0,XVoffset+x1ht]);
% side faces
x2Side1 = fill3(Xx*xSidex+s2X2,Xy*xSidey+Xy/2,Xz*xSidez+Xz/2+XVoffset+x1ht,...
myBlue,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(x2Side1,[0 1 0],-thBragg-2*mirTilt,[s2X2,0,XVoffset+x1ht]);
x2Side2 = fill3(Xx*xSidex+s2X2,Xy*xSidey-Xy/2,Xz*xSidez+Xz/2+XVoffset+x1ht,...
myBlue,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(x2Side2,[0 1 0],-thBragg-2*mirTilt,[s2X2,0,XVoffset+x1ht]);
% end faces
x2End1 = fill3(Xx*xEndx+s2X2+Xx/2,Xy*xEndy,XVoffset+Xz*xEndz+Xz/2+x1ht,...
myBlue,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(x2End1,[0 1 0],-thBragg-2*mirTilt,[s2X2,0,XVoffset+x1ht]);
x2End2 = fill3(Xx*xEndx+s2X2-Xx/2,Xy*xEndy,XVoffset+Xz*xEndz+Xz/2+x1ht,...
myBlue,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(x2End2,[0 1 0],-thBragg-2*mirTilt,[s2X2,0,XVoffset+x1ht]);
% Create 2nd rotation stage
gonioHt = 0.08;
r = 0.125;
tt = 0:pi/45:2*pi;
iir = ii*pi/180;
xgon = r*cos(tt-iir/18)+s2X2;
zgon = r*sin(tt-iir/18)+x2ht;
ygon = 0.0*xgon+Xy/2;
xinn = 0.88*r*cos(tt-iir/18)+s2X2;
zinn = 0.88*r*sin(tt-iir/18)+x2ht;
yinn = 0.0*xgon+Xy/2-0.0001;
fill3(xgon,ygon,zgon,lightGrey,'LineStyle','none',...
'FaceLighting','gouraud','DiffuseStrength',1);
for tt = 1:2:89
plot3([xgon(tt) xinn(tt)],[ygon(tt) yinn(tt)],[zgon(tt) zinn(tt)],...
'color','w','LineWidth', 1);
end
[x,y,z] = cylinder(r,50);
gonioSides2 = surf(x+s2X2,y+Xy/2,gonioHt*z+x2ht,'FaceAlpha',1,'FaceColor',...
midGrey,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(gonioSides2,[1 0 0],-90,[s2X2,Xy/2,x2ht]);
% -------------------------------------------------------------
% Second focussing mirror
% Radiation from X2 to M2
parBeamLength = s2M2-s2X2;
m2ht = x2ht+parBeamLength*tand(2*mirTilt);
r = 2.5/100;
[x,y,z] = cylinder(r,50);
contRad5 = surf(x/2+s2X2,1.25*y,1.25*parBeamLength*z-0.05*parBeamLength...
+x2ht,'FaceAlpha',0.4,'FaceColor',IDcol(arrowcolInd,:),'LineStyle',...
'none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(contRad5,[0 1 0],90-2*mirTilt,[s2X2,0,x2ht]);
% Curved mirror surface
theta = linspace(-pi/36, pi/36, 256) ; % Poloidal angle
phi = linspace(-pi/720, pi/720, 256) ; % Toroidal angle
[p, t] = meshgrid(phi, theta);
x = (Rmajor + aminor.*cos(t)) .* cos(p);
z = (Rmajor + aminor.*cos(t)) .* sin(p);
y = aminor.*sin(t);
mirrorSurf = surf(10*(x-max(max((x))))+s2M2, y, z+m2ht,'FaceAlpha',1,'FaceColor',...
lightGrey,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(mirrorSurf,[0 1 0],270-mirTilt,[s2M2,0,m2ht]);
% Side walls
curvedx = 10*((Rmajor+aminor*cos(pi/36))*cos(phi)-max(max((x))))+s2M2;
curvedx0 = 10*((Rmajor+aminor*cos(0/36))*cos(phi)-max(max((x))))+s2M2;
curvedSidex = [curvedx max(curvedx0)+0.1 max(curvedx0)+0.1 curvedx(1)];
mirL = max((Rmajor + aminor*cos(pi/36))*sin(phi));
curvedz = (Rmajor + aminor*cos(pi/36))*sin(phi)+m2ht;
curvedSidez = [curvedz max(curvedz) min(curvedz) min(curvedz)];
curvedSidey = 0.0*curvedSidex+aminor*sin(pi/36);
curvedSide1 = fill3(curvedSidex,curvedSidey,curvedSidez,midGrey,'LineStyle','none',...
'FaceLighting','gouraud','DiffuseStrength',1);
curvedSide2 = fill3(curvedSidex,curvedSidey-2*aminor*sin(pi/36),curvedSidez,...
midGrey,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(curvedSide1,[0 1 0],270-mirTilt,[s2M2,aminor*sin(pi/36),m2ht]);
rotate(curvedSide2,[0 1 0],270-mirTilt,[s2M2,-aminor*sin(pi/36),m2ht]);
% End faces
endy = aminor*sin(theta);
endFacey = [endy max(endy) min(endy) min(endy)];
endx = 10*((Rmajor+aminor*cos(theta))*cos(pi/720)-max(max((x))))+s2M2;
endFacex = [endx max(curvedx0)+0.1 max(curvedx0)+0.1 endx(1)];
endFacez = 0.0*endFacex;
endFace1 = fill3(endFacex,endFacey,endFacez+m2ht+mirL,midGrey,'LineStyle','none',...
'FaceLighting','gouraud','DiffuseStrength',1);
rotate(endFace1,[0 1 0],270-mirTilt,[s2M2,0,m2ht]);
endFace2 = fill3(endFacex,endFacey,endFacez+m2ht-mirL,midGrey,'LineStyle','none',...
'FaceLighting','gouraud','DiffuseStrength',1);
rotate(endFace2,[0 1 0],270-mirTilt,[s2M2,0,m2ht]);
bottomFacex = [s2M2+0.1 s2M2+0.1 s2M2+0.1 s2M2+0.1];
bottomFacey = [max(endy) max(endy) min(endy) min(endy)];
bottomFacez = [m2ht-mirL m2ht+mirL m2ht+mirL m2ht-mirL];
bottomFace = fill3(bottomFacex,bottomFacey,bottomFacez,midGrey,'LineStyle','none',...
'FaceLighting','gouraud','DiffuseStrength',1);
rotate(bottomFace,[0 1 0],270-mirTilt,[s2M2,0,m2ht]);
% -------------------------------------------------------------
% Pinhole
% Create a cone of radiation focussing through a pinhole
beamLength = s2FZP-s2M2; % Horizontal distance M2 to FZP
t = s2M2:0.05:s2FZP;
r = (2.5/100)*((s2pin-t)/(s2pin-s2M2));
[x,y,z] = cylinder(r,50);
contRad6 = surf(1.5*x+s2M2,1.5*y,m2ht+1.05*beamLength*z-0.05*beamLength,...
'FaceAlpha',0.4,'FaceColor',IDcol(arrowcolInd,:),'LineStyle',...
'none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(contRad6,[0 1 0],90,[s2M2,0,m2ht]);
% Now create pinhole
pinRad = 0.005;
pinSize = 0.08;
pinFramey = [pinSize/2 -pinSize/2 -pinSize/2 pinSize/2 pinSize/2];
pinFramex = [0 0 0 0 0];
pinFramez = [pinSize/2 pinSize/2 -pinSize/2 -pinSize/2 pinSize/2];
pintheta = 0:pi/20:2*pi;
pinholey = pinRad*cos(pintheta);
pinholez = pinRad*sin(pintheta);
pinholex = 0.0*pinholey;
pinx = [pinFramex pinholex]+s2pin;
piny = [pinFramey pinholey];
pinz = [pinFramez pinholez]+m2ht;
fill3(pinx,piny,pinz,midGrey,'LineStyle','none',...
'FaceLighting','gouraud','DiffuseStrength',1);
fill3(pinx+0.005,piny,pinz,darkGrey,'LineStyle','none',...
'FaceLighting','gouraud','DiffuseStrength',1);
% -------------------------------------------------------------
% FZP
FZPth = 0:pi/100:2*pi;
imax = 11;
for jj = imax:-1:0
hold on
r = 0.014*(jj)^0.5;
z_circle = r*cos(FZPth)+m2ht;
y_circle = r*sin(FZPth);
x_circle = 0.0*y_circle+s2FZP+jj/10000;
if (round(jj/2) == jj/2)
FZPcol = [1 1 1];
else
FZPcol = myBlue;
end
zone = fill3(x_circle, y_circle, z_circle, FZPcol,...
'FaceAlpha','1','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);
end
% -------------------------------------------------------------
% BSB - not included for time being - it only confuses novices
% But do your worst...
% -------------------------------------------------------------
% I0
% Top plate
plotcube([0.04 0.04 0.005],[s2I0-0.02 -0.02 m2ht+0.02],1,lightGrey,0);
% Bottom plate
plotcube([0.04 0.04 0.005],[s2I0-0.02 -0.02 m2ht-0.025],1,lightGrey,0);
% HV connectors
plot3([s2I0 s2I0],[0 0],[m2ht+0.02 m2ht+0.07],'color','k', 'LineWidth', 3);
plot3([s2I0 s2I0],[0 0],[m2ht-0.02 m2ht-0.07],'color','k', 'LineWidth', 3);
% -------------------------------------------------------------
% Crystal sample
% Create a cone of radiation focussing on sample
beamLength = s2sample-s2FZP; % Horizontal distance M2 to FZP
t = 0:0.05:beamLength;
r = (s2FZP-s2pin)/(s2pin-s2M2)*(2.5/100)*((beamLength-t)/beamLength);
[x,y,z] = cylinder(r,50);
contRad7 = surf(1.5*x+s2FZP,1.5*y,m2ht+beamLength*z,...
'FaceAlpha',0.4,'FaceColor',IDcol(arrowcolInd,:),'LineStyle',...
'none','FaceLighting','gouraud','DiffuseStrength',1);
rotate(contRad7,[0 1 0],90,[s2FZP,0,m2ht]);
[V,F] = platonic_solid(3,0.04);
V(:,1) = V(:,1)+s2sample;
V(:,3) = V(:,3)+m2ht;
ps = patch('Faces',F,'Vertices',V,'FaceColor',myBlue,'FaceAlpha',0.8, ...
'EdgeColor','none','FaceLighting','gouraud','DiffuseStrength',1,'SpecularStrength',0.7);
direction = [0 1 0];
rotate(ps,direction,ii,[s2sample 0 m2ht]);
% -------------------------------------------------------------
% Detector
detSize = 0.34;
sensSize = 0.3;
plotcube([detSize detSize detSize],[s2detector+0.02 -detSize/2 m2ht-detSize/2],1,darkGrey,0);
plotcube([0.1 sensSize sensSize],[s2detector+0.0001 -sensSize/2 m2ht-sensSize/2],1,midGrey,0);
% Create rotating diffraction pattern
numBP = 10;
phi = ii;
rBP = 0.005;
kV = 1; % Radius of Ewald sphere
for hh = -0.5:0.5/numBP:0.5
for kk = -0.5:0.5/numBP:0.5
for ll = -0.5:0.5/numBP:0.5
hkl = (hh^2 + kk^2 + ll^2)^0.5;
if (hkl <= 0.7) % Plot out Bragg peaks within radius of 0.7 r.l.u.
hnew = hh*cosd(phi) + ll*sind(phi);
knew = kk; % in y-direction
lnew = ll*cosd(phi) - hh*sind(phi);
% Calculate Delk, absolute distance from centre of Ewald sphere to hkl Bragg peak
Delk = ((kV + hnew)^2 + (knew)^2 + (lnew)^2)^0.5;
Del = abs(kV - Delk); % Difference between centre of Bragg peak and Ewald sphere surface
if (Del < 2.0*rBP) % Detected Bragg peaks
% Plot Bragg peaks detected on detector
apos = s2detector;
bpos = 0.2*knew*((1+s2detector-s2sample)/(1+hnew));
cpos = m2ht+0.2*lnew*((1+s2detector-s2sample)/(1+hnew));
if (abs(bpos) < sensSize/2) && (abs(abs(cpos)-m2ht) < sensSize/2)
BP = scatter3(apos,bpos,cpos,4,'y','filled');
end % End of finding out if BP within detector
end
end
end % ll loop
end % kk loop
end % hh loop
% Store the frame
set(gca,'view',[-59,10]);
axis tight
frame = getframe(gcf);
writeVideo(vid,frame);
hold off
end
% Output the movie as an mpg file
close(vid);