Mercurial > octave
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;