changeset 29193:b2d5ee958d7f stable

ode_event_handler.m: Fix mishandling of event edge types and multiple events (bug #59709). * ode_event_handler.m: Correctly find events where either the direction does not matter (dir == 0) or where the edge type matches that specified in direction (dir == sign (evt)). When multiple events occur, just choose the first valid one and return the index and time of the event based on that. This simplifies the calculation of tnew and is Matlab-compatible.
author Rik <rik@octave.org>
date Thu, 17 Dec 2020 15:57:51 -0800
parents ae5d758c10e1
children 01e5c00a8609 900246a8bb37
files scripts/ode/private/ode_event_handler.m
diffstat 1 files changed, 8 insertions(+), 15 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/ode/private/ode_event_handler.m	Thu Dec 17 14:23:19 2020 -0800
+++ b/scripts/ode/private/ode_event_handler.m	Thu Dec 17 15:57:51 2020 -0800
@@ -101,18 +101,10 @@
     ## Check if one or more signs of the event has changed
     signum = (sign (evtold) != sign (evt));
     if (any (signum))         # One or more values have changed
-      idx = find (signum);    # Get the index of the changed values
-
-      if (any (dir(idx) == 0))
-        ## Rising or falling (both are possible)
-        ## Don't change anything, keep the index
-      elseif (any (dir(idx) == sign (evt(idx))))
-        ## Detected rising or falling, need a new index
-        idx = find (dir == sign (evt));
-      else
-        ## Found a zero crossing but must not be notified
-        idx = [];
-      endif
+      ## Find events where either rising and falling edges are counted (dir==0)
+      ## or where the specified edge type matches the event edge type.
+      idx = signum & (dir == 0 | dir == sign (evt));
+      idx = find (idx);  # convert logical to numeric index or []
 
       ## Create new output values if a valid index has been found
       if (! isempty (idx))
@@ -123,11 +115,12 @@
         else
           retcell{1} = any (term(idx));     # Stop integration or not
         endif
-        retcell{2}(evtcnt,1) = idx(1,1);  # Take first event found
+        idx = idx(1);  # Use first event found if there are multiple.
+        retcell{2}(evtcnt,1) = idx;
         ## Calculate the time stamp when the event function returned 0 and
         ## calculate new values for the integration results, we do both by
-        ## a linear interpolation
-        tnew = t - evt(1,idx) * (t - told) / (evt(1,idx) - evtold(1,idx));
+        ## a linear interpolation.
+        tnew = t - evt(idx) * (t - told) / (evt(idx) - evtold(idx));
         ynew = (y - (t - tnew) * (y - yold) / (t - told)).';
         retcell{3}(evtcnt,1) = tnew;
         retcell{4}(evtcnt,:) = ynew;