YetAnotherCoupler 3.2.0
Loading...
Searching...
No Matches
geometry.h
Go to the documentation of this file.
1// Copyright (c) 2024 The YAC Authors
2//
3// SPDX-License-Identifier: BSD-3-Clause
4
13#ifndef GEOMETRY_H
14#define GEOMETRY_H
15
16#ifdef HAVE_CONFIG_H
17// We include 'config.h' in this header file (which is a bad practice)
18// because we need the definition of the 'restrict' keyword for the
19// inlined functions.
20#include "config.h"
21#endif
22
23#include <math.h>
24#include <float.h>
25
26#include "basic_grid.h"
27#include "grid_cell.h"
28#include "utils_core.h"
29
30#define YAC_RAD (0.01745329251994329576923690768489) // M_PI / 180
31
32// angle tolerance
33
34#define yac_angle_tol (1e-9)
35#define yac_sq_angle_tol (1e-9 * 1e-9)
36#define yac_cos_angle_tol (0.9999999999999999995)
37#define yac_angle_low_tol (1e-11)
38#define yac_cos_angle_low_tol (1.0)
39
41 double sin, cos;
42};
43
44static const struct sin_cos_angle SIN_COS_ZERO = {0.0, 1.0};
47
48static const struct sin_cos_angle SIN_COS_M_PI_2 = {1.0, 0.0}; /* PI/2 */
49static const struct sin_cos_angle SIN_COS_M_PI = {0.0, -1.0}; /* PI */
50static const struct sin_cos_angle SIN_COS_7_M_PI_4 = {-0.707106781186547524401, +0.707106781186547524401}; /* (7 * PI)/4 */
51
56
59 double base_vector[3];
64 double sq_crd; // sq_crd = (chord)^2
65};
66
78
79// each circle partitions the sphere into an inside and an outside part
80// (the point is the exception; everything except the point is outside)
81struct yac_circle {
83 union {
84 struct {
85 // the norm vector points into the inside direction
86 double norm_vector[3];
87 } gc, lon;
88 struct {
90 double z;
91 } lat;
92 struct {
93 double vec[3];
94 } p;
96};
97
108int yac_point_in_cell (double point_coords[3], struct yac_grid_cell cell);
109
117int yac_point_in_cell2 (double point_coords[3], struct yac_grid_cell cell,
118 struct bounding_circle bnd_circle);
119
127static inline double get_angle (double a_lon, double b_lon) {
128 double diff = a_lon - b_lon;
129#if defined(CDO)
130 return diff - lround(diff / (2.0 * M_PI)) * (2.0 * M_PI);
131#else
132 return diff - round(diff / (2.0 * M_PI)) * (2.0 * M_PI);
133#endif
134}
135
143static inline double get_angle_deg (double a_lon, double b_lon) {
144 double diff = a_lon - b_lon;
145#if defined(CDO)
146 return diff - lround(diff / 360.0) * (360.0);
147#else
148 return diff - round(diff / 360.0) * (360.0);
149#endif
150}
151
175 enum yac_edge_type edge_type_a, double const a[3], double const b[3],
176 enum yac_edge_type edge_type_b, double const c[3], double const d[3],
177 double p[3], double q[3]);
178
189 struct bounding_circle * bnd_circle);
190
200 double a[3], double b[3], double c[3], struct bounding_circle * bnd_circle);
201
211 double a[3], double b[3], double c[3], struct bounding_circle * bnd_circle);
212
223 double a[3], double b[3], double c[3],
224 struct bounding_circle * bnd_circle);
225
236 double a[3], double b[3], double c[3],
237 struct bounding_circle * bnd_circle);
238
245int yac_extents_overlap(struct bounding_circle * extent_a,
246 struct bounding_circle * extent_b);
247
249static inline double sq_len_diff_vec(double const a[3], double const b[3]) {
250 double ab[3] = {a[0]-b[0], a[1]-b[1], a[2]-b[2]};
251 return ab[0] * ab[0] + ab[1] * ab[1] + ab[2] * ab[2];
252}
253
254static inline double compute_sq_crd(struct sin_cos_angle angle) {
255
256 double temp = 1.0 - angle.cos;
257 return temp * temp + angle.sin * angle.sin;
258}
259
267 double point_vector[3], struct bounding_circle * bnd_circle) {
268
269 if (bnd_circle->sq_crd == DBL_MAX)
270 bnd_circle->sq_crd = compute_sq_crd(bnd_circle->inc_angle);
271
272 return
273 bnd_circle->sq_crd >=
274 sq_len_diff_vec(bnd_circle->base_vector, point_vector);
275}
276
287static inline void LLtoXYZ(double lon, double lat, double p_out[]) {
288
289 while (lon < -M_PI) lon += 2.0 * M_PI;
290 while (lon >= M_PI) lon -= 2.0 * M_PI;
291
292 double cos_lat = cos(lat);
293 p_out[0] = cos_lat * cos(lon);
294 p_out[1] = cos_lat * sin(lon);
295 p_out[2] = sin(lat);
296}
297
304static inline void LLtoXYZ_deg(double lon, double lat, double p_out[]) {
305 LLtoXYZ(lon*YAC_RAD, lat*YAC_RAD, p_out);
306}
307
318static inline void XYZtoLL (double const p_in[], double * lon, double * lat) {
319
320 *lon = atan2(p_in[1] , p_in[0]);
321 *lat = M_PI_2 - acos(p_in[2]);
322}
323
324// computation of the determinant of a 2x2 matrix using Kahan summation
325static inline double det2_kahan(double a, double b, double c, double d) {
326 double bc = b*c;
327 double e_bc = fma(b,c,-bc); // the rounding error of the multiplication
328 double det = fma(a,d,-bc);
329 return det + e_bc;
330}
331
332static inline void crossproduct_kahan (
333 double const a[], double const b[], double cross[]) {
334
335 /* crossproduct in Cartesian coordinates */
336
337 cross[0] = det2_kahan(a[1], a[2], b[1], b[2]);
338 cross[1] = det2_kahan(a[2], a[0], b[2], b[0]);
339 cross[2] = det2_kahan(a[0], a[1], b[0], b[1]);
340
341}
342
347static inline void crossproduct_d (
348 const double a[], const double b[], double cross[]) {
349
350/* crossproduct in Cartesian coordinates */
351
352 cross[0] = a[1] * b[2] - a[2] * b[1];
353 cross[1] = a[2] * b[0] - a[0] * b[2];
354 cross[2] = a[0] * b[1] - a[1] * b[0];
355}
356
364static inline double get_vector_angle(double const a[3], double const b[3]) {
365
366 double dot_product = a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
367
368 // the acos most accurate in the range [-M_SQRT1_2;M_SQRT1_2]
369 if (fabs(dot_product) > M_SQRT1_2) {
370
371 double cross_ab[3];
372
373 crossproduct_kahan(a, b, cross_ab);
374
375 double asin_tmp = asin(sqrt(cross_ab[0]*cross_ab[0]+
376 cross_ab[1]*cross_ab[1]+
377 cross_ab[2]*cross_ab[2]));
378
379 if (dot_product > 0.0) // if the angle is smaller than (PI / 2)
380 return MAX(asin_tmp,0.0);
381 else
382 return MIN(M_PI - asin_tmp, M_PI);
383
384 } else {
385
386 return acos(dot_product);
387 }
388
389 /*
390 // this solution is simpler, but has a much worse performance
391 double cross[3], dot, cross_abs;
392
393 crossproduct_kahan(a, b, cross);
394 dot = a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
395 cross_abs = sqrt(cross[0]*cross[0] + cross[1]*cross[1] + cross[2]*cross[2]);
396
397 return fabs(atan2(cross_abs, dot));
398 */
399}
400
401static inline double clamp_abs_one(double val) {
402 return val > -1.0 ? (val < 1.0 ? val : 1.0) : -1.0;
403}
404
405static inline struct sin_cos_angle sin_cos_angle_new(double sin, double cos) {
406
407 struct sin_cos_angle angle;
408
409 angle.sin = clamp_abs_one(sin);
410 angle.cos = clamp_abs_one(cos);
411
412 return angle;
413}
414
416 double const a[3], double const b[3]) {
417
418 double cross_ab[3];
419 crossproduct_kahan(a, b, cross_ab);
420
421 return sin_cos_angle_new(sqrt(cross_ab[0]*cross_ab[0] +
422 cross_ab[1]*cross_ab[1] +
423 cross_ab[2]*cross_ab[2]),
424 a[0]*b[0] + a[1]*b[1] + a[2]*b[2]);
425}
426
427// works for angles in the range [0;2*PI[
428static inline int compare_angles(
429 struct sin_cos_angle a, struct sin_cos_angle b) {
430
431 // there are 5 sections:
432 // 0: 0 <= angle < PI/4
433 // 1: PI/4 <= angle < 3*PI/4
434 // 2: 3*PI/4 <= angle < 5*PI/4
435 // 3: 5*PI/4 <= angle < 7*PI/4
436 // 4: 7*PI/4 <= angle < 2*PI
437 int t_a = fabs(a.cos) <= M_SQRT1_2;
438 int t_b = fabs(b.cos) <= M_SQRT1_2;
439 int a_section = t_a | ((((a.sin < 0.0) & t_a) |
440 ((a.cos < 0.0) & (fabs(a.sin) < M_SQRT1_2))) << 1);
441 int b_section = t_b | ((((b.sin < 0.0) & t_b) |
442 ((b.cos < 0.0) & (fabs(b.sin) < M_SQRT1_2))) << 1);
443 // if current section is 0, then it could be actually section 4
444 if (!a_section) a_section = (a.sin < 0.0) << 2;
445 if (!b_section) b_section = (b.sin < 0.0) << 2;
446
447 if (a_section != b_section)
448 return (a_section > b_section) - (a_section < b_section);
449
450 int ret;
451
452 switch (a_section) {
453 case(0):
454 case(4):
455 default:
456 ret = (b.sin < a.sin + yac_angle_low_tol) -
457 (a.sin < b.sin + yac_angle_low_tol);
458 if (ret) return ret;
459 else {
460 ret = (a.cos < b.cos + yac_angle_low_tol) -
461 (b.cos < a.cos + yac_angle_low_tol);
462 if (a.sin >= 0.0) return ret;
463 else return -ret;
464 }
465 case(1):
466 ret = (a.cos < b.cos + yac_angle_low_tol) -
467 (b.cos < a.cos + yac_angle_low_tol);
468 if (ret) return ret;
469 else {
470 ret = (b.sin < a.sin + yac_angle_low_tol) -
471 (a.sin < b.sin + yac_angle_low_tol);
472 if (a.cos >= 0.0) return ret;
473 else return -ret;
474 }
475 case(2):
476 ret = (a.sin < b.sin + yac_angle_low_tol) -
477 (b.sin < a.sin + yac_angle_low_tol);
478 if (ret) return ret;
479 else {
480 ret = (a.cos < b.cos + yac_angle_low_tol) -
481 (b.cos < a.cos + yac_angle_low_tol);
482 if (a.sin >= 0.0) return ret;
483 else return -ret;
484 }
485 case(3):
486 ret = (b.cos < a.cos + yac_angle_low_tol) -
487 (a.cos < b.cos + yac_angle_low_tol);
488 if (ret) return ret;
489 else {
490 ret = (a.sin < b.sin + yac_angle_low_tol) -
491 (b.sin < a.sin + yac_angle_low_tol);
492 if (a.cos <= 0.0) return ret;
493 else return -ret;
494 }
495 }
496}
497
505static inline double sin_cos_angle_to_dble(struct sin_cos_angle angle) {
506
507 int t = fabs(angle.cos) <= M_SQRT1_2;
508 int section = t | ((((angle.sin < 0.0) & t) |
509 ((angle.cos < 0.0) & (fabs(angle.sin) < M_SQRT1_2))) << 1);
510 if (!section) section = (angle.sin < 0.0) << 2;
511
512 switch (section) {
513 default:
514 case(0): return 0.0 * M_SQRT1_2 + angle.sin;
515 case(1): return 2.0 * M_SQRT1_2 - angle.cos;
516 case(2): return 4.0 * M_SQRT1_2 - angle.sin;
517 case(3): return 6.0 * M_SQRT1_2 + angle.cos;
518 case(4): return 8.0 * M_SQRT1_2 + angle.sin;
519 }
520}
521
525 struct sin_cos_angle a, struct sin_cos_angle b) {
526
527 return sin_cos_angle_new(a.sin * b.cos + a.cos * b.sin,
528 a.cos * b.cos - a.sin * b.sin);
529}
530
534static inline int sum_angles(
535 struct sin_cos_angle a, struct sin_cos_angle b,
536 struct sin_cos_angle * restrict sum) {
537
538 struct sin_cos_angle sum_ = sum_angles_no_check(a, b);
539
540 *sum = sum_;
541
542 // if a or b is smaller than the result
543 return (compare_angles(sum_, a) < 0) || (compare_angles(sum_, b) < 0);
544}
545
549 struct sin_cos_angle a, struct sin_cos_angle b) {
550
551 return sin_cos_angle_new(a.sin * b.cos - a.cos * b.sin,
552 a.cos * b.cos + a.sin * b.sin);
553}
554
558static inline int sub_angles(
559 struct sin_cos_angle a, struct sin_cos_angle b,
560 struct sin_cos_angle * restrict sub) {
561
562 int compare_result = compare_angles(a, b);
563
564 // return sin=0.0 and cos=1.0 if the angles are equal to each other,
565 // i.e. compare_result == 0
566 *sub = compare_result ? sub_angles_no_check(a, b) : SIN_COS_ZERO;
567
568 return compare_result < 0;
569}
570
572static inline double compute_angle(struct sin_cos_angle angle) {
573
574 // the acos and asin are most accurate in the range [-M_SQRT1_2;M_SQRT1_2]
575 if (angle.cos > M_SQRT1_2) {
576
577 double angle_ = asin(angle.sin);
578
579 if (angle_ < 0.0) angle_ += 2.0 * M_PI;
580 return angle_;
581
582 } else if (angle.cos > - M_SQRT1_2) {
583
584 double angle_ = acos(angle.cos);
585
586 if (angle.sin > 0.0) return angle_;
587 else return 2.0 * M_PI - angle_;
588
589 } else {
590 return M_PI - asin(angle.sin);
591 }
592}
593
600static inline struct sin_cos_angle half_angle(struct sin_cos_angle angle) {
601
602 double x = (1.0 + fabs(angle.cos));
603
604 double scale = 1.0 / sqrt(x * x + angle.sin * angle.sin);
605
606 // first or fourth quadrant
607 if (angle.cos >= 0) {
608 scale = copysign(scale, angle.sin);
609 return sin_cos_angle_new(angle.sin * scale, x * scale);
610
611 // second and third quadrant
612 } else return sin_cos_angle_new(x * scale, angle.sin * scale);
613}
614
617static inline struct sin_cos_angle quarter_angle(struct sin_cos_angle angle) {
618
619 double tan = fabs(angle.sin / (1.0 + fabs(angle.cos)));
620 double one_plus_sq_tan = 1.0 + tan * tan;
621 double sqrt_one_plus_sq_tan = sqrt(one_plus_sq_tan);
622
623 double a = 1.0;
624 double b = tan;
625
626 // second and third quadrant
627 if (angle.cos < 0.0) {
628 a = tan;
629 b = 1.0;
630 }
631
632 double scale = M_SQRT1_2 / sqrt(one_plus_sq_tan + a * sqrt_one_plus_sq_tan);
633 double x = b * scale;
634 double y = (a + sqrt_one_plus_sq_tan) * scale;
635
636 // first and second quadrant
637 if (angle.sin >= 0.0) return sin_cos_angle_new(x, y);
638 // third and fourth quadrant
639 else return sin_cos_angle_new(y, x);
640}
641
648static inline int points_are_identically(double const * a, double const * b) {
649
650 // for small angles the Euclidean distance is nearly identical to
651 // great circle distance between vectors a and b
652 return sq_len_diff_vec(a, b) <= yac_sq_angle_tol;
653}
654
662static inline int compare_coords(double const * a, double const * b) {
663
664 double dot_product = a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
665
666 // if the angle is smaller than ~0.81 degree
667 // (in this range the acos is still rather accurate)
668 if (dot_product > 0.9999) { // (acos(0.9999) = ~0.81 degree)
669
670 // both points are close to each other -> use cross product for higher
671 // accuracy
672
673 double cross_ab[3];
674 crossproduct_kahan(a, b, cross_ab);
675
676 // for very small angles: asin(alpha) = ~alpha (alpha in rad)
677 if (sqrt(cross_ab[0]*cross_ab[0] +
678 cross_ab[1]*cross_ab[1] +
679 cross_ab[2]*cross_ab[2]) < yac_angle_tol) return 0;
680 }
681
682 int ret;
683 if ((ret = (a[0] > b[0]) - (a[0] < b[0]))) return ret;
684 if ((ret = (a[1] > b[1]) - (a[1] < b[1]))) return ret;
685 return (a[2] > b[2]) - (a[2] < b[2]);
686}
687
689static inline void normalise_vector(double v[]) {
690
691 double norm = 1.0 / sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
692
693 v[0] *= norm;
694 v[1] *= norm;
695 v[2] *= norm;
696}
697
704static inline void rotate_vector2(
705 double axis[], struct sin_cos_angle angle, double v_in[], double v_out[]) {
706
707 // using Rodrigues' rotation formula
708 // v_out = v_in * cos(angle) +
709 // (axis x v_in) * sin(angle) +
710 // axis * (axis * v_in) * (1 - cos(angle))
711
712 double cross_axis_v_in[3];
713 crossproduct_d(axis, v_in, cross_axis_v_in);
714
715 double dot_axis_v_in = axis[0]*v_in[0] + axis[1]*v_in[1] + axis[2]*v_in[2];
716 double temp = dot_axis_v_in * (1.0 - angle.cos);
717
718 v_out[0] =
719 v_in[0] * angle.cos + cross_axis_v_in[0] * angle.sin + axis[0] * temp;
720 v_out[1] =
721 v_in[1] * angle.cos + cross_axis_v_in[1] * angle.sin + axis[1] * temp;
722 v_out[2] =
723 v_in[2] * angle.cos + cross_axis_v_in[2] * angle.sin + axis[2] * temp;
724}
725
732static inline void rotate_vector(
733 double axis[], double angle, double v_in[], double v_out[]) {
734
735 double sin_angle = sin(angle);
736 double cos_angle = cos(angle);
738 axis, sin_cos_angle_new(sin_angle, cos_angle), v_in, v_out);
739}
740
755 struct yac_grid_cell cell, size_t start_corner, struct yac_grid_cell * triangles);
756
770 size_t const * corner_indices, size_t num_corners, size_t start_corner,
771 size_t (*triangle_indices)[3]);
772
784 double * vertices, size_t num_vertices, double triangle[][3],
785 size_t num_tests);
786
788 struct yac_circle a, struct yac_circle b, double p[3], double q[3]);
789
791 double p[3], double const a[3], double const b[3],
792 enum yac_circle_type circle_type);
793
794#endif // GEOMETRY_H
795
static struct sin_cos_angle get_vector_angle_2(double const a[3], double const b[3])
Definition geometry.h:415
static int compare_coords(double const *a, double const *b)
Definition geometry.h:662
void yac_get_cell_circumscribe_circle_reg_quad(double a[3], double b[3], double c[3], struct bounding_circle *bnd_circle)
Definition bnd_circle.c:27
static double compute_sq_crd(struct sin_cos_angle angle)
Definition geometry.h:254
static const struct sin_cos_angle SIN_COS_LOW_TOL
Definition geometry.h:46
static double clamp_abs_one(double val)
Definition geometry.h:401
static const struct sin_cos_angle SIN_COS_7_M_PI_4
Definition geometry.h:50
void yac_get_cell_circumscribe_circle_unstruct_triangle(double a[3], double b[3], double c[3], struct bounding_circle *bnd_circle)
Definition bnd_circle.c:45
static const struct sin_cos_angle SIN_COS_ZERO
Definition geometry.h:44
int yac_circle_intersect(struct yac_circle a, struct yac_circle b, double p[3], double q[3])
static double det2_kahan(double a, double b, double c, double d)
Definition geometry.h:325
void yac_triangulate_cell_indices(size_t const *corner_indices, size_t num_corners, size_t start_corner, size_t(*triangle_indices)[3])
Definition grid_cell.c:144
static int points_are_identically(double const *a, double const *b)
Definition geometry.h:648
static const struct sin_cos_angle SIN_COS_M_PI
Definition geometry.h:49
void yac_get_cell_bounding_circle_unstruct_triangle(double a[3], double b[3], double c[3], struct bounding_circle *bnd_circle)
Definition bnd_circle.c:93
void yac_compute_bnd_triangle(double *vertices, size_t num_vertices, double triangle[][3], size_t num_tests)
void yac_get_cell_bounding_circle_reg_quad(double a[3], double b[3], double c[3], struct bounding_circle *bnd_circle)
Definition bnd_circle.c:36
static int yac_point_in_bounding_circle_vec(double point_vector[3], struct bounding_circle *bnd_circle)
Definition geometry.h:266
static double sq_len_diff_vec(double const a[3], double const b[3])
computes square of the lenght of the vector ab
Definition geometry.h:249
static void crossproduct_d(const double a[], const double b[], double cross[])
Definition geometry.h:347
static struct sin_cos_angle sum_angles_no_check(struct sin_cos_angle a, struct sin_cos_angle b)
Definition geometry.h:524
static double sin_cos_angle_to_dble(struct sin_cos_angle angle)
Definition geometry.h:505
static void crossproduct_kahan(double const a[], double const b[], double cross[])
Definition geometry.h:332
static void rotate_vector(double axis[], double angle, double v_in[], double v_out[])
Definition geometry.h:732
static struct sin_cos_angle quarter_angle(struct sin_cos_angle angle)
Definition geometry.h:617
static const struct sin_cos_angle SIN_COS_TOL
Definition geometry.h:45
static const struct sin_cos_angle SIN_COS_M_PI_2
Definition geometry.h:48
#define yac_sq_angle_tol
Definition geometry.h:35
#define yac_cos_angle_low_tol
Definition geometry.h:38
int yac_point_in_cell(double point_coords[3], struct yac_grid_cell cell)
static void rotate_vector2(double axis[], struct sin_cos_angle angle, double v_in[], double v_out[])
Definition geometry.h:704
#define yac_cos_angle_tol
Definition geometry.h:36
void yac_triangulate_cell(struct yac_grid_cell cell, size_t start_corner, struct yac_grid_cell *triangles)
Definition grid_cell.c:66
#define YAC_RAD
Definition geometry.h:30
static void LLtoXYZ_deg(double lon, double lat, double p_out[])
Definition geometry.h:304
int yac_point_in_cell2(double point_coords[3], struct yac_grid_cell cell, struct bounding_circle bnd_circle)
static int compare_angles(struct sin_cos_angle a, struct sin_cos_angle b)
Definition geometry.h:428
int yac_intersect_vec(enum yac_edge_type edge_type_a, double const a[3], double const b[3], enum yac_edge_type edge_type_b, double const c[3], double const d[3], double p[3], double q[3])
static int sub_angles(struct sin_cos_angle a, struct sin_cos_angle b, struct sin_cos_angle *restrict sub)
Definition geometry.h:558
int yac_point_on_edge(double p[3], double const a[3], double const b[3], enum yac_circle_type circle_type)
static struct sin_cos_angle sub_angles_no_check(struct sin_cos_angle a, struct sin_cos_angle b)
Definition geometry.h:548
void yac_get_cell_bounding_circle(struct yac_grid_cell cell, struct bounding_circle *bnd_circle)
Definition bnd_circle.c:160
int yac_extents_overlap(struct bounding_circle *extent_a, struct bounding_circle *extent_b)
Definition bnd_circle.c:204
static double compute_angle(struct sin_cos_angle angle)
return angles in the range of [0;2PI[
Definition geometry.h:572
#define yac_angle_low_tol
Definition geometry.h:37
#define yac_angle_tol
Definition geometry.h:34
static void normalise_vector(double v[])
Definition geometry.h:689
static double get_angle(double a_lon, double b_lon)
Definition geometry.h:127
static void LLtoXYZ(double lon, double lat, double p_out[])
Definition geometry.h:287
static double get_vector_angle(double const a[3], double const b[3])
Definition geometry.h:364
static int sum_angles(struct sin_cos_angle a, struct sin_cos_angle b, struct sin_cos_angle *restrict sum)
Definition geometry.h:534
static double get_angle_deg(double a_lon, double b_lon)
Definition geometry.h:143
yac_circle_type
Definition geometry.h:72
@ POINT
Definition geometry.h:76
@ LAT_CIRCLE
Definition geometry.h:74
@ GREAT_CIRCLE
Definition geometry.h:73
@ LON_CIRCLE
Definition geometry.h:75
static struct sin_cos_angle sin_cos_angle_new(double sin, double cos)
Definition geometry.h:405
static struct sin_cos_angle half_angle(struct sin_cos_angle angle)
Definition geometry.h:600
static void XYZtoLL(double const p_in[], double *lon, double *lat)
Definition geometry.h:318
yac_edge_type
Definition grid_cell.h:12
@ YAC_GREAT_CIRCLE_EDGE
great circle
Definition grid_cell.h:13
@ YAC_LAT_CIRCLE_EDGE
latitude circle
Definition grid_cell.h:14
@ YAC_LON_CIRCLE_EDGE
longitude circle
Definition grid_cell.h:15
struct sin_cos_angle inc_angle
angle between the middle point and the boundary of the spherical cap
Definition geometry.h:61
double base_vector[3]
Definition geometry.h:59
double sq_crd
Definition geometry.h:64
double sin
Definition geometry.h:41
double cos
Definition geometry.h:41
struct yac_circle::@3::@4 gc
union yac_circle::@3 data
double z
Definition geometry.h:90
enum yac_circle_type type
Definition geometry.h:82
struct yac_circle::@3::@6 p
int north_is_out
Definition geometry.h:89
struct yac_circle::@3::@5 lat
struct yac_circle::@3::@4 lon
double norm_vector[3]
Definition geometry.h:86
double vec[3]
Definition geometry.h:93
#define MIN(a, b)
#define MAX(a, b)