// Copyright (c) CNES 2008
//
// This software is part of CelestLab, a CNES toolbox for Scilab
//
// This software is governed by the CeCILL license under French law and
// abiding by the rules of distribution of free software. You can use,
// modify and/ or redistribute the software under the terms of the CeCILL
// license as circulated by CEA, CNRS and INRIA at the following URL
// 'http://www.cecill.info'.
function [kep,jacob] = CL__oe_equin2kep(equin)
// Equinoctial to Keplerian orbital elements
//
// Calling Sequence
// [kep, jacob] = CL__oe_equin2kep(equin)
//
// Description
//
//
Converts equinoctial orbital elements to classical Keplerian orbital elements.
//
The transformation jacobian is optionally computed.
//
See Orbital elements for more details
//
//
//
// Parameters
// equin: Equinoctial orbital elements [sma;ex;ey;hx;hy;L] [m,rad] (6xN)
// kep: Classical Keplerian orbital elements [sma;ecc;inc;pom;gom;anm] [m,rad] (6xN)
// jacob: (optional) Transformation jacobian (See Orbital elements for more details) (6x6xN)
//
// Authors
// CNES - DCT/SB
//
// Examples
// // Example 1
// equin = [7000.e3; 0.1; 0.2; 0.3; 0.4; 1];
// kep = CL__oe_equin2kep(equin);
//
// // Example 2
// equin = [7000.e3; 0.1; 0.2; 0.3; 0.4; 1];
// [kep, jacob1] = CL__oe_equin2kep(equin);
// [equin2, jacob2] = CL__oe_kep2equin(kep);
// equin2 - equin // => zero
// jacob2 * jacob1 // => identity
// Declarations:
EPS_ORB = CL__dataGetEnv("epsOrb", internal=%t);
// Code:
// Handle [] cases
if (equin == [])
kep = [];
jacob = [];
return;
end
// Check validity of input
[isvalid,type_orbit] = CL__oe_isValid("equin",equin);
if (~isvalid); CL__error("Invalid orbital elements"); end;
if (find(type_orbit <> 1) <> []); CL__error("Invalid orbital elements (parabolic or hyperbolic orbit)"); end;
kep = zeros(equin);
// Notes on particular cases handled:
// 1) If the orbit is nearly circular, pom is not defined
// (Circular treshold : ecc < EPS_ORB.cir)
// In that case pom is set to 0 (arbitrary choice)
//
// 1) If the orbit is nearly equatorial, gom is not defined
// (Equatorial treshold : sin(inc) < EPS_ORB.equa)
// In that case gom is set to 0 (arbitrary choice)
// Conversion formulas:
// a(kep) = a(equin)
// e = sqrt(ex^2+ey^2)
// i = 2*atan(sqrt(hx^2+hy^2))
// gom = atan(hy,hx)
// pom = atan(ey,ex) - gom
// M = L - pom - gom
kep(1,:) = equin(1,:);
kep(2,:) = sqrt(equin(2,:).^2 + equin(3,:).^2);
kep(3,:) = 2 * atan(sqrt(equin(4,:).^2 + equin(5,:).^2));
kep(5,:) = atan(equin(5,:),equin(4,:));
Ieq = find(sin(kep(3,:)) < EPS_ORB.equa); // equatorial orbit
kep(5,Ieq) = 0;
kep(4,:) = atan(equin(3,:),equin(2,:)) - kep(5,:);
Icir = find(kep(2,:) < EPS_ORB.cir); // circular orbit
kep(4,Icir) = 0;
kep(6,:) = equin(6,:) - kep(4,:) - kep(5,:);
// Reduce angles to 0 -> 2pi : necessary?
// kep(3:6,:) = CL_rMod(kep(3:6,:),2*%pi)
// Jacobian computation (dkep/dequin)
jacob = [];
if (argn(1) == 2)
// jacob(i,j) = d(kep_i)/d(equin_j)
//
// Formulas used:
// da/da = 1
// de/dex = cos(pom+gom)
// de/dey = sin(pom+gom)
// di/dhx = 2hx/(sqrt(hx^2+hy^2)(1+hx^2+hy^2))
// di/dhy = 2hy/(sqrt(hx^2+hy^2)(1+hx^2+hy^2))
// dpom/dex = -ey/e2
// dpom/dey = ex/e2
// dpom/dhx = hy / (hx^2+hy^2)
// dpom/dhy = -hx / (hx^2+hy^2)
// dgom/dhx = -hy / (hx^2+hy^2)
// dgom/dhy = hx / (hx^2+hy^2)
// dM/dex = ey / e2
// dM/dey = -ex / e2
// dM/dL = 1
// Circular orbit (not computable results --> set to %nan)
ecir = kep(2,:);
ecir(Icir) = %nan;
// Equatorial orbit (not computable results --> set to %nan)
hx2_hy2 = equin(4,:).^2 + equin(5,:).^2;
hx2_hy2(Ieq) = %nan;
N = size(equin,2);
jacob = zeros(6,6,N);
jacob(1,1,:) = 1; // da/da
jacob(2,2,:) = equin(2,:) ./ ecir; // de/dex
jacob(2,3,:) = equin(3,:) ./ ecir; // de/dey
jacob(3,4,:) = 2*equin(4,:) ./ (sqrt(hx2_hy2) .* (1+hx2_hy2)); // di/dhx
jacob(3,5,:) = 2*equin(5,:) ./ (sqrt(hx2_hy2) .* (1+hx2_hy2)); // di/dhy
jacob(4,2,:) = -equin(3,:) ./ ecir.^2; // dpom/dex
jacob(4,3,:) = equin(2,:) ./ ecir.^2; // dpom/dey
jacob(4,4,:) = equin(5,:) ./ hx2_hy2; // dpom/dhx
jacob(4,5,:) = -equin(4,:) ./ hx2_hy2; // dpom/dhy
jacob(5,4,:) = -equin(5,:) ./ hx2_hy2; // dgom/dhx
jacob(5,5,:) = equin(4,:) ./ hx2_hy2; // dgom/dhy
jacob(6,2,:) = equin(3,:) ./ ecir.^2; // dM/dex
jacob(6,3,:) = -equin(2,:) ./ ecir.^2; // dM/dey
jacob(6,6,:) = 1; // dM/dL
end
endfunction