% % 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)