function [xn]=quasi
%
% MATLAB function from:
% "Elements of Modern X-ray Physics" by Jens Als-Nielsen and Des McMorrow
%
% Calculates: Positionso of atoms xn in a Fibonacci chain from the strip projection method,
%             and calculates the scattered intensity
% Calls to: pline, isinpoly, arrow

figure; axes('position',[0.15 0.15 0.8 0.8],'visible','off'); axis equal

tau=(1+sqrt(5))/2;         % golden mean
latp=sqrt(1+tau^2);        % lattice parameter of 2D lattice
Nx=10; Ny=10;              % number of lattice points
angle=atan(1/tau)*180/pi;  % angle of strip
Delta=1+tau;               % width of strip

%----- Draw strip and rotate
h=patch([0 Nx*latp/cos(angle*pi/180) Nx*latp/cos(angle*pi/180) 0],...
        [0 0 Delta Delta],[0.7 0.7 1]);
rotate(h,[0 0 1],angle,[0 0 0])
vp=get(h,'Vertices');

%----- Draw lattice
x=[]; for i=0:Nx; for j=0:Ny; x=[x;i j]; end; end
x=x*latp;
line(x(:,1),x(:,2),'linestyle','none','marker','o','markerfacecolor','g','markersize',6);

%----- Find lattice points that lie in strip
isp=isinpoly(x(:,1),x(:,2),vp(:,1),vp(:,2)); x(find(isp~=1),:)=[];

b=vp(2,1:2);                % end point of line xn
%----- Draw perpendicular lines from points in strip to xn
xn=[0]; yn=[0];
for ix=1:length(x)
    [intx,inty]=pline(b,x(ix,:));   
    xn=[xn; intx]; yn=[yn; inty];
end
line(x(:,1),x(:,2),'linestyle','none','marker','o','markerfacecolor','w','markersize',6);

%------ Label the graph
xnd=diff(xn);L=max(xnd); S=min(xnd);
for id=1:length(xnd)
if abs(xnd(id)-L)< 0.02 ,col=[0.6 0.6 0.6]; lab='L'; else col=[1 0 0]; lab='S'; end   
line([xn(id) xn(id+1)],[yn(id) yn(id+1)],'color',col,'linewidth',2.0)
text(0.5*(xn(id)+xn(id+1)),-2,lab,'color',col,'horizontalalignment','center',...
     'Fontsize',18,'FontName','Times')
end
arrow([-2 5*latp],[-2 6*latp],8,'ends','both')
text(-7,5.5*latp,'\surd(1+\tau^2)','Fontsize',24,'FontName','Times')
arrow(b,b*1.10,10)
text(b(1)*1.10,b(2)*1.0-0.2,'{\it x}_n','Fontsize',24,'FontName','Times',...
     'horizontalalignment','center')
h=arrow(-0.1*b,-0.1*b+[0 1+tau],10,'ends','both'); rotate(h,[0 0 1],angle,[-0.1*b 0])
text(-0.1*b(1),-0.1*b(1)+2.5,'\Delta','Fontsize',24,'FontName','Times')
circ=4.5*latp; ax=circ*cos(angle*pi/180):0.005:circ; ay=sqrt(circ^2-ax.^2); line(ax,ay); 

xn=xn./cos(angle*pi/180); 

%----- Calculate scttareing from chain

figure; axes('position',[0.15 0.15 0.8 0.8])
Q=[0:0.01:20];
F=sum(exp(sqrt(-1)*xn*Q));
plot(Q,F.*conj(F))
set(gca,'FontName','Times','FontSize',24,'Xtick',[0 5 10 15 20])
xlabel('Wavevector transfer (Å^{-1})'); ylabel('Intensity')

function [intx,inty]=pline(b,c)
%
% MATLAB function from:
% "Elements of Modern X-ray Physics" by Jens Als-Nielsen and Des McMorrow
%
% Calculates: Draws a perpendicular line from point c(x,y) 
%             to line that starts at origin and ends at point b(x,y)

if norm(c)==0
   tc=0;
elseif c(1)==0 & c(2)~=0
   tc=pi/2;
else
   tc=atan(c(2)./c(1));
end   
tb=atan(b(2)./b(1)); dt=tc-tb;

intx=norm(c)*cos(dt)*cos(tb); inty=norm(c)*cos(dt)*sin(tb);
line([c(1) intx],[c(2) inty],'color','w','linewidth',1.5)

function  isin = isinpoly(x,y,xp,yp)
% ISIN = ISINPOLY(X,Y,XP,YP)   Finds whether points with coordinates X and Y are inside 
%        or outside of a polygon with vertices XP, YP. Returns matrix ISIN of the same 
%        size as X and Y with 0 for points outside a polygon, 1 for inside points and  
%        0.5 for points belonging to a polygon XP, YP itself.
%  Copyright (c) 1995  by Kirill K. Pankratov
%	kirill@plume.mit.edu, 4/10/94, 8/26/94.

%----- Handle input
if nargin<4
  fprintf('\n  Error: not enough input arguments.\n\n')
  return
end
%----- Make the contour closed and get the sizes
xp = [xp(:); xp(1)]; yp = [yp(:); yp(1)];
sz = size(x); x = x(:); y = y(:);
lp = length(xp); l = length(x);
ep = ones(1,lp); e = ones(1,l);
%----- Calculate cumulative change in azimuth from points x,y to all vertices
A = diff(atan2(yp(:,e)-y(:,ep)',xp(:,e)-x(:,ep)'))/pi;
A = A+2*((A<-1)-(A>1));
isin = any(A==1)-any(A==-1);
isin = (abs(sum(A))-isin)/2;
%----- Check for boundary points
A = (yp(:,e)==y(:,ep)')&(xp(:,e)==x(:,ep)');
fnd = find(any(A));
isin(fnd) = .5*ones(size(fnd));
isin = round(isin*2)/2;
%----- Reshape output to the input size
isin = reshape(isin,sz(1),sz(2));
