26 double a[3],
double b[3]) {
31 return sqrt(cross_ab[0]*cross_ab[0] +
32 cross_ab[1]*cross_ab[1] +
33 cross_ab[2]*cross_ab[2]);
39 double a[3],
double b[3],
double c[3],
42 double middle_point[3];
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];
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]};
64 if (cos_angles[0] < cos_angles[1]) {
65 if (cos_angles[0] < cos_angles[2]) {
73 if (cos_angles[1] < cos_angles[2]) {
86 bnd_circle->
sq_crd = DBL_MAX;
92 double a[3],
double b[3],
double c[3],
double d[3],
95 double middle_point[3];
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];
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];
120 double * temp_x, * temp_y;
121 double temp_cos_angle_x, temp_cos_angle_y;
123 if (cos_angle_a < cos_angle_b) {
125 temp_cos_angle_x = cos_angle_a;
128 temp_cos_angle_x = cos_angle_b;
131 if (cos_angle_c < cos_angle_d) {
133 temp_cos_angle_y = cos_angle_c;
136 temp_cos_angle_y = cos_angle_d;
139 if (temp_cos_angle_x < temp_cos_angle_y) {
153 bnd_circle->
sq_crd = DBL_MAX;
159 double a[3],
double b[3],
double c[3],
double d[3],
double e[3],
162 double middle_point[3];
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];
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];
190 double * temp_x, * temp_y;
191 double temp_cos_angle_x, temp_cos_angle_y;
193 if (cos_angle_a < cos_angle_b) {
195 temp_cos_angle_x = cos_angle_a;
198 temp_cos_angle_x = cos_angle_b;
201 if (cos_angle_c < cos_angle_d) {
203 temp_cos_angle_y = cos_angle_c;
206 temp_cos_angle_y = cos_angle_d;
209 if (cos_angle_e < temp_cos_angle_y) {
211 temp_cos_angle_y = cos_angle_e;
214 if (temp_cos_angle_x < temp_cos_angle_y) {
228 bnd_circle->
sq_crd = DBL_MAX;
234 double a[3],
double b[3],
double c[3],
double d[3],
double e[3],
double f[3],
237 double middle_point[3];
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];
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];
268 double * temp_x, * temp_y;
269 double temp_cos_angle_x, temp_cos_angle_y;
271 if (cos_angle_a < cos_angle_b) {
273 temp_cos_angle_x = cos_angle_a;
276 temp_cos_angle_x = cos_angle_b;
279 if (cos_angle_c < cos_angle_d) {
281 temp_cos_angle_y = cos_angle_c;
284 temp_cos_angle_y = cos_angle_d;
287 if (cos_angle_e < cos_angle_f) {
288 if (cos_angle_e < temp_cos_angle_y) {
290 temp_cos_angle_y = cos_angle_e;
293 if (cos_angle_f < temp_cos_angle_y) {
295 temp_cos_angle_y = cos_angle_f;
299 if (temp_cos_angle_x < temp_cos_angle_y) {
313 bnd_circle->
sq_crd = DBL_MAX;
324 int second_lat_edge_index = first_lat_edge_index + 2;
325 int equator_lat_edge_start_index =
328 first_lat_edge_index:second_lat_edge_index;
330 int equator_lat_edge_end_index =
331 (equator_lat_edge_start_index == 3)?0:(equator_lat_edge_start_index+1);
333 double edge_middle_point[3] =
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;
349 edge_middle_point, bnd_circle);
355 size_t num_corners,
double (* restrict coordinates_xyz)[3],
358 if (num_corners == 0) {
364 bnd_circle->
sq_crd = DBL_MAX;
368 double middle_point[3] = {0.0, 0.0, 0.0};
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];
378 double min_cos_angle = DBL_MAX;
379 double * coordinate_xyz = NULL;
383 for (
size_t i = 0; i < num_corners; ++i) {
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;
403 bnd_circle->
sq_crd = DBL_MAX;
411 double * restrict a,
double * restrict b,
double * restrict middle_point) {
413 double edge_middle_point[3] = {a[0] + b[0], a[1] + b[1], a[2] + b[2]};
425 int gc_edge_flag = 1;
491 double middle_point[3];
493 middle_point[0] = 0.0;
494 middle_point[1] = 0.0;
495 middle_point[2] = 0.0;
500 for (
size_t i = 0; i < num_corners; ++i) {
515 for (
size_t i = 0; i < num_corners-1; ++i) {
521 if (
compare_angles(max_angle, curr_angle) < 0) max_angle = curr_angle;
529 bnd_circle->
sq_crd = DBL_MAX;
static void yac_get_cell_bounding_circle_unstruct_triangle(double a[3], double b[3], double c[3], struct bounding_circle *bnd_circle)
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)
static void yac_get_cell_bounding_circle_reg_quad(struct yac_grid_cell quad, struct bounding_circle *bnd_circle)
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)
static void yac_get_cell_bounding_circle_unstruct(size_t num_corners, double(*restrict coordinates_xyz)[3], struct bounding_circle *bnd_circle)
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)
void yac_get_cell_bounding_circle(struct yac_grid_cell cell, struct bounding_circle *bnd_circle)
int yac_extents_overlap(struct bounding_circle *extent_a, struct bounding_circle *extent_b)
static struct sin_cos_angle compute_edge_inc_angle(double *restrict a, double *restrict b, double *restrict middle_point)
static double get_sin_vector_angle(double a[3], double b[3])
static struct sin_cos_angle get_vector_angle_2(double const a[3], double const b[3])
static const struct sin_cos_angle SIN_COS_ZERO
static struct sin_cos_angle sum_angles_no_check(struct sin_cos_angle a, struct sin_cos_angle b)
static void crossproduct_kahan(double const a[], double const b[], double cross[])
static const struct sin_cos_angle SIN_COS_TOL
static int compare_angles(struct sin_cos_angle a, struct sin_cos_angle b)
static void normalise_vector(double v[])
static int sum_angles(struct sin_cos_angle a, struct sin_cos_angle b, struct sin_cos_angle *restrict sum)
static struct sin_cos_angle sin_cos_angle_new(double sin, double cos)
@ YAC_GREAT_CIRCLE_EDGE
great circle
@ YAC_LAT_CIRCLE_EDGE
latitude circle
@ YAC_LON_CIRCLE_EDGE
longitude circle
struct sin_cos_angle inc_angle
angle between the middle point and the boundary of the spherical cap
enum yac_edge_type * edge_type
double(* coordinates_xyz)[3]