Due to the discussion in this #20054 (comment), it was pointed out that the rotation models unmask the arrays they output. This behavior existed before #18570, meaning that it is outside the scope of #20054. The larger question is, should array masking be preserved by rotations (or even models in general)?

A simple code snippet demonstrating the question is:

import numpy as np
from astropy.modeling.models import RotationSequence3D

x = np.ma.array([1.0, 2.0, 3.0, 4.0], mask=[False, True, False, True])
y = np.ma.array([0.0, 5.0, 6.0, 7.0], mask=[False, False, True, True])
z = np.ma.array([8.0, 9.0, 0.0, 1.0], mask=[False, False, False, True])

print(type(x))
# <class 'numpy.ma.MaskedArray'>

mdl = model([0, 0, 0], axes_order=["x", "y", "z"])

out = mdl(x, y, z)
print(type(out[0])
# <class 'numpy.ndarray'>

Since we input masked arrays into mdl should we expect out to be a tuple of masked arrays? Currently, this is not the case.