function demo_conic_section
% graphical input of a control polygon
% for a conic section

% Initialize data
global AMB_Data;
AMB_Data=[];
AMB_Data.plotpoints=300;
AMB_Data.c=zeros(0,2);
AMB_Data.log_w=0;                 
AMB_Data.w=1;                 
AMB_Data.handles=[0 0 0 0 0];
AMB_Data.movetype=0;
AMB_Data.moveind=0;
AMB_Data.uicontrols=[];
AMB_Data.axis=0;

% Initialize GUI
clf reset;
set(gcf,'WindowButtonDownFcn',@myDownfct);
set(gcf,'WindowButtonMotionFcn',@myMovefct);
set(gcf,'WindowButtonUpFcn',@myUpfct);

% Text for title
uicontrol('Units','normalized', ...
  	'BackgroundColor',get(gcf,'Color'), ...
	'HorizontalAlignment','center', ...
        'FontWeight','bold',...
        'FontSize',20,...
	'Position',[.05,.91,.9,.07],...
	'String','Conic Section',...
	'Style','text');

% Text for w
AMB_Data.uicontrols(1)=uicontrol('Units','normalized', ...
	'BackgroundColor',get(gcf,'Color'), ...
	'HorizontalAlignment','left', ...
	'Position',[.5,.84,.3,.05],...
	'String','w=1.0000, parabola', ...
	'Style','text');

% Slider for w
AMB_Data.uicontrols(2)=uicontrol('Units','normalized', ...
	'BackgroundColor',get(gcf,'Color'), ...
	'Position',[.05,.85,.4,.05], ...
	'Min',-9,'Max',9,'Value',AMB_Data.log_w,...
	'Sliderstep',[.05,.01],...
	'Callback',@slider_update,...
	'Style','slider');

% Info text 
AMB_Data.uicontrols(3)=uicontrol('Units','normalized', ...
	'BackgroundColor',get(gcf,'Color'), ...
	'HorizontalAlignment','right', ...
	'FontWeight','bold', ...
	'ListboxTop',0, ...
	'Position',[.05,.01,.9,.05],...
	'String','Select controlpolygon (3 points) ', ...
	'Style','text');

% Button for new curve
uicontrol('Units','normalized', ...
	'BackgroundColor',get(gcf,'Color'), ...
	'String','Restart', ...
	'Position',[.83,.85,.12,.05], ...
	'Callback',@restart,...
	'Style','pushbutton');

% Axis for curve
AMB_Data.axis(1)=subplot('position',[.05,.07,.9,.75]);
axis([-.01,1.01,-.01,1.01]);
axis manual;
axis off
hold on;

% frame 
fill([-.1,1.1,1.1,-.1],[-.1,-.1,1.1,1.1],[1,1,1],'edgecolor','none');
plot([0,1,1,0,0],[0,0,1,1,0],'-k');

%%% Callback Functions %%%

% Reset button pressed
function restart(varargin) 
global AMB_Data;
AMB_Data.c=zeros(0,2);
AMB_Data.step=0;
AMB_Data.movetype=0;
AMB_Data.moveind=0;
for k=AMB_Data.handles
  if k,delete(k);end
end
AMB_Data.handles=[0,0,0,0,0];
set(AMB_Data.uicontrols(3),'Visible','on');  
set(AMB_Data.uicontrols(2),'Value',0);
AMB_Data.w=1;
AMB_Data.log_w=1;
set(AMB_Data.uicontrols(1),'String','w=1.0000, parabola');

function slider_update(varargin) 
global AMB_Data;
% Read t value
AMB_Data.log_w = get(AMB_Data.uicontrols(2),'Value');
AMB_Data.w = 10^(AMB_Data.log_w/10);

% display w and type
if AMB_Data.log_w == 0
  set(AMB_Data.uicontrols(1),'String','w=1.0000, parabola');
end
if AMB_Data.log_w < 0
  set(AMB_Data.uicontrols(1),'String',sprintf('w=%06.4f, ellipse',AMB_Data.w));
end
if AMB_Data.log_w > 0
  set(AMB_Data.uicontrols(1),'String',sprintf('w=%06.4f, hyperbola',AMB_Data.w));
end
myplot;

% Mouse button pressed
function myDownfct(varargin)
global AMB_Data;

if strcmp(get(gcf,'SelectionType'),'normal') % Left mouse: new point 
  if size(AMB_Data.c,1)<3
    p=get(AMB_Data.axis(1),'CurrentPoint');
    p=p(1,1:2);    
    if min(p) >= 0 & max(p) <=1
      AMB_Data.c=[AMB_Data.c;p];
      n=size(AMB_Data.c,1);
      AMB_Data.movetype=1;
      AMB_Data.moveind=n;
    end
    myplot;
  end
end

% Mouse button released
function myUpfct(varargin)
global AMB_Data;
AMB_Data.movetype=0;
AMB_Data.moveind=0;
myplot;

% Mouse movement
function myMovefct(varargin)
global AMB_Data;
if AMB_Data.movetype
  p=get(AMB_Data.axis(1),'CurrentPoint');
  p=p(1,1:2);
  if min(p) >= 0 & max(p) <=1 
    AMB_Data.c(AMB_Data.moveind,:)=p;
    myplot;
  end 
end     
  
%%% Internal Funcions %%%

function myplot
global AMB_Data;

subplot(AMB_Data.axis(1));
panz = AMB_Data.plotpoints;
n=size(AMB_Data.c,1);

% delete curve and control polygon and c-steps 
if n < 2 
  if AMB_Data.handles(2)~=0  
    delete(AMB_Data.handles(2));
    AMB_Data.handles(2)=0;
  end
end
if n < 1 
  if AMB_Data.handles(1)~=0  
    delete(AMB_Data.handles(1));
    AMB_Data.handles(1)=0;
  end
end

if n>0 
  % draw polygon
  if AMB_Data.handles(1)==0  
    AMB_Data.handles(1)=plot(AMB_Data.c(:,1),AMB_Data.c(:,2),'-ok');
  else
    set(AMB_Data.handles(1),'Xdata',AMB_Data.c(:,1));
    set(AMB_Data.handles(1),'Ydata',AMB_Data.c(:,2));
  end
end
if n==3 & AMB_Data.moveind==0
  % draw closed polygon
  if AMB_Data.handles(1)==0  
    AMB_Data.handles(1)=plot(AMB_Data.c([1,2,3,1],1),AMB_Data.c([1 2 3 1],2),'-ok');
  else
    set(AMB_Data.handles(1),'Xdata',AMB_Data.c([1,2,3,1],1));
    set(AMB_Data.handles(1),'Ydata',AMB_Data.c([1,2,3,1],2));
  end
  % draw median line
  m1=mean(AMB_Data.c([1,3],:));
  m2=(m1+AMB_Data.w*AMB_Data.c(2,:))/(1+AMB_Data.w);
  if AMB_Data.handles(4)==0  
    AMB_Data.handles(4)=plot([m1(1),m2(1),AMB_Data.c(2,1)],...
			    [m1(2),m2(2),AMB_Data.c(2,2)],...
			    '.-','color',[0,.6,0],'Markersize',20);
  else
    set(AMB_Data.handles(4),'Xdata',[m1(1),m2(1),AMB_Data.c(2,1)]);
    set(AMB_Data.handles(4),'Ydata',[m1(2),m2(2),AMB_Data.c(2,2)]);
  end

  % draw parabola
  t=linspace(0,1,panz);
  p=rational_bezier_curve(AMB_Data.c,[1;1;1],t');  
  if AMB_Data.handles(3)==0  
    AMB_Data.handles(3)=plot(p(:,1),p(:,2),':k','Linewidth',1);  
  else
    set(AMB_Data.handles(3),'Xdata',p(:,1));
    set(AMB_Data.handles(3),'Ydata',p(:,2));
  end
 
  % draw curve
  t1=5.^linspace(0.1,2.2,panz);
  t=[1+t1(end),-fliplr(t1),linspace(-1,2,panz),1+t1];
   s=rational_bezier_curve(AMB_Data.c,[1;AMB_Data.w;1],t');  
  s(s<0)=NaN; s(s>1)=NaN;
        
  if AMB_Data.handles(2)==0  
    AMB_Data.handles(2)=plot(s(:,1),s(:,2),'-b','Linewidth',3);  
  else
    set(AMB_Data.handles(2),'Xdata',s(:,1));
    set(AMB_Data.handles(2),'Ydata',s(:,2));
  end  
  set(AMB_Data.uicontrols(3),'Visible','off');  
end

%%% Functions from AMB_Programs (included for standalone use) %%%

function [r,d] = rational_bezier_curve(c,w,t) 
% RATIONAL_BEZIER_CURVE: evaluation of a rational Bezier curve 
% [r,d] = rational_bezier_curve(c,w,t) 
% points r(j,:) at t(j) and weight points d(i,:) 
% for a rational Bezier curve 
% with control points c(k,:) and weights w

% convert to column vector (if necessary)
w = w(:);

% degree and dimension
[np1,dim] = size(c);

% homogeneous control points
W = repmat(w,1,dim);
cw = [W.*c w];

% evaluation
r = casteljau(cw,t);
r= r(:,1:end-1)./repmat(r(:,end),1,dim);

% weight points
K = [1:np1-1];
d = (W(K,:).*c(K,:)+W(K+1,:).*c(K+1,:))./(W(K,:)+W(K+1,:));

function p = casteljau(c,t)
% CASTELJAU: de Casteljau's algorithm for Bezier curves 
% p = casteljau(c,t)
% points p(j,:) on a Bezier curve with control points c(k,:)
% for parameters t(j) 

% convert to column vector (if necessary)
t = t(:); 
 
[np1,dim] = size(c);
% de Casteljau's recursion 
if np1 == 1
   p = (1+0*t)*c;
else
   p0 = casteljau(c(1:end-1,:),t); 
   p1 = casteljau(c(2:end,:),t); 
   t = repmat(t,1,dim); 
   p = (1-t).*p0 + t.*p1;
end

