1 /*--------------------------------------------------------------------
2 *
3 * Copyright (c) 1991-2021 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
4 * See LICENSE.TXT file for copying and redistribution conditions.
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; version 3 or any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU Lesser General Public License for more details.
14 *
15 * Contact info: www.generic-mapping-tools.org
16 *--------------------------------------------------------------------*/
17 /*
18 * Brief synopsis: psrose reads a file [or standard input] with azimuth and length information
19 * and draws a sector or rose diagram. Several options for plot layout are available.
20 * 2 diagrams are possible: Full circle (360) or half circle (180). In the
21 * latter case azimuths > 180 are reversed (-= 180).
22 *
23 * To be compatible with GMT, I assume radial distance to be "x"
24 * and azimuth to be "y". Hence, west = 0.0 and east = max_radius
25 * south/north is -90,90 for halfcircle and 0,360 for full circle
26 *
27 * Update Nov 2017 PW: psrose map machinery is ancient and did not support -J. For modern
28 * mode we do need -J support. To remain backwards compatible we now take these steps:
29 * a) Before parsing of arguments we check if -J was given. If not we find the radius
30 * from -S (if given, else default to 3i) and create a -JX option that we add. If
31 * normalization was set via -S then we retain a plain -S option that now only means
32 * do normalization.
33 * b) psrose then processes -JX and optionally -S. Because -JX is then parsed we must be
34 * on the look-out for any -: as these are inert. We reactivate that swap if needed.
35 * c) We then extract the radius from the -JX string, decode it, and reestablish the old
36 * -Jx1 scaling using the +/- radius as w/e/s/n.
37 * d) The old syntax is undocumented but support by backwards compatibility through GMT 6.
38 *
39 * Author: Paul Wessel
40 * Date: 1-JAN-2010
41 * Version: 6 API
42 */
43
44 #include "gmt_dev.h"
45
46 #define THIS_MODULE_CLASSIC_NAME "psrose"
47 #define THIS_MODULE_MODERN_NAME "rose"
48 #define THIS_MODULE_LIB "core"
49 #define THIS_MODULE_PURPOSE "Plot a polar histogram (rose, sector, windrose diagrams)"
50 #define THIS_MODULE_KEYS "<D{,CC(,ED(,>X},>D),>DI,ID)"
51 #define THIS_MODULE_NEEDS "JR"
52 #define THIS_MODULE_OPTIONS "-:>BJKOPRUVXYbdehiopqstwxy" GMT_OPT("c")
53
54 struct PSROSE_CTRL { /* All control options for this program (except common args) */
55 /* active is true if the option has been activated */
56 struct PSROSE_A { /* -A<sector_angle>[+r] */
57 bool active;
58 bool rose;
59 double inc;
60 } A;
61 struct PSROSE_C { /* -C<cpt> */
62 bool active;
63 char *file;
64 } C;
65 struct PSROSE_D { /* -D */
66 bool active;
67 } D;
68 struct PSROSE_E { /* -Em|[+w]<modefile> */
69 bool active;
70 bool mode;
71 bool mean;
72 char *file;
73 } E;
74 struct PSROSE_F { /* -F */
75 bool active;
76 } F;
77 struct PSROSE_G { /* -G<fill> */
78 bool active;
79 struct GMT_FILL fill;
80 } G;
81 struct PSROSE_I { /* -I */
82 bool active;
83 } I;
84 struct PSROSE_L { /* -L */
85 bool active;
86 char *w, *e, *s, *n;
87 } L;
88 struct PSROSE_M { /* -M[<size>][<modifiers>] */
89 bool active;
90 struct GMT_SYMBOL S;
91 } M;
92 struct PSROSE_N { /* -N[<kind>]+p<pen>, <kind = 0 (for now) OR -N [Deprecated to -S+a] */
93 bool active;
94 bool selected;
95 struct GMT_PEN pen;
96 } N;
97 struct PSROSE_Q { /* -Q<alpha> */
98 bool active;
99 double value;
100 } Q;
101 struct PSROSE_S { /* -S[+a] */
102 bool active;
103 bool normalize;
104 bool area_normalize;
105 double scale; /* Get this via -JX */
106 } S;
107 struct PSROSE_T { /* -T */
108 bool active;
109 } T;
110 struct PSROSE_W { /* -W[v]<pen> */
111 bool active[2];
112 struct GMT_PEN pen[2];
113 } W;
114 struct PSROSE_Z { /* -Zu|<scale> */
115 bool active;
116 unsigned int mode;
117 double scale;
118 } Z;
119 };
120
New_Ctrl(struct GMT_CTRL * GMT)121 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
122 struct PSROSE_CTRL *C = NULL;
123
124 C = gmt_M_memory (GMT, NULL, 1, struct PSROSE_CTRL);
125
126 /* Initialize values whose defaults are not 0/false/NULL */
127 gmt_init_fill (GMT, &C->G.fill, -1.0, -1.0, -1.0);
128 C->N.pen = GMT->current.setting.map_default_pen;
129 C->M.S.symbol = PSL_VECTOR;
130 C->W.pen[0] = C->W.pen[1] = GMT->current.setting.map_default_pen;
131 C->Q.value = 0.05;
132 C->Z.scale = 1.0;
133 return (C);
134 }
135
Free_Ctrl(struct GMT_CTRL * GMT,struct PSROSE_CTRL * C)136 static void Free_Ctrl (struct GMT_CTRL *GMT, struct PSROSE_CTRL *C) { /* Deallocate control structure */
137 if (!C) return;
138 gmt_M_str_free (C->E.file);
139 gmt_M_str_free (C->L.w);
140 gmt_M_str_free (C->L.e);
141 gmt_M_str_free (C->L.s);
142 gmt_M_str_free (C->L.n);
143 gmt_M_free (GMT, C);
144 }
145
psrose_critical_resultant(double alpha,int n)146 GMT_LOCAL double psrose_critical_resultant (double alpha, int n) {
147 /* Return critical resultant for given alpha and sample size.
148 * Based on Rayleigh test for uniformity as approximated by Zaar [1999]
149 * and reported by Berens [2009] in CircStat (MATLAB). Valid for
150 * n >= 10 and for first 3 decimals (gets better with n). */
151 double Rn;
152 Rn = 0.5 * sqrt ((1.0 + 4.0 * n * (1.0 + n) - pow (log (alpha) + 1.0 + 2.0 * n, 2.0))) / n;
153 return (Rn);
154 }
155
usage(struct GMTAPI_CTRL * API,int level)156 static int usage (struct GMTAPI_CTRL *API, int level) {
157 char *choice[2] = {"OFF", "ON"};
158
159 /* This displays the psrose synopsis and optionally full usage information */
160
161 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
162 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
163 GMT_Usage (API, 0, "usage: %s [<table>] [-A<sector_angle>[+r]] [%s] [-C<cpt>] [-D] [-E[m|[+w]<modefile>]] [-F] [-G<fill>] "
164 "[-I] [-JX<diameter>] %s[-L[<wlab>,<elab>,<slab>,<nlab>]] [-M[<size>][<modifiers>]] [-N<mode>[+p<pen>]] %s%s[-Q<alpha>] "
165 "[-R<r0>/<r1>/<theta0>/<theta1>] [-S[+a]] [-T] [%s] [%s] [-W[v]<pen>] [%s] [%s] [-Zu|<scale>] [%s] %s[%s] [%s] [%s] "
166 "[%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s]\n",
167 name, GMT_B_OPT, API->K_OPT, API->O_OPT, API->P_OPT, GMT_U_OPT, GMT_V_OPT, GMT_X_OPT, GMT_Y_OPT, GMT_bi_OPT, API->c_OPT,
168 GMT_di_OPT, GMT_e_OPT, GMT_h_OPT, GMT_i_OPT, GMT_o_OPT, GMT_p_OPT, GMT_qi_OPT, GMT_s_OPT, GMT_t_OPT, GMT_w_OPT, GMT_colon_OPT, GMT_PAR_OPT);
169
170 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
171
172 GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
173 GMT_Option (API, "<");
174 GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
175 GMT_Usage (API, 1, "\n-A<sector_angle>[+r]");
176 GMT_Usage (API, -2, "Set sector width in degrees for sector diagram [Default is windrose].");
177 GMT_Usage (API, 3, "+r Select rose diagram.");
178 GMT_Option (API, "B-");
179 if (gmt_M_showusage (API)) {
180 GMT_Usage (API, -2, "Note: The scale bar length is set to the radial gridline spacing. "
181 "Remember: radial is x-direction, azimuthal is y-direction.");
182 }
183 GMT_Usage (API, 1, "\n-C<cpt>");
184 GMT_Usage (API, -2, "Use CPT to assign fill to sectors based on the r-value. Requires -A (sector diagram).");
185 GMT_Usage (API, 1, "\n-D Shift sectors so that they are centered on the bin interval (e.g., first sector is centered on 0 degrees).");
186 GMT_Usage (API, 1, "\n-E[m|[+w]<modefile>]");
187 GMT_Usage (API, -2, "Plot vectors listed in the <modefile>. For calculated mean direction instead, choose -Em, with optional modifier:");
188 GMT_Usage (API, 3, "+w Write the calculated mean direction to <modefile>.");
189 GMT_Usage (API, 1, "\n-F Do not draw the scale length bar [Default plots scale in lower right corner].");
190 gmt_fill_syntax (API->GMT, 'G', NULL, "Specify color for diagram [Default is no fill].");
191 GMT_Usage (API, 1, "\n-I Inquire mode; only compute and report statistics - no plot is created.");
192 GMT_Usage (API, 1, "\n-J Use -JX<diameter> to set the plot diameter [7.5c].");
193 GMT_Option (API, "K");
194 GMT_Usage (API, 1, "\n-L[<wlab>,<elab>,<slab>,<nlab>]");
195 GMT_Usage (API, -2, "Override default labels [West,East,South,North (depending on GMT_LANGUAGE) "
196 "for full circle and 90W,90E,-,0 for half-circle]. If no argument "
197 "is given then labels will be disabled. Give - to disable an individual label.");
198 GMT_Usage (API, 1, "\n-M Specify arrow attributes. If -E is used then the attributes apply to the -E vector(s). "
199 "Otherwise, if windrose mode is selected we apply vector attributes to individual directions.");
200 gmt_vector_syntax (API->GMT, 15, 3);
201 GMT_Usage (API, -2, "Default is %gp+gblack+p1p.", VECTOR_HEAD_LENGTH);
202 GMT_Usage (API, 1, "\n-N<mode>[+p<pen>]");
203 GMT_Usage (API, -2, "Append <mode> to draw the equivalent von Mises distribution; optionally append desired pen via modifier +p [0.25p,black]. "
204 "<mode> selects which central location and scale to use:");
205 GMT_Usage (API, 3, "0: mean and standard deviation [Default].");
206 GMT_Usage (API, -2, "(Other modes may be added later).");
207 GMT_Option (API, "O,P");
208 GMT_Usage (API, 1, "\n-Q<alpha>");
209 GMT_Usage (API, -2, "Set confidence level for Rayleigh test for uniformity [0.05].");
210 GMT_Usage (API, 1, "\n-R<r0>/<r1>/<theta0>/<theta1>");
211 GMT_Usage (API, -2, "Specify the region (<r0> = 0, <r1> = max_radius). For azimuth: "
212 "Specify <theta0>/<theta1> = -90/90 or 0/180 (half-circles) or 0/360 only). "
213 "If <r0> = <r1> = 0, psrose will compute a reasonable <r1> value.");
214 GMT_Usage (API, 1, "\n-S[+a]");
215 GMT_Usage (API, -2, "Normalize data, i.e., divide all radii (or bin counts) by the maximum radius (or count). Optional modifier:");
216 GMT_Usage (API, 3, "+a Normalize rose plots for area, i.e., take sqrt(r) before plotting [no area normalization].");
217 GMT_Usage (API, 1, "\n-T Indicate that the vectors are oriented (two-headed), not directed [Default]. "
218 "This implies both <azimuth> and <azimuth> + 180 will be counted as inputs. "
219 "Ignored if -R sets a half-circle domain.");
220 GMT_Option (API, "U,V");
221 gmt_pen_syntax (API->GMT, 'W', NULL, "Set pen attributes for outline of rose [Default is no outline].", NULL, 0);
222 GMT_Message (API, GMT_TIME_NONE, "\t Use -Wv<pen> to set a different pen for the vector (requires -E) [Same as rose outline].\n");
223 GMT_Option (API, "X");
224 GMT_Usage (API, 1, "\n-Zu|<scale>");
225 GMT_Usage (API, -2, "Multiply the radii by <scale> before plotting; or use -Zu to set input radii to 1.");
226 GMT_Usage (API, 1, "\n-: Expect (azimuth,radius) input rather than (radius,azimuth) [%s].", choice[API->GMT->current.setting.io_lonlat_toggle[GMT_IN]]);
227 GMT_Option (API, "bi2,c,di,e,h,i,o,p,qi,s,t,w,.");
228
229 return (GMT_MODULE_USAGE);
230 }
231
parse(struct GMT_CTRL * GMT,struct PSROSE_CTRL * Ctrl,struct GMT_OPTION * options)232 static int parse (struct GMT_CTRL *GMT, struct PSROSE_CTRL *Ctrl, struct GMT_OPTION *options) {
233 /* This parses the options provided to psrose and sets parameters in Ctrl.
234 * Note Ctrl has already been initialized and non-zero default values set.
235 * Any GMT common options will override values set previously by other commands.
236 * It also replaces any file names specified as input or output with the data ID
237 * returned when registering these sources/destinations with the API.
238 */
239
240 int n;
241 unsigned int n_errors = 0, k;
242 double range;
243 char txt_a[GMT_LEN256] = {""}, txt_b[GMT_LEN256] = {""}, txt_c[GMT_LEN256] = {""}, txt_d[GMT_LEN256] = {""}, *c = NULL;
244 struct GMT_OPTION *opt = NULL;
245 struct GMTAPI_CTRL *API = GMT->parent;
246
247 for (opt = options; opt; opt = opt->next) { /* Process all the options given */
248
249 switch (opt->option) {
250
251 case '<': /* Input files */
252 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(opt->arg))) n_errors++;;
253 break;
254
255 /* Processes program-specific parameters */
256
257 case 'A': /* Get Sector angle in degrees -A<inc>[+r]*/
258 n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active);
259 Ctrl->A.active = true;
260 k = 0;
261 if ((c = strstr (opt->arg, "+r"))) {
262 Ctrl->A.rose = true;
263 c[0] = '\0'; /* Chop off modifier */
264 }
265 else if (strchr (opt->arg, 'r')) { /* Old syntax -A[r]<inc> */
266 Ctrl->A.rose = true;
267 if (opt->arg[0] == 'r') k = 1;
268 }
269 Ctrl->A.inc = atof (&opt->arg[k]);
270 if (c) c[0] = '+'; /* Restore modifier */
271 break;
272 case 'C':
273 n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active);
274 if (gmt_M_compat_check (GMT, 5)) { /* Need to check for deprecated -Cm|[+w]<modefile> option */
275 if (((c = strstr (opt->arg, "+w")) || (opt->arg[0] == 'm' && opt->arg[1] == '\0') || opt->arg[0] == '\0') && strstr (opt->arg, GMT_CPT_EXTENSION) == NULL) {
276 GMT_Report (API, GMT_MSG_COMPAT, "Option -C for mode-vector(s) is deprecated; use -E instead.\n");
277 Ctrl->E.active = true;
278 if ((c = strstr (opt->arg, "+w"))) { /* Wants to write out mean direction */
279 gmt_M_str_free (Ctrl->E.file);
280 if (c[2]) Ctrl->E.file = strdup (&c[2]);
281 Ctrl->E.mode = GMT_OUT;
282 c[0] = '\0'; /* Chop off temporarily */
283 }
284 if ((opt->arg[0] == 'm' && opt->arg[1] == '\0') || (API->external == 0 && opt->arg[0] == '\0'))
285 Ctrl->E.mean = true;
286 else if (Ctrl->E.mode == GMT_IN) {
287 gmt_M_str_free (Ctrl->E.file);
288 if (opt->arg[0]) Ctrl->E.file = strdup (opt->arg);
289 }
290 break;
291 }
292 }
293 Ctrl->C.active = true;
294 gmt_M_str_free (Ctrl->C.file);
295 if (opt->arg[0]) Ctrl->C.file = strdup (opt->arg);
296 break;
297 case 'D': /* Center the bins */
298 n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active);
299 Ctrl->D.active = true;
300 break;
301 case 'E': /* Read mode file and plot mean directions */
302 n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active);
303 Ctrl->E.active = true;
304 if ((c = strstr (opt->arg, "+w"))) { /* Wants to write out mean direction */
305 gmt_M_str_free (Ctrl->E.file);
306 if (c[2]) Ctrl->E.file = strdup (&c[2]);
307 Ctrl->E.mode = GMT_OUT;
308 c[0] = '\0'; /* Chop off temporarily */
309 }
310 if ((opt->arg[0] == 'm' && opt->arg[1] == '\0') || (API->external == 0 && opt->arg[0] == '\0'))
311 Ctrl->E.mean = true;
312 else if (Ctrl->E.mode == GMT_IN) {
313 gmt_M_str_free (Ctrl->E.file);
314 if (opt->arg[0]) Ctrl->E.file = strdup (opt->arg);
315 }
316 break;
317 case 'F': /* Disable scalebar plotting */
318 n_errors += gmt_M_repeated_module_option (API, Ctrl->F.active);
319 Ctrl->F.active = true;
320 break;
321 case 'G': /* Set Gray shade */
322 n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active);
323 Ctrl->G.active = true;
324 if (gmt_getfill (GMT, opt->arg, &Ctrl->G.fill)) {
325 gmt_fill_syntax (GMT, 'G', NULL, " ");
326 n_errors++;
327 }
328 break;
329 case 'I': /* Compute statistics only - no plot */
330 n_errors += gmt_M_repeated_module_option (API, Ctrl->I.active);
331 Ctrl->I.active = true;
332 break;
333 case 'L': /* Override default labeling */
334 n_errors += gmt_M_repeated_module_option (API, Ctrl->L.active);
335 Ctrl->L.active = true;
336 if (opt->arg[0]) {
337 unsigned int n_comma = 0;
338 for (k = 0; k < strlen (opt->arg); k++) if (opt->arg[k] == ',') n_comma++;
339 if (n_comma == 3) /* New, comma-separated labels */
340 n_errors += gmt_M_check_condition (GMT, sscanf (opt->arg, "%[^,],%[^,],%[^,],%s", txt_a, txt_b, txt_c, txt_d) != 4, "Option -L: Expected\n\t-L<westlabel>,<eastlabel>,<southlabel>,<northlabel>\n");
341 else /* Old slash-separated labels */
342 n_errors += gmt_M_check_condition (GMT, sscanf (opt->arg, "%[^/]/%[^/]/%[^/]/%s", txt_a, txt_b, txt_c, txt_d) != 4, "Option -L: Expected\n\t-L<westlabel>/<eastlabel>/<southlabel>/<northlabel>\n");
343 Ctrl->L.w = strdup (txt_a); Ctrl->L.e = strdup (txt_b);
344 Ctrl->L.s = strdup (txt_c); Ctrl->L.n = strdup (txt_d);
345 }
346 else { /* Turn off all 4 labels */
347 Ctrl->L.w = strdup ("-"); Ctrl->L.e = strdup ("-");
348 Ctrl->L.s = strdup ("-"); Ctrl->L.n = strdup ("-");
349 }
350 break;
351 case 'M': /* Get arrow parameters */
352 n_errors += gmt_M_repeated_module_option (API, Ctrl->M.active);
353 Ctrl->M.active = true;
354 if (gmt_M_compat_check (GMT, 4) && (strchr (opt->arg, '/') && !strchr (opt->arg, '+'))) { /* Old-style args */
355 n = sscanf (opt->arg, "%[^/]/%[^/]/%[^/]/%s", txt_a, txt_b, txt_c, txt_d);
356 if (n != 4 || gmt_getrgb (GMT, txt_d, Ctrl->M.S.v.fill.rgb)) {
357 GMT_Report (API, GMT_MSG_ERROR, "Option -M: Expected\n\t-M<tailwidth/headlength/headwidth/<color>\n");
358 n_errors++;
359 }
360 else { /* Turn the old args into new +a<angle> and pen width */
361 Ctrl->M.S.v.status = PSL_VEC_END + PSL_VEC_FILL + PSL_VEC_OUTLINE;
362 Ctrl->M.S.size_x = VECTOR_HEAD_LENGTH * GMT->session.u2u[GMT_PT][GMT_INCH]; /* 9p */
363 Ctrl->M.S.v.h_length = (float)Ctrl->M.S.size_x; /* 9p */
364 Ctrl->M.S.v.v_angle = 60.0f;
365 Ctrl->M.S.v.pen = GMT->current.setting.map_default_pen;
366 Ctrl->W.active[1] = true;
367 //Ctrl->W.pen[1].width = gmt_M_to_points (GMT, txt_a);
368 Ctrl->M.S.v.v_width = (float)gmt_M_to_inch (GMT, txt_a);
369 Ctrl->M.S.v.h_length = (float)gmt_M_to_inch (GMT, txt_b);
370 Ctrl->M.S.v.h_width = (float)gmt_M_to_inch (GMT, txt_c);
371 Ctrl->M.S.v.v_angle = (float)atand (0.5 * Ctrl->M.S.v.h_width / Ctrl->M.S.v.h_length);
372 Ctrl->M.S.v.status |= (PSL_VEC_OUTLINE + PSL_VEC_FILL);
373 Ctrl->M.S.v.status |= PSL_VEC_FILL2;
374 }
375 Ctrl->M.S.symbol = GMT_SYMBOL_VECTOR_V4;
376 }
377 else {
378 if (opt->arg[0] == '+' || opt->arg[0] == '\0') { /* No size argument (use default), just attributes */
379 n_errors += gmt_parse_vector (GMT, 'v', opt->arg, &Ctrl->M.S);
380 }
381 else { /* Size, plus possible attributes */
382 n = sscanf (opt->arg, "%[^+]%s", txt_a, txt_b); /* txt_a should be symbols size with any +<modifiers> in txt_b */
383 if (n == 1) txt_b[0] = 0; /* No modifiers present, set txt_b to empty */
384 Ctrl->M.S.size_x = gmt_M_to_inch (GMT, txt_a); /* Length of vector */
385 n_errors += gmt_parse_vector (GMT, 'v', txt_b, &Ctrl->M.S);
386 }
387 Ctrl->M.S.v.status |= PSL_VEC_OUTLINE;
388 }
389 break;
390 case 'N': /* Make sectors area be proportional to frequency instead of radius [DEPRECATED in 6.2] */
391 if (opt->arg[0] == '\0') /* Only plain -N is accepted to be backwards compatible */
392 Ctrl->S.area_normalize = true;
393 else { /* Modern -N to draw VPDF with a pen */
394 Ctrl->N.active = true;
395 switch (opt->arg[0]) { /* See which distribution to draw */
396 case '0': break; /* Only allowed mode for now */
397 default:
398 GMT_Report (API, GMT_MSG_ERROR, "Option -N: mode %c unrecognized.\n", opt->arg[0]);
399 n_errors++;
400 break;
401 }
402 Ctrl->N.selected = true;
403 if ((c = strstr (opt->arg, "+p")) != NULL) {
404 if (gmt_getpen (GMT, &c[2], &Ctrl->N.pen)) {
405 gmt_pen_syntax (GMT, 'N', NULL, " ", NULL, 0);
406 n_errors++;
407 }
408 }
409 }
410 break;
411 case 'Q': /* Set critical value [0.05] */
412 n_errors += gmt_M_repeated_module_option (API, Ctrl->Q.active);
413 Ctrl->Q.active = true;
414 if (opt->arg[0]) Ctrl->Q.value = atof (opt->arg);
415 break;
416 case 'S': /* Normalization */
417 n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active);
418 Ctrl->S.active = true;
419 Ctrl->S.normalize = true;
420 if (strstr (opt->arg, "+a"))
421 Ctrl->S.area_normalize = true;
422 break;
423 case 'T': /* Oriented instead of directed data */
424 n_errors += gmt_M_repeated_module_option (API, Ctrl->T.active);
425 Ctrl->T.active = true;
426 break;
427 case 'W': /* Get pen width for outline */
428 n = (opt->arg[0] == 'v') ? 1 : 0;
429 n_errors += gmt_M_repeated_module_option (API, Ctrl->W.active[n]);
430 Ctrl->W.active[n] = true;
431 if (gmt_getpen (GMT, &opt->arg[n], &Ctrl->W.pen[n])) {
432 gmt_pen_syntax (GMT, 'W', NULL, " ", NULL, 0);
433 n_errors++;
434 }
435 break;
436 case 'Z': /* Scale radii before using data */
437 n_errors += gmt_M_repeated_module_option (API, Ctrl->Z.active);
438 Ctrl->Z.active = true;
439 if (opt->arg[0] == 'u')
440 Ctrl->Z.mode = 1;
441 else
442 Ctrl->Z.scale = atof (opt->arg);
443 break;
444
445 default: /* Report bad options */
446 n_errors += gmt_default_error (GMT, opt->option);
447 break;
448 }
449 }
450
451 gmt_consider_current_cpt (API, &Ctrl->C.active, &(Ctrl->C.file));
452
453 /* Check that the options selected are mutually consistent */
454
455 range = GMT->common.R.wesn[YHI] - GMT->common.R.wesn[YLO];
456 if (doubleAlmostEqual (range, 180.0) && Ctrl->T.active) {
457 GMT_Report (API, GMT_MSG_WARNING, "-T only needed for 0-360 range data (ignored)");
458 Ctrl->T.active = false;
459 }
460 n_errors += gmt_M_check_condition (GMT, Ctrl->E.active && Ctrl->E.file && Ctrl->E.mode == GMT_IN && gmt_access (GMT, Ctrl->E.file, R_OK),
461 "Option -E: Cannot read file %s!\n", Ctrl->E.file);
462 n_errors += gmt_M_check_condition (GMT, gmt_M_is_zero (Ctrl->Z.scale), "Option -Z: factor must be nonzero\n");
463 n_errors += gmt_M_check_condition (GMT, Ctrl->A.inc < 0.0, "Option -A: sector width must be positive\n");
464 n_errors += gmt_M_check_condition (GMT, Ctrl->Q.value <= 0.0 || Ctrl->Q.value >= 1.0, "Option -Q: confidence level must be in 0-1 range\n");
465 n_errors += gmt_M_check_condition (GMT, Ctrl->C.active && Ctrl->G.active, "Option -C: Cannot give both -C and -G\n");
466 n_errors += gmt_M_check_condition (GMT, Ctrl->C.active && Ctrl->A.rose, "Option -C: Cannot be used with -A+r\n");
467 n_errors += gmt_M_check_condition (GMT, Ctrl->C.active && !Ctrl->A.active, "Option -C: Requires -A\n");
468 n_errors += gmt_M_check_condition (GMT, Ctrl->N.active && !Ctrl->A.active, "Option -N: Requires -A\n");
469 n_errors += gmt_M_check_condition (GMT, Ctrl->N.active && !doubleAlmostEqual (range, 360.0), "Option -N: Requires the full circle in -R\n");
470 if (GMT->common.J.active) { /* Impose our conditions on -JX */
471 n_errors += gmt_M_check_condition (GMT, GMT->common.J.string[0] != 'X', "Option -J: Must specify -JX<diameter>\n");
472 n_errors += gmt_M_check_condition (GMT, strchr (GMT->common.J.string, '/'), "Option -J: Must specify -JX<diameter>\n");
473 }
474 if (!Ctrl->I.active) {
475 n_errors += gmt_M_check_condition (GMT, !GMT->common.J.active, "Must specify -JX option\n");
476 n_errors += gmt_M_check_condition (GMT, !GMT->common.R.active[RSET], "Must specify -R option\n");
477 n_errors += gmt_M_check_condition (GMT, !((GMT->common.R.wesn[YLO] == -90.0 && GMT->common.R.wesn[YHI] == 90.0) \
478 || (GMT->common.R.wesn[YLO] == 0.0 && GMT->common.R.wesn[YHI] == 180.0)
479 || (GMT->common.R.wesn[YLO] == 0.0 && GMT->common.R.wesn[YHI] == 360.0)),
480 "Option -R: theta0/theta1 must be either -90/90, 0/180 or 0/360\n");
481 n_errors += gmt_M_check_condition (GMT, GMT->common.R.wesn[XLO] != 0.0, "Option -R: r0/r1 must have r0 == 0 and r1 must be positive\n");
482 }
483 n_errors += gmt_check_binary_io (GMT, 2);
484
485 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
486 }
487
488 #define bailout(code) {gmt_M_free_options (mode); return (code);}
489 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
490
GMT_psrose(void * V_API,int mode,void * args)491 EXTERN_MSC int GMT_psrose (void *V_API, int mode, void *args) {
492 bool do_fill = false, automatic = false, sector_plot = false, windrose = true, do_labels = true;
493 unsigned int n_bins, n_modes = 0, form, n_in, half_only = 0, bin, save;
494 int error = 0, k, n_annot, sbin, significant, index;
495
496 uint64_t n = 0, i;
497
498 size_t n_alloc = GMT_CHUNK;
499
500 char text[GMT_BUFSIZ] = {""}, format[GMT_BUFSIZ] = {""};
501
502 double max = 0.0, radius, az, x_origin, y_origin, tmp, one_or_two = 1.0, s, c, f;
503 double angle1, angle2, angle, x, y, mean_theta, mean_radius, xr = 0.0, yr = 0.0, area = 0.0;
504 double x1, x2, y1, y2, total = 0.0, total_arc, off, max_radius, az_offset, start_angle;
505 double asize, lsize, this_az, half_bin_width, diameter, wesn[4], mean_vector, mean_resultant;
506 double *xx = NULL, *yy = NULL, *in = NULL, *sum = NULL, *azimuth = NULL, critical_resultant;
507 double *length = NULL, *mode_direction = NULL, *mode_length = NULL, dim[PSL_MAX_DIMS], rgb[4];
508
509 struct PSROSE_CTRL *Ctrl = NULL;
510 struct GMT_FILL *F = NULL;
511 struct GMT_DATASET *Cin = NULL;
512 struct GMT_DATATABLE *T = NULL;
513 struct GMT_PALETTE *P = NULL;
514 struct GMT_RECORD *In = NULL;
515 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL; /* General GMT internal parameters */
516 struct GMT_OPTION *options = NULL;
517 struct PSL_CTRL *PSL = NULL; /* General PSL internal parameters */
518 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
519
520 /*----------------------- Standard module initialization and parsing ----------------------*/
521
522 if (API == NULL) return (GMT_NOT_A_SESSION);
523 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
524 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
525
526 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
527
528 /* Parse the command-line arguments; return if errors are encountered */
529
530 if ((GMT = gmt_init_module (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_KEYS, THIS_MODULE_NEEDS, NULL, &options, &GMT_cpy)) == NULL) bailout (API->error); /* Save current state */
531 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
532 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
533 if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);
534
535 /*---------------------------- This is the psrose main code ----------------------------*/
536
537 GMT_Report (API, GMT_MSG_INFORMATION, "Processing input table data\n");
538 gmt_M_memset (dim, PSL_MAX_DIMS, double);
539
540 max_radius = GMT->common.R.wesn[XHI];
541 if (doubleAlmostEqual (GMT->common.R.wesn[YLO], -90.0))
542 half_only = 1;
543 else if (doubleAlmostEqual (GMT->common.R.wesn[YHI], 180.0))
544 half_only = 2;
545 if (Ctrl->A.rose) windrose = false;
546 sector_plot = (Ctrl->A.inc > 0.0);
547 if (sector_plot) windrose = false; /* Draw rose diagram instead of sector diagram */
548 if (!Ctrl->S.normalize) Ctrl->S.area_normalize = false; /* Only do this if data is normalized for length also */
549 if (!Ctrl->I.active && !GMT->common.R.active[RSET]) automatic = true;
550 if (Ctrl->T.active) one_or_two = 2.0;
551 half_bin_width = Ctrl->D.active * Ctrl->A.inc * 0.5;
552 if (half_only == 1) {
553 total_arc = 180.0;
554 az_offset = 90.0;
555 start_angle = 90.0;
556 }
557 else if (half_only == 2) {
558 total_arc = 180.0;
559 az_offset = 0.0;
560 start_angle = 180.0;
561 }
562 else {
563 total_arc = 360.0;
564 az_offset = 0.0;
565 start_angle = 90.0;
566 }
567 n_bins = (Ctrl->A.inc <= 0.0) ? 1U : urint (total_arc / Ctrl->A.inc);
568
569 /* Read data and do some stats */
570
571 if (Ctrl->C.active && (P = GMT_Read_Data (API, GMT_IS_PALETTE, GMT_IS_FILE, GMT_IS_NONE, GMT_READ_NORMAL, NULL, Ctrl->C.file, NULL)) == NULL) {
572 Return (GMT_DATA_READ_ERROR);
573 }
574
575 n = 0;
576 n_in = (GMT->common.i.select && GMT->common.i.n_cols == 1) ? 1 : 2;
577
578 if ((error = GMT_Set_Columns (API, GMT_IN, n_in, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) {
579 Return (error);
580 }
581 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_IN, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Register data input */
582 Return (API->error);
583 }
584 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_IN, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data input and sets access mode */
585 Return (API->error);
586 }
587
588 /* Allocate arrays */
589
590 sum = gmt_M_memory (GMT, NULL, n_bins, double);
591 xx = gmt_M_memory (GMT, NULL, n_bins+2, double);
592 yy = gmt_M_memory (GMT, NULL, n_bins+2, double);
593 azimuth = gmt_M_memory (GMT, NULL, n_alloc, double);
594 length = gmt_M_memory (GMT, NULL, n_alloc, double);
595
596 /* Because of -JX being parsed already, any -: will have no effect. For backwards compatibility we
597 * check if -: was given and turn that on again here before reading */
598 if (GMT->common.colon.active)
599 GMT->current.setting.io_lonlat_toggle[GMT_IN] = true;
600
601 do { /* Keep returning records until we reach EOF */
602 if ((In = GMT_Get_Record (API, GMT_READ_DATA, NULL)) == NULL) { /* Read next record, get NULL if special case */
603 if (gmt_M_rec_is_error (GMT)) { /* Bail if there are any read errors */
604 gmt_M_free (GMT, length); gmt_M_free (GMT, xx); gmt_M_free (GMT, sum);
605 gmt_M_free (GMT, azimuth); gmt_M_free (GMT, yy);
606 Return (GMT_RUNTIME_ERROR);
607 }
608 else if (gmt_M_rec_is_eof (GMT)) /* Reached end of file */
609 break;
610 continue; /* Go back and read the next record */
611 }
612
613 if (In->data == NULL) {
614 gmt_quit_bad_record (API, In);
615 Return (API->error);
616 }
617
618 /* Data record to process */
619 in = In->data; /* Only need to process numerical part here */
620
621 if (n_in == 2) { /* Read azimuth and length */
622 length[n] = in[GMT_X];
623 azimuth[n] = in[GMT_Y];
624 if (Ctrl->Z.active) {
625 if (Ctrl->Z.mode) length[n] = 1.0;
626 else if (Ctrl->Z.scale != 1.0) length[n] *= Ctrl->Z.scale;
627 }
628 }
629 else { /* Only read azimuth; set length = weight = 1 */
630 length[n] = 1.0;
631 azimuth[n] = in[GMT_X];
632 }
633
634 /* Make sure azimuth is in 0 <= az < 360 range */
635
636 while (azimuth[n] < 0.0) azimuth[n] += 360.0;
637 while (azimuth[n] >= 360.0) azimuth[n] -= 360.0;
638
639 if (half_only == 1) { /* Flip azimuths about E-W line i.e. -90 < az <= 90 */
640 if (azimuth[n] > 90.0 && azimuth[n] <= 270.0) azimuth[n] -= 180.0;
641 if (azimuth[n] > 270.0) azimuth[n] -= 360.0;
642 }
643 else if (half_only == 2) { /* Flip azimuths about N-S line i.e. 0 < az <= 180 */
644 if (azimuth[n] > 180.0) azimuth[n] -= 180.0;
645 }
646 else if (Ctrl->T.active) {
647 azimuth[n] = 0.5 * fmod (2.0 * azimuth[n], 360.0);
648 }
649
650 /* Double angle to find mean azimuth */
651
652 sincosd (one_or_two * azimuth[n], &s, &c);
653 xr += length[n] * c;
654 yr += length[n] * s;
655
656 total += length[n];
657
658 n++;
659 if (n == n_alloc) { /* Get more memory */
660 n_alloc <<= 1;
661 azimuth = gmt_M_memory (GMT, azimuth, n_alloc, double);
662 length = gmt_M_memory (GMT, length, n_alloc, double);
663 }
664 } while (true);
665
666 if (Ctrl->A.inc > 0.0) { /* Sum up sector diagram info */
667 for (i = 0; i < n; i++) {
668 if (Ctrl->D.active) { /* Center bin by removing half bin width here */
669 this_az = azimuth[i] - half_bin_width;
670 if (!half_only && this_az < 0.0) this_az += 360.0;
671 if (half_only == 1 && this_az < -90.0) this_az += 180.0;
672 if (half_only == 2 && this_az < 0.0) this_az += 180.0;
673 }
674 else
675 this_az = azimuth[i];
676 sbin = irint (floor ((this_az + az_offset) / Ctrl->A.inc));
677 assert (sbin >= 0);
678 bin = sbin;
679 if (bin == n_bins) {
680 bin = 0;
681 }
682 assert (bin < n_bins);
683 sum[bin] += length[i];
684 if (Ctrl->T.active) { /* Also count the other end of the orientation */
685 this_az += 180.0; if (this_az >= 360.0) this_az -= 360.0;
686 bin = irint (floor ((this_az + az_offset) / Ctrl->A.inc));
687 sum[bin] += length[i];
688 }
689 }
690 }
691
692 mean_theta = d_atan2d (yr, xr) / one_or_two;
693 if (mean_theta < 0.0) mean_theta += 360.0;
694 mean_vector = hypot (xr, yr) / n;
695 mean_resultant = mean_radius = hypot (xr, yr) / total;
696 critical_resultant = psrose_critical_resultant (Ctrl->Q.value, (int)n);
697 significant = (mean_resultant > critical_resultant);
698 if (!Ctrl->S.normalize) mean_radius *= max_radius;
699
700 if (Ctrl->A.inc > 0.0) { /* Find max of the bins */
701 for (bin = 0; bin < n_bins; bin++) if (sum[bin] > max) max = sum[bin];
702 if (Ctrl->S.normalize) for (bin = 0; bin < n_bins; bin++) sum[bin] /= max;
703 if (Ctrl->S.area_normalize) for (bin = 0; bin < n_bins; bin++) sum[bin] = sqrt (sum[bin]);
704 for (bin = 0; bin < n_bins; bin++) area += sum[bin];
705 }
706 else { /* Find max length of individual vectors */
707 for (i = 0; i < n; i++) if (length[i] > max) max = length[i];
708 if (Ctrl->S.normalize)
709 for (i = 0; i < n; i++) length[i] /= max;
710 }
711
712 if (Ctrl->I.active || gmt_M_is_verbose (GMT, GMT_MSG_INFORMATION)) {
713 char *kind[2] = {"r", "bin sum"};
714 sprintf (format, "Info for data: n = %% " PRIu64 " mean az = %s mean r = %s mean resultant length = %s max %s = %s scaled mean r = %s linear length sum = %s sign@%%.2f = %%d\n",
715 GMT->current.setting.format_float_out, GMT->current.setting.format_float_out, GMT->current.setting.format_float_out, kind[Ctrl->A.active],
716 GMT->current.setting.format_float_out, GMT->current.setting.format_float_out, GMT->current.setting.format_float_out);
717 GMT_Report (API, GMT_MSG_INFORMATION, format, n, mean_theta, mean_vector, mean_resultant, max, mean_radius, total, Ctrl->Q.value, significant);
718 if (Ctrl->I.active) { /* That was all we needed to do, wrap up */
719 double out[7];
720 unsigned int col_type[2];
721 struct GMT_RECORD *Rec = gmt_new_record (GMT, out, NULL);
722 gmt_M_memcpy (col_type, GMT->current.io.col_type[GMT_OUT], 2U, unsigned int); /* Save first 2 current output col types */
723 gmt_set_column_type (GMT, GMT_OUT, GMT_X, GMT_IS_FLOAT);
724 gmt_set_column_type (GMT, GMT_OUT, GMT_Y, GMT_IS_FLOAT);
725 if ((error = GMT_Set_Columns (API, GMT_OUT, 7U, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) {
726 gmt_M_free (GMT, sum);
727 gmt_M_free (GMT, xx);
728 gmt_M_free (GMT, yy);
729 gmt_M_free (GMT, azimuth);
730 gmt_M_free (GMT, length);
731 gmt_M_free (GMT, Rec);
732 Return (error);
733 }
734 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_NONE, GMT_OUT, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Establishes data output */
735 gmt_M_free (GMT, sum);
736 gmt_M_free (GMT, xx);
737 gmt_M_free (GMT, yy);
738 gmt_M_free (GMT, azimuth);
739 gmt_M_free (GMT, length);
740 gmt_M_free (GMT, Rec);
741 Return (API->error);
742 }
743 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_ON) != GMT_NOERROR) {
744 gmt_M_free (GMT, sum);
745 gmt_M_free (GMT, xx);
746 gmt_M_free (GMT, yy);
747 gmt_M_free (GMT, azimuth);
748 gmt_M_free (GMT, length);
749 gmt_M_free (GMT, Rec);
750 Return (API->error); /* Enables data output and sets access mode */
751 }
752 if (GMT_Set_Geometry (API, GMT_OUT, GMT_IS_NONE) != GMT_NOERROR) { /* Sets output geometry */
753 gmt_M_free (GMT, sum);
754 gmt_M_free (GMT, xx);
755 gmt_M_free (GMT, yy);
756 gmt_M_free (GMT, azimuth);
757 gmt_M_free (GMT, length);
758 gmt_M_free (GMT, Rec);
759 Return (API->error);
760 }
761 sprintf (format, "n\tmean_az\tmean_r\tmean_resultant_length\tmax\tscaled_mean_r\tlinear_length_sum");
762 out[0] = (double)n; out[1] = mean_theta; out[2] = mean_vector; out[3] = mean_resultant;
763 out[4] = max; out[5] = mean_radius; out[6] = total;
764 GMT_Put_Record (API, GMT_WRITE_TABLE_HEADER, format); /* Write this to output if -ho */
765 GMT_Put_Record (API, GMT_WRITE_DATA, Rec);
766 gmt_M_free (GMT, Rec);
767 gmt_M_free (GMT, sum);
768 gmt_M_free (GMT, xx);
769 gmt_M_free (GMT, yy);
770 gmt_M_free (GMT, azimuth);
771 gmt_M_free (GMT, length);
772 if (GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) { /* Disables further data output */
773 Return (API->error);
774 }
775 gmt_M_memcpy (GMT->current.io.col_type[GMT_OUT], col_type, 2U, unsigned int); /* Restore 2 current output col types */
776 Return (GMT_NOERROR);
777 }
778 }
779
780 if (automatic) {
781 if (gmt_M_is_zero (GMT->current.map.frame.axis[GMT_X].item[GMT_ANNOT_UPPER].interval)) {
782 tmp = pow (10.0, floor (d_log10 (GMT, max)));
783 if ((max / tmp) < 3.0) tmp *= 0.5;
784 }
785 else
786 tmp = GMT->current.map.frame.axis[GMT_X].item[GMT_ANNOT_UPPER].interval;
787 max_radius = ceil (max / tmp) * tmp;
788 }
789
790 /* Ready to plot. So set up GMT projections (not used by psrose), we set region to actual plot width and scale to 1 */
791
792 Ctrl->S.scale = 0.5 * gmt_M_to_inch (GMT, &GMT->common.J.string[1]); /* Get radius from full width */
793 GMT->common.J.active = false; /* Reset projection machinery */
794 gmt_parse_common_options (GMT, "J", 'J', "x1i");
795 GMT->common.R.active[RSET] = GMT->common.J.active = true;
796 wesn[XLO] = wesn[YLO] = -Ctrl->S.scale; wesn[XHI] = wesn[YHI] = Ctrl->S.scale;
797 if (gmt_map_setup (GMT, wesn)) {
798 gmt_M_free (GMT, length); gmt_M_free (GMT, xx); gmt_M_free (GMT, sum);
799 gmt_M_free (GMT, azimuth); gmt_M_free (GMT, yy);
800 Return (GMT_PROJECTION_ERROR);
801 }
802
803 if (GMT->current.map.frame.paint[GMT_Z]) { /* Until psrose uses a polar projection we must bypass the basemap fill and do it ourself here */
804 GMT->current.map.frame.paint[GMT_Z] = false; /* Turn off so gmt_plotinit won't fill */
805 do_fill = true;
806 }
807 if ((PSL = gmt_plotinit (GMT, options)) == NULL) Return (GMT_RUNTIME_ERROR);
808
809 x_origin = Ctrl->S.scale; y_origin = ((half_only) ? 0.0 : Ctrl->S.scale);
810 diameter = 2.0 * Ctrl->S.scale;
811 PSL_setorigin (PSL, x_origin, y_origin, 0.0, PSL_FWD);
812 gmt_plane_perspective (GMT, GMT->current.proj.z_project.view_plane, GMT->current.proj.z_level);
813 gmt_plotcanvas (GMT); /* Fill canvas if requested */
814
815 if (!Ctrl->S.normalize) Ctrl->S.scale /= max_radius;
816
817 /* Must redo any automatically set intervals via -Bag is junk, so the wesn passed to gmt_map_setup had no angles, we must rerun them */
818 /* The factor of 2 is a bit ad-hoc but yields graticules that are more square than otherwise */
819 f = (GMT->current.map.frame.axis[GMT_X].item[GMT_ANNOT_UPPER].generated) ? 1.0 : 2.0;
820 GMT->common.R.wesn[XLO] = -f*max_radius; GMT->common.R.wesn[XHI] = f*max_radius;
821 GMT->common.R.wesn[YLO] = 0.0; GMT->common.R.wesn[YHI] = 360.0;
822 for (k = GMT_X; k < GMT_Z; k++) { /* If an interval was generated, rest to 0 */
823 if (GMT->current.map.frame.axis[k].item[GMT_ANNOT_UPPER].generated) GMT->current.map.frame.axis[k].item[GMT_ANNOT_UPPER].interval = 0.0;
824 if (GMT->current.map.frame.axis[k].item[GMT_GRID_UPPER].generated) GMT->current.map.frame.axis[k].item[GMT_GRID_UPPER].interval = 0.0;
825 }
826 save = gmt_get_column_type (GMT, GMT_IN, GMT_Y);
827 GMT->current.io.col_type[GMT_IN][GMT_Y] = GMT_IS_GEO; /* Let y be geographic to get division fo 90 etc */
828 /* Update, if generated previously */
829 gmt_auto_frame_interval (GMT, GMT_X, GMT_ANNOT_UPPER);
830 gmt_auto_frame_interval (GMT, GMT_Y, GMT_ANNOT_UPPER);
831 /* Reset to what it was, i.e. Cartesian square box */
832 gmt_M_memcpy (GMT->common.R.wesn, wesn, 4U, double);
833 gmt_set_column_type (GMT, GMT_IN, GMT_Y, save); /* Reset */
834
835 if (GMT->current.map.frame.draw && !GMT->current.map.frame.no_frame && gmt_M_is_zero (GMT->current.map.frame.axis[GMT_Y].item[GMT_ANNOT_UPPER].interval) && gmt_M_is_zero (GMT->current.map.frame.axis[GMT_Y].item[GMT_GRID_UPPER].interval)) do_labels = false;
836
837 if (do_fill) { /* Until psrose uses a polar projection we must bypass the basemap fill and do it ourself here */
838 double dim = 2.0 * Ctrl->S.scale;
839 GMT->current.map.frame.paint[GMT_Z] = true; /* Restore original setting */
840 if (half_only) { /* Clip the bottom half of the circle */
841 double xc[4], yc[4];
842 xc[0] = xc[3] = -Ctrl->S.scale; xc[1] = xc[2] = Ctrl->S.scale;
843 yc[0] = yc[1] = 0.0; yc[2] = yc[3] = Ctrl->S.scale;
844 PSL_beginclipping (PSL, xc, yc, 4, GMT->session.no_rgb, 3);
845 }
846 gmt_setfill (GMT, &GMT->current.map.frame.fill[GMT_Z], 0);
847 PSL_plotsymbol (PSL, 0.0, 0.0, &dim, PSL_CIRCLE);
848 if (half_only) PSL_endclipping (PSL, 1); /* Reduce polygon clipping by one level */
849 }
850 if (GMT->common.B.active[0] && !GMT->current.map.frame.no_frame ) { /* Draw frame */
851 int n_alpha, n_radii;
852
853 /* Lay down gridlines before histogram */
854 gmt_setpen (GMT, &GMT->current.setting.map_grid_pen[GMT_PRIMARY]);
855 off = max_radius * Ctrl->S.scale;
856 n_alpha = (GMT->current.map.frame.axis[GMT_Y].item[GMT_GRID_UPPER].interval > 0.0) ? irint (total_arc / GMT->current.map.frame.axis[GMT_Y].item[GMT_GRID_UPPER].interval) : -1;
857 for (k = 0; k <= n_alpha; k++) {
858 angle = k * GMT->current.map.frame.axis[GMT_Y].item[GMT_GRID_UPPER].interval;
859 sincosd (angle, &s, &c);
860 x = max_radius * Ctrl->S.scale * c;
861 y = max_radius * Ctrl->S.scale * s;
862 PSL_plotsegment (PSL, 0.0, 0.0, x, y);
863 }
864
865 if (GMT->current.map.frame.axis[GMT_X].item[GMT_GRID_UPPER].interval > 0.0) {
866 n_radii = urint (max_radius / GMT->current.map.frame.axis[GMT_X].item[GMT_GRID_UPPER].interval);
867 for (k = 1; k <= n_radii; k++)
868 PSL_plotarc (PSL, 0.0, 0.0, k * GMT->current.map.frame.axis[GMT_X].item[GMT_GRID_UPPER].interval * Ctrl->S.scale, 0.0, total_arc, PSL_MOVE|PSL_STROKE);
869 }
870 }
871
872 gmt_setpen (GMT, &Ctrl->W.pen[0]);
873 if (Ctrl->M.S.v.status & PSL_VEC_OUTLINE2) { /* Gave specific head outline pen */
874 PSL_defpen (GMT->PSL, "PSL_vecheadpen", Ctrl->M.S.v.pen.width, Ctrl->M.S.v.pen.style, Ctrl->M.S.v.pen.offset, Ctrl->M.S.v.pen.rgb);
875 dim[PSL_VEC_HEAD_PENWIDTH] = Ctrl->M.S.v.pen.width;
876 }
877 else if (Ctrl->M.active) {
878 PSL_defpen (GMT->PSL, "PSL_vecheadpen", 0.5 * Ctrl->W.pen[1].width, Ctrl->W.pen[1].style, Ctrl->W.pen[1].offset, Ctrl->W.pen[1].rgb);
879 dim[PSL_VEC_HEAD_PENWIDTH] = 0.5 * Ctrl->W.pen[1].width;
880 }
881 if (windrose) { /* Here we draw individual vectors */
882 if (Ctrl->M.active) { /* Initialize vector head settings */
883 gmt_init_vector_param (GMT, &Ctrl->M.S, false, false, NULL, false, NULL);
884 if (Ctrl->M.S.symbol == PSL_VECTOR) Ctrl->M.S.v.v_width = (float)(Ctrl->W.pen[1].width * GMT->session.u2u[GMT_PT][GMT_INCH]);
885 dim[PSL_VEC_HEAD_SHAPE] = Ctrl->M.S.v.v_shape;
886 dim[PSL_VEC_STATUS] = (double)Ctrl->M.S.v.status;
887 dim[PSL_VEC_HEAD_TYPE_BEGIN] = (double)Ctrl->M.S.v.v_kind[0];
888 dim[PSL_VEC_HEAD_TYPE_END] = (double)Ctrl->M.S.v.v_kind[1];
889 if (Ctrl->M.S.v.status & PSL_VEC_OUTLINE) gmt_setpen (GMT, &Ctrl->W.pen[1]);
890 if (Ctrl->M.S.v.status & PSL_VEC_FILL2) gmt_setfill (GMT, &Ctrl->M.S.v.fill, 1); /* Use fill structure */
891 }
892 for (i = 0; i < n; i++) {
893 sincosd (start_angle - azimuth[i], &s, &c);
894 radius = length[i] * Ctrl->S.scale;
895 if (Ctrl->M.active) { /* Set end point of vector */
896 dim[PSL_VEC_XTIP] = radius * c;
897 dim[PSL_VEC_YTIP] = radius * s;
898 f = (radius < Ctrl->M.S.v.v_norm) ? radius / Ctrl->M.S.v.v_norm : 1.0;
899 dim[PSL_VEC_TAIL_WIDTH] = f * Ctrl->M.S.v.v_width;
900 dim[PSL_VEC_HEAD_LENGTH] = f * Ctrl->M.S.v.h_length;
901 dim[PSL_VEC_HEAD_WIDTH] = f * Ctrl->M.S.v.h_width;
902 }
903 if (Ctrl->T.active) {
904 if (Ctrl->M.active) { /* Draw two-headed vectors */
905 if (Ctrl->M.S.symbol == GMT_SYMBOL_VECTOR_V4) {
906 int v4_outline = Ctrl->W.active[1];
907 double *this_rgb = NULL;
908 if (Ctrl->M.S.v.status & PSL_VEC_FILL2)
909 this_rgb = Ctrl->M.S.v.fill.rgb;
910 else
911 this_rgb = GMT->session.no_rgb;
912 if (v4_outline) gmt_setpen (GMT, &Ctrl->W.pen[1]);
913 v4_outline += 8; /* Double-headed */
914 dim[PSL_VEC_HEAD_SHAPE] = GMT->current.setting.map_vector_shape;
915 psl_vector_v4 (PSL, -radius * c, -radius * s, dim, this_rgb, v4_outline);
916 }
917 else
918 PSL_plotsymbol (PSL, -radius * c, -radius * s, dim, PSL_VECTOR);
919 }
920 else
921 PSL_plotsegment (PSL, -radius * c, -radius * s, radius * c, radius * s);
922 }
923 else {
924 if (Ctrl->M.active) { /* Draw one-headed vectors */
925 if (Ctrl->M.S.symbol == GMT_SYMBOL_VECTOR_V4) {
926 int v4_outline = Ctrl->W.active[1];
927 double *this_rgb = NULL;
928 if (Ctrl->M.S.v.status & PSL_VEC_FILL2)
929 this_rgb = Ctrl->M.S.v.fill.rgb;
930 else
931 this_rgb = GMT->session.no_rgb;
932 if (v4_outline) gmt_setpen (GMT, &Ctrl->W.pen[1]);
933 dim[PSL_VEC_HEAD_SHAPE] = GMT->current.setting.map_vector_shape;
934 psl_vector_v4 (PSL, 0.0, 0.0, dim, this_rgb, v4_outline);
935 }
936 else
937 PSL_plotsymbol (PSL, 0.0, 0.0, dim, PSL_VECTOR);
938 }
939 else
940 PSL_plotsegment (PSL, 0.0, 0.0, radius * c, radius * s);
941 }
942 }
943 }
944
945 if (sector_plot && !Ctrl->A.rose && (Ctrl->C.active || Ctrl->G.active)) { /* Draw pie slices for sector plot if fill is requested */
946
947 if (Ctrl->G.active)
948 gmt_setfill (GMT, &(Ctrl->G.fill), 0);
949 dim[7] = 0;
950 if (Ctrl->G.active) dim[7] = 1;
951 if (Ctrl->W.active[0]) dim[7] += 2;
952 for (bin = 0; bin < n_bins; bin++) {
953 if (Ctrl->C.active) {
954 index = gmt_get_rgb_from_z (GMT, P, sum[bin] * Ctrl->S.scale, rgb);
955 F = gmt_M_get_cptslice_pattern (P, index);
956 if (F) /* Pattern */
957 gmt_setfill (GMT, F, 0);
958 else
959 PSL_setfill (PSL, rgb, 0);
960 }
961 az = bin * Ctrl->A.inc - az_offset + half_bin_width;
962 dim[1] = (start_angle - az - Ctrl->A.inc);
963 dim[2] = dim[1] + Ctrl->A.inc;
964 dim[0] = sum[bin] * Ctrl->S.scale;
965 PSL_plotsymbol (PSL, 0.0, 0.0, dim, PSL_WEDGE);
966 }
967 }
968 else if (Ctrl->A.rose) { /* Draw rose diagram */
969
970 for (i = bin = 0; bin < n_bins; bin++, i++) {
971 az = (bin - 0.5) * Ctrl->A.inc - az_offset + half_bin_width;
972 sincosd (start_angle - az - Ctrl->A.inc, &s, &c);
973 xx[i] = Ctrl->S.scale * sum[bin] * c;
974 yy[i] = Ctrl->S.scale * sum[bin] * s;
975 }
976 if (half_only) {
977 xx[i] = Ctrl->S.scale * 0.5 * (sum[0] + sum[n_bins-1]);
978 yy[i++] = 0.0;
979 xx[i] = -xx[i-1];
980 yy[i++] = 0.0;
981 }
982 PSL_setfill (PSL, Ctrl->G.fill.rgb, Ctrl->W.active[0]);
983 PSL_plotpolygon (PSL, xx, yy, (int)i);
984 }
985
986 if (sector_plot && Ctrl->W.active[0] && !Ctrl->A.rose) { /* Draw a line outlining the pie slices */
987 angle1 = ((half_only) ? 180.0 : 90.0) - half_bin_width;
988 angle2 = ((half_only) ? 0.0 : 90.0) - half_bin_width;
989 sincosd (angle1, &s, &c);
990 x1 = (sum[0] * Ctrl->S.scale) * c;
991 y1 = (sum[0] * Ctrl->S.scale) * s;
992 sincosd (angle2, &s, &c);
993 x2 = (sum[n_bins-1] * Ctrl->S.scale) * c;
994 y2 = (sum[n_bins-1] * Ctrl->S.scale) * s;
995 PSL_plotpoint (PSL, x1, y1, PSL_MOVE);
996 PSL_plotpoint (PSL, x2, y2, PSL_DRAW);
997 for (bin = n_bins; bin > 0; bin--) {
998 k = bin - 1;
999 az = k * Ctrl->A.inc - az_offset + half_bin_width;
1000 angle1 = 90.0 - az - Ctrl->A.inc;
1001 angle2 = angle1 + Ctrl->A.inc;
1002 PSL_plotarc (PSL, 0.0, 0.0, sum[k] * Ctrl->S.scale, angle1, angle2, (k == 0) ? PSL_STROKE : PSL_DRAW);
1003 }
1004 }
1005
1006 if (Ctrl->N.active) { /* Draw best-fitting probability distribution curve */
1007 unsigned int k, NP = 361;
1008 double a, mu, kappa, pdf, inc = 360.0 / (NP - 1), scl = 2.0 * area * Ctrl->A.inc * D2R;
1009 double *xp = gmt_M_memory (GMT, NULL, NP, double);
1010 double *yp = gmt_M_memory (GMT, NULL, NP, double);
1011
1012 mu = gmt_von_mises_mu_and_kappa (GMT, azimuth, length, n, &kappa);
1013 for (k = 0; k < NP; k++) {
1014 a = inc * k;
1015 pdf = scl * gmt_vonmises_pdf (GMT, a, mu, kappa);
1016 sincosd (start_angle - a, &s, &c);
1017 xp[k] = c * pdf;
1018 yp[k] = s * pdf;
1019 }
1020 gmt_setpen (GMT, &Ctrl->N.pen);
1021 PSL_plotline (PSL, xp, yp, NP, PSL_MOVE|PSL_STROKE);
1022 gmt_M_free (GMT, xp);
1023 gmt_M_free (GMT, yp);
1024 }
1025
1026 if (Ctrl->E.active) {
1027 unsigned int this_mode;
1028 int v4_outline = Ctrl->W.active[1];
1029 double *this_rgb = NULL;
1030 if (Ctrl->E.mode == GMT_OUT) { /* Write mean vector and statistics to file */
1031 uint64_t dim[GMT_DIM_SIZE] = {1, 1, 1, 8}; /* One record with 8 columns */
1032 struct GMT_DATASET *V = NULL;
1033 struct GMT_DATASEGMENT *S = NULL;
1034 if ((V = GMT_Create_Data (GMT->parent, GMT_IS_DATASET, GMT_IS_POINT, 0, dim, NULL, NULL, 0, 0, NULL)) == NULL) {
1035 Return (API->error);
1036 }
1037 S = V->table[0]->segment[0]; /* The only segment */
1038 sprintf (format, "mean_az\tmean_r\tmean_resultant\tmax\tscaled_mean_r\tlength_sum\tn\tsign@%.2f", Ctrl->Q.value);
1039 S->data[GMT_X][0] = mean_theta;
1040 S->data[GMT_Y][0] = mean_radius;
1041 S->data[GMT_Z][0] = mean_resultant;
1042 S->data[3][0] = max;
1043 S->data[4][0] = mean_radius;
1044 S->data[5][0] = total;
1045 S->data[6][0] = (double)n;
1046 S->data[7][0] = significant;
1047 S->n_rows = 1;
1048 if (GMT_Set_Comment (API, GMT_IS_DATASET, GMT_COMMENT_IS_COLNAMES, format, V)) {
1049 Return (API->error);
1050 }
1051 if (GMT_Write_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_WRITE_SET, NULL, Ctrl->E.file, V) != GMT_NOERROR) {
1052 Return (API->error);
1053 }
1054 }
1055 if (!Ctrl->W.active[1]) Ctrl->W.pen[1] = Ctrl->W.pen[0]; /* No separate pen specified; use same as for rose outline */
1056 if (Ctrl->E.mean) { /* Not given, calculate and use mean direction only */
1057 n_modes = 1;
1058 mode_direction = gmt_M_memory (GMT, NULL, 1, double);
1059 mode_length = gmt_M_memory (GMT, NULL, 1, double);
1060 mode_direction[0] = mean_theta;
1061 mode_length[0] = mean_radius;
1062 }
1063 else if (Ctrl->E.mode == GMT_IN) { /* Get mode parameters from separate file */
1064 if ((Cin = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_NONE, GMT_READ_NORMAL, NULL, Ctrl->E.file, NULL)) == NULL) {
1065 Return (API->error);
1066 }
1067 if (Cin->n_columns < 2) {
1068 GMT_Report (API, GMT_MSG_ERROR, "Input file %s has %d column(s) but at least 2 are needed\n", Ctrl->E.file, (int)Cin->n_columns);
1069 Return (GMT_DIM_TOO_SMALL);
1070 }
1071 T = Cin->table[0]; /* Can only be one table since we read a single file; We also only use the first segment */
1072 n_modes = (unsigned int)T->n_records;
1073 mode_direction = T->segment[0]->data[GMT_X];
1074 mode_length = T->segment[0]->data[GMT_Y];
1075 }
1076 if (!Ctrl->M.active) { /* Must supply defaults for the vector attributes */
1077 Ctrl->M.S.size_x = VECTOR_HEAD_LENGTH * GMT->session.u2u[GMT_PT][GMT_INCH]; /* 9p */
1078 Ctrl->M.S.v.v_width = (float)(VECTOR_LINE_WIDTH * GMT->session.u2u[GMT_PT][GMT_INCH]); /* 9p */
1079 Ctrl->M.S.v.v_angle = 30.0f;
1080 Ctrl->M.S.v.status |= (PSL_VEC_OUTLINE + PSL_VEC_OUTLINE2 + PSL_VEC_FILL + PSL_VEC_FILL2 + PSL_VEC_END);
1081 gmt_init_pen (GMT, &Ctrl->M.S.v.pen, VECTOR_LINE_WIDTH);
1082 gmt_init_fill (GMT, &Ctrl->M.S.v.fill, 0.0, 0.0, 0.0); /* Default vector fill = black */
1083 }
1084 gmt_init_vector_param (GMT, &Ctrl->M.S, false, false, NULL, false, NULL);
1085 if (Ctrl->M.S.symbol == PSL_VECTOR) Ctrl->M.S.v.v_width = (float)(Ctrl->W.pen[1].width * GMT->session.u2u[GMT_PT][GMT_INCH]);
1086 dim[PSL_VEC_TAIL_WIDTH] = Ctrl->M.S.v.v_width;
1087 dim[PSL_VEC_HEAD_LENGTH] = Ctrl->M.S.v.h_length;
1088 dim[PSL_VEC_HEAD_WIDTH] = Ctrl->M.S.v.h_width;
1089 dim[PSL_VEC_HEAD_SHAPE] = Ctrl->M.S.v.v_shape;
1090 dim[PSL_VEC_STATUS] = (double)Ctrl->M.S.v.status;
1091 dim[PSL_VEC_HEAD_TYPE_BEGIN] = (double)Ctrl->M.S.v.v_kind[0];
1092 dim[PSL_VEC_HEAD_TYPE_END] = (double)Ctrl->M.S.v.v_kind[1];
1093 if (Ctrl->M.S.v.status & PSL_VEC_OUTLINE) gmt_setpen (GMT, &Ctrl->W.pen[1]);
1094 if (Ctrl->M.S.v.status & PSL_VEC_FILL2) gmt_setfill (GMT, &Ctrl->M.S.v.fill, 1); /* Use fill structure */
1095 if (Ctrl->M.S.symbol == GMT_SYMBOL_VECTOR_V4) {
1096 dim[PSL_VEC_HEAD_SHAPE] = GMT->current.setting.map_vector_shape;
1097 if (Ctrl->M.S.v.status & PSL_VEC_FILL2)
1098 this_rgb = Ctrl->M.S.v.fill.rgb;
1099 else
1100 this_rgb = GMT->session.no_rgb;
1101 if (v4_outline) gmt_setpen (GMT, &Ctrl->W.pen[1]);
1102 }
1103 else {
1104 if (Ctrl->M.S.v.status & PSL_VEC_OUTLINE2) { /* Gave specific head outline pen */
1105 PSL_defpen (GMT->PSL, "PSL_vecheadpen", Ctrl->M.S.v.pen.width, Ctrl->M.S.v.pen.style, Ctrl->M.S.v.pen.offset, Ctrl->M.S.v.pen.rgb);
1106 dim[PSL_VEC_HEAD_PENWIDTH] = Ctrl->M.S.v.pen.width;
1107 }
1108 else if (Ctrl->M.active) {
1109 PSL_defpen (GMT->PSL, "PSL_vecheadpen", 0.5 * Ctrl->W.pen[1].width, Ctrl->W.pen[1].style, Ctrl->W.pen[1].offset, Ctrl->W.pen[1].rgb);
1110 dim[PSL_VEC_HEAD_PENWIDTH] = 0.5 * Ctrl->W.pen[1].width;
1111 }
1112 }
1113 for (this_mode = 0; this_mode < n_modes; this_mode++) {
1114 if (Ctrl->S.area_normalize) mode_length[this_mode] = sqrt (mode_length[this_mode]);
1115 if (half_only && mode_direction[this_mode] > 90.0 && mode_direction[this_mode] <= 270.0) mode_direction[this_mode] -= 180.0;
1116 angle = start_angle - mode_direction[this_mode];
1117 sincosd (angle, &s, &c);
1118 xr = Ctrl->S.scale * mode_length[this_mode] * c;
1119 yr = Ctrl->S.scale * mode_length[this_mode] * s;
1120 dim[PSL_VEC_XTIP] = xr, dim[PSL_VEC_YTIP] = yr;
1121 if (Ctrl->M.S.symbol == GMT_SYMBOL_VECTOR_V4)
1122 psl_vector_v4 (PSL, 0.0, 0.0, dim, this_rgb, v4_outline);
1123 else
1124 PSL_plotsymbol (PSL, 0.0, 0.0, dim, PSL_VECTOR);
1125 }
1126
1127 }
1128 if (GMT_End_IO (API, GMT_IN, 0) != GMT_NOERROR) { /* Disables further data input */
1129 Return (API->error);
1130 }
1131
1132 asize = GMT->current.setting.font_annot[GMT_PRIMARY].size * GMT->session.u2u[GMT_PT][GMT_INCH];
1133 lsize = GMT->current.setting.font_annot[GMT_PRIMARY].size * GMT->session.u2u[GMT_PT][GMT_INCH];
1134
1135 if (Ctrl->L.active) { /* Deactivate those with - */
1136 if (Ctrl->L.w[0] == '-' && Ctrl->L.w[1] == '\0') Ctrl->L.w[0] = '\0';
1137 if (Ctrl->L.e[0] == '-' && Ctrl->L.e[1] == '\0') Ctrl->L.e[0] = '\0';
1138 if (Ctrl->L.s[0] == '-' && Ctrl->L.s[1] == '\0') Ctrl->L.s[0] = '\0';
1139 if (Ctrl->L.n[0] == '-' && Ctrl->L.n[1] == '\0') Ctrl->L.n[0] = '\0';
1140 }
1141 else { /* Use default labels */
1142 Ctrl->L.w = strdup (GMT->current.language.cardinal_name[0][0]);
1143 Ctrl->L.e = strdup (GMT->current.language.cardinal_name[0][1]);
1144 Ctrl->L.s = strdup (GMT->current.language.cardinal_name[0][2]);
1145 Ctrl->L.n = strdup (GMT->current.language.cardinal_name[0][3]);
1146 }
1147 if ((GMT->common.B.active[GMT_PRIMARY] || GMT->common.B.active[GMT_SECONDARY]) && !GMT->current.map.frame.no_frame) {
1148 PSL_setcolor (PSL, GMT->current.setting.map_frame_pen.rgb, PSL_IS_STROKE);
1149 y = lsize + 6.0 * GMT->current.setting.map_annot_offset[GMT_PRIMARY];
1150 gmt_map_title (GMT, 0.0, off + y);
1151
1152 gmt_get_format (GMT, GMT->current.map.frame.axis[GMT_X].item[GMT_ANNOT_UPPER].interval, GMT->current.map.frame.axis[GMT_X].unit, GMT->current.map.frame.axis[GMT_X].prefix, format);
1153 if (do_labels) {
1154 if (half_only) {
1155 char text[GMT_LEN64] = {""};
1156 if (!Ctrl->L.active) { /* Use default labels */
1157 gmt_M_str_free (Ctrl->L.w);
1158 gmt_M_str_free (Ctrl->L.e);
1159 gmt_M_str_free (Ctrl->L.n);
1160 if (GMT->current.setting.map_degree_symbol == gmt_none) {
1161 if (half_only == 1) {
1162 sprintf (text, "90%s", GMT->current.language.cardinal_name[2][0]); Ctrl->L.w = strdup (text);
1163 sprintf (text, "90%s", GMT->current.language.cardinal_name[2][1]); Ctrl->L.e = strdup (text);
1164 sprintf (text, "0"); Ctrl->L.n = strdup (text);
1165 }
1166 else {
1167 sprintf (text, "0%s", GMT->current.language.cardinal_name[2][3]); Ctrl->L.w = strdup (text);
1168 sprintf (text, "180%s", GMT->current.language.cardinal_name[2][2]); Ctrl->L.e = strdup (text);
1169 sprintf (text, "90%s", GMT->current.language.cardinal_name[2][1]); Ctrl->L.n = strdup (text);
1170 }
1171 }
1172 else {
1173 if (half_only == 1) {
1174 sprintf (text, "90%c%s", (int)GMT->current.setting.ps_encoding.code[GMT->current.setting.map_degree_symbol], GMT->current.language.cardinal_name[2][0]); Ctrl->L.w = strdup (text);
1175 sprintf (text, "90%c%s", (int)GMT->current.setting.ps_encoding.code[GMT->current.setting.map_degree_symbol], GMT->current.language.cardinal_name[2][1]); Ctrl->L.e = strdup (text);
1176 sprintf (text, "0%c", (int)GMT->current.setting.ps_encoding.code[GMT->current.setting.map_degree_symbol]); Ctrl->L.n = strdup (text);
1177 }
1178 else {
1179 sprintf (text, "0%c%s", (int)GMT->current.setting.ps_encoding.code[GMT->current.setting.map_degree_symbol], GMT->current.language.cardinal_name[2][3]); Ctrl->L.w = strdup (text);
1180 sprintf (text, "180%c%s", (int)GMT->current.setting.ps_encoding.code[GMT->current.setting.map_degree_symbol], GMT->current.language.cardinal_name[2][2]); Ctrl->L.e = strdup (text);
1181 sprintf (text, "90%c%s", (int)GMT->current.setting.ps_encoding.code[GMT->current.setting.map_degree_symbol], GMT->current.language.cardinal_name[2][1]); Ctrl->L.n = strdup (text);
1182 }
1183 }
1184 }
1185 form = gmt_setfont (GMT, &GMT->current.setting.font_label);
1186 y = -(3.0 * GMT->current.setting.map_annot_offset[GMT_PRIMARY] + GMT->session.font[GMT->current.setting.font_annot[GMT_PRIMARY].id].height * asize);
1187 if (GMT->current.map.frame.axis[GMT_X].label[0]) PSL_plottext (PSL, 0.0, y, GMT->current.setting.font_label.size, GMT->current.map.frame.axis[GMT_X].label, 0.0, PSL_TC, form);
1188 y = -(5.0 * GMT->current.setting.map_annot_offset[GMT_PRIMARY] + GMT->session.font[GMT->current.setting.font_annot[GMT_PRIMARY].id].height * lsize + GMT->session.font[GMT->current.setting.font_label.id].height * lsize);
1189 if (GMT->current.map.frame.axis[GMT_Y].label[0]) PSL_plottext (PSL, 0.0, y, GMT->current.setting.font_label.size, GMT->current.map.frame.axis[GMT_Y].label, 0.0, PSL_TC, form);
1190 form = gmt_setfont (GMT, &GMT->current.setting.font_annot[GMT_PRIMARY]);
1191 n_annot = (GMT->current.map.frame.axis[GMT_X].item[GMT_ANNOT_UPPER].interval > 0.0) ? irint (max_radius / GMT->current.map.frame.axis[GMT_X].item[GMT_ANNOT_UPPER].interval) : -1;
1192 if (n_annot > 0) PSL_plottext (PSL, 0.0, -GMT->current.setting.map_annot_offset[GMT_PRIMARY], GMT->current.setting.font_annot[GMT_PRIMARY].size, "0", 0.0, PSL_TC, form);
1193 for (k = 1; n_annot > 0 && k <= n_annot; k++) {
1194 x = k * GMT->current.map.frame.axis[GMT_X].item[GMT_ANNOT_UPPER].interval;
1195 sprintf (text, format, x);
1196 x *= Ctrl->S.scale;
1197 PSL_plottext (PSL, x, -GMT->current.setting.map_annot_offset[GMT_PRIMARY], GMT->current.setting.font_annot[GMT_PRIMARY].size, text, 0.0, PSL_TC, form);
1198 PSL_plottext (PSL, -x, -GMT->current.setting.map_annot_offset[GMT_PRIMARY], GMT->current.setting.font_annot[GMT_PRIMARY].size, text, 0.0, PSL_TC, form);
1199 }
1200 }
1201 else {
1202 form = gmt_setfont (GMT, &GMT->current.setting.font_label);
1203 PSL_plottext (PSL, 0.0, -off - 2.0 * GMT->current.setting.map_annot_offset[GMT_PRIMARY], GMT->current.setting.font_label.size, Ctrl->L.s, 0.0, PSL_TC, form);
1204 if (!Ctrl->F.active && GMT->current.map.frame.axis[GMT_X].item[GMT_GRID_UPPER].interval > 0.0) { /* Draw scale bar but only if x-grid interval is set */
1205 double xp[4], yp[4];
1206 xp[0] = xp[1] = off; xp[2] = xp[3] = (max_radius - GMT->current.map.frame.axis[GMT_X].item[GMT_GRID_UPPER].interval) * Ctrl->S.scale;
1207 yp[0] = yp[3] = GMT->current.setting.map_tick_length[GMT_PRIMARY] - off; yp[1] = yp[2] = -off;
1208 gmt_setpen (GMT, &GMT->current.setting.map_tick_pen[GMT_PRIMARY]);
1209 PSL_plotline (PSL, xp, yp, 4, PSL_MOVE|PSL_STROKE);
1210 //PSL_plotsegment (PSL, off, -off, (max_radius - GMT->current.map.frame.axis[GMT_X].item[GMT_GRID_UPPER].interval) * Ctrl->S.scale, -off);
1211 //PSL_plotsegment (PSL, off, -off, off, GMT->current.setting.map_tick_length[GMT_PRIMARY] - off);
1212 //PSL_plotsegment (PSL, (max_radius - GMT->current.map.frame.axis[GMT_X].item[GMT_GRID_UPPER].interval) * Ctrl->S.scale, -off, (max_radius - GMT->current.map.frame.axis[GMT_X].item[GMT_GRID_UPPER].interval) * Ctrl->S.scale, GMT->current.setting.map_tick_length[GMT_PRIMARY] - off);
1213 if (GMT->current.map.frame.axis[GMT_X].label[0]) {
1214 strcat (format, " %s");
1215 sprintf (text, format, GMT->current.map.frame.axis[GMT_X].item[GMT_GRID_UPPER].interval, GMT->current.map.frame.axis[GMT_X].label);
1216 }
1217 else
1218 sprintf (text, format, GMT->current.map.frame.axis[GMT_X].item[GMT_GRID_UPPER].interval);
1219 form = gmt_setfont (GMT, &GMT->current.setting.font_annot[GMT_PRIMARY]);
1220 PSL_plottext (PSL, (max_radius - 0.5 * GMT->current.map.frame.axis[GMT_X].item[GMT_GRID_UPPER].interval) * Ctrl->S.scale, -(off + GMT->current.setting.map_annot_offset[GMT_PRIMARY]), GMT->current.setting.font_annot[GMT_PRIMARY].size, text, 0.0, PSL_TC, form);
1221 }
1222 y = -(off + 5.0 * GMT->current.setting.map_annot_offset[GMT_PRIMARY] + GMT->session.font[GMT->current.setting.font_annot[GMT_PRIMARY].id].height * lsize + GMT->session.font[GMT->current.setting.font_label.id].height * lsize);
1223 form = gmt_setfont (GMT, &GMT->current.setting.font_label);
1224 if (GMT->current.map.frame.axis[GMT_Y].label[0]) PSL_plottext (PSL, 0.0, y, GMT->current.setting.font_label.size, GMT->current.map.frame.axis[GMT_Y].label, 0.0, PSL_TC, form);
1225 }
1226 form = gmt_setfont (GMT, &GMT->current.setting.font_label);
1227 PSL_plottext (PSL, off + 2.0 * GMT->current.setting.map_annot_offset[GMT_PRIMARY], 0.0, GMT->current.setting.font_label.size, Ctrl->L.e, 0.0, PSL_ML, form);
1228 PSL_plottext (PSL, -off - 2.0 * GMT->current.setting.map_annot_offset[GMT_PRIMARY], 0.0, GMT->current.setting.font_label.size, Ctrl->L.w, 0.0, PSL_MR, form);
1229 PSL_plottext (PSL, 0.0, off + 2.0 * GMT->current.setting.map_annot_offset[GMT_PRIMARY], GMT->current.setting.font_label.size, Ctrl->L.n, 0.0, PSL_BC, form);
1230 PSL_setcolor (PSL, GMT->current.setting.map_default_pen.rgb, PSL_IS_STROKE);
1231 }
1232 }
1233
1234 if (GMT->common.B.active[0] && !GMT->current.map.frame.no_frame ) { /* Finally draw frame */
1235 int symbol = (half_only) ? PSL_WEDGE : PSL_CIRCLE;
1236 double dim[PSL_MAX_DIMS];
1237 struct GMT_FILL no_fill;
1238
1239 gmt_init_fill (GMT, &no_fill, -1.0, -1.0, -1.0);
1240 gmt_M_memset (dim, PSL_MAX_DIMS, double);
1241 dim[0] = (half_only) ? 0.5 * diameter : diameter;
1242 dim[1] = 0.0;
1243 dim[2] = (half_only) ? 180.0 : 360.0;
1244 dim[7] = 2; /* Do draw line */
1245 gmt_setpen (GMT, &GMT->current.setting.map_frame_pen);
1246 gmt_setfill (GMT, &no_fill, 1);
1247 PSL_plotsymbol (PSL, 0.0, 0.0, dim, symbol);
1248 }
1249
1250 gmt_plane_perspective (GMT, -1, 0.0);
1251 PSL_setorigin (PSL, -x_origin, -y_origin, 0.0, PSL_INV);
1252 gmt_plotend (GMT);
1253
1254 gmt_M_free (GMT, sum);
1255 gmt_M_free (GMT, xx);
1256 gmt_M_free (GMT, yy);
1257 gmt_M_free (GMT, azimuth);
1258 gmt_M_free (GMT, length);
1259 if (Ctrl->E.active) {
1260 if (Ctrl->E.mean) {
1261 gmt_M_free (GMT, mode_length);
1262 gmt_M_free (GMT, mode_direction);
1263 }
1264 }
1265
1266 Return (GMT_NOERROR);
1267 }
1268
GMT_rose(void * V_API,int mode,void * args)1269 EXTERN_MSC int GMT_rose (void *V_API, int mode, void *args) {
1270 /* This is the GMT6 modern mode name */
1271 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
1272 if (API->GMT->current.setting.run_mode == GMT_CLASSIC && !API->usage) {
1273 GMT_Report (API, GMT_MSG_ERROR, "Shared GMT module not found: rose\n");
1274 return (GMT_NOT_A_VALID_MODULE);
1275 }
1276 return GMT_psrose (V_API, mode, args);
1277 }
1278