YetAnotherCoupler 3.2.0
Loading...
Searching...
No Matches
check_overlap.c
Go to the documentation of this file.
1// Copyright (c) 2024 The YAC Authors
2//
3// SPDX-License-Identifier: BSD-3-Clause
4
5#include <stdlib.h>
6#include <math.h>
7#include "utils_core.h"
8#include "geometry.h"
9#include "grid_cell.h"
10
11static double const tol = 1.0e-12;
12
13// returns the direction of the given edge along the equator
14// switching the points a and b of the edge give the following result:
15// a->b == -(b-a)
16// (exeption: if the edge is a meridian the result in always 0)
17static inline int edge_direction(double * a, double * b) {
18
19 double cross_ab_2 = a[0] * b[1] - a[1] * b[0];
20
21 return (cross_ab_2 > - tol) - (cross_ab_2 < tol);
22}
23
25 double point_coords[3], struct yac_grid_cell cell) {
26
27 double second_point[3];
28
29 // if the point is on the pole
30 if (fabs(fabs(point_coords[2]) - 1.0) < tol) {
31 second_point[0] = 1;
32 second_point[1] = 0;
33 second_point[2] = 0;
34 } else if (point_coords[2] > 0) {
35 second_point[0] = 0;
36 second_point[1] = 0;
37 second_point[2] = -1;
38 } else {
39 second_point[0] = 0;
40 second_point[1] = 0;
41 second_point[2] = 1;
42 }
43
44 int through_cell_corner;
45 int edge_crossings;
46 through_cell_corner = 0;
47 edge_crossings = 0;
48
49 // for all edges of cell
50 for (size_t i = 0; i < cell.num_corners; ++i) {
51
52 double * a = cell.coordinates_xyz[i];
53 double * b = cell.coordinates_xyz[(i+1)%cell.num_corners];
54
55 int ret_value;
56 double p[3], q[3];
57
58 ret_value = yac_intersect_vec(cell.edge_type[i], a, b, YAC_LON_CIRCLE_EDGE,
59 point_coords, second_point, p, q);
60
61 // if both edges do not intersect
62 if ((ret_value == -1) || (ret_value == 0)) continue;
63
64 // if p is not the intersection point of both edges
65 if ((ret_value & ((1 << 0) + (1 << 2))) != (1 << 0) + (1 << 2)) {
66
67 // if q is the intersection point of both edges
68 if ((ret_value & ((1 << 1) + (1 << 3))) == (1 << 1) + (1 << 3))
69 p[0] = q[0], p[1] = q[1], p[2] = q[2];
70 else
71 continue;
72 }
73
74 // if the intersection is the point itself
75 if (points_are_identically(point_coords, p))
76 return 1;
77
78 //----------
79 // in case the point edge goes through a endpoint of a cell edge
80 // there are three special cases that need to be taken care of:
81 // 1: the cell is concave an the point edge only touched an inner corner of the cell -> does not switch in/outside
82 // 2: the point edge only touched an outer corner of the cell -> does not switch in/outside
83 // 3: the point edge entered/left the cell through a corner -> switch in/outside
84 //---------
85 if ((points_are_identically(p, a)) ||
86 (points_are_identically(p, b))) {
87
88 // check the direction of the cell edge (edge can have a positive or
89 // negative direction in terms of longitudes)
90 // each cell endpoint is checked twice (for each adjacent edge)
91 // If the direction for two adjacent edges sharing an endpoint that
92 // is crossed by the point edge, then the related endpoint
93 // intersection can be classified as case 3.
94 // case 1 and 2 result in no change to through_cell_corner
95 // case 3 result in an in/decrease of through_cell_corner by 2
96
97 through_cell_corner += edge_direction(a, b);
98 } else {
99
100 //standard case -> we crossed an edge
101 edge_crossings++;
102 }
103 } // (j = 0; j < cell.num_corners; ++j)
104
105 edge_crossings += through_cell_corner / 2;
106
107 // if we have an odd number of edge crossings -> point was inside cell
108 // or the point was directly on an longitude edge
109 return ( (edge_crossings & 1) || (through_cell_corner & 1) );
110}
111
113 double point_coords[3], struct yac_grid_cell cell) {
114
115 double * a = cell.coordinates_xyz[0];
116 double * b = cell.coordinates_xyz[1];
117 double * c = cell.coordinates_xyz[2];
118
119 double cross_edges[3][3];
120
121 crossproduct_kahan(a, b, &(cross_edges[0][0]));
122 crossproduct_kahan(b, c, &(cross_edges[1][0]));
123 crossproduct_kahan(c, a, &(cross_edges[2][0]));
124
125 double sq_sin_edge[3] = {
126 cross_edges[0][0] * cross_edges[0][0] +
127 cross_edges[0][1] * cross_edges[0][1] +
128 cross_edges[0][2] * cross_edges[0][2],
129 cross_edges[1][0] * cross_edges[1][0] +
130 cross_edges[1][1] * cross_edges[1][1] +
131 cross_edges[1][2] * cross_edges[1][2],
132 cross_edges[2][0] * cross_edges[2][0] +
133 cross_edges[2][1] * cross_edges[2][1] +
134 cross_edges[2][2] * cross_edges[2][2]};
135
136 // if all edges of the triangle have a length of at least yac_angle_tol
137 if ((sq_sin_edge[0] > (yac_angle_tol * yac_angle_tol)) &&
138 (sq_sin_edge[1] > (yac_angle_tol * yac_angle_tol)) &&
139 (sq_sin_edge[2] > (yac_angle_tol * yac_angle_tol))) {
140
141 normalise_vector(&(cross_edges[0][0]));
142 normalise_vector(&(cross_edges[1][0]));
143 normalise_vector(&(cross_edges[2][0]));
144
145 double sin_point_edge[3] = {
146 cross_edges[0][0] * point_coords[0] +
147 cross_edges[0][1] * point_coords[1] +
148 cross_edges[0][2] * point_coords[2],
149 cross_edges[1][0] * point_coords[0] +
150 cross_edges[1][1] * point_coords[1] +
151 cross_edges[1][2] * point_coords[2],
152 cross_edges[2][0] * point_coords[0] +
153 cross_edges[2][1] * point_coords[1] +
154 cross_edges[2][2] * point_coords[2]};
155
156 int cell_corner_order =
157 cross_edges[0][0] * c[0] + cross_edges[0][1] * c[1] + cross_edges[0][2] * c[2] >= 0;
158
159 if (cell_corner_order)
160 return (sin_point_edge[0] > - yac_angle_tol) &&
161 (sin_point_edge[1] > - yac_angle_tol) &&
162 (sin_point_edge[2] > - yac_angle_tol);
163 else
164 return (sin_point_edge[0] < yac_angle_tol) &&
165 (sin_point_edge[1] < yac_angle_tol) &&
166 (sin_point_edge[2] < yac_angle_tol);
167
168 // if the length of the first edge is shorter than yac_angle_tol
169 } else if (sq_sin_edge[0] <= (yac_angle_tol * yac_angle_tol)) {
170
171 // if the length of the second edge is also shorter than yac_angle_tol
172 // -> triangle is a point
173 if (sq_sin_edge[1] <= (yac_angle_tol * yac_angle_tol))
174 return
175 compare_angles(get_vector_angle_2(a, point_coords), SIN_COS_TOL) <= 0;
176 // -> triangle is an edge
177 else {
178
179 normalise_vector(&(cross_edges[1][0]));
180 double sin_point_edge_1 =
181 cross_edges[1][0] * point_coords[0] +
182 cross_edges[1][1] * point_coords[1] +
183 cross_edges[1][2] * point_coords[2];
184
185 struct sin_cos_angle edge_angle =
187 struct sin_cos_angle ap_angle = get_vector_angle_2(a, point_coords);
188 struct sin_cos_angle cp_angle = get_vector_angle_2(c, point_coords);
189 // point is on the plane of the edge and between b and c
190 return (fabs(sin_point_edge_1) <= yac_angle_tol) &&
191 (compare_angles(ap_angle, edge_angle) <= 0) &&
192 (compare_angles(cp_angle, edge_angle) <= 0);
193 }
194
195 // -> triangle is an edge
196 } else {
197
198 normalise_vector(&(cross_edges[0][0]));
199
200 double sin_point_edge_0 =
201 cross_edges[0][0] * point_coords[0] +
202 cross_edges[0][1] * point_coords[1] +
203 cross_edges[0][2] * point_coords[2];
204
205 struct sin_cos_angle edge_angle =
207 struct sin_cos_angle ap_angle = get_vector_angle_2(a, point_coords);
208 struct sin_cos_angle bp_angle = get_vector_angle_2(b, point_coords);
209 // point is on the plane of the edge and between b and c
210 return (fabs(sin_point_edge_0) <= yac_angle_tol) &&
211 (compare_angles(ap_angle, edge_angle) <= 0) &&
212 (compare_angles(bp_angle, edge_angle) <= 0);
213 }
214}
215
216// static int any_point_in_cell_unstruct_triangle(
217 // double * point_coords, size_t num_points, struct yac_grid_cell cell) {
218
219 // double * a = cell.coordinates_xyz[0];
220 // double * b = cell.coordinates_xyz[1];
221 // double * c = cell.coordinates_xyz[2];
222
223 // double cross_edges[3][3];
224
225 // crossproduct_kahan(a, b, &(cross_edges[0][0]));
226 // crossproduct_kahan(b, c, &(cross_edges[1][0]));
227 // crossproduct_kahan(c, a, &(cross_edges[2][0]));
228
229 // double sq_sin_edge[3] = {
230 // cross_edges[0][0] * cross_edges[0][0] +
231 // cross_edges[0][1] * cross_edges[0][1] +
232 // cross_edges[0][2] * cross_edges[0][2],
233 // cross_edges[1][0] * cross_edges[1][0] +
234 // cross_edges[1][1] * cross_edges[1][1] +
235 // cross_edges[1][2] * cross_edges[1][2],
236 // cross_edges[2][0] * cross_edges[2][0] +
237 // cross_edges[2][1] * cross_edges[2][1] +
238 // cross_edges[2][2] * cross_edges[2][2]};
239
240 // normalise_vector(&(cross_edges[0][0]));
241 // normalise_vector(&(cross_edges[1][0]));
242 // normalise_vector(&(cross_edges[2][0]));
243
244 // // if all edges of the triangle have a length of at least yac_angle_tol
245 // if ((sq_sin_edge[0] > (yac_angle_tol * yac_angle_tol)) &&
246 // (sq_sin_edge[1] > (yac_angle_tol * yac_angle_tol)) &&
247 // (sq_sin_edge[2] > (yac_angle_tol * yac_angle_tol))) {
248
249 // int cell_corner_order =
250 // cross_edges[0][0] * c[0] + cross_edges[0][1] * c[1] + cross_edges[0][2] * c[2] >= 0;
251
252 // if (cell_corner_order) {
253 // for (size_t i = 0; i < num_points; ++i) {
254 // double sin_point_edge[3] = {
255 // cross_edges[0][0] * point_coords[0+3*i] +
256 // cross_edges[0][1] * point_coords[1+3*i] +
257 // cross_edges[0][2] * point_coords[2+3*i],
258 // cross_edges[1][0] * point_coords[0+3*i] +
259 // cross_edges[1][1] * point_coords[1+3*i] +
260 // cross_edges[1][2] * point_coords[2+3*i],
261 // cross_edges[2][0] * point_coords[0+3*i] +
262 // cross_edges[2][1] * point_coords[1+3*i] +
263 // cross_edges[2][2] * point_coords[2+3*i]};
264 // if ((sin_point_edge[0] > - yac_angle_tol) &&
265 // (sin_point_edge[1] > - yac_angle_tol) &&
266 // (sin_point_edge[2] > - yac_angle_tol)) return 1;
267 // }
268 // } else {
269 // for (size_t i = 0; i < num_points; ++i) {
270 // double sin_point_edge[3] = {
271 // cross_edges[0][0] * point_coords[0+3*i] +
272 // cross_edges[0][1] * point_coords[1+3*i] +
273 // cross_edges[0][2] * point_coords[2+3*i],
274 // cross_edges[1][0] * point_coords[0+3*i] +
275 // cross_edges[1][1] * point_coords[1+3*i] +
276 // cross_edges[1][2] * point_coords[2+3*i],
277 // cross_edges[2][0] * point_coords[0+3*i] +
278 // cross_edges[2][1] * point_coords[1+3*i] +
279 // cross_edges[2][2] * point_coords[2+3*i]};
280 // if ((sin_point_edge[0] < yac_angle_tol) &&
281 // (sin_point_edge[1] < yac_angle_tol) &&
282 // (sin_point_edge[2] < yac_angle_tol)) return 1;
283 // }
284 // }
285
286 // // if the length of the first edge is shorter than yac_angle_tol
287 // } else if (sq_sin_edge[0] <= (yac_angle_tol * yac_angle_tol)) {
288
289 // // if the length of the second edge is also shorter than yac_angle_tol
290 // // -> triangle is a point
291 // if (sq_sin_edge[1] <= (yac_angle_tol * yac_angle_tol))
292 // return
293 // compare_angles(get_vector_angle_2(a, point_coords), SIN_COS_TOL) <= 0;
294 // // -> triangle is an edge
295 // else {
296 // struct sin_cos_angle edge_angle =
297 // sum_angles_no_check(get_vector_angle_2(a, c), SIN_COS_TOL);
298
299 // for (size_t i = 0; i < num_points; ++i) {
300 // struct sin_cos_angle ap_angle = get_vector_angle_2(a, point_coords+3*i);
301 // struct sin_cos_angle cp_angle = get_vector_angle_2(c, point_coords+3*i);
302 // double sin_point_edge = cross_edges[1][0] * point_coords[0+3*i] +
303 // cross_edges[1][1] * point_coords[1+3*i] +
304 // cross_edges[1][2] * point_coords[2+3*i];
305 // // point is on the plane of the edge and between b and c
306 // if ((fabs(sin_point_edge) <= yac_angle_tol) &&
307 // (compare_angles(ap_angle, edge_angle) <= 0) &&
308 // (compare_angles(cp_angle, edge_angle) <= 0)) return 1;
309 // }
310 // }
311
312 // // -> triangle is an edge
313 // } else {
314 // struct sin_cos_angle edge_angle =
315 // sum_angles_no_check(get_vector_angle_2(a, b), SIN_COS_TOL);
316
317 // for (size_t i = 0; i < num_points; ++i) {
318 // struct sin_cos_angle ap_angle = get_vector_angle_2(a, point_coords+3*i);
319 // struct sin_cos_angle bp_angle = get_vector_angle_2(b, point_coords+3*i);
320 // double sin_point_edge = cross_edges[0][0] * point_coords[0+3*i] +
321 // cross_edges[0][1] * point_coords[1+3*i] +
322 // cross_edges[0][2] * point_coords[2+3*i];
323 // // point is on the plane of the edge and between b and c
324 // if ((fabs(sin_point_edge) <= yac_angle_tol) &&
325 // (compare_angles(ap_angle, edge_angle) <= 0) &&
326 // (compare_angles(bp_angle, edge_angle) <= 0)) return 1;
327 // }
328 // }
329
330 // return 0;
331// }
332
333static inline double dotproduct(double a[], double b[]) {
334 return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
335}
336
338 double point_coords[3], struct yac_grid_cell cell) {
339
340 double * a = cell.coordinates_xyz[0];
341 double * b = cell.coordinates_xyz[1];
342 double * c = cell.coordinates_xyz[2];
343 double * d = cell.coordinates_xyz[3];
344
345 double cross[4][3];
346
347 crossproduct_kahan(a, b, &(cross[0][0]));
348 crossproduct_kahan(b, c, &(cross[1][0]));
349 crossproduct_kahan(c, d, &(cross[2][0]));
350 crossproduct_kahan(d, a, &(cross[3][0]));
351
352 double sq_sin_edge[4] = {
353 cross[0][0]*cross[0][0]+cross[0][1]*cross[0][1]+cross[0][2]*cross[0][2],
354 cross[1][0]*cross[1][0]+cross[1][1]*cross[1][1]+cross[1][2]*cross[1][2],
355 cross[2][0]*cross[2][0]+cross[2][1]*cross[2][1]+cross[2][2]*cross[2][2],
356 cross[3][0]*cross[3][0]+cross[3][1]*cross[3][1]+cross[3][2]*cross[3][2]};
357
358 // if all edges of the quad have a length of at least yac_angle_tol
359 if ((sq_sin_edge[0] > (yac_angle_tol * yac_angle_tol)) &&
360 (sq_sin_edge[1] > (yac_angle_tol * yac_angle_tol)) &&
361 (sq_sin_edge[2] > (yac_angle_tol * yac_angle_tol)) &&
362 (sq_sin_edge[3] > (yac_angle_tol * yac_angle_tol))) {
363
364 normalise_vector(&(cross[0][0]));
365 normalise_vector(&(cross[1][0]));
366 normalise_vector(&(cross[2][0]));
367 normalise_vector(&(cross[3][0]));
368
369 int tmp = dotproduct(c, cross[0]) >= 0.0;
370 tmp += dotproduct(d, cross[0]) >= 0.0;
371 tmp += dotproduct(a, cross[1]) >= 0.0;
372 tmp += dotproduct(d, cross[1]) >= 0.0;
373 tmp += dotproduct(a, cross[2]) >= 0.0;
374 tmp += dotproduct(b, cross[2]) >= 0.0;
375 tmp += dotproduct(b, cross[3]) >= 0.0;
376 tmp += dotproduct(c, cross[3]) >= 0.0;
377 // check whether the quad is convex
378 if (!(tmp & 7)) {
379 if (tmp)
380 return
381 (dotproduct(point_coords, cross[0]) + yac_angle_tol >= 0.0) &&
382 (dotproduct(point_coords, cross[1]) + yac_angle_tol >= 0.0) &&
383 (dotproduct(point_coords, cross[2]) + yac_angle_tol >= 0.0) &&
384 (dotproduct(point_coords, cross[3]) + yac_angle_tol >= 0.0);
385 else
386 return
387 (dotproduct(point_coords, cross[0]) - yac_angle_tol <= -0.0) &&
388 (dotproduct(point_coords, cross[1]) - yac_angle_tol <= -0.0) &&
389 (dotproduct(point_coords, cross[2]) - yac_angle_tol <= -0.0) &&
390 (dotproduct(point_coords, cross[3]) - yac_angle_tol <= -0.0);
391 }
392 }
393
394 return point_in_cell_generic(point_coords, cell);
395}
396
397// static int any_point_in_cell_unstruct_quad(
398 // double * point_coords, size_t num_points, struct yac_grid_cell cell) {
399
400 // double * a = cell.coordinates_xyz[0];
401 // double * b = cell.coordinates_xyz[1];
402 // double * c = cell.coordinates_xyz[2];
403 // double * d = cell.coordinates_xyz[3];
404
405 // double cross[4][3];
406
407 // crossproduct_kahan(a, b, &(cross[0][0]));
408 // crossproduct_kahan(b, c, &(cross[1][0]));
409 // crossproduct_kahan(c, d, &(cross[2][0]));
410 // crossproduct_kahan(d, a, &(cross[3][0]));
411
412 // double sq_sin_edge[4] = {
413 // cross[0][0]*cross[0][0]+cross[0][1]*cross[0][1]+cross[0][2]*cross[0][2],
414 // cross[1][0]*cross[1][0]+cross[1][1]*cross[1][1]+cross[1][2]*cross[1][2],
415 // cross[2][0]*cross[2][0]+cross[2][1]*cross[2][1]+cross[2][2]*cross[2][2],
416 // cross[3][0]*cross[3][0]+cross[3][1]*cross[3][1]+cross[3][2]*cross[3][2]};
417
418 // normalise_vector(&(cross[0][0]));
419 // normalise_vector(&(cross[1][0]));
420 // normalise_vector(&(cross[2][0]));
421 // normalise_vector(&(cross[3][0]));
422
423 // // if all edges of the quad have a length of at least yac_angle_tol
424 // if ((sq_sin_edge[0] > (yac_angle_tol * yac_angle_tol)) &&
425 // (sq_sin_edge[1] > (yac_angle_tol * yac_angle_tol)) &&
426 // (sq_sin_edge[2] > (yac_angle_tol * yac_angle_tol)) &&
427 // (sq_sin_edge[3] > (yac_angle_tol * yac_angle_tol))) {
428
429 // int tmp = dotproduct(c, cross[0]) >= 0.0;
430 // tmp += dotproduct(d, cross[0]) >= 0.0;
431 // tmp += dotproduct(a, cross[1]) >= 0.0;
432 // tmp += dotproduct(d, cross[1]) >= 0.0;
433 // tmp += dotproduct(a, cross[2]) >= 0.0;
434 // tmp += dotproduct(b, cross[2]) >= 0.0;
435 // tmp += dotproduct(b, cross[3]) >= 0.0;
436 // tmp += dotproduct(c, cross[3]) >= 0.0;
437 // // check whether the quad is convex
438 // if (!(tmp & 7)) {
439 // if (tmp) {
440 // for (size_t i = 0; i < num_points; ++i) {
441 // if ((dotproduct(point_coords+3*i, cross[0]) + yac_angle_tol >= 0.0) &&
442 // (dotproduct(point_coords+3*i, cross[1]) + yac_angle_tol >= 0.0) &&
443 // (dotproduct(point_coords+3*i, cross[2]) + yac_angle_tol >= 0.0) &&
444 // (dotproduct(point_coords+3*i, cross[3]) + yac_angle_tol >= 0.0))
445 // return 1;
446 // }
447 // } else {
448 // for (size_t i = 0; i < num_points; ++i) {
449 // if ((dotproduct(point_coords+3*i, cross[0]) - yac_angle_tol <= -0.0) &&
450 // (dotproduct(point_coords+3*i, cross[1]) - yac_angle_tol <= -0.0) &&
451 // (dotproduct(point_coords+3*i, cross[2]) - yac_angle_tol <= -0.0) &&
452 // (dotproduct(point_coords+3*i, cross[3]) - yac_angle_tol <= -0.0))
453 // return 1;
454 // }
455 // }
456 // }
457 // } else {
458 // for (size_t i = 0; i < num_points; ++i)
459 // return point_in_cell_generic(point_coords+3*i, cell);
460 // }
461 // return 0;
462// }
463
464// checks whether a point is within a cell
465int yac_point_in_cell (double point_coords[3], struct yac_grid_cell cell) {
466
467 struct bounding_circle bnd_circle;
468 yac_get_cell_bounding_circle(cell, &bnd_circle);
469
470 return yac_point_in_cell2(point_coords, cell, bnd_circle);
471}
472
473// checks whether a point is within a cell
474int yac_point_in_cell2 (double point_coords[3], struct yac_grid_cell cell,
475 struct bounding_circle bnd_circle) {
476
477 // check whether the point is within the bounding circle of the cell;
478
479 if (!yac_point_in_bounding_circle_vec(point_coords, &bnd_circle))
480 return 1 == 0;
481
482 // check if the point matches any of cell corners
483 for (size_t i = 0; i < cell.num_corners; ++i)
484 if (points_are_identically(point_coords, cell.coordinates_xyz[i]))
485 return 1;
486
487 if ((cell.num_corners == 3) &&
488 (cell.edge_type[0] == YAC_GREAT_CIRCLE_EDGE) &&
489 (cell.edge_type[1] == YAC_GREAT_CIRCLE_EDGE) &&
490 (cell.edge_type[2] == YAC_GREAT_CIRCLE_EDGE))
491 return point_in_cell_unstruct_triangle(point_coords, cell);
492 else if ((cell.num_corners == 4) &&
493 (cell.edge_type[0] == YAC_GREAT_CIRCLE_EDGE) &&
494 (cell.edge_type[1] == YAC_GREAT_CIRCLE_EDGE) &&
495 (cell.edge_type[2] == YAC_GREAT_CIRCLE_EDGE) &&
496 (cell.edge_type[3] == YAC_GREAT_CIRCLE_EDGE))
497 return point_in_cell_unstruct_quad(point_coords, cell);
498 else
499 return point_in_cell_generic(point_coords, cell);
500}
void yac_get_cell_bounding_circle(struct yac_grid_cell cell, struct bounding_circle *bnd_circle)
Definition bnd_circle.c:160
static int edge_direction(double *a, double *b)
int yac_point_in_cell(double point_coords[3], struct yac_grid_cell cell)
static double dotproduct(double a[], double b[])
int yac_point_in_cell2(double point_coords[3], struct yac_grid_cell cell, struct bounding_circle bnd_circle)
static int point_in_cell_unstruct_quad(double point_coords[3], struct yac_grid_cell cell)
static int point_in_cell_unstruct_triangle(double point_coords[3], struct yac_grid_cell cell)
static double const tol
static int point_in_cell_generic(double point_coords[3], struct yac_grid_cell cell)
static struct sin_cos_angle get_vector_angle_2(double const a[3], double const b[3])
Definition geometry.h:415
static int points_are_identically(double const *a, double const *b)
Definition geometry.h:648
static int yac_point_in_bounding_circle_vec(double point_vector[3], struct bounding_circle *bnd_circle)
Definition geometry.h:266
static struct sin_cos_angle sum_angles_no_check(struct sin_cos_angle a, struct sin_cos_angle b)
Definition geometry.h:524
static void crossproduct_kahan(double const a[], double const b[], double cross[])
Definition geometry.h:332
static const struct sin_cos_angle SIN_COS_TOL
Definition geometry.h:45
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])
#define yac_angle_tol
Definition geometry.h:34
static void normalise_vector(double v[])
Definition geometry.h:689
@ YAC_GREAT_CIRCLE_EDGE
great circle
Definition grid_cell.h:13
@ YAC_LON_CIRCLE_EDGE
longitude circle
Definition grid_cell.h:15
size_t num_corners
Definition grid_cell.h:21
enum yac_edge_type * edge_type
Definition grid_cell.h:20
double(* coordinates_xyz)[3]
Definition grid_cell.h:19