13static double const tol = 1.0e-12;
40 double const norm_vector_a[3],
double const norm_vector_b[3],
41 double p[3],
double q[3]) {
50 double sq_sin_plane_angle =
51 temp_p[0] * temp_p[0] + temp_p[1] * temp_p[1] + temp_p[2] * temp_p[2];
56 double scale = 1.0 / sqrt(sq_sin_plane_angle);
57 p[0] = (temp_p[0] * scale);
58 p[1] = (temp_p[1] * scale);
59 p[2] = (temp_p[2] * scale);
60 q[0] = -p[0], q[1] = -p[1], q[2] = -p[2];
65#if defined __INTEL_COMPILER
66#pragma intel optimization_level 0
72 double const gc_norm_vector[3],
double const z,
73 double p[3],
double q[3]) {
80 double b = sqrt(
MAX(1.0 - sq_z,0.0));
81 double b_n0 = b * gc_norm_vector[0];
82 double b_n1 = b * gc_norm_vector[1];
95 double const gc_norm_vector[3],
double const z,
96 double p[3],
double q[3]) {
131 return loncxlatc(gc_norm_vector, z, p, q);
133 double sq_n[3] = {gc_norm_vector[0] * gc_norm_vector[0],
134 gc_norm_vector[1] * gc_norm_vector[1],
135 gc_norm_vector[2] * gc_norm_vector[2]};
138 double max_sq_gc_z = sq_n[0] + sq_n[1];
147 int num_intersections =
148 ((max_sq_gc_z > (sq_z -
tol)) - (max_sq_gc_z < (sq_z +
tol))) + 1;
150 switch (num_intersections) {
159 double a = -(gc_norm_vector[2] * z) / max_sq_gc_z;
161 p[0] = ((q[0] = gc_norm_vector[0] * a));
162 p[1] = ((q[1] = gc_norm_vector[1] * a));
170 double a = -(gc_norm_vector[2] * z) / max_sq_gc_z;
171 double b = sqrt(max_sq_gc_z - sq_z)/max_sq_gc_z;
172 double a_n0 = a * gc_norm_vector[0];
173 double a_n1 = a * gc_norm_vector[1];
174 double b_n0 = b * gc_norm_vector[0];
175 double b_n1 = b * gc_norm_vector[1];
187 return num_intersections;
195 double const norm_vector_a[3],
double const norm_vector_b[3],
196 double p[3],
double q[3]) {
198 double sq_plane_diff =
199 MIN((norm_vector_a[0] - norm_vector_b[0]) *
200 (norm_vector_a[0] - norm_vector_b[0]) +
201 (norm_vector_a[1] - norm_vector_b[1]) *
202 (norm_vector_a[1] - norm_vector_b[1]),
203 (norm_vector_a[0] + norm_vector_b[0]) *
204 (norm_vector_a[0] + norm_vector_b[0]) +
205 (norm_vector_a[1] + norm_vector_b[1]) *
206 (norm_vector_a[1] + norm_vector_b[1]));
212 p[0] = 0, p[1] = 0; p[2] = 1;
213 q[0] = 0, q[1] = 0; q[2] = -1;
223 double const vec[3],
double const norm_vector[3],
double p[3],
double q[3]) {
227 double dot = vec[0] * norm_vector[0];
228 double err = fma(vec[0], norm_vector[0], -dot);
229 dot = fma(vec[1], norm_vector[1], dot);
230 err += fma(vec[1], norm_vector[1], -vec[1]*norm_vector[1]);
231 dot = fma(vec[2], norm_vector[2], dot);
232 err += fma(vec[2], norm_vector[2], -vec[2]*norm_vector[2]);
238 q[0] = -((p[0] = vec[0]));
239 q[1] = -((p[1] = vec[1]));
240 q[2] = -((p[2] = vec[2]));
245 double const vec[3],
double z,
double p[3],
double q[3]) {
248 if (fabs(vec[2] - z) >
tol)
return 0;
250 q[0] = -((p[0] = vec[0]));
251 q[1] = -((p[1] = vec[1]));
252 q[2] = ((p[2] = vec[2]));
257 double const vec_a[3],
double const vec_b[3],
double p[3],
double q[3]) {
261 q[0] = -((p[0] = vec_a[0]));
262 q[1] = -((p[1] = vec_a[1]));
263 q[2] = -((p[2] = vec_a[2]));
282 "ERROR(yac_circle_intersect): unknown circle type.")
320 double const a[],
double const b[],
double const p[],
321 double sq_len_diff_ab) {
333 return sq_len_diff_ab +
tol >= sq_len_diff_vec_ap + sq_len_diff_vec_bp;
337 double const a[],
double const b[],
double const p[]) {
347 if (fabs(fabs(a[2]) - 1.0) <
tol)
return 1;
349 double fabs_a[2] = {fabs(a[0]), fabs(a[1])};
350 int flag = fabs_a[0] > fabs_a[1];
351 double a_0 = a[flag], a_1 = a[flag^1];
352 double b_0 = b[flag], b_1 = b[flag^1];
353 double p_0 = p[flag], p_1 = p[flag^1];
355 double temp = b_0 - (b_1 * a_0) / a_1;
359 "ERROR(vector_is_between_lat): "
360 "routine does not support zero length edges")
365 double beta = (p_0 - (p_1 * a_0) / a_1) / temp;
367 if (beta < -
tol)
return 0;
369 double alpha = (p_1 - b_1 * beta) / a_1;
375 enum yac_edge_type edge_type,
double const a[3],
double const b[3],
382 "ERROR(compute_edge_middle_point_vec): invalide edge_type")
385 middle[0] = a[0]+b[0], middle[1] = a[1]+b[1], middle[2] = a[2];
386 double temp = middle[0]*middle[0] + middle[1]*middle[1];
387 double scale = sqrt((1.0 - middle[2] * middle[2])/temp);
391 middle[0] = a[0]+b[0], middle[1] = a[1]+b[1], middle[2] = a[2]+b[2];
397 enum yac_edge_type edge_type,
double const a[3],
double const b[3],
398 double const c[3],
double const d[3],
double p[3],
double q[3],
399 int intersection_type) {
417 switch (intersection_type) {
433 p[0] = a[0], p[1] = a[1], p[2] = a[2];
434 q[0] = b[0], q[1] = b[1], q[2] = b[2];
441 p[0] = c[0], p[1] = c[1], p[2] = c[2];
442 q[0] = d[0], q[1] = d[1], q[2] = d[2];
447 p[0] = a[0], p[1] = a[1], p[2] = a[2];
448 q[0] = c[0], q[1] = c[1], q[2] = c[2];
451 p[0] = a[0], p[1] = a[1], p[2] = a[2];
452 q[0] = d[0], q[1] = d[1], q[2] = d[2];
455 p[0] = b[0], p[1] = b[1], p[2] = b[2];
456 q[0] = c[0], q[1] = c[1], q[2] = c[2];
459 p[0] = b[0], p[1] = b[1], p[2] = b[2];
460 q[0] = d[0], q[1] = d[1], q[2] = d[2];
469 double const a[3],
double const b[3],
double const c[3],
double const d[3],
470 double p[3],
double q[3]) {
486 double const a[3],
double const b[3],
double const c[3],
double const d[3],
487 double const p[3],
double const q[3]) {
500 double const a[3],
double const b[3],
double norm_vector[3]) {
502 norm_vector[0] = a[1] * b[2] - a[2] * b[1];
503 norm_vector[1] = a[2] * b[0] - a[0] * b[2];
504 norm_vector[2] = a[0] * b[1] - a[1] * b[0];
506 double scale = 1.0 / sqrt(norm_vector[0] * norm_vector[0] +
507 norm_vector[1] * norm_vector[1] +
508 norm_vector[2] * norm_vector[2]);
509 norm_vector[0] *= scale;
510 norm_vector[1] *= scale;
511 norm_vector[2] *= scale;
536 double const a[3],
double const b[3],
double const c[3],
double const d[3],
537 double p[3],
double q[3]) {
539 double norm_vector_ab[3], norm_vector_cd[3];
544 switch (
gcxgc(norm_vector_ab, norm_vector_cd, p, q)) {
571 double const a[3],
double const b[3],
double const c[3],
double const d[3],
572 double p[3],
double q[3]) {
575 if (fabs(a[2] - c[2]) >
tol)
return -1;
587 double const a[3],
double const b[3],
double norm_vector[3]) {
589 norm_vector[0] = a[1] * b[2] - a[2] * b[1];
590 norm_vector[1] = a[2] * b[0] - a[0] * b[2];
591 norm_vector[2] = 0.0;
593 double scale = 1.0 / sqrt(norm_vector[0] * norm_vector[0] +
594 norm_vector[1] * norm_vector[1]);
595 norm_vector[0] *= scale;
596 norm_vector[1] *= scale;
617 double const a[3],
double const b[3],
double const c[3],
double const d[3],
618 double p[3],
double q[3]) {
620 double norm_vector_ab[3], norm_vector_cd[3];
624 switch(
loncxlonc(norm_vector_ab, norm_vector_cd, p, q)) {
636 double const a[3],
double const b[3],
double const c[3],
double const d[3],
637 double const p[3],
double const q[3]) {
667 double const a[3],
double const b[3],
double const c[3],
double const d[3],
668 double p[3],
double q[3]) {
670 double norm_vector_ab[3];
673 int num_intersections =
loncxlatc(norm_vector_ab, c[2], p, q);
674 YAC_ASSERT(num_intersections == 2,
"ERROR(yac_loncxlatc_vec): internal error")
702#if defined __INTEL_COMPILER
703#pragma intel optimization_level 0
709 double const a[3],
double const b[3],
double const c[3],
double const d[3],
710 double p[3],
double q[3]) {
712 double norm_vector_ab[3];
715 int num_intersections =
gcxlatc(norm_vector_ab, c[2], p, q);
717 switch (num_intersections) {
734 double const a[3],
double const b[3],
double p[3],
double q[3]) {
740 double const point[3],
double const a[3],
double const b[3],
741 double p[3],
double q[3]) {
748 if (
pxgc(point, norm_ab, p, q)) {
764 double const point[3],
double const a[3],
double const b[3],
765 double p[3],
double q[3]) {
767 if (
pxlatc(point, a[2], p, q)) {
781 return (ret_value & (~(1 + 2 + 4 + 8))) +
782 ((ret_value & (1 + 2)) << 2) +
783 ((ret_value & (4 + 8)) >> 2);
787 enum yac_edge_type edge_type_a,
double const a[3],
double const b[3],
788 enum yac_edge_type edge_type_b,
double const c[3],
double const d[3],
789 double p[3],
double q[3]) {
803 if (edge_a_is_point && edge_b_is_point) {
808 }
else if (edge_a_is_point || edge_b_is_point) {
810 enum yac_edge_type edge_type[2] = {edge_type_a, edge_type_b};
811 double const * edge[2][2] = {{a,b}, {c,d}};
812 double const * point[2] = {c,a};
815 switch (edge_type[edge_a_is_point]) {
819 point[edge_a_is_point],
820 edge[edge_a_is_point][0], edge[edge_a_is_point][1], p, q);
827 point[edge_a_is_point],
828 edge[edge_a_is_point][0], edge[edge_a_is_point][1], p, q);
832 if (edge_a_is_point)
return ret_value;
884 double p[3],
double const a[3],
double const b[3],
891 (circle_type ==
POINT),
892 "ERROR(yac_point_on_edge): unknown edge type.")
894 switch (circle_type) {
static int points_are_identically(double const *a, double const *b)
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_kahan(double const a[], double const b[], double cross[])
#define yac_angle_low_tol
static void normalise_vector(double v[])
@ YAC_GREAT_CIRCLE_EDGE
great circle
@ YAC_LAT_CIRCLE_EDGE
latitude circle
@ YAC_LON_CIRCLE_EDGE
longitude circle
static int yac_pxlat_vec(double const point[3], double const a[3], double const b[3], double p[3], double q[3])
static int yac_check_pq_gcxlatc(double const a[3], double const b[3], double const c[3], double const d[3], double const p[3], double const q[3])
static int yac_latcxlatc_vec(double const a[3], double const b[3], double const c[3], double const d[3], double p[3], double q[3])
compute the intersection point two circles of latitude
static int gcxgc(double const norm_vector_a[3], double const norm_vector_b[3], double p[3], double q[3])
int yac_circle_intersect(struct yac_circle a, struct yac_circle b, double p[3], double q[3])
static int yac_identical_gcxgc_vec(double const a[3], double const b[3], double const c[3], double const d[3], double p[3], double q[3])
static int vector_is_between(double const a[], double const b[], double const p[], double sq_len_diff_ab)
static int yac_pxgc_vec(double const point[3], double const a[3], double const b[3], double p[3], double q[3])
static void compute_edge_middle_point_vec(enum yac_edge_type edge_type, double const a[3], double const b[3], double middle[3])
static int gcxlatc(double const gc_norm_vector[3], double const z, double p[3], double q[3])
static int vector_is_between_lat(double const a[], double const b[], double const p[])
static int yac_check_pq_gcxgc(double const a[3], double const b[3], double const c[3], double const d[3], double const p[3], double const q[3])
static int pxgc(double const vec[3], double const norm_vector[3], double p[3], double q[3])
static int adjust_ret_value(int ret_value)
static void compute_norm_vector(double const a[3], double const b[3], double norm_vector[3])
static int yac_gcxgc_vec(double const a[3], double const b[3], double const c[3], double const d[3], double p[3], double q[3])
static int loncxlonc(double const norm_vector_a[3], double const norm_vector_b[3], double p[3], double q[3])
static int yac_loncxlatc_vec(double const a[3], double const b[3], double const c[3], double const d[3], double p[3], double q[3])
compute the intersection point of a meridian and a parallel
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])
int yac_point_on_edge(double p[3], double const a[3], double const b[3], enum yac_circle_type circle_type)
static int pxlatc(double const vec[3], double z, double p[3], double q[3])
static int yac_gcxlatc_vec(double const a[3], double const b[3], double const c[3], double const d[3], double p[3], double q[3])
compute the intersection of a great circle with the parallel
static int yac_pxp_vec(double const a[3], double const b[3], double p[3], double q[3])
static int yac_loncxlonc_vec(double const a[3], double const b[3], double const c[3], double const d[3], double p[3], double q[3])
compute the intersection point two circles of longitude
static int yac_identical_circles_vec(enum yac_edge_type edge_type, double const a[3], double const b[3], double const c[3], double const d[3], double p[3], double q[3], int intersection_type)
static void compute_lon_norm_vector(double const a[3], double const b[3], double norm_vector[3])
static int loncxlatc(double const gc_norm_vector[3], double const z, double p[3], double q[3])
static int pxp(double const vec_a[3], double const vec_b[3], double p[3], double q[3])
struct yac_circle::@3::@4 gc
union yac_circle::@3 data
enum yac_circle_type type
struct yac_circle::@3::@6 p
struct yac_circle::@3::@5 lat
struct yac_circle::@3::@4 lon
#define YAC_ASSERT(exp, msg)