30#define YAC_RAD (0.01745329251994329576923690768489)
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)
127static inline double get_angle (
double a_lon,
double b_lon) {
128 double diff = a_lon - b_lon;
130 return diff - lround(diff / (2.0 * M_PI)) * (2.0 * M_PI);
132 return diff - round(diff / (2.0 * M_PI)) * (2.0 * M_PI);
144 double diff = a_lon - b_lon;
146 return diff - lround(diff / 360.0) * (360.0);
148 return diff - round(diff / 360.0) * (360.0);
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]);
202 double ab[3] = {a[0]-b[0], a[1]-b[1], a[2]-b[2]};
203 return ab[0] * ab[0] + ab[1] * ab[1] + ab[2] * ab[2];
208 double temp = 1.0 - angle.
cos;
209 return temp * temp + angle.
sin * angle.
sin;
221 if (bnd_circle->
sq_crd == DBL_MAX)
230 return val > -1.0 ? (val < 1.0 ? val : 1.0) : -1.0;
242 double angle,
double * sin_value,
double * cos_value) {
244 double abs_angle = fabs(angle);
252 *sin_value = copysign(1.0, angle);
270static inline void LLtoXYZ(
double lon,
double lat,
double p_out[]) {
272 while (lon < -M_PI) lon += 2.0 * M_PI;
273 while (lon >= M_PI) lon -= 2.0 * M_PI;
275 double sin_lon, cos_lon, sin_lat, cos_lat;
279 p_out[0] = cos_lat * cos_lon;
280 p_out[1] = cos_lat * sin_lon;
290static inline void LLtoXYZ_deg(
double lon,
double lat,
double p_out[]) {
304static inline void XYZtoLL (
double const p_in[],
double * lon,
double * lat) {
306 *lon = atan2(p_in[1] , p_in[0]);
307 *lat = M_PI_2 - acos(p_in[2]);
311static inline double det2_kahan(
double a,
double b,
double c,
double d) {
313 double e_bc = fma(b,c,-bc);
314 double det = fma(a,d,-bc);
319 double const a[],
double const b[],
double cross[]) {
323 cross[0] =
det2_kahan(a[1], a[2], b[1], b[2]);
324 cross[1] =
det2_kahan(a[2], a[0], b[2], b[0]);
325 cross[2] =
det2_kahan(a[0], a[1], b[0], b[1]);
334 const double a[],
const double b[],
double cross[]) {
338 cross[0] = a[1] * b[2] - a[2] * b[1];
339 cross[1] = a[2] * b[0] - a[0] * b[2];
340 cross[2] = a[0] * b[1] - a[1] * b[0];
352 double dot_product = a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
355 if (fabs(dot_product) > M_SQRT1_2) {
361 double asin_tmp = asin(sqrt(cross_ab[0]*cross_ab[0]+
362 cross_ab[1]*cross_ab[1]+
363 cross_ab[2]*cross_ab[2]));
365 if (dot_product > 0.0)
366 return MAX(asin_tmp,0.0);
368 return MIN(M_PI - asin_tmp, M_PI);
372 return acos(dot_product);
398 double const a[3], double const b[3]) {
404 cross_ab[1]*cross_ab[1] +
405 cross_ab[2]*cross_ab[2]),
406 a[0]*b[0] + a[1]*b[1] + a[2]*b[2]);
419 int t_a = fabs(a.
cos) <= M_SQRT1_2;
420 int t_b = fabs(b.
cos) <= M_SQRT1_2;
421 int a_section = t_a | ((((a.
sin < 0.0) & t_a) |
422 ((a.
cos < 0.0) & (fabs(a.
sin) < M_SQRT1_2))) << 1);
423 int b_section = t_b | ((((b.
sin < 0.0) & t_b) |
424 ((b.
cos < 0.0) & (fabs(b.
sin) < M_SQRT1_2))) << 1);
426 if (!a_section) a_section = (a.
sin < 0.0) << 2;
427 if (!b_section) b_section = (b.
sin < 0.0) << 2;
429 if (a_section != b_section)
430 return (a_section > b_section) - (a_section < b_section);
444 if (a.
sin >= 0.0)
return ret;
454 if (a.
cos >= 0.0)
return ret;
464 if (a.
sin >= 0.0)
return ret;
474 if (a.
cos <= 0.0)
return ret;
489 int t = fabs(angle.
cos) <= M_SQRT1_2;
490 int section = t | ((((angle.
sin < 0.0) & t) |
491 ((angle.
cos < 0.0) & (fabs(angle.
sin) < M_SQRT1_2))) << 1);
492 if (!section) section = (angle.
sin < 0.0) << 2;
496 case(0):
return 0.0 * M_SQRT1_2 + angle.
sin;
497 case(1):
return 2.0 * M_SQRT1_2 - angle.
cos;
498 case(2):
return 4.0 * M_SQRT1_2 - angle.
sin;
499 case(3):
return 6.0 * M_SQRT1_2 + angle.
cos;
500 case(4):
return 8.0 * M_SQRT1_2 + angle.
sin;
510 a.cos * b.cos - a.sin * b.sin);
534 a.cos * b.cos + a.sin * b.sin);
550 return compare_result < 0;
557 if (angle.
cos > M_SQRT1_2) {
559 double angle_ = asin(angle.
sin);
561 if (angle_ < 0.0) angle_ += 2.0 * M_PI;
564 }
else if (angle.
cos > - M_SQRT1_2) {
566 double angle_ = acos(angle.
cos);
568 if (angle.
sin > 0.0)
return angle_;
569 else return 2.0 * M_PI - angle_;
572 return M_PI - asin(angle.
sin);
584 double x = (1.0 + fabs(angle.cos));
586 double scale = 1.0 / sqrt(x * x + angle.sin * angle.sin);
589 if (angle.cos >= 0) {
590 scale = copysign(scale, angle.sin);
601 double tan = fabs(angle.sin / (1.0 + fabs(angle.cos)));
602 double one_plus_sq_tan = 1.0 + tan * tan;
603 double sqrt_one_plus_sq_tan = sqrt(one_plus_sq_tan);
609 if (angle.cos < 0.0) {
614 double scale = M_SQRT1_2 / sqrt(one_plus_sq_tan + a * sqrt_one_plus_sq_tan);
615 double x = b * scale;
616 double y = (a + sqrt_one_plus_sq_tan) * scale;
646 double dot_product = a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
650 if (dot_product > 0.9999) {
659 if (sqrt(cross_ab[0]*cross_ab[0] +
660 cross_ab[1]*cross_ab[1] +
665 if ((ret = (a[0] > b[0]) - (a[0] < b[0])))
return ret;
666 if ((ret = (a[1] > b[1]) - (a[1] < b[1])))
return ret;
667 return (a[2] > b[2]) - (a[2] < b[2]);
673 double norm = 1.0 / sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
687 double axis[],
struct sin_cos_angle angle,
double v_in[],
double v_out[]) {
694 double cross_axis_v_in[3];
697 double dot_axis_v_in = axis[0]*v_in[0] + axis[1]*v_in[1] + axis[2]*v_in[2];
698 double temp = dot_axis_v_in * (1.0 - angle.
cos);
701 v_in[0] * angle.
cos + cross_axis_v_in[0] * angle.
sin + axis[0] * temp;
703 v_in[1] * angle.
cos + cross_axis_v_in[1] * angle.
sin + axis[1] * temp;
705 v_in[2] * angle.
cos + cross_axis_v_in[2] * angle.
sin + axis[2] * temp;
715 double axis[],
double angle,
double v_in[],
double v_out[]) {
717 double sin_angle = sin(angle);
718 double cos_angle = cos(angle);
752 size_t const * corner_indices,
size_t num_corners,
size_t start_corner,
753 size_t (*triangle_indices)[3]);
766 double * vertices,
size_t num_vertices,
double triangle[][3],
773 double p[3],
double const a[3],
double const b[3],
static struct sin_cos_angle get_vector_angle_2(double const a[3], double const b[3])
static int compare_coords(double const *a, double const *b)
static double compute_sq_crd(struct sin_cos_angle angle)
static const struct sin_cos_angle SIN_COS_LOW_TOL
static double clamp_abs_one(double val)
static const struct sin_cos_angle SIN_COS_7_M_PI_4
static const struct sin_cos_angle SIN_COS_ZERO
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)
void yac_triangulate_cell_indices(size_t const *corner_indices, size_t num_corners, size_t start_corner, size_t(*triangle_indices)[3])
static int points_are_identically(double const *a, double const *b)
static const struct sin_cos_angle SIN_COS_M_PI
static void compute_sin_cos(double angle, double *sin_value, double *cos_value)
void yac_compute_bnd_triangle(double *vertices, size_t num_vertices, double triangle[][3], size_t num_tests)
static int yac_point_in_bounding_circle_vec(double point_vector[3], struct bounding_circle *bnd_circle)
static double sq_len_diff_vec(double const a[3], double const b[3])
computes square of the lenght of the vector ab
static void crossproduct_d(const double a[], const double b[], double cross[])
static struct sin_cos_angle sum_angles_no_check(struct sin_cos_angle a, struct sin_cos_angle b)
static double sin_cos_angle_to_dble(struct sin_cos_angle angle)
static void crossproduct_kahan(double const a[], double const b[], double cross[])
static void rotate_vector(double axis[], double angle, double v_in[], double v_out[])
static struct sin_cos_angle quarter_angle(struct sin_cos_angle angle)
static const struct sin_cos_angle SIN_COS_TOL
static const struct sin_cos_angle SIN_COS_M_PI_2
#define yac_cos_angle_low_tol
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[])
#define yac_cos_angle_tol
void yac_triangulate_cell(struct yac_grid_cell cell, size_t start_corner, struct yac_grid_cell *triangles)
static void LLtoXYZ_deg(double lon, double lat, double p_out[])
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)
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)
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)
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 double compute_angle(struct sin_cos_angle angle)
return angles in the range of [0;2PI[
#define yac_angle_low_tol
static void normalise_vector(double v[])
static double get_angle(double a_lon, double b_lon)
static void LLtoXYZ(double lon, double lat, double p_out[])
static double get_vector_angle(double const a[3], double const b[3])
static int sum_angles(struct sin_cos_angle a, struct sin_cos_angle b, struct sin_cos_angle *restrict sum)
static double get_angle_deg(double a_lon, double b_lon)
static struct sin_cos_angle sin_cos_angle_new(double sin, double cos)
static struct sin_cos_angle half_angle(struct sin_cos_angle angle)
static void XYZtoLL(double const p_in[], double *lon, double *lat)
@ 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
struct yac_circle::@2::@4 lat
struct yac_circle::@2::@3 lon
enum yac_circle_type type
struct yac_circle::@2::@3 gc
struct yac_circle::@2::@5 p
union yac_circle::@2 data