Function gmtplot_circle_pen_poly

src/gmt_plot.c:4883–4927

}

GMT_LOCAL void gmtplot_circle_pen_poly (struct GMT_CTRL *GMT, double *A, double *B, bool longway, double rot, struct GMT_PEN *pen, struct GMT_SYMBOL *S, struct GMT_CIRCLE *C, double scale) {
	/* Given vectors A and B, return a small circle polygon path sampled every step that approximates a pen of given width
	 * drawn on the map.  Use small circle pole and its opening rot */
	uint64_t k, n, n2;
	double Ai[3], Ao[3], Px[3], X[3], R[3][3], R0[3][3], w, step;
	struct GMT_DATASEGMENT *L = NULL;
	gmt_M_unused(B);

	gmt_cross3v (GMT, A, C->P, Px);	/* Px is Pole to plane through A and P  */
	gmt_normalize3v (GMT, Px);			/* Rotation pole unit vector */
	/* Rotate A back/fore by rotation angle of +/- pen halfwidth about Px */
	w = 0.5 * scale * pen->width * GMT->session.u2u[GMT_PT][GMT_INCH] * S->v.scale;		/* Half-width of pen in degrees */
	gmt_make_rot_matrix2 (GMT, Px, +w, R);		/* Rotation of rot_v degrees about pole P */
	gmt_matrix_vect_mult (GMT, 3U, R, A, Ai);	/* Get Ai = R * A */
	gmt_make_rot_matrix2 (GMT, Px, -w, R);		/* Rotation of rot_v degrees about pole P */
	gmt_matrix_vect_mult (GMT, 3U, R, A, Ao);	/* Get Ao = R * A */
	if (longway) rot = 360.0 - rot;

	step = GMT->current.map.path_step;		/* Use default map-step if given as 0 */
	n = lrint (ceil (fabs (rot) / step)) + 1;	/* Number of segments needed for smooth curve from A to B inclusive */
	step = D2R * rot / (n - 1);			/* Adjust step for exact fit, convert to radians */
	L = GMT_Alloc_Segment (GMT->parent, GMT_NO_STRINGS, 2*n+1, 2, NULL, NULL);	/* Allocate polygon to draw filled path */
	n2 = 2*n-1;
	gmtlib_init_rot_matrix (R0, C->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, 9U, double);			/* Get a copy of the "0-angle" rotation matrix */
		gmtlib_load_rot_matrix (w, R, C->P);			/* Build the actual rotation matrix for this angle */
		gmt_matrix_vect_mult (GMT, 3U, R, Ai, X);		/* Rotate point Ai towards B and get X */
		gmt_cart_to_geo (GMT, &L->data[GMT_Y][k], &L->data[GMT_X][k], X, true);	/* Get lon/lat of this point along arc */
		gmt_matrix_vect_mult (GMT, 3U, R, Ao, X);		/* Rotate point Ai towards B and get X */
		gmt_cart_to_geo (GMT, &L->data[GMT_Y][n2-k], &L->data[GMT_X][n2-k], X, true);	/* Get lon/lat of this point along arc */
	}
	L->data[GMT_X][2*n] = L->data[GMT_X][0];	/* Explicitly close the polygon */
	L->data[GMT_Y][2*n] = L->data[GMT_Y][0];

	/* Plot pen as a closed filled polygon without outline */
	PSL_command (GMT->PSL, "V\n");
	gmt_setpen (GMT, pen);		/* Set pen width just so later setpen's will work */
	PSL_setfill (GMT->PSL, pen->rgb, 0);	/* Fill, no outline */
	gmt_geo_polygons (GMT, L);	/* "Draw" the line */
	PSL_command (GMT->PSL, "U\n");

	gmt_free_segment (GMT, &L);
}

GMT_LOCAL void gmtplot_gcircle_sub (struct GMT_CTRL *GMT, double lon0, double lat0, double azimuth, double length, struct GMT_SYMBOL *S, struct GMT_CIRCLE *C) {
	/* We must determine points A and B, whose great-circle connector is the arc we seek to draw */