1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.hipparchus.ode;
19
20 import java.util.Arrays;
21 import java.util.List;
22
23 import org.hipparchus.UnitTestUtils;
24 import org.hipparchus.exception.MathIllegalArgumentException;
25 import org.hipparchus.exception.MathIllegalStateException;
26 import org.hipparchus.ode.VariationalEquation.MismatchedEquations;
27 import org.hipparchus.ode.nonstiff.DormandPrince54Integrator;
28 import org.hipparchus.util.FastMath;
29 import org.junit.Assert;
30 import org.junit.Test;
31
32 public class VariationalEquationTest {
33
34 @Test
35 public void testLowAccuracyExternalDifferentiation()
36 throws MathIllegalArgumentException, MathIllegalStateException {
37
38
39
40
41
42
43
44
45
46 ODEIntegrator integ =
47 new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-4, 1.0e-4 }, new double[] { 1.0e-4, 1.0e-4 });
48 double hP = 1.0e-12;
49 UnitTestUtils.SimpleStatistics residualsP0 = new UnitTestUtils.SimpleStatistics();
50 UnitTestUtils.SimpleStatistics residualsP1 = new UnitTestUtils.SimpleStatistics();
51 for (double b = 2.88; b < 3.08; b += 0.001) {
52 Brusselator brusselator = new Brusselator(b);
53 double[] y = { 1.3, b };
54 y = integ.integrate(brusselator, new ODEState(0, y), 20.0).getPrimaryState();
55 double[] yP = { 1.3, b + hP };
56 yP = integ.integrate(brusselator, new ODEState(0, yP), 20.0).getPrimaryState();
57 residualsP0.addValue((yP[0] - y[0]) / hP - brusselator.dYdP0());
58 residualsP1.addValue((yP[1] - y[1]) / hP - brusselator.dYdP1());
59 }
60 Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) > 500);
61 Assert.assertTrue(residualsP0.getStandardDeviation() > 30);
62 Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) > 700);
63 Assert.assertTrue(residualsP1.getStandardDeviation() > 40);
64 }
65
66 @Test
67 public void testHighAccuracyExternalDifferentiation()
68 throws MathIllegalArgumentException, MathIllegalStateException {
69 ODEIntegrator integ =
70 new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-10, 1.0e-10 }, new double[] { 1.0e-10, 1.0e-10 });
71 double hP = 1.0e-12;
72 UnitTestUtils.SimpleStatistics residualsP0 = new UnitTestUtils.SimpleStatistics();
73 UnitTestUtils.SimpleStatistics residualsP1 = new UnitTestUtils.SimpleStatistics();
74 for (double b = 2.88; b < 3.08; b += 0.001) {
75 ParamBrusselator brusselator = new ParamBrusselator(b);
76 double[] y = { 1.3, b };
77 y = integ.integrate(brusselator, new ODEState(0, y), 20.0).getPrimaryState();
78 double[] yP = { 1.3, b + hP };
79 brusselator.setParameter("b", b + hP);
80 yP = integ.integrate(brusselator, new ODEState(0, yP), 20.0).getPrimaryState();
81 residualsP0.addValue((yP[0] - y[0]) / hP - brusselator.dYdP0());
82 residualsP1.addValue((yP[1] - y[1]) / hP - brusselator.dYdP1());
83 }
84 Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) > 0.02);
85 Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) < 0.03);
86 Assert.assertTrue(residualsP0.getStandardDeviation() > 0.003);
87 Assert.assertTrue(residualsP0.getStandardDeviation() < 0.004);
88 Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) > 0.04);
89 Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) < 0.05);
90 Assert.assertTrue(residualsP1.getStandardDeviation() > 0.007);
91 Assert.assertTrue(residualsP1.getStandardDeviation() < 0.008);
92 }
93
94 @Test
95 public void testWrongParameterName() {
96 final String name = "an-unknown-parameter";
97 try {
98 ParamBrusselator brusselator = new ParamBrusselator(2.9);
99 brusselator.setParameter(name, 3.0);
100 Assert.fail("an exception should have been thrown");
101 } catch (MathIllegalArgumentException upe) {
102 Assert.assertEquals(LocalizedODEFormats.UNKNOWN_PARAMETER, upe.getSpecifier());
103 Assert.assertEquals(name, upe.getParts()[0]);
104 }
105 }
106
107 @Test
108 public void testMismatchedEquations() {
109 try {
110 ExpandableODE efode = new ExpandableODE(new ParamBrusselator(2.9));
111 ODEJacobiansProvider jode = new Brusselator(2.9);
112 new VariationalEquation(efode, jode);
113 Assert.fail("an exception should have been thrown");
114 } catch (VariationalEquation.MismatchedEquations upe) {
115 Assert.assertEquals(LocalizedODEFormats.UNMATCHED_ODE_IN_EXPANDED_SET,
116 upe.getSpecifier());
117 }
118 }
119
120 @Test
121 public void testDefaultMethods() {
122 final String name = "name";
123 ODEJacobiansProvider jode = new ODEJacobiansProvider() {
124 public int getDimension() {
125 return 1;
126 }
127 public double[] computeDerivatives(double t, double[] y) {
128 return y;
129 }
130 public double[][] computeMainStateJacobian(double t, double[] y, double[] yDot) {
131 return null;
132 }
133 };
134 Assert.assertFalse(jode.isSupported(name));
135 Assert.assertTrue(jode.getParametersNames().isEmpty());
136 try {
137 double t = 0.0;
138 double[] y = new double[] { 0.0 };
139 jode.computeParameterJacobian(t, y, jode.computeDerivatives(t, y), name);
140 Assert.fail("an exception should have been thrown");
141 } catch (MathIllegalArgumentException miae) {
142 Assert.assertEquals(LocalizedODEFormats.UNKNOWN_PARAMETER, miae.getSpecifier());
143 Assert.assertSame(name, miae.getParts()[0]);
144 }
145 }
146
147 @Test
148 public void testInternalDifferentiation()
149 throws MathIllegalArgumentException, MathIllegalStateException, MismatchedEquations {
150 AbstractIntegrator integ =
151 new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-4, 1.0e-4 }, new double[] { 1.0e-4, 1.0e-4 });
152 double hP = 1.0e-12;
153 double hY = 1.0e-12;
154 UnitTestUtils.SimpleStatistics residualsP0 = new UnitTestUtils.SimpleStatistics();
155 UnitTestUtils.SimpleStatistics residualsP1 = new UnitTestUtils.SimpleStatistics();
156 for (double b = 2.88; b < 3.08; b += 0.001) {
157 ParamBrusselator brusselator = new ParamBrusselator(b);
158 brusselator.setParameter(ParamBrusselator.B, b);
159 double[] z = { 1.3, b };
160
161 ExpandableODE efode = new ExpandableODE(brusselator);
162 VariationalEquation jacob = new VariationalEquation(efode, brusselator, new double[] { hY, hY },
163 brusselator,
164 new ParameterConfiguration(ParamBrusselator.B, hP));
165 jacob.setInitialParameterJacobian(ParamBrusselator.B, new double[] { 0.0, 1.0 });
166
167 integ.setMaxEvaluations(5000);
168 final ODEState initialState = jacob.setUpInitialState(new ODEState(0, z));
169 final ODEStateAndDerivative finalState = integ.integrate(efode, initialState, 20.0);
170 final double[] dZdP = jacob.extractParameterJacobian(finalState, ParamBrusselator.B);
171
172
173
174
175 residualsP0.addValue(dZdP[0] - brusselator.dYdP0());
176 residualsP1.addValue(dZdP[1] - brusselator.dYdP1());
177 }
178 Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) < 0.02);
179 Assert.assertTrue(residualsP0.getStandardDeviation() < 0.003);
180 Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) < 0.05);
181 Assert.assertTrue(residualsP1.getStandardDeviation() < 0.01);
182 }
183
184 @Test
185 public void testAnalyticalDifferentiation()
186 throws MathIllegalArgumentException, MathIllegalStateException, MismatchedEquations {
187 AbstractIntegrator integ =
188 new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-4, 1.0e-4 }, new double[] { 1.0e-4, 1.0e-4 });
189 UnitTestUtils.SimpleStatistics residualsP0 = new UnitTestUtils.SimpleStatistics();
190 UnitTestUtils.SimpleStatistics residualsP1 = new UnitTestUtils.SimpleStatistics();
191 for (double b = 2.88; b < 3.08; b += 0.001) {
192 Brusselator brusselator = new Brusselator(b);
193 double[] z = { 1.3, b };
194
195 ExpandableODE efode = new ExpandableODE(brusselator);
196 VariationalEquation jacob = new VariationalEquation(efode, brusselator);
197 jacob.setInitialParameterJacobian(Brusselator.B, new double[] { 0.0, 1.0 });
198
199 integ.setMaxEvaluations(5000);
200 final ODEState initialState = jacob.setUpInitialState(new ODEState(0, z));
201 final ODEStateAndDerivative finalState = integ.integrate(efode, initialState, 20.0);
202 final double[] dZdP = jacob.extractParameterJacobian(finalState, Brusselator.B);
203
204
205
206 residualsP0.addValue(dZdP[0] - brusselator.dYdP0());
207 residualsP1.addValue(dZdP[1] - brusselator.dYdP1());
208 }
209 Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) < 0.014);
210 Assert.assertTrue(residualsP0.getStandardDeviation() < 0.003);
211 Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) < 0.05);
212 Assert.assertTrue(residualsP1.getStandardDeviation() < 0.01);
213 }
214
215 @Test
216 public void testFinalResult()
217 throws MathIllegalArgumentException, MathIllegalStateException, MismatchedEquations {
218
219 AbstractIntegrator integ =
220 new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-10, 1.0e-10 }, new double[] { 1.0e-10, 1.0e-10 });
221 double[] y = new double[] { 0.0, 1.0 };
222 Circle circle = new Circle(y, 1.0, 1.0, 0.1);
223
224 ExpandableODE efode = new ExpandableODE(circle);
225 VariationalEquation jacob = new VariationalEquation(efode, circle);
226 jacob.setInitialMainStateJacobian(circle.exactDyDy0(0));
227 jacob.setInitialParameterJacobian(Circle.CX, circle.exactDyDcx(0));
228 jacob.setInitialParameterJacobian(Circle.CY, circle.exactDyDcy(0));
229 jacob.setInitialParameterJacobian(Circle.OMEGA, circle.exactDyDom(0));
230
231 integ.setMaxEvaluations(5000);
232
233 double t = 18 * FastMath.PI;
234 final ODEState initialState = jacob.setUpInitialState(new ODEState(0, y));
235 final ODEStateAndDerivative finalState = integ.integrate(efode, initialState, t);
236 y = finalState.getPrimaryState();
237 for (int i = 0; i < y.length; ++i) {
238 Assert.assertEquals(circle.exactY(t)[i], y[i], 1.0e-9);
239 }
240
241 double[][] dydy0 = jacob.extractMainSetJacobian(finalState);
242 for (int i = 0; i < dydy0.length; ++i) {
243 for (int j = 0; j < dydy0[i].length; ++j) {
244 Assert.assertEquals(circle.exactDyDy0(t)[i][j], dydy0[i][j], 1.0e-9);
245 }
246 }
247 double[] dydcx = jacob.extractParameterJacobian(finalState, Circle.CX);
248 for (int i = 0; i < dydcx.length; ++i) {
249 Assert.assertEquals(circle.exactDyDcx(t)[i], dydcx[i], 1.0e-7);
250 }
251 double[] dydcy = jacob.extractParameterJacobian(finalState, Circle.CY);
252 for (int i = 0; i < dydcy.length; ++i) {
253 Assert.assertEquals(circle.exactDyDcy(t)[i], dydcy[i], 1.0e-7);
254 }
255 double[] dydom = jacob.extractParameterJacobian(finalState, Circle.OMEGA);
256 for (int i = 0; i < dydom.length; ++i) {
257 Assert.assertEquals(circle.exactDyDom(t)[i], dydom[i], 1.0e-7);
258 }
259 }
260
261 @Test
262 public void testParameterizable()
263 throws MathIllegalArgumentException, MathIllegalStateException, MismatchedEquations {
264
265 AbstractIntegrator integ =
266 new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-10, 1.0e-10 }, new double[] { 1.0e-10, 1.0e-10 });
267 double[] y = new double[] { 0.0, 1.0 };
268 ParameterizedCircle pcircle = new ParameterizedCircle(y, 1.0, 1.0, 0.1);
269
270 double hP = 1.0e-12;
271 double hY = 1.0e-12;
272
273 ExpandableODE efode = new ExpandableODE(pcircle);
274 VariationalEquation jacob = new VariationalEquation(efode, pcircle, new double[] { hY, hY },
275 pcircle,
276 new ParameterConfiguration(ParameterizedCircle.CX, hP),
277 new ParameterConfiguration(ParameterizedCircle.CY, hP),
278 new ParameterConfiguration(ParameterizedCircle.OMEGA, hP));
279 jacob.setInitialMainStateJacobian(pcircle.exactDyDy0(0));
280 jacob.setInitialParameterJacobian(ParameterizedCircle.CX, pcircle.exactDyDcx(0));
281 jacob.setInitialParameterJacobian(ParameterizedCircle.CY, pcircle.exactDyDcy(0));
282 jacob.setInitialParameterJacobian(ParameterizedCircle.OMEGA, pcircle.exactDyDom(0));
283
284 integ.setMaxEvaluations(50000);
285
286 double t = 18 * FastMath.PI;
287 final ODEState initialState = jacob.setUpInitialState(new ODEState(0, y));
288 final ODEStateAndDerivative finalState = integ.integrate(efode, initialState, t);
289 y = finalState.getPrimaryState();
290 for (int i = 0; i < y.length; ++i) {
291 Assert.assertEquals(pcircle.exactY(t)[i], y[i], 1.0e-9);
292 }
293
294 double[][] dydy0 = jacob.extractMainSetJacobian(finalState);
295 for (int i = 0; i < dydy0.length; ++i) {
296 for (int j = 0; j < dydy0[i].length; ++j) {
297 Assert.assertEquals(pcircle.exactDyDy0(t)[i][j], dydy0[i][j], 5.0e-4);
298 }
299 }
300
301 double[] dydp0 = jacob.extractParameterJacobian(finalState, ParameterizedCircle.CX);
302 for (int i = 0; i < dydp0.length; ++i) {
303 Assert.assertEquals(pcircle.exactDyDcx(t)[i], dydp0[i], 5.0e-4);
304 }
305
306 double[] dydp1 = jacob.extractParameterJacobian(finalState, ParameterizedCircle.OMEGA);
307 for (int i = 0; i < dydp1.length; ++i) {
308 Assert.assertEquals(pcircle.exactDyDom(t)[i], dydp1[i], 1.0e-2);
309 }
310 }
311
312 private static class Brusselator implements ODEJacobiansProvider {
313
314 public static final String B = "b";
315
316 private double b;
317
318 public Brusselator(double b) {
319 this.b = b;
320 }
321
322 @Override
323 public int getDimension() {
324 return 2;
325 }
326
327 @Override
328 public double[] computeDerivatives(double t, double[] y) {
329 double prod = y[0] * y[0] * y[1];
330 return new double[] {
331 1 + prod - (b + 1) * y[0],
332 b * y[0] - prod
333 };
334 }
335
336 @Override
337 public double[][] computeMainStateJacobian(double t, double[] y, double[] yDot) {
338 double p = 2 * y[0] * y[1];
339 double y02 = y[0] * y[0];
340 return new double[][] {
341 { p - (1 + b), y02 },
342 { b - p, -y02}
343 };
344 }
345
346 @Override
347 public List<String> getParametersNames() {
348 return Arrays.asList(B);
349 }
350
351 @Override
352 public boolean isSupported(String name) {
353 return B.equals(name);
354 }
355
356 @Override
357 public double[] computeParameterJacobian(double t, double[] y, double[] yDot,
358 String paramName) {
359 if (isSupported(paramName)) {
360 return new double[] { -y[0], y[0] };
361 } else {
362 throw new MathIllegalArgumentException(LocalizedODEFormats.UNKNOWN_PARAMETER,
363 paramName);
364 }
365 }
366
367 public double dYdP0() {
368 return -1088.232716447743 + (1050.775747149553 + (-339.012934631828 + 36.52917025056327 * b) * b) * b;
369 }
370
371 public double dYdP1() {
372 return 1502.824469929139 + (-1438.6974831849952 + (460.959476642384 - 49.43847385647082 * b) * b) * b;
373 }
374
375 }
376
377 private static class ParamBrusselator extends AbstractParameterizable
378 implements OrdinaryDifferentialEquation, ParametersController {
379
380 public static final String B = "b";
381
382 private double b;
383
384 public ParamBrusselator(double b) {
385 super(B);
386 this.b = b;
387 }
388
389 @Override
390 public int getDimension() {
391 return 2;
392 }
393
394
395 @Override
396 public double getParameter(final String name)
397 throws MathIllegalArgumentException {
398 complainIfNotSupported(name);
399 return b;
400 }
401
402
403 @Override
404 public void setParameter(final String name, final double value)
405 throws MathIllegalArgumentException {
406 complainIfNotSupported(name);
407 b = value;
408 }
409
410 @Override
411 public double[] computeDerivatives(double t, double[] y) {
412 double prod = y[0] * y[0] * y[1];
413 return new double[] {
414 1 + prod - (b + 1) * y[0],
415 b * y[0] - prod
416 };
417 }
418
419 public double dYdP0() {
420 return -1088.232716447743 + (1050.775747149553 + (-339.012934631828 + 36.52917025056327 * b) * b) * b;
421 }
422
423 public double dYdP1() {
424 return 1502.824469929139 + (-1438.6974831849952 + (460.959476642384 - 49.43847385647082 * b) * b) * b;
425 }
426
427 }
428
429
430 private static class Circle implements ODEJacobiansProvider {
431
432 public static final String CX = "cx";
433 public static final String CY = "cy";
434 public static final String OMEGA = "omega";
435
436 private final double[] y0;
437 private double cx;
438 private double cy;
439 private double omega;
440
441 public Circle(double[] y0, double cx, double cy, double omega) {
442 this.y0 = y0.clone();
443 this.cx = cx;
444 this.cy = cy;
445 this.omega = omega;
446 }
447
448 @Override
449 public int getDimension() {
450 return 2;
451 }
452
453 @Override
454 public double[] computeDerivatives(double t, double[] y) {
455 return new double[] {
456 omega * (cy - y[1]),
457 omega * (y[0] - cx)
458 };
459 }
460
461 @Override
462 public double[][] computeMainStateJacobian(double t, double[] y, double[] yDot) {
463 return new double[][] {
464 { 0, -omega },
465 { omega, 0 }
466 };
467 }
468
469 @Override
470 public List<String> getParametersNames() {
471 return Arrays.asList(CX, CY, OMEGA);
472 }
473
474 @Override
475 public boolean isSupported(String name) {
476 return CX.equals(name) || CY.equals(name) || OMEGA.equals(name);
477 }
478
479 @Override
480 public double[] computeParameterJacobian(double t, double[] y, double[] yDot, String paramName)
481 throws MathIllegalArgumentException {
482 if (CX.equals(paramName)) {
483 return new double[] { 0, -omega };
484 } else if (CY.equals(paramName)) {
485 return new double[] { omega, 0 };
486 } else if (OMEGA.equals(paramName)) {
487 return new double[] { cy - y[1], y[0] - cx };
488 } else {
489 throw new MathIllegalArgumentException(LocalizedODEFormats.UNKNOWN_PARAMETER,
490 paramName);
491 }
492 }
493
494 public double[] exactY(double t) {
495 double cos = FastMath.cos(omega * t);
496 double sin = FastMath.sin(omega * t);
497 double dx0 = y0[0] - cx;
498 double dy0 = y0[1] - cy;
499 return new double[] {
500 cx + cos * dx0 - sin * dy0,
501 cy + sin * dx0 + cos * dy0
502 };
503 }
504
505 public double[][] exactDyDy0(double t) {
506 double cos = FastMath.cos(omega * t);
507 double sin = FastMath.sin(omega * t);
508 return new double[][] {
509 { cos, -sin },
510 { sin, cos }
511 };
512 }
513
514 public double[] exactDyDcx(double t) {
515 double cos = FastMath.cos(omega * t);
516 double sin = FastMath.sin(omega * t);
517 return new double[] {1 - cos, -sin};
518 }
519
520 public double[] exactDyDcy(double t) {
521 double cos = FastMath.cos(omega * t);
522 double sin = FastMath.sin(omega * t);
523 return new double[] {sin, 1 - cos};
524 }
525
526 public double[] exactDyDom(double t) {
527 double cos = FastMath.cos(omega * t);
528 double sin = FastMath.sin(omega * t);
529 double dx0 = y0[0] - cx;
530 double dy0 = y0[1] - cy;
531 return new double[] { -t * (sin * dx0 + cos * dy0) , t * (cos * dx0 - sin * dy0) };
532 }
533
534 }
535
536
537 private static class ParameterizedCircle extends AbstractParameterizable
538 implements OrdinaryDifferentialEquation, ParametersController {
539
540 public static final String CX = "cx";
541 public static final String CY = "cy";
542 public static final String OMEGA = "omega";
543
544 private final double[] y0;
545 private double cx;
546 private double cy;
547 private double omega;
548
549 public ParameterizedCircle(double[] y0, double cx, double cy, double omega) {
550 super(CX,CY,OMEGA);
551 this.y0 = y0.clone();
552 this.cx = cx;
553 this.cy = cy;
554 this.omega = omega;
555 }
556
557 @Override
558 public int getDimension() {
559 return 2;
560 }
561
562 @Override
563 public double[] computeDerivatives(double t, double[] y) {
564 return new double[] {
565 omega * (cy - y[1]),
566 omega * (y[0] - cx)
567 };
568 }
569
570 @Override
571 public double getParameter(final String name)
572 throws MathIllegalArgumentException {
573 if (name.equals(CX)) {
574 return cx;
575 } else if (name.equals(CY)) {
576 return cy;
577 } else if (name.equals(OMEGA)) {
578 return omega;
579 } else {
580 throw new MathIllegalArgumentException(LocalizedODEFormats.UNKNOWN_PARAMETER, name);
581 }
582 }
583
584 @Override
585 public void setParameter(final String name, final double value)
586 throws MathIllegalArgumentException {
587 if (name.equals(CX)) {
588 cx = value;
589 } else if (name.equals(CY)) {
590 cy = value;
591 } else if (name.equals(OMEGA)) {
592 omega = value;
593 } else {
594 throw new MathIllegalArgumentException(LocalizedODEFormats.UNKNOWN_PARAMETER, name);
595 }
596 }
597
598 public double[] exactY(double t) {
599 double cos = FastMath.cos(omega * t);
600 double sin = FastMath.sin(omega * t);
601 double dx0 = y0[0] - cx;
602 double dy0 = y0[1] - cy;
603 return new double[] {
604 cx + cos * dx0 - sin * dy0,
605 cy + sin * dx0 + cos * dy0
606 };
607 }
608
609 public double[][] exactDyDy0(double t) {
610 double cos = FastMath.cos(omega * t);
611 double sin = FastMath.sin(omega * t);
612 return new double[][] {
613 { cos, -sin },
614 { sin, cos }
615 };
616 }
617
618 public double[] exactDyDcx(double t) {
619 double cos = FastMath.cos(omega * t);
620 double sin = FastMath.sin(omega * t);
621 return new double[] {1 - cos, -sin};
622 }
623
624 public double[] exactDyDcy(double t) {
625 double cos = FastMath.cos(omega * t);
626 double sin = FastMath.sin(omega * t);
627 return new double[] {sin, 1 - cos};
628 }
629
630 public double[] exactDyDom(double t) {
631 double cos = FastMath.cos(omega * t);
632 double sin = FastMath.sin(omega * t);
633 double dx0 = y0[0] - cx;
634 double dy0 = y0[1] - cy;
635 return new double[] { -t * (sin * dx0 + cos * dy0) , t * (cos * dx0 - sin * dy0) };
636 }
637
638 }
639
640 }