diff --git a/.github/workflows/oas.yml b/.github/workflows/oas.yml
index 199536f74..c11171576 100644
--- a/.github/workflows/oas.yml
+++ b/.github/workflows/oas.yml
@@ -51,6 +51,7 @@ jobs:
     # we need OpenVSP to run vsp tests.
     - name: Install OpenVSP
       run: |
+        sudo apt-get update
         export PYTHON_INCLUDE_DIR=$(python -c "from distutils.sysconfig import get_python_inc; print(get_python_inc())")
         export PYTHON_LIBRARY=$(python -c "import distutils.sysconfig as sysconfig; print(sysconfig.get_config_var('LIBDIR'))")
         export INST_PREFIX=$(python -c "import distutils.sysconfig as sysconfig; print(sysconfig.get_config_var('prefix'))")
diff --git a/openaerostruct/__init__.py b/openaerostruct/__init__.py
index b482efeb4..2614ce9d9 100644
--- a/openaerostruct/__init__.py
+++ b/openaerostruct/__init__.py
@@ -1 +1 @@
-__version__ = "2.6.2"
+__version__ = "2.7.0"
diff --git a/openaerostruct/docs/user_reference/mesh_surface_dict.rst b/openaerostruct/docs/user_reference/mesh_surface_dict.rst
index d1af3aa25..795ad1e91 100644
--- a/openaerostruct/docs/user_reference/mesh_surface_dict.rst
+++ b/openaerostruct/docs/user_reference/mesh_surface_dict.rst
@@ -95,10 +95,10 @@ The surface dict will be provided to Groups, including ``Geometry``, ``AeroPoint
       - np.array([0.1, 5])
       - m
       - B-spline control points for chord distribution. Array convention is the same than ``twist_cp``.
-    * - chord_scaling_pos
+    * - ref_axis_pos
       - 0.25
       - 
-      - Chord position at which the chord scaling factor is applied. 1 is the trailing edge, 0 is the leading edge.
+      - Position of reference axis along the chord about which to apply twist, chord, taper, and span geometry transformations. 1 is the trailing edge, 0 is the leading edge.
 
 .. list-table:: Aerodynamics definitions
     :widths: 20 20 5 55
diff --git a/openaerostruct/examples/rectangular_wing/mphys_opt_chord.py b/openaerostruct/examples/rectangular_wing/mphys_opt_chord.py
index 0b1a73494..87f37ef13 100644
--- a/openaerostruct/examples/rectangular_wing/mphys_opt_chord.py
+++ b/openaerostruct/examples/rectangular_wing/mphys_opt_chord.py
@@ -31,6 +31,7 @@ def setup(self):
             "S_ref_type": "wetted",  # how we compute the wing area,
             # can be 'wetted' or 'projected'
             "twist_cp": np.zeros(3),  # Define twist using 3 B-spline cp's
+            "ref_axis_pos": 0.25,  # Define the reference axis position. 0 is the leading edge, 1 is the trailing edge.
             # Aerodynamic performance of the lifting surface at
             # an angle of attack of 0 (alpha=0).
             # These CL0 and CD0 values are added to the CL and CD
diff --git a/openaerostruct/examples/rectangular_wing/mphys_opt_twist.py b/openaerostruct/examples/rectangular_wing/mphys_opt_twist.py
index a0523e252..63aebc595 100644
--- a/openaerostruct/examples/rectangular_wing/mphys_opt_twist.py
+++ b/openaerostruct/examples/rectangular_wing/mphys_opt_twist.py
@@ -31,6 +31,7 @@ def setup(self):
             "S_ref_type": "wetted",  # how we compute the wing area,
             # can be 'wetted' or 'projected'
             "twist_cp": np.zeros(3),  # Define twist using 3 B-spline cp's
+            "ref_axis_pos": 0.25,  # Define the reference axis position. 0 is the leading edge, 1 is the trailing edge.
             # Aerodynamic performance of the lifting surface at
             # an angle of attack of 0 (alpha=0).
             # These CL0 and CD0 values are added to the CL and CD
diff --git a/openaerostruct/examples/rectangular_wing/opt_chord.py b/openaerostruct/examples/rectangular_wing/opt_chord.py
index 59038a152..f32432e46 100644
--- a/openaerostruct/examples/rectangular_wing/opt_chord.py
+++ b/openaerostruct/examples/rectangular_wing/opt_chord.py
@@ -57,7 +57,7 @@
     "S_ref_type": "projected",  # how we compute the wing area,
     # can be 'wetted' or 'projected'
     "chord_cp": np.ones(3),  # Define chord using 3 B-spline cp's
-    "chord_scaling_pos": 0.25,  # Define the chord scaling position. 0 is the leading edge, 1 is the trailing edge.
+    "ref_axis_pos": 0.25,  # Define the reference axis position. 0 is the leading edge, 1 is the trailing edge.
     # distributed along span
     "mesh": mesh,
     # Aerodynamic performance of the lifting surface at
diff --git a/openaerostruct/examples/rectangular_wing/opt_twist.py b/openaerostruct/examples/rectangular_wing/opt_twist.py
index 983a97bd9..d781cef4e 100644
--- a/openaerostruct/examples/rectangular_wing/opt_twist.py
+++ b/openaerostruct/examples/rectangular_wing/opt_twist.py
@@ -56,6 +56,7 @@
     "S_ref_type": "projected",  # how we compute the wing area,
     # can be 'wetted' or 'projected'
     "twist_cp": np.zeros(3),  # Define twist using 3 B-spline cp's
+    "ref_axis_pos": 0.25,  # Define the reference axis position. 0 is the leading edge, 1 is the trailing edge.
     # distributed along span
     "mesh": mesh,
     # Aerodynamic performance of the lifting surface at
diff --git a/openaerostruct/geometry/geometry_mesh.py b/openaerostruct/geometry/geometry_mesh.py
index 03dcfd821..940861eda 100644
--- a/openaerostruct/geometry/geometry_mesh.py
+++ b/openaerostruct/geometry/geometry_mesh.py
@@ -15,7 +15,6 @@
     ShearZ,
     Rotate,
 )
-import warnings
 
 
 class GeometryMesh(om.Group):
@@ -55,6 +54,11 @@ def initialize(self):
     def setup(self):
         surface = self.options["surface"]
 
+        if "ref_axis_pos" in surface:
+            ref_axis_pos = surface["ref_axis_pos"]
+        else:
+            ref_axis_pos = 0.25  # if no reference axis line is specified : it is the quarter-chord
+
         mesh = surface["mesh"]
         ny = mesh.shape[1]
         mesh_shape = mesh.shape
@@ -73,26 +77,21 @@ def setup(self):
             val = 1.0
             promotes = []
 
-        self.add_subsystem("taper", Taper(val=val, mesh=mesh, symmetry=symmetry), promotes_inputs=promotes)
+        self.add_subsystem(
+            "taper", Taper(val=val, mesh=mesh, symmetry=symmetry, ref_axis_pos=ref_axis_pos), promotes_inputs=promotes
+        )
 
         # 2. Scale X
 
         val = np.ones(ny)
-        chord_scaling_pos = 0.25  # if no scaling position is specified : chord scaling w.r.t quarter of chord
         if "chord_cp" in surface:
             promotes = ["chord"]
-            if "chord_scaling_pos" in surface:
-                chord_scaling_pos = surface["chord_scaling_pos"]
         else:
-            if "chord_scaling_pos" in surface:
-                warnings.warn(
-                    "Chord_scaling_pos has been specified but no chord design variable available", stacklevel=2
-                )
             promotes = []
 
         self.add_subsystem(
             "scale_x",
-            ScaleX(val=val, mesh_shape=mesh_shape, chord_scaling_pos=chord_scaling_pos),
+            ScaleX(val=val, mesh_shape=mesh_shape, ref_axis_pos=ref_axis_pos),
             promotes_inputs=promotes,
         )
 
@@ -124,15 +123,17 @@ def setup(self):
             val = surface["span"]
         else:
             # Compute span. We need .real to make span to avoid OpenMDAO warnings.
-            quarter_chord = 0.25 * mesh[-1, :, :] + 0.75 * mesh[0, :, :]
-            span = max(quarter_chord[:, 1]).real - min(quarter_chord[:, 1]).real
+            ref_axis = ref_axis_pos * mesh[-1, :, :] + (1 - ref_axis_pos) * mesh[0, :, :]
+            span = max(ref_axis[:, 1]).real - min(ref_axis[:, 1]).real
             if symmetry:
                 span *= 2.0
             val = span
             promotes = []
 
         self.add_subsystem(
-            "stretch", Stretch(val=val, mesh_shape=mesh_shape, symmetry=symmetry), promotes_inputs=promotes
+            "stretch",
+            Stretch(val=val, mesh_shape=mesh_shape, symmetry=symmetry, ref_axis_pos=ref_axis_pos),
+            promotes_inputs=promotes,
         )
 
         # 6. Shear Y
@@ -179,7 +180,7 @@ def setup(self):
 
         self.add_subsystem(
             "rotate",
-            Rotate(val=val, mesh_shape=mesh_shape, symmetry=symmetry),
+            Rotate(val=val, mesh_shape=mesh_shape, symmetry=symmetry, ref_axis_pos=ref_axis_pos),
             promotes_inputs=promotes,
             promotes_outputs=["mesh"],
         )
diff --git a/openaerostruct/geometry/geometry_mesh_transformations.py b/openaerostruct/geometry/geometry_mesh_transformations.py
index ea5f8ba5e..514adebe3 100644
--- a/openaerostruct/geometry/geometry_mesh_transformations.py
+++ b/openaerostruct/geometry/geometry_mesh_transformations.py
@@ -9,7 +9,7 @@
 class Taper(om.ExplicitComponent):
     """
     OpenMDAO component that manipulates the mesh by altering the spanwise chord linearly to produce
-    a tapered wing. Note that we apply taper around the quarter-chord line.
+    a tapered wing. Note that we apply taper around the reference axis line which is the quarter-chord by default.
 
     Parameters
     ----------
@@ -31,10 +31,16 @@ def initialize(self):
         self.options.declare(
             "symmetry", default=False, desc="Flag set to true if surface is reflected about y=0 plane."
         )
+        self.options.declare(
+            "ref_axis_pos",
+            default=0.25,
+            desc="Fraction of the chord to use as the reference axis",
+        )
 
     def setup(self):
         mesh = self.options["mesh"]
         val = self.options["val"]
+        self.ref_axis_pos = self.options["ref_axis_pos"]
 
         self.add_input("taper", val=val)
 
@@ -51,8 +57,8 @@ def compute(self, inputs, outputs):
         le = mesh[0]
         te = mesh[-1]
         num_x, num_y, _ = mesh.shape
-        quarter_chord = 0.25 * te + 0.75 * le
-        x = quarter_chord[:, 1]
+        ref_axis = self.ref_axis_pos * te + (1 - self.ref_axis_pos) * le
+        x = ref_axis[:, 1]
         span = x[-1] - x[0]
 
         # If symmetric, solve for the correct taper ratio, which is a linear
@@ -70,7 +76,7 @@ def compute(self, inputs, outputs):
         taper = np.interp(x.real, xp.real, fp.real)
 
         # Modify the mesh based on the taper amount computed per spanwise section
-        outputs["mesh"] = np.einsum("ijk,j->ijk", mesh - quarter_chord, taper) + quarter_chord
+        outputs["mesh"] = np.einsum("ijk,j->ijk", mesh - ref_axis, taper) + ref_axis
 
     def compute_partials(self, inputs, partials):
         mesh = self.options["mesh"]
@@ -81,8 +87,8 @@ def compute_partials(self, inputs, partials):
         le = mesh[0]
         te = mesh[-1]
         num_x, num_y, _ = mesh.shape
-        quarter_chord = 0.25 * te + 0.75 * le
-        x = quarter_chord[:, 1]
+        ref_axis = self.ref_axis_pos * te + (1 - self.ref_axis_pos) * le
+        x = ref_axis[:, 1]
         span = x[-1] - x[0]
 
         # If symmetric, solve for the correct taper ratio, which is a linear
@@ -104,7 +110,7 @@ def compute_partials(self, inputs, partials):
         else:
             dtaper = (1.0 - taper) / (1.0 - taper_ratio)
 
-        partials["mesh", "taper"] = np.einsum("ijk, j->ijk", mesh - quarter_chord, dtaper)
+        partials["mesh", "taper"] = np.einsum("ijk, j->ijk", mesh - ref_axis, dtaper)
 
 
 class ScaleX(om.ExplicitComponent):
@@ -132,15 +138,15 @@ def initialize(self):
         self.options.declare("val", desc="Initial value for chord lengths")
         self.options.declare("mesh_shape", desc="Tuple containing mesh shape (nx, ny).")
         self.options.declare(
-            "chord_scaling_pos",
+            "ref_axis_pos",
             default=0.25,
-            desc="float which indicates the chord fraction where to do the chord scaling. 1. is trailing edge, 0. is leading edge",
+            desc="Fraction of the chord to use as the reference axis",
         )
 
     def setup(self):
         mesh_shape = self.options["mesh_shape"]
         val = self.options["val"]
-        self.chord_scaling_pos = self.options["chord_scaling_pos"]
+        self.ref_axis_pos = self.options["ref_axis_pos"]
         self.add_input("chord", units="m", val=val)
         self.add_input("in_mesh", shape=mesh_shape, units="m")
 
@@ -171,9 +177,9 @@ def compute(self, inputs, outputs):
 
         te = mesh[-1]
         le = mesh[0]
-        chord_pos = self.chord_scaling_pos * te + (1 - self.chord_scaling_pos) * le
+        ref_axis = self.ref_axis_pos * te + (1 - self.ref_axis_pos) * le
 
-        outputs["mesh"] = np.einsum("ijk,j->ijk", mesh - chord_pos, chord_dist) + chord_pos
+        outputs["mesh"] = np.einsum("ijk,j->ijk", mesh - ref_axis, chord_dist) + ref_axis
 
     def compute_partials(self, inputs, partials):
         mesh = inputs["in_mesh"]
@@ -181,9 +187,9 @@ def compute_partials(self, inputs, partials):
 
         te = mesh[-1]
         le = mesh[0]
-        chord_pos = self.chord_scaling_pos * te + (1 - self.chord_scaling_pos) * le
+        ref_axis = self.ref_axis_pos * te + (1 - self.ref_axis_pos) * le
 
-        partials["mesh", "chord"] = (mesh - chord_pos).flatten()
+        partials["mesh", "chord"] = (mesh - ref_axis).flatten()
 
         nx, ny, _ = mesh.shape
         nn = nx * ny * 3
@@ -192,12 +198,12 @@ def compute_partials(self, inputs, partials):
 
         d_qc = (np.einsum("ij,i->ij", np.ones((ny, 3)), 1.0 - chord_dist)).flatten()
         nnq = (nx - 1) * ny * 3
-        partials["mesh", "in_mesh"][nn : nn + nnq] = np.tile(self.chord_scaling_pos * d_qc, nx - 1)
-        partials["mesh", "in_mesh"][nn + nnq :] = np.tile((1 - self.chord_scaling_pos) * d_qc, nx - 1)
+        partials["mesh", "in_mesh"][nn : nn + nnq] = np.tile(self.ref_axis_pos * d_qc, nx - 1)
+        partials["mesh", "in_mesh"][nn + nnq :] = np.tile((1 - self.ref_axis_pos) * d_qc, nx - 1)
 
         nnq = ny * 3
-        partials["mesh", "in_mesh"][nn - nnq : nn] += self.chord_scaling_pos * d_qc
-        partials["mesh", "in_mesh"][:nnq] += (1 - self.chord_scaling_pos) * d_qc
+        partials["mesh", "in_mesh"][nn - nnq : nn] += self.ref_axis_pos * d_qc
+        partials["mesh", "in_mesh"][:nnq] += (1 - self.ref_axis_pos) * d_qc
 
 
 class Sweep(om.ExplicitComponent):
@@ -425,10 +431,16 @@ def initialize(self):
         self.options.declare(
             "symmetry", default=False, desc="Flag set to true if surface is reflected about y=0 plane."
         )
+        self.options.declare(
+            "ref_axis_pos",
+            default=0.25,
+            desc="Fraction of the chord to use as the reference axis",
+        )
 
     def setup(self):
         mesh_shape = self.options["mesh_shape"]
         val = self.options["val"]
+        self.ref_axis_pos = self.options["ref_axis_pos"]
 
         self.add_input("span", val=val, units="m")
         self.add_input("in_mesh", shape=mesh_shape, units="m")
@@ -475,7 +487,7 @@ def compute(self, inputs, outputs):
         # Set the span along the quarter-chord line
         le = mesh[0]
         te = mesh[-1]
-        quarter_chord = 0.25 * te + 0.75 * le
+        ref_axis = self.ref_axis_pos * te + (1 - self.ref_axis_pos) * le
 
         # The user always deals with the full span, so if they input a specific
         # span value and have symmetry enabled, we divide this value by 2.
@@ -484,8 +496,8 @@ def compute(self, inputs, outputs):
 
         # Compute the previous span and determine the scalar needed to reach the
         # desired span
-        prev_span = quarter_chord[-1, 1] - quarter_chord[0, 1]
-        s = quarter_chord[:, 1] / prev_span
+        prev_span = ref_axis[-1, 1] - ref_axis[0, 1]
+        s = ref_axis[:, 1] / prev_span
 
         outputs["mesh"][:] = mesh
         outputs["mesh"][:, :, 1] = s * span
@@ -499,7 +511,7 @@ def compute_partials(self, inputs, partials):
         # Set the span along the quarter-chord line
         le = mesh[0]
         te = mesh[-1]
-        quarter_chord = 0.25 * te + 0.75 * le
+        ref_axis = self.ref_axis_pos * te + (1 - self.ref_axis_pos) * le
 
         # The user always deals with the full span, so if they input a specific
         # span value and have symmetry enabled, we divide this value by 2.
@@ -508,10 +520,10 @@ def compute_partials(self, inputs, partials):
 
         # Compute the previous span and determine the scalar needed to reach the
         # desired span
-        prev_span = quarter_chord[-1, 1] - quarter_chord[0, 1]
-        s = quarter_chord[:, 1] / prev_span
+        prev_span = ref_axis[-1, 1] - ref_axis[0, 1]
+        s = ref_axis[:, 1] / prev_span
 
-        d_prev_span = -quarter_chord[:, 1] / prev_span**2
+        d_prev_span = -ref_axis[:, 1] / prev_span**2
         d_prev_span_qc0 = np.zeros((ny,))
         d_prev_span_qc1 = np.zeros((ny,))
         d_prev_span_qc0[0] = d_prev_span_qc1[-1] = 1.0 / prev_span
@@ -525,17 +537,21 @@ def compute_partials(self, inputs, partials):
         partials["mesh", "in_mesh"][:nn] = 1.0
 
         nn2 = nx * ny
-        partials["mesh", "in_mesh"][nn : nn + nn2] = np.tile(-0.75 * span * (d_prev_span - d_prev_span_qc0), nx)
+        partials["mesh", "in_mesh"][nn : nn + nn2] = np.tile(
+            -(1 - self.ref_axis_pos) * span * (d_prev_span - d_prev_span_qc0), nx
+        )
         nn3 = nn + nn2 * 2
-        partials["mesh", "in_mesh"][nn + nn2 : nn3] = np.tile(0.75 * span * (d_prev_span + d_prev_span_qc1), nx)
+        partials["mesh", "in_mesh"][nn + nn2 : nn3] = np.tile(
+            (1 - self.ref_axis_pos) * span * (d_prev_span + d_prev_span_qc1), nx
+        )
         nn4 = nn3 + nn2
-        partials["mesh", "in_mesh"][nn3:nn4] = np.tile(-0.25 * span * (d_prev_span - d_prev_span_qc0), nx)
+        partials["mesh", "in_mesh"][nn3:nn4] = np.tile(-self.ref_axis_pos * span * (d_prev_span - d_prev_span_qc0), nx)
         nn5 = nn4 + nn2
-        partials["mesh", "in_mesh"][nn4:nn5] = np.tile(0.25 * span * (d_prev_span + d_prev_span_qc1), nx)
+        partials["mesh", "in_mesh"][nn4:nn5] = np.tile(self.ref_axis_pos * span * (d_prev_span + d_prev_span_qc1), nx)
 
         nn6 = nn5 + nx * (ny - 2)
-        partials["mesh", "in_mesh"][nn5:nn6] = 0.75 * span / prev_span
-        partials["mesh", "in_mesh"][nn6:] = 0.25 * span / prev_span
+        partials["mesh", "in_mesh"][nn5:nn6] = (1 - self.ref_axis_pos) * span / prev_span
+        partials["mesh", "in_mesh"][nn6:] = self.ref_axis_pos * span / prev_span
 
 
 class ShearY(om.ExplicitComponent):
@@ -825,10 +841,16 @@ def initialize(self):
             "always be applied perpendicular to the wing (say, in the case of "
             "a winglet).",
         )
+        self.options.declare(
+            "ref_axis_pos",
+            default=0.25,
+            desc="Fraction of the chord to use as the reference axis",
+        )
 
     def setup(self):
         mesh_shape = self.options["mesh_shape"]
         val = self.options["val"]
+        self.ref_axis_pos = self.options["ref_axis_pos"]
 
         self.add_input("twist", val=val, units="deg")
         self.add_input("in_mesh", shape=mesh_shape, units="m")
@@ -904,26 +926,26 @@ def compute(self, inputs, outputs):
 
         te = mesh[-1]
         le = mesh[0]
-        quarter_chord = 0.25 * te + 0.75 * le
+        ref_axis = self.ref_axis_pos * te + (1 - self.ref_axis_pos) * le
 
         _, ny, _ = mesh.shape
 
         if rotate_x:
             # Compute spanwise z displacements along quarter chord
             if symmetry:
-                dz_qc = quarter_chord[:-1, 2] - quarter_chord[1:, 2]
-                dy_qc = quarter_chord[:-1, 1] - quarter_chord[1:, 1]
+                dz_qc = ref_axis[:-1, 2] - ref_axis[1:, 2]
+                dy_qc = ref_axis[:-1, 1] - ref_axis[1:, 1]
                 theta_x = np.arctan(dz_qc / dy_qc)
 
                 # Prepend with 0 so that root is not rotated
                 rad_theta_x = np.append(theta_x, 0.0)
             else:
                 root_index = int((ny - 1) / 2)
-                dz_qc_left = quarter_chord[:root_index, 2] - quarter_chord[1 : root_index + 1, 2]
-                dy_qc_left = quarter_chord[:root_index, 1] - quarter_chord[1 : root_index + 1, 1]
+                dz_qc_left = ref_axis[:root_index, 2] - ref_axis[1 : root_index + 1, 2]
+                dy_qc_left = ref_axis[:root_index, 1] - ref_axis[1 : root_index + 1, 1]
                 theta_x_left = np.arctan(dz_qc_left / dy_qc_left)
-                dz_qc_right = quarter_chord[root_index + 1 :, 2] - quarter_chord[root_index:-1, 2]
-                dy_qc_right = quarter_chord[root_index + 1 :, 1] - quarter_chord[root_index:-1, 1]
+                dz_qc_right = ref_axis[root_index + 1 :, 2] - ref_axis[root_index:-1, 2]
+                dy_qc_right = ref_axis[root_index + 1 :, 1] - ref_axis[root_index:-1, 1]
                 theta_x_right = np.arctan(dz_qc_right / dy_qc_right)
 
                 # Concatenate thetas
@@ -950,7 +972,7 @@ def compute(self, inputs, outputs):
         mats[:, 2, 1] = sin_rtx
         mats[:, 2, 2] = cos_rtx * cos_rty
 
-        outputs["mesh"] = np.einsum("ikj, mij -> mik", mats, mesh - quarter_chord) + quarter_chord
+        outputs["mesh"] = np.einsum("ikj, mij -> mik", mats, mesh - ref_axis) + ref_axis
 
     def compute_partials(self, inputs, partials):
         symmetry = self.options["symmetry"]
@@ -960,15 +982,15 @@ def compute_partials(self, inputs, partials):
 
         te = mesh[-1]
         le = mesh[0]
-        quarter_chord = 0.25 * te + 0.75 * le
+        ref_axis = self.ref_axis_pos * te + (1 - self.ref_axis_pos) * le
 
         nx, ny, _ = mesh.shape
 
         if rotate_x:
             # Compute spanwise z displacements along quarter chord
             if symmetry:
-                dz_qc = quarter_chord[:-1, 2] - quarter_chord[1:, 2]
-                dy_qc = quarter_chord[:-1, 1] - quarter_chord[1:, 1]
+                dz_qc = ref_axis[:-1, 2] - ref_axis[1:, 2]
+                dy_qc = ref_axis[:-1, 1] - ref_axis[1:, 1]
                 theta_x = np.arctan(dz_qc / dy_qc)
 
                 # Prepend with 0 so that root is not rotated
@@ -982,11 +1004,11 @@ def compute_partials(self, inputs, partials):
 
             else:
                 root_index = int((ny - 1) / 2)
-                dz_qc_left = quarter_chord[:root_index, 2] - quarter_chord[1 : root_index + 1, 2]
-                dy_qc_left = quarter_chord[:root_index, 1] - quarter_chord[1 : root_index + 1, 1]
+                dz_qc_left = ref_axis[:root_index, 2] - ref_axis[1 : root_index + 1, 2]
+                dy_qc_left = ref_axis[:root_index, 1] - ref_axis[1 : root_index + 1, 1]
                 theta_x_left = np.arctan(dz_qc_left / dy_qc_left)
-                dz_qc_right = quarter_chord[root_index + 1 :, 2] - quarter_chord[root_index:-1, 2]
-                dy_qc_right = quarter_chord[root_index + 1 :, 1] - quarter_chord[root_index:-1, 1]
+                dz_qc_right = ref_axis[root_index + 1 :, 2] - ref_axis[root_index:-1, 2]
+                dy_qc_right = ref_axis[root_index + 1 :, 1] - ref_axis[root_index:-1, 1]
                 theta_x_right = np.arctan(dz_qc_right / dy_qc_right)
 
                 # Concatenate thetas
@@ -1031,7 +1053,7 @@ def compute_partials(self, inputs, partials):
         dmats_dthy[:, 2, 0] = -cos_rtx * cos_rty * deg2rad
         dmats_dthy[:, 2, 2] = -cos_rtx * sin_rty * deg2rad
 
-        d_dthetay = np.einsum("ikj, mij -> mik", dmats_dthy, mesh - quarter_chord)
+        d_dthetay = np.einsum("ikj, mij -> mik", dmats_dthy, mesh - ref_axis)
         partials["mesh", "twist"] = d_dthetay.flatten()
 
         nn = nx * ny * 9
@@ -1042,8 +1064,8 @@ def compute_partials(self, inputs, partials):
         d_qch = (eye - mats).flatten()
 
         nqc = ny * 9
-        partials["mesh", "in_mesh"][:nqc] += 0.75 * d_qch
-        partials["mesh", "in_mesh"][nn - nqc : nn] += 0.25 * d_qch
+        partials["mesh", "in_mesh"][:nqc] += (1 - self.ref_axis_pos) * d_qch
+        partials["mesh", "in_mesh"][nn - nqc : nn] += self.ref_axis_pos * d_qch
 
         if rotate_x:
             dmats_dthx = np.zeros((ny, 3, 3))
@@ -1054,7 +1076,7 @@ def compute_partials(self, inputs, partials):
             dmats_dthx[:, 2, 1] = cos_rtx
             dmats_dthx[:, 2, 2] = -sin_rtx * cos_rty
 
-            d_dthetax = np.einsum("ikj, mij -> mik", dmats_dthx, mesh - quarter_chord)
+            d_dthetax = np.einsum("ikj, mij -> mik", dmats_dthx, mesh - ref_axis)
             d_dq = np.einsum("ijk, jm -> ijkm", d_dthetax, dthx_dq)
 
             d_dq_flat = d_dq.flatten()
@@ -1062,18 +1084,18 @@ def compute_partials(self, inputs, partials):
             del_n = nn - 9 * ny
             nn2 = nn + del_n
             nn3 = nn2 + del_n
-            partials["mesh", "in_mesh"][nn:nn2] = 0.75 * d_dq_flat[-del_n:]
-            partials["mesh", "in_mesh"][nn2:nn3] = 0.25 * d_dq_flat[:del_n]
+            partials["mesh", "in_mesh"][nn:nn2] = (1 - self.ref_axis_pos) * d_dq_flat[-del_n:]
+            partials["mesh", "in_mesh"][nn2:nn3] = self.ref_axis_pos * d_dq_flat[:del_n]
 
             # Contribution back to main diagonal.
             del_n = 9 * ny
-            partials["mesh", "in_mesh"][:nqc] += 0.75 * d_dq_flat[:del_n]
-            partials["mesh", "in_mesh"][nn - nqc : nn] += 0.25 * d_dq_flat[-del_n:]
+            partials["mesh", "in_mesh"][:nqc] += (1 - self.ref_axis_pos) * d_dq_flat[:del_n]
+            partials["mesh", "in_mesh"][nn - nqc : nn] += self.ref_axis_pos * d_dq_flat[-del_n:]
 
             # Quarter chord direct contribution.
             d_qch_od = np.tile(d_qch.flatten(), nx - 1)
-            partials["mesh", "in_mesh"][nn:nn2] += 0.75 * d_qch_od
-            partials["mesh", "in_mesh"][nn2:nn3] += 0.25 * d_qch_od
+            partials["mesh", "in_mesh"][nn:nn2] += (1 - self.ref_axis_pos) * d_qch_od
+            partials["mesh", "in_mesh"][nn2:nn3] += self.ref_axis_pos * d_qch_od
 
             # off-off diagonal pieces
             if symmetry:
@@ -1081,9 +1103,9 @@ def compute_partials(self, inputs, partials):
 
                 del_n = nn - 9 * nx
                 nn4 = nn3 + del_n
-                partials["mesh", "in_mesh"][nn3:nn4] = -0.75 * d_dq_flat
+                partials["mesh", "in_mesh"][nn3:nn4] = -(1 - self.ref_axis_pos) * d_dq_flat
                 nn5 = nn4 + del_n
-                partials["mesh", "in_mesh"][nn4:nn5] = -0.25 * d_dq_flat
+                partials["mesh", "in_mesh"][nn4:nn5] = -self.ref_axis_pos * d_dq_flat
 
             else:
                 d_dq_flat1 = d_dq[:, :root_index, :, :].flatten()
@@ -1091,10 +1113,10 @@ def compute_partials(self, inputs, partials):
 
                 del_n = nx * root_index * 9
                 nn4 = nn3 + del_n
-                partials["mesh", "in_mesh"][nn3:nn4] = -0.75 * d_dq_flat1
+                partials["mesh", "in_mesh"][nn3:nn4] = -(1 - self.ref_axis_pos) * d_dq_flat1
                 nn5 = nn4 + del_n
-                partials["mesh", "in_mesh"][nn4:nn5] = -0.75 * d_dq_flat2
+                partials["mesh", "in_mesh"][nn4:nn5] = -(1 - self.ref_axis_pos) * d_dq_flat2
                 nn6 = nn5 + del_n
-                partials["mesh", "in_mesh"][nn5:nn6] = -0.25 * d_dq_flat1
+                partials["mesh", "in_mesh"][nn5:nn6] = -self.ref_axis_pos * d_dq_flat1
                 nn7 = nn6 + del_n
-                partials["mesh", "in_mesh"][nn6:nn7] = -0.25 * d_dq_flat2
+                partials["mesh", "in_mesh"][nn6:nn7] = -self.ref_axis_pos * d_dq_flat2
diff --git a/openaerostruct/geometry/tests/test_geometry_mesh_transformations.py b/openaerostruct/geometry/tests/test_geometry_mesh_transformations.py
index e75b1633b..392163e4f 100644
--- a/openaerostruct/geometry/tests/test_geometry_mesh_transformations.py
+++ b/openaerostruct/geometry/tests/test_geometry_mesh_transformations.py
@@ -68,8 +68,9 @@ def test_taper(self):
         group = prob.model
 
         val = self.rng.random(1)
+        ref_axis_pos = self.rng.random(1)
 
-        comp = Taper(val=val, mesh=mesh, symmetry=symmetry)
+        comp = Taper(val=val, mesh=mesh, symmetry=symmetry, ref_axis_pos=ref_axis_pos)
         group.add_subsystem("comp", comp)
 
         prob.setup()
@@ -86,8 +87,9 @@ def test_taper_symmetry(self):
         group = prob.model
 
         val = self.rng.random(1)
+        ref_axis_pos = self.rng.random(1)
 
-        comp = Taper(val=val, mesh=mesh, symmetry=symmetry)
+        comp = Taper(val=val, mesh=mesh, symmetry=symmetry, ref_axis_pos=ref_axis_pos)
         group.add_subsystem("comp", comp)
 
         prob.setup()
@@ -104,8 +106,9 @@ def test_scalex(self):
         group = prob.model
 
         val = self.rng.random(NY)
+        ref_axis_pos = self.rng.random(1)
 
-        comp = ScaleX(val=val, mesh_shape=mesh.shape)
+        comp = ScaleX(val=val, mesh_shape=mesh.shape, ref_axis_pos=ref_axis_pos)
         group.add_subsystem("comp", comp)
 
         prob.setup()
@@ -125,8 +128,9 @@ def test_scalex_symmetry(self):
         group = prob.model
 
         val = self.rng.random(NY)
+        ref_axis_pos = self.rng.random(1)
 
-        comp = ScaleX(val=val, mesh_shape=mesh.shape)
+        comp = ScaleX(val=val, mesh_shape=mesh.shape, ref_axis_pos=ref_axis_pos)
         group.add_subsystem("comp", comp)
 
         prob.setup()
@@ -138,30 +142,7 @@ def test_scalex_symmetry(self):
         check = prob.check_partials(compact_print=True, abs_err_tol=1e-5, rel_err_tol=1e-5)
         assert_check_partials(check, atol=1e-6, rtol=1e-6)
 
-    def test_scalex_chord_scaling_pos_random(self):
-        symmetry = False
-        mesh = get_mesh(symmetry)
-
-        # Test for random values of chord_scaling_pos, check derivatives
-        prob = om.Problem()
-        group = prob.model
-
-        val = self.rng.random(NY)
-        chord_scaling_pos = self.rng.random(1)
-
-        comp = ScaleX(val=val, mesh_shape=mesh.shape, chord_scaling_pos=chord_scaling_pos)
-        group.add_subsystem("comp", comp)
-
-        prob.setup()
-
-        prob["comp.in_mesh"] = mesh
-
-        prob.run_model()
-
-        check = prob.check_partials(compact_print=True, abs_err_tol=1e-5, rel_err_tol=1e-5)
-        assert_check_partials(check, atol=1e-6, rtol=1e-6)
-
-    def test_scalex_chord_scaling_pos_trailing_edge(self):
+    def test_scalex_ref_axis_trailing_edge(self):
         symmetry = True
         mesh = get_mesh(symmetry)
 
@@ -171,7 +152,7 @@ def test_scalex_chord_scaling_pos_trailing_edge(self):
 
         val = self.rng.random(NY)
 
-        comp = ScaleX(val=val, mesh_shape=mesh.shape, chord_scaling_pos=1)
+        comp = ScaleX(val=val, mesh_shape=mesh.shape, ref_axis_pos=1)
         group.add_subsystem("comp", comp)
 
         prob.setup()
@@ -254,8 +235,9 @@ def test_stretch(self):
         group = prob.model
 
         val = self.rng.random(1)
+        ref_axis_pos = self.rng.random(1)
 
-        comp = Stretch(val=val, mesh_shape=mesh.shape, symmetry=symmetry)
+        comp = Stretch(val=val, mesh_shape=mesh.shape, symmetry=symmetry, ref_axis_pos=ref_axis_pos)
         group.add_subsystem("comp", comp)
 
         prob.setup()
@@ -275,8 +257,9 @@ def test_stretch_symmetry(self):
         group = prob.model
 
         val = self.rng.random(1)
+        ref_axis_pos = self.rng.random(1)
 
-        comp = Stretch(val=val, mesh_shape=mesh.shape, symmetry=symmetry)
+        comp = Stretch(val=val, mesh_shape=mesh.shape, symmetry=symmetry, ref_axis_pos=ref_axis_pos)
         group.add_subsystem("comp", comp)
 
         prob.setup()
@@ -380,8 +363,9 @@ def test_rotate(self):
         group = prob.model
 
         val = self.rng.random(NY)
+        ref_axis_pos = self.rng.random(1)
 
-        comp = Rotate(val=val, mesh_shape=mesh.shape, symmetry=symmetry)
+        comp = Rotate(val=val, mesh_shape=mesh.shape, symmetry=symmetry, ref_axis_pos=ref_axis_pos)
         group.add_subsystem("comp", comp)
 
         prob.setup()
@@ -401,8 +385,9 @@ def test_rotate_symmetry(self):
         group = prob.model
 
         val = self.rng.random(NY)
+        ref_axis_pos = self.rng.random(1)
 
-        comp = Rotate(val=val, mesh_shape=mesh.shape, symmetry=symmetry)
+        comp = Rotate(val=val, mesh_shape=mesh.shape, symmetry=symmetry, ref_axis_pos=ref_axis_pos)
         group.add_subsystem("comp", comp)
 
         prob.setup()
@@ -414,6 +399,27 @@ def test_rotate_symmetry(self):
         check = prob.check_partials(compact_print=True, abs_err_tol=1e-5, rel_err_tol=1e-5)
         assert_check_partials(check, atol=1e-6, rtol=1e-6)
 
+    def test_rotate_trailing_edge(self):
+        symmetry = False
+        mesh = get_mesh(symmetry)
+
+        prob = om.Problem()
+        group = prob.model
+
+        val = self.rng.random(NY)
+
+        comp = Rotate(val=val, mesh_shape=mesh.shape, symmetry=symmetry, ref_axis_pos=1)
+        group.add_subsystem("comp", comp)
+
+        prob.setup()
+
+        prob["comp.in_mesh"] = mesh
+
+        prob.run_model()
+
+        # If chord_scaling_pos = 1, TE should not move
+        assert_near_equal(mesh[-1, :, :], prob["comp.mesh"][-1, :, :], tolerance=1e-10)
+
 
 if __name__ == "__main__":
     unittest.main()