%
% file:      	shells.m, (c) Matthew Roughan, Fri Aug 20 2010
% created: 	Fri Aug 20 2010 
% author:  	Matthew Roughan 
% email:   	matthew.roughan@adelaide.edu.au
% 
% Set the parameters for a shell, and draw it:
%
% Note that the shell follows a log-spiral, which defines the centre of mass curve by 
%    x = a*exp(|alpha|*phi).*cos(-sign(alpha)*phi);
%    y = a*exp(|alpha|*phi).*sin(-sign(alpha)*phi);
%    z = -h*exp(|alpha|*phi);
% i.e., the centre of mass is a log-spiral on a cone, oriented along the z-axis
%
% The shell surface is obtained by moving an ellipse along this path
%
% Parameters:
%
% log-spiral parameters
%    alpha          rate of decrease of spiral sections (default -0.06)
%                       negative for anticlockwise spiral (dextral)
% 		        positive for clockwise spiral (sinistral)
%    width          horizontal width at base of cone (default = 1)
%                       a = width/2
%    h              vertical height    (default = 1)
%    N              number of spirals  (default is such that largest spiral > 100x smallest)
%
% opening parameters
%    r_1            axis of ellipse in z-direction           (default = 0.5)
%    r_2            axis of ellipse along radii from centre  (default = 0.5)
%    beta           the ellipse is rotated at this angle with respect to vertical
%
%       this could be more general -- we could just specify a set of points
%
% decorations
%    f              vector giving frequencies for periodic growth rings
%    A              vector giving amplitudes (e.g. obtain with Fourier transform)
%   
%
% surface parameters
%    N_phi = number of steps, per revolution of spiral in surface
%    N_theta = number of steps, around the ellipse in the surface
%

% log-spiral parameters 
%    x = a*exp(|alpha|*phi).*cos(-sign(alpha)*phi);
%    y = a*exp(|alpha|*phi).*sin(-sign(alpha)*phi);
%    z = -h*exp(alpha*phi);
if (~exist('alpha','var'))
  alpha = -0.06;
end
if (alpha < 0)
  sign_alpha = -1;
  if (~exist('N','var'))
    N = ceil( (1/alpha) * log(1/100)/(2*pi) );
  end
elseif (alpha > 0)
  alpha = -alpha;
  sign_alpha = 1;
  if (~exist('N','var'))
    N = ceil( (1/alpha) * log(1/100)/(2*pi) );
  end
else
  alpha = 0;
  sign_alpha = 1;
  if (~exist('N','var'))
    N = 10; % spiral doesn't change size, so just choose a value
  end
end
if (~exist('width','var'))
  width = 1;
end
a = abs(width)/2;
if (~exist('h','var'))
  h = 1;
end
  

% opening parameters
%    r_1            axis of ellipse in z-direction           (default = 0.5)
%    r_2            axis of ellipse along radii from centre  (default = 0.5)
%    beta           the ellipse is rotated at this angle with respect to vertical
if (~exist('r1','var'))
  r1 = 0.5;
end
if (~exist('r2','var'))
  r2 = 0.5;
end
if (~exist('beta','var'))
  beta = 0;
end

% decoration parameters
if (~exist('f','var') | ~exist('A','var'))
  f = 0;
  A = 0;
end

% surface parameters
%    N_phi = number of steps, per revolution of spiral in surface
%    N_theta = number of steps, around the ellipse in the surface
if (~exist('N_phi','var'))
  N_phi = 40;
end
if (~exist('N_theta','var'))
  N_theta = 40;
end

%%%%
%%%%  Now draw the shell and its bits
%%%%

% center of mass of spiral
step = pi/15;
phi = 0:step:N*2*pi;
x = a*exp(alpha*phi).*cos(-sign_alpha*phi);
y = a*exp(alpha*phi).*sin(-sign_alpha*phi);
z = -h*exp(alpha*phi);
 
% cone
step = a/10;
r_cone = 0:step:a;
theta_cone = 0:pi/50:2*pi;
X_cone = r_cone'*cos(theta_cone);
Y_cone = r_cone'*sin(theta_cone);
Z_cone = -h * sqrt(X_cone.^2 + Y_cone.^2) / a - h/50;

figure(1)
hold off
plot3(x,y,z, 'k', 'linewidth', 3)
hold on
mesh(X_cone, Y_cone, Z_cone);
axis equal
title('log-spiral centre of mass of shell');


% wireframe of surface
step_phi = 2*pi/N_phi;
step_theta = 2*pi/N_theta;
phi = 0:step_phi:N*2*pi;
theta = pi + (0:step_theta:2*pi);
[Theta, Phi] = meshgrid(theta, phi);
x = a*exp(alpha*phi).*cos(phi);
y = a*exp(alpha*phi).*sin(phi);
z = -h*exp(alpha*phi);
Xd = r1*cos(Theta);
Yd = r2*sin(Theta);
R = exp(alpha*Phi);
X = a*R.*cos(Phi);
Y = a*R.*sin(Phi);
Z = -h*R;
surf_X = X + cos(Phi) .* Xd .* R + Yd .* R .* sin(Phi) * sin(beta);
surf_Y = Y + sin(Phi) .* Xd .* R - Yd .* R .* cos(Phi) * sin(beta);
surf_Z = Z + Yd .* R * cos(beta);
figure(2)
hold off
g = mesh(surf_X, surf_Y, surf_Z,'EdgeColor','k');
axis equal vis3d
hold on
plot3(X(1,:) + Xd(1,:), Y(1,:) - Yd(1,:)*sin(beta), Z(1,:) + Yd(1,:)*cos(beta), 'linewidth', 3);
title('shell surface wireframe');

% a picture of the opening
figure(21)
phi = 0;
hold off
plot3(X(1,:) + Xd(1,:), Y(1,:) - Yd(1,:)*sin(beta), Z(1,:) + Yd(1,:)*cos(beta), 'linewidth', 3);
axis equal vis3d
title('shell aperture');

% surface features and color
step_phi = 2*pi/(N_phi);
step_theta = 2*pi/(N_theta);
phi = 0:step_phi:N*2*pi;
theta = pi + (0:step_theta:2*pi);
x = a*exp(alpha*phi).*cos(phi);
y = a*exp(alpha*phi).*sin(phi);
z = -h*exp(alpha*phi);
[Theta, Phi] = meshgrid(theta, phi);
Xd = r1*cos(Theta);
Yd = r2*sin(Theta);
R = exp(alpha*Phi);
X = a*R.*cos(Phi);
Y = a*R.*sin(Phi);
Z = -h*R;

pattern = 1;
for i=1:length(f)
  pattern = pattern + A(i)*sin(f(i)*Phi);
end
R2 = exp(alpha*Phi) .* pattern; 
surf_X = X + cos(Phi) .* Xd .* R2 + Yd .* R2 .* sin(Phi) * sin(beta);
surf_Y = Y + sin(Phi) .* Xd .* R2 - Yd .* R2 .* cos(Phi) * sin(beta);
surf_Z = Z + Yd .* R2 * cos(beta);

figure(3)
hold off
plot(Phi/(2*pi), pattern, 'b');
hold on
plot(Phi/(2*pi), R2(:,1), 'r-');
xlabel('rotations');


figure(4)
hold off
g = surf(surf_X, surf_Y, surf_Z,'EdgeColor','none');
set(g, 'FaceLighting','phong');
sigma = 0.05;
m = 5;
kernal = ones(m,m)/m^2;
clear colors
colors(:,:,1) = 0.7 + sigma*conv2(randn(size(surf_X)), kernal, 'same');
colors(:,:,2) = 0.6 + sigma*conv2(randn(size(surf_X)), kernal, 'same');
colors(:,:,3) = 0.5 + sigma*conv2(randn(size(surf_X)), kernal, 'same');
set(g, 'CData', colors);
set(g, 'AmbientStrength', 0.3);
set(g, 'DiffuseStrength', 0.7);
set(g, 'SpecularStrength', 0.1);
set(g, 'SpecularExponent', 10);
set(g, 'SpecularColorReflectance', 0);

camlight right
colormap('gray')
axis off equal
az = 45;
el = 15;
view(az, el);



% filename = sprintf('Shell_Plots/shell_%06.3f_%06.3f_%06.3f_%06.3f_%06.3f_%06.3f.eps', ...
% 		   a, h, -alpha, r1, r2, beta);
% print('-depsc', filename)

% az = 0;
% el = 90;
% view(az, el);
% filename = sprintf('Shell_Plots/shell2_%06.3f_%06.3f_%06.3f_%06.3f_%06.3f_%06.3f.eps', ...
% 		   a, h, -alpha, r1, r2, beta);
% print('-depsc', filename)
