Function gmtplot_scircle_sub

src/gmt_plot.c:4984–5045

}

GMT_LOCAL void gmtplot_scircle_sub (struct GMT_CTRL *GMT, double lon0, double lat0, double angle_1, double angle_2, struct GMT_SYMBOL *S, struct GMT_CIRCLE *C) {
	/* We must determine points A and B, whose small-circle connector about pole P is the arc we seek to draw */

	int justify = PSL_vec_justify (S->v.status);	/* Return justification as 0-3 */
	double R[3][3], M[3];
	gmt_M_memset (C, 1, struct GMT_CIRCLE);	/* Set all to zero */
	/* Requires the rotation matrix for pole S->v.pole */

	/* Here angle_1, angle_2 are not necessarily that, depending on S->v.status:
	 * S->v.pole & PSL_VEC_ANGLES : angle_1 is opening angle1 and angle_2 is opening angle2 about the pole.
	 * Otherwise:	angle_2 is the length of the arc in km */
	gmt_geo_to_cart (GMT, lat0, lon0, M, true);	/* Given input point */
	gmt_geo_to_cart (GMT, S->v.pole[GMT_Y], S->v.pole[GMT_X], C->P, true);
	C->colat = d_acosd (gmt_dot3v (GMT, M, C->P));	/* Colatitude of input point relative to pole, in degrees */

	if (S->v.status & PSL_VEC_ANGLES) {
		/* Was given the two opening angles; compute A and B accordingly */
		gmt_make_rot_matrix (GMT, S->v.pole[GMT_X], S->v.pole[GMT_Y], angle_1, R);
		gmt_matrix_vect_mult (GMT, 3U, R, M, C->A);	/* Get A */
		gmt_cart_to_geo (GMT, &C->lat[0], &C->lon[0], C->A, true);
		gmt_make_rot_matrix (GMT, S->v.pole[GMT_X], S->v.pole[GMT_Y], angle_2, R);
		gmt_matrix_vect_mult (GMT, 3U, R, M, C->B);	/* Get B */
		gmt_cart_to_geo (GMT, &C->lat[1], &C->lon[1], C->B, true);
		C->rot = C->r0 = C->r = angle_2 - angle_1;
	}
	else {
		/* Here A, B, or midpoint was given, + the arc length via angle_2 */
		/* Determine co-latitude for this point */
		C->rot = C->r0 = C->r = (angle_1 / GMT->current.proj.DIST_KM_PR_DEG) / sind (C->colat);	/* Opening angle in spherical degrees */
		switch (justify) {	/* A and B depends on chosen justification */
			case 0: /* Was given coordinates of A; determine B */
				gmt_M_memcpy (C->A, M, 3, double);
				C->lon[0] = lon0;	C->lat[0] = lat0;
				if (C->r > 180.0) {C->longway = true; C->r -= 180.0;}	/* Temporarily adjust if arcs > 180 degrees are chosen */
				/* Rotate A by C->r0 degrees about P to get B */
				gmt_make_rot_matrix (GMT, S->v.pole[GMT_X], S->v.pole[GMT_Y], C->r0, R);
				gmt_matrix_vect_mult (GMT, 3U, R, C->A, C->B);	/* Get B */
				gmt_cart_to_geo (GMT, &C->lat[1], &C->lon[1], C->B, true);
				break;
			case 1: /* Was given coordinates of halfway point; determine A and B */
				if (C->r > 180.0) C->longway = true;	/* Temporarily adjust if arcs > 180 degrees are chosen */
				/* Rotate M by -C->r0/2 degrees about P to get A */
				gmt_make_rot_matrix (GMT, S->v.pole[GMT_X], S->v.pole[GMT_Y], -0.5 * C->r0, R);
				gmt_matrix_vect_mult (GMT, 3U, R, M, C->A);	/* Get A */
				gmt_cart_to_geo (GMT, &C->lat[0], &C->lon[0], C->A, true);
				/* Rotate M by +C->r0/2 degrees about P to get B */
				gmt_make_rot_matrix (GMT, S->v.pole[GMT_X], S->v.pole[GMT_Y], +0.5 * C->r0, R);
				gmt_matrix_vect_mult (GMT, 3U, R, M, C->B);	/* Get B */
				gmt_cart_to_geo (GMT, &C->lat[1], &C->lon[1], C->B, true);
				break;
			case 2: /* Was given coordinates of B point; determine A */
				gmt_M_memcpy (C->B, M, 3, double);
				C->lon[1] = lon0;	C->lat[1] = lat0;
				if (C->r > 180.0) {C->longway = true; C->r -= 180.0;}	/* Temporarily adjust if arcs > 180 degrees are chosen */
				/* Rotate B by -C->r0 degrees about P to get A */
				gmt_make_rot_matrix (GMT, S->v.pole[GMT_X], S->v.pole[GMT_Y], -C->r0, R);
				gmt_matrix_vect_mult (GMT, 3U, R, C->B, C->A);	/* Get A */
				gmt_cart_to_geo (GMT, &C->lat[0], &C->lon[0], C->A, true);