function conicdemo
%CONICDEMO Interactive demo of the conic sections generated by
% intersecting the positive and negative unit cones with an arbitrary
% plane.
% Tim Farajian 7/2005
% timfarajian@verizon.net
% Initialize data
hAx = []; hPlane = []; hCont = [];
[cx cy] = meshgrid(linspace(-1, 1, 50)); %Generate x and y data for interpolation
zcone = sqrt(cx.^2 + cy.^2); %Generate cone defined at same x and y as plane
AxisBigPos = [0.25 0.05 0.7 0.9];
AxisSmallPos = [0.02 0.75 0.2 0.2];
offset = -1/2; %Set default plane offset
normal = [1; 1; 1]; %Set default plane normal
% Create figure
hFig = figure('units', 'norm', 'pos', [0.05 0.1 0.9 0.8],'resize','off',...
'menu', 'none', 'name', 'Conic Sections', 'NumberTitle', 'off',...
'doublebuffer','off');
figColor = get(hFig, 'color');
% Create slider text boxes
hTxt(4) = uicontrol('position',[85 22 38 17],'string',offset);
hTxt(3) = uicontrol('position',[85 52 38 17],'string',normal(3));
hTxt(2) = uicontrol('position',[85 82 38 17],'string',normal(2));
hTxt(1) = uicontrol('position',[85 112 38 17],'string',normal(1));
set(hTxt,'style','text','backgroundcolor','w','fontsize',10);
% Create A,B,C,D text boxes
tmp(1)=uicontrol('position',[125 22 20 17], 'string','D');
tmp(2)=uicontrol('position',[125 52 20 17], 'string','C');
tmp(3)=uicontrol('position',[125 82 20 17], 'string','B');
tmp(4)=uicontrol('position',[125 112 20 17],'string','A');
tmp(5)=uicontrol('position',[0 140 170 17],'string','Ax + By + Cz + D = 0');
set(tmp,'style','Text','Back',figColor,'fontsize',10);
% Create sliders
hSld(4) = uicontrol('position',[20 20 60 20],'value',offset,'userdata',4);
hSld(3) = uicontrol('position',[20 50 60 20],'value',normal(3),'userdata',3);
hSld(2) = uicontrol('position',[20 80 60 20],'value',normal(2),'userdata',2);
hSld(1) = uicontrol('position',[20 110 60 20],'value',normal(1),'userdata',1);
set(hSld,'style','slider','min',-1,'max',1,'callback',@slidercall);
% Create Switch button
uicontrol('style','pushbutton','position',[70 400 60 30],...
'String','Switch','Callback',@SwitchAxes);
% Initialize axes
Init3dAxis;
Init2dAxis;
slidercall(hSld(1));
function Init3dAxis
hAx(1) = axes('position',AxisBigPos,'Xticklabel',[],...
'Yticklabel',[],'Zticklabel',[]);
view(3)
hold on
shading interp
axis([-1 1 -1 1 -1 1])
grid on
[x, y, z] = cylinder(linspace(0, 1), 50); %Generate cylinder
surf(x,y,z,'linestyle','none'); % Plot upper cone
surf(x,y,-z,'linestyle','none'); % Plot lower cone
hPlane = patch(1,1,1,1,'FaceColor','interp');
end
function Init2dAxis
hAx(2) = axes('pos',AxisSmallPos,'Xtickla',[],'Ytickla',[]);
grid on, hold on
axis([-1 1 -1 1])
end
%%%%%%%%%%%%%%%%%%%
% Slider callback %
%%%%%%%%%%%%%%%%%%%
function slidercall(hcb, ignore)
h = get(hcb, 'userdata'); %Get handle to corresponding text box
if h < 4 %If normal coefficient slider
normal(h) = get(hcb, 'value'); %Get value of slider
val = normal(h); %Assign to v
else
offset = get(hcb, 'value'); %Get value of offset
val = offset; %Assign to v
end
if abs(val) < 1e-4
str = '0'; %Ensure near 0 is shown as 0 in text box
else
str = sprintf('%0.3g', val); %Generate text box string
end
set(hTxt(h), 'string', str); %Set string of corresponding text box
axes(hAx(1));
Conic; %Plot 3-D plot
end
function Conic
% Plots the 3-D intersection of the plane with the unit cones
% by finding the 0-contour of the difference between the z data for the cone
% and the z data for the plane defined at the same x and y values, then
% interpolating to find the z value of the intersection.
%Ensure normal(3) ~= 0
if normal(3)>=0 && normal(3)<1e-4
normal(3) = 1e-4;
elseif normal(3)<0 && normal(3)>1e-4
normal(3) = -1e-4;
end
% Determine corners of plane
xx = [-1 1 1 -1]; yy = [-1 -1 1 1]; %Define 4 corners
zz = -1/normal(3)*(normal(1)*xx + normal(2)*yy + offset); %Generate plane
c = zz; c(c>1) = 1; c(c<-1) = -1; %Define colorscale
set(hPlane,'Xdata',xx,'Ydata',yy,'Zdata',zz,'Cdata',c);
% Generate intersection
zplane = -1/normal(3)*(normal(1)*cx + normal(2)*cy + offset); %Generate plane z-data
dz = zplane - zcone; %Compute difference between surfaces
%%%Determine rotation matrices%%%
[tht phi] = cart2sph(normal(1), normal(2), normal(3)); %Determine cartesian coords of normal
R1 = [cos(tht) sin(tht) 0; -sin(tht) cos(tht) 0; 0 0 1]; %Determine first rot matrix
R2 = [sin(phi) 0 -cos(phi); 0 1 0; cos(phi) 0 sin(phi)]; %Determine second rot matrix
M = -offset/sqrt(dot(normal(1:3),normal(1:3))); %Determine scaling factor
%%%First time through is for positive cone, second for negative%%%%
delete(hCont); hCont = [];
for k = 1:2
Z = []; %Initialize Z
try
C = contours(cx, cy, dz, [0 0]); %Compute 0 contour
catch
C = [];
end
if ~isempty(C)
[X Y v] = contoursconvert(C); %Convert to usable form
%%%Interpolate to find z values of intersection
for j = 1:length(v) %For each line
tv = 1:v(j); %Create indexing vector
Z(tv, j) = interp2(cx, cy, zplane, X(tv, j), Y(tv, j));
end
X(X>1|X<-1) = NaN; Y(Y>1|Y<-1) = NaN; Z(Z>1|Z<-1) = NaN; %Restrict display of line
th = plot3(hAx(1),X, Y, Z, 'k', 'linewidth', 2); % Plot intersection
hCont = [hCont; th]; % Append handles
% Create 2-D conic section
XYZ = [X(:) Y(:) Z(:)].'; %Create [x; y; z] matrix
xyz = R2*R1*XYZ; %Perform rotation
X = reshape(xyz(1, :), size(X)); %%
Y = reshape(xyz(2, :), size(X)); %%Convert back to mesh
Z = reshape(xyz(3, :), size(X)); %%
th = plot3(hAx(2),X, Y, Z, 'k'); % Plot intersection
hCont = [hCont; th]; % Append handle
axis(hAx(2),[-1 1 -1 1])
end
dz = zplane + zcone; %Now set to compute for negative cone
end
end
% Switch the position of the two axes
function SwitchAxes(ignore,ignore2)
pos1 = get(hAx(1),'Position');
set(hAx(1),'Position',get(hAx(2),'Position'));
set(hAx(2),'Position',pos1);
end
end
function [x, y, v] = contoursconvert(C)
%Converts the output of CONTOURS to a form usable by PLOT
%
% C = [info1 line1 info2 line2 ...]
% where info# are 2-by-1 vectors of information about line# in the form
% [ignore; N#], and line# is the 2-by-N# vector of [x; y] data
%
%Converts to the form:
%
% [x y] where the columns of x and y each define a different line. For
% lines where N# is less than the max(N#), NaN values fill the remainder of
% the columns.
%
% v is the number of data points for each line.
limit = size(C,2);
i = 1;
n = 0;
while(i < limit)
n = n + 1;
I(n) = i;
v(n) = C(2,i);
i = i + 1 + C(2, i);
end
x = NaN*ones(n, max(v));
y = x;
for n = 1:length(v)
x(n,1:v(n)) = C(1,I(n)+1:I(n)+v(n));
y(n,1:v(n)) = C(2,I(n)+1:I(n)+v(n));
end
x = x.'; y = y.';
end
