function pt=comp_points(n)
syms a s d f;
a;
pt=[a s; d f];
pt(:,:)=[-1 1
         0 0];
for i=1:n    
    pt = erase_dubels(simplify([pt comp_ints(pt)]));    
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function int=comp_ints(pt);
n=size(pt,2);
int=[];
for i=1:n;
for j=i:n;
for k=max(j,i+1):n;
for l=max(k,j+1):n;
    if (~(i==j))&(~(k==l));
        int=[int int_ll(pt(:,[i j k l])) ...
            int_lc(pt(:,[i j k l])) int_lc(pt(:,[i j l k])) int_lc(pt(:,[k l i j ])) int_lc(pt(:,[k l j i ]))...
            int_cc(pt(:,[i j k l])) int_cc(pt(:,[j i k l])) int_cc(pt(:,[i j l k])) int_cc(pt(:,[j i l k])) ];
    end;
    if (~(i==k))&(~(j==l));
        int=[int int_ll(pt(:,[i k j l]))...
            int_lc(pt(:,[i k j l])) int_lc(pt(:,[i k l j])) int_lc(pt(:,[j l i k])) int_lc(pt(:,[j l k i]))...
            int_cc(pt(:,[i k j l])) int_cc(pt(:,[k i j l])) int_cc(pt(:,[i k l j])) int_cc(pt(:,[k i l j])) ];
    end;
    if (~(i==l))&(~(j==k));
        int=[int int_ll(pt(:,[i l j k]))...
            int_lc(pt(:,[i l j k])) int_lc(pt(:,[i l k j])) int_lc(pt(:,[j k i l])) int_lc(pt(:,[j k l i]))...
            int_cc(pt(:,[i l j k])) int_cc(pt(:,[l i j k])) int_cc(pt(:,[i l k j])) int_cc(pt(:,[l i k j])) ];
    end;
end
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function m = erase_dubels(m);
i=1;
n=size(m,2);
while i<size(m,2)
    m=[m(:,1:(i-1)) erase_firsts_like(m(:,i:end))];
    i=i+1;
end    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function m=erase_firsts_like(m)
n=size(m,2);
m(:,find(all(simplify(m(:,1)*ones(1,(n-1))-m(:,2:end))==0))+1)=[];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function pt=int_ll(pts)
v=pts(:,1)-pts(:,3);
d1=pts(:,2)-pts(:,1);
d2=pts(:,4)-pts(:,3);
den= det([d1 d2]);
if simplify(den)==0;
    pt=[];
    return
end
t=-det([v d2])/den;
pt=pts(:,1)+d1*t;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function pt=int_lc(pts)
v=pts(:,1)-pts(:,3);
d=pts(:,2)-pts(:,1);
d2=(pts(:,4)-pts(:,3));
r2=d2'*d2;
a=d'*d;
b=2*v'*d;
c=v'*v-r2;
del=sqrt(b*b-4*a*c);

if ~(del==del');
    pt=[];
    return
end
tow_a=2*a;
if (simplify(del)==0);
    pt=pts(:,1)-d*b/tow_a;
    return
end
t1=(-b+del)/tow_a;
t2=(-b-del)/tow_a;
pt(:,1)=pts(:,1)+d*t1;
pt(:,2)=pts(:,1)+d*t2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function pt=int_cc(pts)
v=pts(:,1)-pts(:,3);
b=sqrt(v'*v);
if (b==0);
    pt=[];
    return
end
d1=pts(:,2)-pts(:,1);
r1=sqrt(d1'*d1);
d2=(pts(:,4)-pts(:,3));
r2=sqrt(d2'*d2);

p=(b+r1+r2)/2;
s=sqrt(p*(p-b)*(p-r1)*(p-r2));
if ~(s==s');
    pt=[];
    return
end
h=2*s/b;

p1=sqrt(r1*r1-h*h);
v=-v/b;
pt0=pts(:,1)+v*p1;
if (simplify(h)==0);
    pt=pt0;
    return
end

v=[v(2); -v(1)];
pt=[pt0+v*h,pt0-v*h];