function merry_christmas
% Merry Christmas
% vand@dtu.dk 2013

[X,Y,Z] = sphere(0.05); % sphere for decoration
dh = 0.75; % height of each segment
N = 4; % number of segments
R0 = 1.25; % radius of the lowest segment
RN = 1; % radius of the highest segment
D = 28; % number of decorations per segment

h = 0:dh:(N-1)*dh; % vector containing segment height position
r = linspace(R0,RN,length(h)); % vector containing segment radii
figure
%set(gcf,'units','normalized','position',[0 0 1 1]); % full screen
set(gca,'position',[0 0 1 1]);
hold on, axis vis3d off, campos([1 -2 1]*20), camlight('right')

P = zeros(N*D,3); % placeholder for decorations
for i = 1:length(h) % drawing the tree
    [Xs,Ys,Zs,p] = segment(r(i),h(i),D); % creating a segment
    surf(Xs,Ys,Zs,'EdgeColor','none','FaceColor',[0 0.9 0.3],...
        'FaceLighting','phong') % drawing a segment
    P((i-1)*D+1:i*D,:) = p;
end
zlim([0 dh*N+dh])
pause(0.5)

for j = 1:size(P,1) % adding the decoration
    surf(X+P(j,1),Y+P(j,2),Z+P(j,3),'EdgeColor','none','FaceColor','r',...
        'FaceLighting','phong')
    pause(0.05)
end
pause(0.5)

[xs,ys,zs] = star;
h = tubeplot(0.3*xs,0.3*ys,0.3*zs+N*dh+0.5,0.01); % adding the star
set(h,'EdgeColor','none','FaceColor','y','FaceLighting','phong')

function [X,Y,Z] = sphere(a)
% creates a sphere of radius a
t = linspace(0,pi,8); % 8 height points
[X,Y,Z] = rotate_profil(a*sin(t),a*cos(t),8); % 8 rotation angles

function [X,Y,Z,p] = segment(a,h,D)
% creates a segment of radius a, at the height h, with D decorations
x = linspace(0,1,32); % 32 height points
z = ((abs(x)-1.5).^2-1/4)/2;
[X,Y,Z] = rotate_profil([x,0]*a,[z,0]+h,32); % 32 rotation angles
i = randi(numel(X),[D,1]); % 24 decorations per segment
p = [X(i),Y(i),Z(i)];

function [X,Y,Z] = rotate_profil(x,z,n)
% rotates an (x,z) profile around z=0 axis

a = linspace(0,2*pi,n+1);
X = x(:)*sin(a);
Y = x(:)*cos(a);
Z = z(:)*ones(size(a));

function [x,y,z] = star
% starshape curve

a = 1;
b = 2/5*a;
t = linspace(0,4*pi);

x = (a-b)*sin(t) - b*sin((a/b-1)*t);
y = zeros(size(t));
z = (a-b)*cos(t) + b*cos((a/b-1)*t);

function h = tubeplot(x,y,z,r)
% plots a curve as a tube

[~,n,b]=frenet(x(:),y(:),z(:));

theta=linspace(0,2*pi,7);
rep = ones(1,7);

X = x(:)*rep + r*(n(:,1)*cos(theta) + b(:,2)*sin(theta));
Y = y(:)*rep + r*(n(:,2)*cos(theta) + b(:,2)*sin(theta));
Z = z(:)*rep + r*(n(:,3)*cos(theta) + b(:,3)*sin(theta));

h = surf(X,Y,Z);

function [t,n,b]=frenet(x,y,z)
% Frenet–Serret frame of a curve: tangent, normal and binormal vector

p=[x(:),y(:),z(:)];
t = normalize(gradient(p));
n = normalize(gradient(t));
b = normalize(cross(t,n,2));

function x = normalize(x)
m = (sum((x.^2),2)).^0.5;
n = m~=0;
x(n,:) = x(n,:)./(m(n)*[1 1 1]);