Function gmtplot_great_circle_arc

src/gmt_plot.c:4819–4847

}

GMT_LOCAL uint64_t gmtplot_great_circle_arc (struct GMT_CTRL *GMT, double *A, double *B, double step, bool longway, double **xp, double **yp) {
	/* Given vectors A and B, return great circle path sampled every step.  Shorest path is selected unless longway is true */
	/* Determine unit vector pole of great circle or use the one given by small circle pole and its opening rot */
	uint64_t k, n;
	double P[3], X[3], R[3][3], R0[3][3], c, w, *xx = NULL, *yy = NULL;

	gmt_cross3v (GMT, A, B, P);	/* Parallel to rotation pole */
	gmt_normalize3v (GMT, P);		/* Rotation pole unit vector */
	c = d_acosd (gmt_dot3v (GMT, A, B));	/* opening angle in degrees */

	if (longway) {	/* Want to go the long way */
		c = 360.0 - c;
		P[0] = -P[0], P[1] = -P[1], P[2] = -P[2];
	}
	if (gmt_M_is_zero (step)) step = GMT->current.map.path_step;	/* Use default map-step if given as 0 */
	n = lrint (ceil (c / step)) + 1;	/* Number of segments needed for smooth curve from A to B inclusive */
	step = D2R * c / (n - 1);	/* Adjust step for exact fit, convert to radians */
	gmt_M_malloc2 (GMT, xx, yy, n, NULL, double);	/* Allocate space for arrays */
	gmtlib_init_rot_matrix (R0, P);			/* Get partial rotation matrix since no actual angle is applied yet */
	for (k = 0; k < n; k++) {	/* March along the arc */
		w = k * step;					/* Opening angle from A to this point X */
		gmt_M_memcpy (R, R0, 9, double);			/* Get a copy of the "0-angle" rotation matrix */
		gmtlib_load_rot_matrix (w, R, P);			/* Build the actual rotation matrix for this angle */
		gmt_matrix_vect_mult (GMT, 3U, R, A, X);			/* Rotate point A towards B and get X */
		gmt_cart_to_geo (GMT, &yy[k], &xx[k], X, true);	/* Get lon/lat of this point along arc */
	}
	*xp = xx;	*yp = yy;
	return (n);
}

GMT_LOCAL uint64_t gmtplot_small_circle_arc (struct GMT_CTRL *GMT, double *A, double step, double P[], double rot, double **xp, double **yp) {
	/* Given vectors A and B, return small circle path sampled every step. */