1########################################################################
2##
3## Copyright (C) 2006-2021 The Octave Project Developers
4##
5## See the file COPYRIGHT.md in the top-level directory of this
6## distribution or <https://octave.org/copyright/>.
7##
8## This file is part of Octave.
9##
10## Octave is free software: you can redistribute it and/or modify it
11## under the terms of the GNU General Public License as published by
12## the Free Software Foundation, either version 3 of the License, or
13## (at your option) any later version.
14##
15## Octave is distributed in the hope that it will be useful, but
16## WITHOUT ANY WARRANTY; without even the implied warranty of
17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18## GNU General Public License for more details.
19##
20## You should have received a copy of the GNU General Public License
21## along with Octave; see the file COPYING.  If not, see
22## <https://www.gnu.org/licenses/>.
23##
24########################################################################
25
26## -*- texinfo -*-
27## @deftypefn {} {@var{retval} =} ode_event_handler (@var{@@evtfun}, @var{t}, @var{y}, @var{flag}, @var{par1}, @var{par2}, @dots{})
28##
29## Return the solution of the event function (@var{@@evtfun}) which is
30## specified in the form of a function handle.
31##
32## The second input argument @var{t} is a scalar double and specifies the time
33## of the event evaluation.
34##
35## The third input argument @var{y} may be a column vector of type double
36## (for ODEs and DAEs) which specifies the solutions.  Alternatives, @var{y}
37## may be a cell array (for IDEs and DDEs) which specifies the solutions and
38## derivatives.
39##
40## The fourth input argument @var{flag} is of type string.  Valid values are:
41##
42## @table @option
43## @item  @qcode{"init"}
44## Initialize internal persistent variables of the function
45## @code{ode_event_handler} and return an empty cell array of size 4.
46##
47## @item  @qcode{"calc"}
48## Evaluate the event function and return the solution @var{retval} as a cell
49## array of size 4.
50##
51## @item  @qcode{"done"}
52## Clean up internal variables of the function @code{ode_event_handler} and
53## return an empty cell array of size 4.
54## @end table
55##
56## If additional input arguments @var{par1}, @var{par2}, @dots{} are given
57## these parameters are passed directly to the event function.
58##
59## This function is an ODE internal helper function and it should never be
60## necessary to call it directly.
61## @end deftypefn
62
63function retval = ode_event_handler (evtfun, t, y, flag = "", varargin)
64
65  ## No error handling has been implemented in this function to achieve
66  ## the highest performance possible.
67
68  ## retval{1} is true (to terminate) or false (to continue)
69  ## retval{2} is the index information for which an event occurred
70  ## retval{3} is the time information column vector
71  ## retval{4} is the line by line result information matrix
72
73  ## These persistent variables store the results and time value from the
74  ## processing in the previous time stamp.
75  ## evtold  the results from the event function
76  ## told    the time stamp
77  ## yold    the ODE result
78  ## retcell the return values cell array
79  ## evtcnt  the counter for how often this function has been called
80  persistent evtold told yold retcell;
81  persistent evtcnt = 1;   # Don't remove.  Required for Octave parser.
82  persistent firstrun = true;
83
84  if (isempty (flag))
85    ## Process the event, i.e.,
86    ## find the zero crossings for either a rising or falling edge
87    if (! iscell (y))
88      inpargs = {evtfun, t, y};
89    else
90      inpargs = {evtfun, t, y{1}, y{2}};
91      y = y{1};  # Delete cell element 2
92    endif
93    if (nargin > 4)
94      inpargs = {inpargs{:}, varargin{:}};
95    endif
96    [evt, term, dir] = feval (inpargs{:});
97
98    ## We require that all return values be row vectors
99    evt = evt(:).'; term = term(:).'; dir = dir(:).';
100
101    ## Check if one or more signs of the event has changed
102    signum = (sign (evtold) != sign (evt));
103    if (any (signum))         # One or more values have changed
104      ## Find events where either rising and falling edges are counted (dir==0)
105      ## or where the specified edge type matches the event edge type.
106      idx = signum & (dir == 0 | dir == sign (evt));
107      idx = find (idx);  # convert logical to numeric index or []
108
109      ## Create new output values if a valid index has been found
110      if (! isempty (idx))
111        ## Change the persistent result cell array
112        if (firstrun)
113          ## Matlab compatibility requires ignoring condition on first run.
114          retcell{1} = false;
115        else
116          retcell{1} = any (term(idx));     # Stop integration or not
117        endif
118        idx = idx(1);  # Use first event found if there are multiple.
119        retcell{2}(evtcnt,1) = idx;
120        ## Calculate the time stamp when the event function returned 0 and
121        ## calculate new values for the integration results, we do both by
122        ## a linear interpolation.
123        tnew = t - evt(idx) * (t - told) / (evt(idx) - evtold(idx));
124        ynew = (y - (t - tnew) * (y - yold) / (t - told)).';
125        retcell{3}(evtcnt,1) = tnew;
126        retcell{4}(evtcnt,:) = ynew;
127        evtcnt += 1;
128      endif
129
130    endif
131
132    firstrun = false;
133    evtold = evt; told = t; yold = y;
134    retval = retcell;
135
136  elseif (strcmp (flag, "init"))
137    ## Call the event function if an event function has been defined to
138    ## initialize the internal variables of the event function and to get
139    ## a value for evtold.
140
141    firstrun = true;
142
143    if (! iscell (y))
144      inpargs = {evtfun, t, y};
145    else
146      inpargs = {evtfun, t, y{1}, y{2}};
147      y = y{1};  # Delete cell element 2
148    endif
149    if (nargin > 4)
150      inpargs = {inpargs{:}, varargin{:}};
151    endif
152    [evtold, ~, ~] = feval (inpargs{:});
153
154    ## We require that all return values be row vectors
155    evtold = evtold(:).'; told = t; yold = y;
156    evtcnt = 1;
157    retval = retcell = cell (1,4);
158
159  elseif (strcmp (flag, "done"))
160    ## Clear this event handling function
161    firstrun = true;
162    evtold = told = yold = evtcnt = [];
163    retval = retcell = cell (1,4);
164
165  endif
166
167endfunction
168