add radar measurement

Author: AndreyAPopov

toolbox/+otp/+cr3bp/+presets/CircularEarthOrbit.m

@@ -10,10 +10,12 @@
             % ----------
             % OrbitalRadius : numeric(1, 1)
             %    The radius of the orbit above Earth's surface in km.
+            %    The default radius is 340 km, which is approximately where
+            %    the ISS is/was located.
 
             p = inputParser();
             p.KeepUnmatched = true;
-            p.addParameter('OrbitalRadius', 340, @isnumeric);
+            p.addParameter('OrbitalRadius', 340, @(x) isnumeric(x) && isscalar(x));
             p.parse(varargin{:});
             results = p.Results;
             varargin = [fieldnames(p.Unmatched), struct2cell(p.Unmatched)].';
@@ -22,6 +24,7 @@
             earthMoonDist = 385000;
             orbitalradius = results.OrbitalRadius;
 
+            % set the relative distance from earth to the moon
             delta = (equatorialRadius + orbitalradius)/earthMoonDist;
 
             mE = otp.utils.PhysicalConstants.EarthMass;

toolbox/+otp/+cr3bp/CR3BPProblem.m

@@ -77,7 +77,7 @@
     properties (SetAccess = private)
         JacobiConstant
         JacobiConstantJacobian
-        RHSStab
+        RadarMeasurement
     end
 
     methods (Access = protected)
@@ -88,6 +88,8 @@ function onSettingsChanged(obj)
             obj.JacobiConstant = @(y) otp.cr3bp.jacobiconstant(y, mu, soft);
             obj.JacobiConstantJacobian = @(y) otp.cr3bp.jacobiconstantjacobian(y, mu, soft);
 
+            obj.RadarMeasurement = @(y, radary) otp.cr3bp.radarMeasurement(t, y, mu, soft, radary);
+
             obj.RHS = otp.RHS(@(t, y) otp.cr3bp.f(t, y, mu, soft), ...
                 'Jacobian', @(t, y) otp.cr3bp.jacobian(t, y, mu, soft), ...
                 'JacobianVectorProduct', @(t, y, v) otp.cr3bp.jacobianVectorProduct(t, y, v, mu, soft), ...

toolbox/+otp/+cr3bp/radarMeasurement.m

@@ -0,0 +1,16 @@
+function rmeas = radarMeasurement(~, y, ~, ~, radary)
+
+x = y(1:3, :);
+v = y(4:6, :);
+
+relpos = x - radary;
+
+r = sqrt(relpos(1, :).^2 + relpos(2, :).^2 + relpos(3, :).^2);
+drdt = sum(relpos.*v, 1)/r;
+theta = atan2(relpos(2, :), relpos(1, :));
+phi = atan2(relpos(3, :), r);
+
+
+rmeas = [r; drdt; theta; phi];
+
+end