geo/algorithm/geodesic_area.rs
1use crate::geometry::*;
2use geographiclib_rs::{Geodesic, PolygonArea, Winding};
3
4/// Determine the perimeter and area of a geometry on an ellipsoidal model of the earth.
5///
6/// This uses the geodesic measurement methods given by [Karney (2013)].
7///
8/// [Karney (2013)]: https://arxiv.org/pdf/1109.4448.pdf
9pub trait GeodesicArea<T> {
10 /// Determine the area of a geometry on an ellipsoidal model of the earth.
11 ///
12 /// This uses the geodesic measurement methods given by [Karney (2013)].
13 ///
14 /// # Assumptions
15 /// - Polygons are assumed to be wound in a counter-clockwise direction
16 /// for the exterior ring and a clockwise direction for interior rings.
17 /// This is the standard winding for geometries that follow the Simple Feature standard.
18 /// Alternative windings may result in a negative area. See "Interpreting negative area values" below.
19 /// - Polygons are assumed to be smaller than half the size of the earth. If you expect to be dealing
20 /// with polygons larger than this, please use the `unsigned` methods.
21 ///
22 /// # Units
23 ///
24 /// - return value: meter²
25 ///
26 /// # Interpreting negative area values
27 ///
28 /// A negative value can mean one of two things:
29 /// 1. The winding of the polygon is in the clockwise direction (reverse winding). If this is the case, and you know the polygon is smaller than half the area of earth, you can take the absolute value of the reported area to get the correct area.
30 /// 2. The polygon is larger than half the planet. In this case, the returned area of the polygon is not correct. If you expect to be dealing with very large polygons, please use the `unsigned` methods.
31 ///
32 /// # Examples
33 /// ```rust
34 /// use geo::prelude::*;
35 /// use geo::polygon;
36 /// use geo::Polygon;
37 ///
38 /// // The O2 in London
39 /// let mut polygon: Polygon<f64> = polygon![
40 /// (x: 0.00388383, y: 51.501574),
41 /// (x: 0.00538587, y: 51.502278),
42 /// (x: 0.00553607, y: 51.503299),
43 /// (x: 0.00467777, y: 51.504181),
44 /// (x: 0.00327229, y: 51.504435),
45 /// (x: 0.00187754, y: 51.504168),
46 /// (x: 0.00087976, y: 51.503380),
47 /// (x: 0.00107288, y: 51.502324),
48 /// (x: 0.00185608, y: 51.501770),
49 /// (x: 0.00388383, y: 51.501574),
50 /// ];
51 ///
52 /// let area = polygon.geodesic_area_unsigned();
53 ///
54 /// assert_eq!(
55 /// 78_596., // meters
56 /// area.round()
57 /// );
58 /// ```
59 /// [Karney (2013)]: https://arxiv.org/pdf/1109.4448.pdf
60 fn geodesic_area_signed(&self) -> T;
61
62 /// Determine the area of a geometry on an ellipsoidal model of the earth. Supports very large geometries that cover a significant portion of the earth.
63 ///
64 /// This uses the geodesic measurement methods given by [Karney (2013)].
65 ///
66 /// # Assumptions
67 /// - Polygons are assumed to be wound in a counter-clockwise direction
68 /// for the exterior ring and a clockwise direction for interior rings.
69 /// This is the standard winding for geometries that follow the Simple Features standard.
70 /// Using alternative windings will result in incorrect results.
71 ///
72 /// # Units
73 ///
74 /// - return value: meter²
75 ///
76 /// # Examples
77 /// ```rust
78 /// use geo::prelude::*;
79 /// use geo::polygon;
80 /// use geo::Polygon;
81 ///
82 /// // Describe a polygon that covers all of the earth EXCEPT this small square.
83 /// // The outside of the polygon is in this square, the inside of the polygon is the rest of the earth.
84 /// let mut polygon: Polygon<f64> = polygon![
85 /// (x: 0.0, y: 0.0),
86 /// (x: 0.0, y: 1.0),
87 /// (x: 1.0, y: 1.0),
88 /// (x: 1.0, y: 0.0),
89 /// ];
90 ///
91 /// let area = polygon.geodesic_area_unsigned();
92 ///
93 /// // Over 5 trillion square meters!
94 /// assert_eq!(area, 510053312945726.94);
95 /// ```
96 /// [Karney (2013)]: https://arxiv.org/pdf/1109.4448.pdf
97 fn geodesic_area_unsigned(&self) -> T;
98
99 /// Determine the perimeter of a geometry on an ellipsoidal model of the earth.
100 ///
101 /// This uses the geodesic measurement methods given by [Karney (2013)].
102 ///
103 /// For a polygon this returns the sum of the perimeter of the exterior ring and interior rings.
104 /// To get the perimeter of just the exterior ring of a polygon, do `polygon.exterior().geodesic_length()`.
105 ///
106 /// # Units
107 ///
108 /// - return value: meter
109 ///
110 /// [Karney (2013)]: https://arxiv.org/pdf/1109.4448.pdf
111 fn geodesic_perimeter(&self) -> T;
112
113 /// Determine the perimeter and area of a geometry on an ellipsoidal model of the earth, all in one operation.
114 ///
115 /// This returns the perimeter and area in a `(perimeter, area)` tuple and uses the geodesic measurement methods given by [Karney (2013)].
116 ///
117 /// # Area Assumptions
118 /// - Polygons are assumed to be wound in a counter-clockwise direction
119 /// for the exterior ring and a clockwise direction for interior rings.
120 /// This is the standard winding for Geometries that follow the Simple Features standard.
121 /// Alternative windings may result in a negative area. See "Interpreting negative area values" below.
122 /// - Polygons are assumed to be smaller than half the size of the earth. If you expect to be dealing
123 /// with polygons larger than this, please use the 'unsigned' methods.
124 ///
125 /// # Perimeter
126 /// For a polygon this returns the sum of the perimeter of the exterior ring and interior rings.
127 /// To get the perimeter of just the exterior ring of a polygon, do `polygon.exterior().geodesic_length()`.
128 ///
129 /// # Units
130 ///
131 /// - return value: (meter, meter²)
132 ///
133 /// # Interpreting negative area values
134 ///
135 /// A negative area value can mean one of two things:
136 /// 1. The winding of the polygon is in the clockwise direction (reverse winding). If this is the case, and you know the polygon is smaller than half the area of earth, you can take the absolute value of the reported area to get the correct area.
137 /// 2. The polygon is larger than half the planet. In this case, the returned area of the polygon is not correct. If you expect to be dealing with very large polygons, please use the 'unsigned' methods.
138 ///
139 /// [Karney (2013)]: https://arxiv.org/pdf/1109.4448.pdf
140 fn geodesic_perimeter_area_signed(&self) -> (T, T);
141
142 /// Determine the perimeter and area of a geometry on an ellipsoidal model of the earth, all in one operation. Supports very large geometries that cover a significant portion of the earth.
143 ///
144 /// This returns the perimeter and area in a `(perimeter, area)` tuple and uses the geodesic measurement methods given by [Karney (2013)].
145 ///
146 /// # Area Assumptions
147 /// - Polygons are assumed to be wound in a counter-clockwise direction
148 /// for the exterior ring and a clockwise direction for interior rings.
149 /// This is the standard winding for Geometries that follow the Simple Features standard.
150 /// Using alternative windings will result in incorrect results.
151 ///
152 /// # Perimeter
153 /// For a polygon this returns the perimeter of the exterior ring and interior rings.
154 /// To get the perimeter of just the exterior ring of a polygon, do `polygon.exterior().geodesic_length()`.
155 ///
156 /// # Units
157 ///
158 /// - return value: (meter, meter²)
159 ///
160 /// [Karney (2013)]: https://arxiv.org/pdf/1109.4448.pdf
161 fn geodesic_perimeter_area_unsigned(&self) -> (T, T);
162}
163
164impl GeodesicArea<f64> for Polygon {
165 fn geodesic_perimeter(&self) -> f64 {
166 let (perimeter, _area) = geodesic_area(self, true, false, false);
167 perimeter
168 }
169
170 fn geodesic_area_signed(&self) -> f64 {
171 let (_perimeter, area) = geodesic_area(self, true, false, false);
172 area
173 }
174
175 fn geodesic_area_unsigned(&self) -> f64 {
176 let (_perimeter, area) = geodesic_area(self, false, false, false);
177 area
178 }
179
180 fn geodesic_perimeter_area_signed(&self) -> (f64, f64) {
181 geodesic_area(self, true, false, false)
182 }
183
184 fn geodesic_perimeter_area_unsigned(&self) -> (f64, f64) {
185 geodesic_area(self, false, false, false)
186 }
187}
188
189fn geodesic_area(poly: &Polygon, sign: bool, reverse: bool, exterior_only: bool) -> (f64, f64) {
190 let g = Geodesic::wgs84();
191
192 let (exterior_winding, interior_winding) = if reverse {
193 (Winding::Clockwise, Winding::CounterClockwise)
194 } else {
195 (Winding::CounterClockwise, Winding::Clockwise)
196 };
197
198 // Add the exterior ring
199 let (outer_perimeter, outer_area) = {
200 let mut pa = PolygonArea::new(&g, exterior_winding);
201 poly.exterior().points().for_each(|p| {
202 pa.add_point(p.y(), p.x());
203 });
204 let (perimeter, area, _) = pa.compute(sign);
205 (perimeter, area)
206 };
207
208 // Add the interior rings
209 let (interior_perimeter, mut inner_area) = if exterior_only {
210 (0.0, 0.0)
211 } else {
212 let mut inner_area = 0.;
213 let mut inner_perimeter = 0.;
214 poly.interiors().iter().for_each(|ring| {
215 let mut pa = PolygonArea::new(&g, interior_winding);
216 ring.points().for_each(|p| {
217 pa.add_point(p.y(), p.x());
218 });
219 let (perimeter, area, _) = pa.compute(sign);
220 inner_area += area.abs();
221 inner_perimeter += perimeter;
222 });
223 (inner_perimeter, inner_area)
224 };
225
226 if outer_area < 0.0 && inner_area > 0.0 {
227 inner_area = -inner_area;
228 }
229
230 (
231 outer_perimeter + interior_perimeter,
232 outer_area - inner_area,
233 )
234}
235
236/// Generate a `GeodesicArea` implementation where the result is zero.
237macro_rules! zero_impl {
238 ($type:ident) => {
239 impl GeodesicArea<f64> for $type {
240 fn geodesic_perimeter(&self) -> f64 {
241 0.0
242 }
243
244 fn geodesic_area_signed(&self) -> f64 {
245 0.0
246 }
247
248 fn geodesic_area_unsigned(&self) -> f64 {
249 0.0
250 }
251
252 fn geodesic_perimeter_area_signed(&self) -> (f64, f64) {
253 (0.0, 0.0)
254 }
255
256 fn geodesic_perimeter_area_unsigned(&self) -> (f64, f64) {
257 (0.0, 0.0)
258 }
259 }
260 };
261}
262
263/// Generate a `GeodesicArea` implementation which delegates to the `Polygon`
264/// implementation.
265macro_rules! to_polygon_impl {
266 ($type:ident) => {
267 impl GeodesicArea<f64> for $type {
268 fn geodesic_perimeter(&self) -> f64 {
269 self.to_polygon().geodesic_perimeter()
270 }
271
272 fn geodesic_area_signed(&self) -> f64 {
273 self.to_polygon().geodesic_area_signed()
274 }
275
276 fn geodesic_area_unsigned(&self) -> f64 {
277 self.to_polygon().geodesic_area_unsigned()
278 }
279
280 fn geodesic_perimeter_area_signed(&self) -> (f64, f64) {
281 self.to_polygon().geodesic_perimeter_area_signed()
282 }
283
284 fn geodesic_perimeter_area_unsigned(&self) -> (f64, f64) {
285 self.to_polygon().geodesic_perimeter_area_unsigned()
286 }
287 }
288 };
289}
290
291/// Generate a `GeodesicArea` implementation which calculates the area for each of its
292/// sub-components and sums them up.
293macro_rules! sum_impl {
294 ($type:ident) => {
295 impl GeodesicArea<f64> for $type {
296 fn geodesic_perimeter(&self) -> f64 {
297 self.iter()
298 .fold(0.0, |total, next| total + next.geodesic_perimeter())
299 }
300
301 fn geodesic_area_signed(&self) -> f64 {
302 self.iter()
303 .fold(0.0, |total, next| total + next.geodesic_area_signed())
304 }
305
306 fn geodesic_area_unsigned(&self) -> f64 {
307 self.iter()
308 .fold(0.0, |total, next| total + next.geodesic_area_unsigned())
309 }
310
311 fn geodesic_perimeter_area_signed(&self) -> (f64, f64) {
312 self.iter()
313 .fold((0.0, 0.0), |(total_perimeter, total_area), next| {
314 let (perimeter, area) = next.geodesic_perimeter_area_signed();
315 (total_perimeter + perimeter, total_area + area)
316 })
317 }
318
319 fn geodesic_perimeter_area_unsigned(&self) -> (f64, f64) {
320 self.iter()
321 .fold((0.0, 0.0), |(total_perimeter, total_area), next| {
322 let (perimeter, area) = next.geodesic_perimeter_area_unsigned();
323 (total_perimeter + perimeter, total_area + area)
324 })
325 }
326 }
327 };
328}
329
330zero_impl!(Point);
331zero_impl!(Line);
332zero_impl!(LineString);
333zero_impl!(MultiPoint);
334zero_impl!(MultiLineString);
335to_polygon_impl!(Rect);
336to_polygon_impl!(Triangle);
337sum_impl!(GeometryCollection);
338sum_impl!(MultiPolygon);
339
340impl GeodesicArea<f64> for Geometry<f64> {
341 crate::geometry_delegate_impl! {
342 fn geodesic_perimeter(&self) -> f64;
343 fn geodesic_area_signed(&self) -> f64;
344 fn geodesic_area_unsigned(&self) -> f64;
345 fn geodesic_perimeter_area_signed(&self) -> (f64, f64);
346 fn geodesic_perimeter_area_unsigned(&self) -> (f64, f64);
347 }
348}
349
350#[cfg(test)]
351mod test {
352 use super::*;
353 use crate::algorithm::geodesic_length::GeodesicLength;
354 use crate::polygon;
355
356 #[test]
357 fn test_negative() {
358 let polygon = polygon![
359 (x: 125., y: -15.),
360 (x: 144., y: -15.),
361 (x: 154., y: -27.),
362 (x: 148., y: -39.),
363 (x: 130., y: -33.),
364 (x: 117., y: -37.),
365 (x: 113., y: -22.),
366 (x: 125., y: -15.),
367 ];
368 assert_relative_eq!(
369 -7786102826806.07,
370 polygon.geodesic_area_signed(),
371 epsilon = 0.01
372 );
373
374 let geoid = geographiclib_rs::Geodesic::wgs84();
375 assert_relative_eq!(
376 geoid.area() - 7786102826806.07,
377 polygon.geodesic_area_unsigned(),
378 epsilon = 0.01
379 );
380
381 // Confirm that the exterior ring geodesic_length is the same as the perimeter
382 assert_relative_eq!(
383 polygon.exterior().geodesic_length(),
384 polygon.geodesic_perimeter()
385 );
386 }
387
388 #[test]
389 fn test_positive() {
390 let polygon = polygon![
391 (x: 125., y: -15.),
392 (x: 113., y: -22.),
393 (x: 117., y: -37.),
394 (x: 130., y: -33.),
395 (x: 148., y: -39.),
396 (x: 154., y: -27.),
397 (x: 144., y: -15.),
398 (x: 125., y: -15.),
399 ];
400 assert_relative_eq!(
401 7786102826806.07,
402 polygon.geodesic_area_signed(),
403 epsilon = 0.01
404 );
405 assert_relative_eq!(
406 7786102826806.07,
407 polygon.geodesic_area_unsigned(),
408 epsilon = 0.01
409 );
410
411 // Confirm that the exterior ring geodesic_length is the same as the perimeter
412 assert_relative_eq!(
413 polygon.exterior().geodesic_length(),
414 polygon.geodesic_perimeter()
415 );
416 }
417
418 #[test]
419 fn test_missing_endpoint() {
420 let polygon = polygon![
421 (x: 125., y: -15.),
422 (x: 113., y: -22.),
423 (x: 117., y: -37.),
424 (x: 130., y: -33.),
425 (x: 148., y: -39.),
426 (x: 154., y: -27.),
427 (x: 144., y: -15.),
428 // (x: 125., y: -15.), <-- missing endpoint
429 ];
430 assert_relative_eq!(
431 7786102826806.07,
432 polygon.geodesic_area_signed(),
433 epsilon = 0.01
434 );
435 assert_relative_eq!(
436 7786102826806.07,
437 polygon.geodesic_area_unsigned(),
438 epsilon = 0.01
439 );
440
441 // Confirm that the exterior ring geodesic_length is the same as the perimeter
442 assert_relative_eq!(
443 polygon.exterior().geodesic_length(),
444 polygon.geodesic_perimeter()
445 );
446 }
447
448 #[test]
449 fn test_holes() {
450 let mut poly = polygon![
451 exterior: [
452 (x: 0., y: 0.),
453 (x: 10., y: 0.),
454 (x: 10., y: 10.),
455 (x: 0., y: 10.),
456 (x: 0., y: 0.)
457 ],
458 interiors: [
459 [
460 (x: 1., y: 1.),
461 (x: 1., y: 2.),
462 (x: 2., y: 2.),
463 (x: 2., y: 1.),
464 (x: 1., y: 1.),
465 ],
466 [
467 (x: 5., y: 5.),
468 (x: 5., y: 6.),
469 (x: 6., y: 6.),
470 (x: 6., y: 5.),
471 (x: 5., y: 5.)
472 ],
473 ],
474 ];
475
476 assert_relative_eq!(
477 1203317999173.7063,
478 poly.geodesic_area_signed(),
479 epsilon = 0.01
480 );
481 assert_relative_eq!(
482 1203317999173.7063,
483 poly.geodesic_area_unsigned(),
484 epsilon = 0.01
485 );
486 assert_relative_eq!(5307742.446635911, poly.geodesic_perimeter(), epsilon = 0.01);
487
488 let (perimeter, area) = poly.geodesic_perimeter_area_signed();
489
490 assert_relative_eq!(5307742.446635911, perimeter, epsilon = 0.01);
491 assert_relative_eq!(1203317999173.7063, area, epsilon = 0.01);
492
493 let (perimeter, area) = poly.geodesic_perimeter_area_unsigned();
494
495 assert_relative_eq!(5307742.446635911, perimeter, epsilon = 0.01);
496 assert_relative_eq!(1203317999173.7063, area, epsilon = 0.01);
497
498 // Test with exterior and interior both with CW winding
499 use crate::algorithm::winding_order::Winding;
500 poly.exterior_mut(|exterior| {
501 exterior.make_cw_winding();
502 });
503
504 let (perimeter, area) = poly.geodesic_perimeter_area_signed();
505 assert_relative_eq!(-1203317999173.7063, area, epsilon = 0.01);
506 assert_relative_eq!(5307742.446635911, perimeter, epsilon = 0.01);
507
508 // Test with exterior CW and interior CCW winding
509 poly.interiors_mut(|interiors| {
510 for interior in interiors {
511 interior.make_ccw_winding();
512 }
513 });
514
515 let (perimeter, area) = poly.geodesic_perimeter_area_signed();
516 assert_relative_eq!(-1203317999173.7063, area, epsilon = 0.01);
517 assert_relative_eq!(5307742.446635911, perimeter, epsilon = 0.01);
518
519 // Test with exterior and interior both with CCW winding
520 poly.exterior_mut(|exterior| {
521 exterior.make_ccw_winding();
522 });
523
524 let (perimeter, area) = poly.geodesic_perimeter_area_signed();
525 assert_relative_eq!(1203317999173.7063, area, epsilon = 0.01);
526 assert_relative_eq!(5307742.446635911, perimeter, epsilon = 0.01);
527 }
528
529 #[test]
530 fn test_bad_interior_winding() {
531 let poly = polygon![
532 exterior: [
533 (x: 0., y: 0.),
534 (x: 10., y: 0.),
535 (x: 10., y: 10.),
536 (x: 0., y: 10.),
537 (x: 0., y: 0.)
538 ],
539 interiors: [
540 [
541 (x: 1., y: 1.),
542 (x: 2., y: 1.),
543 (x: 2., y: 2.),
544 (x: 1., y: 2.),
545 (x: 1., y: 1.),
546 ],
547 [
548 (x: 5., y: 5.),
549 (x: 6., y: 5.),
550 (x: 6., y: 6.),
551 (x: 5., y: 6.),
552 (x: 5., y: 5.)
553 ],
554 ],
555 ];
556
557 assert_relative_eq!(1203317999173.7063, poly.geodesic_area_signed());
558 }
559
560 #[test]
561 fn test_diamond() {
562 // a diamond shape
563 let mut diamond = polygon![
564 // exterior oriented counter-clockwise
565 exterior: [
566 (x: 1.0, y: 0.0),
567 (x: 2.0, y: 1.0),
568 (x: 1.0, y: 2.0),
569 (x: 0.0, y: 1.0),
570 (x: 1.0, y: 0.0),
571 ],
572 // interior oriented clockwise
573 interiors: [
574 [
575 (x: 1.0, y: 0.5),
576 (x: 0.5, y: 1.0),
577 (x: 1.0, y: 1.5),
578 (x: 1.5, y: 1.0),
579 (x: 1.0, y: 0.5),
580 ],
581 ],
582 ];
583 assert_relative_eq!(18462065880.09138, diamond.geodesic_area_unsigned());
584 assert_relative_eq!(18462065880.09138, diamond.geodesic_area_signed());
585 assert_relative_eq!(941333.0085011568, diamond.geodesic_perimeter());
586
587 let (perimeter, area) = diamond.geodesic_perimeter_area_signed();
588 assert_relative_eq!(941333.0085011568, perimeter);
589 assert_relative_eq!(18462065880.09138, area);
590
591 let (perimeter, area) = diamond.geodesic_perimeter_area_unsigned();
592 assert_relative_eq!(941333.0085011568, perimeter);
593 assert_relative_eq!(18462065880.09138, area);
594
595 // Test with exterior and interior both with CW winding
596 use crate::algorithm::winding_order::Winding;
597 diamond.exterior_mut(|exterior| {
598 exterior.make_cw_winding();
599 });
600
601 let (perimeter, area) = diamond.geodesic_perimeter_area_signed();
602 assert_relative_eq!(-18462065880.09138, area);
603 assert_relative_eq!(941333.0085011568, perimeter);
604
605 // Test with exterior CW and interior CCW winding
606 diamond.interiors_mut(|interiors| {
607 for interior in interiors {
608 interior.make_ccw_winding();
609 }
610 });
611
612 let (perimeter, area) = diamond.geodesic_perimeter_area_signed();
613 assert_relative_eq!(-18462065880.09138, area);
614 assert_relative_eq!(941333.0085011568, perimeter);
615
616 // Test with exterior and interior both with CCW winding
617 diamond.exterior_mut(|exterior| {
618 exterior.make_ccw_winding();
619 });
620
621 let (perimeter, area) = diamond.geodesic_perimeter_area_signed();
622 assert_relative_eq!(18462065880.09138, area);
623 assert_relative_eq!(941333.0085011568, perimeter);
624 }
625
626 #[test]
627 fn test_very_large_polygon() {
628 // Describe a polygon that covers all of the earth EXCEPT this small square.
629 // The outside of the polygon is in this square, the inside of the polygon is the rest of the earth.
630 let polygon_large: Polygon<f64> = polygon![
631 (x: 0.0, y: 0.0),
632 (x: 0.0, y: 1.0),
633 (x: 1.0, y: 1.0),
634 (x: 1.0, y: 0.0),
635 ];
636
637 let area = polygon_large.geodesic_area_unsigned();
638 assert_eq!(area, 510053312945726.94);
639
640 // A very large polygon that covers nearly all the earth, and then a hole that also covers nearly all the earth as well.
641 // This is a neat polygon because signed and unsigned areas are the same, regardless of the winding order.
642 let polygon_large_with_hole: Polygon<f64> = polygon![
643 exterior: [
644 (x: 0.5, y: 0.5),
645 (x: 0.5, y: 1.0),
646 (x: 1.0, y: 1.0),
647 (x: 1.0, y: 0.5),
648 (x: 0.5, y: 0.5),
649 ],
650 interiors: [
651 [
652 (x: 0.0, y: 0.0),
653 (x: 2.0, y: 0.0),
654 (x: 2.0, y: 2.0),
655 (x: 0.0, y: 2.0),
656 (x: 0.0, y: 0.0),
657 ],
658 ],
659 ];
660
661 let area = polygon_large_with_hole.geodesic_area_signed();
662 assert_relative_eq!(area, 46154562709.8, epsilon = 0.1);
663
664 let area = polygon_large_with_hole.geodesic_area_unsigned();
665 assert_relative_eq!(area, 46154562709.8, epsilon = 0.1);
666 }
667}