diff --git a/.github/workflows/docs.yaml b/.github/workflows/docs.yaml
index b363031..0ac2884 100644
--- a/.github/workflows/docs.yaml
+++ b/.github/workflows/docs.yaml
@@ -10,7 +10,7 @@ jobs:
     runs-on: ${{ matrix.os }}
     strategy:
       matrix:
-        julia-version: [1.6.0]
+        julia-version: [1.11.0]
         julia-arch: [x86]
         os: [ubuntu-latest]
     steps:
diff --git a/docs/Project.toml b/docs/Project.toml
index 00559b3..383cee9 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -11,4 +11,4 @@ SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
-Documenter = "0.25"
\ No newline at end of file
+Documenter = "1.0"
\ No newline at end of file
diff --git a/docs/make.jl b/docs/make.jl
index da21ae5..493086a 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -13,6 +13,7 @@ makedocs(
     repo="https://github.com/byuflowlab/CCBlade.jl/blob/{commit}{path}#L{line}",
     sitename="CCBlade.jl",
     authors="Andrew Ning <aning@byu.edu>",
+    warnonly = Documenter.except(:linkcheck, :footnote),
 )
 
 deploydocs(
diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md
index 3c9b0fe..20d3914 100644
--- a/docs/src/tutorial.md
+++ b/docs/src/tutorial.md
@@ -2,12 +2,16 @@
 
 This starter tutorial walks through the mechanics of running a analysis.  This is designing as a starting point and does not explain every option or consideration.  More specific and advanced usage are described in the [how to guide](howto.md).
 
-We will simulate the APC thin electric 10 x 5 propeller.  The geometry, and wind tunnel data for this propeller is available from [UIUC](https://m-selig.ae.illinois.edu/props/volume-1/propDB-volume-1.html#APC). Let's load CCBlade and a plotting package (I chose to use PyPlot in this example).  
+We will simulate the APC thin electric 10 x 5 propeller.  The geometry, and wind tunnel data for this propeller is available from [UIUC](https://m-selig.ae.illinois.edu/props/volume-1/propDB-volume-1.html#APC). Let's load CCBlade and a plotting package (I chose to use PyPlot in this example).
 
 ```@setup prop
 using CCBlade
 using PyPlot
 ```
+```julia
+using CCBlade
+using PyPlot
+```
 There are two parts of the geometry we need to define: the first are quantities for the whole rotor, and the second part defines properties that vary along the blade radius.  For the rotor, we only need to define the hub radius, tip radius, and the number of blades.
 
 This is a two-bladed propeller, and the 10 in the name (10 x 5) means it has a diameter of 10 inches (and so the tip radius is half of that).  I prefer to do all calculations in metric units so I'll convert it.  From the [geometry table](https://m-selig.ae.illinois.edu/props/volume-1/data/apce_10x5_geom.txt) for this propeller we see that the hub radius is less than ``0.15 R_{tip}``.  It isn't defined exactly, and is less critical, we'll assume ``0.10 R_{tip}`` for this example.
@@ -151,7 +155,7 @@ J &= \frac{V}{n D}
 ```
 
 !!! note
-    Efficiency is set to zero if the thrust is negative (producing drag).  
+    Efficiency is set to zero if the thrust is negative (producing drag).
 
 The code below performs this analysis then plots thrust coefficient, power coefficient, and efficiency as a function of advance ratio as compared to the experimental data.
 
diff --git a/src/CCBlade.jl b/src/CCBlade.jl
index e33e098..bb6c95f 100644
--- a/src/CCBlade.jl
+++ b/src/CCBlade.jl
@@ -26,10 +26,10 @@ include("airfoils.jl")  # all the code related to airfoil data
 # --------- structs -------------
 
 """
-    Rotor(Rhub, Rtip, B; precone=0.0, turbine=false, 
+    Rotor(Rhub, Rtip, B; precone=0.0, turbine=false,
         mach=nothing, re=nothing, rotation=nothing, tip=PrandtlTipHub())
 
-Parameters defining the rotor (apply to all sections).  
+Parameters defining the rotor (apply to all sections).
 
 **Arguments**
 - `Rhub::Float64`: hub radius (along blade length)
@@ -42,8 +42,8 @@ Parameters defining the rotor (apply to all sections).
 - `rotation::RotationCorrection`: correction method for blade rotation
 - `tip::TipCorrection`: correction method for hub/tip loss
 """
-struct Rotor{TF, TI, TB, 
-        T1 <: Union{Nothing, MachCorrection}, T2 <: Union{Nothing, ReCorrection}, 
+struct Rotor{TF, TI, TB,
+        T1 <: Union{Nothing, MachCorrection}, T2 <: Union{Nothing, ReCorrection},
         T3 <: Union{Nothing, RotationCorrection}, T4 <: Union{Nothing, TipCorrection}}
     Rhub::TF
     Rtip::TF
@@ -65,7 +65,7 @@ Rotor(Rhub, Rtip, B; precone=0.0, turbine=false, mach=nothing, re=nothing, rotat
     Section(r, chord, theta, af)
 
 Define sectional properties for one station along rotor
-    
+
 **Arguments**
 - `r::Float64`: radial location along blade
 - `chord::Float64`: corresponding local chord length
@@ -77,7 +77,7 @@ struct Section{TF, TAF}
     chord::TF
     theta::TF
     af::TAF
-end  
+end
 
 # promote to same type, e.g., duals
 function Section(r, chord, theta, af)
@@ -87,15 +87,19 @@ end
 
 
 # convenience function to access fields within an array of structs
-function Base.getproperty(obj::AbstractVector{<:Section}, sym::Symbol)
-    return getfield.(obj, sym)
+function Base.getproperty(obj::Vector{<:Section}, sym::Symbol)
+    if sym in (:ref, :size)
+        return getfield(obj, sym)
+    else
+        return getfield.(obj, sym)
+    end
 end # This is not always type stable b/c we don't know if the return type will be float or af function.
 
 """
     OperatingPoint(Vx, Vy, rho; pitch=0.0, mu=1.0, asound=1.0)
 
-Operation point for a rotor.  
-The x direction is the axial direction, and y direction is the tangential direction in the rotor plane.  
+Operation point for a rotor.
+The x direction is the axial direction, and y direction is the tangential direction in the rotor plane.
 See Documentation for more detail on coordinate systems.
 `Vx` and `Vy` vary radially at same locations as `r` in the rotor definition.
 
@@ -111,7 +115,7 @@ struct OperatingPoint{TF}
     Vx::TF
     Vy::TF
     rho::TF
-    pitch::TF  
+    pitch::TF
     mu::TF
     asound::TF
 end
@@ -124,8 +128,12 @@ OperatingPoint(Vx, Vy, rho, pitch, mu, asound) = OperatingPoint(promote(Vx, Vy,
 OperatingPoint(Vx, Vy, rho; pitch=zero(rho), mu=one(rho), asound=one(rho)) = OperatingPoint(Vx, Vy, rho, pitch, mu, asound)
 
 # convenience function to access fields within an array of structs
-function Base.getproperty(obj::AbstractVector{<:OperatingPoint}, sym::Symbol)
-    return getfield.(obj, sym)
+function Base.getproperty(obj::Vector{<:OperatingPoint}, sym::Symbol)
+    if sym in (:ref, :size)
+        return getfield(obj, sym)
+    else
+        return getfield.(obj, sym)
+    end
 end
 
 
@@ -176,8 +184,12 @@ Outputs(Np, Tp, a, ap, u, v, phi, alpha, W, cl, cd, cn, ct, F, G) = Outputs(prom
 Outputs() = Outputs(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
 
 # convenience function to access fields within an array of structs
-function Base.getproperty(obj::AbstractVector{<:Outputs}, sym::Symbol)
-    return getfield.(obj, sym)
+function Base.getproperty(obj::Vector{<:Outputs}, sym::Symbol)
+    if sym in (:ref, :size)
+        return getfield(obj, sym)
+    else
+        return getfield.(obj, sym)
+    end
 end
 
 # -------------------------------
@@ -195,7 +207,7 @@ function residual_and_outputs(phi, x, p)  #rotor, section, op)
     # unpack inputs
     r, chord, theta, Rhub, Rtip, Vx, Vy, rho, pitch, mu, asound = x  # variables
     af, B, turbine, re_corr, mach_corr, rotation_corr, tip_corr = p  # parameters
-    
+
     # constants
     sigma_p = B*chord/(2.0*pi*r)
     sphi = sin(phi)
@@ -235,7 +247,7 @@ function residual_and_outputs(phi, x, p)  #rotor, section, op)
     # hub/tip loss
     F = 1.0
     if !isnothing(tip_corr)
-        F = tip_correction(tip_corr, r, Rhub, Rtip, phi, B)   
+        F = tip_correction(tip_corr, r, Rhub, Rtip, phi, B)
     end
 
     # sec parameters
@@ -258,7 +270,7 @@ function residual_and_outputs(phi, x, p)  #rotor, section, op)
         a = zero(phi)
         ap = zero(phi)
         R = sign(Vx) + kp
-    
+
     else
 
         if real(phi) < 0
@@ -272,8 +284,8 @@ function residual_and_outputs(phi, x, p)  #rotor, section, op)
         if real(k) >= -2.0/3  # momentum region
             a = k/(1 - k)
 
-        else  # empirical region. Not Buhl's correction but instead uses Buhl with F = 1 then multiplied by F.  
-            # (Buhl(F = 1)*F).  The original method does not force CT -> 0 as f->0.  This can be problematic if 
+        else  # empirical region. Not Buhl's correction but instead uses Buhl with F = 1 then multiplied by F.
+            # (Buhl(F = 1)*F).  The original method does not force CT -> 0 as f->0.  This can be problematic if
             # using CT/CP directly as a design variables.  Suggestion courtesy of Kenneth Lønbæk.
             g1 = 2*k + 1.0/9
             g2 = -2*k - 1.0/3
@@ -305,8 +317,8 @@ function residual_and_outputs(phi, x, p)  #rotor, section, op)
     Np = cn*0.5*rho*W^2*chord
     Tp = ct*0.5*rho*W^2*chord
 
-    # The BEM methodology applies hub/tip losses to the loads rather than to the velocities.  
-    # This is the most common way to implement a BEM, but it means that the raw velocities are misleading 
+    # The BEM methodology applies hub/tip losses to the loads rather than to the velocities.
+    # This is the most common way to implement a BEM, but it means that the raw velocities are misleading
     # as they do not contain any hub/tip loss corrections.
     # To fix this we compute the effective hub/tip losses that would produce the same thrust/torque.
     # In other words:
@@ -519,13 +531,13 @@ function solve(rotor, section, op; npts=10, forcebackwardsearch=false, epsilon_e
             end
             _, outputs = residual_and_outputs(phistar, xv, pv)
             return outputs
-        end    
-    end    
+        end
+    end
 
     # it shouldn't get to this point.  if it does it means no solution was found
     # it will return empty outputs
     # alternatively, one could increase npts and try again
-    
+
     @warn "Invalid data (likely) for this section.  Zero loading assumed."
     return Outputs()
 end
@@ -558,7 +570,7 @@ function simple_op(Vinf, Omega, r, rho; pitch=zero(rho), mu=one(rho), asound=one
         error("You passed in an vector for r, but this function does not accept an vector.\nProbably you intended to use broadcasting")
     end
 
-    Vx = Vinf * cos(precone) 
+    Vx = Vinf * cos(precone)
     Vy = Omega * r * cos(precone)
 
     return OperatingPoint(Vx, Vy, rho, pitch, mu, asound)
@@ -633,7 +645,7 @@ end
 """
     thrusttorque(rotor, sections, outputs::AbstractVector{TO}) where TO
 
-integrate the thrust/torque across the blade, 
+integrate the thrust/torque across the blade,
 including 0 loads at hub/tip, using a trapezoidal rule.
 
 **Arguments**
diff --git a/test/Project.toml b/test/Project.toml
index b3d7195..621909d 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,4 +1,6 @@
 [deps]
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/runtests.jl b/test/runtests.jl
index 3778bb4..6723eef 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -37,7 +37,7 @@ function affunc(alpha, Re, M)
     cl = 0.084*alpha*180/pi
 
     return cl, 0.0
-end 
+end
 
 sections = Section.(r, chord, theta, Ref(affunc))
 
@@ -79,8 +79,8 @@ betavec = 90 .- out.phi*180/pi
 @test isapprox(betavec[6], 84.7113, atol=1e-3)  # using my more converged solution
 
 
-# 
-# 
+#
+#
 
 # idx = 6
 
@@ -143,7 +143,7 @@ function affunc2(alpha, Re, M)
     cd = 0.008 - 0.003*cl + 0.01*cl*cl
 
     return cl, cd
-end 
+end
 
 sections = Section.(r, chord, theta, Ref(affunc2))
 
@@ -159,7 +159,7 @@ qsim = 1e2*[0.803638686218187, 0.806984572453978, 0.809709290183008, 0.811743686
 for i = 1:60
 
     Vinf = float(i)
-    Omega = RPM * pi/30 
+    Omega = RPM * pi/30
 
     ops = simple_op.(Vinf, Omega, r, rho)
 
@@ -209,13 +209,13 @@ function affunc3(alpha, Re, M)
     cd = 0.008 - 0.003*cl + 0.01*cl*cl
 
     return cl, cd
-end 
+end
 
 
 sections = Section.(r, chord, theta, Ref(affunc3))
 
 Vinf = 5.0
-Omega = RPM * pi/30 
+Omega = RPM * pi/30
 ops = simple_op.(Vinf, Omega, r, rho)
 
 out = solve.(Ref(rotor_no_F), sections, ops)
@@ -266,7 +266,7 @@ end
 # empirical data.  Rather they are from the figures used in the documentaiton.  Qualitatively
 # the output is about right.  Minor changes in, for example, the airfoil interpolation method
 # coudl slightly change the outputs.  The main purpose of these tests is to alert us if something
-# significant changes.  
+# significant changes.
 
 Rhub = 1.5
 Rtip = 63.0
@@ -298,7 +298,7 @@ aftypes[8] = AlphaAF("airfoils/NACA64_A17.dat", radians=false)
 # indices correspond to which airfoil is used at which station
 af_idx = [1, 1, 2, 3, 4, 4, 5, 6, 6, 7, 7, 8, 8, 8, 8, 8, 8]
 
-# create airfoil array 
+# create airfoil array
 airfoils = aftypes[af_idx]
 
 sections = Section.(r, chord, theta, airfoils)
@@ -491,7 +491,7 @@ CQ = zeros(nP)
 
 
 for i = 1:nP
-    
+
     op = simple_op.(Vinf, Omega, r, rho, pitch=pitch[i])
     outputs = solve.(Ref(rotor), sections, op)
     T, Q = thrusttorque(rotor, sections, outputs)
@@ -674,7 +674,7 @@ function affunc(alpha, Re, M)
     cd = 0.008 - 0.003*cl + 0.01*cl*cl
 
     return cl, cd
-end 
+end
 
 n = length(r)
 airfoils = fill(affunc, n)
@@ -687,15 +687,54 @@ precone = 0.0
 rho = 1.225
 Vinf = 30.0
 RPM = 2100
-Omega = RPM * pi/30 
+Omega = RPM * pi/30
+
+function ccbladewrapper(x)
+
+    # unpack
+    nall = length(x)
+    nvec = nall - 7
+    n = nvec ÷ 3
+
+    rp = x[1:n]
+    chordp = x[n+1:2*n]
+    thetap = x[2*n+1:3*n]
+    Rhubp = x[3*n+1]
+    Rtipp = x[3*n+2]
+    pitchp = x[3*n+3]
+    preconep = x[3*n+4]
+    Vinfp = x[3*n+5]
+    Omegap = x[3*n+6]
+    rhop = x[3*n+7]
+
+    rotor = Rotor(Rhubp, Rtipp, B; turbine=turbine, precone=preconep)
+    sections = Section.(rp, chordp, thetap, airfoils)
+    ops = simple_op.(Vinfp, Omegap, rp, rhop; pitch=pitchp)
+
+    outputs = solve.(Ref(rotor), sections, ops)
+
+    T, Q = thrusttorque(rotor, sections, outputs)
+
+    return [T; Q]
+end
 
 import ForwardDiff
+
+x = [r; chord; theta; Rhub; Rtip; pitch; precone; Vinf; Omega; rho]
+
+J = ForwardDiff.jacobian(ccbladewrapper, x)
+
+# using BenchmarkTools
+# @btime ForwardDiff.jacobian($ccbladewrapper, $x)
+# original: 584.041 μs (9910 allocations: 1.15 MiB)
+# with ImplicitAD: 323.208 μs (11862 allocations: 747.50 KiB)
+
 import FiniteDiff
 
 J_no_implicitad = nothing
 for implicitad_option in (false, true)
 
-    function ccbladewrapper(x)
+    function ccbladewrapper_fdcheck(x)
         
         # unpack
         nall = length(x)
@@ -726,7 +765,7 @@ for implicitad_option in (false, true)
 
     x = [r; chord; theta; Rhub; Rtip; pitch; precone; Vinf; Omega; rho]
 
-    J = ForwardDiff.jacobian(ccbladewrapper, x)
+    J = ForwardDiff.jacobian(ccbladewrapper_fdcheck, x)
     if !implicitad_option
         J_no_implicitad = J
     # else
@@ -739,11 +778,11 @@ for implicitad_option in (false, true)
     # original: 584.041 μs (9910 allocations: 1.15 MiB) 
     # with ImplicitAD: 323.208 μs (11862 allocations: 747.50 KiB)
 
-    J2 = FiniteDiff.finite_difference_jacobian(ccbladewrapper, x, Val{:central})
+    J2 = FiniteDiff.finite_difference_jacobian(ccbladewrapper_fdcheck, x, Val{:central})
 
     @test maximum(abs.(J - J2)) < 1e-6
 
-    J3 = FiniteDiff.finite_difference_jacobian(ccbladewrapper, x, Val{:complex})
+    J3 = FiniteDiff.finite_difference_jacobian(ccbladewrapper_fdcheck, x, Val{:complex})
 
     @test maximum(abs.(J - J3)) < 3e-11
 
@@ -763,10 +802,10 @@ function ccbladewrapper(x)
 
         cl = 6.2*alpha
         cd = 0.008 - 0.003*cl + 0.01*cl*cl
-    
+
         return cl, cd
-    end 
-    
+    end
+
     # unpack
     nall = length(x)
     nvec = nall - 7
@@ -811,7 +850,7 @@ precone = 0.0
 rho = 1.225
 Vinf = 30.0
 RPM = 2100
-Omega = RPM * pi/30 
+Omega = RPM * pi/30
 
 xvec = [r; chord; theta; Rhub; Rtip; pitch; precone; Vinf; Omega; rho]
 
@@ -831,3 +870,64 @@ end
 
 
 end
+
+using OffsetArrays
+using FillArrays
+
+@testset "esoteric arrays" begin
+    radii = rand(10)
+    chord = rand(length(radii))
+    theta = rand(length(radii))
+    sections = Section.(radii, chord, theta, alpha->(1.0, 1.0))
+
+    Vx = rand(length(radii))
+    Vy = rand(length(radii))
+    rho = rand(length(radii))
+    pitch = rand(length(radii))
+    mu = rand(length(radii))
+    asound = rand(length(radii))
+    ops = OperatingPoint.(Vx, Vy, rho, pitch, mu, asound)
+
+    Np = rand(length(radii))
+    Tp = rand(length(radii))
+    a = rand(length(radii))
+    ap = rand(length(radii))
+    u = rand(length(radii))
+    v = rand(length(radii))
+    phi = rand(length(radii))
+    alpha = rand(length(radii))
+    W = rand(length(radii))
+    cl = rand(length(radii))
+    cd = rand(length(radii))
+    cn = rand(length(radii))
+    ct = rand(length(radii))
+    F = rand(length(radii))
+    G = rand(length(radii))
+    outs = Outputs.(Np, Tp, a, ap, u, v, phi, alpha, W, cl, cd, cn, ct, F, G)
+
+    @testset "OffsetArrays" begin
+        sections_oa = OffsetArray(sections, 0:length(sections)-1)
+        @test sections[1].r ≈ sections_oa[0].r
+
+        ops_oa = OffsetArray(ops, 0:length(ops)-1)
+        @test ops[1].Vx ≈ ops_oa[0].Vx
+
+        outs_oa = OffsetArray(outs, 0:length(outs)-1)
+        @test outs[1].Np ≈ outs_oa[0].Np
+    end
+
+    @testset "FillArrays" begin
+        sections_fill = Fill(sections[1], 3)
+        @test length(sections_fill) == 3
+        sections_fill[3].r ≈ sections[1].r
+
+        ops_fill = Fill(ops[1], 3)
+        @test length(ops_fill) == 3
+        ops_fill[3].Vx ≈ ops[1].Vx
+
+        outs_fill = Fill(outs[1], 3)
+        @test length(outs_fill) == 3
+        outs_fill[3].Np ≈ outs[1].Np
+    end
+
+end