1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math.ode.events;
19
20 import org.apache.commons.math.ConvergenceException;
21 import org.apache.commons.math.FunctionEvaluationException;
22 import org.apache.commons.math.analysis.UnivariateRealFunction;
23 import org.apache.commons.math.analysis.solvers.BrentSolver;
24 import org.apache.commons.math.ode.DerivativeException;
25 import org.apache.commons.math.ode.sampling.StepInterpolator;
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41 public class EventState {
42
43
44 private final EventHandler handler;
45
46
47 private final double maxCheckInterval;
48
49
50 private final double convergence;
51
52
53 private final int maxIterationCount;
54
55
56 private double t0;
57
58
59 private double g0;
60
61
62 private boolean g0Positive;
63
64
65 private boolean pendingEvent;
66
67
68 private double pendingEventTime;
69
70
71 private double previousEventTime;
72
73
74 private boolean forward;
75
76
77
78
79 private boolean increasing;
80
81
82 private int nextAction;
83
84
85
86
87
88
89
90
91
92
93 public EventState(final EventHandler handler, final double maxCheckInterval,
94 final double convergence, final int maxIterationCount) {
95 this.handler = handler;
96 this.maxCheckInterval = maxCheckInterval;
97 this.convergence = Math.abs(convergence);
98 this.maxIterationCount = maxIterationCount;
99
100
101 t0 = Double.NaN;
102 g0 = Double.NaN;
103 g0Positive = true;
104 pendingEvent = false;
105 pendingEventTime = Double.NaN;
106 previousEventTime = Double.NaN;
107 increasing = true;
108 nextAction = EventHandler.CONTINUE;
109
110 }
111
112
113
114
115 public EventHandler getEventHandler() {
116 return handler;
117 }
118
119
120
121
122 public double getMaxCheckInterval() {
123 return maxCheckInterval;
124 }
125
126
127
128
129 public double getConvergence() {
130 return convergence;
131 }
132
133
134
135
136 public int getMaxIterationCount() {
137 return maxIterationCount;
138 }
139
140
141
142
143
144
145
146
147
148 public void reinitializeBegin(final double t0, final double[] y0)
149 throws EventException {
150 this.t0 = t0;
151 g0 = handler.g(t0, y0);
152 g0Positive = (g0 >= 0);
153 }
154
155
156
157
158
159
160
161
162
163
164
165
166 public boolean evaluateStep(final StepInterpolator interpolator)
167 throws DerivativeException, EventException, ConvergenceException {
168
169 try {
170
171 forward = interpolator.isForward();
172 final double t1 = interpolator.getCurrentTime();
173 final int n = Math.max(1, (int) Math.ceil(Math.abs(t1 - t0) / maxCheckInterval));
174 final double h = (t1 - t0) / n;
175
176 double ta = t0;
177 double ga = g0;
178 double tb = t0 + (interpolator.isForward() ? convergence : -convergence);
179 for (int i = 0; i < n; ++i) {
180
181
182 tb += h;
183 interpolator.setInterpolatedTime(tb);
184 final double gb = handler.g(tb, interpolator.getInterpolatedState());
185
186
187 if (g0Positive ^ (gb >= 0)) {
188
189
190
191 increasing = (gb >= ga);
192
193 final UnivariateRealFunction f = new UnivariateRealFunction() {
194 public double value(final double t) throws FunctionEvaluationException {
195 try {
196 interpolator.setInterpolatedTime(t);
197 return handler.g(t, interpolator.getInterpolatedState());
198 } catch (DerivativeException e) {
199 throw new FunctionEvaluationException(e, t);
200 } catch (EventException e) {
201 throw new FunctionEvaluationException(e, t);
202 }
203 }
204 };
205 final BrentSolver solver = new BrentSolver();
206 solver.setAbsoluteAccuracy(convergence);
207 solver.setMaximalIterationCount(maxIterationCount);
208 double root;
209 try {
210 root = (ta <= tb) ? solver.solve(f, ta, tb) : solver.solve(f, tb, ta);
211 } catch (IllegalArgumentException iae) {
212
213 root = Double.NaN;
214 }
215 if (Double.isNaN(root) ||
216 ((Math.abs(root - ta) <= convergence) &&
217 (Math.abs(root - previousEventTime) <= convergence))) {
218
219 ta = tb;
220 ga = gb;
221 } else if (Double.isNaN(previousEventTime) ||
222 (Math.abs(previousEventTime - root) > convergence)) {
223 pendingEventTime = root;
224 if (pendingEvent && (Math.abs(t1 - pendingEventTime) <= convergence)) {
225
226
227
228
229 return false;
230 }
231
232
233 pendingEvent = true;
234 return true;
235 }
236
237 } else {
238
239 ta = tb;
240 ga = gb;
241 }
242
243 }
244
245
246 pendingEvent = false;
247 pendingEventTime = Double.NaN;
248 return false;
249
250 } catch (FunctionEvaluationException e) {
251 final Throwable cause = e.getCause();
252 if ((cause != null) && (cause instanceof DerivativeException)) {
253 throw (DerivativeException) cause;
254 } else if ((cause != null) && (cause instanceof EventException)) {
255 throw (EventException) cause;
256 }
257 throw new EventException(e);
258 }
259
260 }
261
262
263
264
265
266
267 public double getEventTime() {
268 return pendingEventTime;
269 }
270
271
272
273
274
275
276
277
278
279 public void stepAccepted(final double t, final double[] y)
280 throws EventException {
281
282 t0 = t;
283 g0 = handler.g(t, y);
284
285 if (pendingEvent) {
286
287 previousEventTime = t;
288 g0Positive = increasing;
289 nextAction = handler.eventOccurred(t, y, !(increasing ^ forward));
290 } else {
291 g0Positive = (g0 >= 0);
292 nextAction = EventHandler.CONTINUE;
293 }
294 }
295
296
297
298
299
300 public boolean stop() {
301 return nextAction == EventHandler.STOP;
302 }
303
304
305
306
307
308
309
310
311
312
313 public boolean reset(final double t, final double[] y)
314 throws EventException {
315
316 if (! pendingEvent) {
317 return false;
318 }
319
320 if (nextAction == EventHandler.RESET_STATE) {
321 handler.resetState(t, y);
322 }
323 pendingEvent = false;
324 pendingEventTime = Double.NaN;
325
326 return (nextAction == EventHandler.RESET_STATE) ||
327 (nextAction == EventHandler.RESET_DERIVATIVES);
328
329 }
330
331 }