// The Klein model
function limitCircle()
	t=0:0.01:2*%pi;
	plot(cos(t),sin(t),'r:');
endfunction

function a = angle(p1,q1,p2,q2)
	if ( p1 == p2 ) then
		a = %pi/2;
	else
		a = atan((q2-q1)/(p2-p1));
	end
endfunction

function q = rotate90(p)
	q = [ -p(2), p(1)];
endfunction

function n = normal(p, q)
	// unit vector normal to p2-p1
	if p == q then
		error("Distint points are needed")
	end
	n = rotate90(q-p);
	n = n/norm(n);
endfunction

function r = dotprod(p,q)
	r = sum( p .* q);
endfunction

// n is a unit vector
function y = vectorProj(p, n)
	y = dotprod(p,n) * n;
endfunction

// n is a unit vector
function y = scalarProj(p, n)
	y = dotprod(p,n);
endfunction

function [e1, e2] = endPoints(n,d)
	m = rotate90(n);
	lambda = sqrt(1 - d*d);
	e1 = d*n + lambda*m;
	e2 = d*n - lambda*m;
endfunction

function plotLineSeq(p,q,c)
	plot([p(1),q(1)],[p(2),q(2)],c);
endfunction

function plotPoint(p, name)
	plot(p(1),p(2),'ko');
	xstring(p(1),p(2),name);
endfunction

function p=pole(n,d)
	p = n / d;
endfunction

p=[0.5,0.3];
q=[0.1,0.6];

clf
limitCircle()
n=normal(p,q);
d=scalarProj(p,n);
[e1,e2] = endPoints(n,d);
plotLineSeq(e1,e2,'bo-')
plotPoint(p,"P");
plotPoint(q,"Q");
npole = pole(n,d);
plotPoint(npole,"Pole");
plotLineSeq(e1,npole,'m:')
plotLineSeq(e2,npole,'m:')

a = gca();
a.isoview="on";
a.box="off";
a.axes_visible=["off","off","off"];
//a.data_bounds=[-1.1, -1.1; 1.1, 1.1];
