[Orekit Developers] recursive orbital estimation

[<brianh4321@hotmail.com>]

Hello/Bonjour

I've been testing a recursive orbital estimator using by importing the orekit
libraries directly into matlab. Being somewhat new to the field I've had many
questions and Luc and his team have been very gracious with their help. To
further these discussions in a larger domain I wish to post the following
questions.

Not sure if anyone works much with matlab, but I've been testing an recursive
orbital estimator using the new NumericalPropagatorWithJacobians (which imports
directly into matlab). 

1. The orbital state in numerical propagator is equinoctial expressed in
"true latitude", however most books seem to express equinoctial in terms of
longitude. I seem to find that while the equations match most reference as
longitude, the description of your functions seems to replace latitude for
longitude. Is there something I'm missing?

2. Also, most reference expresses equinoctial in terms of mean
latitude/longitude as oppose to true latitude/longitude. I was hoping to
find a reference that expresses the measurement Jacobian as a function of
your state. Do anyone know of such a reference? (A Jacobian relating your state
and a cartesian state would also work). I'm currently computing the
measurement Jacobian with finite differencing, but I'd love to have a closed
form.

3. Here is a sub-section of the Kalman estimator code. It works but the answers
don't seem to converge with the accuracy I'd expect. My initial state
covariance is somewhat of a guess, I don't have much experience with
equinoctial and they are somewhat of a guess. The problem is derived from an
example in Satellite Orbits, pg 291 (by Oliver MontenBruck)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Pest = diag([10e5,1e-2,1e-2,1e-4,1e-3,1e-3].^2);
Xest= Xref0;
Xref= Xref0;
b=zeros(3,N_obs);
sigma_angle= .01*pi/180;
sigma_range= 10;
Q= zeros(6);
linear= 1;
for ii=1:N_obs
    date= propagator.getInitialState().getDate();
    dt= Obs{ii,1}.durationFrom(date);
    
    % Propagate
    predictedState = propagator.propagate(date.shiftedBy(dt), PhiEq, dFdP);
    predictedPos= predictedState.getPVCoordinates().getPosition();
    % Predicted Measurements
    Elev= staF.getElevation(predictedPos,inertialFrame,Obs{ii,1});
    Azim= staF.getAzimuth(predictedPos,inertialFrame,Obs{ii,1});
    Range= staF.getRange(predictedPos,inertialFrame,Obs{ii,1});
    % Store predicted reference and compute Phi, Inverse Measurement Cov(R)
    % and Measurement Jacobian H (finite differencing)
    oXref= Xref;
    Xref= [ predictedState.getA();
	predictedState.getEquinoctialEx();
	predictedState.getEquinoctialEy();
	predictedState.getHx();
	predictedState.getHy();
	predictedState.getLv()];
    mPhiEq= ja2mat(PhiEq);
    mPhiEq= mPhiEq(1:6,1:6);
    Rinv= diag([sigma_angle*cos(Elev),sigma_angle,sigma_range].^2);
    G= computeH(staF,predictedState);
    
    % Time Update
    Pbar= mPhiEq*Pest*mPhiEq.'+Q;
    if(linear)
	% Linearized Kalman Update
	Xbar= Xref+ mPhiEq*(Xest-oXref);
	% Measurement Update
	b(:,ii)= [Obs{ii,2};Obs{ii,3};Obs{ii,4}]-[Azim;Elev;Range];
	K= Pbar*G'*inv(G*Pbar*G'+Rinv);
	Xest= Xbar+K*(b(:,ii)-G*(Xbar-Xref));
	Pest= Pbar-K*G*Pbar;
	propagator.setInitialState(predictedState);
    else
	% Extended Kalman Update
	Xbar= Xref;
	b(:,ii)= [Obs{ii,2};Obs{ii,3};Obs{ii,4}]-[Azim;Elev;Range];
	K= Pbar*G'*inv(G*Pbar*G'+Rinv);
	Pest= Pbar-K*G*Pbar;
	Xest= Xbar+K*b(:,ii);
	stateEst= EquinoctialOrbit(Xest(1),Xest(2), Xest(3), Xest(4), ...
	    Xest(5), Xest(6),EquinoctialOrbit.TRUE_LATITUDE_ARGUMENT,...
	    FramesFactory.getEME2000(), Obs{ii,1}, mu);
	propagator.setInitialState(SpacecraftState(stateEst));
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Any help or comment are appreciated. 

Brian