YAC 3.12.0
Yet Another Coupler
Loading...
Searching...
No Matches
bnd_circle.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#ifdef HAVE_CONFIG_H
6// Get the definition of the 'restrict' keyword.
7#include "config.h"
8#endif
9
10#include <math.h>
11#include <stdlib.h>
12#include <string.h>
13#include <float.h>
14
15#include "geometry.h"
16#include "utils_core.h"
17#include "yac_types.h"
18#include "ensure_array_size.h"
19
25static inline double get_sin_vector_angle(
26 double a[3], double b[3]) {
27
28 double cross_ab[3];
29 crossproduct_kahan(a, b, cross_ab);
30
31 return sqrt(cross_ab[0]*cross_ab[0] +
32 cross_ab[1]*cross_ab[1] +
33 cross_ab[2]*cross_ab[2]);
34}
35
36// computes the bounding circle of a triangle on the sphere
37// it is assumed that all edges are great circles
39 double a[3], double b[3], double c[3],
40 struct bounding_circle * bnd_circle) {
41
42 double middle_point[3];
43
44 middle_point[0] = a[0] + b[0] + c[0];
45 middle_point[1] = a[1] + b[1] + c[1];
46 middle_point[2] = a[2] + b[2] + c[2];
47
48 normalise_vector(middle_point);
49
50 double cos_angles[3] = {middle_point[0] * a[0] +
51 middle_point[1] * a[1] +
52 middle_point[2] * a[2],
53 middle_point[0] * b[0] +
54 middle_point[1] * b[1] +
55 middle_point[2] * b[2],
56 middle_point[0] * c[0] +
57 middle_point[1] * c[1] +
58 middle_point[2] * c[2]};
59
60 struct sin_cos_angle inc_angle;
61
62 // find the biggest angle
63
64 if (cos_angles[0] < cos_angles[1]) {
65 if (cos_angles[0] < cos_angles[2]) {
66 inc_angle =
67 sin_cos_angle_new(get_sin_vector_angle(middle_point, a), cos_angles[0]);
68 } else {
69 inc_angle =
70 sin_cos_angle_new(get_sin_vector_angle(middle_point, c), cos_angles[2]);
71 }
72 } else {
73 if (cos_angles[1] < cos_angles[2]) {
74 inc_angle =
75 sin_cos_angle_new(get_sin_vector_angle(middle_point, b), cos_angles[1]);
76 } else {
77 inc_angle =
78 sin_cos_angle_new(get_sin_vector_angle(middle_point, c), cos_angles[2]);
79 }
80 }
81
82 bnd_circle->base_vector[0] = middle_point[0];
83 bnd_circle->base_vector[1] = middle_point[1];
84 bnd_circle->base_vector[2] = middle_point[2];
85 bnd_circle->inc_angle = sum_angles_no_check(inc_angle, SIN_COS_TOL);
86 bnd_circle->sq_crd = DBL_MAX;
87}
88
89// computes the bounding circle of a quad on the sphere
90// it is assumed that all edges are great circles
92 double a[3], double b[3], double c[3], double d[3],
93 struct bounding_circle * bnd_circle) {
94
95 double middle_point[3];
96
97 middle_point[0] = a[0] + b[0] + c[0] + d[0];
98 middle_point[1] = a[1] + b[1] + c[1] + d[1];
99 middle_point[2] = a[2] + b[2] + c[2] + d[2];
100
101 normalise_vector(middle_point);
102
103 double cos_angle_a = middle_point[0] * a[0] +
104 middle_point[1] * a[1] +
105 middle_point[2] * a[2];
106 double cos_angle_b = middle_point[0] * b[0] +
107 middle_point[1] * b[1] +
108 middle_point[2] * b[2];
109 double cos_angle_c = middle_point[0] * c[0] +
110 middle_point[1] * c[1] +
111 middle_point[2] * c[2];
112 double cos_angle_d = middle_point[0] * d[0] +
113 middle_point[1] * d[1] +
114 middle_point[2] * d[2];
115
116 struct sin_cos_angle inc_angle;
117
118 // find the biggest angle
119
120 double * temp_x, * temp_y;
121 double temp_cos_angle_x, temp_cos_angle_y;
122
123 if (cos_angle_a < cos_angle_b) {
124 temp_x = a;
125 temp_cos_angle_x = cos_angle_a;
126 } else {
127 temp_x = b;
128 temp_cos_angle_x = cos_angle_b;
129 }
130
131 if (cos_angle_c < cos_angle_d) {
132 temp_y = c;
133 temp_cos_angle_y = cos_angle_c;
134 } else {
135 temp_y = d;
136 temp_cos_angle_y = cos_angle_d;
137 }
138
139 if (temp_cos_angle_x < temp_cos_angle_y) {
140 inc_angle =
142 get_sin_vector_angle(middle_point, temp_x), temp_cos_angle_x);
143 } else {
144 inc_angle =
146 get_sin_vector_angle(middle_point, temp_y), temp_cos_angle_y);
147 }
148
149 bnd_circle->base_vector[0] = middle_point[0];
150 bnd_circle->base_vector[1] = middle_point[1];
151 bnd_circle->base_vector[2] = middle_point[2];
152 bnd_circle->inc_angle = sum_angles_no_check(inc_angle, SIN_COS_TOL);
153 bnd_circle->sq_crd = DBL_MAX;
154}
155
156// computes the bounding circle of a pentagon on the sphere
157// it is assumed that all edges are great circles
159 double a[3], double b[3], double c[3], double d[3], double e[3],
160 struct bounding_circle * bnd_circle) {
161
162 double middle_point[3];
163
164 middle_point[0] = a[0] + b[0] + c[0] + d[0] + e[0];
165 middle_point[1] = a[1] + b[1] + c[1] + d[1] + e[1];
166 middle_point[2] = a[2] + b[2] + c[2] + d[2] + e[2];
167
168 normalise_vector(middle_point);
169
170 double cos_angle_a = middle_point[0] * a[0] +
171 middle_point[1] * a[1] +
172 middle_point[2] * a[2];
173 double cos_angle_b = middle_point[0] * b[0] +
174 middle_point[1] * b[1] +
175 middle_point[2] * b[2];
176 double cos_angle_c = middle_point[0] * c[0] +
177 middle_point[1] * c[1] +
178 middle_point[2] * c[2];
179 double cos_angle_d = middle_point[0] * d[0] +
180 middle_point[1] * d[1] +
181 middle_point[2] * d[2];
182 double cos_angle_e = middle_point[0] * e[0] +
183 middle_point[1] * e[1] +
184 middle_point[2] * e[2];
185
186 struct sin_cos_angle inc_angle;
187
188 // find the biggest angle
189
190 double * temp_x, * temp_y;
191 double temp_cos_angle_x, temp_cos_angle_y;
192
193 if (cos_angle_a < cos_angle_b) {
194 temp_x = a;
195 temp_cos_angle_x = cos_angle_a;
196 } else {
197 temp_x = b;
198 temp_cos_angle_x = cos_angle_b;
199 }
200
201 if (cos_angle_c < cos_angle_d) {
202 temp_y = c;
203 temp_cos_angle_y = cos_angle_c;
204 } else {
205 temp_y = d;
206 temp_cos_angle_y = cos_angle_d;
207 }
208
209 if (cos_angle_e < temp_cos_angle_y) {
210 temp_y = e;
211 temp_cos_angle_y = cos_angle_e;
212 }
213
214 if (temp_cos_angle_x < temp_cos_angle_y) {
215 inc_angle =
217 get_sin_vector_angle(middle_point, temp_x), temp_cos_angle_x);
218 } else {
219 inc_angle =
221 get_sin_vector_angle(middle_point, temp_y), temp_cos_angle_y);
222 }
223
224 bnd_circle->base_vector[0] = middle_point[0];
225 bnd_circle->base_vector[1] = middle_point[1];
226 bnd_circle->base_vector[2] = middle_point[2];
227 bnd_circle->inc_angle = sum_angles_no_check(inc_angle, SIN_COS_TOL);
228 bnd_circle->sq_crd = DBL_MAX;
229}
230
231// computes the bounding circle of a hexagon on the sphere
232// it is assumed that all edges are great circles
234 double a[3], double b[3], double c[3], double d[3], double e[3], double f[3],
235 struct bounding_circle * bnd_circle) {
236
237 double middle_point[3];
238
239 middle_point[0] = a[0] + b[0] + c[0] + d[0] + e[0] + f[0];
240 middle_point[1] = a[1] + b[1] + c[1] + d[1] + e[1] + f[1];
241 middle_point[2] = a[2] + b[2] + c[2] + d[2] + e[2] + f[2];
242
243 normalise_vector(middle_point);
244
245 double cos_angle_a = middle_point[0] * a[0] +
246 middle_point[1] * a[1] +
247 middle_point[2] * a[2];
248 double cos_angle_b = middle_point[0] * b[0] +
249 middle_point[1] * b[1] +
250 middle_point[2] * b[2];
251 double cos_angle_c = middle_point[0] * c[0] +
252 middle_point[1] * c[1] +
253 middle_point[2] * c[2];
254 double cos_angle_d = middle_point[0] * d[0] +
255 middle_point[1] * d[1] +
256 middle_point[2] * d[2];
257 double cos_angle_e = middle_point[0] * e[0] +
258 middle_point[1] * e[1] +
259 middle_point[2] * e[2];
260 double cos_angle_f = middle_point[0] * f[0] +
261 middle_point[1] * f[1] +
262 middle_point[2] * f[2];
263
264 struct sin_cos_angle inc_angle;
265
266 // find the biggest angle
267
268 double * temp_x, * temp_y;
269 double temp_cos_angle_x, temp_cos_angle_y;
270
271 if (cos_angle_a < cos_angle_b) {
272 temp_x = a;
273 temp_cos_angle_x = cos_angle_a;
274 } else {
275 temp_x = b;
276 temp_cos_angle_x = cos_angle_b;
277 }
278
279 if (cos_angle_c < cos_angle_d) {
280 temp_y = c;
281 temp_cos_angle_y = cos_angle_c;
282 } else {
283 temp_y = d;
284 temp_cos_angle_y = cos_angle_d;
285 }
286
287 if (cos_angle_e < cos_angle_f) {
288 if (cos_angle_e < temp_cos_angle_y) {
289 temp_y = e;
290 temp_cos_angle_y = cos_angle_e;
291 }
292 } else {
293 if (cos_angle_f < temp_cos_angle_y) {
294 temp_y = f;
295 temp_cos_angle_y = cos_angle_f;
296 }
297 }
298
299 if (temp_cos_angle_x < temp_cos_angle_y) {
300 inc_angle =
302 get_sin_vector_angle(middle_point, temp_x), temp_cos_angle_x);
303 } else {
304 inc_angle =
306 get_sin_vector_angle(middle_point, temp_y), temp_cos_angle_y);
307 }
308
309 bnd_circle->base_vector[0] = middle_point[0];
310 bnd_circle->base_vector[1] = middle_point[1];
311 bnd_circle->base_vector[2] = middle_point[2];
312 bnd_circle->inc_angle = sum_angles_no_check(inc_angle, SIN_COS_TOL);
313 bnd_circle->sq_crd = DBL_MAX;
314}
315
316// computes the bounding circle of a quad on the sphere
317// it is assumed that the edges are circles of longitude and latitude
319 struct yac_grid_cell quad, struct bounding_circle * bnd_circle) {
320
321 // compute edge middle point of the outwards pointing edge
322 // (the one closest to the equator)
323 int first_lat_edge_index = quad.edge_type[0] != YAC_LAT_CIRCLE_EDGE;
324 int second_lat_edge_index = first_lat_edge_index + 2;
325 int equator_lat_edge_start_index =
326 (fabs(quad.coordinates_xyz[first_lat_edge_index][2]) <
327 fabs(quad.coordinates_xyz[second_lat_edge_index][2]))?
328 first_lat_edge_index:second_lat_edge_index;
329
330 int equator_lat_edge_end_index =
331 (equator_lat_edge_start_index == 3)?0:(equator_lat_edge_start_index+1);
332
333 double edge_middle_point[3] =
334 {quad.coordinates_xyz[equator_lat_edge_start_index][0] +
335 quad.coordinates_xyz[equator_lat_edge_end_index][0],
336 quad.coordinates_xyz[equator_lat_edge_start_index][1] +
337 quad.coordinates_xyz[equator_lat_edge_end_index][1],
338 quad.coordinates_xyz[equator_lat_edge_start_index][2]};
339 double norm =
340 sqrt((1.0 - edge_middle_point[2] * edge_middle_point[2]) /
341 (edge_middle_point[0] * edge_middle_point[0] +
342 edge_middle_point[1] * edge_middle_point[1]));
343 edge_middle_point[0] *= norm;
344 edge_middle_point[1] *= norm;
345
347 quad.coordinates_xyz[0], quad.coordinates_xyz[1],
348 quad.coordinates_xyz[2], quad.coordinates_xyz[3],
349 edge_middle_point, bnd_circle);
350}
351
352// computes the bounding circle of a polygon on the sphere
353// it is assumed that all edges are great circles
355 size_t num_corners, double (* restrict coordinates_xyz)[3],
356 struct bounding_circle * bnd_circle) {
357
358 if (num_corners == 0) {
359
360 bnd_circle->base_vector[0] = 1.0;
361 bnd_circle->base_vector[1] = 0.0;
362 bnd_circle->base_vector[2] = 0.0;
363 bnd_circle->inc_angle = SIN_COS_ZERO;
364 bnd_circle->sq_crd = DBL_MAX;
365 return;
366 }
367
368 double middle_point[3] = {0.0, 0.0, 0.0};
369
370 for (size_t i = 0; i < num_corners; ++i) {
371 middle_point[0] += coordinates_xyz[i][0];
372 middle_point[1] += coordinates_xyz[i][1];
373 middle_point[2] += coordinates_xyz[i][2];
374 }
375
376 normalise_vector(middle_point);
377
378 double min_cos_angle = DBL_MAX;
379 double * coordinate_xyz = NULL;
380
381 // find the biggest angle
382
383 for (size_t i = 0; i < num_corners; ++i) {
384
385 double * restrict temp_coordinate_xyz = coordinates_xyz[i];
386 double temp_cos_angle = middle_point[0] * temp_coordinate_xyz[0] +
387 middle_point[1] * temp_coordinate_xyz[1] +
388 middle_point[2] * temp_coordinate_xyz[2];
389 if (temp_cos_angle < min_cos_angle) {
390 min_cos_angle = temp_cos_angle;
391 coordinate_xyz = temp_coordinate_xyz;
392 }
393 }
394
395 struct sin_cos_angle inc_angle =
397 get_sin_vector_angle(middle_point, coordinate_xyz), min_cos_angle);
398
399 bnd_circle->base_vector[0] = middle_point[0];
400 bnd_circle->base_vector[1] = middle_point[1];
401 bnd_circle->base_vector[2] = middle_point[2];
402 bnd_circle->inc_angle = sum_angles_no_check(inc_angle, SIN_COS_TOL);
403 bnd_circle->sq_crd = DBL_MAX;
404}
405
406// computes a) the angle between the middle point of the edge with the middle
407// point of the polygon
408// b) half of the angle between between the two vertices of the edge
409// returns the sum of both angles
411 double * restrict a, double * restrict b, double * restrict middle_point) {
412
413 double edge_middle_point[3] = {a[0] + b[0], a[1] + b[1], a[2] + b[2]};
414 normalise_vector(edge_middle_point);
415
416 struct sin_cos_angle t1 = get_vector_angle_2(edge_middle_point, a);
417 struct sin_cos_angle t2 = get_vector_angle_2(edge_middle_point, middle_point);
418
419 return sum_angles_no_check(t1, t2);
420}
421
423 struct bounding_circle * bnd_circle) {
424
425 int gc_edge_flag = 1;
426 for (size_t i = 0; i < cell.num_corners; ++i)
427 gc_edge_flag &= (cell.edge_type[i] == YAC_GREAT_CIRCLE_EDGE);
428
429 // check for cells that only have gc edges
430 if (gc_edge_flag) {
431
432 // check for triangle
433 if (cell.num_corners == 3)
434
436 cell.coordinates_xyz[0],
437 cell.coordinates_xyz[1],
438 cell.coordinates_xyz[2], bnd_circle);
439
440 // check for quad
441 else if (cell.num_corners == 4)
442
444 cell.coordinates_xyz[0],
445 cell.coordinates_xyz[1],
446 cell.coordinates_xyz[2],
447 cell.coordinates_xyz[3], bnd_circle);
448
449 // check for pentagon
450 else if (cell.num_corners == 5)
451
453 cell.coordinates_xyz[0],
454 cell.coordinates_xyz[1],
455 cell.coordinates_xyz[2],
456 cell.coordinates_xyz[3],
457 cell.coordinates_xyz[4], bnd_circle);
458
459 // check for hexagon
460 else if (cell.num_corners == 6)
461
463 cell.coordinates_xyz[0],
464 cell.coordinates_xyz[1],
465 cell.coordinates_xyz[2],
466 cell.coordinates_xyz[3],
467 cell.coordinates_xyz[4],
468 cell.coordinates_xyz[5], bnd_circle);
469
470 else
471
473 cell.num_corners, cell.coordinates_xyz, bnd_circle);
474
475 // check for regular quad
476 } else if ((cell.num_corners == 4) &&
477 (((cell.edge_type[0] == YAC_LAT_CIRCLE_EDGE) &&
478 (cell.edge_type[1] == YAC_LON_CIRCLE_EDGE) &&
479 (cell.edge_type[2] == YAC_LAT_CIRCLE_EDGE) &&
480 (cell.edge_type[3] == YAC_LON_CIRCLE_EDGE)) ||
481 ((cell.edge_type[0] == YAC_LON_CIRCLE_EDGE) &&
482 (cell.edge_type[1] == YAC_LAT_CIRCLE_EDGE) &&
483 (cell.edge_type[2] == YAC_LON_CIRCLE_EDGE) &&
484 (cell.edge_type[3] == YAC_LAT_CIRCLE_EDGE)))) {
485
487
488 // general case application to all cell types
489 } else {
490
491 double middle_point[3];
492
493 middle_point[0] = 0.0;
494 middle_point[1] = 0.0;
495 middle_point[2] = 0.0;
496
497 size_t num_corners = cell.num_corners;
498
499 // compute the coordinates in rad and 3d
500 for (size_t i = 0; i < num_corners; ++i) {
501
502 middle_point[0] += cell.coordinates_xyz[i][0];
503 middle_point[1] += cell.coordinates_xyz[i][1];
504 middle_point[2] += cell.coordinates_xyz[i][2];
505 }
506
507 normalise_vector(middle_point);
508
509 // compute the angle required for the bounding circle
510 struct sin_cos_angle max_angle =
512 cell.coordinates_xyz[num_corners-1],
513 cell.coordinates_xyz[0], middle_point);
514
515 for (size_t i = 0; i < num_corners-1; ++i) {
516
517 struct sin_cos_angle curr_angle =
519 cell.coordinates_xyz[i], cell.coordinates_xyz[i+1], middle_point);
520
521 if (compare_angles(max_angle, curr_angle) < 0) max_angle = curr_angle;
522 }
523
524 bnd_circle->base_vector[0] = middle_point[0];
525 bnd_circle->base_vector[1] = middle_point[1];
526 bnd_circle->base_vector[2] = middle_point[2];
527
528 bnd_circle->inc_angle = sum_angles_no_check(max_angle, SIN_COS_TOL);
529 bnd_circle->sq_crd = DBL_MAX;
530 }
531}
532
534 struct bounding_circle * extent_b) {
535
536 struct sin_cos_angle base_vector_angle =
537 get_vector_angle_2(extent_a->base_vector, extent_b->base_vector);
538
539 struct sin_cos_angle tmp_angle, inc_angle_sum;
540 int big_sum =
541 sum_angles(extent_a->inc_angle, extent_b->inc_angle, &tmp_angle);
542
543 if (big_sum ||
544 sum_angles(tmp_angle, SIN_COS_TOL, &inc_angle_sum))
545 return 1;
546
547 return compare_angles(base_vector_angle, inc_angle_sum) <= 0;
548}
static 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:38
static void yac_get_cell_bounding_circle_unstruct_hexa(double a[3], double b[3], double c[3], double d[3], double e[3], double f[3], struct bounding_circle *bnd_circle)
Definition bnd_circle.c:233
static void yac_get_cell_bounding_circle_reg_quad(struct yac_grid_cell quad, struct bounding_circle *bnd_circle)
Definition bnd_circle.c:318
static void yac_get_cell_bounding_circle_unstruct_penta(double a[3], double b[3], double c[3], double d[3], double e[3], struct bounding_circle *bnd_circle)
Definition bnd_circle.c:158
static void yac_get_cell_bounding_circle_unstruct(size_t num_corners, double(*restrict coordinates_xyz)[3], struct bounding_circle *bnd_circle)
Definition bnd_circle.c:354
static void yac_get_cell_bounding_circle_unstruct_quad(double a[3], double b[3], double c[3], double d[3], struct bounding_circle *bnd_circle)
Definition bnd_circle.c:91
void yac_get_cell_bounding_circle(struct yac_grid_cell cell, struct bounding_circle *bnd_circle)
Definition bnd_circle.c:422
int yac_extents_overlap(struct bounding_circle *extent_a, struct bounding_circle *extent_b)
Definition bnd_circle.c:533
static struct sin_cos_angle compute_edge_inc_angle(double *restrict a, double *restrict b, double *restrict middle_point)
Definition bnd_circle.c:410
static double get_sin_vector_angle(double a[3], double b[3])
Definition bnd_circle.c:25
static struct sin_cos_angle get_vector_angle_2(double const a[3], double const b[3])
Definition geometry.h:361
static const struct sin_cos_angle SIN_COS_ZERO
Definition geometry.h:36
static struct sin_cos_angle sum_angles_no_check(struct sin_cos_angle a, struct sin_cos_angle b)
Definition geometry.h:470
static void crossproduct_kahan(double const a[], double const b[], double cross[])
Definition geometry.h:297
static const struct sin_cos_angle SIN_COS_TOL
Definition geometry.h:37
static int compare_angles(struct sin_cos_angle a, struct sin_cos_angle b)
Definition geometry.h:374
static void normalise_vector(double v[])
Definition geometry.h:635
static int sum_angles(struct sin_cos_angle a, struct sin_cos_angle b, struct sin_cos_angle *restrict sum)
Definition geometry.h:480
static struct sin_cos_angle sin_cos_angle_new(double sin, double cos)
Definition geometry.h:351
@ 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:53
double base_vector[3]
Definition geometry.h:51
double sq_crd
Definition geometry.h:56
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