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