diff --git a/docs/source/morphpy.rst b/docs/source/morphpy.rst
index 8031b3b..1fd4f1c 100644
--- a/docs/source/morphpy.rst
+++ b/docs/source/morphpy.rst
@@ -106,19 +106,26 @@ exclude: list of str
     Exclude a manipulations from refinement by name
     (e.g. exclude=["scale", "stretch"] excludes the scale and stretch morphs).
 scale: float
-    Apply scale factor. This multiplies the function ordinate by scale.
+    Apply scale factor.
+
+    This multiplies the function ordinate by scale.
 stretch: float
-    Stretch function grid by a fraction stretch. Specifically, this multiplies the function grid by 1+stretch.
+    Stretch function grid by a fraction stretch.
+
+    This multiplies the function grid by 1+stretch.
 squeeze: list of float
     Squeeze function grid given a polynomial
-    p(x) = squeeze[0]+squeeze[1]*x+...+squeeze[n]*x^n. n is dependent on the number
+    p(x) = squeeze[0]+squeeze[1]*x+...+squeeze[n]*x^n.
+
+    n is dependent on the number
     of values in the user-inputted comma-separated list.
     The morph transforms the function grid from x to x+p(x).
     When this parameter is given, hshift is disabled.
     When n>1, stretch is disabled.
 smear: float
-    Smear the peaks with a Gaussian of width smear. This
-    is done by convolving the function with a Gaussian
+    Smear the peaks with a Gaussian of width smear.
+
+    This is done by convolving the function with a Gaussian
     with standard deviation smear. If both smear and
     smear_pdf are used, only smear_pdf will be
     applied.
@@ -129,6 +136,7 @@ smear_pdf: float
     applied.
 slope: float
     Slope of the baseline used in converting from PDF to RDF.
+
     This is used with the option smear_pdf. The slope will
     be estimated if not provided.
 hshift: float
@@ -139,39 +147,99 @@ qdamp: float
     Dampen PDF by a factor qdamp.
 radius: float
     Apply characteristic function of sphere with radius
-    given by parameter radius. If pradius is also specified, instead apply
+    given by parameter radius.
+
+    If pradius is also specified, instead apply
     characteristic function of spheroid with equatorial
     radius radius and polar radius pradius.
 pradius: float
     Apply characteristic function of spheroid with
     equatorial radius given by above parameter radius and polar radius pradius.
+
     If only pradius is specified, instead apply
     characteristic function of sphere with radius pradius.
 iradius: float
     Apply inverse characteristic function of sphere with
-    radius iradius. If ipradius is also specified, instead
+    radius iradius.
+
+    If ipradius is also specified, instead
     apply inverse characteristic function of spheroid with
     equatorial radius iradius and polar radius ipradius.
 ipradius: float
     Apply inverse characteristic function of spheroid with
     equatorial radius iradius and polar radius ipradius.
+
     If only ipradius is specified, instead apply inverse
     characteristic function of sphere with radius ipradius.
 funcy: tuple (function, dict)
+    Apply a function to the y-axis of the (two-column) data.
+
     This morph applies the function funcy[0] with parameters given in funcy[1].
-    The function funcy[0] must be a function of both the abscissa and ordinate
-    (e.g. take in at least two inputs with as many additional parameters as needed).
+    The function funcy[0] take in as parameters both the abscissa and ordinate
+    (i.e. take in at least two inputs with as many additional parameters as needed).
+    The y-axis values of the data are then replaced by the return value of funcy[0].
+
     For example, let's start with a two-column table with abscissa x and ordinate y.
     let us say we want to apply the function ::
 
         def linear(x, y, a, b, c):
             return a * x + b * y + c
 
-    This function takes in both the abscissa and ordinate on top of three additional
-    parameters a, b, and c. To use the funcy parameter with initial guesses
-    a=1.0, b=2.0, c=3.0, we would pass ``funcy=(linear, {a: 1.0, b: 2.0, c: 3.0})``.
-    For an example use-case, see the Python-Specific Morphs section below.
+    This example function above takes in both the abscissa and ordinate on top of
+    three additional parameters a, b, and c.
+    To use the funcy parameter with parameter values a=1.0, b=2.0, and c=3.0,
+    we would pass ``funcy=(linear, {"a": 1.0, "b": 2.0, "c": 3.0})``.
+    For an explicit example, see the Python-Specific Morphs section below.
+funcx: tuple (function, dict)
+    Apply a function to the x-axis of the (two-column) data.
+
+    This morph works fundamentally differently from the other grid morphs
+    (e.g. stretch and squeeze) as it directly modifies the grid of the
+    morph function.
+    The other morphs maintain the original grid and apply the morphs by interpolating
+    the function ***.
+
+    This morph applies the function funcx[0] with parameters given in funcx[1].
+    The function funcx[0] take in as parameters both the abscissa and ordinate
+    (i.e. take in at least two inputs with as many additional parameters as needed).
+    The x-axis values of the data are then replaced by the return value of funcx[0].
+    Note that diffpy.morph requires the x-axis be monotonic increasing
+    (i.e. for i < j, x[i] < x[j]): as such,
+    if funcx[0] is not a monotonic increasing function of the provided x-axis data,
+    the error ``x must be a strictly increasing sequence`` will be thrown.
+
+    For example, let's start with a two-column table with abscissa x and ordinate y.
+    let us say we want to apply the function ::
+
+        def exponential(x, y, amp, decay):
+            return abs(amp) * (1 - 2**(-decay * x))
 
+    This example function above takes in both the abscissa and ordinate on top of
+    three additional parameters amp and decay.
+    (Even though the ordinate is not used in the function,
+    it is still required that the function take in both acscissa and ordinate.)
+    To use the funcx parameter with parameter values amp=1.0 and decay=2.0,
+    we would pass ``funcx=(exponential, {"amp": 1.0, "decay:: 2.0})``.
+    For an explicit example, see the Python-Specific Morphs section below.
+funcxy: tuple (function, dict)
+    Apply a function the (two-column) data.
+
+    This morph applies the function funcxy[0] with parameters given in funcxy[1].
+    The function funcxy[0] take in as parameters both the abscissa and ordinate
+    (i.e. take in at least two inputs with as many additional parameters as needed).
+    The two columns of the data are then replaced by the two return values of funcxy[0].
+
+    For example, let's start with a two-column table with abscissa x and ordinate y.
+    let us say we want to apply the function ::
+
+        def shift(x, y, hshift, vshift):
+            return x + hshift, y + vshift
+
+    This example function above takes in both the abscissa and ordinate on top of
+    two additional parameters hshift and vshift.
+    To use the funcy parameter with parameter values hshift=1.0 and vshift=2.0,
+    we would pass ``funcy=(shift, {"hshift": 1.0, "vshift": 1.0})``.
+    For an example use-case, see the Python-Specific Morphs section below.
 
 Python-Specific Morphs
 ======================
@@ -179,17 +247,25 @@ Python-Specific Morphs
 Some morphs in ``diffpy.morph`` are supported only in Python. Here, we detail
 how they are used and how to call them.
 
-MorphFuncy: Applying custom functions
+MorphFunc: Applying custom functions
 -------------------------------------
 
+In these tutorial, we walk through how to use the ``MorphFunc`` morphs
+(``MorphFuncy``, ``MorphFuncx``, ``MorphFuncxy``)
+with some example transformations.
+
+Unlike other morphs that can be run from the command line,
+``MorphFunc`` moprhs require a Python function and is therefore
+intended to be used through Python scripting.
+
+MorphFuncy:
+^^^^^^^^^^^
+
 The ``MorphFuncy`` morph allows users to apply a custom Python function
 to the y-axis values of a dataset, enabling flexible and user-defined
 transformations.
 
-In this tutorial, we walk through how to use ``MorphFuncy`` with an example
-transformation. Unlike other morphs that can be run from the command line,
-``MorphFuncy`` requires a Python function and is therefore intended to be used
-through Python scripting.
+Let's try out this morph!
 
     1. Import the necessary modules into your Python script:
 
@@ -230,7 +306,7 @@ through Python scripting.
 
        .. code-block:: python
 
-            morph_params, morph_table = morph_arrays(np.array([x_morph, y_morph]).T,np.array([x_target, y_target]).T,
+            morph_params, morph_table = morph_arrays(np.array([x_morph, y_morph]).T, np.array([x_target, y_target]).T,
             funcy=(linear_function,{'scale': 1.2, 'offset': 0.1}))
 
     5. Extract the fitted parameters from the result:
@@ -246,3 +322,218 @@ to generate the target (scale=20 & offset=0.8). This example shows how
 ``MorphFuncy`` can be used to fit and apply custom transformations. Now
 it's your turn to experiment with other custom functions that may be useful
 for analyzing your data.
+
+MorphFuncx:
+^^^^^^^^^^^
+
+The ``MorphFuncx`` morph allows users to apply a custom Python function
+to the x-axis values of a dataset, similar to the ``MorphFuncy`` morph.
+
+One caveat to this morph is that the x-axis values must remain monotonic
+increasing, so it is possible to run into errors when applying this morph.
+For example, if your initial grid is ``[-1, 0, 1]``, and your function is
+``lambda x, y: x**2``, the grid after the function is applied will be
+``[1, 0, 1]``, which is no longer monotonic increasing.
+In this case, the error ``x must be a strictly increasing sequence``
+will be thrown.
+
+Let's try out this morph!
+
+    1. Import the necessary modules into your Python script:
+
+       .. code-block:: python
+
+            from diffpy.morph.morphpy import morph_arrays
+            import numpy as np
+
+    2. Define a custom Python function to apply a transformation to the data.
+       The function must take ``x`` and ``y`` (1D arrays of the same length)
+       along with named parameters, and return a transformed ``x`` array of the
+       same length. Recall that this function must maintain the monotonic
+       increasing nature of the ``x`` array.
+
+       For this example, we will use a simple exponential function transformation that
+       greatly modifies the input:
+
+       .. code-block:: python
+
+            def exp_function(x, y, scale, rate):
+                return np.abs(scale) * np.exp(np.abs(rate) * x)
+
+       Notice that, though the function only uses the ``x`` input,
+       the function signature takes in both ``x`` and ``y``.
+
+    3. Like in the previous example, we will use a sine function for the morph
+       data and generate the target data by applying the decay transfomration
+       with a known scale and rate:
+
+       .. code-block:: python
+
+            x_morph = np.linspace(0, 10, 1001)
+            y_morph = np.sin(x_morph)
+            x_target = x_target = 20 * np.exp(0.8 * x_morph)
+            y_target = y_morph.copy()
+
+    4. Setup and run the morph using the ``morph_arrays(...)``.
+       ``morph_arrays`` expects the morph and target data as **2D arrays** in
+       *two-column* format ``[[x0, y0], [x1, y1], ...]``. This will apply
+       the user-defined function and refine the parameters to best align the
+       morph data with the target data. This includes both the transformation
+       parameters (our initial guess) and the transformation function itself:
+
+       .. code-block:: python
+
+            morph_params, morph_table = morph_arrays(np.array([x_morph, y_morph]).T, np.array([x_target, y_target]).T,
+            funcx=(decay_function, {'scale': 1.2, 'rate': 1.0}))
+
+    5. Extract the fitted parameters from the result:
+
+       .. code-block:: python
+
+            fitted_params = morph_params["funcx"]
+            print(f"Fitted scale: {fitted_params['scale']}")
+            print(f"Fitted rate: {fitted_params['rate']}")
+
+Again, we should see that the fitted scale and offset values match the ones used
+to generate the target (scale=20 & rate=0.8).
+
+For fun, you can plot the original function to the morphed function to see
+how much the
+
+MorphFuncxy:
+^^^^^^^^^^^^
+The ``MorphFuncxy`` morph allows users to apply a custom Python function
+to a dataset that modifies both the ``x`` and ``y`` column values.
+This is equivalent to applying a ``MorphFuncx`` and ``MorphFuncy``
+simultaneously.
+
+This morph is useful when you want to apply operations that modify both
+the grid and function value. A PDF-specific example includes computing
+PDFs from 1D diffraction data (see paragraph at the end of this section).
+
+For this tutorial, we will go through two examples. One simple one
+involving shifting a function in the ``x`` and ``y`` directions, and
+another involving a Fourier transform.
+
+    1. Let's start by taking a simple ``sine`` function.
+
+       .. code-block:: python
+
+            import numpy as np
+            morph_x = np.linspace(0, 10, 101)
+            morph_y = np.sin(morph_x)
+            morph_table = np.array([morph_x, morph_y]).T
+
+    2. Then, let our target function be that same ``sine`` function shifted
+       to the right by ``0.3`` and up by ``0.7``.
+
+       .. code-block:: python
+
+            target_x = morph_x + 0.3
+            target_y = morph_y + 0.7
+            target_table = np.array([target_x, target_y]).T
+
+    3. While we could use the ``hshift`` and ``vshift`` morphs,
+       this would require us to refine over two separate morph
+       operations. We can instead perform these morphs simultaneously
+       by defining a function:
+
+       .. code-block:: python
+
+            def shift(x, y, hshift, vshift):
+                return x + hshift, y + vshift
+
+    4. Now, let's try finding the optimal shift parameters using the ``MorphFuncxy`` morph.
+       We can try an initial guess of ``hshift=0.0`` and ``vshift=0.0``.
+
+       .. code-block:: python
+
+            from diffpy.morph.morphpy import morph_arrays
+            initial_guesses = {"hshift": 0.0, "vshift": 0.0}
+            info, table = morph_arrays(morph_table, target_table, funcxy=(shift, initial_guesses))
+
+    5. Finally, to see the refined ``hshift`` and ``vshift`` parameters, we extract them from ``info``.
+
+       .. code-block:: python
+
+            print(f"Refined hshift: {info["funcxy"]["hshift"]}")
+            print(f"Refined vshift: {info["funcxy"]["vshift"]}")
+
+Now for an example involving a Fourier transform.
+
+    1. Let's say you measured a signal of the form :math:`f(x)=\exp\{\cos(\pi x)\}`.
+       Unfortunately, your measurement was taken against a noisy sinusoidal
+       background of the form :math:`n(x)=A\sin(Bx)`, where ``A``, ``B`` are unknown.
+       For our example, let's say (unknown to us) that ``A=2`` and ``B=1.7``.
+
+       .. code-block:: python
+
+            import numpy as np
+            n = 201
+            dx = 0.01
+            measured_x = np.linspace(0, 2, n)
+
+            def signal(x):
+                return np.exp(np.cos(np.pi * x))
+
+            def noise(x, A, B):
+                return A * np.sin(B * x)
+
+            measured_f = signal(measured_x) + noise(measured_x, 2, 1.7)
+            morph_table = np.array([measured_x, measured_f]).T
+
+    2. Your colleague remembers they previously computed the Fourier transform
+       of the function and has sent that to you.
+
+       .. code-block:: python
+
+            # We only consider the region where the grid is positive for simplicity
+            target_x = np.fft.fftfreq(n, dx)[:n//2]
+            target_f = np.real(np.fft.fft(signal(measured_x))[:n//2])
+            target_table = np.array([target_x, target_f]).T
+
+    3. We can now write a noise subtraction function that takes in our measured
+       signal and guesses for parameters ``A``, ``B``, and computes the Fourier
+       transform post-noise-subtraction.
+
+       .. code-block:: python
+
+            def noise_subtracted_ft(x, y, A, B):
+                n = 201
+                dx = 0.01
+                background_subtracted_y = y - noise(x, A, B)
+
+                ft_x = np.fft.fftfreq(n, dx)[:n//2]
+                ft_f = np.real(np.fft.fft(background_subtracted_y)[:n//2])
+
+                return ft_x, ft_f
+
+    4. Finally, we can provide initial guesses of ``A=0`` and ``B=1`` to the
+       ``MorphFuncxy`` morph and see what refined values we get.
+
+       .. code-block:: python
+
+            from diffpy.morph.morphpy import morph_arrays
+            initial_guesses = {"A": 0, "B": 1}
+            info, table = morph_arrays(morph_table, target_table, funcxy=(background_subtracted_ft, initial_guesses))
+
+    5. Print these values to see if they match with the true values of
+       of ``A=2.0`` and ``B=1.7``!
+
+       .. code-block:: python
+
+            print(f"Refined A: {info["funcxy"]["A"]}")
+            print(f"Refined B: {info["funcxy"]["B"]}")
+
+You can also use this morph to help find optimal parameters
+(e.g. ``rpoly``, ``qmin``, ``qmax``, ``bgscale``) for computing
+PDFs of materials with known structures.
+One does this by setting the ``MorphFuncxy`` function to a PDF
+computing function such as
+`PDFgetx3 <https://www.diffpy.org/products/pdfgetx.html>`_.
+The input (morphed) 1D function should be the 1D diffraction data
+one wishes to compute the PDF of and the target 1D function
+can be the PDF of a target material with similar geometry.
+More information about this will be released in the ``diffpy.morph``
+manuscript, and we plan to integrate this feature automatically into
+``PDFgetx3`` soon.
diff --git a/news/morphfuncxy.rst b/news/morphfuncxy.rst
new file mode 100644
index 0000000..ef89e30
--- /dev/null
+++ b/news/morphfuncxy.rst
@@ -0,0 +1,24 @@
+**Added:**
+
+* morphfuncx added: apply a function to the grid of your morphed function; this function should maintain the monotonic increasing nature of the grid
+* morphfuncxy added: apply a general function which can modify both the ordinate and abscissa; useful when applying fourier transform or integration functions
+
+**Changed:**
+
+* <news item>
+
+**Deprecated:**
+
+* <news item>
+
+**Removed:**
+
+* <news item>
+
+**Fixed:**
+
+* <news item>
+
+**Security:**
+
+* <news item>
diff --git a/src/diffpy/morph/morph_api.py b/src/diffpy/morph/morph_api.py
index f4b06a7..e13f4d0 100644
--- a/src/diffpy/morph/morph_api.py
+++ b/src/diffpy/morph/morph_api.py
@@ -209,7 +209,7 @@ def morph(
             refpars.append("baselineslope")
         elif k == "funcy":
             morph_inst = morph_cls()
-            morph_inst.function = rv_cfg.get("function", None)
+            morph_inst.function = rv_cfg.get("funcy_function", None)
             if morph_inst.function is None:
                 raise ValueError(
                     "Must provide a 'function' when using 'parameters'"
diff --git a/src/diffpy/morph/morph_io.py b/src/diffpy/morph/morph_io.py
index 0c35a3f..fd016af 100644
--- a/src/diffpy/morph/morph_io.py
+++ b/src/diffpy/morph/morph_io.py
@@ -79,32 +79,41 @@ def single_morph_output(
                 rw_pos + idx, (f"squeeze a{idx}", sq_dict[f"a{idx}"])
             )
         mr_copy = dict(morph_results_list)
-    funcy_function = None
-    if "function" in mr_copy:
-        funcy_function = mr_copy.pop("function")
-        print(funcy_function)
-    if "funcy" in mr_copy:
-        fy_dict = mr_copy.pop("funcy")
-        rw_pos = list(mr_copy.keys()).index("Rw")
-        morph_results_list = list(mr_copy.items())
-        for idx, key in enumerate(fy_dict):
-            morph_results_list.insert(
-                rw_pos + idx, (f"funcy {key}", fy_dict[key])
-            )
-        mr_copy = dict(morph_results_list)
+
+    # Handle special inputs (functional remove)
+    func_dicts = {
+        "funcxy": [None, None],
+        "funcx": [None, None],
+        "funcy": [None, None],
+    }
+    for func in func_dicts.keys():
+        if f"{func}_function" in mr_copy:
+            func_dicts[func][0] = mr_copy.pop(f"{func}_function")
+        if func in mr_copy:
+            func_dicts[func][1] = mr_copy.pop(func)
+            rw_pos = list(mr_copy.keys()).index("Rw")
+            morph_results_list = list(mr_copy.items())
+            for idx, key in enumerate(func_dicts[func][1]):
+                morph_results_list.insert(
+                    rw_pos + idx, (f"{func} {key}", func_dicts[func][1][key])
+                )
+            mr_copy = dict(morph_results_list)
+
     # Normal inputs
     morphs_out += "\n".join(
         f"# {key} = {mr_copy[key]:.6f}" for key in mr_copy.keys()
     )
-    # Special inputs (functional)
-    if funcy_function is not None:
-        morphs_in += '# funcy function =\n"""\n'
-        f_code, _ = inspect.getsourcelines(funcy_function)
-        n_leading = len(f_code[0]) - len(f_code[0].lstrip())
-        for idx, f_line in enumerate(f_code):
-            f_code[idx] = f_line[n_leading:]
-        morphs_in += "".join(f_code)
-        morphs_in += '"""\n'
+
+    # Handle special inputs (functional add)
+    for func in func_dicts.keys():
+        if func_dicts[func][0] is not None:
+            morphs_in += f'# {func} function =\n"""\n'
+            f_code, _ = inspect.getsourcelines(func_dicts[func][0])
+            n_leading = len(f_code[0]) - len(f_code[0].lstrip())
+            for idx, f_line in enumerate(f_code):
+                f_code[idx] = f_line[n_leading:]
+            morphs_in += "".join(f_code)
+            morphs_in += '"""\n'
 
     # Printing to terminal
     if stdout_flag:
diff --git a/src/diffpy/morph/morphapp.py b/src/diffpy/morph/morphapp.py
index 3647614..9981e3a 100755
--- a/src/diffpy/morph/morphapp.py
+++ b/src/diffpy/morph/morphapp.py
@@ -520,12 +520,26 @@ def single_morph(
 
     # Python-Specific Morphs
     if pymorphs is not None:
-        # funcy value is a tuple (function,{param_dict})
+        # funcxy/funcx/funcy value is a tuple (function,{param_dict})
+        if "funcxy" in pymorphs:
+            mfxy_function = pymorphs["funcxy"][0]
+            mfxy_params = pymorphs["funcxy"][1]
+            chain.append(morphs.MorphFuncxy())
+            config["funcxy_function"] = mfxy_function
+            config["funcxy"] = mfxy_params
+            refpars.append("funcxy")
+        if "funcx" in pymorphs:
+            mfx_function = pymorphs["funcx"][0]
+            mfx_params = pymorphs["funcx"][1]
+            chain.append(morphs.MorphFuncx())
+            config["funcx_function"] = mfx_function
+            config["funcx"] = mfx_params
+            refpars.append("funcx")
         if "funcy" in pymorphs:
             mfy_function = pymorphs["funcy"][0]
             mfy_params = pymorphs["funcy"][1]
             chain.append(morphs.MorphFuncy())
-            config["function"] = mfy_function
+            config["funcy_function"] = mfy_function
             config["funcy"] = mfy_params
             refpars.append("funcy")
 
@@ -692,6 +706,9 @@ def single_morph(
 
     # FOR FUTURE MAINTAINERS
     # Any new morph should have their input morph parameters updated here
+    # You should also update the IO in morph_io
+    # if you think there requires special handling
+
     # Input morph parameters
     morph_inputs = {
         "scale": scale_in,
@@ -705,11 +722,25 @@ def single_morph(
         for idx, _ in enumerate(squeeze_dict):
             morph_inputs.update({f"squeeze a{idx}": squeeze_dict[f"a{idx}"]})
     if pymorphs is not None:
+        if "funcxy" in pymorphs:
+            for funcxy_param in pymorphs["funcxy"][1].keys():
+                morph_inputs.update(
+                    {
+                        f"funcxy {funcxy_param}": pymorphs["funcxy"][1][
+                            funcxy_param
+                        ]
+                    }
+                )
         if "funcy" in pymorphs:
             for funcy_param in pymorphs["funcy"][1].keys():
                 morph_inputs.update(
                     {f"funcy {funcy_param}": pymorphs["funcy"][1][funcy_param]}
                 )
+        if "funcx" in pymorphs:
+            for funcy_param in pymorphs["funcx"][1].keys():
+                morph_inputs.update(
+                    {f"funcx {funcy_param}": pymorphs["funcx"][1][funcy_param]}
+                )
 
     # Output morph parameters
     morph_results = dict(config.items())
diff --git a/src/diffpy/morph/morphpy.py b/src/diffpy/morph/morphpy.py
index f623815..bac6eee 100644
--- a/src/diffpy/morph/morphpy.py
+++ b/src/diffpy/morph/morphpy.py
@@ -26,7 +26,7 @@ def get_args(parser, params, kwargs):
 
 def __get_morph_opts__(parser, scale, stretch, smear, plot, **kwargs):
     # Check for Python-specific options
-    python_morphs = ["funcy"]
+    python_morphs = ["funcy", "funcx", "funcxy"]
     pymorphs = {}
     for pmorph in python_morphs:
         if pmorph in kwargs:
diff --git a/src/diffpy/morph/morphs/__init__.py b/src/diffpy/morph/morphs/__init__.py
index a1cbdc8..e66f8e3 100644
--- a/src/diffpy/morph/morphs/__init__.py
+++ b/src/diffpy/morph/morphs/__init__.py
@@ -17,6 +17,8 @@
 
 from diffpy.morph.morphs.morph import Morph  # noqa: F401
 from diffpy.morph.morphs.morphchain import MorphChain  # noqa: F401
+from diffpy.morph.morphs.morphfuncx import MorphFuncx
+from diffpy.morph.morphs.morphfuncxy import MorphFuncxy
 from diffpy.morph.morphs.morphfuncy import MorphFuncy
 from diffpy.morph.morphs.morphishape import MorphISphere, MorphISpheroid
 from diffpy.morph.morphs.morphresolution import MorphResolutionDamping
@@ -42,6 +44,8 @@
     MorphShift,
     MorphSqueeze,
     MorphFuncy,
+    MorphFuncx,
+    MorphFuncxy,
 ]
 
 # End of file
diff --git a/src/diffpy/morph/morphs/morphfuncx.py b/src/diffpy/morph/morphs/morphfuncx.py
new file mode 100644
index 0000000..d2c6edc
--- /dev/null
+++ b/src/diffpy/morph/morphs/morphfuncx.py
@@ -0,0 +1,87 @@
+"""Class MorphFuncx -- apply a user-supplied python function to the
+x-axis."""
+
+from diffpy.morph.morphs.morph import LABEL_GR, LABEL_RA, Morph
+
+
+class MorphFuncx(Morph):
+    """Apply a custom function to the x-axis (grid) of the morph
+    function.
+
+    General morph function that applies a user-supplied function to the
+    x-coordinates of morph data to make it align with a target.
+
+    Notice: the function provided must preserve the monotonic
+    increase of the grid.
+    I.e. the function f applied on the grid x must ensure for all
+    indices i<j, f(x[i]) < f(x[j]).
+
+    Configuration Variables
+    -----------------------
+    function: callable
+        The user-supplied function that applies a transformation to the
+        x-coordinates of the data.
+
+    parameters: dict
+        A dictionary of parameters to pass to the function.
+
+    Returns
+    -------
+        A tuple (x_morph_out, y_morph_out, x_target_out, y_target_out)
+        where the target values remain the same and the morph data is
+        transformed according to the user-specified function and parameters
+        The morphed data is returned on the same grid as the unmorphed data
+
+    Example
+    -------
+    Import the funcx morph function:
+
+        >>> from diffpy.morph.morphs.morphfuncx import MorphFuncx
+
+    Define or import the user-supplied transformation function:
+
+        >>> import numpy as np
+        >>> def exp_function(x, y, scale, rate):
+        >>>     return abs(scale) * np.exp(rate * x)
+
+    Note that this transformation is monotonic increasing, so will preserve
+    the monotonic increasing nature of the provided grid.
+
+    Provide initial guess for parameters:
+
+        >>> parameters = {'scale': 1, 'rate': 1}
+
+    Run the funcy morph given input morph array (x_morph, y_morph)and target
+    array (x_target, y_target):
+
+        >>> morph = MorphFuncx()
+        >>> morph.funcx_function = exp_function
+        >>> morph.funcx = parameters
+        >>> x_morph_out, y_morph_out, x_target_out, y_target_out =
+        ... morph.morph(x_morph, y_morph, x_target, y_target)
+
+    To access parameters from the morph instance:
+
+        >>> x_morph_in = morph.x_morph_in
+        >>> y_morph_in = morph.y_morph_in
+        >>> x_target_in = morph.x_target_in
+        >>> y_target_in = morph.y_target_in
+        >>> parameters_out = morph.funcx
+    """
+
+    # Define input output types
+    summary = "Apply a Python function to the x-axis data"
+    xinlabel = LABEL_RA
+    yinlabel = LABEL_GR
+    xoutlabel = LABEL_RA
+    youtlabel = LABEL_GR
+    parnames = ["funcx_function", "funcx"]
+
+    def morph(self, x_morph, y_morph, x_target, y_target):
+        """Apply the user-supplied Python function to the x-coordinates
+        of the morph data."""
+        Morph.morph(self, x_morph, y_morph, x_target, y_target)
+        self.x_morph_out = self.funcx_function(
+            self.x_morph_in, self.y_morph_in, **self.funcx
+        )
+        return self.xyallout
diff --git a/src/diffpy/morph/morphs/morphfuncxy.py b/src/diffpy/morph/morphs/morphfuncxy.py
new file mode 100644
index 0000000..1a802bc
--- /dev/null
+++ b/src/diffpy/morph/morphs/morphfuncxy.py
@@ -0,0 +1,85 @@
+"""Class MorphFuncxy -- apply a user-supplied python function to both
+the x and y axes."""
+
+from diffpy.morph.morphs.morph import LABEL_GR, LABEL_RA, Morph
+
+
+class MorphFuncxy(Morph):
+    """Apply a custom function to the morph function.
+
+    General morph function that applies a user-supplied function to the
+    morph data to make it align with a target.
+
+    This function may modify both the grid (x-axis) and function (y-axis)
+    of the morph data.
+
+    The user-provided function must return a two-column 1D function.
+
+    Configuration Variables
+    -----------------------
+    function: callable
+        The user-supplied function that applies a transformation to the
+        grid (x-axis) and morph function (y-axis).
+
+    parameters: dict
+        A dictionary of parameters to pass to the function.
+
+    Returns
+    -------
+        A tuple (x_morph_out, y_morph_out, x_target_out, y_target_out)
+        where the target values remain the same and the morph data is
+        transformed according to the user-specified function and parameters
+        The morphed data is returned on the same grid as the unmorphed data
+
+    Example (EDIT)
+    -------
+    Import the funcxy morph function:
+
+        >>> from diffpy.morph.morphs.morphfuncxy import MorphFuncxy
+
+    Define or import the user-supplied transformation function:
+
+        >>> import numpy as np
+        >>> def shift_function(x, y, hshift, vshift):
+        >>>     return x + hshift, y + vshift
+
+    Provide initial guess for parameters:
+
+        >>> parameters = {'hshift': 1, 'vshift': 1}
+
+    Run the funcy morph given input morph array (x_morph, y_morph)and target
+    array (x_target, y_target):
+
+        >>> morph = MorphFuncxy()
+        >>> morph.function = shift_function
+        >>> morph.funcy = parameters
+        >>> x_morph_out, y_morph_out, x_target_out, y_target_out =
+        ... morph.morph(x_morph, y_morph, x_target, y_target)
+
+    To access parameters from the morph instance:
+
+        >>> x_morph_in = morph.x_morph_in
+        >>> y_morph_in = morph.y_morph_in
+        >>> x_target_in = morph.x_target_in
+        >>> y_target_in = morph.y_target_in
+        >>> parameters_out = morph.funcxy
+    """
+
+    # Define input output types
+    summary = (
+        "Apply a Python function to the data (y-axis) and data grid (x-axis)"
+    )
+    xinlabel = LABEL_RA
+    yinlabel = LABEL_GR
+    xoutlabel = LABEL_RA
+    youtlabel = LABEL_GR
+    parnames = ["funcxy_function", "funcxy"]
+
+    def morph(self, x_morph, y_morph, x_target, y_target):
+        """Apply the user-supplied Python function to the y-coordinates
+        of the morph data."""
+        Morph.morph(self, x_morph, y_morph, x_target, y_target)
+        self.x_morph_out, self.y_morph_out = self.funcxy_function(
+            self.x_morph_in, self.y_morph_in, **self.funcxy
+        )
+        return self.xyallout
diff --git a/src/diffpy/morph/morphs/morphfuncy.py b/src/diffpy/morph/morphs/morphfuncy.py
index 41716e7..176a024 100644
--- a/src/diffpy/morph/morphs/morphfuncy.py
+++ b/src/diffpy/morph/morphs/morphfuncy.py
@@ -34,6 +34,7 @@ class MorphFuncy(Morph):
 
     Define or import the user-supplied transformation function:
 
+        >>> import numpy as np
         >>> def sine_function(x, y, amplitude, frequency):
         >>>     return amplitude * np.sin(frequency * x) * y
 
@@ -45,7 +46,7 @@ class MorphFuncy(Morph):
     array (x_target, y_target):
 
         >>> morph = MorphFuncy()
-        >>> morph.function = sine_function
+        >>> morph.funcy_function = sine_function
         >>> morph.funcy = parameters
         >>> x_morph_out, y_morph_out, x_target_out, y_target_out =
         ... morph.morph(x_morph, y_morph, x_target, y_target)
@@ -65,13 +66,13 @@ class MorphFuncy(Morph):
     yinlabel = LABEL_GR
     xoutlabel = LABEL_RA
     youtlabel = LABEL_GR
-    parnames = ["function", "funcy"]
+    parnames = ["funcy_function", "funcy"]
 
     def morph(self, x_morph, y_morph, x_target, y_target):
         """Apply the user-supplied Python function to the y-coordinates
         of the morph data."""
         Morph.morph(self, x_morph, y_morph, x_target, y_target)
-        self.y_morph_out = self.function(
+        self.y_morph_out = self.funcy_function(
             self.x_morph_in, self.y_morph_in, **self.funcy
         )
         return self.xyallout
diff --git a/src/diffpy/morph/morphs/morphrgrid.py b/src/diffpy/morph/morphs/morphrgrid.py
index b6f3e54..d0c6ce1 100644
--- a/src/diffpy/morph/morphs/morphrgrid.py
+++ b/src/diffpy/morph/morphs/morphrgrid.py
@@ -51,22 +51,39 @@ class MorphRGrid(Morph):
     youtlabel = LABEL_GR
     parnames = ["rmin", "rmax", "rstep"]
 
+    # Define rmin rmax holders for adaptive x-grid refinement
+    # Without these, the program r-grid can only decrease in interval size
+    rmin_origin = None
+    rmax_origin = None
+    rstep_origin = None
+
     def morph(self, x_morph, y_morph, x_target, y_target):
         """Resample arrays onto specified grid."""
+        if self.rmin is not None:
+            self.rmin_origin = self.rmin
+        if self.rmax is not None:
+            self.rmax_origin = self.rmax
+        if self.rstep is not None:
+            self.rstep_origin = self.rstep
+
         Morph.morph(self, x_morph, y_morph, x_target, y_target)
-        rmininc = max(self.x_target_in[0], self.x_morph_in[0])
-        r_step_target = self.x_target_in[1] - self.x_target_in[0]
-        r_step_morph = self.x_morph_in[1] - self.x_morph_in[0]
+        rmininc = max(min(self.x_target_in), min(self.x_morph_in))
+        r_step_target = (max(self.x_target_in) - min(self.x_target_in)) / (
+            len(self.x_target_in) - 1
+        )
+        r_step_morph = (max(self.x_morph_in) - min(self.x_morph_in)) / (
+            len(self.x_morph_in) - 1
+        )
         rstepinc = max(r_step_target, r_step_morph)
         rmaxinc = min(
-            self.x_target_in[-1] + r_step_target,
-            self.x_morph_in[-1] + r_step_morph,
+            max(self.x_target_in) + r_step_target,
+            max(self.x_morph_in) + r_step_morph,
         )
-        if self.rmin is None or self.rmin < rmininc:
+        if self.rmin_origin is None or self.rmin_origin < rmininc:
             self.rmin = rmininc
-        if self.rmax is None or self.rmax > rmaxinc:
+        if self.rmax_origin is None or self.rmax_origin > rmaxinc:
             self.rmax = rmaxinc
-        if self.rstep is None or self.rstep < rstepinc:
+        if self.rstep_origin is None or self.rstep_origin < rstepinc:
             self.rstep = rstepinc
         # roundoff tolerance for selecting bounds on arrays.
         epsilon = self.rstep / 2
diff --git a/src/diffpy/morph/refine.py b/src/diffpy/morph/refine.py
index c0c5eb6..688e06e 100644
--- a/src/diffpy/morph/refine.py
+++ b/src/diffpy/morph/refine.py
@@ -15,7 +15,7 @@
 """refine -- Refine a morph or morph chain
 """
 
-from numpy import concatenate, dot, exp, ones_like
+from numpy import array, concatenate, dot, exp, ones_like
 from scipy.optimize import leastsq
 from scipy.stats import pearsonr
 
@@ -54,6 +54,10 @@ def __init__(
         self.pars = []
         self.residual = self._residual
         self.flat_to_grouped = {}
+
+        # Padding required for the residual vector to ensure constant length
+        # across the entire morph process
+        self.res_length = None
         return
 
     def _update_chain(self, pvals):
@@ -79,6 +83,26 @@ def _residual(self, pvals):
             self.x_morph, self.y_morph, self.x_target, self.y_target
         )
         rvec = _y_target - _y_morph
+
+        # If first time computing residual
+        if self.res_length is None:
+            self.res_length = len(rvec)
+        # Ensure residual length is constant
+        else:
+            # Padding
+            if len(rvec) < self.res_length:
+                diff_length = self.res_length - len(rvec)
+                rvec = list(rvec)
+                rvec.extend([0] * diff_length)
+                rvec = array(rvec)
+            # Removal
+            # For removal, pass the average RMS
+            # This is fast and easy to compute
+            # For sufficiently functions, this approximation becomes exact
+            elif len(rvec) > self.res_length:
+                avg_rms = sum(rvec**2) / len(rvec)
+                rvec = array([avg_rms for _ in range(self.res_length)])
+
         return rvec
 
     def _pearson(self, pvals):
@@ -147,8 +171,8 @@ def refine(self, *args, **kw):
 
         sol, cov_sol, infodict, emesg, ier = leastsq(
             self.residual,
-            initial,
-            full_output=1,
+            array(initial),
+            full_output=True,
             ftol=self.tolerance,
             xtol=self.tolerance,
         )
diff --git a/tests/test_morph_func.py b/tests/test_morph_func.py
index cb0442d..2fba0a3 100644
--- a/tests/test_morph_func.py
+++ b/tests/test_morph_func.py
@@ -145,7 +145,7 @@ def linear_function(x, y, scale, offset):
     x_target = x_morph.copy()
     y_target = np.sin(x_target) * 2 * x_target + 0.4
     cfg = morph_default_config(funcy={"scale": 1.2, "offset": 0.1})
-    cfg["function"] = linear_function
+    cfg["funcy_function"] = linear_function
     morph_rv = morph(x_morph, y_morph, x_target, y_target, **cfg)
     morphed_cfg = morph_rv["morphed_config"]
     x_morph_out, y_morph_out, x_target_out, y_target_out = morph_rv[
diff --git a/tests/test_morphfuncx.py b/tests/test_morphfuncx.py
new file mode 100644
index 0000000..48444db
--- /dev/null
+++ b/tests/test_morphfuncx.py
@@ -0,0 +1,68 @@
+import numpy as np
+import pytest
+
+from diffpy.morph.morphs.morphfuncx import MorphFuncx
+
+
+def x_exponential_function(x, y, x_amplitude, x_rate):
+    return x_amplitude * np.exp(x_rate * x)
+
+
+def x_linear_function(x, y, x_slope, x_intercept):
+    return x_slope * x + x_intercept
+
+
+def x_cubic_function(x, y, x_amplitude, x_shift):
+    return x_amplitude * (x - x_shift) ** 3
+
+
+def x_arctan_function(x, y, x_amplitude, x_frequency):
+    return x_amplitude * np.arctan(x_frequency * x)
+
+
+funcx_test_suite = [
+    (
+        x_exponential_function,
+        {"x_amplitude": 2, "x_rate": 5},
+        lambda x, y: 2 * np.exp(5 * x),
+    ),
+    (
+        x_linear_function,
+        {"x_slope": 5, "x_intercept": 0.1},
+        lambda x, y: 5 * x + 0.1,
+    ),
+    (
+        x_cubic_function,
+        {"x_amplitude": 2, "x_shift": 5},
+        lambda x, y: 2 * (x - 5) ** 3,
+    ),
+    (
+        x_arctan_function,
+        {"x_amplitude": 4, "x_frequency": 2},
+        lambda x, y: 4 * np.arctan(2 * x),
+    ),
+]
+
+
+@pytest.mark.parametrize(
+    "function, parameters, expected_function",
+    funcx_test_suite,
+)
+def test_funcy(function, parameters, expected_function):
+    x_morph = np.linspace(0, 10, 101)
+    y_morph = np.sin(x_morph)
+    x_target = x_morph.copy()
+    y_target = y_morph.copy()
+    x_morph_expected = expected_function(x_morph, y_morph)
+    y_morph_expected = y_morph
+    morph = MorphFuncx()
+    morph.funcx_function = function
+    morph.funcx = parameters
+    x_morph_actual, y_morph_actual, x_target_actual, y_target_actual = (
+        morph.morph(x_morph, y_morph, x_target, y_target)
+    )
+
+    assert np.allclose(y_morph_actual, y_morph_expected)
+    assert np.allclose(x_morph_actual, x_morph_expected)
+    assert np.allclose(x_target_actual, x_target)
+    assert np.allclose(y_target_actual, y_target)
diff --git a/tests/test_morphfuncxy.py b/tests/test_morphfuncxy.py
new file mode 100644
index 0000000..1de8687
--- /dev/null
+++ b/tests/test_morphfuncxy.py
@@ -0,0 +1,71 @@
+import numpy as np
+import pytest
+
+from diffpy.morph.morphs.morphfuncxy import MorphFuncxy
+
+from .test_morphfuncx import funcx_test_suite
+from .test_morphfuncy import funcy_test_suite
+
+funcxy_test_suite = []
+for entry_y in funcy_test_suite:
+    for entry_x in funcx_test_suite:
+        funcxy_test_suite.append(
+            (
+                entry_x[0],
+                entry_y[0],
+                entry_x[1],
+                entry_y[1],
+                entry_x[2],
+                entry_y[2],
+            )
+        )
+
+
+# FIXME:
+@pytest.mark.parametrize(
+    "funcx_func, funcy_func, funcx_params, funcy_params, "
+    "funcx_lambda, funcy_lambda",
+    funcxy_test_suite,
+)
+def test_funcy(
+    funcx_func,
+    funcy_func,
+    funcx_params,
+    funcy_params,
+    funcx_lambda,
+    funcy_lambda,
+):
+    x_morph = np.linspace(0, 10, 101)
+    y_morph = np.sin(x_morph)
+    x_target = x_morph.copy()
+    y_target = y_morph.copy()
+    x_morph_expected = funcx_lambda(x_morph, y_morph)
+    y_morph_expected = funcy_lambda(x_morph, y_morph)
+
+    funcxy_params = {}
+    funcxy_params.update(funcx_params)
+    funcxy_params.update(funcy_params)
+
+    def funcxy_func(x, y, **funcxy_params):
+        funcx_params = {}
+        funcy_params = {}
+        for param in funcxy_params.keys():
+            if param[:2] == "x_":
+                funcx_params.update({param: funcxy_params[param]})
+            elif param[:2] == "y_":
+                funcy_params.update({param: funcxy_params[param]})
+        return funcx_func(x, y, **funcx_params), funcy_func(
+            x, y, **funcy_params
+        )
+
+    morph = MorphFuncxy()
+    morph.funcxy_function = funcxy_func
+    morph.funcxy = funcxy_params
+    x_morph_actual, y_morph_actual, x_target_actual, y_target_actual = (
+        morph.morph(x_morph, y_morph, x_target, y_target)
+    )
+
+    assert np.allclose(y_morph_actual, y_morph_expected)
+    assert np.allclose(x_morph_actual, x_morph_expected)
+    assert np.allclose(x_target_actual, x_target)
+    assert np.allclose(y_target_actual, y_target)
diff --git a/tests/test_morphfuncy.py b/tests/test_morphfuncy.py
index 31749fd..a775326 100644
--- a/tests/test_morphfuncy.py
+++ b/tests/test_morphfuncy.py
@@ -4,55 +4,58 @@
 from diffpy.morph.morphs.morphfuncy import MorphFuncy
 
 
-def sine_function(x, y, amplitude, frequency):
-    return amplitude * np.sin(frequency * x) * y
+def y_sine_function(x, y, y_amplitude, y_frequency):
+    return y_amplitude * np.sin(y_frequency * x) * y
 
 
-def exponential_decay_function(x, y, amplitude, decay_rate):
-    return amplitude * np.exp(-decay_rate * x) * y
+def y_exponential_decay_function(x, y, y_amplitude, y_decay_rate):
+    return y_amplitude * np.exp(-y_decay_rate * x) * y
 
 
-def gaussian_function(x, y, amplitude, mean, sigma):
-    return amplitude * np.exp(-((x - mean) ** 2) / (2 * sigma**2)) * y
+def y_gaussian_function(x, y, y_amplitude, y_mean, y_sigma):
+    return y_amplitude * np.exp(-((x - y_mean) ** 2) / (2 * y_sigma**2)) * y
 
 
-def polynomial_function(x, y, a, b, c):
-    return (a * x**2 + b * x + c) * y
+def y_polynomial_function(x, y, y_a, y_b, y_c):
+    return (y_a * x**2 + y_b * x + y_c) * y
 
 
-def logarithmic_function(x, y, scale):
-    return scale * np.log(1 + x) * y
+def y_logarithmic_function(x, y, y_scale):
+    return y_scale * np.log(1 + x) * y
+
+
+funcy_test_suite = [
+    (
+        y_sine_function,
+        {"y_amplitude": 2, "y_frequency": 5},
+        lambda x, y: 2 * np.sin(5 * x) * y,
+    ),
+    (
+        y_exponential_decay_function,
+        {"y_amplitude": 5, "y_decay_rate": 0.1},
+        lambda x, y: 5 * np.exp(-0.1 * x) * y,
+    ),
+    (
+        y_gaussian_function,
+        {"y_amplitude": 1, "y_mean": 5, "y_sigma": 1},
+        lambda x, y: np.exp(-((x - 5) ** 2) / (2 * 1**2)) * y,
+    ),
+    (
+        y_polynomial_function,
+        {"y_a": 1, "y_b": 2, "y_c": 0},
+        lambda x, y: (x**2 + 2 * x) * y,
+    ),
+    (
+        y_logarithmic_function,
+        {"y_scale": 0.5},
+        lambda x, y: 0.5 * np.log(1 + x) * y,
+    ),
+]
 
 
 @pytest.mark.parametrize(
     "function, parameters, expected_function",
-    [
-        (
-            sine_function,
-            {"amplitude": 2, "frequency": 5},
-            lambda x, y: 2 * np.sin(5 * x) * y,
-        ),
-        (
-            exponential_decay_function,
-            {"amplitude": 5, "decay_rate": 0.1},
-            lambda x, y: 5 * np.exp(-0.1 * x) * y,
-        ),
-        (
-            gaussian_function,
-            {"amplitude": 1, "mean": 5, "sigma": 1},
-            lambda x, y: np.exp(-((x - 5) ** 2) / (2 * 1**2)) * y,
-        ),
-        (
-            polynomial_function,
-            {"a": 1, "b": 2, "c": 0},
-            lambda x, y: (x**2 + 2 * x) * y,
-        ),
-        (
-            logarithmic_function,
-            {"scale": 0.5},
-            lambda x, y: 0.5 * np.log(1 + x) * y,
-        ),
-    ],
+    funcy_test_suite,
 )
 def test_funcy(function, parameters, expected_function):
     x_morph = np.linspace(0, 10, 101)
@@ -62,7 +65,7 @@ def test_funcy(function, parameters, expected_function):
     x_morph_expected = x_morph
     y_morph_expected = expected_function(x_morph, y_morph)
     morph = MorphFuncy()
-    morph.function = function
+    morph.funcy_function = function
     morph.funcy = parameters
     x_morph_actual, y_morph_actual, x_target_actual, y_target_actual = (
         morph.morph(x_morph, y_morph, x_target, y_target)
diff --git a/tests/test_morphpy.py b/tests/test_morphpy.py
index 4308ca8..621f844 100644
--- a/tests/test_morphpy.py
+++ b/tests/test_morphpy.py
@@ -249,6 +249,68 @@ def gaussian_like_function(x, y, mu):
         assert pytest.approx(abs(morph_info["smear"])) == 4.0
         assert pytest.approx(morph_info["funcy"]["mu"]) == 50.0
 
+    # FIXME:
+    def test_morphfuncx(self, setup_morph):
+        def gaussian(x, mu, sigma):
+            return np.exp(-((x - mu) ** 2) / (2 * sigma**2)) / (
+                sigma * np.sqrt(2 * np.pi)
+            )
+
+        def gaussian_like_function(x, y, mu):
+            return gaussian((x + y) / 2, mu, 3)
+
+        morph_r = np.linspace(0, 100, 1001)
+        morph_gr = np.linspace(0, 100, 1001)
+
+        target_r = np.linspace(0, 100, 1001)
+        target_gr = 0.5 * gaussian(target_r, 50, 5) + 0.05
+
+        morph_info, _ = morph_arrays(
+            np.array([morph_r, morph_gr]).T,
+            np.array([target_r, target_gr]).T,
+            scale=1,
+            smear=3.75,
+            vshift=0.01,
+            funcy=(gaussian_like_function, {"mu": 47.5}),
+            tolerance=1e-12,
+        )
+
+        assert pytest.approx(morph_info["scale"]) == 0.5
+        assert pytest.approx(morph_info["vshift"]) == 0.05
+        assert pytest.approx(abs(morph_info["smear"])) == 4.0
+        assert pytest.approx(morph_info["funcy"]["mu"]) == 50.0
+
+    # FIXME:
+    def test_morphfuncxy(self, setup_morph):
+        def gaussian(x, mu, sigma):
+            return np.exp(-((x - mu) ** 2) / (2 * sigma**2)) / (
+                sigma * np.sqrt(2 * np.pi)
+            )
+
+        def gaussian_like_function(x, y, mu):
+            return gaussian((x + y) / 2, mu, 3)
+
+        morph_r = np.linspace(0, 100, 1001)
+        morph_gr = np.linspace(0, 100, 1001)
+
+        target_r = np.linspace(0, 100, 1001)
+        target_gr = 0.5 * gaussian(target_r, 50, 5) + 0.05
+
+        morph_info, _ = morph_arrays(
+            np.array([morph_r, morph_gr]).T,
+            np.array([target_r, target_gr]).T,
+            scale=1,
+            smear=3.75,
+            vshift=0.01,
+            funcy=(gaussian_like_function, {"mu": 47.5}),
+            tolerance=1e-12,
+        )
+
+        assert pytest.approx(morph_info["scale"]) == 0.5
+        assert pytest.approx(morph_info["vshift"]) == 0.05
+        assert pytest.approx(abs(morph_info["smear"])) == 4.0
+        assert pytest.approx(morph_info["funcy"]["mu"]) == 50.0
+
     def test_morphpy_outputs(self, tmp_path):
         r = np.linspace(0, 1, 11)
         gr = np.linspace(0, 1, 11)
diff --git a/tests/test_refine.py b/tests/test_refine.py
index d05a3db..4cf1bc7 100644
--- a/tests/test_refine.py
+++ b/tests/test_refine.py
@@ -9,6 +9,8 @@
 from diffpy.morph.morph_helpers.transformpdftordf import TransformXtalPDFtoRDF
 from diffpy.morph.morph_helpers.transformrdftopdf import TransformXtalRDFtoPDF
 from diffpy.morph.morphs.morphchain import MorphChain
+from diffpy.morph.morphs.morphfuncx import MorphFuncx
+from diffpy.morph.morphs.morphrgrid import MorphRGrid
 from diffpy.morph.morphs.morphscale import MorphScale
 from diffpy.morph.morphs.morphsmear import MorphSmear
 from diffpy.morph.morphs.morphstretch import MorphStretch
@@ -114,6 +116,71 @@ def test_refine_tolerance(self, setup):
         assert not pytest.approx(config["scale"], stol, stol) == 3.0
         return
 
+    def test_refine_grid_change(self):
+        err = 1e-08
+
+        # First test what occurs when the grid overlap increases
+        # As we shift, the overlap number increases
+        # In this case, overlap goes from 41 -> 51
+        exp_hshift = 1
+        grid1 = numpy.linspace(0, 5, 51)
+        grid2 = numpy.linspace(0 + exp_hshift, 5 + exp_hshift, 51)
+        func1 = numpy.zeros(grid1.shape)
+        func1[(1 < grid1) & (grid1 < 4)] = 1
+        func2 = numpy.zeros(grid2.shape)
+        func2[(1 + exp_hshift < grid2) & (grid2 < 4 + exp_hshift)] = 1
+
+        def shift(x, y, hshift):
+            return x + hshift
+
+        config = {
+            "funcx_function": shift,
+            "funcx": {"hshift": 0},
+            "rmin": 0,
+            "rmax": 7,
+            "rstep": 0.01,
+        }
+
+        mfuncx = MorphFuncx(config)
+        mrgrid = MorphRGrid(config)
+        chain = MorphChain(config, mfuncx, mrgrid)
+        refiner = Refiner(chain, grid1, func1, grid2, func2)
+        refpars = ["funcx"]
+        res = refiner.refine(*refpars)
+
+        assert res < err
+
+        # Second test when the grid overlap decreases
+        # As we stretch, the grid spacing increases
+        # Thus, the overlap number decreases
+        # For this test, overlap goes from 12 -> 10
+        grid1 = numpy.linspace(0, 4, 41)
+        grid2 = numpy.linspace(2, 4, 21)
+        func1 = numpy.zeros(grid1.shape)
+        func1[grid1 <= 2] = 1
+        func1[2 < grid1] = 2
+        func2 = numpy.zeros(grid2.shape) + 1
+
+        def stretch(x, y, stretch):
+            return x * (1 + stretch)
+
+        config = {
+            "funcx_function": stretch,
+            "funcx": {"stretch": 0.7},
+            "rmin": 0,
+            "rmax": 4,
+            "rstep": 0.01,
+        }
+
+        mfuncx = MorphFuncx(config)
+        mrgrid = MorphRGrid(config)
+        chain = MorphChain(config, mfuncx, mrgrid)
+        refiner = Refiner(chain, grid1, func1, grid2, func2)
+        refpars = ["funcx"]
+        res = refiner.refine(*refpars)
+
+        assert res < err
+
 
 # End of class TestRefine