function  [alpha,sigma,in_out, validez] = DME_DME_Upgrade_Area (lat_a, lon_a,  lat_D1, lon_D1, R1,  lat_D2, lon_D2, R2, PLOT, fh)
%% -----------------------------------------------------------------------------------------------------------
% Version:      22/09/2015    
% (c):          Vicente Calabuig Simón
%               vicasim@etsid.upv.es
%
% 1.- Utilidad:
%
%     Análisis de la figura III - 1-3-1 del documento OACI 8168 - V2
%               --> Área maxima de actualización de dos estaciones DME
%
% 2.- Parámetros de entrada:
% 
%       lat_a  : latitud de la aeronave            [deg]           lon_a  : longitud de la aeronave            [deg]   
%       
%       lat_D1 : latitud de la estación DME #1     [deg]           lon_D1 : longitud de la estación DME #1     [deg]
%          
%       lat_D2 : latitud de la estación DME #2     [deg]           lon_D2 : longitud de la estación DME #2     [deg] 
%
%       R1 = Rango de funcionamiento del DME #1    [NM]            R1 = Rango de funcionamiento del DME #2      [NM]
%
%       PLOT   : Puede valer 'P' ó 0.           P --> Con plot del escenario.
%                                               0 --> Sin plot del escenario.
%
%       fh   :  Figure handler                 0  --> La función crea su propia figura 
%                                              fh --> La función plotea sobre la figura especificada por fh
% 3.- Parámetros de salida:
%
%
%       alpha : Ángulo de intersección entre la aeronave y la pareja de DME's 
%
%       Sigma : Precisión de la pareja de DME's
%
%       in_out :  Puede valer '1' o '0'          1  --> Si la aeronave se encuentra DENTRO del área máxima de actualización
%                                                         2 --> Si la aeronave se encuentra FUERA del área máxima de actualización.
%
%       validez :          Puede valer '1' o '0'          1  --> Se verifican todos los requisitos (30º<alpha<150, sigma < 5MN, in_out = 1)
%                                                         0  --> Alguno de los requisitos no se cumple.
%                                                  
% 4.- Ejemplo de uso:    
%
%      [alpha,sigma,in_out, validez] = DME_DME_Upgrade_Area (10, 11,  12, 13, 120,  14, 13, 130, 'P', 0)       Cobertura Invalida
%
%      [alpha,sigma,in_out, validez] = DME_DME_Upgrade_Area(13, 12.5,  12, 13, 120,  14, 13, 130, 'P', 0)       Cobertura Valida
%
%      
%      [alpha,sigma,in_out, validez] = DME_DME_Upgrade_Area (39.489, -0.4816  ,39.484, -0.47031, 18,  39.707, -0.98643 , 130, 'P', 0)
%      
%      Estaciones DME muy cerca entre si
%      [alpha,sigma,in_out, validez] = DME_DME_Upgrade_Area (39.489, -0.4816,  39.484, -0.47031, 18,  39.486 , -0.48304 , 130, 'P', 0)
%      
%                                           
%
%% -----------------------------------------------------------------------------------------------------------------------------------
% Creacion propia figura // Dejando posibilidad a que sea ploteada en otra.
if strcmp(PLOT, 'P') 
    if  (fh == 0)
       update_area = figure('Name', 'Check Cobertura');  figure(update_area);
    else
        figure(fh);
    end;    
end

% Tiempo de espera y numero de puntos [DEPURACION]
t=0;    np = 200;        na = 200;

% Manejo de varias figuras:

% Conversion de unidades( NM -> Km )a los rangos:
NM2Km = 1.852;    Km2NM = 1/NM2Km;      R1 = R1 * NM2Km;       R2 = R2 * NM2Km;  


% Asumiendo la simplificacion de 'Tierra Plana' : 1 grado latitud  ~ 110 Km | 1 grado longitud ~ 111 Km  
latDeg2Km = 110;             Km2LatDeg = 1/latDeg2Km;       lonDeg2Km = 111;      Km2LonDeg = 1/lonDeg2Km; 
          
% Trabaremos ahora en unidades de distancia (Km)
x_a  = lat_a  * latDeg2Km ;     y_a  = lon_a  * lonDeg2Km ;            % Aeronave
x_D1 = lat_D1 * latDeg2Km ;     y_D1 = lon_D1 * lonDeg2Km ;            % DME 1
x_D2 = lat_D2 * latDeg2Km ;     y_D2 = lon_D2 * lonDeg2Km ;            % DME 2


%% PLOT RANGOS ESTACIONES DME
% Vectores genereales X e Y con las coordenadas de los puntos de interes
X0 = [x_a, x_D1, x_D2];        Y0 = [y_a, y_D1, y_D2];    

% Circunferencias con centro en cada estación DME y radio igual a su rango de funcionamiento:
a = linspace(0, 360, np);  NM2Km = 1.852;
% Circunferencia sobre el DME 1
x1 = R1*cosd(a) + x_D1;   y1 = R1*sind(a) + y_D1;  if strcmp('P', PLOT)  plot(x1, y1, '-black', 'LineWidth',2); hold on; end
% Circunferencia sobre el DME 2
x2 = R2*cosd(a) + x_D2;   y2 = R2*sind(a) + y_D2;  if strcmp('P', PLOT)  plot(x2, y2, '-black', 'LineWidth',2); hold on; end  

% Circunferencias de intercepción con radio 'D': D es la distancia entre los DME en Km
D = sqrt ( (x_D1 - x_D2 )^2 + (y_D1 - y_D2 )^2 );
%% PLOT AERONAVE Y DME's antes de regiones de intersección.
% Posiciones puntuales:
if strcmp('P', PLOT)
    for i = 1:3
        if      i == 1             % Posición Aeronave
            plot( X0(i), Y0(i),     '^red','LineWidth',4);        hold on;
        else                       % Posiciones DME 1 y DME 2
            plot( X0(i), Y0(i),     '^blue', 'LineWidth',4);      hold on;
        end
    end
end

% Segmento Aeronave - DME 1 - Prolongar un 10% de D en direccion Aeronave -> DME
ms1  = (y_D1 - y_a)/(x_D1 - x_a); 
if      x_a > x_D1   % Hay que incrementar la x hacia la izquierda (negativamente)
    s1 = -1;
elseif  x_a == x_D1  % Solo hay que incrementar la y
    s1 = 0;
elseif x_a < x_D1    % Hay que incrementar la x hacia la derecha, (positivamente)
    s1 = +1;
end
% Comprobacion pendiente infinita:
if ms1 ~= Inf && ms1 ~= -Inf
    b1 = y_a - ms1*x_a;    x_f1 = x_D1 + s1*(0.1)*D;        y_f1 = ms1*(x_f1) + b1;
elseif ms1 == -Inf
    b1 = 0;                x_f1 = x_D1;                     y_f1 = y_D1 - (0.1)*D;
elseif ms1 == Inf
    b1 = 0;                x_f1 = x_D1;                     y_f1 = y_D1 + (0.1)*D;
end
if strcmp('P', PLOT) plot([x_a, x_f1],[y_a, y_f1], '--black', 'LineWidth',2);      hold on; end
% Segmento Aeronave - DME 2
ms2  = (y_D2 - y_a)/(x_D2 - x_a);
if      x_a > x_D2   % Hay que incrementar la x hacia la izquierda (negativamente)
    s2 = -1;
elseif  x_a == x_D2  % Solo hay que incrementar la y
    s2= 0;
elseif x_a < x_D2    % Hay que incrementar la x hacia la derecha, (positivamente)
    s2 = +1;
end
% Comprobacion pendiente infinita:
if ms2 ~= Inf && ms2 ~= -Inf
    b2 = y_a - ms2*x_a;    x_f2 = x_D2 + s2*(0.1)*D;        y_f2 = ms2*(x_f2) + b2;
elseif ms2 == -Inf
    b2 = 0;                x_f2 = x_D2;                     y_f2 = y_D2 - (0.1)*D;
elseif ms2 == Inf
    b2 = 0;                x_f2 = x_D2;                     y_f2 = y_D2 + (0.1)*D;
end
if strcmp('P', PLOT) plot([x_a, x_f2],[y_a, y_f2], '--black', 'LineWidth',2);      hold on; end
if D>= R1+R2
    fprintf('Pareja de DMEs no válida, no hay zona de cobertura común\n');
    fprintf('Distancia entre las estaciones: D = %g Km\n', D);
    fprintf('Rango de la estación 1: R1 = %g Km\n', R1);       fprintf('Rango de la estación 2: R2 = %g Km\n', R2);
    fprintf('La suma de ambos rangos cubre una distancia R1+R2 = %g Km\n', R1+R2);
    fprintf('Quedando un espacio sin cobertura de dimensión: D-(R1+R2) = %g Km\n', D -(R1+R2));
    alpha = 0; sigma = 0; in_out = 0; validez = 0; if strcmp(PLOT, 'P');    figure(gcf);  hold off;    end
    return;
elseif (x_D1 == x_D2 && y_D1== y_D2) 
    fprintf('Error: Ambas estaciones DME no pueden tener las mismas coordenadas\n');
    return;
    
elseif  (x_a == x_D1 && y_a == y_D1) ||  (x_a == x_D2 && y_a == y_D2)
    fprintf('Error: Aeronave situada en las mismas coordenadas que una de las estaciones DME\n');
    if update_area;  close gcf; end
    return;
end

%% Circunferencias de intersección con radio D
D = sqrt ( (x_D1 - x_D2 )^2 + (y_D1 - y_D2 )^2 );
syms x0s y0s
ec1 = (x_D1 - x0s)^2 + (y_D1 - y0s)^2 - D^2;             ec2 = (x_D2 - x0s)^2 + (y_D2 - y0s)^2 - D^2;
[x0s, y0s] = solve(ec1, ec2);                            x0n = double(x0s);         y0n = double(y0s);
xi1 = D*cosd(a) + x0n(1);  yi1 = D*sind(a) + y0n(1);    if strcmp('P', PLOT) plot(xi1, yi1, '-blue', 'LineWidth',2); hold on;  end
xi2 = D*cosd(a) + x0n(2);  yi2 = D*sind(a) + y0n(2);    if strcmp('P', PLOT)plot(xi2, yi2, '-blue', 'LineWidth',2); hold on;   end
%% Circunferencias de actualizacion de radio mínimo:
Rmin = 1*NM2Km;
xm1 = Rmin*cosd(a) + x_D1;  ym1 = Rmin*sind(a) + y_D1;     if strcmp('P', PLOT)  plot(xm1, ym1, '-red', 'LineWidth',2); hold on; end
xm2 = Rmin*cosd(a) + x_D2;  ym2 = Rmin*sind(a) + y_D2;     if strcmp('P', PLOT)plot(xm2, ym2, '-red', 'LineWidth',2); hold on;   end

%% CALCULO del Ángulo de Interseccion alpha y de la precision sigma de la pareja:

% Vector director unitario del segmento que va desde la aeronave hasta el DME_1
vd1(1) = x_D1 - x_a;        vd1(2) = y_D1 - y_a;       u1 = vd1/norm(vd1); 
% Vector director unitario del segmento que va desde la aeronave hasta el DME_2
vd2(1) = x_D2 - x_a;        vd2(2) = y_D2 - y_a;       u2 = vd2/norm(vd2); 

% Alpha como el ángulo formado entre ambos vectores directores
CosenoAlpha = ( abs( u1(1)*u2(1) +  u1(2)*u2(2)) ) / ( sqrt(u1(1)^2 + u1(2)^2)* sqrt(u2(1)^2 + u2(2)^2) );
alpha = acosd(CosenoAlpha);   alpha = real(alpha);   % [º]   % La formula siempre devuelve alpha 

% Queremos el ángulo INTERIOR a los segmentos
% Necesitamos la distancia de la aeronave a cada estacion DME: en Km
d1 = sqrt ( (x_D1 - x_a)^2 + (y_D1 - y_a)^2 );  d2 = sqrt ( (x_D2 - x_a)^2 + (y_D2 - y_a)^2 );   %Km
if D> max(d1, d2)
    alpha = 180 - alpha;
end

% Calculo de la precisión sigma de la pareja de DMEs:  OACI 8168 V2 pág 537
s_sis = 0.05 * NM2Km; 
s1_air = max (0.085*NM2Km, 0.00125*d1);          s2_air = max (0.085*NM2Km, 0.00125*d2);             %Km

sigma = sqrt( (s1_air^2 + s_sis^2) + (s2_air^2 + s_sis^2) )/( sind (alpha) );  
sigma = sigma*Km2NM;

%% REGIONES DE INTERSECCION:
% Intersección entre las circunferencias de rango máximo de cada estación DME
from_2_into_1 = inpolygon(x2,y2,x1,y1);         from_1_into_2 = inpolygon(x1,y1,x2,y2);        
xri = [];   yri =[]; 
for i = 1 : length(x2)
    if      from_2_into_1(i)
                xri = [xri; x2(i)];     yri = [yri; y2(i)];
                %plot(x2(i), y2(i), 'oblack', 'LineWidth',2);                  pause(t);
    elseif  from_1_into_2(i)
                xri = [xri; x1(i)];     yri = [yri; y1(i)];
                %plot(x1(i), y1(i), 'oblack', 'LineWidth',2);                  pause(t);
    end
end  
[xri, yri] = Arrange(xri, yri);
% Region anterior que además intersecta con la Circunferencia de interseccion 1   : Región de actualizacion 1
xri1 = [];   yri1 =[];
check = inpolygon(xri,yri,xi1,yi1) ;              
for i = 1 : length(xri)
    if      check(i)
                xri1 = [xri1; xri(i)];     yri1 = [yri1; yri(i)];
                %plot(xri(i), yri(i), 'ogreen', 'LineWidth',2);                  pause(t);
    end
end
check = inpolygon(xi1,yi1,xri,yri) ;
for i = 1:length(xi1);
    if  check(i)
                xri1 = [xri1; xi1(i)];     yri1 = [yri1; yi1(i)]; 
                %plot(xi1(i), yi1(i), 'ogreen', 'LineWidth',2);                  pause(t);
    end
end
[xri1, yri1] = Arrange(xri1, yri1);
if strcmp('P', PLOT)
    fill(xri1, yri1, 'green') 
end
% Region anterior que además intersecta con la Circunferencia de interseccion 2    : Región de actualización 2
xri2 = [];   yri2 =[];
check = inpolygon(xri,yri,xi2,yi2) ;              
for i = 1 : length(xri)
    if      check(i)
                xri2 = [xri2; xri(i)];     yri2 = [yri2; yri(i)]; 
                %plot(xri(i), yri(i), 'ogreen', 'LineWidth',2);                  pause(t);
    end
end
check = inpolygon(xi2,yi2,xri,yri) ;
for i = 1:length(xi2);
    if  check(i)
                xri2 = [xri2; xi2(i)];     yri2 = [yri2; yi2(i)];
                %plot(xi2(i), yi2(i), 'ogreen', 'LineWidth',2);                     pause(t);
    end
end  
[xri2, yri2] = Arrange(xri2, yri2);
if strcmp('P', PLOT)
    fill(xri2, yri2, 'green')  
end

% Region de interseccion entre ambas circunferencias de interseccion: NO actualizable 
xri3 = [];   yri3 =[];
check = inpolygon(xi1,yi1,xi2,yi2) ;              
for i = 1 : length(xi1)
    if      check(i)
                xri3 = [xri3; xi1(i)];     yri3 = [yri3; yi1(i)];  
    end
end
check = inpolygon(xi2,yi2,xi1,yi1) ;
for i = 1:length(xi2);
    if  check(i)
                xri3 = [xri3; xi2(i)];     yri3 = [yri3; yi2(i)];
    end
end
[xri3, yri3] = Arrange(xri3, yri3);
if strcmp('P', PLOT)
    fill(xri3, yri3, 'white')
end
% Regiones definidas por las circunferencias de radio mínimo:
if strcmp('P', PLOT)
    fill(xm1, ym1, 'white');            fill(xm2, ym2, 'white');        plot( x_a, y_a, '^red','LineWidth',4);    
end
[xri3, yri3] = Arrange(xri3, yri3);
% Comprobamos si la aeronave está dentro o fuera de la región de actualización:
check1 = inpolygon(x_a, y_a, xri1, yri1);       
check2 = inpolygon(x_a, y_a, xri2, yri2);        check3 = inpolygon(x_a, y_a, xri3, yri3);
check4 = inpolygon(x_a, y_a, xm1,  ym1);         check5 = inpolygon(x_a, y_a, xm2,  ym2);

if check3 || check4 || check5
    in_out = 0;
elseif  check1 || check2
    in_out = 1; 
else 
    in_out = 0;
end

%% Validez final de la cobertura, una vez se conocen alpha, sigma e in_out:
if ((30<alpha && alpha<150)  && (in_out) && (sigma<5))   validez = 1;      else        validez = 0;       end 
%% PLOT AERONAVE Y DME's
% Posiciones puntuales:
if strcmp('P', PLOT)
    for i = 1:3
        if      i == 1             % Posición Aeronave
            plot( X0(i), Y0(i),     '^red','LineWidth',4);        hold on;
        else                       % Posiciones DME 1 y DME 2
            plot( X0(i), Y0(i),     '^blue', 'LineWidth',4);      hold on;
        end
    end
end

% Segmento Aeronave - DME 1 - Prolongar un 40% de D en direccion Aeronave -> DME
ms1  = (y_D1 - y_a)/(x_D1 - x_a); 
if      x_a > x_D1   % Hay que incrementar la x hacia la izquierda (negativamente)
    s1 = -1;
elseif  x_a == x_D1  % Solo hay que incrementar la y
    s1 = 0;
elseif x_a < x_D1    % Hay que incrementar la x hacia la derecha, (positivamente)
    s1 = +1;
end
% Comprobacion pendiente infinita:
if ms1 ~= Inf && ms1 ~= -Inf
    b1 = y_a - ms1*x_a;    x_f1 = x_D1 + s1*(0.4)*D;        y_f1 = ms1*(x_f1) + b1;
elseif ms1 == -Inf
    b1 = 0;                x_f1 = x_D1;                     y_f1 = y_D1 - (0.4)*D;
elseif ms1 == Inf
    b1 = 0;                x_f1 = x_D1;                     y_f1 = y_D1 + (0.4)*D;
end
if strcmp('P', PLOT) plot([x_a, x_f1],[y_a, y_f1], '--black', 'LineWidth',2);      hold on; end
% Segmento Aeronave - DME 2
ms2  = (y_D2 - y_a)/(x_D2 - x_a);
if      x_a > x_D2   % Hay que incrementar la x hacia la izquierda (negativamente)
    s2 = -1;
elseif  x_a == x_D2  % Solo hay que incrementar la y
    s2= 0;
elseif x_a < x_D2    % Hay que incrementar la x hacia la derecha, (positivamente)
    s2 = +1;
end
% Comprobacion pendiente infinita:
if ms2 ~= Inf && ms2 ~= -Inf
    b2 = y_a - ms2*x_a;    x_f2 = x_D2 + s2*(0.4)*D;        y_f2 = ms2*(x_f2) + b2;
elseif ms2 == -Inf
    b2 = 0;                x_f2 = x_D2;                     y_f2 = y_D2 - (0.4)*D;
elseif ms2 == Inf
    b2 = 0;                x_f2 = x_D2;                     y_f2 = y_D2 + (0.4)*D;
end
if strcmp('P', PLOT) plot([x_a, x_f2],[y_a, y_f2], '--black', 'LineWidth',2);      hold on; end

%% PLOT DEL ANGULO
if strcmp('P', PLOT)
       Beta1 = atand(ms1);         Beta2 = atand(ms2);
    % Ángulo inicial a utilizar:   9 casos posible
    X = [min(x_D1,x_D2), max(x_D1, x_D2)];          Y = [min(y_D1,y_D2), max(y_D1, y_D2)];
    
   % Ecuación recta que pasa por ambas estaciones DME
   mDD  = (Y(2) - Y(1) )/(X(2) - X(1));   bDD = Y(1) - mDD*X(1); 
   
   % Centroide del triangulo formado por ambas estaciones y la aeronave:
    x_ct = (x_a + x_D1 + x_D2) /3;       y_ct = (y_a + y_D1 + y_D2) /3;
    D_ct_a = sqrt(   (x_ct - x_a)^2   +   (y_ct - y_a)^2   ); 
   % Distancias a las estaciones DME
    d = [d1, d2];                  
    if      X(1) < x_a && x_a < X(2)  % RECTA CENTRAL
        RA = D/2;
        if X(1) == x_a && x_a == X(2)
                    B  = atand(mDD);                        a0 = B;            sa = 1;
        end
        if      (y_a == mDD*x_a + bDD) % sobre la recta:   alpha = 180º                       
                    B  = atand(mDD);                        a0 = B;            sa = 1;                       
        elseif y_a > mDD*x_a + bDD  % Por encima
                    B  = max(atand(ms1), atand(ms2));       a0 = B +180;       sa = 1; 

        elseif y_a < mDD*x_a + bDD    %Por debajo
                    B  = max(atand(ms1), atand(ms2));       a0 = B;            sa = 1; 
        end
    elseif      X(2)<= x_a    % A la DERECHA de la recta                        EN ELLO ESTAMOS
        if      min(d) == d(1)   RA = d1/2;  elseif  min(d) == d(2)   RA = d2/2;  end
        if      (y_a == mDD*x_a + bDD) % Sobre la recta:   alpha = 0º       
                    B  = atand(mDD);                        a0 = B;            sa = 1;                       
        elseif y_a > mDD*x_a + bDD  % Por encima
                    B  = min(atand(ms1), atand(ms2));    
               if x_a ~= X(2)         a0 = B +180;       sa = 1;      else    a0 = B;            sa = -1;        end
                    
        elseif y_a < mDD*x_a + bDD    %Por debajo
               B  = min(atand(ms1), atand(ms2));
               if x_a ~= X(2)         a0 = B + 180;      sa = 1;  else  a0 = B + 180;            sa = -1;        end
        end
    elseif      x_a <= X(1)   % A la IZQUIERDA
        if      min(d) == d(1)   RA = d1/2;  elseif  min(d) == d(2)   RA = d2/2;  end
        if      (y_a == mDD*x_a + bDD) % Sobre la recta:   alpha = 0º
                    B  = atand(mDD);                        a0 = B;            sa = 1;                       
        elseif y_a > mDD*x_a + bDD  % Por encima
                    B  = min(atand(ms1), atand(ms2));       
                    if x_a~= X(1)     a0 = B;            sa = 1; else   a0 = B;            sa = 1;        end

        elseif y_a < mDD*x_a + bDD    %Por debajo
                    B  = min(atand(ms1), atand(ms2));        
                    if x_a~= X(1)       a0 = B;             sa = 1;     else      a0 = B;             sa = 1;   end
        end
    end
    af = a0 + (sa)*alpha;
    % Incremento angular en coordenadas polares
     a0r = a0 * pi/180;       afr = af * pi/180;              theta = linspace(a0r, afr, na);       rho  = ones(1,na) * RA;
    
    % Transformación a coordenadas cartesianas
    [x,y] = pol2cart(theta,rho);        x = x + x_a;      y = y + y_a;        plot(x,y, '-red', 'LineWidth', 2);        hold on
    
    % Ángulo numerico como texto
    L  = length(x);     xt = x(L/2);               yt = y(L/2);
    s1 = '\alpha = ';   s2 = num2str(alpha);    string = [s1, s2, 'º'];
    text(xt,yt, string);     axis equal; 
   
    %% ENFOQUE de los elementos : Tiene preferencia siempre el trianguo formado por los DME y la aeronave
    X0 = [x_a, x_D1, x_D2];        Y0 = [y_a, y_D1, y_D2];          dprefs = [d1, d2, D];    dpref = max(dprefs);
    X = X0;                        Y = Y0;
    if R1 < 3.5*dpref         % Orden de magnitud de R1 es tal que pueda graficarse sin afectar a la visualizacion del triangulo? SI, lo añadimos
        X = [X, x1];    Y = [Y, y1];
    end
    if R2 < 3.5*dpref         
        X = [X, x2];    Y = [Y, y2];
    end
    if D < 3.5*dpref         
        X = [X, xi1, xi2];    Y = [Y, yi1, yi2];
    end
        
    x_min = min(X);     x_max = max(X);          y_min = min(Y);     y_max = max(Y);        axis([x_min x_max y_min y_max])
    %axis([(x_a- 3*dpref ) (x_a+ 3*dpref) (y_a- 3*dpref) (y_a+ 3*dpref)]);
   
                       
    grid off;    axis equal;     xlabel('Eje Y: Latitud en Km (1º~ 110 Km)','FontSize',20);    ylabel('Eje X: Longitud en Km (1º~ 111 Km)','FontSize',20);  hold off;

    % TEXTOS en gráfica: Cuadro resumen en la esquina superior derecha             
    % alpha entre 30º y 150º:
    sa1 = '30º < \alpha < 150º : ';    if (30<alpha && alpha<150)  sa2 = ['\color{black}SI ', '\color{green} \surd'];  else       sa2 = ['\color{black}NO ', '\color{red} X'] ;   end
    sa = [sa1, sa2];                                       
    % in_out 
    si_o1 = '\color{black}Dentro-Fuera :  ';       if (in_out)  si_o2 = ['\color{black}DENTRO ', '\color{green} \surd'];  else si_o2 = ['\color{black}FUERA ','\color{red} X'];       end    
    si_o = [si_o1, si_o2];                        
    % Precision sigma 
    sig1 = '\color{black}Precisión \sigma :  ';     sig2 = num2str(sigma); sig3 = ' NM';   if sigma<5 sig4 =['\color{green} \surd']; else sig4 =['\color{red} X']; end
    sig = [sig1, sig2, sig3, sig4];               
    % Evaluacion final de la cobertura:
    fin1 = '\color{black}\bfCobertura Válida: ';
    if validez   fin2= '\color{black}SI';      else        fin2 = '\color{black}NO';       end  
    fin = [fin1, fin2]; 
    % Escritura de las lineas en el cuadro de texto
    annotation(gcf,'textbox','String',{sa,'',si_o,'', sig,'' ,fin}, ...
   'FontSize',14,'units','pix', 'Position', [1500 700 400 200],'Interpreter','Tex', 'BackgroundColor', 'white');

    % LONGITUD AL EJE X  // LATITUD AL EJE Y
    view(-90,90) 
    set(gca,'ydir','reverse')
2+2;
end
end

function [x_ord, y_ord] = Arrange(x, y)
for i = 1 : (length(x)-1)        % Pasamos por todos los puntos menos el ultimo
        d = [];
    for j = (i+1): length(x)     % Anteriores ya ordenados
        d(j) = sqrt( (x(i) - x(j) )^2 + ( y(i) - y(j) )^2 );
    end
    % El punto siguiente más cercano tiene como índice imin:
    ipos = find(d>0);
    if isempty(ipos)
       break;
       fprintf('ipos está vacio\n');
       return;
    end
    [value,imin] = min(d(ipos)); 
    imin = imin + i;
    % Colócalo justo siguiendo al punto actual: intercambia posciones
    save  = x(i+1);   
    x(i+1)= x(imin);  
    x(imin) = save;
   
    save  = y(i+1);   
    y(i+1)= y(imin);   
    y(imin) = save;
end
% Por si el punto final fuera problemático
x (end) = x(1);                  y(end) = y(1);

% Finalmente:
x_ord = x;      y_ord = y;
end