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