Add the so wished option to plot spheres in psxyz. (#8844)

* Add the so wished option to plot spheres in psxyz.

They aren't really spheres but circles in 3D plots.

This is mostly a Claude.ai creation but a one that needed a lot of baby-sitting development.

* Change -SQ to -SP for sphere plots. Add a sphere.sh test

* Add the baseline image sphere.ps to dvc

* Update src/psxyz.c

Co-authored-by: Dongdong Tian <seisman.info@gmail.com>

* Fixes and updates to this test.

* Make the -SP+e & +a work better.

Oddly the PS seem right but conversion to raster flips spheres upside-down.

---------

Co-authored-by: Dongdong Tian <seisman.info@gmail.com>

Author: joa-quim

doc/rst/source/explain_3D_symbols.rst_

@@ -31,3 +31,17 @@
         **q** if *size* is a quantity in x-units [Default is plot-distance units].
         The facet colors will be modified to simulate shading.
         Use **-SU** to disable 3-D illumination.
+
+    **-SP**\ *size*\ [**c**\|\ **i**\|\ **p**][**+a**\ *azimuth*][**+e**\ *elevation*][**+f**][**+n**]
+        Sphere (3-D) with diameter *size* [Default unit is **c** (cm)].
+        The sphere is rendered with a radial gradient from white at the light source
+        to the fill color (**-G**) at the opposite side, creating a 3-D appearance.
+        Modifiers:
+
+        - **+a**\ *azimuth* - Set light source azimuth [default: 0, from the right].
+        - **+e**\ *elevation* - Set light source elevation [default: 90, perpendicular to viewing plane].
+        - **+f** - Use flat/constant fill color (no gradient shading).
+        - **+n** - Draw outline only (no fill).
+
+        The outline pen is controlled by **-W**. If no fill color is specified via **-G**,
+        the sphere defaults to black.

src/PSL_prologue.ps

@@ -1038,6 +1038,35 @@ static char *PSL_prologue_str =
 "/Sp {N 3 -1 roll 0 360 arc fs N}!\n"
 "% Patch fill: x1 y1 ... xn yn n\n"
 "/SP {M {D} repeat FO}!\n"
+"% Sphere with 3D shading: radius xc yc\n"
+"/SPhere {/yc_SP edef /xc_SP edef /radius_SP edef\n"
+"  /SP_flat where {pop} {/SP_flat false def} ifelse\n"
+"  /SP_no_fill where {pop} {/SP_no_fill false def} ifelse\n"
+"  gsave\n"
+"  xc_SP yc_SP matrix currentmatrix transform /yc_SP_user edef /xc_SP_user edef\n"
+"  radius_SP 0 matrix currentmatrix dtransform dup mul exch dup mul add sqrt /radius_SP_user edef\n"
+"  /SP_lx where {pop /SP_lx_orig SP_lx def /SP_ly_orig SP_ly def} {/SP_lx_orig 0.0 def /SP_ly_orig 0.0 def} ifelse\n"
+"  radius_SP_user SP_lx_orig mul /SP_lx_user edef  radius_SP_user SP_ly_orig mul /SP_ly_user edef\n"
+"  matrix setmatrix\n"
+"  SP_no_fill {\n"
+"    % No fill (outline only, drawn separately in C code)\n"
+"  } {\n"
+"    fc\n"
+"    SP_flat {\n"
+"      % Flat/constant color (no gradient)\n"
+"      xc_SP_user yc_SP_user T N 0 0 radius_SP_user 0 360 arc fill\n"
+"    } {\n"
+"      % 3D gradient shading\n"
+"      currentrgbcolor /B_SP edef /G_SP edef /R_SP edef\n"
+"      xc_SP_user yc_SP_user T N 0 0 radius_SP_user 0 360 arc\n"
+"      /Pattern setcolorspace\n"
+"      << /PatternType 2 /Shading << /ShadingType 3 /ColorSpace /DeviceRGB\n"
+"      /Coords [SP_lx_user SP_ly_user 0 0 0 radius_SP_user] /Function <<\n"
+"      /FunctionType 2 /Domain [0 1] /C0 [1 1 1] /C1 [R_SP G_SP B_SP] /N 1 >> >>\n"
+"      >> matrix makepattern setcolor fill\n"
+"    } ifelse\n"
+"  } ifelse\n"
+"  grestore}!\n"
 "% Rectangle: height width xc yc\n"
 "/Sr {M dup -2 div 2 index -2 div G dup 0 D exch 0 exch D neg 0 D FO}!\n"
 "% Rounded rectangle: height width radius xc yc\n"

src/PSL_strings.h

@@ -17118,7 +17118,7 @@ int gmt_parse_symbol_option (struct GMT_CTRL *GMT, char *text, struct GMT_SYMBOL
 	unsigned int ju, col;
 	char symbol_type, txt_a[GMT_LEN256] = {""}, txt_b[GMT_LEN256] = {""}, txt_c[GMT_LEN256] = {""}, txt_d[GMT_LEN256] = {""};
 	char text_cp[GMT_LEN256] = {""}, diameter[GMT_LEN32] = {""}, *c = NULL;
-	static char *allowed_symbols[2] = {"~=-+AaBbCcDdEefGgHhIiJjMmNnpqRrSsTtVvWwxy", "=-+AabCcDdEefGgHhIiJjMmNnOopqRrSsTtUuVvWwxy"};
+	static char *allowed_symbols[2] = {"~=-+AaBbCcDdEefGgHhIiJjMmNnpQqRrSsTtVvWwxy", "=-+AabCcDdEefGgHhIiJjMmNnOopQqRrSsTtUuVvWwxy"};
 	static char *bar_symbols[2] = {"Bb", "-BbOoUu"};
 	if (cmd) {
 		p->base = GMT->session.d_NaN;
@@ -17213,6 +17213,28 @@ int gmt_parse_symbol_option (struct GMT_CTRL *GMT, char *text, struct GMT_SYMBOL
 			}
 		}
 	}
+	else if (text[0] == 'P') {	/* Sphere symbol with optional modifiers for light position, flat color, or no fill */
+		char arg[GMT_LEN64] = {""};
+		n = sscanf (text, "%c%[^+]", &symbol_type, arg);	/* arg should be symbol size with no +<modifiers> at the end */
+		if (n == 1) {	/* No modifiers or no size given */
+			if (text[1] && text[1] != '+') {
+				/* Gave size without modifiers */
+				strncpy (arg, &text[1], GMT_LEN64-1);
+			}
+		}
+		if (arg[0] && arg[0] != '+') {	/* Need to get size */
+			if (cmd) p->read_size_cmd = false;
+			p->size_x = p->given_size_x = gmt_M_to_inch (GMT, arg);
+			check = false;
+		}
+		else if (!text[1] || text[1] == '+') {	/* No size given */
+			if (p->size_x == 0.0) p->size_x = p->given_size_x;
+			if (p->size_y == 0.0) p->size_y = p->given_size_y;
+			if (p->size_x == 0.0)
+				col_off++;
+			if (cmd) p->read_size_cmd = true;
+		}
+	}
 	else if (strchr (GMT_VECTOR_CODES, text[0])) {
 		/* Vectors gets separate treatment because of optional modifiers [+j<just>+b+e+s+l+r+a<angle>+n<norm>] */
 		int one;
@@ -17823,14 +17845,67 @@ int gmt_parse_symbol_option (struct GMT_CTRL *GMT, char *text, struct GMT_SYMBOL
 				GMT_Report (GMT->parent, GMT_MSG_ERROR, "Option -S: Symbol type %c is 3-D only\n", symbol_type);
 			}
 			break;
-		case 'P':
 		case 'p':
 			p->symbol = PSL_DOT;
 			if (p->size_x == 0.0 && !p->read_size) {	/* User forgot to set size */
 				p->size_x = GMT_DOT_SIZE;
 				check = false;
 			}
 			break;
+		case 'P':	/* Sphere symbol: -SP<size>[+a<azim>][+e<elev>][+f][+n] */
+			p->symbol = PSL_SPHERE;
+			/* Set default light position: center (perpendicular to viewing plane) */
+			p->SP_lx = 0.0;
+			p->SP_ly = 0.0;
+			p->SP_light_set = false;
+			p->SP_flat = false;
+			p->SP_no_fill = false;
+			/* Process +a, +e, +f, and +n modifiers */
+			{
+				double azimuth = 0.0, elevation = 90.0;	/* Default values: centered light */
+				bool got_azim = false, got_elev = false;
+				char mod[GMT_LEN64] = {""};
+				unsigned int pos = 0, error = 0;
+				if ((c = strchr(text, '+'))) {	/* Got modifiers */
+					while (gmt_getmodopt(GMT, 'S', c, "aefn", &pos, mod, &error) && error == 0) {
+						switch (mod[0]) {
+							case 'a':	/* Azimuth */
+								if (mod[1]) {
+									azimuth = atof(&mod[1]);
+									got_azim = true;
+								}
+								else {
+									GMT_Report(GMT->parent, GMT_MSG_ERROR, "Option -SP: +a modifier requires azimuth value\n");
+									decode_error++;
+								}
+								break;
+							case 'e':	/* Elevation */
+								if (mod[1]) {
+									elevation = atof(&mod[1]);
+									got_elev = true;
+								}
+								else {
+									GMT_Report(GMT->parent, GMT_MSG_ERROR, "Option -SP: +e modifier requires elevation value\n");
+									decode_error++;
+								}
+								break;
+							case 'f':	/* Flat/constant color (no gradient) */
+								p->SP_flat = true;
+								break;
+							case 'n':	/* No fill (outline only) */
+								p->SP_no_fill = true;
+								break;
+						}
+					}
+					if (got_azim || got_elev) {	/* Store light azimuth/elevation for later projection */
+						/* Store the light direction - will be projected to 2D at drawing time */
+						p->SP_light_az = azimuth;
+						p->SP_light_el = elevation;
+						p->SP_light_set = true;
+					}
+				}
+			}
+			break;
 		case 'q':	/* Quoted lines: -Sq[d|n|l|s|x]<info>[:<labelinfo>] */
 			p->symbol = GMT_SYMBOL_QUOTED_LINE;
 			check = false;

src/gmt_init.c

@@ -190,6 +190,15 @@ struct GMT_SYMBOL {
 	struct GMT_FRONTLINE f;	/* parameters needed for a front */
 	struct GMT_CUSTOM_SYMBOL *custom;	/* pointer to a custom symbol */
 
+	/* These apply to sphere symbols */
+	double SP_lx;		/* Light source x position for sphere symbol [0.0] */
+	double SP_ly;		/* Light source y position for sphere symbol [0.0] */
+	double SP_light_az;	/* Light source azimuth for sphere symbol [0.0] */
+	double SP_light_el;	/* Light source elevation for sphere symbol [0.0] */
+	bool SP_light_set;	/* true if +a or +e modifier was used to set light position */
+	bool SP_flat;		/* true if +f modifier was used to disable gradient (flat/solid color) */
+	bool SP_no_fill;	/* true if +n modifier was used to draw outline only (no fill) */
+
 	struct GMT_CONTOUR G;	/* For quoted lines */
 	struct GMT_DECORATE D;	/* For decorated lines */
 };

src/gmt_plot.h

@@ -3995,6 +3995,10 @@ int PSL_plotsymbol (struct PSL_CTRL *PSL, double x, double y, double size[], int
 			PSL_command (PSL, "%d %d %d S%c\n", psl_iz (PSL, 0.5 * size[0]), psl_ix (PSL, x), psl_iy (PSL, y), (char)symbol);
 			break;
 
+	case PSL_SPHERE:	/* Sphere */
+		PSL_command (PSL, "%d %d %d SPhere\n", psl_iz (PSL, 0.5 * size[0]), psl_ix (PSL, x), psl_iy (PSL, y));
+		break;
+
 		/* Multi-parameter fillable symbols */
 
 		case PSL_WEDGE:		/* A wedge or pie-slice. size[0] = radius, size[1..2] = azimuth range of arc */

src/postscriptlight.c

@@ -69,6 +69,7 @@ extern "C" {
 #define PSL_MARC		((int)'m')
 #define PSL_PENTAGON		((int)'n')
 #define PSL_DOT			((int)'p')
+#define PSL_SPHERE		((int)'P')
 #define PSL_RECT		((int)'r')
 #define PSL_RNDRECT		((int)'R')
 #define PSL_SQUARE		((int)'s')

src/postscriptlight.h

@@ -616,7 +616,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
 	GMT_Usage (API, 2, "\n%s Basic geometric symbol. Append one:", GMT_LINE_BULLET);
 	GMT_Usage (API, -3, "-(xdash), +(plus), st(a)r, (b|B)ar, (c)ircle, (d)iamond, (e)llipse, "
 		"(f)ront, octa(g)on, (h)exagon, (i)nvtriangle, (j)rotated rectangle, "
-		"(k)ustom, (l)etter, (m)athangle, pe(n)tagon, (p)oint, (q)uoted line, (r)ectangle, "
+		"(k)ustom, (l)etter, (m)athangle, pe(n)tagon, (p)oint, s(Q)here, (q)uoted line, (r)ectangle, "
 		"(R)ounded rectangle, (s)quare, (t)riangle, (v)ector, (w)edge, (x)cross, (y)dash, or "
 		"=(geovector, i.e., great or small circle vectors) or ~(decorated line).");
 	GMT_Usage (API, -3, "If no size is specified, then the 3rd column must have sizes. "

src/psxy.c

@@ -238,7 +238,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
 	GMT_Usage (API, -3, "-(xdash), +(plus), st(a)r, (b|B)ar, (c)ircle, (d)iamond, (e)llipse, "
 		"(f)ront, octa(g)on, (h)exagon (i)nvtriangle, (j)rotated rectangle, "
 		"(k)ustom, (l)etter, (m)athangle, pe(n)tagon, c(o)lumn, (p)oint, "
-		"(q)uoted line, (r)ectangle, (R)ounded rectangle, (s)quare, (t)riangle, "
+		"s(P)here, (q)uoted line, (r)ectangle, (R)ounded rectangle, (s)quare, (t)riangle, "
 		"c(u)be, (v)ector, (w)edge, (x)cross, (y)dash, (z)dash, or "
 		"=(geovector, i.e., great or small circle vectors).");
 
@@ -276,6 +276,13 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
 	GMT_Usage (API, 2, "\n%s 3-D Cube: Give <size> as the length of all sides; append q if <size> "
 		"is a quantity in x-units.", GMT_LINE_BULLET);
 
+	GMT_Usage (API, 2, "\n%s 3-D Sphere: Give <size> as sphere diameter [Default unit is cm]. "
+		"Sphere is rendered with radial gradient shading from white at light source to fill color. Modifiers:", GMT_LINE_BULLET);
+	GMT_Usage (API, 3, "+a Set light source azimuth [0, from the right].");
+	GMT_Usage (API, 3, "+e Set light source elevation [90, perpendicular to viewing plane].");
+	GMT_Usage (API, 3, "+f Use flat/constant fill color (no gradient shading).");
+	GMT_Usage (API, 3, "+n Draw outline only (no fill). Outline color from -W.");
+
 	GMT_Usage (API, 2, "\n%s Ellipse: If not given, we read direction, major, and minor axis from columns 4-6. "
 		"If -SE rather than -Se is selected, %s will expect azimuth, and "
 		"axes [in km], and convert azimuths based on map projection. "
@@ -1048,6 +1055,10 @@ EXTERN_MSC int GMT_psxyz (void *V_API, int mode, void *args) {
 	fill_active = Ctrl->G.active;	/* Make copies because we will change the values */
 	outline_active =  Ctrl->W.active;
 	if (not_line && !outline_active && S.symbol != PSL_WEDGE && !fill_active && !get_rgb && !QR_symbol) outline_active = true;	/* If no fill nor outline for symbols then turn outline on */
+	if (S.symbol == PSL_SPHERE && !fill_active && !S.SP_no_fill) {	/* Sphere needs a default fill for 3D shading */
+		current_fill.rgb[0] = current_fill.rgb[1] = current_fill.rgb[2] = 0.5;	/* Black */
+		fill_active = true;
+	}
 
 	if (Ctrl->D.active) {
 		/* Shift the plot a bit. This is a bit frustrating, since the only way to do this
@@ -1892,6 +1903,31 @@ EXTERN_MSC int GMT_psxyz (void *V_API, int mode, void *args) {
 						gmt_plane_perspective (GMT, GMT_Z, data[i].z);
 						PSL_plotsymbol (PSL, xpos[item], data[i].y, data[i].dim, data[i].symbol);
 						break;
+					case PSL_SPHERE:	/* Case created by Claude.ai */
+						gmt_plane_perspective(GMT, GMT_Z, data[i].z);
+						if (S.SP_light_set) {	/* Calculate and output custom light position for sphere */
+							/* Simple model: azimuth controls horizontal, elevation controls vertical */
+							/* Relative to viewing direction */
+							double dazim = GMT->current.proj.z_project.view_azimuth - S.SP_light_az;
+							if (dazim > 90) dazim = 90.0;		/* Do let illum from the hidden hemisphere */
+							if (dazim < -90) dazim = -90.0;
+							double SP_lx_proj = sind(dazim);
+							double SP_ly_proj = sind(S.SP_light_el - GMT->current.proj.z_project.view_elevation);
+							PSL_command(PSL, "/SP_lx %.12g def /SP_ly %.12g def\n", SP_lx_proj, SP_ly_proj);
+						}
+						if (S.SP_flat)	/* Output flat/constant color flag for sphere */
+							PSL_command(PSL, "/SP_flat true def\n");
+						if (S.SP_no_fill)	/* Output no-fill flag for sphere */
+							PSL_command(PSL, "/SP_no_fill true def\n");
+						PSL_plotsymbol(PSL, xpos[item], data[i].y, data[i].dim, data[i].symbol);
+						/* Draw outline circle in user space using pen color */
+						if (data[i].outline) {
+							/* The 1/4 factor was obtained by trial-and-error. Couldn't find the true logic. (JL) */
+							PSL_command(PSL, "V /DeviceRGB setcolorspace matrix setmatrix %s "
+								"xc_SP_user yc_SP_user T N 0 0 radius_SP_user 0 360 arc S U\n",
+								PSL_makepen(PSL, (1.0/4.0) * data[i].p.width, data[i].p.rgb, data[i].p.style, data[i].p.offset));
+						}
+						break;
 					case PSL_ELLIPSE:
 						gmt_plane_perspective (GMT, GMT_Z, data[i].z);
 						if (data[i].flag & 2)

src/psxyz.c

@@ -1,5 +1,6 @@
 outs:
-- md5: e8300f3b8e45690deb50c78d03d3e548.dir
-  nfiles: 28
+- md5: d294c350ed35533a2b8baf4729953c9d.dir
+  nfiles: 29
   path: psxyz
   hash: md5
+  size: 1205108