Mercurial > octave
view scripts/ode/private/odepkg_event_handle.m @ 20552:eb9e2d187ed2
maint: Use Octave coding conventions in scripts/ode/private dir.
* AbsRel_Norm.m, fuzzy_compare.m, hermite_quartic_interpolation.m,
integrate_adaptive.m, integrate_const.m, integrate_n_steps.m, kahan.m,
ode_struct_value_check.m, odepkg_event_handle.m, odepkg_structure_check.m,
runge_kutta_45_dorpri.m, starting_stepsize.m:
Wrap long lines to < 80 chars.
Use double quotes rather than single quotes where possible.
Use ';' at end of keywords "return;" and "break;"
Use '##" for stand-alone comments and '#' for end-of-line comments.
Use two spaces after period before starting new sentence.
Use '!' instead of '~' for logical negation.
Use specific form of end (endif, endfor, etc.).
Don't use line continuation marker '...' unless necessary.
author | Rik <rik@octave.org> |
---|---|
date | Sun, 04 Oct 2015 22:18:54 -0700 |
parents | 25623ef2ff4f |
children | b7ac1e94266e |
line wrap: on
line source
## Copyright (C) 2006-2012, Thomas Treichl <treichl@users.sourceforge.net> ## ## This file is part of Octave. ## ## Octave is free software; you can redistribute it and/or modify it ## under the terms of the GNU General Public License as published by ## the Free Software Foundation; either version 3 of the License, or (at ## your option) any later version. ## ## Octave is distributed in the hope that it will be useful, but ## WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ## General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with Octave; see the file COPYING. If not, see ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {Function File} {@var{sol} =} odepkg_event_handle (@var{@@fun}, @var{time}, @var{y}, @var{flag}, @var{par1}, @var{par2}, @dots{}) ## ## Return the solution of the event function that is specified as the first ## input argument @var{@@fun} in the form of a function handle. ## ## The second input argument @var{time} is of type double scalar and ## specifies the time of the event evaluation, the third input argument ## @var{y} either is of type double column vector (for ODEs and DAEs) and ## specifies the solutions or is of type cell array (for IDEs and DDEs) and ## specifies the derivatives or the history values, the third input argument ## @var{flag} is of type string and can be of the form ## ## @table @option ## @item @qcode{"init"} ## then initialize internal persistent variables of the function ## @command{odepkg_event_handle} and return an empty cell array of size 4, ## ## @item @qcode{"calc"} ## then do the evaluation of the event function and return the solution ## @var{sol} as type cell array of size 4, ## ## @item @qcode{"done"} ## then cleanup internal variables of the function ## @command{odepkg_event_handle} and return an empty cell array of size 4. ## @end table ## ## Optionally if further input arguments @var{par1}, @var{par2}, @dots{} of ## any type are given then pass these parameters through ## @command{odepkg_event_handle} to the event function. ## ## This function is an OdePkg internal helper function therefore it should ## never be necessary that this function is called directly by a user. There ## is only little error detection implemented in this function file to ## achieve the highest performance. ## @end deftypefn ## ## @seealso{odepkg} function vretval = odepkg_event_handle (vevefun, vt, vy, vflag, varargin) ## No error handling has been implemented in this function to achieve ## the highest performance available. ## vretval{1} is true or false; either to terminate or to continue ## vretval{2} is the index information for which event occured ## vretval{3} is the time information column vector ## vretval{4} is the line by line result information matrix ## These persistent variables are needed to store the results and the ## time value from the processing in the time stamp before, veveold ## are the results from the event function, vtold the time stamp, ## vretcell the return values cell array, vyold the result of the ode ## and vevecnt the counter for how often this event handling ## has been called persistent veveold; persistent vtold; persistent vretcell; persistent vyold; persistent vevecnt; ## Call the event function if an event function has been defined to ## initialize the internal variables of the event function an to get ## a value for veveold if (strcmp (vflag, "init")) if (! iscell (vy)) vinpargs = {vevefun, vt, vy}; else vinpargs = {vevefun, vt, vy{1}, vy{2}}; vy = vy{1}; # Delete cell element 2 endif if (nargin > 4) vinpargs = {vinpargs{:}, varargin{:}}; endif [veveold, vterm, vdir] = feval (vinpargs{:}); ## We assume that all return values must be column vectors veveold = veveold(:)'; vterm = vterm(:)'; vdir = vdir(:)'; vtold = vt; vyold = vy; vevecnt = 1; vretcell = cell (1,4); ## Process the event, find the zero crossings either for a rising ## or for a falling edge elseif (isempty (vflag)) if (! iscell (vy)) vinpargs = {vevefun, vt, vy}; else vinpargs = {vevefun, vt, vy{1}, vy{2}}; vy = vy{1}; # Delete cell element 2 endif if (nargin > 4) vinpargs = {vinpargs{:}, varargin{:}}; endif [veve, vterm, vdir] = feval (vinpargs{:}); ## We assume that all return values must be column vectors veve = veve(:)'; vterm = vterm(:)'; vdir = vdir(:)'; ## Check if one or more signs of the event has changed vsignum = (sign (veveold) != sign (veve)); if (any (vsignum)) # One or more values have changed vindex = find (vsignum); # Get the index of the changed values if (any (vdir(vindex) == 0)) ## Rising or falling (both are possible) ## Don't change anything, keep the index elseif (any (vdir(vindex) == sign (veve(vindex)))) ## Detected rising or falling, need a new index vindex = find (vdir == sign (veve)); else ## Found a zero crossing but must not be notified vindex = []; endif ## Create new output values if a valid index has been found if (! isempty (vindex)) ## Change the persistent result cell array vretcell{1} = any (vterm(vindex)); # Stop integration or not vretcell{2}(vevecnt,1) = vindex(1,1); # Take first event found ## 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 vtnew = vt - veve(1,vindex) * (vt - vtold) / ... (veve(1,vindex) - veveold(1,vindex)); vynew = (vy - (vt - vtnew) * (vy - vyold) / (vt - vtold))'; vretcell{3}(vevecnt,1) = vtnew; vretcell{4}(vevecnt,:) = vynew; vevecnt = vevecnt + 1; endif endif # Check for one or more signs ... veveold = veve; vtold = vt; vretval = vretcell; vyold = vy; elseif (strcmp (vflag, "done")) # Clear this event handling function clear ("veveold", "vtold", "vretcell", "vyold", "vevecnt"); vretcell = cell (1,4); endif endfunction