// quadrilateral
t=[0:0.1:2*%pi]';
x=2*cos(t);
y=sin(t);
plot2d(x,y);

u=2*%pi*gsort(rand(1,4,'uniform'),'g','i');
v=[u,u(1),u(3),u(2),u(4)];
w=[2*cos(v);sin(v)]';
plot(w(:,1),w(:,2));

z(1,1)= 2*(sin(u(1))-sin(u(2)))/sin(u(1)-u(2));
z(1,2)=-1*(cos(u(1))-cos(u(2)))/sin(u(1)-u(2));
z(2,1)= 2*(sin(u(2))-sin(u(3)))/sin(u(2)-u(3));
z(2,2)=-1*(cos(u(2))-cos(u(3)))/sin(u(2)-u(3));
z(3,1)= 2*(sin(u(3))-sin(u(4)))/sin(u(3)-u(4));
z(3,2)=-1*(cos(u(3))-cos(u(4)))/sin(u(3)-u(4));
z(4,1)= 2*(sin(u(4))-sin(u(1)))/sin(u(4)-u(1));
z(4,2)=-1*(cos(u(4))-cos(u(1)))/sin(u(4)-u(1));
z(5,:)=z(1,:); z(6,:)=z(3,:); z(7,:)=z(2,:); z(8,:)=z(4,:);
plot(z(:,1),z(:,2),'r');
