Error in ode45 and 'Events' of a Satellite Tracking Radar Model

I'm trying to simulate a satellite tracker. With ode45 I have de state of the satellite for each time and then convert position to Azimuth-Elevation local coordinates. If the object enters in the field of view (expressed as a mask of Az-El), the integration stops and changes the step to a smaller one in order to obtain more points inside. I'm trying to use an 'Events' function, which compares de Az-El of the satellite with the radar's mask. The problem is ode45 detects an event when it's not supposed to, and vice versa. I discarded the transformation to angular coordinates as the origin of the problem.

load("maskmin.mat")
load("maskmax.mat")
Azmaskmin = maskmin(:,1);
Elmaskmin = maskmin(:,2);
Azmaskmax = maskmax(:,1);
Elmaskmax = maskmax(:,2);
t_total = [];    y_total = [];
az_total = [];   el_total = [];
t_actual = t0;   y_actual = y0;
options = odeset('RelTol',3e-14,'AbsTol',3e-14,'Refine',10,'NormControl','on',...
    'Events', @(t,y) eventos_radar(t, y, lat_radar,...
    lon_radar, alt_radar,Azmaskmin,Azmaskmax,Elmaskmax,Elmaskmin));
aze = [];
ele = [];
yee = [];
tee = [];
% Principal loop
r1 = 100; % Big step
rt = 10; % Small step
dt = r1; % Initial condition isn't inside FOV
while t_actual < tf - dt
    [t, y, te, ye, ie] = ode45(@(t,y)ecuacion_movimiento(t, y, mu),...
        [t_actual:dt:tf], y_actual, options);
    t_total = [t_total; t];
    y_total = [y_total; y];
    % Update next integration
    t_actual = t(end);
    y_actual = y(end,:)';
    % Event check
    if ~isempty(ie)
        if dt == r1
            dt = rt;
            fprintf('Entrance in t = %.2f segundos\n', t_actual);
        elseif dt == rt
            dt = r1;
            fprintf('Exit in t = %.2f segundos\n', t_actual);
        end
        t_actual = te(end);
        y_actual = ye(end, :)';
    end
end
% Procesing results
procesar_resultados(t_total, y_total, lat_radar, lon_radar, alt_radar,...
    Azmaskmin,Azmaskmax,Elmaskmax,Elmaskmin);

function [value, isterminal, direction] = eventos_radar(t, y, lat_radar,...
    lon_radar, alt_radar,Azmaskmin,Azmaskmax,Elmaskmax,Elmaskmin)
[az, el] = calcular_az_el(y(1:3)', lat_radar, lon_radar, alt_radar, t);
if az>=Azmaskmax(1) && az<=Azmaskmax(end)
    el_min = interp1(Azmaskmin, Elmaskmin, az);
    el_max = interp1(Azmaskmax, Elmaskmax, az);
    value = [el_max - el, el - el_min];
    isterminal = [1 1]; % Detect changes of sign
else
    value = [Elmaskmax(1) - el, el - Elmaskmax(1)];
    % Elmaskmax(1) is the El value from the FOV corner in order to get continuous values of 'value' when the object is             % outside
    isterminal = [0 0]; % Ignores changes of sign
end
direction = [0 0];

The left plot is using the y_total values once the simulation stops, and the rigth plot is using the values of the y that the event function receives from the ode 45.

Anyone knows what is going on? Any help is welcome!!

Answer 1

[Torsten]

The "tspan" vector you provide for "ode45" has no influence on the internal time steps taken by the integrator. In zones with steep gradients of the solution variables, "ode45" automatically chooses smaller time steps that will be able to resolve the solution curve. Thus I suggest using no event function, defining "tspan" as a two-element vector [ts,tf] and starting with tolerances (RelTol, AbsTol) in the order of 1e-8. If necessary, you can make them even more stringent, but values of 1e-14 are too low in my opinion.

You didn't show us the function "ecuacion_movimiento" yet. If you have defined time points where the derivative function becomes discontinuous, you should integrate within such a time interval and restart the solver at the time when a discontinuity happens. If the time points where the derivative function becomes discontinuous are not known in advance and are implicitly defined, you should choose these implicit conditions as events.