clc,
clear all;
close all
tic

%Poloha letadla
Px1=245;
Py1=100;
P=[Px1,Py1];
%%Polohy stanic 
x1=0;
y1=0;
x2=400;
y2=500;
x3=600;
y3=100;
Stanice1=[x1,y1];
Stanice2=[x2,y2];
Stanice3=[x3,y3];


%Vypocet casu za ktery dorazi signal z letadla do stanice
v=300000000;%rychlost sireni 
SP1=((Px1-x1)^2+(Py1-y1)^2)^(1/2);
SP2=((Px1-x2)^2+(Py1-y2)^2)^(1/2);
SP3=((Px1-x3)^2+(Py1-y3)^2)^(1/2);

t1=SP1/v;
t2=SP2/v;
t3=SP3/v;

%Nastaveni zacatku a konce intervalu, ve kterem se hleda prusecik
%PS: Konec intervalu nastavte na polovinu pozadovane hodnoty(v algoritmu je
%skok po dvou)Priklad: chci od 100 do 500, Nastavim: Zacatek=100;Konec=250
%PS2:cim vetsi interval tim delsi vypocetni cas
Zacatek=480;
Konec=75;

%Polomer je hodnota vzdalenosti, pod kterou se bod vyhodnoti jako
%prusecik, viz BP podkapitola 3.5.2

Polomer = 35;
%Poznamka: kdyz spatne vyhodnoti program plochu blizko teziste je nutne
%zapnout podminku na radku 4156
% Nastaveni velikosti casove odchylky
dt=50e-09;
t1p=t1+dt;
t2p=t2+dt;
t3p=t3+dt;
t1m=abs(t1-dt);
t2m=abs(t2-dt);
t3m=abs(t3-dt);

%%Vypocet e
e12=((x2-x1)^2+(y2-y1)^2)^(1/2)/2;
e13=((x3-x1)^2+(y3-y1)^2)^(1/2)/2;
e23=((x3-x2)^2+(y3-y2)^2)^(1/2)/2;

%%Vypocet x

x12=(2*e12-abs(t1-t2)*v)/2;
x13=(2*e13-abs(t3-t1)*v)/2;
x23=(2*e23-abs(t3-t2)*v)/2;

x12A=(2*e12-abs(t1m-t2p)*v)/2;
x13A=(2*e13-abs(t3m-t1p)*v)/2;
x23A=(2*e23-abs(t3m-t2p)*v)/2;

x12B=(2*e12-abs(t1p-t2m)*v)/2;
x13B=(2*e13-abs(t3p-t1m)*v)/2;
x23B=(2*e23-abs(t3p-t2m)*v)/2;

% vypocet a
a12=e12-x12;
a13=e13-x13;
a23=e23-x23;

a12A=e12-abs(x12A);
a13A=e13-abs(x13A);
a23A=e23-abs(x23A);

a12B=e12-abs(x12B);
a13B=e13-abs(x13B);
a23B=e23-abs(x23B);

% vypocet b

b12=((e12)^2-(a12)^2)^(1/2);
b13=((e13)^2-(a13)^2)^(1/2);
b23=((e23)^2-(a23)^2)^(1/2);

b12A=((e12)^2-(a12A)^2)^(1/2);
b13A=((e13)^2-(a13A)^2)^(1/2);
b23A=((e23)^2-(a23A)^2)^(1/2);

b12B=((e12)^2-(a12A)^2)^(1/2);
b13B=((e13)^2-(a13A)^2)^(1/2);
b23B=((e23)^2-(a23A)^2)^(1/2);


%%Vypocet polohy stredu stanic S[m,n]
m12=(x2+x1)/2;
n12=(y2+y1)/2;
m13=(x3+x1)/2;
n13=(y3+y1)/2;
m23=(x3+x2)/2;
n23=(y3+y2)/2;
Stred12=[m12,n12];
Stred13=[m13,n13];
Stred23=[m23,n23];

%%vypocet uhlu q (mezi osou x a primkou protinajici body)

q12=atand(y2/x2);
q13=atand(y3/x3);
q23=atand((y2-y3)/(x2-x3));

%%parametry hyperboly mezi stanicemi 1 a 2
m=m12;
n=n12;
a=a12;
b=b12;

sinq12=sind(q12);
cosq12=cosd(q12);
% vypocet hyperboly (Viz tøetí kapitola BP)
% zakladni rovnice(osove rovnobezna s osami x=0 a y=0):
% (x^2/a^2)-(y^2/b^2)-1=0
% posunuti do stredu spojnice stanic 1 a 2 S12[m,n]:
% x=(x-n) 
% y=(y-m)
% otoceni o uhel(q) svirajici spojnice bodu a osa x=0:
% x=(x*cos(q)+y*sin(q))
% y=(-x*sin(q)+y*cos(q))

%hyperbola f12
syms xg yg
f12 = @(x,y) (((x-m)*cosq12+(y-n12)*sinq12).^2)/a^2 -(((-(x-m)*sinq12+(y-n)*cosq12).^2)/b^2)-1;
hold on
eq12= (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12^2)==1  ;
h12=ezplot(f12,[-100 700]);
set(h12, 'Color','r' )
%parametry hyperboly mezi stanicemi 1 a 3
m=m13;
n=n13;
a=a13;
b=b13;

sinq13=sind(q13);
cosq13=cosd(q13);
%hyperbola f13
f13 = @(x,y) (((x-m)*cosq13+(y-n)*sinq13).^2)/a^2 -(((-(x-m)*sinq13+(y-n)*cosq13).^2)/b^2)-1;
eq13 = (((xg-m13)*cosq13+(yg-n13)*sinq13).^2)/a13^2 -(((-(xg-m13)*sinq13+(yg-n13)*cosq13).^2)/b13^2)==1;
h13 =ezplot(f13,[-100 700]);
set(h13, 'Color','r' )
%parametry hyperboly mezi stanicemi 2 a 3
m=m23;
n=n23;
a=a23;
b=b23;

sinq=sind(q23);
cosq=cosd(q23);
%%hyperbola f23
f23 = @(x,y) abs(((x-m)*cosq+(y-n)*sinq).^2)/a^2 -(((-(x-m)*sinq+(y-n)*cosq).^2)/b^2)-1;
eq23=(((xg-m)*cosq+(yg-n)*sinq).^2)/a^2 -(((-(xg-m)*sinq+(yg-n)*cosq).^2)/b^2)==1;
h =ezplot(f23,[-100 700]);
set(h, 'Color','r' )

%%vypocet 4 pruseciku dvou puvodnich hyperbol
[Xg,Yg] = solve(eq12,eq13);%2 hyperboly maji 4 pruseciky o x-ove a y-ove souradnici
xg1=double(round(vpa(Xg(1,1))));
xg2=double(round(vpa(Xg(2,1))));
xg3=double(round(vpa(Xg(3,1))));
xg4=double(round(vpa(Xg(4,1))));
yg1=double(round(vpa(Yg(1,1))));
yg2=double(round(vpa(Yg(2,1))));
yg3=double(round(vpa(Yg(3,1))));
yg4=double(round(vpa(Yg(4,1))));

 x1zkouska = imag(xg1);
 x2zkouska = imag(xg2); 
% Zde zacina overeni jestli funkce solve nevyhodila nesmysl, teorie
% v kapitole 3.5.
 if (x1zkouska>0 || x2zkouska>0 || xg1 >3000||xg1 <-3000||yg1 <-3000||yg1 <-3000 || xg2 >3000||xg2 <-3000||yg2 <-3000||yg2 <-3000)   
  xg=Zacatek;
  IntersectionPoint2=[];
  pom2=0;
  pom3=1;
  IntersectionPoint3=[];
  pomA = 0;
  pomB = 0;
  pomC = 0;
  pomD = 0;  
  for i=1:Konec

eq12G= (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12^2)==1  ;
GH1=round(vpa(solve(eq12G,yg)));

Point1=[xg;GH1(1,1)];
Point2=[xg;GH1(2,1)];
eq13G = (((xg-m13)*cosq13+(yg-n13)*sinq13).^2)/a13^2 -(((-(xg-m13)*sinq13+(yg-n13)*cosq13).^2)/b13^2)==1;
GH2=round(vpa(solve(eq13G,yg)));
Point3=[xg;GH2(1,1)];
Point4=[xg;GH2(2,1)];

d13=sqrt((xg-xg).^2+(GH1(1,1)-GH2(1,1)).^2);
d14=sqrt((xg-xg).^2+(GH1(1,1)-GH2(2,1)).^2);
d23=sqrt((xg-xg).^2+(GH1(2,1)-GH2(1,1)).^2);
d24=sqrt((xg-xg).^2+(GH1(2,1)-GH2(2,1)).^2);
try

if d13<Polomer
    Podezrely = [xg;GH1(1,1)];
end
if d14<Polomer
    Podezrely = [xg;GH1(1,1)];
end
if d23<Polomer
    Podezrely = [xg;GH1(2,1)];
end
if d24<Polomer
    Podezrely = [xg;GH1(2,1)];
end
if Podezrely(1,1)==xg %kdyz je nejaky Podezrely

try %zkus
    
[dm,dn] = size(IntersectionPoint2);%zjisti kolik je uz pruseciku v mnozine
dl=dm/2; %pomocny vypocet, pruseiku je v mnozine tolik kolik je radku deleno dvema, protoze kazdy podezrely ma na 1 radku x-ovou a na druhem y-ovou souradnici

if dl==0% kdyz v mnozine zadny jeste neni tak priradine prvniho podezrelyho jako prusecik
    IntersectionPoint2 = [Podezrely];
    Podezrely=null; 
else %kdyz tam zadny neni
    
for j=1:dl % cyklus, ktery se opakuje tolikrat kolik je pruseciku v mnozine
    
    blabla = IntersectionPoint2(pom3,1); %vyber spravneho pruseciku
    
   
    if    (blabla+2+pom2)==xg 
      switch pom3
        case 1
        pomA=pomA+1; 
        case 3
        pomB=pomB+1; 
        case 5
        pomC=pomC+1; 
        case 7
        pomD=pomD+1; 
        end
        pom2=pom2+2;
        Podezrely=null;
  % kdyz je tento podezrely(xg je jeho x-va souradnice) priliz blizko pruseciku..
  %Pozn: jedna se o to, ze pro jeden prusecik, nalezne algoritmus vice vhodnych bodu, nam staci 1, ostatni nevybereme    
  %      
        
        
    else
      
       IntersectionPoint2 = [IntersectionPoint2;Podezrely];
       pom3=pom3+2;
       pom2=0;
       Podezrely=null;
    end
    end

end
end
end
end
xg=xg+2;
  end
%zde zacina korekce pruseciku( zjisti se kolik bylo podezrelych u kazdeho
%pruseciku, pricte se polovicni hodnota k x a pro tu se dopocita y..Nachazi
%se zde osetreni kdy jsou nalezeny 4,2 a 1 prusecik

  try
MBx1=IntersectionPoint2(1,1);
MBy1=IntersectionPoint2(2,1);
MBx2=IntersectionPoint2(3,1);
MBy2=IntersectionPoint2(4,1);

%MB1
xg =MBx1+pomA; 
eq12G= (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12^2)==1  ;
GH1=round(vpa(solve(eq12G,yg)));

d1G=sqrt((MBy1-GH1(1,1)).^2);
d2G=sqrt((MBy1-GH1(2,1)).^2);
if d1G<d2G
IntersectionPoint3 = [xg;GH1(1,1)];

end
if d1G>d2G
IntersectionPoint3 = [xg;GH1(2,1)];    
end
%MB2
xg =MBx2+pomB; 
eq12G= (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12^2)==1  ;
GH1=round(vpa(solve(eq12G,yg)));

d1G=sqrt((MBy2-GH1(1,1)).^2);
d2G=sqrt((MBy2-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
try
MBx3=IntersectionPoint2(5,1);
MBy3=IntersectionPoint2(6,1);
MBx4=IntersectionPoint2(7,1);
MBy4=IntersectionPoint2(8,1);
%MB3
xg =MBx3+pomC; 
eq12G= (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12^2)==1  ;
GH1=round(vpa(solve(eq12G,yg)));

d1G=sqrt((MBy3-GH1(1,1)).^2);
d2G=sqrt((MBy3-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
%MB4
xg =MBx4+pomD; 
eq12G= (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12^2)==1  ;
GH1=round(vpa(solve(eq12G,yg)));

d1G=sqrt((MBy4-GH1(1,1)).^2);
d2G=sqrt((MBy4-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
end
end
 try 
xg1=round(vpa(IntersectionPoint3(1,1)));
xg2=round(vpa(IntersectionPoint3(3,1)));
xg3=round(vpa(IntersectionPoint3(5,1)));
xg4=round(vpa(IntersectionPoint3(7,1)));
yg1=round(vpa(IntersectionPoint3(2,1)));
yg2=round(vpa(IntersectionPoint3(4,1)));
yg3=round(vpa(IntersectionPoint3(6,1)));
yg4=round(vpa(IntersectionPoint3(8,1)));

 catch
     try
xg1=round(vpa(IntersectionPoint3(1,1)));
xg2=round(vpa(IntersectionPoint3(3,1)));
yg1=round(vpa(IntersectionPoint3(2,1)));
yg2=round(vpa(IntersectionPoint3(4,1)));
     catch
MBx1=round(vpa(IntersectionPoint2(1,1)));
MBy1=round(vpa(IntersectionPoint2(2,1)));
          
         
xg =MBx1+pomA; 
eq12G= (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12^2)==1  ;
GH1=round(vpa(solve(eq12G,yg)));
d1G=sqrt((MBy1-GH1(1,1)).^2);
d2G=sqrt((MBy1-GH1(2,1)).^2);
if d1G<d2G

xg1=xg;

yg1=GH1(1,1);

end
if d1G>d2G

xg1=xg;

yg1=GH1(2,1);
end
         
      
     end
    
 end
 syms xg;
 end
 

P1=[xg1,yg1];
P2=[xg2,yg2];
P3=[xg3,yg3];
P4=[xg4,yg4];
  
B = [P1 ; P2 ; P3 ; P4 ]; %Mnozina pruseciku


%Pomocna metoda kdyz jsou jen 2 pruseciky (vyradi zbyle dva pro dalsi
%vypocty)
%kdybych tam nechal komplexní èíslo, shodí to program v dalších èástech
x1zkouska = imag(xg3);
x2zkouska = imag(xg4); 
if (x1zkouska>0 || x2zkouska>0)

dP1S1=double(sqrt((x1-xg1).^2+(y1-yg1).^2));
dP2S1=double(sqrt((x1-xg2).^2+(y1-yg2).^2));
dP1S2=double(sqrt((x2-xg1).^2+(y2-yg1).^2));
dP2S2=double(sqrt((x2-xg2).^2+(y2-yg2).^2));
%*******************************Tohle jeste poresit*********500,100;
if t1<t2
    if dP1S1<dP2S1
        
        IntercestionalPoint= [xg1, yg1];
        
        
    else
        IntercestionalPoint= [xg2, yg2];
    end
else
    
    if dP1S1<dP2S1
        if dP1S2 > dP2S2
        IntercestionalPoint= [xg2, yg2];
        else
        IntercestionalPoint= [xg1, yg1];    
        end
    else
       if dP1S2 < dP2S2
        IntercestionalPoint= [xg2, yg2];
        else
        IntercestionalPoint= [xg1, yg1];    
        end
    end
end
C= IntercestionalPoint;
else
dP1S1=double(sqrt((x1-xg1).^2+(y1-yg1).^2));
dP1S2=double(sqrt((x2-xg1).^2+(y2-yg1).^2));
dP1S3=double(sqrt((x3-xg1).^2+(y3-yg1).^2));
dP2S1=double(sqrt((x1-xg2).^2+(y1-yg2).^2));
dP2S2=double(sqrt((x2-xg2).^2+(y2-yg2).^2));
dP2S3=double(sqrt((x3-xg2).^2+(y3-yg2).^2));
dP3S1=double(sqrt((x1-xg3).^2+(y1-yg3).^2));
dP3S2=double(sqrt((x2-xg3).^2+(y2-yg3).^2));
dP3S3=double(sqrt((x3-xg3).^2+(y3-yg3).^2));
dP4S1=double(sqrt((x1-xg4).^2+(y1-yg4).^2));
dP4S2=double(sqrt((x2-xg4).^2+(y2-yg4).^2));
dP4S3=double(sqrt((x3-xg4).^2+(y3-yg4).^2));
Points=[];
Points2=[];
if t1<t2  %kdyz signal z letadla dorazi na prvni stanici drive nez na druhou...
    % 1. podminka) a tento prusecik je blize prvni stanici nez druhe...  
    % jde do "uzsiho vyberu", ve kterem zbyvaji 2 pruseciky
    if dP1S1<dP1S2  
      Points = [Points; P1]; 

    end  
    if dP2S1<dP2S2     
      Points = [Points; P2];
  
    end 
    if dP3S1<dP3S2        
      Points = [Points; P3];
   
    end 
    if dP4S1<dP4S2  
      Points = [Points; P4];

    end
    
    if t1<t3 % 2. podminka) a tento prusecik je blize prvni stanici nez treti
        
       if dP1S1<dP1S3 
    Points2 = [Points2; P1]; % jde do "uzsiho vyberu", ve kterem zbyvaji 2 pruseciky(jina skupina)
       end
       if dP2S1<dP2S3
    Points2 = [Points2; P2];
       end
       if dP3S1<dP3S3   
    Points2 = [Points2; P3];
       end
       if dP4S1<dP4S3    
    Points2 = [Points2; P4];
       end
    else     % opacny pripad druhe podminky
        
       if dP1S1>dP1S3 
    Points2 = [Points2; P1]; % jde do "uzsiho vyberu", ve kterem zbyvaji 2 pruseciky(jina skupina)
       end
       if dP2S1>dP2S3
    Points2 = [Points2; P2];
       end
       if dP3S1>dP3S3   
    Points2 = [Points2; P3];
       end
       if dP4S1>dP4S3    
    Points2 = [Points2; P4];
       end 
    end
else         %opacny pripad prvni podminky
    if dP1S1>dP1S2       
      Points = [Points; P1];
    end
    if dP2S1>dP2S2
      Points = [Points; P2];
    end
    if dP3S1>dP3S2  
      Points = [Points; P3];
    end
    if dP4S1>dP4S2  
      Points = [Points; P4];
    end
    
    if t1<t3  % opet 2. podminka) a tento prusecik je blize prvni stanici nez treti
       if dP1S1<dP1S3   
         Points2 = [Points2; P1];
       end
       if dP2S1<dP2S3   
          Points2 = [Points2; P2];
       end
       if dP3S1<dP3S3
          Points2 = [Points2; P3];
       end
       if dP4S1<dP4S3  
          Points2 = [Points2; P4];
       end
    else     % opacny pripad druhe podminky
        
       if dP1S1>dP1S3 
    Points2 = [Points2; P1]; % jde do "uzsiho vyberu", ve kterem zbyvaji 2 pruseciky(jina skupina)
       end
       if dP2S1>dP2S3
    Points2 = [Points2; P2];
       end
       if dP3S1>dP3S3   
    Points2 = [Points2; P3];
       end
       if dP4S1>dP4S3    
    Points2 = [Points2; P4];
       end 
    end 
end

%Zde se trochu nesikovne porovnavaj obe dve mnoziny a ten body ktery je v
%obou mnozinach(splnuje obe 2 podminky) je ten prusecik ktery hledame

try
Pointa11 = Points(1,1);
Pointa12 = Points(1,2);
Pointa21 = Points(2,1);
Pointa22 = Points(2,2);

Pointb11 = Points2(1,1);
Pointb12 = Points2(1,2);
Pointb21 = Points2(2,1);
Pointb22 = Points2(2,2);
if (Pointa11 == Pointb11 && Pointa12 == Pointb12)
    IntercestionalPoint= [Pointa11, Pointa12];
end
if (Pointa11 == Pointb21 && Pointa12 == Pointb22)
    IntercestionalPoint= [Pointa11, Pointa12];
end
if (Pointa21 == Pointb11 && Pointa22 == Pointb12)
    IntercestionalPoint= [Pointa21, Pointa22];
end
if (Pointa21 == Pointb21 && Pointa22 == Pointb22)
    IntercestionalPoint= [Pointa21, Pointa22];
end
%catch kdyz jsou jenom 2 pruseciky
catch 
Pointa11 = Points(1,1);
Pointa12 = Points(1,2);

    IntercestionalPoint= [Pointa11, Pointa12];
end
%Hledany prusecik:
C= [IntercestionalPoint ];
end



% V této èásti vytváøíme hyperboly, které odpovídají maximální èasové odchylce 
% pùvodních tøí hyperbol
%***************************************************************************





%%----parametry hyperboly A  mezi stanicemi 2 a 3----
m=m23;
n=n23;
a=a23A;
b=b23A;

sinq23=sind(q23);
cosq23=cosd(q23);
%hyperbola f23A
f23A = @(x,y) (((x-m)*cosq23+(y-n)*sinq23).^2)/a^2 -(((-(x-m)*sinq23+(y-n)*cosq23).^2)/b^2)-1;
eq23a = (((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23A^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23A^2)==1;


%%----parametry hyperboly B  mezi stanicemi 2 a 3----
m=m23;
n=n23;
a=a23B;
b=b23B;

sinq23=sind(q23);
cosq23=cosd(q23);
%hyperbola f23B
f23B = @(x,y) (((x-m)*cosq23+(y-n)*sinq23).^2)/a^2 -(((-(x-m)*sinq23+(y-n)*cosq23).^2)/b^2)-1;
eq23b=(((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23B^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23B^2)==1; 


%%----parametry hyperboly A  mezi stanicemi 1 a 2----
m=m12;
n=n12;
a=a12A;
b=b12A;

sinq12=sind(q12);
cosq12=cosd(q12);
%hyperbola f12A
f12A = @(x,y) (((x-m)*cosq12+(y-n)*sinq12).^2)/a^2 -(((-(x-m)*sinq12+(y-n)*cosq12).^2)/b^2)-1;
eq12a = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12A^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12A^2)==1;


%%----parametry hyperboly B  mezi stanicemi 1 a 2----
m=m12;
n=n12;
a=a12B;
b=b12B;

sinq12=sind(q12);
cosq12=cosd(q12);
%hyperbola f12B
f12B = @(x,y) (((x-m)*cosq12+(y-n)*sinq12).^2)/a^2 -(((-(x-m)*sinq12+(y-n)*cosq12).^2)/b^2)-1;
eq12b = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12B^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12B^2)==1;



%%----parametry hyperboly A  mezi stanicemi 1 a 3----
m=m13;
n=n13;
a=a13A;
b=b13A;

sinq13=sind(q13);
cosq13=cosd(q13);
%hyperbola f13A
f13A = @(x,y) (((x-m)*cosq13+(y-n)*sinq13).^2)/a^2 -(((-(x-m)*sinq13+(y-n)*cosq13).^2)/b^2)-1;
eq13a = (((xg-m13)*cosq13+(yg-n13)*sinq13).^2)/a13A^2 -(((-(xg-m13)*sinq13+(yg-n13)*cosq13).^2)/b13A^2)==1;


%%----parametry hyperboly B  mezi stanicemi 1 a 3----
m=m13;
n=n13;
a=a13B;
b=b13B;

sinq13=sind(q13);
cosq13=cosd(q13);
%hyperbola f13B
f13B = @(x,y) (((x-m)*cosq13+(y-n)*sinq13).^2)/a^2 -(((-(x-m)*sinq13+(y-n)*cosq13).^2)/b^2)-1;
eq13b = (((xg-m13)*cosq13+(yg-n13)*sinq13).^2)/a13B^2 -(((-(xg-m13)*sinq13+(yg-n13)*cosq13).^2)/b13B^2)==1;

%vykresleni techto 6-ti hyperbol do grafu
ezplot(f12A,[-100 700])
ezplot(f12B,[-100 700])
ezplot(f13A,[-100 700])
ezplot(f13B,[-100 700])
ezplot(f23A,[-100 700])
ezplot(f23B,[-100 700])

%Vypocty pruseciku pasem  tvorenych hyperbolama 
%**************************************************************************
%okomentovany pouze pro prvni bod, u zbylych 11 je princip stejny
IntercestionalPoint
ipx=IntercestionalPoint(1,1);
ipy=IntercestionalPoint(1,2);
disp('1)')
%1
%reseni pomoci soustavy  rovnic (BP kapitola 3.5.1)
[Xga,Yga] = solve(eq23b,eq13b);
xg1a=round(vpa(Xga(1,1)));
xg2a=round(vpa(Xga(2,1)));
xg3a=round(vpa(Xga(3,1)));
xg4a=round(vpa(Xga(4,1)));
yg1a=round(vpa(Yga(1,1)));
yg2a=round(vpa(Yga(2,1)));
yg3a=round(vpa(Yga(3,1)));
yg4a=round(vpa(Yga(4,1)));


% vyhodnoceni vzdalenosti pruseciku od letadla
d1aIP=double(sqrt((ipx-xg1a).^2+(ipy-yg1a).^2));
d2aIP=double(sqrt((ipx-xg2a).^2+(ipy-yg2a).^2));
d3aIP=double(sqrt((ipx-xg3a).^2+(ipy-yg3a).^2));
d4aIP=double(sqrt((ipx-xg4a).^2+(ipy-yg4a).^2));


try
    if (xg1a >3000||xg1a <-3000||yg1a <-3000||yg1a <-3000 || xg2a >3000||xg2a <-3000||yg2a <-3000||yg2a <-3000)
        block  %kdyz funkce solve vyhodila "nesmyslne" vysledky shodi vypocet pomoci solve a pouzije se druha metoda
    end
    %porovnani vzdalenosti,ten nejbližší k letounu je ten hledany
  if (d1aIP <d2aIP && d1aIP <d3aIP && d1aIP <d4aIP)
    Point12B23B = [xg1a, yg1a];
    d1 = d1aIP;
  end
  if (d2aIP <d1aIP && d2aIP <d3aIP && d2aIP <d4aIP)
    Point12B23B = [xg2a, yg2a];
     d1 = d2aIP;
  end
  if (d3aIP <d2aIP && d3aIP <d1aIP && d3aIP <d4aIP)
    Point12B23B = [xg3a, yg3a];
     d1 = d3aIP;
  end
  if (d4aIP <d2aIP && d4aIP <d3aIP && d4aIP <d1aIP)
    Point12B23B = [xg4a, yg4a];
     d1 = d4aIP;
  end
  %vysledny hledany prusecik
D= [Point12B23B];
catch
  %druha metoda hledani pruseciku (BP 3.5.2)
  xg=Zacatek;
  IntersectionPoint2=[];
  pom2=0;
  pom3=1;
  IntersectionPoint3=[];
  pomA = 0;
  pomB = 0;
  pomC = 0;
  pomD = 0;  
  for i=1:Konec

eq23bG=(((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23B^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23B^2)==1; 
GH1=round(vpa(solve(eq23bG,yg)));

Point1=[xg;GH1(1,1)];
Point2=[xg;GH1(2,1)];
eq13bG = (((xg-m13)*cosq13+(yg-n13)*sinq13).^2)/a13B^2 -(((-(xg-m13)*sinq13+(yg-n13)*cosq13).^2)/b13B^2)==1;
GH2=round(vpa(solve(eq13bG,yg)));
Point3=[xg;GH2(1,1)];
Point4=[xg;GH2(2,1)];

d13=sqrt((xg-xg).^2+(GH1(1,1)-GH2(1,1)).^2);
d14=sqrt((xg-xg).^2+(GH1(1,1)-GH2(2,1)).^2);
d23=sqrt((xg-xg).^2+(GH1(2,1)-GH2(1,1)).^2);
d24=sqrt((xg-xg).^2+(GH1(2,1)-GH2(2,1)).^2);
try

if d13<Polomer
    Podezrely = [xg;GH1(1,1)];
end
if d14<Polomer
    Podezrely = [xg;GH1(1,1)];
end
if d23<Polomer
    Podezrely = [xg;GH1(2,1)];
end
if d24<Polomer
    Podezrely = [xg;GH1(2,1)];
end
if Podezrely(1,1)==xg %kdyz je nejaky Podezrely

try%zkus
    
[dm,dn] = size(IntersectionPoint2);%zjisti kolik je uz pruseciku v mnozine
dl=dm/2; %pomocny vypocet, pruseciku je v mnozine tolik kolik je radku deleno dvema, protoze kazdy podezrely ma na 1 radku x-ovou a na druhem y-ovou souradnici

if dl==0% kdyz v mnozine zadny jeste neni tak priradine prvniho podezrelyho jako prusecik
    IntersectionPoint2 = [Podezrely];
    Podezrely=null; 
else %kdyz tam zadny neni
    
for j=1:dl % cyklus, ktery se opakuje tolikrat kolik je pruseciku v mnozine
    
    blabla = IntersectionPoint2(pom3,1); %vyber spravneho pruseciku
    
   
    if    (blabla+2+pom2)==xg 
      switch pom3
        case 1
        pomA=pomA+1; 
        case 3
        pomB=pomB+1; 
        case 5
        pomC=pomC+1; 
        case 7
        pomD=pomD+1; 
        end
        pom2=pom2+2;
        Podezrely=null;
  % kdyz je tento podezrely(xg je jeho x-va souradnice) priliz blizko pruseciku..
  %Pozn: jedna se o to, ze pro jeden prusecik, nalezne algoritmus vice vhodnych bodu, nam staci 1, ostatni nevybereme    
        
        
        
    else
      
       IntersectionPoint2 = [IntersectionPoint2;Podezrely];
       pom3=pom3+2;
       pom2=0;
       Podezrely=null;
    end
    end

end
end
end
end
xg=xg+2;
  end
  %zde dochazi ke korekci pruseciku ( mnozina podezrelych se omezi na 1
  %nejblizsi ( nejdrive znam hodnotu x a z rovnice dopocitam y)
try
MBx1=IntersectionPoint2(1,1);
MBy1=IntersectionPoint2(2,1);
MBx2=IntersectionPoint2(3,1);
MBy2=IntersectionPoint2(4,1);

%MB1
xg =MBx1+pomA; 
eq23bG=(((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23B^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23B^2)==1; 
GH1=round(vpa(solve(eq23bG,yg)));
d1G=sqrt((MBy1-GH1(1,1)).^2);
d2G=sqrt((MBy1-GH1(2,1)).^2);
if d1G<d2G
IntersectionPoint3 = [xg;GH1(1,1)];

end
if d1G>d2G
IntersectionPoint3 = [xg;GH1(2,1)];    
end
%MB2
xg =MBx2+pomB; 
eq23bG=(((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23B^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23B^2)==1; 
GH1=round(vpa(solve(eq23bG,yg)));
d1G=sqrt((MBy2-GH1(1,1)).^2);
d2G=sqrt((MBy2-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
try
MBx3=IntersectionPoint2(5,1);
MBy3=IntersectionPoint2(6,1);
MBx4=IntersectionPoint2(7,1);
MBy4=IntersectionPoint2(8,1);
%MB3
xg =MBx3+pomC; 
eq23bG=(((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23B^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23B^2)==1; 
GH1=round(vpa(solve(eq23bG,yg)));
d1G=sqrt((MBy3-GH1(1,1)).^2);
d2G=sqrt((MBy3-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
%MB4
xg =MBx4+pomD; 
eq23bG=(((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23B^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23B^2)==1; 
GH1=round(vpa(solve(eq23bG,yg)));
d1G=sqrt((MBy4-GH1(1,1)).^2);
d2G=sqrt((MBy4-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
end
end  
  % zde dochazi k vybrani z mnoziny pruseciku ten hledany(nejblizsi k
  % letounu) 
  %nejdrive zkouska pro 4 pruseciky
 try 
xg1a=round(vpa(IntersectionPoint3(1,1)));
xg2a=round(vpa(IntersectionPoint3(3,1)));
xg3a=round(vpa(IntersectionPoint3(5,1)));
xg4a=round(vpa(IntersectionPoint3(7,1)));
yg1a=round(vpa(IntersectionPoint3(2,1)));
yg2a=round(vpa(IntersectionPoint3(4,1)));
yg3a=round(vpa(IntersectionPoint3(6,1)));
yg4a=round(vpa(IntersectionPoint3(8,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);
d3aIP=sqrt((ipx-xg3a).^2+(ipy-yg3a).^2);
d4aIP=sqrt((ipx-xg4a).^2+(ipy-yg4a).^2);

  if (d1aIP <d2aIP && d1aIP <d3aIP && d1aIP <d4aIP)
    Point12B23B = [xg1a, yg1a];
    d1 = d1aIP;
  end
  if (d2aIP <d1aIP && d2aIP <d3aIP && d2aIP <d4aIP)
    Point12B23B = [xg2a, yg2a];
     d1 = d2aIP;
  end
  if (d3aIP <d2aIP && d3aIP <d1aIP && d3aIP <d4aIP)
    Point12B23B = [xg3a, yg3a];
     d1 = d3aIP;
  end
  if (d4aIP <d2aIP && d4aIP <d3aIP && d4aIP <d1aIP)
    Point12B23B = [xg4a, yg4a];
     d1 = d4aIP;
  end
D= [Point12B23B];
 catch
     %pote kdyz nejsou 4, zkusi jestli jsou 2
     try
xg1a=round(vpa(IntersectionPoint3(1,1)));
xg2a=round(vpa(IntersectionPoint3(3,1)));
yg1a=round(vpa(IntersectionPoint3(2,1)));
yg2a=round(vpa(IntersectionPoint3(4,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);


  if (d1aIP <d2aIP )
    Point12B23B = [xg1a, yg1a];
    d1 = d1aIP;
  end
  if (d2aIP <d1aIP)
    Point12B23B = [xg2a, yg2a];
     d1 = d2aIP;
  end
  
D= [Point12B23B];
     catch
    %kdyz metoda nasla pouze 1 prusecik
         MBx1=round(vpa(IntersectionPoint2(1,1)));
         MBy1=round(vpa(IntersectionPoint2(2,1)));
          
         
xg =MBx1+pomA; 
eq23bG=(((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23B^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23B^2)==1; 
GH1=round(vpa(solve(eq23bG,yg)));
d1G=sqrt((MBy1-GH1(1,1)).^2);
d2G=sqrt((MBy1-GH1(2,1)).^2);
if d1G<d2G
Point12B23B= [xg,GH1(1,1)];
d1=sqrt((ipx-xg).^2+(ipy-GH1(1,1)).^2);

end
if d1G>d2G
Point12B23B = [xg,GH1(2,1)];   
d1=sqrt((ipx-xg).^2+(ipy-GH1(2,1)).^2);
end
         
        D= [Point12B23B]; 
     end
 end
 syms xg;
end
disp('2)')
%2
[Xga,Yga] = solve(eq23a,eq13b);
xg1a=round(vpa(Xga(1,1)));
xg2a=round(vpa(Xga(2,1)));
xg3a=round(vpa(Xga(3,1)));
xg4a=round(vpa(Xga(4,1)));
yg1a=round(vpa(Yga(1,1)));
yg2a=round(vpa(Yga(2,1)));
yg3a=round(vpa(Yga(3,1)));
yg4a=round(vpa(Yga(4,1)));



d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);
d3aIP=sqrt((ipx-xg3a).^2+(ipy-yg3a).^2);
d4aIP=sqrt((ipx-xg4a).^2+(ipy-yg4a).^2);


try
    
  if (xg1a >3000||xg1a <-3000||yg1a <-3000||yg1a <-3000 || xg2a >3000||xg2a <-3000||yg2a <-3000||yg2a <-3000)
        block 
   end
  if (d1aIP <d2aIP && d1aIP <d3aIP && d1aIP <d4aIP)
    Point12B23B = [xg1a, yg1a];
    d2 = d1aIP;
  end
  if (d2aIP <d1aIP && d2aIP <d3aIP && d2aIP <d4aIP)
    Point12B23B = [xg2a, yg2a];
     d2 = d2aIP;
  end
  if (d3aIP <d2aIP && d3aIP <d1aIP && d3aIP <d4aIP)
    Point12B23B = [xg3a, yg3a];
     d2 = d3aIP;
  end
  if (d4aIP <d2aIP && d4aIP <d3aIP && d4aIP <d1aIP)
    Point12B23B = [xg4a, yg4a];
     d2 = d4aIP;
  end
D= [D ; Point12B23B];



catch
  xg=Zacatek;
  IntersectionPoint2=[];
  pom2=0;
  pom3=1;

  IntersectionPoint3=[];
  pomA = 0;
  pomB = 0;
  pomC = 0;
  pomD = 0;  
  for i=1:Konec

eq23aG = (((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23A^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23A^2)==1;
GH1=round(vpa(solve(eq23aG,yg)));

Point1=[xg;GH1(1,1)];
Point2=[xg;GH1(2,1)];
eq13bG = (((xg-m13)*cosq13+(yg-n13)*sinq13).^2)/a13B^2 -(((-(xg-m13)*sinq13+(yg-n13)*cosq13).^2)/b13B^2)==1;
GH2=round(vpa(solve(eq13bG,yg)));
Point3=[xg;GH2(1,1)];
Point4=[xg;GH2(2,1)];

d13=sqrt((xg-xg).^2+(GH1(1,1)-GH2(1,1)).^2);
d14=sqrt((xg-xg).^2+(GH1(1,1)-GH2(2,1)).^2);
d23=sqrt((xg-xg).^2+(GH1(2,1)-GH2(1,1)).^2);
d24=sqrt((xg-xg).^2+(GH1(2,1)-GH2(2,1)).^2);
try

if d13<Polomer
    Podezrely = [xg;GH1(1,1)];
end
if d14<Polomer
    Podezrely = [xg;GH1(1,1)];
end
if d23<Polomer
    Podezrely = [xg;GH1(2,1)];
end
if d24<Polomer
    Podezrely = [xg;GH1(2,1)];
end
if Podezrely(1,1)==xg %kdyz je nejaky Podezrely

try%zkus
    
[dm,dn] = size(IntersectionPoint2);%zjisti kolik je uz pruseciku v mnozine
dl=dm/2; %pomocny vypocet, pruseiku je v mnozine tolik kolik je radku deleno dvema, protoze kazdy podezrely ma na 1 radku x-ovou a na druhem y-ovou souradnici

if dl==0% kdyz v mnozine zadny jeste neni tak priradine prvniho podezrelyho jako prusecik
    IntersectionPoint2 = [Podezrely];
    Podezrely=null; 
else %kdyz tam zadny neni
    
for j=1:dl % cyklus, ktery se opakuje tolikrat kolik je pruseciku v mnozine
    
    blabla = IntersectionPoint2(pom3,1); %vyber spravneho pruseciku
    
   
    if    (blabla+2+pom2)==xg 
      switch pom3
        case 1
        pomA=pomA+1; 
        case 3
        pomB=pomB+1; 
        case 5
        pomC=pomC+1; 
        case 7
        pomD=pomD+1; 
        end
        pom2=pom2+2;
        Podezrely=null;
  % kdyz je tento podezrely(xg je jeho x-va souradnice) priliz blizko pruseciku..
  %Pozn: jedna se o to, ze pro jeden prusecik, nalezne algoritmus vice vhodnych bodu, nam staci 1, ostatni nevybereme    
        
        
        
    else
      
       IntersectionPoint2 = [IntersectionPoint2;Podezrely];
       pom3=pom3+2;
       pom2=0;
       Podezrely=null;
    end
    end

end
end
end
end
xg=xg+2;
  end
 
  try
MBx1=IntersectionPoint2(1,1);
MBy1=IntersectionPoint2(2,1);
MBx2=IntersectionPoint2(3,1);
MBy2=IntersectionPoint2(4,1);

%MB1

xg =MBx1+pomA; 
eq23aG = (((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23A^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23A^2)==1;
GH1=round(vpa(solve(eq23aG,yg)));
d1G=sqrt((MBy1-GH1(1,1)).^2);
d2G=sqrt((MBy1-GH1(2,1)).^2);
if d1G<d2G
IntersectionPoint3 = [xg;GH1(1,1)];

end
if d1G>d2G
IntersectionPoint3 = [xg;GH1(2,1)];    
end
%MB2
xg =MBx2+pomB; 
eq23aG = (((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23A^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23A^2)==1;
GH1=round(vpa(solve(eq23aG,yg)));
d1G=sqrt((MBy2-GH1(1,1)).^2);
d2G=sqrt((MBy2-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
try
MBx3=IntersectionPoint2(5,1);
MBy3=IntersectionPoint2(6,1);
MBx4=IntersectionPoint2(7,1);
MBy4=IntersectionPoint2(8,1);
%MB3
xg =MBx3+pomC; 
eq23aG = (((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23A^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23A^2)==1;
GH1=round(vpa(solve(eq23aG,yg)));
d1G=sqrt((MBy3-GH1(1,1)).^2);
d2G=sqrt((MBy3-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
%MB4
xg =MBx4+pomD; 
eq23aG = (((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23A^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23A^2)==1;
GH1=round(vpa(solve(eq23aG,yg)));
d1G=sqrt((MBy4-GH1(1,1)).^2);
d2G=sqrt((MBy4-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
end
  end

 try 
xg1a=round(vpa(IntersectionPoint3(1,1)));
xg2a=round(vpa(IntersectionPoint3(3,1)));
xg3a=round(vpa(IntersectionPoint3(5,1)));
xg4a=round(vpa(IntersectionPoint3(7,1)));
yg1a=round(vpa(IntersectionPoint3(2,1)));
yg2a=round(vpa(IntersectionPoint3(4,1)));
yg3a=round(vpa(IntersectionPoint3(6,1)));
yg4a=round(vpa(IntersectionPoint3(8,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);
d3aIP=sqrt((ipx-xg3a).^2+(ipy-yg3a).^2);
d4aIP=sqrt((ipx-xg4a).^2+(ipy-yg4a).^2);

  if (d1aIP <d2aIP && d1aIP <d3aIP && d1aIP <d4aIP)
    Point12B23B = [xg1a, yg1a];
    d2 = d1aIP;
  end
  if (d2aIP <d1aIP && d2aIP <d3aIP && d2aIP <d4aIP)
    Point12B23B = [xg2a, yg2a];
     d2 = d2aIP;
  end
  if (d3aIP <d2aIP && d3aIP <d1aIP && d3aIP <d4aIP)
    Point12B23B = [xg3a, yg3a];
     d2 = d3aIP;
  end
  if (d4aIP <d2aIP && d4aIP <d3aIP && d4aIP <d1aIP)
    Point12B23B = [xg4a, yg4a];
     d2 = d4aIP;
  end
D= [D ; Point12B23B];
 catch
     try
xg1a=round(vpa(IntersectionPoint3(1,1)));
xg2a=round(vpa(IntersectionPoint3(3,1)));
yg1a=round(vpa(IntersectionPoint3(2,1)));
yg2a=round(vpa(IntersectionPoint3(4,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);


  if (d1aIP <d2aIP )
    Point12B23B = [xg1a, yg1a];
    d2 = d1aIP;
  end
  if (d2aIP <d1aIP)
    Point12B23B = [xg2a, yg2a];
     d2 = d2aIP;
  end
  
D= [D ; Point12B23B];
     catch
          MBx1=round(vpa(IntersectionPoint2(1,1)));
         MBy1=round(vpa(IntersectionPoint2(2,1)));
          
         
xg =MBx1+pomA; 
eq23aG = (((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23A^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23A^2)==1;
GH1=round(vpa(solve(eq23aG,yg)));
d1G=sqrt((MBy1-GH1(1,1)).^2);
d2G=sqrt((MBy1-GH1(2,1)).^2);
if d1G<d2G
Point12B23B= [xg,GH1(1,1)];
d2=sqrt((ipx-xg).^2+(ipy-GH1(1,1)).^2);

end
if d1G>d2G
Point12B23B = [xg,GH1(2,1)];   
d2=sqrt((ipx-xg).^2+(ipy-GH1(2,1)).^2);
end
         
        D= [D ;Point12B23B]; 
     end
 end
 syms xg;
end
disp('3)')

%3
[Xga,Yga] = solve(eq12b,eq13b);
xg1a=round(vpa(Xga(1,1)));
xg2a=round(vpa(Xga(2,1)));
xg3a=round(vpa(Xga(3,1)));
xg4a=round(vpa(Xga(4,1)));
yg1a=round(vpa(Yga(1,1)));
yg2a=round(vpa(Yga(2,1)));
yg3a=round(vpa(Yga(3,1)));
yg4a=round(vpa(Yga(4,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);
d3aIP=sqrt((ipx-xg3a).^2+(ipy-yg3a).^2);
d4aIP=sqrt((ipx-xg4a).^2+(ipy-yg4a).^2);


try
    
    if (xg1a >3000||xg1a <-3000||yg1a <-3000||yg1a <-3000 || xg2a >3000||xg2a <-3000||yg2a <-3000||yg2a <-3000)
        block 
    end
    
  if (d1aIP <d2aIP && d1aIP <d3aIP && d1aIP <d4aIP)
    Point12B23B = [xg1a, yg1a];
    d3 = d1aIP;
  end
  if (d2aIP <d1aIP && d2aIP <d3aIP && d2aIP <d4aIP)
    Point12B23B = [xg2a, yg2a];
     d3 = d2aIP;
  end
  if (d3aIP <d2aIP && d3aIP <d1aIP && d3aIP <d4aIP)
    Point12B23B = [xg3a, yg3a];
     d3 = d3aIP;
  end
  if (d4aIP <d2aIP && d4aIP <d3aIP && d4aIP <d1aIP)
    Point12B23B = [xg4a, yg4a];
     d3 = d4aIP;
  end
D= [D ; Point12B23B];
catch
  xg=Zacatek;
  IntersectionPoint2=[];
  pom2=0;
  pom3=1;
  IntersectionPoint3=[];
  pomA = 0;
  pomB = 0;
  pomC = 0;
  pomD = 0;  
  for i=1:Konec

eq12bG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12B^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12B^2)==1;
GH1=round(vpa(solve(eq12bG,yg)));

Point1=[xg;GH1(1,1)];
Point2=[xg;GH1(2,1)];
eq13bG = (((xg-m13)*cosq13+(yg-n13)*sinq13).^2)/a13B^2 -(((-(xg-m13)*sinq13+(yg-n13)*cosq13).^2)/b13B^2)==1;
GH2=round(vpa(solve(eq13bG,yg)));
Point3=[xg;GH2(1,1)];
Point4=[xg;GH2(2,1)];

d13=sqrt((xg-xg).^2+(GH1(1,1)-GH2(1,1)).^2);
d14=sqrt((xg-xg).^2+(GH1(1,1)-GH2(2,1)).^2);
d23=sqrt((xg-xg).^2+(GH1(2,1)-GH2(1,1)).^2);
d24=sqrt((xg-xg).^2+(GH1(2,1)-GH2(2,1)).^2);
try

if d13<Polomer
    Podezrely = [xg;GH1(1,1)];
end
if d14<Polomer
    Podezrely = [xg;GH1(1,1)];
end
if d23<Polomer
    Podezrely = [xg;GH1(2,1)];
end
if d24<Polomer
    Podezrely = [xg;GH1(2,1)];
end
if Podezrely(1,1)==xg %kdyz je nejaky Podezrely

try%zkus
    
[dm,dn] = size(IntersectionPoint2);%zjisti kolik je uz pruseciku v mnozine
dl=dm/2; %pomocny vypocet, pruseiku je v mnozine tolik kolik je radku deleno dvema, protoze kazdy podezrely ma na 1 radku x-ovou a na druhem y-ovou souradnici

if dl==0% kdyz v mnozine zadny jeste neni tak priradine prvniho podezrelyho jako prusecik
    IntersectionPoint2 = [Podezrely];
    Podezrely=null; 
else %kdyz tam zadny neni
    
for j=1:dl % cyklus, ktery se opakuje tolikrat kolik je pruseciku v mnozine
    
    blabla = IntersectionPoint2(pom3,1); %vyber spravneho pruseciku
    
   
    if    (blabla+2+pom2)==xg 
      switch pom3
        case 1
        pomA=pomA+1; 
        case 3
        pomB=pomB+1; 
        case 5
        pomC=pomC+1; 
        case 7
        pomD=pomD+1; 
       end
        pom2=pom2+2;
        Podezrely=null;
        
  % kdyz je tento podezrely(xg je jeho x-va souradnice) priliz blizko pruseciku..
  %Pozn: jedna se o to, ze pro jeden prusecik, nalezne algoritmus vice vhodnych bodu, nam staci 1, ostatni nevybereme    
        
        
        
    else
      
       IntersectionPoint2 = [IntersectionPoint2;Podezrely];
       pom3=pom3+2;
       pom2=0;
       Podezrely=null;
    end
    end

end
end
end
end
xg=xg+2;
  end
 
  try
MBx1=IntersectionPoint2(1,1);
MBy1=IntersectionPoint2(2,1);
MBx2=IntersectionPoint2(3,1);
MBy2=IntersectionPoint2(4,1);

%MB1
xg =MBx1+pomA; 
eq12bG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12B^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12B^2)==1;
GH1=round(vpa(solve(eq12bG,yg)));
d1G=sqrt((MBy1-GH1(1,1)).^2);
d2G=sqrt((MBy1-GH1(2,1)).^2);
if d1G<d2G
IntersectionPoint3 = [xg;GH1(1,1)];

end
if d1G>d2G
IntersectionPoint3 = [xg;GH1(2,1)];    
end
%MB2
xg =MBx2+pomB; 
eq12bG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12B^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12B^2)==1;
GH1=round(vpa(solve(eq12bG,yg)));
d1G=sqrt((MBy2-GH1(1,1)).^2);
d2G=sqrt((MBy2-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
try
MBx3=IntersectionPoint2(5,1);
MBy3=IntersectionPoint2(6,1);
MBx4=IntersectionPoint2(7,1);
MBy4=IntersectionPoint2(8,1);
%MB3
xg =MBx3+pomC; 
eq12bG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12B^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12B^2)==1;
GH1=round(vpa(solve(eq12bG,yg)));
d1G=sqrt((MBy3-GH1(1,1)).^2);
d2G=sqrt((MBy3-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
%MB4
xg =MBx4+pomD; 
eq12bG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12B^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12B^2)==1;
GH1=round(vpa(solve(eq12bG,yg)));
d1G=sqrt((MBy4-GH1(1,1)).^2);
d2G=sqrt((MBy4-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
end
  end

 try 
xg1a=round(vpa(IntersectionPoint3(1,1)));
xg2a=round(vpa(IntersectionPoint3(3,1)));
xg3a=round(vpa(IntersectionPoint3(5,1)));
xg4a=round(vpa(IntersectionPoint3(7,1)));
yg1a=round(vpa(IntersectionPoint3(2,1)));
yg2a=round(vpa(IntersectionPoint3(4,1)));
yg3a=round(vpa(IntersectionPoint3(6,1)));
yg4a=round(vpa(IntersectionPoint3(8,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);
d3aIP=sqrt((ipx-xg3a).^2+(ipy-yg3a).^2);
d4aIP=sqrt((ipx-xg4a).^2+(ipy-yg4a).^2);

  if (d1aIP <d2aIP && d1aIP <d3aIP && d1aIP <d4aIP)
    Point12B23B = [xg1a, yg1a];
    d3 = d1aIP;
  end
  if (d2aIP <d1aIP && d2aIP <d3aIP && d2aIP <d4aIP)
    Point12B23B = [xg2a, yg2a];
     d3 = d2aIP;
  end
  if (d3aIP <d2aIP && d3aIP <d1aIP && d3aIP <d4aIP)
    Point12B23B = [xg3a, yg3a];
     d3 = d3aIP;
  end
  if (d4aIP <d2aIP && d4aIP <d3aIP && d4aIP <d1aIP)
    Point12B23B = [xg4a, yg4a];
     d3 = d4aIP;
  end
D= [D ; Point12B23B];
 catch
     try
xg1a=round(vpa(IntersectionPoint3(1,1)));
xg2a=round(vpa(IntersectionPoint3(3,1)));
yg1a=round(vpa(IntersectionPoint3(2,1)));
yg2a=round(vpa(IntersectionPoint3(4,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);


  if (d1aIP <d2aIP )
    Point12B23B = [xg1a, yg1a];
    d3 = d1aIP;
  end
  if (d2aIP <d1aIP)
    Point12B23B = [xg2a, yg2a];
     d3 = d2aIP;
  end
  
D= [D ; Point12B23B];
catch
         MBx1=round(vpa(IntersectionPoint2(1,1)));
         MBy1=round(vpa(IntersectionPoint2(2,1)));
          
         
xg =MBx1+pomA; 
eq12bG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12B^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12B^2)==1;
GH1=round(vpa(solve(eq12bG,yg)));
d1G=sqrt((MBy1-GH1(1,1)).^2);
d2G=sqrt((MBy1-GH1(2,1)).^2);
if d1G<d2G
Point12B23B= [xg,GH1(1,1)];
d3=sqrt((ipx-xg).^2+(ipy-GH1(1,1)).^2);

end
if d1G>d2G
Point12B23B = [xg,GH1(2,1)];   
d3=sqrt((ipx-xg).^2+(ipy-GH1(2,1)).^2);
end
 
D= [D ; Point12B23B]; 

 end
 end
 syms xg;
end

%4
disp('4)')

[Xga,Yga] = solve(eq12a,eq13b);
xg1a=round(vpa(Xga(1,1)));
xg2a=round(vpa(Xga(2,1)));
xg3a=round(vpa(Xga(3,1)));
xg4a=round(vpa(Xga(4,1)));
yg1a=round(vpa(Yga(1,1)));
yg2a=round(vpa(Yga(2,1)));
yg3a=round(vpa(Yga(3,1)));
yg4a=round(vpa(Yga(4,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);
d3aIP=sqrt((ipx-xg3a).^2+(ipy-yg3a).^2);
d4aIP=sqrt((ipx-xg4a).^2+(ipy-yg4a).^2);

try
    if (xg1a >3000||xg1a <-3000||yg1a <-3000||yg1a <-3000 || xg2a >3000||xg2a <-3000||yg2a <-3000||yg2a <-3000)
        block 
    end
    
  if (d1aIP <d2aIP && d1aIP <d3aIP && d1aIP <d4aIP)
    Point12B23B = [xg1a, yg1a];
    d4 = d1aIP;
  end
  if (d2aIP <d1aIP && d2aIP <d3aIP && d2aIP <d4aIP)
    Point12B23B = [xg2a, yg2a];
     d4 = d2aIP;
  end
  if (d3aIP <d2aIP && d3aIP <d1aIP && d3aIP <d4aIP)
    Point12B23B = [xg3a, yg3a];
     d4 = d3aIP;
  end
  if (d4aIP <d2aIP && d4aIP <d3aIP && d4aIP <d1aIP)
    Point12B23B = [xg4a, yg4a];
     d4 = d4aIP;
  end
D= [D ; Point12B23B];
catch
  xg=Zacatek;
  IntersectionPoint2=[];
  pom2=0;
  pom3=1;
 IntersectionPoint3=[];
  pomA = 0;
  pomB = 0;
  pomC = 0;
  pomD = 0;   
  for i=1:Konec

eq12aG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12A^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12A^2)==1;
GH1=round(vpa(solve(eq12aG,yg)));

Point1=[xg;GH1(1,1)];
Point2=[xg;GH1(2,1)];
eq13bG = (((xg-m13)*cosq13+(yg-n13)*sinq13).^2)/a13B^2 -(((-(xg-m13)*sinq13+(yg-n13)*cosq13).^2)/b13B^2)==1;
GH2=round(vpa(solve(eq13bG,yg)));
Point3=[xg;GH2(1,1)];
Point4=[xg;GH2(2,1)];

d13=sqrt((xg-xg).^2+(GH1(1,1)-GH2(1,1)).^2);
d14=sqrt((xg-xg).^2+(GH1(1,1)-GH2(2,1)).^2);
d23=sqrt((xg-xg).^2+(GH1(2,1)-GH2(1,1)).^2);
d24=sqrt((xg-xg).^2+(GH1(2,1)-GH2(2,1)).^2);
try

if d13<Polomer
    Podezrely = [xg;GH1(1,1)];
end
if d14<Polomer
    Podezrely = [xg;GH1(1,1)];
end
if d23<Polomer
    Podezrely = [xg;GH1(2,1)];
end
if d24<Polomer
    Podezrely = [xg;GH1(2,1)];
end
if Podezrely(1,1)==xg %kdyz je nejaky Podezrely

try%zkus
    
[dm,dn] = size(IntersectionPoint2);%zjisti kolik je uz pruseciku v mnozine
dl=dm/2; %pomocny vypocet, pruseiku je v mnozine tolik kolik je radku deleno dvema, protoze kazdy podezrely ma na 1 radku x-ovou a na druhem y-ovou souradnici

if dl==0% kdyz v mnozine zadny jeste neni tak priradine prvniho podezrelyho jako prusecik
    IntersectionPoint2 = [Podezrely];
    Podezrely=null; 
else %kdyz tam zadny neni
    
for j=1:dl % cyklus, ktery se opakuje tolikrat kolik je pruseciku v mnozine
    
    blabla = IntersectionPoint2(pom3,1); %vyber spravneho pruseciku
    
   
    if    (blabla+2+pom2)==xg 
      switch pom3
        case 1
        pomA=pomA+1; 
        case 3
        pomB=pomB+1; 
        case 5
        pomC=pomC+1; 
        case 7
        pomD=pomD+1; 
        end
        pom2=pom2+2;
        Podezrely=null;
  % kdyz je tento podezrely(xg je jeho x-va souradnice) priliz blizko pruseciku..
  %Pozn: jedna se o to, ze pro jeden prusecik, nalezne algoritmus vice vhodnych bodu, nam staci 1, ostatni nevybereme    
        
        
        
    else
      
       IntersectionPoint2 = [IntersectionPoint2;Podezrely];
       pom3=pom3+2;
       pom2=0;
       Podezrely=null;
    end
    end

end
end
end
end
xg=xg+2;
  end
  try
MBx1=IntersectionPoint2(1,1);
MBy1=IntersectionPoint2(2,1);
MBx2=IntersectionPoint2(3,1);
MBy2=IntersectionPoint2(4,1);

%MB1
xg =MBx1+pomA; 
eq12aG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12A^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12A^2)==1;
GH1=round(vpa(solve(eq12aG,yg)));
d1G=sqrt((MBy1-GH1(1,1)).^2);
d2G=sqrt((MBy1-GH1(2,1)).^2);
if d1G<d2G
IntersectionPoint3 = [xg;GH1(1,1)];

end
if d1G>d2G
IntersectionPoint3 = [xg;GH1(2,1)];    
end
%MB2
xg =MBx2+pomB; 
eq12aG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12A^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12A^2)==1;
GH1=round(vpa(solve(eq12aG,yg)));
d1G=sqrt((MBy2-GH1(1,1)).^2);
d2G=sqrt((MBy2-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
try
MBx3=IntersectionPoint2(5,1);
MBy3=IntersectionPoint2(6,1);
MBx4=IntersectionPoint2(7,1);
MBy4=IntersectionPoint2(8,1);
%MB3
xg =MBx3+pomC; 
eq12aG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12A^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12A^2)==1;
GH1=round(vpa(solve(eq12aG,yg)));
d1G=sqrt((MBy3-GH1(1,1)).^2);
d2G=sqrt((MBy3-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
%MB4
xg =MBx4+pomD; 
eq12aG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12A^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12A^2)==1;
GH1=round(vpa(solve(eq12aG,yg)));
d1G=sqrt((MBy4-GH1(1,1)).^2);
d2G=sqrt((MBy4-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
end
end
 try 
xg1a=round(vpa(IntersectionPoint3(1,1)));
xg2a=round(vpa(IntersectionPoint3(3,1)));
xg3a=round(vpa(IntersectionPoint3(5,1)));
xg4a=round(vpa(IntersectionPoint3(7,1)));
yg1a=round(vpa(IntersectionPoint3(2,1)));
yg2a=round(vpa(IntersectionPoint3(4,1)));
yg3a=round(vpa(IntersectionPoint3(6,1)));
yg4a=round(vpa(IntersectionPoint3(8,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);
d3aIP=sqrt((ipx-xg3a).^2+(ipy-yg3a).^2);
d4aIP=sqrt((ipx-xg4a).^2+(ipy-yg4a).^2);

  if (d1aIP <d2aIP && d1aIP <d3aIP && d1aIP <d4aIP)
    Point12B23B = [xg1a, yg1a];
    d4 = d1aIP;
  end
  if (d2aIP <d1aIP && d2aIP <d3aIP && d2aIP <d4aIP)
    Point12B23B = [xg2a, yg2a];
     d4 = d2aIP;
  end
  if (d3aIP <d2aIP && d3aIP <d1aIP && d3aIP <d4aIP)
    Point12B23B = [xg3a, yg3a];
     d4 = d3aIP;
  end
  if (d4aIP <d2aIP && d4aIP <d3aIP && d4aIP <d1aIP)
    Point12B23B = [xg4a, yg4a];
     d4 = d4aIP;
  end
D= [D ; Point12B23B];
 catch
     try
xg1a=round(vpa(IntersectionPoint3(1,1)));
xg2a=round(vpa(IntersectionPoint3(3,1)));
yg1a=round(vpa(IntersectionPoint3(2,1)));
yg2a=round(vpa(IntersectionPoint3(4,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);


  if (d1aIP <d2aIP )
    Point12B23B = [xg1a, yg1a];
    d4 = d1aIP;
  end
  if (d2aIP <d1aIP)
    Point12B23B = [xg2a, yg2a];
     d4 = d2aIP;
  end
  
D= [D ; Point12B23B];
catch
          MBx1=round(vpa(IntersectionPoint2(1,1)));
         MBy1=round(vpa(IntersectionPoint2(2,1)));
          
         
xg =MBx1+pomA; 
eq12aG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12A^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12A^2)==1;
GH1=round(vpa(solve(eq12aG,yg)));
d1G=sqrt((MBy1-GH1(1,1)).^2);
d2G=sqrt((MBy1-GH1(2,1)).^2);
if d1G<d2G
Point12B23B= [xg,GH1(1,1)];
d4=sqrt((ipx-xg).^2+(ipy-GH1(1,1)).^2);

end
if d1G>d2G
Point12B23B = [xg,GH1(2,1)];   
d4=sqrt((ipx-xg).^2+(ipy-GH1(2,1)).^2);
end
         
        D= [D ;Point12B23B]; 
     end
 end
 syms xg;
end

disp('5)')
%5
[Xga,Yga] = solve(eq13a,eq23b);
xg1a=round(vpa(Xga(1,1)));
xg2a=round(vpa(Xga(2,1)));
xg3a=round(vpa(Xga(3,1)));
xg4a=round(vpa(Xga(4,1)));
yg1a=round(vpa(Yga(1,1)));
yg2a=round(vpa(Yga(2,1)));
yg3a=round(vpa(Yga(3,1)));
yg4a=round(vpa(Yga(4,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);
d3aIP=sqrt((ipx-xg3a).^2+(ipy-yg3a).^2);
d4aIP=sqrt((ipx-xg4a).^2+(ipy-yg4a).^2);


try
    
    if (xg1a >3000||xg1a <-3000||yg1a <-3000||yg1a <-3000 || xg2a >3000||xg2a <-3000||yg2a <-3000||yg2a <-3000)
        block 
    end
    
  if (d1aIP <d2aIP && d1aIP <d3aIP && d1aIP <d4aIP)
    Point12B23B = [xg1a, yg1a];
    d5 = d1aIP;
  end
  if (d2aIP <d1aIP && d2aIP <d3aIP && d2aIP <d4aIP)
    Point12B23B = [xg2a, yg2a];
     d5 = d2aIP;
  end
  if (d3aIP <d2aIP && d3aIP <d1aIP && d3aIP <d4aIP)
    Point12B23B = [xg3a, yg3a];
     d5 = d3aIP;
  end
  if (d4aIP <d2aIP && d4aIP <d3aIP && d4aIP <d1aIP)
    Point12B23B = [xg4a, yg4a];
     d5 = d4aIP;
  end
D= [D ; Point12B23B];
catch
  xg=Zacatek;
  IntersectionPoint2=[];
  pom2=0;
  pom3=1;
  IntersectionPoint3=[];
  pomA = 0;
  pomB = 0;
  pomC = 0;
  pomD = 0;  
  for i=1:Konec

eq13aG = (((xg-m13)*cosq13+(yg-n13)*sinq13).^2)/a13A^2 -(((-(xg-m13)*sinq13+(yg-n13)*cosq13).^2)/b13A^2)==1;
GH1=round(vpa(solve(eq13aG,yg)));

Point1=[xg;GH1(1,1)];
Point2=[xg;GH1(2,1)];
eq23bG=(((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23B^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23B^2)==1; 
GH2=round(vpa(solve(eq23bG,yg)));
Point3=[xg;GH2(1,1)];
Point4=[xg;GH2(2,1)];

d13=sqrt((xg-xg).^2+(GH1(1,1)-GH2(1,1)).^2);
d14=sqrt((xg-xg).^2+(GH1(1,1)-GH2(2,1)).^2);
d23=sqrt((xg-xg).^2+(GH1(2,1)-GH2(1,1)).^2);
d24=sqrt((xg-xg).^2+(GH1(2,1)-GH2(2,1)).^2);
try

if d13<Polomer
    Podezrely = [xg;GH1(1,1)];
end
if d14<Polomer
    Podezrely = [xg;GH1(1,1)];
end
if d23<Polomer
    Podezrely = [xg;GH1(2,1)];
end
if d24<Polomer
    Podezrely = [xg;GH1(2,1)];
end
if Podezrely(1,1)==xg %kdyz je nejaky Podezrely

try%zkus
    
[dm,dn] = size(IntersectionPoint2);%zjisti kolik je uz pruseciku v mnozine
dl=dm/2; %pomocny vypocet, pruseiku je v mnozine tolik kolik je radku deleno dvema, protoze kazdy podezrely ma na 1 radku x-ovou a na druhem y-ovou souradnici

if dl==0% kdyz v mnozine zadny jeste neni tak priradine prvniho podezrelyho jako prusecik
    IntersectionPoint2 = [Podezrely];
    Podezrely=null; 
else %kdyz tam zadny neni
    
for j=1:dl % cyklus, ktery se opakuje tolikrat kolik je pruseciku v mnozine
    
    blabla = IntersectionPoint2(pom3,1); %vyber spravneho pruseciku
    
   
    if    (blabla+2+pom2)==xg 
      switch pom3
        case 1
        pomA=pomA+1; 
        case 3
        pomB=pomB+1; 
        case 5
        pomC=pomC+1; 
        case 7
        pomD=pomD+1; 
        end
        pom2=pom2+2;
        Podezrely=null;
  % kdyz je tento podezrely(xg je jeho x-va souradnice) priliz blizko pruseciku..
  %Pozn: jedna se o to, ze pro jeden prusecik, nalezne algoritmus vice vhodnych bodu, nam staci 1, ostatni nevybereme    
        
        
        
    else
      
       IntersectionPoint2 = [IntersectionPoint2;Podezrely];
       pom3=pom3+2;
       pom2=0;
       Podezrely=null;
    end
    end

end
end
end
end
xg=xg+2;
  end
  IntersectionPoint2
  pomA
  pomB
  pomC
  pomD  
  try
MBx1=IntersectionPoint2(1,1);
MBy1=IntersectionPoint2(2,1);
MBx2=IntersectionPoint2(3,1);
MBy2=IntersectionPoint2(4,1);

%MB1
xg =MBx1+pomA; 
eq23bG=(((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23B^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23B^2)==1; 
GH1=round(vpa(solve(eq23bG,yg)));
d1G=sqrt((MBy1-GH1(1,1)).^2);
d2G=sqrt((MBy1-GH1(2,1)).^2);
if d1G<d2G
IntersectionPoint3 = [xg;GH1(1,1)];

end
if d1G>d2G
IntersectionPoint3 = [xg;GH1(2,1)];    
end
%MB2
xg =MBx2+pomB; 
eq23bG=(((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23B^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23B^2)==1; 
GH1=round(vpa(solve(eq23bG,yg)));
d1G=sqrt((MBy2-GH1(1,1)).^2);
d2G=sqrt((MBy2-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
try
MBx3=IntersectionPoint2(5,1);
MBy3=IntersectionPoint2(6,1);
MBx4=IntersectionPoint2(7,1);
MBy4=IntersectionPoint2(8,1);
%MB3
xg =MBx3+pomC; 
eq23bG=(((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23B^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23B^2)==1; 
GH1=round(vpa(solve(eq23bG,yg)));
d1G=sqrt((MBy3-GH1(1,1)).^2);
d2G=sqrt((MBy3-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
%MB4
xg =MBx4+pomD; 
eq23bG=(((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23B^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23B^2)==1; 
GH1=round(vpa(solve(eq23bG,yg)));
d1G=sqrt((MBy4-GH1(1,1)).^2);
d2G=sqrt((MBy4-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
end
end
 try 
xg1a=round(vpa(IntersectionPoint3(1,1)));
xg2a=round(vpa(IntersectionPoint3(3,1)));
xg3a=round(vpa(IntersectionPoint3(5,1)));
xg4a=round(vpa(IntersectionPoint3(7,1)));
yg1a=round(vpa(IntersectionPoint3(2,1)));
yg2a=round(vpa(IntersectionPoint3(4,1)));
yg3a=round(vpa(IntersectionPoint3(6,1)));
yg4a=round(vpa(IntersectionPoint3(8,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);
d3aIP=sqrt((ipx-xg3a).^2+(ipy-yg3a).^2);
d4aIP=sqrt((ipx-xg4a).^2+(ipy-yg4a).^2);

  if (d1aIP <d2aIP && d1aIP <d3aIP && d1aIP <d4aIP)
    Point12B23B = [xg1a, yg1a];
    d5 = d1aIP;
  end
  if (d2aIP <d1aIP && d2aIP <d3aIP && d2aIP <d4aIP)
    Point12B23B = [xg2a, yg2a];
     d5 = d2aIP;
  end
  if (d3aIP <d2aIP && d3aIP <d1aIP && d3aIP <d4aIP)
    Point12B23B = [xg3a, yg3a];
     d5 = d3aIP;
  end
  if (d4aIP <d2aIP && d4aIP <d3aIP && d4aIP <d1aIP)
    Point12B23B = [xg4a, yg4a];
     d5 = d4aIP;
  end
D= [D ; Point12B23B];
 catch
     try
xg1a=round(vpa(IntersectionPoint3(1,1)));
xg2a=round(vpa(IntersectionPoint3(3,1)));
yg1a=round(vpa(IntersectionPoint3(2,1)));
yg2a=round(vpa(IntersectionPoint3(4,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);


  if (d1aIP <d2aIP )
    Point12B23B = [xg1a, yg1a];
    d5 = d1aIP;
  end
  if (d2aIP <d1aIP)
    Point12B23B = [xg2a, yg2a];
     d5 = d2aIP;
  end
  
D= [D ; Point12B23B];
catch
                  MBx1=round(vpa(IntersectionPoint2(1,1)));
         MBy1=round(vpa(IntersectionPoint2(2,1)));
          
         
xg =MBx1+pomA; 
eq23bG=(((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23B^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23B^2)==1; 
GH1=round(vpa(solve(eq23bG,yg)));
d1G=sqrt((MBy1-GH1(1,1)).^2);
d2G=sqrt((MBy1-GH1(2,1)).^2);
if d1G<d2G
Point12B23B= [xg,GH1(1,1)];
d5=sqrt((ipx-xg).^2+(ipy-GH1(1,1)).^2);

end
if d1G>d2G
Point12B23B = [xg,GH1(2,1)];   
d5=sqrt((ipx-xg).^2+(ipy-GH1(2,1)).^2);
end
         
        D= [D ;Point12B23B]; 
     end
 end
 syms xg;
end

disp('6)')
%6

[Xga,Yga] = solve(eq12b,eq23b);

xg1a=round(vpa(Xga(1,1)));
xg2a=round(vpa(Xga(2,1)));
xg3a=round(vpa(Xga(3,1)));
xg4a=round(vpa(Xga(4,1)));
yg1a=round(vpa(Yga(1,1)));
yg2a=round(vpa(Yga(2,1)));
yg3a=round(vpa(Yga(3,1)));
yg4a=round(vpa(Yga(4,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);
d3aIP=sqrt((ipx-xg3a).^2+(ipy-yg3a).^2);
d4aIP=sqrt((ipx-xg4a).^2+(ipy-yg4a).^2);
try
    if (xg1a >3000||xg1a <-3000||yg1a <-3000||yg1a <-3000 || xg2a >3000||xg2a <-3000||yg2a <-3000||yg2a <-3000)
        block 
    end
  if (d1aIP <d2aIP && d1aIP <d3aIP && d1aIP <d4aIP)
    Point12B23B = [xg1a, yg1a];
    d6 = d1aIP;
  end
  if (d2aIP <d1aIP && d2aIP <d3aIP && d2aIP <d4aIP)
    Point12B23B = [xg2a, yg2a];
     d6 = d2aIP;
  end
  if (d3aIP <d2aIP && d3aIP <d1aIP && d3aIP <d4aIP)
    Point12B23B = [xg3a, yg3a];
     d6 = d3aIP;
  end
  if (d4aIP <d2aIP && d4aIP <d3aIP && d4aIP <d1aIP)
    Point12B23B = [xg4a, yg4a];
     d6 = d4aIP;
  end
D= [D ; Point12B23B];
catch
  xg=Zacatek;
  IntersectionPoint2=[];
  pom2=0;
  pom3=1;
  IntersectionPoint3=[];
  pomA = 0;
  pomB = 0;
  pomC = 0;
  pomD = 0;  
  for i=1:Konec

eq12bG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12B^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12B^2)==1;
GH1=round(vpa(solve(eq12bG,yg)));

Point1=[xg;GH1(1,1)];
Point2=[xg;GH1(2,1)];
eq23bG=(((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23B^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23B^2)==1; 
GH2=round(vpa(solve(eq23bG,yg)));
Point3=[xg;GH2(1,1)];
Point4=[xg;GH2(2,1)];

d13=sqrt((xg-xg).^2+(GH1(1,1)-GH2(1,1)).^2);
d14=sqrt((xg-xg).^2+(GH1(1,1)-GH2(2,1)).^2);
d23=sqrt((xg-xg).^2+(GH1(2,1)-GH2(1,1)).^2);
d24=sqrt((xg-xg).^2+(GH1(2,1)-GH2(2,1)).^2);
try

if d13<Polomer
    Podezrely = [xg;GH1(1,1)];
end
if d14<Polomer
    Podezrely = [xg;GH1(1,1)];
end
if d23<Polomer
    Podezrely = [xg;GH1(2,1)];
end
if d24<Polomer
    Podezrely = [xg;GH1(2,1)];
end
if Podezrely(1,1)==xg %kdyz je nejaky Podezrely

try%zkus
    
[dm,dn] = size(IntersectionPoint2);%zjisti kolik je uz pruseciku v mnozine
dl=dm/2; %pomocny vypocet, pruseiku je v mnozine tolik kolik je radku deleno dvema, protoze kazdy podezrely ma na 1 radku x-ovou a na druhem y-ovou souradnici

if dl==0% kdyz v mnozine zadny jeste neni tak priradine prvniho podezrelyho jako prusecik
    IntersectionPoint2 = [Podezrely];
    Podezrely=null; 
else %kdyz tam zadny neni
    
for j=1:dl % cyklus, ktery se opakuje tolikrat kolik je pruseciku v mnozine
    
    blabla = IntersectionPoint2(pom3,1); %vyber spravneho pruseciku
    
   
    if    (blabla+2+pom2)==xg 
      switch pom3
        case 1
        pomA=pomA+1; 
        case 3
        pomB=pomB+1; 
        case 5
        pomC=pomC+1; 
        case 7
        pomD=pomD+1; 
        end
        pom2=pom2+2;
        Podezrely=null;
  % kdyz je tento podezrely(xg je jeho x-va souradnice) priliz blizko pruseciku..
  %Pozn: jedna se o to, ze pro jeden prusecik, nalezne algoritmus vice vhodnych bodu, nam staci 1, ostatni nevybereme    
        
        
        
    else
      
       IntersectionPoint2 = [IntersectionPoint2;Podezrely];
       pom3=pom3+2;
       pom2=0;
       Podezrely=null;
    end
    end

end
end
end
end
xg=xg+2;
  end
  
 try
MBx1=IntersectionPoint2(1,1);
MBy1=IntersectionPoint2(2,1);
MBx2=IntersectionPoint2(3,1);
MBy2=IntersectionPoint2(4,1);

%MB1
xg =MBx1+pomA; 
eq23bG=(((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23B^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23B^2)==1; 
GH1=round(vpa(solve(eq23bG,yg)));
d1G=sqrt((MBy1-GH1(1,1)).^2);
d2G=sqrt((MBy1-GH1(2,1)).^2);
if d1G<d2G
IntersectionPoint3 = [xg;GH1(1,1)];

end
if d1G>d2G
IntersectionPoint3 = [xg;GH1(2,1)];    
end
%MB2
xg =MBx2+pomB; 
eq23bG=(((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23B^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23B^2)==1; 
GH1=round(vpa(solve(eq23bG,yg)));
d1G=sqrt((MBy2-GH1(1,1)).^2);
d2G=sqrt((MBy2-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
try
MBx3=IntersectionPoint2(5,1);
MBy3=IntersectionPoint2(6,1);
MBx4=IntersectionPoint2(7,1);
MBy4=IntersectionPoint2(8,1);
%MB3
xg =MBx3+pomC; 
eq23bG=(((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23B^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23B^2)==1; 
GH1=round(vpa(solve(eq23bG,yg)));
d1G=sqrt((MBy3-GH1(1,1)).^2);
d2G=sqrt((MBy3-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
%MB4
xg =MBx4+pomD; 
eq23bG=(((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23B^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23B^2)==1; 
GH1=round(vpa(solve(eq23bG,yg)));
d1G=sqrt((MBy4-GH1(1,1)).^2);
d2G=sqrt((MBy4-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
end
end
 try 
xg1a=round(vpa(IntersectionPoint3(1,1)));
xg2a=round(vpa(IntersectionPoint3(3,1)));
xg3a=round(vpa(IntersectionPoint3(5,1)));
xg4a=round(vpa(IntersectionPoint3(7,1)));
yg1a=round(vpa(IntersectionPoint3(2,1)));
yg2a=round(vpa(IntersectionPoint3(4,1)));
yg3a=round(vpa(IntersectionPoint3(6,1)));
yg4a=round(vpa(IntersectionPoint3(8,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);
d3aIP=sqrt((ipx-xg3a).^2+(ipy-yg3a).^2);
d4aIP=sqrt((ipx-xg4a).^2+(ipy-yg4a).^2);

  if (d1aIP <d2aIP && d1aIP <d3aIP && d1aIP <d4aIP)
    Point12B23B = [xg1a, yg1a];
    d6 = d1aIP;
  end
  if (d2aIP <d1aIP && d2aIP <d3aIP && d2aIP <d4aIP)
    Point12B23B = [xg2a, yg2a];
     d6 = d2aIP;
  end
  if (d3aIP <d2aIP && d3aIP <d1aIP && d3aIP <d4aIP)
    Point12B23B = [xg3a, yg3a];
     d6 = d3aIP;
  end
  if (d4aIP <d2aIP && d4aIP <d3aIP && d4aIP <d1aIP)
    Point12B23B = [xg4a, yg4a];
     d6 = d4aIP;
  end
D= [D ; Point12B23B];
 catch
     try
xg1a=round(vpa(IntersectionPoint3(1,1)));
xg2a=round(vpa(IntersectionPoint3(3,1)));
yg1a=round(vpa(IntersectionPoint3(2,1)));
yg2a=round(vpa(IntersectionPoint3(4,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);


  if (d1aIP <d2aIP )
    Point12B23B = [xg1a, yg1a];
    d6 = d1aIP;
  end
  if (d2aIP <d1aIP)
    Point12B23B = [xg2a, yg2a];
     d6 = d2aIP;
  end
  
D= [D ; Point12B23B];
catch
MBx1=round(vpa(IntersectionPoint2(1,1)));
MBy1=round(vpa(IntersectionPoint2(2,1)));
          
         
xg =MBx1+pomA; 
eq23bG=(((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23B^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23B^2)==1; 
GH1=round(vpa(solve(eq23bG,yg)));
d1G=sqrt((MBy1-GH1(1,1)).^2);
d2G=sqrt((MBy1-GH1(2,1)).^2);
if d1G<d2G
Point12B23B= [xg,GH1(1,1)];
d6=sqrt((ipx-xg).^2+(ipy-GH1(1,1)).^2);

end
if d1G>d2G
Point12B23B = [xg,GH1(2,1)];   
d6=sqrt((ipx-xg).^2+(ipy-GH1(2,1)).^2);
end  
        D= [D ;Point12B23B]; 
     end
 end
 syms xg;
end

disp('7)')
%7
[Xga,Yga] = solve(eq12a,eq23b);
xg1a=round(vpa(Xga(1,1)));
xg2a=round(vpa(Xga(2,1)));
xg3a=round(vpa(Xga(3,1)));
xg4a=round(vpa(Xga(4,1)));
yg1a=round(vpa(Yga(1,1)));
yg2a=round(vpa(Yga(2,1)));
yg3a=round(vpa(Yga(3,1)));
yg4a=round(vpa(Yga(4,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);
d3aIP=sqrt((ipx-xg3a).^2+(ipy-yg3a).^2);
d4aIP=sqrt((ipx-xg4a).^2+(ipy-yg4a).^2);

try
    if (xg1a >3000||xg1a <-3000||yg1a <-3000||yg1a <-3000 || xg2a >3000||xg2a <-3000||yg2a <-3000||yg2a <-3000)
        block 
    end
    
  if (d1aIP <d2aIP && d1aIP <d3aIP && d1aIP <d4aIP)
    Point12B23B = [xg1a, yg1a];
    d7 = d1aIP;
  end
  if (d2aIP <d1aIP && d2aIP <d3aIP && d2aIP <d4aIP)
    Point12B23B = [xg2a, yg2a];
     d7 = d2aIP;
  end
  if (d3aIP <d2aIP && d3aIP <d1aIP && d3aIP <d4aIP)
    Point12B23B = [xg3a, yg3a];
     d7 = d3aIP;
  end
  if (d4aIP <d2aIP && d4aIP <d3aIP && d4aIP <d1aIP)
    Point12B23B = [xg4a, yg4a];
     d7 = d4aIP;
  end
D= [D ; Point12B23B];
catch
  xg=Zacatek;
  IntersectionPoint2=[];
  pom2=0;
  pom3=1;
  IntersectionPoint3=[];
  pomA = 0;
  pomB = 0;
  pomC = 0;
  pomD = 0;  
  for i=1:Konec

eq12aG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12A^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12A^2)==1;
GH1=round(vpa(solve(eq12aG,yg)));

Point1=[xg;GH1(1,1)];
Point2=[xg;GH1(2,1)];
eq23bG=(((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23B^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23B^2)==1; 
GH2=round(vpa(solve(eq23bG,yg)));
Point3=[xg;GH2(1,1)];
Point4=[xg;GH2(2,1)];

d13=sqrt((xg-xg).^2+(GH1(1,1)-GH2(1,1)).^2);
d14=sqrt((xg-xg).^2+(GH1(1,1)-GH2(2,1)).^2);
d23=sqrt((xg-xg).^2+(GH1(2,1)-GH2(1,1)).^2);
d24=sqrt((xg-xg).^2+(GH1(2,1)-GH2(2,1)).^2);
try

if d13<Polomer
    Podezrely = [xg;GH1(1,1)];
end
if d14<Polomer
    Podezrely = [xg;GH1(1,1)];
end
if d23<Polomer
    Podezrely = [xg;GH1(2,1)];
end
if d24<Polomer
    Podezrely = [xg;GH1(2,1)];
end
if Podezrely(1,1)==xg %kdyz je nejaky Podezrely

try %zkus
    
[dm,dn] = size(IntersectionPoint2);%zjisti kolik je uz pruseciku v mnozine
dl=dm/2; %pomocny vypocet, pruseiku je v mnozine tolik kolik je radku deleno dvema, protoze kazdy podezrely ma na 1 radku x-ovou a na druhem y-ovou souradnici

if dl==0% kdyz v mnozine zadny jeste neni tak priradine prvniho podezrelyho jako prusecik
    IntersectionPoint2 = [Podezrely];
    Podezrely=null; 
else %kdyz tam zadny neni
    
for j=1:dl % cyklus, ktery se opakuje tolikrat kolik je pruseciku v mnozine
    
    blabla = IntersectionPoint2(pom3,1); %vyber spravneho pruseciku
    
   
    if    (blabla+2+pom2)==xg 
      switch pom3
        case 1
        pomA=pomA+1; 
        case 3
        pomB=pomB+1; 
        case 5
        pomC=pomC+1; 
        case 7
        pomD=pomD+1; 
        end
        pom2=pom2+2;
        Podezrely=null;
  % kdyz je tento podezrely(xg je jeho x-va souradnice) priliz blizko pruseciku..
  %Pozn: jedna se o to, ze pro jeden prusecik, nalezne algoritmus vice vhodnych bodu, nam staci 1, ostatni nevybereme    
        
        
        
    else
      
       IntersectionPoint2 = [IntersectionPoint2;Podezrely];
       pom3=pom3+2;
       pom2=0;
       Podezrely=null;
    end
    end

end
end
end
end
xg=xg+2;
  end
  try
MBx1=IntersectionPoint2(1,1);
MBy1=IntersectionPoint2(2,1);
MBx2=IntersectionPoint2(3,1);
MBy2=IntersectionPoint2(4,1);

%MB1
xg =MBx1+pomA; 
eq23bG=(((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23B^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23B^2)==1; 
GH1=round(vpa(solve(eq23bG,yg)));
d1G=sqrt((MBy1-GH1(1,1)).^2);
d2G=sqrt((MBy1-GH1(2,1)).^2);
if d1G<d2G
IntersectionPoint3 = [xg;GH1(1,1)];

end
if d1G>d2G
IntersectionPoint3 = [xg;GH1(2,1)];    
end
%MB2
xg =MBx2+pomB; 
eq23bG=(((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23B^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23B^2)==1; 
GH1=round(vpa(solve(eq23bG,yg)));
d1G=sqrt((MBy2-GH1(1,1)).^2);
d2G=sqrt((MBy2-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
try
MBx3=IntersectionPoint2(5,1);
MBy3=IntersectionPoint2(6,1);
MBx4=IntersectionPoint2(7,1);
MBy4=IntersectionPoint2(8,1);
%MB3
xg =MBx3+pomC; 
eq23bG=(((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23B^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23B^2)==1; 
GH1=round(vpa(solve(eq23bG,yg)));
d1G=sqrt((MBy3-GH1(1,1)).^2);
d2G=sqrt((MBy3-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
%MB4
xg =MBx4+pomD; 
eq23bG=(((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23B^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23B^2)==1; 
GH1=round(vpa(solve(eq23bG,yg)));
d1G=sqrt((MBy4-GH1(1,1)).^2);
d2G=sqrt((MBy4-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
end
end
 try 
xg1a=round(vpa(IntersectionPoint3(1,1)));
xg2a=round(vpa(IntersectionPoint3(3,1)));
xg3a=round(vpa(IntersectionPoint3(5,1)));
xg4a=round(vpa(IntersectionPoint3(7,1)));
yg1a=round(vpa(IntersectionPoint3(2,1)));
yg2a=round(vpa(IntersectionPoint3(4,1)));
yg3a=round(vpa(IntersectionPoint3(6,1)));
yg4a=round(vpa(IntersectionPoint3(8,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);
d3aIP=sqrt((ipx-xg3a).^2+(ipy-yg3a).^2);
d4aIP=sqrt((ipx-xg4a).^2+(ipy-yg4a).^2);

  if (d1aIP <d2aIP && d1aIP <d3aIP && d1aIP <d4aIP)
    Point12B23B = [xg1a, yg1a];
    d7 = d1aIP;
  end
  if (d2aIP <d1aIP && d2aIP <d3aIP && d2aIP <d4aIP)
    Point12B23B = [xg2a, yg2a];
     d7 = d2aIP;
  end
  if (d3aIP <d2aIP && d3aIP <d1aIP && d3aIP <d4aIP)
    Point12B23B = [xg3a, yg3a];
     d7 = d3aIP;
  end
  if (d4aIP <d2aIP && d4aIP <d3aIP && d4aIP <d1aIP)
    Point12B23B = [xg4a, yg4a];
     d7 = d4aIP;
  end
D= [D ; Point12B23B];
 catch
     try
xg1a=round(vpa(IntersectionPoint3(1,1)));
xg2a=round(vpa(IntersectionPoint3(3,1)));
yg1a=round(vpa(IntersectionPoint3(2,1)));
yg2a=round(vpa(IntersectionPoint3(4,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);


  if (d1aIP <d2aIP )
    Point12B23B = [xg1a, yg1a];
    d7 = d1aIP;
  end
  if (d2aIP <d1aIP)
    Point12B23B = [xg2a, yg2a];
     d7 = d2aIP;
  end
  
D= [D ; Point12B23B];
catch
MBx1=round(vpa(IntersectionPoint2(1,1)));
MBy1=round(vpa(IntersectionPoint2(2,1)));
          
         
xg =MBx1+pomA; 
eq23bG=(((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23B^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23B^2)==1; 
GH1=round(vpa(solve(eq23bG,yg)));
d1G=sqrt((MBy1-GH1(1,1)).^2);
d2G=sqrt((MBy1-GH1(2,1)).^2);
if d1G<d2G
Point12B23B= [xg,GH1(1,1)];
d7=sqrt((ipx-xg).^2+(ipy-GH1(1,1)).^2);

end
if d1G>d2G
Point12B23B = [xg,GH1(2,1)];   
d7=sqrt((ipx-xg).^2+(ipy-GH1(2,1)).^2);
end
         
         
        D= [D ;Point12B23B]; 
     end
 end
 syms xg;
end
disp('8)')
%8
[Xga,Yga] = solve(eq12a,eq13a);
xg1a=round(vpa(Xga(1,1)));
xg2a=round(vpa(Xga(2,1)));
xg3a=round(vpa(Xga(3,1)));
xg4a=round(vpa(Xga(4,1)));
yg1a=round(vpa(Yga(1,1)));
yg2a=round(vpa(Yga(2,1)));
yg3a=round(vpa(Yga(3,1)));
yg4a=round(vpa(Yga(4,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);
d3aIP=sqrt((ipx-xg3a).^2+(ipy-yg3a).^2);
d4aIP=sqrt((ipx-xg4a).^2+(ipy-yg4a).^2);

try
    if (xg1a >3000||xg1a <-3000||yg1a <-3000||yg1a <-3000 || xg2a >3000||xg2a <-3000||yg2a <-3000||yg2a <-3000)
        block 
    end
    
  if (d1aIP <d2aIP && d1aIP <d3aIP && d1aIP <d4aIP)
    Point12B23B = [xg1a, yg1a];
    d8 = d1aIP;
  end
  if (d2aIP <d1aIP && d2aIP <d3aIP && d2aIP <d4aIP)
    Point12B23B = [xg2a, yg2a];
     d8 = d2aIP;
  end
  if (d3aIP <d2aIP && d3aIP <d1aIP && d3aIP <d4aIP)
    Point12B23B = [xg3a, yg3a];
     d8 = d3aIP;
  end
  if (d4aIP <d2aIP && d4aIP <d3aIP && d4aIP <d1aIP)
    Point12B23B = [xg4a, yg4a];
     d8 = d4aIP;
  end
D= [D ; Point12B23B];
catch
  xg=Zacatek;
  IntersectionPoint2=[];
  pom2=0;
  pom3=1;
  IntersectionPoint3=[];
  pomA = 0;
  pomB = 0;
  pomC = 0;
  pomD = 0;  
  for i=1:Konec

eq12aG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12A^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12A^2)==1;
GH1=round(vpa(solve(eq12aG,yg)));

Point1=[xg;GH1(1,1)];
Point2=[xg;GH1(2,1)];
eq13aG = (((xg-m13)*cosq13+(yg-n13)*sinq13).^2)/a13A^2 -(((-(xg-m13)*sinq13+(yg-n13)*cosq13).^2)/b13A^2)==1;
GH2=round(vpa(solve(eq13aG,yg)));
Point3=[xg;GH2(1,1)];
Point4=[xg;GH2(2,1)];

d13=sqrt((xg-xg).^2+(GH1(1,1)-GH2(1,1)).^2);
d14=sqrt((xg-xg).^2+(GH1(1,1)-GH2(2,1)).^2);
d23=sqrt((xg-xg).^2+(GH1(2,1)-GH2(1,1)).^2);
d24=sqrt((xg-xg).^2+(GH1(2,1)-GH2(2,1)).^2);
try

if d13<Polomer
    Podezrely = [xg;GH1(1,1)];
end
if d14<Polomer
    Podezrely = [xg;GH1(1,1)];
end
if d23<Polomer
    Podezrely = [xg;GH1(2,1)];
end
if d24<Polomer
    Podezrely = [xg;GH1(2,1)];
end
if Podezrely(1,1)==xg %kdyz je nejaky Podezrely

try %zkus
    
[dm,dn] = size(IntersectionPoint2);%zjisti kolik je uz pruseciku v mnozine
dl=dm/2; %pomocny vypocet, pruseiku je v mnozine tolik kolik je radku deleno dvema, protoze kazdy podezrely ma na 1 radku x-ovou a na druhem y-ovou souradnici

if dl==0% kdyz v mnozine zadny jeste neni tak priradine prvniho podezrelyho jako prusecik
    IntersectionPoint2 = [Podezrely];
    Podezrely=null; 
else %kdyz tam zadny neni
    
for j=1:dl % cyklus, ktery se opakuje tolikrat kolik je pruseciku v mnozine
    
    blabla = IntersectionPoint2(pom3,1); %vyber spravneho pruseciku
    
   
    if    (blabla+2+pom2)==xg 
      switch pom3
        case 1
        pomA=pomA+1; 
        case 3
        pomB=pomB+1; 
        case 5
        pomC=pomC+1; 
        case 7
        pomD=pomD+1; 
        end
        pom2=pom2+2;
        Podezrely=null;
  % kdyz je tento podezrely(xg je jeho x-va souradnice) priliz blizko pruseciku..
  %Pozn: jedna se o to, ze pro jeden prusecik, nalezne algoritmus vice vhodnych bodu, nam staci 1, ostatni nevybereme    
        
        
        
    else
      
       IntersectionPoint2 = [IntersectionPoint2;Podezrely];
       pom3=pom3+2;
       pom2=0;
       Podezrely=null;
    end
    end

end
end
end
end
xg=xg+2;
  end
  try
MBx1=IntersectionPoint2(1,1);
MBy1=IntersectionPoint2(2,1);
MBx2=IntersectionPoint2(3,1);
MBy2=IntersectionPoint2(4,1);

%MB1
xg =MBx1+pomA; 
eq12aG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12A^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12A^2)==1;
GH1=round(vpa(solve(eq12aG,yg)));
d1G=sqrt((MBy1-GH1(1,1)).^2);
d2G=sqrt((MBy1-GH1(2,1)).^2);
if d1G<d2G
IntersectionPoint3 = [xg;GH1(1,1)];

end
if d1G>d2G
IntersectionPoint3 = [xg;GH1(2,1)];    
end
%MB2
xg =MBx2+pomB; 
eq12aG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12A^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12A^2)==1;
GH1=round(vpa(solve(eq12aG,yg)));
d1G=sqrt((MBy2-GH1(1,1)).^2);
d2G=sqrt((MBy2-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
try
MBx3=IntersectionPoint2(5,1);
MBy3=IntersectionPoint2(6,1);
MBx4=IntersectionPoint2(7,1);
MBy4=IntersectionPoint2(8,1);
%MB3
xg =MBx3+pomC; 
eq12aG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12A^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12A^2)==1;
GH1=round(vpa(solve(eq12aG,yg)));
d1G=sqrt((MBy3-GH1(1,1)).^2);
d2G=sqrt((MBy3-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
%MB4
xg =MBx4+pomD; 
eq12aG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12A^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12A^2)==1;
GH1=round(vpa(solve(eq12aG,yg)));
d1G=sqrt((MBy4-GH1(1,1)).^2);
d2G=sqrt((MBy4-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
end
end
  
  
 try 
xg1a=round(vpa(IntersectionPoint3(1,1)));
xg2a=round(vpa(IntersectionPoint3(3,1)));
xg3a=round(vpa(IntersectionPoint3(5,1)));
xg4a=round(vpa(IntersectionPoint3(7,1)));
yg1a=round(vpa(IntersectionPoint3(2,1)));
yg2a=round(vpa(IntersectionPoint3(4,1)));
yg3a=round(vpa(IntersectionPoint3(6,1)));
yg4a=round(vpa(IntersectionPoint3(8,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);
d3aIP=sqrt((ipx-xg3a).^2+(ipy-yg3a).^2);
d4aIP=sqrt((ipx-xg4a).^2+(ipy-yg4a).^2);

  if (d1aIP <d2aIP && d1aIP <d3aIP && d1aIP <d4aIP)
    Point12B23B = [xg1a, yg1a];
    d8 = d1aIP;
  end
  if (d2aIP <d1aIP && d2aIP <d3aIP && d2aIP <d4aIP)
    Point12B23B = [xg2a, yg2a];
     d8 = d2aIP;
  end
  if (d3aIP <d2aIP && d3aIP <d1aIP && d3aIP <d4aIP)
    Point12B23B = [xg3a, yg3a];
     d8 = d3aIP;
  end
  if (d4aIP <d2aIP && d4aIP <d3aIP && d4aIP <d1aIP)
    Point12B23B = [xg4a, yg4a];
     d8 = d4aIP;
  end
D= [D ; Point12B23B];
 catch
     try
xg1a=round(vpa(IntersectionPoint3(1,1)));
xg2a=round(vpa(IntersectionPoint3(3,1)));
yg1a=round(vpa(IntersectionPoint3(2,1)));
yg2a=round(vpa(IntersectionPoint3(4,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);


  if (d1aIP <d2aIP )
    Point12B23B = [xg1a, yg1a];
    d8 = d1aIP;
  end
  if (d2aIP <d1aIP)
    Point12B23B = [xg2a, yg2a];
     d8 = d2aIP;
  end
  
D= [D ; Point12B23B];
catch
MBx1=round(vpa(IntersectionPoint2(1,1)));
MBy1=round(vpa(IntersectionPoint2(2,1)));
          
         
xg =MBx1+pomA; 
eq12aG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12A^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12A^2)==1;
GH1=round(vpa(solve(eq12aG,yg)));
d1G=sqrt((MBy1-GH1(1,1)).^2);
d2G=sqrt((MBy1-GH1(2,1)).^2);
if d1G<d2G
Point12B23B= [xg,GH1(1,1)];
d8=sqrt((ipx-xg).^2+(ipy-GH1(1,1)).^2);

end
if d1G>d2G
Point12B23B = [xg,GH1(2,1)];   
d8=sqrt((ipx-xg).^2+(ipy-GH1(2,1)).^2);
end         
        D= [D ;Point12B23B]; 
     end
 end
 syms xg;
end

%9
disp('9)')
[Xga,Yga] = solve(eq12b,eq13a);
xg1a=round(vpa(Xga(1,1)));
xg2a=round(vpa(Xga(2,1)));
xg3a=round(vpa(Xga(3,1)));
xg4a=round(vpa(Xga(4,1)));
yg1a=round(vpa(Yga(1,1)));
yg2a=round(vpa(Yga(2,1)));
yg3a=round(vpa(Yga(3,1)));
yg4a=round(vpa(Yga(4,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);
d3aIP=sqrt((ipx-xg3a).^2+(ipy-yg3a).^2);
d4aIP=sqrt((ipx-xg4a).^2+(ipy-yg4a).^2);

try
    if (xg1a >3000||xg1a <-3000||yg1a <-3000||yg1a <-3000 || xg2a >3000||xg2a <-3000||yg2a <-3000||yg2a <-3000)
        block 
    end
    
  if (d1aIP <d2aIP && d1aIP <d3aIP && d1aIP <d4aIP)
    Point12B23B = [xg1a, yg1a];
    d9 = d1aIP;
  end
  if (d2aIP <d1aIP && d2aIP <d3aIP && d2aIP <d4aIP)
    Point12B23B = [xg2a, yg2a];
     d9= d2aIP;
  end
  if (d3aIP <d2aIP && d3aIP <d1aIP && d3aIP <d4aIP)
    Point12B23B = [xg3a, yg3a];
     d9 = d3aIP;
  end
  if (d4aIP <d2aIP && d4aIP <d3aIP && d4aIP <d1aIP)
    Point12B23B = [xg4a, yg4a];
     d9 = d4aIP;
  end
D= [D ; Point12B23B];
catch
  xg=Zacatek;
  IntersectionPoint2=[];
  pom2=0;
  pom3=1;
  IntersectionPoint3=[];
  pomA = 0;
  pomB = 0;
  pomC = 0;
  pomD = 0;  
  for i=1:Konec

eq13aG = (((xg-m13)*cosq13+(yg-n13)*sinq13).^2)/a13A^2 -(((-(xg-m13)*sinq13+(yg-n13)*cosq13).^2)/b13A^2)==1;
GH1=round(vpa(solve(eq13aG,yg)));

Point1=[xg;GH1(1,1)];
Point2=[xg;GH1(2,1)];
eq12bG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12B^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12B^2)==1;
GH2=round(vpa(solve(eq12bG,yg)));
Point3=[xg;GH2(1,1)];
Point4=[xg;GH2(2,1)];

d13=sqrt((xg-xg).^2+(GH1(1,1)-GH2(1,1)).^2);
d14=sqrt((xg-xg).^2+(GH1(1,1)-GH2(2,1)).^2);
d23=sqrt((xg-xg).^2+(GH1(2,1)-GH2(1,1)).^2);
d24=sqrt((xg-xg).^2+(GH1(2,1)-GH2(2,1)).^2);
try

if d13<Polomer
    Podezrely = [xg;GH1(1,1)];
end
if d14<Polomer
    Podezrely = [xg;GH1(1,1)];
end
if d23<Polomer
    Podezrely = [xg;GH1(2,1)];
end
if d24<Polomer
    Podezrely = [xg;GH1(2,1)];
end
if Podezrely(1,1)==xg %kdyz je nejaky Podezrely

try %zkus
    
[dm,dn] = size(IntersectionPoint2);%zjisti kolik je uz pruseciku v mnozine
dl=dm/2; %pomocny vypocet, pruseiku je v mnozine tolik kolik je radku deleno dvema, protoze kazdy podezrely ma na 1 radku x-ovou a na druhem y-ovou souradnici

if dl==0% kdyz v mnozine zadny jeste neni tak priradine prvniho podezrelyho jako prusecik
    IntersectionPoint2 = [Podezrely];
    Podezrely=null; 
else %kdyz tam zadny neni
    
for j=1:dl % cyklus, ktery se opakuje tolikrat kolik je pruseciku v mnozine
    
    blabla = IntersectionPoint2(pom3,1); %vyber spravneho pruseciku
    
   
    if    (blabla+2+pom2)==xg 
      switch pom3
        case 1
        pomA=pomA+1; 
        case 3
        pomB=pomB+1; 
        case 5
        pomC=pomC+1; 
        case 7
        pomD=pomD+1; 
        end
        pom2=pom2+2;
        Podezrely=null;
  % kdyz je tento podezrely(xg je jeho x-va souradnice) priliz blizko pruseciku..
  %Pozn: jedna se o to, ze pro jeden prusecik, nalezne algoritmus vice vhodnych bodu, nam staci 1, ostatni nevybereme    
        
        
        
    else
      
       IntersectionPoint2 = [IntersectionPoint2;Podezrely];
       pom3=pom3+2;
       pom2=0;
       Podezrely=null;
    end
    end

end
end
end
end
xg=xg+2;
  end
IntersectionPoint2 
pomA 
  pomB 
  pomC 
  pomD 
  try
MBx1=IntersectionPoint2(1,1);
MBy1=IntersectionPoint2(2,1);
MBx2=IntersectionPoint2(3,1);
MBy2=IntersectionPoint2(4,1);

%MB1
xg =MBx1+pomA; 
eq12bG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12B^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12B^2)==1;
GH1=round(vpa(solve(eq12bG,yg)));
d1G=sqrt((MBy1-GH1(1,1)).^2);
d2G=sqrt((MBy1-GH1(2,1)).^2);
if d1G<d2G
IntersectionPoint3 = [xg;GH1(1,1)];

end
if d1G>d2G
IntersectionPoint3 = [xg;GH1(2,1)];    
end
%MB2
xg =MBx2+pomB; 
eq12bG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12B^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12B^2)==1;
GH1=round(vpa(solve(eq12bG,yg)));

d1G=sqrt((MBy2-GH1(1,1)).^2);
d2G=sqrt((MBy2-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
try
MBx3=IntersectionPoint2(5,1);
MBy3=IntersectionPoint2(6,1);
MBx4=IntersectionPoint2(7,1);
MBy4=IntersectionPoint2(8,1);
%MB3
xg =MBx3+pomC; 
eq12bG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12B^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12B^2)==1;
GH1=round(vpa(solve(eq12bG,yg)));

d1G=sqrt((MBy3-GH1(1,1)).^2);
d2G=sqrt((MBy3-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
%MB4
xg =MBx4+pomD; 
eq12bG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12B^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12B^2)==1;
GH1=round(vpa(solve(eq12bG,yg)));

d1G=sqrt((MBy4-GH1(1,1)).^2);
d2G=sqrt((MBy4-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
end
end
 try 
xg1a=round(vpa(IntersectionPoint3(1,1)));
xg2a=round(vpa(IntersectionPoint3(3,1)));
xg3a=round(vpa(IntersectionPoint3(5,1)));
xg4a=round(vpa(IntersectionPoint3(7,1)));
yg1a=round(vpa(IntersectionPoint3(2,1)));
yg2a=round(vpa(IntersectionPoint3(4,1)));
yg3a=round(vpa(IntersectionPoint3(6,1)));
yg4a=round(vpa(IntersectionPoint3(8,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);
d3aIP=sqrt((ipx-xg3a).^2+(ipy-yg3a).^2);
d4aIP=sqrt((ipx-xg4a).^2+(ipy-yg4a).^2);

  if (d1aIP <d2aIP && d1aIP <d3aIP && d1aIP <d4aIP)
    Point12B23B = [xg1a, yg1a];
    d9 = d1aIP;
  end
  if (d2aIP <d1aIP && d2aIP <d3aIP && d2aIP <d4aIP)
    Point12B23B = [xg2a, yg2a];
     d9 = d2aIP;
  end
  if (d3aIP <d2aIP && d3aIP <d1aIP && d3aIP <d4aIP)
    Point12B23B = [xg3a, yg3a];
     d9 = d3aIP;
  end
  if (d4aIP <d2aIP && d4aIP <d3aIP && d4aIP <d1aIP)
    Point12B23B = [xg4a, yg4a];
     d9 = d4aIP;
  end
D= [D ; Point12B23B];
 catch
     try
xg1a=round(vpa(IntersectionPoint3(1,1)));
xg2a=round(vpa(IntersectionPoint3(3,1)));
yg1a=round(vpa(IntersectionPoint3(2,1)));
yg2a=round(vpa(IntersectionPoint3(4,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);


  if (d1aIP <d2aIP )
    Point12B23B = [xg1a, yg1a];
    d9 = d1aIP;
  end
  if (d2aIP <d1aIP)
    Point12B23B = [xg2a, yg2a];
     d9 = d2aIP;
  end
  
D= [D ; Point12B23B];
catch
         MBx1=round(vpa(IntersectionPoint2(1,1)));
MBy1=round(vpa(IntersectionPoint2(2,1)));
          
         
xg =MBx1+pomA; 
eq12bG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12B^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12B^2)==1;
GH1=round(vpa(solve(eq12bG,yg)));
d1G=sqrt((MBy1-GH1(1,1)).^2);
d2G=sqrt((MBy1-GH1(2,1)).^2);
if d1G<d2G
Point12B23B= [xg,GH1(1,1)];
d9=sqrt((ipx-xg).^2+(ipy-GH1(1,1)).^2);

end
if d1G>d2G
Point12B23B = [xg,GH1(2,1)];   
d9=sqrt((ipx-xg).^2+(ipy-GH1(2,1)).^2);
end
         
        D= [D ;Point12B23B]; 
     end
 end
 syms xg;
end

%10
disp('10)')
[Xga,Yga] = solve(eq12b,eq23a);
xg1a=round(vpa(Xga(1,1)));
xg2a=round(vpa(Xga(2,1)));
xg3a=round(vpa(Xga(3,1)));
xg4a=round(vpa(Xga(4,1)));
yg1a=round(vpa(Yga(1,1)));
yg2a=round(vpa(Yga(2,1)));
yg3a=round(vpa(Yga(3,1)));
yg4a=round(vpa(Yga(4,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);
d3aIP=sqrt((ipx-xg3a).^2+(ipy-yg3a).^2);
d4aIP=sqrt((ipx-xg4a).^2+(ipy-yg4a).^2);

try
    if (xg1a >3000||xg1a <-3000||yg1a <-3000||yg1a <-3000 || xg2a >3000||xg2a <-3000||yg2a <-3000||yg2a <-3000)
        block 
    end
    
  if (d1aIP <d2aIP && d1aIP <d3aIP && d1aIP <d4aIP)
    Point12B23B = [xg1a, yg1a];
    d10 = d1aIP;
  end
  if (d2aIP <d1aIP && d2aIP <d3aIP && d2aIP <d4aIP)
    Point12B23B = [xg2a, yg2a];
     d10 = d2aIP;
  end
  if (d3aIP <d2aIP && d3aIP <d1aIP && d3aIP <d4aIP)
    Point12B23B = [xg3a, yg3a];
     d10 = d3aIP;
  end
  if (d4aIP <d2aIP && d4aIP <d3aIP && d4aIP <d1aIP)
    Point12B23B = [xg4a, yg4a];
     d10 = d4aIP;
  end
D= [D ; Point12B23B];
catch
  xg=Zacatek;
  IntersectionPoint2=[];
  pom2=0;
  pom3=1;
  IntersectionPoint3=[];
  pomA = 0;
  pomB = 0;
  pomC = 0;
  pomD = 0;  
  for i=1:Konec

eq12bG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12B^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12B^2)==1;
GH1=round(vpa(solve(eq12bG,yg)));

Point1=[xg;GH1(1,1)];
Point2=[xg;GH1(2,1)];
eq23aG = (((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23A^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23A^2)==1;
GH2=round(vpa(solve(eq23aG,yg)));
Point3=[xg;GH2(1,1)];
Point4=[xg;GH2(2,1)];

d13=sqrt((xg-xg).^2+(GH1(1,1)-GH2(1,1)).^2);
d14=sqrt((xg-xg).^2+(GH1(1,1)-GH2(2,1)).^2);
d23=sqrt((xg-xg).^2+(GH1(2,1)-GH2(1,1)).^2);
d24=sqrt((xg-xg).^2+(GH1(2,1)-GH2(2,1)).^2);
try

if d13<Polomer
    Podezrely = [xg;GH1(1,1)];
end
if d14<Polomer
    Podezrely = [xg;GH1(1,1)];
end
if d23<Polomer
    Podezrely = [xg;GH1(2,1)];
end
if d24<Polomer
    Podezrely = [xg;GH1(2,1)];
end
if Podezrely(1,1)==xg %kdyz je nejaky Podezrely

try %zkus
    
[dm,dn] = size(IntersectionPoint2);%zjisti kolik je uz pruseciku v mnozine
dl=dm/2; %pomocny vypocet, pruseiku je v mnozine tolik kolik je radku deleno dvema, protoze kazdy podezrely ma na 1 radku x-ovou a na druhem y-ovou souradnici

if dl==0% kdyz v mnozine zadny jeste neni tak priradine prvniho podezrelyho jako prusecik
    IntersectionPoint2 = [Podezrely];
    Podezrely=null; 
else %kdyz tam zadny neni
    
for j=1:dl % cyklus, ktery se opakuje tolikrat kolik je pruseciku v mnozine
    
    blabla = IntersectionPoint2(pom3,1); %vyber spravneho pruseciku
    
   
    if    (blabla+2+pom2)==xg 
      switch pom3
        case 1
        pomA=pomA+1; 
        case 3
        pomB=pomB+1; 
        case 5
        pomC=pomC+1; 
        case 7
        pomD=pomD+1; 
        end
        pom2=pom2+2;
        Podezrely=null;
  % kdyz je tento podezrely(xg je jeho x-va souradnice) priliz blizko pruseciku..
  %Pozn: jedna se o to, ze pro jeden prusecik, nalezne algoritmus vice vhodnych bodu, nam staci 1, ostatni nevybereme    
        
        
        
    else
      
       IntersectionPoint2 = [IntersectionPoint2;Podezrely];
       pom3=pom3+2;
       pom2=0;
       Podezrely=null;
    end
    end

end
end
end
end
xg=xg+2;
  end

  try
MBx1=IntersectionPoint2(1,1);
MBy1=IntersectionPoint2(2,1);
MBx2=IntersectionPoint2(3,1);
MBy2=IntersectionPoint2(4,1);

%MB1
xg =MBx1+pomA; 
eq12bG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12B^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12B^2)==1;
GH1=round(vpa(solve(eq12bG,yg)));

d1G=sqrt((MBy1-GH1(1,1)).^2);
d2G=sqrt((MBy1-GH1(2,1)).^2);
if d1G<d2G
IntersectionPoint3 = [xg;GH1(1,1)];

end
if d1G>d2G
IntersectionPoint3 = [xg;GH1(2,1)];    
end
%MB2
xg =MBx2+pomB; 
eq12bG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12B^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12B^2)==1;
GH1=round(vpa(solve(eq12bG,yg)));

d1G=sqrt((MBy2-GH1(1,1)).^2);
d2G=sqrt((MBy2-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
try
MBx3=IntersectionPoint2(5,1);
MBy3=IntersectionPoint2(6,1);
MBx4=IntersectionPoint2(7,1);
MBy4=IntersectionPoint2(8,1);
%MB3
xg =MBx3+pomC; 
eq12bG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12B^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12B^2)==1;
GH1=round(vpa(solve(eq12bG,yg)));

d1G=sqrt((MBy3-GH1(1,1)).^2);
d2G=sqrt((MBy3-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
%MB4
xg =MBx4+pomD; 
eq12bG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12B^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12B^2)==1;
GH1=round(vpa(solve(eq12bG,yg)));

d1G=sqrt((MBy4-GH1(1,1)).^2);
d2G=sqrt((MBy4-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
end
end
 try 
xg1a=round(vpa(IntersectionPoint3(1,1)));
xg2a=round(vpa(IntersectionPoint3(3,1)));
xg3a=round(vpa(IntersectionPoint3(5,1)));
xg4a=round(vpa(IntersectionPoint3(7,1)));
yg1a=round(vpa(IntersectionPoint3(2,1)));
yg2a=round(vpa(IntersectionPoint3(4,1)));
yg3a=round(vpa(IntersectionPoint3(6,1)));
yg4a=round(vpa(IntersectionPoint3(8,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);
d3aIP=sqrt((ipx-xg3a).^2+(ipy-yg3a).^2);
d4aIP=sqrt((ipx-xg4a).^2+(ipy-yg4a).^2);

  if (d1aIP <d2aIP && d1aIP <d3aIP && d1aIP <d4aIP)
    Point12B23B = [xg1a, yg1a];
    d10 = d1aIP;
  end
  if (d2aIP <d1aIP && d2aIP <d3aIP && d2aIP <d4aIP)
    Point12B23B = [xg2a, yg2a];
     d10 = d2aIP;
  end
  if (d3aIP <d2aIP && d3aIP <d1aIP && d3aIP <d4aIP)
    Point12B23B = [xg3a, yg3a];
     d10 = d3aIP;
  end
  if (d4aIP <d2aIP && d4aIP <d3aIP && d4aIP <d1aIP)
    Point12B23B = [xg4a, yg4a];
     d10 = d4aIP;
  end
D= [D ; Point12B23B];
 catch
     try
xg1a=round(vpa(IntersectionPoint3(1,1)));
xg2a=round(vpa(IntersectionPoint3(3,1)));
yg1a=round(vpa(IntersectionPoint3(2,1)));
yg2a=round(vpa(IntersectionPoint3(4,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);


  if (d1aIP <d2aIP )
    Point12B23B = [xg1a, yg1a];
    d10 = d1aIP;
  end
  if (d2aIP <d1aIP)
    Point12B23B = [xg2a, yg2a];
     d10 = d2aIP;
  end
  
D= [D ; Point12B23B];
catch
         MBx1=round(vpa(IntersectionPoint2(1,1)));
MBy1=round(vpa(IntersectionPoint2(2,1)));
          
         
xg =MBx1+pomA; 
eq12bG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12B^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12B^2)==1;
GH1=round(vpa(solve(eq12bG,yg)));
d1G=sqrt((MBy1-GH1(1,1)).^2);
d2G=sqrt((MBy1-GH1(2,1)).^2);
if d1G<d2G
Point12B23B= [xg,GH1(1,1)];
d10=sqrt((ipx-xg).^2+(ipy-GH1(1,1)).^2);

end
if d1G>d2G
Point12B23B = [xg,GH1(2,1)];   
d10=sqrt((ipx-xg).^2+(ipy-GH1(2,1)).^2);
end
         
        D= [D ;Point12B23B]; 
     end
 end
 syms xg;
end


%11
disp('11)')
[Xga,Yga] = solve(eq13a,eq23a);
xg1a=round(vpa(Xga(1,1)));
xg2a=round(vpa(Xga(2,1)));
xg3a=round(vpa(Xga(3,1)));
xg4a=round(vpa(Xga(4,1)));
yg1a=round(vpa(Yga(1,1)));
yg2a=round(vpa(Yga(2,1)));
yg3a=round(vpa(Yga(3,1)));
yg4a=round(vpa(Yga(4,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);
d3aIP=sqrt((ipx-xg3a).^2+(ipy-yg3a).^2);
d4aIP=sqrt((ipx-xg4a).^2+(ipy-yg4a).^2);


try
    if (xg1a >3000||xg1a <-3000||yg1a <-3000||yg1a <-3000 || xg2a >3000||xg2a <-3000||yg2a <-3000||yg2a <-3000)
        block 
    end
    
  if (d1aIP <d2aIP && d1aIP <d3aIP && d1aIP <d4aIP)
    Point12B23B = [xg1a, yg1a];
    d11 = d1aIP;
  end
  if (d2aIP <d1aIP && d2aIP <d3aIP && d2aIP <d4aIP)
    Point12B23B = [xg2a, yg2a];
     d11 = d2aIP;
  end
  if (d3aIP <d2aIP && d3aIP <d1aIP && d3aIP <d4aIP)
    Point12B23B = [xg3a, yg3a];
     d11 = d3aIP;
  end
  if (d4aIP <d2aIP && d4aIP <d3aIP && d4aIP <d1aIP)
    Point12B23B = [xg4a, yg4a];
     d11 = d4aIP;
  end
D= [D ; Point12B23B];
catch
  xg=Zacatek;
  IntersectionPoint2=[];
  pom2=0;
  pom3=1;
  IntersectionPoint3=[];
  pomA = 0;
  pomB = 0;
  pomC = 0;
  pomD = 0;  
  for i=1:Konec

eq13aG = (((xg-m13)*cosq13+(yg-n13)*sinq13).^2)/a13A^2 -(((-(xg-m13)*sinq13+(yg-n13)*cosq13).^2)/b13A^2)==1;
GH1=round(vpa(solve(eq13aG,yg)));

Point1=[xg;GH1(1,1)];
Point2=[xg;GH1(2,1)];
eq23aG = (((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23A^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23A^2)==1;
GH2=round(vpa(solve(eq23aG,yg)));
Point3=[xg;GH2(1,1)];
Point4=[xg;GH2(2,1)];

d13=sqrt((xg-xg).^2+(GH1(1,1)-GH2(1,1)).^2);
d14=sqrt((xg-xg).^2+(GH1(1,1)-GH2(2,1)).^2);
d23=sqrt((xg-xg).^2+(GH1(2,1)-GH2(1,1)).^2);
d24=sqrt((xg-xg).^2+(GH1(2,1)-GH2(2,1)).^2);
try

if d13<Polomer
    Podezrely = [xg;GH1(1,1)];
end
if d14<Polomer
    Podezrely = [xg;GH1(1,1)];
end
if d23<Polomer
    Podezrely = [xg;GH1(2,1)];
end
if d24<Polomer
    Podezrely = [xg;GH1(2,1)];
end
if Podezrely(1,1)==xg %kdyz je nejaky Podezrely

try %zkus
    
[dm,dn] = size(IntersectionPoint2);%zjisti kolik je uz pruseciku v mnozine
dl=dm/2; %pomocny vypocet, pruseiku je v mnozine tolik kolik je radku deleno dvema, protoze kazdy podezrely ma na 1 radku x-ovou a na druhem y-ovou souradnici

if dl==0% kdyz v mnozine zadny jeste neni tak priradine prvniho podezrelyho jako prusecik
    IntersectionPoint2 = [Podezrely];
    Podezrely=null; 
else %kdyz tam zadny neni
    
for j=1:dl % cyklus, ktery se opakuje tolikrat kolik je pruseciku v mnozine
    
    blabla = IntersectionPoint2(pom3,1); %vyber spravneho pruseciku
    
   
    if    (blabla+2+pom2)==xg 
      switch pom3
        case 1
        pomA=pomA+1; 
        case 3
        pomB=pomB+1; 
        case 5
        pomC=pomC+1; 
        case 7
        pomD=pomD+1; 
        end
        pom2=pom2+2;
        Podezrely=null;
  % kdyz je tento podezrely(xg je jeho x-va souradnice) priliz blizko pruseciku..
  %Pozn: jedna se o to, ze pro jeden prusecik, nalezne algoritmus vice vhodnych bodu, nam staci 1, ostatni nevybereme    
        
        
        
    else
      
       IntersectionPoint2 = [IntersectionPoint2;Podezrely];
       pom3=pom3+2;
       pom2=0;
       Podezrely=null;
    end
    end

end
end
end
end
xg=xg+2;
  end

  try
MBx1=IntersectionPoint2(1,1);
MBy1=IntersectionPoint2(2,1);
MBx2=IntersectionPoint2(3,1);
MBy2=IntersectionPoint2(4,1);

%MB1
xg =MBx1+pomA; 
eq23aG = (((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23A^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23A^2)==1;
GH1=round(vpa(solve(eq23aG,yg)));
d1G=sqrt((MBy1-GH1(1,1)).^2);
d2G=sqrt((MBy1-GH1(2,1)).^2);
if d1G<d2G
IntersectionPoint3 = [xg;GH1(1,1)];

end
if d1G>d2G
IntersectionPoint3 = [xg;GH1(2,1)];    
end
%MB2
xg =MBx2+pomB; 
eq23aG = (((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23A^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23A^2)==1;
GH1=round(vpa(solve(eq23aG,yg)));
d1G=sqrt((MBy2-GH1(1,1)).^2);
d2G=sqrt((MBy2-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
try
MBx3=IntersectionPoint2(5,1);
MBy3=IntersectionPoint2(6,1);
MBx4=IntersectionPoint2(7,1);
MBy4=IntersectionPoint2(8,1);
%MB3
xg =MBx3+pomC; 
eq23aG = (((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23A^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23A^2)==1;
GH1=round(vpa(solve(eq23aG,yg)));
d1G=sqrt((MBy3-GH1(1,1)).^2);
d2G=sqrt((MBy3-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
%MB4
xg =MBx4+pomD; 
eq23aG = (((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23A^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23A^2)==1;
GH1=round(vpa(solve(eq23aG,yg)));
d1G=sqrt((MBy4-GH1(1,1)).^2);
d2G=sqrt((MBy4-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
end
end
 try 
xg1a=round(vpa(IntersectionPoint3(1,1)));
xg2a=round(vpa(IntersectionPoint3(3,1)));
xg3a=round(vpa(IntersectionPoint3(5,1)));
xg4a=round(vpa(IntersectionPoint3(7,1)));
yg1a=round(vpa(IntersectionPoint3(2,1)));
yg2a=round(vpa(IntersectionPoint3(4,1)));
yg3a=round(vpa(IntersectionPoint3(6,1)));
yg4a=round(vpa(IntersectionPoint3(8,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);
d3aIP=sqrt((ipx-xg3a).^2+(ipy-yg3a).^2);
d4aIP=sqrt((ipx-xg4a).^2+(ipy-yg4a).^2);

  if (d1aIP <d2aIP && d1aIP <d3aIP && d1aIP <d4aIP)
    Point12B23B = [xg1a, yg1a];
    d11 = d1aIP;
  end
  if (d2aIP <d1aIP && d2aIP <d3aIP && d2aIP <d4aIP)
    Point12B23B = [xg2a, yg2a];
     d11 = d2aIP;
  end
  if (d3aIP <d2aIP && d3aIP <d1aIP && d3aIP <d4aIP)
    Point12B23B = [xg3a, yg3a];
     d11 = d3aIP;
  end
  if (d4aIP <d2aIP && d4aIP <d3aIP && d4aIP <d1aIP)
    Point12B23B = [xg4a, yg4a];
     d11 = d4aIP;
  end
D= [D ; Point12B23B];
 catch
     try
xg1a=round(vpa(IntersectionPoint3(1,1)));
xg2a=round(vpa(IntersectionPoint3(3,1)));
yg1a=round(vpa(IntersectionPoint3(2,1)));
yg2a=round(vpa(IntersectionPoint3(4,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);


  if (d1aIP <d2aIP )
    Point12B23B = [xg1a, yg1a];
    d11 = d1aIP;
  end
  if (d2aIP <d1aIP)
    Point12B23B = [xg2a, yg2a];
     d11 = d2aIP;
  end
  
D= [D ; Point12B23B];
catch
              MBx1=round(vpa(IntersectionPoint2(1,1)));
MBy1=round(vpa(IntersectionPoint2(2,1)));
          
         
xg =MBx1+pomA; 
eq23aG = (((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23A^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23A^2)==1;
GH1=round(vpa(solve(eq23aG,yg)));
d1G=sqrt((MBy1-GH1(1,1)).^2);
d2G=sqrt((MBy1-GH1(2,1)).^2);
if d1G<d2G
Point12B23B= [xg,GH1(1,1)];
d11=sqrt((ipx-xg).^2+(ipy-GH1(1,1)).^2);

end
if d1G>d2G
Point12B23B = [xg,GH1(2,1)];   
d11=sqrt((ipx-xg).^2+(ipy-GH1(2,1)).^2);
end
         
        D= [D ;Point12B23B]; 
     end
 end
 syms xg;
end

%12
disp('12)')
[Xga,Yga] = solve(eq12a,eq23a);
xg1a=round(vpa(Xga(1,1)));
xg2a=round(vpa(Xga(2,1)));
xg3a=round(vpa(Xga(3,1)));
xg4a=round(vpa(Xga(4,1)));
yg1a=round(vpa(Yga(1,1)));
yg2a=round(vpa(Yga(2,1)));
yg3a=round(vpa(Yga(3,1)));
yg4a=round(vpa(Yga(4,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);
d3aIP=sqrt((ipx-xg3a).^2+(ipy-yg3a).^2);
d4aIP=sqrt((ipx-xg4a).^2+(ipy-yg4a).^2);

try
    if (xg1a >3000||xg1a <-3000||yg1a <-3000||yg1a <-3000 || xg2a >3000||xg2a <-3000||yg2a <-3000||yg2a <-3000)
        block 
    end
    
  if (d1aIP <d2aIP && d1aIP <d3aIP && d1aIP <d4aIP)
    Point12B23B = [xg1a, yg1a];
    d12 = d1aIP;
  end
  if (d2aIP <d1aIP && d2aIP <d3aIP && d2aIP <d4aIP)
    Point12B23B = [xg2a, yg2a];
     d12 = d2aIP;
  end
  if (d3aIP <d2aIP && d3aIP <d1aIP && d3aIP <d4aIP)
    Point12B23B = [xg3a, yg3a];
     d12 = d3aIP;
  end
  if (d4aIP <d2aIP && d4aIP <d3aIP && d4aIP <d1aIP)
    Point12B23B = [xg4a, yg4a];
     d12 = d4aIP;
  end
D= [D ; Point12B23B]
catch
  xg=Zacatek;
  IntersectionPoint2=[];
  pom2=0;
  pom3=1;
  IntersectionPoint3=[];
  pomA = 0;
  pomB = 0;
  pomC = 0;
  pomD = 0;
  for i=1:Konec

eq12aG = (((xg-m12)*cosq12+(yg-n12)*sinq12).^2)/a12A^2 -(((-(xg-m12)*sinq12+(yg-n12)*cosq12).^2)/b12A^2)==1;
GH1=round(vpa(solve(eq12aG,yg)));

Point1=[xg;GH1(1,1)];
Point2=[xg;GH1(2,1)];
eq23aG = (((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23A^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23A^2)==1;
GH2=round(vpa(solve(eq23aG,yg)));
Point3=[xg;GH2(1,1)];
Point4=[xg;GH2(2,1)];

d13=sqrt((xg-xg).^2+(GH1(1,1)-GH2(1,1)).^2);
d14=sqrt((xg-xg).^2+(GH1(1,1)-GH2(2,1)).^2);
d23=sqrt((xg-xg).^2+(GH1(2,1)-GH2(1,1)).^2);
d24=sqrt((xg-xg).^2+(GH1(2,1)-GH2(2,1)).^2);
try

if d13<Polomer
    Podezrely = [xg;GH1(1,1)];
end
if d14<Polomer
    Podezrely = [xg;GH1(1,1)];
end
if d23<Polomer
    Podezrely = [xg;GH1(2,1)];
end
if d24<Polomer
    Podezrely = [xg;GH1(2,1)];
end
if Podezrely(1,1)==xg %kdyz je nejaky Podezrely

try %zkus
    
[dm,dn] = size(IntersectionPoint2);%zjisti kolik je uz pruseciku v mnozine
dl=dm/2; %pomocny vypocet, pruseiku je v mnozine tolik kolik je radku deleno dvema, protoze kazdy podezrely ma na 1 radku x-ovou a na druhem y-ovou souradnici

if dl==0% kdyz v mnozine zadny jeste neni tak priradine prvniho podezrelyho jako prusecik
    IntersectionPoint2 = [Podezrely];
    Podezrely=null; 
else %kdyz tam zadny neni
    
for j=1:dl % cyklus, ktery se opakuje tolikrat kolik je pruseciku v mnozine
    
    blabla = IntersectionPoint2(pom3,1); %vyber spravneho pruseciku
    
   
    if    (blabla+2+pom2)==xg 
      switch pom3
        case 1
        pomA=pomA+1; 
        case 3
        pomB=pomB+1; 
        case 5
        pomC=pomC+1; 
        case 7
        pomD=pomD+1; 
        end
        pom2=pom2+2;
        Podezrely=null;
  % kdyz je tento podezrely(xg je jeho x-va souradnice) priliz blizko pruseciku..
  %Pozn: jedna se o to, ze pro jeden prusecik, nalezne algoritmus vice vhodnych bodu, nam staci 1, ostatni nevybereme    
        
        
        
    else
      
       IntersectionPoint2 = [IntersectionPoint2;Podezrely];
       pom3=pom3+2;
       pom2=0;
       Podezrely=null;
    end
    end

end
end
end
end
xg=xg+2;
  end

 
try
MBx1=IntersectionPoint2(1,1);
MBy1=IntersectionPoint2(2,1);
MBx2=IntersectionPoint2(3,1);
MBy2=IntersectionPoint2(4,1);

%MB1
xg =MBx1+pomA; 
eq23aG = (((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23A^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23A^2)==1;
GH1=round(vpa(solve(eq23aG,yg)));
d1G=sqrt((MBy1-GH1(1,1)).^2);
d2G=sqrt((MBy1-GH1(2,1)).^2);
if d1G<d2G
IntersectionPoint3 = [xg;GH1(1,1)];

end
if d1G>d2G
IntersectionPoint3 = [xg;GH1(2,1)];    
end
%MB2
xg =MBx2+pomB; 
eq23aG = (((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23A^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23A^2)==1;
GH1=round(vpa(solve(eq23aG,yg)));
d1G=sqrt((MBy2-GH1(1,1)).^2);
d2G=sqrt((MBy2-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
try
MBx3=IntersectionPoint2(5,1);
MBy3=IntersectionPoint2(6,1);
MBx4=IntersectionPoint2(7,1);
MBy4=IntersectionPoint2(8,1);
%MB3
xg =MBx3+pomC; 
eq23aG = (((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23A^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23A^2)==1;
GH1=round(vpa(solve(eq23aG,yg)));
d1G=sqrt((MBy3-GH1(1,1)).^2);
d2G=sqrt((MBy3-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
%MB4
xg =MBx4+pomD; 
eq23aG = (((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23A^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23A^2)==1;
GH1=round(vpa(solve(eq23aG,yg)));
d1G=sqrt((MBy4-GH1(1,1)).^2);
d2G=sqrt((MBy4-GH1(2,1)).^2);
if d1G<d2G
Dalsibod =[xg;GH1(1,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
if d1G>d2G
Dalsibod =[xg;GH1(2,1)];  
IntersectionPoint3 = [IntersectionPoint3 ;Dalsibod];

end
end
end
 try 
xg1a=round(vpa(IntersectionPoint3(1,1)));
xg2a=round(vpa(IntersectionPoint3(3,1)));
xg3a=round(vpa(IntersectionPoint3(5,1)));
xg4a=round(vpa(IntersectionPoint3(7,1)));
yg1a=round(vpa(IntersectionPoint3(2,1)));
yg2a=round(vpa(IntersectionPoint3(4,1)));
yg3a=round(vpa(IntersectionPoint3(6,1)));
yg4a=round(vpa(IntersectionPoint3(8,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);
d3aIP=sqrt((ipx-xg3a).^2+(ipy-yg3a).^2);
d4aIP=sqrt((ipx-xg4a).^2+(ipy-yg4a).^2);

  if (d1aIP <d2aIP && d1aIP <d3aIP && d1aIP <d4aIP)
    Point12B23B = [xg1a, yg1a];
    d12 = d1aIP;
  end
  if (d2aIP <d1aIP && d2aIP <d3aIP && d2aIP <d4aIP)
    Point12B23B = [xg2a, yg2a];
     d12 = d2aIP;
  end
  if (d3aIP <d2aIP && d3aIP <d1aIP && d3aIP <d4aIP)
    Point12B23B = [xg3a, yg3a];
     d12 = d3aIP;
  end
  if (d4aIP <d2aIP && d4aIP <d3aIP && d4aIP <d1aIP)
    Point12B23B = [xg4a, yg4a];
     d12 = d4aIP;
  end
D= [D ; Point12B23B]
 catch
     try
xg1a=round(vpa(IntersectionPoint3(1,1)));
xg2a=round(vpa(IntersectionPoint3(3,1)));
yg1a=round(vpa(IntersectionPoint3(2,1)));
yg2a=round(vpa(IntersectionPoint3(4,1)));

d1aIP=sqrt((ipx-xg1a).^2+(ipy-yg1a).^2);
d2aIP=sqrt((ipx-xg2a).^2+(ipy-yg2a).^2);


  if (d1aIP <d2aIP )
    Point12B23B = [xg1a, yg1a];
    d12 = d1aIP;
  end
  if (d2aIP <d1aIP)
    Point12B23B = [xg2a, yg2a];
     d12 = d2aIP;
  end
  
D= [D ; Point12B23B]
catch
               MBx1=round(vpa(IntersectionPoint2(1,1)));
MBy1=round(vpa(IntersectionPoint2(2,1)));
          
         
xg =MBx1+pomA; 
eq23aG = (((xg-m23)*cosq23+(yg-n23)*sinq23).^2)/a23A^2 -(((-(xg-m23)*sinq23+(yg-n23)*cosq23).^2)/b23A^2)==1;
GH1=round(vpa(solve(eq23aG,yg)));
d1G=sqrt((MBy1-GH1(1,1)).^2);
d2G=sqrt((MBy1-GH1(2,1)).^2);
if d1G<d2G
Point12B23B= [xg,GH1(1,1)];
d12=sqrt((ipx-xg).^2+(ipy-GH1(1,1)).^2);

end
if d1G>d2G
Point12B23B = [xg,GH1(2,1)];   
d12=sqrt((ipx-xg).^2+(ipy-GH1(2,1)).^2);
end
         
        D= [D ;Point12B23B] 
     end
 end
 syms xg;
end


%****************************************************************************
%HLEDANI PLOCHY
%teorie v BP 3.6 a 3.7
%mnozina vzdalenosti bodu od hlavniho pruseciku
E=[vpa(d1);vpa(d2);vpa(d3);vpa(d4);vpa(d5);vpa(d6);vpa(d7);vpa(d8);vpa(d9);vpa(d10);vpa(d11);vpa(d12)];
%spojeni mnoziny vzdalenosti s jejich body
F=[D , E];
%serazeni podle vzdalenosti, pote odstraneni poslednich 4 bodu a odstraneni
%vzdalenosti
F;
G = sortrows(F,3);
G;
H=G(1:8,:);
dH1=H(1,3);
%podminka, ze kdyz se jedna o priblizne pravidelnej sestiuhelnik, rovnou se
%vyhodnoti 6 prvnich body jako hledane
dH6=H(6,3);
dH=dH6-dH1;

%obcas se stane ze program spatne vyhodnoti
%plochu, kdyz je letoun blizko teziste, v tom pripade zapnout tuto podminku
%Vyhodnoti nejblizsich sest jako vyslednou plochu
%end teto metody je zvyraznen touto sipkou: <------------------ (radek
%4322)
%if dH<5
%    N=H(1:6,:);
%    N(:,3) = [];
%else
H(:,3) = [];





%vyrazene body 

I=G(9:12,:);
I(:,3) = [];

I;
%ted vyberu 1 bod z vyrazenych, a najdu k nemu nejblizsi bod, to bude prvni dvojice
%zbyle 2 body budou dalsi dvojice
i11=I(1,1);
i12=I(1,2);
i21=I(2,1);
i22=I(2,2);
i31=I(3,1);
i32=I(3,2);
i41=I(4,1);
i42=I(4,2);

di12=sqrt((i11-i21).^2+(i12-i22).^2);
di13=sqrt((i11-i31).^2+(i12-i32).^2);
di14=sqrt((i11-i41).^2+(i12-i42).^2);

if (di12 <di13 && di12 <di14 )
    Pair1=[i11, i12; i21, i22];
    Pair2=[i31, i32; i41 ,i42];
end
if (di13 <di12 && di13 <di14 )
    Pair1=[i11, i12; i31, i32];
    Pair2=[i21, i22; i41 ,i42];
end
if (di14 <di13 && di14 <di12 )
    Pair1=[i11, i12; i41, i42];
    Pair2=[i31, i32; i21 ,i22];
end
%nyni najdu bod z 8-mice, ktery je k me dvojici nejblize, ten bude urcite v konecne
%mnozine, udelam pro obe dvojice
h11=H(1,1);
h12=H(1,2);
h21=H(2,1);
h22=H(2,2);
h31=H(3,1);
h32=H(3,2);
h41=H(4,1);
h42=H(4,2);
h51=H(5,1);
h52=H(5,2);
h61=H(6,1);
h62=H(6,2);
h71=H(7,1);
h72=H(7,2);
h81=H(8,1);
h82=H(8,2);
dh1=sqrt((i11-h11).^2+(i12-h12).^2);
dh2=sqrt((i11-h21).^2+(i12-h22).^2);
dh3=sqrt((i11-h31).^2+(i12-h32).^2);
dh4=sqrt((i11-h41).^2+(i12-h42).^2);
dh5=sqrt((i11-h51).^2+(i12-h52).^2);
dh6=sqrt((i11-h61).^2+(i12-h62).^2);
dh7=sqrt((i11-h71).^2+(i12-h72).^2);
dh8=sqrt((i11-h81).^2+(i12-h82).^2);

J=[vpa(dh1);vpa(dh2);vpa(dh3);vpa(dh4);vpa(dh5);vpa(dh6);vpa(dh7);vpa(dh8)];
K=[H , J];
K = sortrows(K,3);
ClosestPoints=[K(1,1) , K(1,2)];

%pro druhou dvojici
i11=Pair2(1,1);
i12=Pair2(1,2);
dh1=sqrt((i11-h11).^2+(i12-h12).^2);
dh2=sqrt((i11-h21).^2+(i12-h22).^2);
dh3=sqrt((i11-h31).^2+(i12-h32).^2);
dh4=sqrt((i11-h41).^2+(i12-h42).^2);
dh5=sqrt((i11-h51).^2+(i12-h52).^2);
dh6=sqrt((i11-h61).^2+(i12-h62).^2);
dh7=sqrt((i11-h71).^2+(i12-h72).^2);
dh8=sqrt((i11-h81).^2+(i12-h82).^2);

J=[vpa(dh1);vpa(dh2);vpa(dh3);vpa(dh4);vpa(dh5);vpa(dh6);vpa(dh7);vpa(dh8)];
K=[H , J];
K = sortrows(K,3);
ClosestPoints=[ClosestPoints; K(1,1) , K(1,2)];
%k temto dvou bodum najdu nejblizsi 2, ktery jsou k nim nejblize a mam
%konecnou mnozinu bodu
i11=ClosestPoints(1,1);
i12=ClosestPoints(1,2);
dh1=sqrt((i11-h11).^2+(i12-h12).^2);
dh2=sqrt((i11-h21).^2+(i12-h22).^2);
dh3=sqrt((i11-h31).^2+(i12-h32).^2);
dh4=sqrt((i11-h41).^2+(i12-h42).^2);
dh5=sqrt((i11-h51).^2+(i12-h52).^2);
dh6=sqrt((i11-h61).^2+(i12-h62).^2);
dh7=sqrt((i11-h71).^2+(i12-h72).^2);
dh8=sqrt((i11-h81).^2+(i12-h82).^2);

J=[vpa(dh1);vpa(dh2);vpa(dh3);vpa(dh4);vpa(dh5);vpa(dh6);vpa(dh7);vpa(dh8)];
K=[H , J];
L = sortrows(K,3);

Lx1=double(L(1,1));
Ly1=double(L(1,2));
Lx2=double(L(2,1));
Ly2=double(L(2,2));
Lx3=double(L(3,1));
Ly3=double(L(3,2));
%pomocna metoda, nekdy se stava ze jiny bod je blize nez ktery hledame, a
%ten lezi na stejne hyperbole jako zbyle dva, tudiz uhel ktery sviraji je
%bud priblizne 0 nebo 180 stupnu a podle toho ho vyradime
% diky nepresnosti druhe metody hledani pruseciku( BP kapitola 3.5.2)
% musime nechat rezervu, zvolena 20 a 160 stupnu
qL=(atand((Ly1-Ly2)/(Lx1-Lx2)))-(atand((Ly3-Ly2)/(Lx3-Lx2)));
qcorectL=abs(qL);

if (qcorectL<20 || qcorectL>160)

L(3,1)=L(4,1);
L(3,2)=L(4,2);
L=L(1:3,:);
else
        L=L(1:3,:);
end

L(:,3) = [];
%druhy bod 
i11=ClosestPoints(2,1);
i12=ClosestPoints(2,2);
dh1=sqrt((i11-h11).^2+(i12-h12).^2);
dh2=sqrt((i11-h21).^2+(i12-h22).^2);
dh3=sqrt((i11-h31).^2+(i12-h32).^2);
dh4=sqrt((i11-h41).^2+(i12-h42).^2);
dh5=sqrt((i11-h51).^2+(i12-h52).^2);
dh6=sqrt((i11-h61).^2+(i12-h62).^2);
dh7=sqrt((i11-h71).^2+(i12-h72).^2);
dh8=sqrt((i11-h81).^2+(i12-h82).^2);

J=[vpa(dh1);vpa(dh2);vpa(dh3);vpa(dh4);vpa(dh5);vpa(dh6);vpa(dh7);vpa(dh8)];
K=[H , J];
M = sortrows(K,3);
M;
Mx1=double(M(1,1));
My1=double(M(1,2));
Mx2=double(M(2,1));
My2=double(M(2,2));
Mx3=double(M(3,1));
My3=double(M(3,2));

qM=(atand((My1-My2)/(Mx1-Mx2)))-(atand((My3-My2)/(Mx3-Mx2)));
qcorectM=abs(qM);

if (qcorectM<20 || qcorectM>160)

M(3,1)=M(4,1);
M(3,2)=M(4,2);      
end

if M(3,1) == L(3,1)
M(3,1)=M(5,1);
M(3,2)=M(5,2);  
end
M=M(1:3,:);
M(:,3) = [];
N=[L ; M]
%end<------------konec pomocne metody<------
%Nyni vyuzijeme funkci convhull, ktera nam vypocita plochu mnohouhelnika z
%mnoziny bodu N
o11=double(N(1,1));
o12=double(N(1,2));
o21=double(N(2,1));
o22=double(N(2,2));
o31=double(N(3,1));
o32=double(N(3,2));
o41=double(N(4,1));
o42=double(N(4,2));
o51=double(N(5,1));
o52=double(N(5,2));
o61=double(N(6,1));
o62=double(N(6,2));
xN=[o11,o21,o31,o41,o51,o61];
yN=[o12,o22,o32,o42,o52,o62];
[CH,	Area]	=	convhull (xN,yN);
Area
%****************************************************************************


A = [Stanice1; Stanice2 ; Stanice3];
 plot(A(:,1), A(:,2), '*')
 
 plot(C(:,1), C(:,2), 'O')

  %plot(N(:,1), N(:,2), 'O')
plot(xN(CH),yN(CH),'b',	xN,yN,'bo')
 plot(G(:,1), G(:,2), '*')
title(['Výpoèet polohy pomocí TDOA '])
grid on
axis equal
N=double(N);
D=double(D);
%zapis hodnot do souboru
fileID = fopen('hodnoty.txt','a');
fprintf(fileID,'%6s\r\n',' ');
fprintf(fileID,'%6s\r\n','*********************');
fprintf(fileID,'%6s\r\n','Zpoždìní:');
fprintf(fileID,'%6.9f\r\n',dt);
fprintf(fileID,'%6s\r\n','Poloha letadla:');
fprintf(fileID,'%6s %12s\r\n','x','y');
fprintf(fileID,'%6.0f %12.0f\r\n',P);
fprintf(fileID,'%6s\r\n','Množina bodù:');
fprintf(fileID,'%6s %12s\r\n','x','y');
fprintf(fileID,'%6.0f %12.0f\r\n',D);
fprintf(fileID,'%6s\r\n','Posledních šest:');
fprintf(fileID,'%6s %12s\r\n','x','y');
fprintf(fileID,'%6.0f %12.0f\r\n',N);
fprintf(fileID,'%6s','Plocha:');
fprintf(fileID,'%6.0f m^2\r\n',Area);
fclose(fileID);
toc