Mercurial > octave-nkf

comparison scripts/ode/private/odepkg_event_handle.m @ 20584:eb9e2d187ed2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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
comparison
equal deleted inserted replaced
20583:d746695bf494 20584:eb9e2d187ed2
15 ## You should have received a copy of the GNU General Public License 15 ## You should have received a copy of the GNU General Public License
16 ## along with Octave; see the file COPYING. If not, see 16 ## along with Octave; see the file COPYING. If not, see
17 ## <http://www.gnu.org/licenses/>. 17 ## <http://www.gnu.org/licenses/>.
18 18
19 ## -*- texinfo -*- 19 ## -*- texinfo -*-
20 ## @deftypefn {Function File} {[@var{sol}] =} odepkg_event_handle (@var{@@fun}, @var{time}, @var{y}, @var{flag}, @var{par1}, @var{par2}, @dots{}) 20 ## @deftypefn {Function File} {@var{sol} =} odepkg_event_handle (@var{@@fun}, @var{time}, @var{y}, @var{flag}, @var{par1}, @var{par2}, @dots{})
21 ## 21 ##
22 ## Return the solution of the event function that is specified as the first 22 ## Return the solution of the event function that is specified as the first
23 ## input argument @var{@@fun} in form of a function handle. 23 ## input argument @var{@@fun} in the form of a function handle.
24 ## 24 ##
25 ## The second input argument @var{time} is of type double scalar and 25 ## The second input argument @var{time} is of type double scalar and
26 ## specifies the time of the event evaluation, the third input argument 26 ## specifies the time of the event evaluation, the third input argument
27 ## @var{y} either is of type double column vector (for ODEs and DAEs) and 27 ## @var{y} either is of type double column vector (for ODEs and DAEs) and
28 ## specifies the solutions or is of type cell array (for IDEs and DDEs) and 28 ## specifies the solutions or is of type cell array (for IDEs and DDEs) and
53 ## achieve the highest performance. 53 ## achieve the highest performance.
54 ## @end deftypefn 54 ## @end deftypefn
55 ## 55 ##
56 ## @seealso{odepkg} 56 ## @seealso{odepkg}
57 57
58 function [vretval] = odepkg_event_handle (vevefun, vt, vy, vflag, varargin) 58 function vretval = odepkg_event_handle (vevefun, vt, vy, vflag, varargin)
59 59
60 ## No error handling has been implemented in this function to achieve 60 ## No error handling has been implemented in this function to achieve
61 ## the highest performance available. 61 ## the highest performance available.
62 62
63 ## vretval{1} is true or false; either to terminate or to continue 63 ## vretval{1} is true or false; either to terminate or to continue
76 persistent vevecnt; 76 persistent vevecnt;
77 77
78 ## Call the event function if an event function has been defined to 78 ## Call the event function if an event function has been defined to
79 ## initialize the internal variables of the event function an to get 79 ## initialize the internal variables of the event function an to get
80 ## a value for veveold 80 ## a value for veveold
81 if (strcmp (vflag, 'init')) 81 if (strcmp (vflag, "init"))
82 82
83 if (~iscell (vy)) 83 if (! iscell (vy))
84 vinpargs = {vevefun, vt, vy}; 84 vinpargs = {vevefun, vt, vy};
85 else 85 else
86 vinpargs = {vevefun, vt, vy{1}, vy{2}}; 86 vinpargs = {vevefun, vt, vy{1}, vy{2}};
87 vy = vy{1}; ## Delete cell element 2 87 vy = vy{1}; # Delete cell element 2
88 end 88 endif
89 if (nargin > 4) 89 if (nargin > 4)
90 vinpargs = {vinpargs{:}, varargin{:}}; 90 vinpargs = {vinpargs{:}, varargin{:}};
91 end 91 endif
92 [veveold, vterm, vdir] = feval (vinpargs{:}); 92 [veveold, vterm, vdir] = feval (vinpargs{:});
93 93
94 ## We assume that all return values must be column vectors 94 ## We assume that all return values must be column vectors
95 veveold = veveold(:)'; vterm = vterm(:)'; vdir = vdir(:)'; 95 veveold = veveold(:)'; vterm = vterm(:)'; vdir = vdir(:)';
96 vtold = vt; vyold = vy; vevecnt = 1; vretcell = cell (1,4); 96 vtold = vt; vyold = vy; vevecnt = 1; vretcell = cell (1,4);
97 97
98 ## Process the event, find the zero crossings either for a rising 98 ## Process the event, find the zero crossings either for a rising
99 ## or for a falling edge 99 ## or for a falling edge
100 elseif (isempty (vflag)) 100 elseif (isempty (vflag))
101 101
102 if (~iscell (vy)) 102 if (! iscell (vy))
103 vinpargs = {vevefun, vt, vy}; 103 vinpargs = {vevefun, vt, vy};
104 else 104 else
105 vinpargs = {vevefun, vt, vy{1}, vy{2}}; 105 vinpargs = {vevefun, vt, vy{1}, vy{2}};
106 vy = vy{1}; ## Delete cell element 2 106 vy = vy{1}; # Delete cell element 2
107 end 107 endif
108 if (nargin > 4) 108 if (nargin > 4)
109 vinpargs = {vinpargs{:}, varargin{:}}; 109 vinpargs = {vinpargs{:}, varargin{:}};
110 end 110 endif
111 [veve, vterm, vdir] = feval (vinpargs{:}); 111 [veve, vterm, vdir] = feval (vinpargs{:});
112 112
113 ## We assume that all return values must be column vectors 113 ## We assume that all return values must be column vectors
114 veve = veve(:)'; vterm = vterm(:)'; vdir = vdir(:)'; 114 veve = veve(:)'; vterm = vterm(:)'; vdir = vdir(:)';
115 115
116 ## Check if one or more signs of the event has changed 116 ## Check if one or more signs of the event has changed
117 vsignum = (sign (veveold) ~= sign (veve)); 117 vsignum = (sign (veveold) != sign (veve));
118 if (any (vsignum)) ## One or more values have changed 118 if (any (vsignum)) # One or more values have changed
119 vindex = find (vsignum); ## Get the index of the changed values 119 vindex = find (vsignum); # Get the index of the changed values
120 120
121 if (any (vdir(vindex) == 0)) 121 if (any (vdir(vindex) == 0))
122 ## Rising or falling (both are possible) 122 ## Rising or falling (both are possible)
123 ## Don't change anything, keep the index 123 ## Don't change anything, keep the index
124 elseif (any (vdir(vindex) == sign (veve(vindex)))) 124 elseif (any (vdir(vindex) == sign (veve(vindex))))
125 ## Detected rising or falling, need a new index 125 ## Detected rising or falling, need a new index
126 vindex = find (vdir == sign (veve)); 126 vindex = find (vdir == sign (veve));
127 else 127 else
128 ## Found a zero crossing but must not be notified 128 ## Found a zero crossing but must not be notified
129 vindex = []; 129 vindex = [];
130 end 130 endif
131 131
132 ## Create new output values if a valid index has been found 132 ## Create new output values if a valid index has been found
133 if (~isempty (vindex)) 133 if (! isempty (vindex))
134 ## Change the persistent result cell array 134 ## Change the persistent result cell array
135 vretcell{1} = any (vterm(vindex)); ## Stop integration or not 135 vretcell{1} = any (vterm(vindex)); # Stop integration or not
136 vretcell{2}(vevecnt,1) = vindex(1,1); ## Take first event found 136 vretcell{2}(vevecnt,1) = vindex(1,1); # Take first event found
137 ## Calculate the time stamp when the event function returned 0 and 137 ## Calculate the time stamp when the event function returned 0 and
138 ## calculate new values for the integration results, we do both by 138 ## calculate new values for the integration results, we do both by
139 ## a linear interpolation 139 ## a linear interpolation
140 vtnew = vt - veve(1,vindex) * (vt - vtold) / (veve(1,vindex) - veveold(1,vindex)); 140 vtnew = vt - veve(1,vindex) * (vt - vtold) / ...
141 (veve(1,vindex) - veveold(1,vindex));
141 vynew = (vy - (vt - vtnew) * (vy - vyold) / (vt - vtold))'; 142 vynew = (vy - (vt - vtnew) * (vy - vyold) / (vt - vtold))';
142 vretcell{3}(vevecnt,1) = vtnew; 143 vretcell{3}(vevecnt,1) = vtnew;
143 vretcell{4}(vevecnt,:) = vynew; 144 vretcell{4}(vevecnt,:) = vynew;
144 vevecnt = vevecnt + 1; 145 vevecnt = vevecnt + 1;
145 end ## if (~isempty (vindex)) 146 endif
146 147
147 end ## Check for one or more signs ... 148 endif # Check for one or more signs ...
148 veveold = veve; vtold = vt; vretval = vretcell; vyold = vy; 149 veveold = veve; vtold = vt; vretval = vretcell; vyold = vy;
149 150
150 elseif (strcmp (vflag, 'done')) ## Clear this event handling function 151 elseif (strcmp (vflag, "done")) # Clear this event handling function
151 clear ('veveold', 'vtold', 'vretcell', 'vyold', 'vevecnt'); 152 clear ("veveold", "vtold", "vretcell", "vyold", "vevecnt");
152 vretcell = cell (1,4); 153 vretcell = cell (1,4);
153 154
154 end 155 endif
156 endfunction
157