function coe = coe_from_sv(R,V,mu)
% –––––––––––––––––––––––––––––––––––––––––––––---------------------------
%% Notas aclaratorias 
% Esta función calcula los elementos orbitales clásicos (COE)
% a partir del vector de estado (R,V)
%
% mu   - parámetro gravitacional (km^3/s^2)
% R    - vector de posición en el marco ecuatorial geocéntrico (km)
% V    - vector de velocidad en el marco ecuatorial geocéntrico (km)
% r, v - las magnitudes de R y V
% vr   - componente de velocidad radial (km/s)
% H    - el vector de momento angular (km^2/s)
% h    - la magnitud de H (km^2/s)
% incl - inclinación de la órbita (rad)
% N    - el vector de la línea de nodos (km^2/s)
% n    - la magnitud de N
% cp   - producto cruzado de N y R
% RA   - ascensión recta del nodo ascendente (rad)
% E    - vector de excentricidad
% e    - excentricidad (magnitud de E)
% eps  - un número pequeño por debajo del cual se considera que
%        la excentricidad es cero
% w    - argumento del perigeo (rad)
% TA   - anomalía verdadera (rad)
% a    - semieje mayor (km)
% pi   - 3.1415926...
% coe  - vector de elementos orbitales [h e RA incl w TA a]
%
% Funciones M del usuario requeridas: Ninguna
% –––––––––––––––––––––––––––––––––––––––––––––---------------------------


eps = 1.e-10;

% Calculo los modulos del vector posicion y velocidad
r = norm(R);
v = norm(V);

% Calculo de velocidad radial
vr = dot(R,V)/r;

% CALCULO DE LOS PARAMETROS ORBITALES
%% Momento angular especifico y su modulo
H = cross(R,V);
h = norm(H);

%% Inclincacion
incl = acos(H(3)/h);
% Recuerde que i debe estar entre 0° y 180°, que es precisamente el rango
% (valores principales) de la función arcocoseno. Por lo tanto, no hay
% ambigüedad de cuadrante con la que lidiar aquí. Si 90° < i ≤ 180°, el
% momento angular h apunta en dirección sur. En ese caso, la órbita es
% retrógrada, lo que significa que el movimiento del satélite alrededor de
% la Tierra es opuesto a la rotación de la Tierra.     

%% Ascensión recta del nodo ascendente (RAAN o Ω),
N = cross([0 0 1],H);
n = norm(N);
if n ~= 0
    if N(2) >= 0
        RA = acos(N(1)/n);
    else 
        RA = 2*pi - acos(N(1)/n);
    end
else
    RA = 0;
end
% La ascensión recta del nodo ascendente, Ω = cos⁻¹(Nx/N). 
% Si Nx > 0, entonces Ω se encuentra en el primer o cuarto cuadrante. 
% Si Nx < 0, entonces Ω se encuentra en el segundo o tercer cuadrante.  
% Por lo tanto, Ny > 0 implica que 0 ≤ Ω < 180°, 
% mientras que Ny < 0 implica que 180° ≤ Ω < 360°. 

%% Excentricidad
E = 1/mu*((v^2 - mu/r)*R - r*vr*V);
e = norm(E);

%% Argumento del perigeo
if n ~= 0
    if e > eps
        if E(3) >= 0
            w = acos(dot(N,E)/n/e);
        else
            w = 2*pi - acos(dot(N,E)/n/e);
        end
    else
        w = 0;
    end
else
    w = 0;
end
% Si N·e > 0, entonces ω se encuentra en el primer o cuarto cuadrante. 
% Si N·e < 0, entonces ω se encuentra en el segundo o tercer cuadrante. 
% Para ubicar ω en el cuadrante adecuado, observa que el perigeo 
% se encuentra por encima del plano ecuatorial (0° ≤ ω < 180°) si e apunta hacia arriba 
% (en la dirección positiva del eje Z) y que el perigeo se encuentra por debajo del plano 
% (180° ≤ ω < 360°) si e apunta hacia abajo. 
% Por lo tanto, eZ ≥ 0 implica que 0° ≤ ω < 180°, mientras que eZ < 0 implica que 
% 180° ≤ ω < 360°. 

%% Anomalia verdadera
if e > eps
    if vr >= 0
        TA = acos(dot(E,R)/e/r);
    else
        TA = 2*pi - acos(dot(E,R)/e/r);
    end
else
    cp = cross(N,R);
    if cp(3) >= 0
        TA = acos(dot(N,R)/n/r);
    else
        TA = 2*pi - acos(dot(N,R)/n/r);
    end
end


%% Semieje mayor
a = h^2/mu/(1 - e^2);

%% Parametros orbitales de salida
coe = [h e RA incl w TA a];
end %coe_from_sv