Fix wrong column selections in rotsmoother (#6795)

* Fix wrong column selections in rotsmoother

Per documentation, input format shall be lon, lat, time, angle, [weight], but code flipped time and angle.  It also did not allow a single time (since rotations usually correspond to an interval).  Now, if -T<single_t> is used we include all rotations in the calculation of a single rotation.

* Update rotsmoother.c

Fix index of last rotation outside group and also give message if not enough rotations found for a group.

doc/rst/source/supplements/spotter/rotsmoother.rst

@@ -39,13 +39,13 @@ Description
 
 **rotsmoother** reads a table of total reconstructions and computes mean
 rotations (and optionally covariance matrices) for sub-groups of rotations
-based on rotation age.
+based on rotation age (or all of the rotations).
 
 Required Arguments
 ------------------
 
 *table*
-    Name of a rotation table containing (lon, lat, time, angle, [weight]) values.
+    Name of a rotation table containing (*lon lat time angle* [*weight*]) values.
 
 Optional Arguments
 ------------------
@@ -80,12 +80,13 @@ Optional Arguments
 .. _-T:
 
 **-T**\ *ages*
-    Sets the desired groups of ages.  For a single time append
-    the desired time.  For an equidistant range of reconstruction ages
-    give **-T**\ *start*\ /\ *stop*\ /\ *inc* or **-T**\ *start*\ /\ *stop*\ /\ *npoints*\ **+n**.
+    Sets the desired groups of ages.  For a single time append the desired time
+    and all rotations are used to compute a single mean rotation.  For an equidistant
+    range of reconstruction ages give **-T**\ *start*\ /\ *stop*\ /\ *inc* or
+    **-T**\ *start*\ /\ *stop*\ /\ *npoints*\ **+n**.
     For an non-equidistant set of reconstruction ages please pass them
-    via the first column in a file, e.g., **-T**\ *agefile*.  The ages indicate
-    read or generated becomes the bin-boundaries and we output the average time of
+    via the first column in a file, e.g., **-T**\ *agefile*.  The ages we
+    read or generate becomes the bin-boundaries and we output the average time of
     all rotations inside each bin.
 
 .. |Add_-V| replace:: |Add_-V_links|
@@ -132,12 +133,14 @@ Optional Arguments
 Examples
 --------
 
-To smooth rotation groups in increments of 3 Myr and ensure northern hemisphere poles, try
-
-   ::
+To smooth rotation groups in increments of 3 Myr and ensure northern hemisphere poles, try::
 
     gmt rotsmoother rotations.txt -N -T3/3/30 -V > rot_means.txt
 
+To smooth all rotations and compute a single mean rotation (assigned to time = 5) with
+corresponding covariance matrix, try::
+
+    gmt rotsmoother rotations.txt -C -T5 -V > rot_means_cov.txt
 
 See Also
 --------

src/spotter/rotsmoother.c

@@ -263,7 +263,7 @@ EXTERN_MSC int GMT_rotsmoother (void *V_API, int mode, void *args) {
 	bool stop;
 	uint64_t n_read = 0, rot, p, first = 0, last, n_use = 0, n_out = 0, n_total_use = 0, n_minimum, n_alloc = GMT_CHUNK;
 	int error = 0, n_fields;
-	unsigned int n_in = 3, k, j, t_col, w_col, t, n_cols = 4, matrix_dim = 3, nrots;
+	unsigned int n_in = 3, k, j, a_col, w_col, start_t, t, n_cols = 4, matrix_dim = 3, nrots;
 	double *in = NULL, min_rot_angle, min_rot_age, max_rot_angle, max_rot_age;
 	double sum_rot_angle, sum_rot_angle2, sum_rot_age, sum_rot_age2, sum_weights;
 	double out[20], khat = 1.0, g = 1.0e-5;	/* Common scale factor for all Covariance terms */
@@ -304,7 +304,7 @@ EXTERN_MSC int GMT_rotsmoother (void *V_API, int mode, void *args) {
 
 	GMT_Report (API, GMT_MSG_INFORMATION, "Processing input table data\n");
 
-	if (!Ctrl->A.active) n_in++;	/* Got time */
+	if (!Ctrl->A.active) n_in++;	/* Got time and angle */
 	if (Ctrl->W.active) n_in++;		/* Got weights */
 	if (Ctrl->C.active) n_cols = 19;	/* Want everything */
 
@@ -319,8 +319,9 @@ EXTERN_MSC int GMT_rotsmoother (void *V_API, int mode, void *args) {
 	if (GMT_Begin_IO (API, GMT_IS_DATASET,  GMT_IN, GMT_HEADER_ON) != GMT_NOERROR) {	/* Enables data input and sets access mode */
 		Return (API->error);
 	}
-	t_col = (Ctrl->A.active) ? GMT_Z : 3;	/* If no time we use angle as proxy for time */
-	w_col = t_col + 1;
+	/* time (or angle proxy) is always in column 2 */
+	a_col = (Ctrl->A.active) ? GMT_Z : 3;	/* Angle is in column 2 or 3 */
+	w_col = a_col + 1;	/* Optional final column */
 	D = (struct ROTSMOOTHER_AGEROT *) gmt_M_memory (GMT, NULL, n_alloc, struct ROTSMOOTHER_AGEROT);
 	Out = gmt_new_record (GMT, out, NULL);	/* Since we only need to worry about numerics in this module */
 
@@ -352,8 +353,8 @@ EXTERN_MSC int GMT_rotsmoother (void *V_API, int mode, void *args) {
 		/* Convert to geocentric, load parameters  */
 		D[n_read].wxyasn[K_LON]    = in[GMT_X];
 		D[n_read].wxyasn[K_LAT]    = gmt_lat_swap (GMT, in[GMT_Y], GMT_LATSWAP_G2O);
-		D[n_read].wxyasn[K_ANGLE]  = in[GMT_Z];
-		D[n_read].wxyasn[K_AGE]    = in[t_col];
+		D[n_read].wxyasn[K_ANGLE]  = in[a_col];
+		D[n_read].wxyasn[K_AGE]    = in[GMT_Z];
 		D[n_read].wxyasn[K_WEIGHT] = (Ctrl->W.active) ? in[w_col] : 1.0;	/* Optionally use weights */
 		n_read++;
 		if (n_read == n_alloc) {	/* Need larger arrays */
@@ -429,17 +430,29 @@ EXTERN_MSC int GMT_rotsmoother (void *V_API, int mode, void *args) {
 	z_unit_vector[0] = z_unit_vector[1] = 0.0;	z_unit_vector[2] = 1.0;	/* The local z unit vector */
 	n_minimum = (Ctrl->C.active) ? 2 : 1;	/* Need at least two rotations to compute covariance, at least one to report the mean */
 
-	for (t = 1; t < Ctrl->T.n_times; t++) {	/* For each desired output time interval */
-		t_lo = Ctrl->T.value[t-1];	t_hi = Ctrl->T.value[t];	/* The current interval */
-		for (rot = first, stop = false; !stop && rot < n_read; rot++)	/* Determine index of first rotation inside this age window */
-			if (D[rot].wxyasn[K_AGE] >= t_lo) stop = true;
-		first = rot - 1;	/* Index to first rotation inside this time interval */
-		for (stop = false; !stop && rot < n_read; rot++)	/* Determine index of last rotation inside this age window */
-			if (D[rot].wxyasn[K_AGE] > t_hi) stop = true;
-		last = rot - 1;	/* Index to first rotation outside this time interval */
-		n_use = last - first;	/* Number of rotations in the interval */
-		GMT_Report (API, GMT_MSG_INFORMATION, "Found %d rots for the time interval %g <= t < %g\n", n_use, t_lo, t_hi);
-		if (n_use < n_minimum) continue;	/* Need at least 1 or 2 poles to do anything useful */
+	start_t = (Ctrl->T.n_times > 1) ? 1 : 0;	/* For a fixed single time there is no interval */
+	for (t = start_t; t < Ctrl->T.n_times; t++) {	/* For each desired output time interval */
+		if (start_t) {	/* We have intervals */
+			t_lo = Ctrl->T.value[t-1];	t_hi = Ctrl->T.value[t];	/* The current interval */
+			for (rot = first, stop = false; !stop && rot < n_read; rot++)	/* Determine index of first rotation inside this age window */
+				if (D[rot].wxyasn[K_AGE] >= t_lo) stop = true;
+			first = rot - 1;	/* Index to first rotation inside this time interval */
+			for (stop = false; !stop && rot < n_read; rot++)	/* Determine index of last rotation inside this age window */
+				if (D[rot].wxyasn[K_AGE] > t_hi) stop = true;
+			last = (stop) ? rot - 1 : n_read;	/* Index to first rotation outside this time interval */
+			n_use = last - first;	/* Number of rotations in the interval */
+			GMT_Report (API, GMT_MSG_INFORMATION, "Found %d rots for the time interval %g <= t < %g\n", n_use, t_lo, t_hi);
+		}
+		else {
+			t_lo = Ctrl->T.value[0];	t_hi = Ctrl->T.value[0];	/* The current time */
+			first = 0;	last = n_read;
+			n_use = n_read;	/* Number of rotations in the interval */
+			GMT_Report (API, GMT_MSG_INFORMATION, "Found %d rots for time = %g\n", n_use, t_lo);
+		}
+		if (n_use < n_minimum) {	/* Need at least 1 or 2 poles to do anything useful */
+			GMT_Report (API, GMT_MSG_INFORMATION, "Not enough rotations to compute anything - skipping this group\n");
+			continue;
+		}
 
 		/* Now estimate the average rotation */