diff --git a/Project.toml b/Project.toml
index f8a8f98..9ce60ff 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "RegisterFit"
 uuid = "36121b08-3789-5198-aff2-59a3443d9b59"
-version = "1.0.0"
+version = "1.0.1"
 authors = ["Tim Holy <tim.holy@gmail.com>"]
 
 [deps]
diff --git a/README.md b/README.md
index 2023059..3f70f50 100644
--- a/README.md
+++ b/README.md
@@ -50,7 +50,7 @@ using RegisterFit
 fixed  = zeros(5, 7); fixed[3, 2:6]  .= 1.0
 moving = zeros(7, 5); moving[2:6, 3] .= 1.0
 
-tfms = pat_rotation(fixed, moving)   # 2 candidate AffineMap transforms in 2D
+tfms = principalaxes_rotation(fixed, moving)   # 2 candidate AffineMap transforms in 2D
 # Evaluate each candidate against mismatch data and pick the best
 ```
 
diff --git a/docs/src/api.md b/docs/src/api.md
index 2990c39..73bc6d0 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -4,7 +4,7 @@
 
 ```@docs
 mismatch2affine
-pat_rotation
+principalaxes_rotation
 optimize_per_aperture
 ```
 
@@ -12,6 +12,7 @@ optimize_per_aperture
 
 ```@docs
 qfit
+mms2fit
 mms2fit!
 qbuild
 ```
@@ -20,6 +21,7 @@ qbuild
 
 ```@docs
 uisvalid
+uclamp
 uclamp!
 principalaxes
 ```
diff --git a/docs/src/index.md b/docs/src/index.md
index 1b3f1a5..ca0120c 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -53,7 +53,7 @@ result is an `AffineMap` from
 
 ### Principal Axes Transform
 
-[`pat_rotation`](@ref) provides a complementary rigid-alignment approach that
+[`principalaxes_rotation`](@ref) provides a complementary rigid-alignment approach that
 matches the intensity-weighted covariance ellipsoids of two images. It is useful
 as an initial guess before running the quadratic optimization. Because an ellipsoid
 is symmetric, there are 2 candidate rotations in 2D and 4 in 3D; you must evaluate
@@ -74,7 +74,7 @@ cs, Qs, mmis = mms2fit!(mms, thresh)
 tform = mismatch2affine(mms, thresh, knots)
 
 # 4. (Optional) rigid pre-alignment with PAT
-candidates = pat_rotation(fixed, moving)
+candidates = principalaxes_rotation(fixed, moving)
 # … pick best candidate, then refine with mismatch2affine / RegisterDeformation
 ```
 
@@ -100,6 +100,6 @@ using RegisterFit
 fixed  = zeros(5, 7); fixed[3, 2:6]  .= 1.0   # horizontal bar
 moving = zeros(7, 5); moving[2:6, 3] .= 1.0   # vertical bar
 
-tfms = pat_rotation(fixed, moving)   # 2 candidate AffineMap transforms
+tfms = principalaxes_rotation(fixed, moving)   # 2 candidate AffineMap transforms
 # Select the candidate with the smallest mismatch
 ```
diff --git a/src/RegisterFit.jl b/src/RegisterFit.jl
index fc483e6..fcc2406 100644
--- a/src/RegisterFit.jl
+++ b/src/RegisterFit.jl
@@ -14,14 +14,16 @@ import Base.Cartesian: @nloops, @nexprs, @nref, @nif
 
 export
     mismatch2affine,
+    mms2fit,
     mms2fit!,
     optimize_per_aperture,
-    pat_rotation,
     principalaxes,
+    principalaxes_rotation,
     qbuild,
     qfit,
-    uisvalid,
-    uclamp!
+    uclamp,
+    uclamp!,
+    uisvalid
 
 """
 RegisterFit provides functions that compute affine transformations minimizing
@@ -30,7 +32,7 @@ image registration mismatch, given per-aperture mismatch data from `RegisterMism
 ## Global optimization
 
 - [`mismatch2affine`](@ref): affine transform from mismatch data by least squares
-- [`pat_rotation`](@ref): rigid alignment via a Principal Axes Transformation
+- [`principalaxes_rotation`](@ref): rigid alignment via a Principal Axes Transformation
 - [`optimize_per_aperture`](@ref): naive per-aperture displacement search
 
 ## Utilities
@@ -225,7 +227,7 @@ julia> r[0, 0]
 5.0
 ```
 """
-function qbuild(E0::Real, umin::Vector, Q::Matrix, maxshift)
+function qbuild(E0::Real, umin::AbstractVector, Q::AbstractMatrix, maxshift::Union{AbstractVector,Tuple})
     d = length(maxshift)
     (size(Q, 1) == d && size(Q, 2) == d && length(umin) == d) || error("Size mismatch")
     szout = ((2 * [maxshift...] .+ 1)...,)
@@ -261,7 +263,7 @@ julia> uisvalid([2.5, 0.5], (3, 3))
 false
 ```
 """
-function uisvalid(u::AbstractArray{T}, maxshift) where {T <: Number}
+function uisvalid(u::AbstractArray{T}, maxshift::Union{AbstractVector,Tuple}) where {T <: Number}
     nd = size(u, 1)
     sztail = size(u)[2:end]
     for j in CartesianIndices(sztail), idim in 1:nd
@@ -292,7 +294,7 @@ julia> uclamp!(u, (3, 3))
  -2.49
 ```
 """
-function uclamp!(u::AbstractArray{T}, maxshift) where {T <: Number}
+function uclamp!(u::AbstractArray{T}, maxshift::Union{AbstractVector,Tuple}) where {T <: Number}
     nd = size(u, 1)
     sztail = size(u)[2:end]
     for j in CartesianIndices(sztail), idim in 1:nd
@@ -301,11 +303,19 @@ function uclamp!(u::AbstractArray{T}, maxshift) where {T <: Number}
     return u
 end
 
-function uclamp!(u::AbstractArray{T}, maxshift) where {T <: StaticVector}
+function uclamp!(u::AbstractArray{T}, maxshift::Union{AbstractVector,Tuple}) where {T <: StaticVector}
     uclamp!(reshape(reinterpret(eltype(T), vec(u)), (length(T), size(u)...)), maxshift)
     return u
 end
 
+"""
+    uclamp(u, maxshift)
+
+Non-mutating counterpart to [`uclamp!`](@ref). Returns a clamped copy of `u`;
+the original is unmodified.
+"""
+uclamp(u::AbstractArray, maxshift::Union{AbstractVector,Tuple}) = uclamp!(copy(u), maxshift)
+
 """
     center, cov = principalaxes(img)
 
@@ -379,10 +389,10 @@ end
 end
 
 """
-    tfms = pat_rotation(fixed, moving)
-    tfms = pat_rotation(fixed, moving, SD)
-    tfms = pat_rotation(fixedpa, moving)
-    tfms = pat_rotation(fixedpa, moving, SD)
+    tfms = principalaxes_rotation(fixed, moving)
+    tfms = principalaxes_rotation(fixed, moving, SD)
+    tfms = principalaxes_rotation(fixedpa, moving)
+    tfms = principalaxes_rotation(fixedpa, moving, SD)
 
 Compute the Principal Axes Transform (PAT) aligning the low-order intensity
 moments of two images. `fixed` is the reference image and `moving` is the image
@@ -406,7 +416,7 @@ julia> fixed = zeros(5, 7); fixed[3, 2:6] .= 1.0;   # horizontal bar
 
 julia> moving = zeros(7, 5); moving[2:6, 3] .= 1.0;  # vertical bar
 
-julia> tfms = pat_rotation(fixed, moving);
+julia> tfms = principalaxes_rotation(fixed, moving);
 
 julia> length(tfms)
 2
@@ -417,7 +427,7 @@ julia> tfms[1].linear   # ≈ 90° rotation
  -1.0  0.0
 ```
 """
-function pat_rotation(
+function principalaxes_rotation(
         fixedmoments::Tuple{Vector, Matrix}, moving::AbstractArray,
         SD = Matrix{Float64}(I, ndims(moving), ndims(moving))
     )
@@ -463,8 +473,8 @@ function pat_rotation(
     return tfms
 end
 
-pat_rotation(fixed::AbstractArray, moving::AbstractArray, SD = Matrix{Float64}(I, ndims(fixed), ndims(fixed))) =
-    pat_rotation(principalaxes(fixed), moving, SD)
+principalaxes_rotation(fixed::AbstractArray, moving::AbstractArray, SD = Matrix{Float64}(I, ndims(fixed), ndims(fixed))) =
+    principalaxes_rotation(principalaxes(fixed), moving, SD)
 
 function pat_at(S, SD, c, fmean, mmean)
     Sp = SD \ (S * SD)
@@ -527,6 +537,8 @@ exist, returns `(zero(T), zeros(T, d), zeros(T, d, d))`.
 `maxsep` restricts the fit to shifts satisfying `|u[d] - u0[d]| ≤ maxsep[d]`.
 Setting `opt=false` uses a fast heuristic for `Q` instead of a full nonlinear
 solve, trading accuracy for speed.
+`solver_kwargs` is a `NamedTuple` of keyword arguments forwarded to `NLsolve.nlsolve`
+when `opt=true` (e.g. `solver_kwargs=(iterations=100, ftol=1e-10)`).
 
 # Returns
 - `E0::T` — mismatch value at the fitted minimum
@@ -555,11 +567,11 @@ julia> Q ≈ [1.0 0.0; 0.0 1.0]
 true
 ```
 """
-function qfit(mm::MismatchArray, thresh::Real; maxsep = size(mm), opt::Bool = true)
-    return qfit(mm, thresh, maxsep, opt)
+function qfit(mm::MismatchArray, thresh::Real; maxsep = size(mm), opt::Bool = true, solver_kwargs = (;))
+    return _qfit(mm, thresh, maxsep, opt; solver_kwargs)
 end
 
-function qfit(mm::MismatchArray, thresh::Real, maxsep, opt::Bool)
+function _qfit(mm::MismatchArray, thresh::Real, maxsep, opt::Bool; solver_kwargs = (;))
     T = eltype(eltype(mm))
     threshT = convert(T, thresh)
     d = ndims(mm)
@@ -621,7 +633,7 @@ function qfit(mm::MismatchArray, thresh::Real, maxsep, opt::Bool)
     end
     local results
     function solveql(C, dE, QL, x)
-        return nlsolve((fx, x) -> QLerr!(x, fx, C, dE, similar(QL)), (gx, x) -> QLjac!(x, gx, C, similar(QL)), x)
+        return nlsolve((fx, x) -> QLerr!(x, fx, C, dE, similar(QL)), (gx, x) -> QLjac!(x, gx, C, similar(QL)), x; solver_kwargs...)
     end
     try
         results = solveql(C, dE, QL, x)
@@ -686,6 +698,16 @@ function mms2fit!(mms::AbstractArray{A, N}, thresh) where {A <: MismatchArray, N
     return cs, Qs, mmis
 end
 
+"""
+    mms2fit(mms, thresh)
+
+Non-mutating counterpart to [`mms2fit!`](@ref). Returns the same `(cs, Qs, mmis)` tuple
+without modifying `mms`. Uses `deepcopy` because `interpolate_mm!` writes B-spline
+coefficients into the underlying data of each `MismatchArray` element in-place.
+"""
+mms2fit(mms::AbstractArray{A, N}, thresh) where {A <: MismatchArray, N} =
+    mms2fit!(deepcopy(mms), thresh)
+
 function unpackL!(QL, x)
     d = size(QL, 1)
     indx = 0
diff --git a/test/runtests.jl b/test/runtests.jl
index a429b47..b8955c2 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,6 +1,6 @@
 import RegisterFit
 using Test, Aqua, ExplicitImports, Documenter, CoordinateTransformations, Interpolations, ImageBase, ImageTransformations, LinearAlgebra
-using RegisterCore
+using RegisterCore, StaticArrays
 
 using RegisterUtilities
 
@@ -74,6 +74,10 @@ end
     v1 = 0.3*((-5:5).+1).^2
     v2 = 0.5*((-5:5).-1).^2
     @test A.data ≈ v1.+v2'.+2
+
+    # accepts StaticArrays (as produced by mms2fit!)
+    A2 = RegisterFit.qbuild(2, SVector(-1.0, 1.0), SMatrix{2,2}(0.3, 0.0, 0.0, 0.5), (5,5))
+    @test A2.data ≈ A.data
 end
 
 @testset "uisvalid and uclamp!" begin
@@ -89,6 +93,18 @@ end
     @test abs(u_clamp[1]) < maxshift[1]
     @test abs(u_clamp[2]) < maxshift[2]
     @test abs(u_clamp[3]) < maxshift[3]
+
+    # scalar maxshift must be rejected (Union{AbstractVector,Tuple} annotation)
+    @test_throws MethodError RegisterFit.uisvalid(u_valid, 3)
+    @test_throws MethodError RegisterFit.uclamp!(copy(u_clamp), 3)
+
+    # non-mutating uclamp returns clamped copy, leaves original intact
+    u_orig = reshape([5.0, -6.0, 0.1], 3, 1)
+    u_snapshot = copy(u_orig)
+    u_result = RegisterFit.uclamp(u_orig, maxshift)
+    @test u_result !== u_orig
+    @test u_orig == u_snapshot
+    @test u_result == RegisterFit.uclamp!(copy(u_snapshot), maxshift)
 end
 
 @testset "qfit edge cases" begin
@@ -99,6 +115,16 @@ end
     @test E0 == 0
     @test all(c .== 0)
     @test all(Q .== 0)
+
+    # solver_kwargs forwarded to nlsolve (iterations=1 forces early exit but must not error)
+    denom2 = ones(11, 11)
+    Qmat = rand(Float64, 2, 2); Qmat = Qmat'*Qmat
+    num2 = quadratic(11, 11, [1, -2], Qmat)
+    @test_nowarn RegisterFit.qfit(MismatchArray(num2, denom2), 1e-3; solver_kwargs=(iterations=1,))
+
+    # positional 4-arg form no longer exists (Breaking: CHUNK-007)
+    mm2 = MismatchArray(num2, denom2)
+    @test_throws MethodError RegisterFit.qfit(mm2, 1e-3, size(mm2), false)
 end
 
 @testset "optimize_per_aperture" begin
@@ -128,6 +154,22 @@ end
     @test cs[2, 1] ≈ [0.0,  1.0] atol=1e-10
 end
 
+@testset "mms2fit" begin
+    Q = [1.0 0; 0 1.0]
+    mm1 = MismatchArray(quadratic(5, 5, [1, -1], Q), ones(5, 5))
+    mm2 = MismatchArray(quadratic(5, 5, [0,  1], Q), ones(5, 5))
+    mm3 = MismatchArray(quadratic(5, 5, [-1, 0], Q), ones(5, 5))
+    mm4 = MismatchArray(quadratic(5, 5, [1,  0], Q), ones(5, 5))
+    mms = reshape([mm1, mm2, mm3, mm4], 2, 2)
+    # snapshot raw data before call; interpolate_mm! writes B-spline coefficients in-place
+    data_before = copy(mms[1, 1].data)
+    cs, Qs, mmis = RegisterFit.mms2fit(mms, 0.5)
+    @test mms[1, 1].data == data_before        # original unmodified
+    @test size(cs) == (2, 2)
+    @test cs[1, 1] ≈ [1.0, -1.0] atol=1e-10   # same results as mms2fit!
+    @test cs[2, 1] ≈ [0.0,  1.0] atol=1e-10
+end
+
 @testset "PAT" begin
     # Principal Axes Transformation
     fixed = zeros(10,11)
@@ -137,7 +179,7 @@ end
     moving[3:7,8] .= 1
     moving[2:8,7] .= 1
     fmean, fvar = RegisterFit.principalaxes(fixed)
-    tfm = RegisterFit.pat_rotation((fmean, fvar), moving)
+    tfm = RegisterFit.principalaxes_rotation((fmean, fvar), moving)
     for i = 1:2
         S = tfm[i].linear
         @test abs(S[1,2]) ≈ 1
@@ -147,7 +189,7 @@ end
     end
 
     # Test the array-input convenience wrapper
-    tfms2 = RegisterFit.pat_rotation(fixed, moving)
+    tfms2 = RegisterFit.principalaxes_rotation(fixed, moving)
     @test length(tfms2) == length(tfm)
     for i in eachindex(tfm)
         @test tfms2[i].linear ≈ tfm[i].linear