skip to primary navigationskip to content
 

ConicDemo.m

Objective-C source code icon ConicDemo.m — Objective-C source code, 8 KB (8593 bytes)

File contents

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
% with modification by Cerian Webb & Atayab 2012/13


% 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.03 0.75 0.2 0.2];

offset = 1/2; %Set default plane offset
normal = [0; 0; 0]; %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(1) = uicontrol('position',[85 32 38 25],'string',normal(1));
hTxt(2) = uicontrol('position',[85 72 38 25],'string',normal(2));
hTxt(3) = uicontrol('position',[85 112 38 25],'string',offset);
set(hTxt,'style','text','backgroundcolor','w','fontsize',12);

% Create A,B,C text boxes
tmp(1)=uicontrol('position',[125 30 200 20], 'string','Rotate about y-axis');
tmp(2)=uicontrol('position',[125 70 200 20], 'string','Rotate about x-axis');
tmp(3)=uicontrol('position',[125 110 200 20],'string','Translate up & down z-axis');
set(tmp,'style','Text','Back',figColor,'fontsize',14);

% Create sliders
hSld2(1) = uicontrol('position',[20 115 60 20],'value',offset,'userdata',3);
hSld(2) = uicontrol('position',[20 75 60 20],'value',normal(2),'userdata',2);
hSld(1) = uicontrol('position',[20 35 60 20],'value',normal(1),'userdata',1);
%CRW:Changed range to allow more flexibility from [-1,1] to [-4 4]
set(hSld,'style','slider','min',-4,'max',4,'callback',@slidercall);
%CRW added different scale so z translation can be between [-1 1]
set(hSld2,'style','slider','min',-1,'max',1,'callback',@slidercall);

% Create Switch button
% CRW: Change background color of switch
uicontrol('style','pushbutton','BackgroundColor','w','FontSize',16,'position',[70 400 160 30],...
    'String','Switch graphs','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');
        
        xlabel('x - axis','FontSize',14);
        xlabh = get(gca,'XLabel');
        set(xlabh,'Position',get(xlabh,'Position') + [0 .5 0.1])
        ylabel('y - axis','FontSize',14);
        ylabh = get(gca,'YLabel');
        set(ylabh,'Position',get(ylabh,'Position') + [0.2 .5 0.2])
        zlabel('z - axis','FontSize',14);
        
    end

    function Init2dAxis
        hAx(2) = axes('pos',AxisSmallPos,'Xtickla',[],'Ytickla',[]);
        grid on, hold on
        axis square
        axis([-1 1 -1 1])
        xlabel('x - axis','FontSize',14);
        ylabel('y - axis','FontSize',14);
        
    end



%%%%%%%%%%%%%%%%%%%
% Slider callback %
%%%%%%%%%%%%%%%%%%%

    function slidercall(hcb, ignore)

        h = get(hcb, 'userdata'); %Get handle to corresponding text box
        if h < 3 %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(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(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), 1); %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:2),normal(1:2))); %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
                % CRW: Change line width where intersection plotted on cone
                th = plot3(hAx(1),X, Y, Z, 'r', 'linewidth', 3); % 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)); %%
                % CRW: Change line width and color on 2d plot
                th = plot3(hAx(2),X, Y, Z, 'r', 'linewidth', 3); % Plot intersection
                hCont = [hCont; th]; % Append handle
                %CRW: change axis limits as currently missing off some
                %circle from [-1 1 -1 1]
                axis(hAx(2),[-1.4 1.4 -1.4 1.4])
            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