Let grdtrend fit model along xx or yy only. (#5496)

* Let grdtrend fit model along xx or yy only.

* Update docs and add a test script

* Finalize

* Polish

* Update grdtrend.rst

* Update grdtrend.rst

* Update decision on trivial solution

The cutoff is different for +x or +y since we are doing a linear fit only.

* Update grdtrend.c

Co-authored-by: Paul Wessel <pwessel@hawaii.edu>

doc/rst/source/grdtrend.rst

@@ -12,7 +12,7 @@ Synopsis
 
 .. include:: common_SYN_OPTs.rst_
 
-**gmt grdtrend** *grdfile* |-N|\ *n\_model*\ [**+r**]
+**gmt grdtrend** *grdfile* |-N|\ *n_model*\ [**+r**][**+x**\|\ **y**]
 [ |-D|\ *diff.nc* ]
 [ |SYN_OPT-R| ]
 [ |-T|\ *trend.nc* ] [ |-W|\ *weight.nc* ]
@@ -30,14 +30,16 @@ is defined by:
 .. math::
     m_1 + m_2x + m_3y + m_4xy + m_5x^2 + m_6y^2 + m_7x^3 + m_8x^2y + m_9xy^2 + m_{10}y^3.
 
-The user must specify **-N**\ *n\_model*, the number of model parameters
+The user must specify **-N**\ *n_model*, the number of model parameters
 to use; thus, **-N**\ *3* fits a bilinear trend, **-N**\ *6* a quadratic
 surface, and so on. Optionally, append **+r** to the **-N** option to
 perform a robust fit. In this case, the program will iteratively
 reweight the data based on a robust scale estimate, in order to converge
 to a solution insensitive to outliers. This may be handy when separating
 a "regional" field from a "residual" which should have non-zero mean,
 such as a local mountain on a regional surface.
+Optionally, you may choose to fit a trend that varies only along the *x* or *y* axis,
+in which case you select an *n_model* from 1 (constant) to 4 (cubic).
 
 If data file has values set to NaN, these will be ignored during
 fitting; if output files are written, these will also have NaN in the
@@ -51,9 +53,13 @@ Required Arguments
 
 .. _-N:
 
-**-N**\ *n\_model*\ [**+r**]
-    *n\_model* sets the number of model parameters to fit.
-    Append **+r** for robust fit.
+**-N**\ *n_model*\ [**+r**][**+x**\|\ **y**]
+    *n_model* sets the ID of the highest model parameters to fit.
+    Append **+r** for robust fit.  As an option, append either **+x** or **+y** to only
+    fit a model that depends on *x* or *y* terms, respectively. This means we either fit
+    :math:`m_1 + m_2x + m_3x^2 + m_4x^3` or :math:`m_1 + m_2y + m_3y^2 + m_4y^3`.
+    Note that *n_model* may only be 1-4 for the one-dimensional fits but may be 1-10
+    for the two-dimensional surface fits.
 
 Optional Arguments
 ------------------
@@ -98,7 +104,7 @@ basis functions are built from Legendre polynomials. These have a
 numerical advantage in the form of the matrix which must be inverted and
 allow more accurate solutions. NOTE: The model parameters listed with
 **-V** are Legendre polynomial coefficients; they are not numerically
-equivalent to the m#s in the equation described above. The description
+equivalent to the :math:`m_j` in the equation described above. The description
 above is to allow the user to match **-N** with the order of the
 polynomial surface. See :doc:`grdmath` if you need to evaluate the trend
 using the reported coefficients.