diff --git a/Project.toml b/Project.toml
index a925491d..e0dda35e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -50,6 +50,7 @@ BenchmarkTools = "1.6"
 CodecXz = "0.7.4"
 Colors = "0"
 ControlPlots = "0.2.4"
+ControlSystemsBase = "1"
 DSP = "0.7.10"
 Dierckx = "0.5.4"
 DiffEqBase = "6.161.0"
@@ -104,6 +105,7 @@ LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
 
 [targets]
-test = ["Test", "BenchmarkTools", "PackageCompiler", "Documenter", "ControlPlots", "Colors", "LaTeXStrings", "StatsBase", "JLD2", "Aqua", "DSP"]
+test = ["Test", "BenchmarkTools", "PackageCompiler", "Documenter", "ControlPlots", "Colors", "LaTeXStrings", "StatsBase", "JLD2", "Aqua", "DSP", "ControlSystemsBase"]
diff --git a/data/prob_1.10_ram_dynamic_3_seg.bin.default.xz b/data/prob_1.10_ram_dynamic_3_seg.bin.default.xz
index f2eca88b..30824cac 100644
Binary files a/data/prob_1.10_ram_dynamic_3_seg.bin.default.xz and b/data/prob_1.10_ram_dynamic_3_seg.bin.default.xz differ
diff --git a/data/prob_1.11_ram_dynamic_3_seg.bin.default.xz b/data/prob_1.11_ram_dynamic_3_seg.bin.default.xz
index f6df56fa..9dfc2f1a 100644
Binary files a/data/prob_1.11_ram_dynamic_3_seg.bin.default.xz and b/data/prob_1.11_ram_dynamic_3_seg.bin.default.xz differ
diff --git a/data/settings_ram.yaml b/data/settings_ram.yaml
index 9c405176..f5d6a9e2 100644
--- a/data/settings_ram.yaml
+++ b/data/settings_ram.yaml
@@ -1,5 +1,5 @@
 system:
-    log_file: "data/log_8700W_8ms" # filename without extension  [replay only]
+    log_file: "data/ram_air_kite"  # filename without extension  [replay only]
                                    #   use / as path delimiter, even on Windows 
     time_lapse:   1.0              # relative replay speed
     sim_time:   100.0              # simulation time             [sim only]
@@ -10,17 +10,13 @@ system:
     fixed_font: ""                 # name or filepath+filename of alternative fixed pitch font
 
 initial:
-    l_tether: 50.0        # initial tether length       [m]
-    elevation: 70.8        # initial elevation angle   [deg]
-    v_reel_out: 0.0        # initial reel out speed    [m/s]
+    l_tethers: [50.0, 50.0, 50.0]  # initial tether length       [m]
+    elevation: 70.8                # initial elevation angle   [deg]
+    v_reel_out: 0.0                # initial reel out speed    [m/s]
 
 solver:
-    abs_tol: 0.01        # absolute tolerance of the DAE solver [m, m/s]
-    rel_tol: 0.01         # relative tolerance of the DAE solver [-]
-    solver: "DFBDF"        # DAE solver, IDA or DImplicitEuler, DFBDF
-    linear_solver: "GMRES" # can be GMRES or Dense or LapackDense (only for IDA)
-    max_order: 4           # maximal order, usually between 3 and 5
-    max_iter:  10000         # max number of iterations of the steady-state-solver
+    abs_tol: 0.01          # absolute tolerance of the DAE solver [m, m/s]
+    rel_tol: 0.01          # relative tolerance of the DAE solver [-]
 
 kite:
     model: "data/ram_air_kite_body.obj"     # 3D model of the kite
@@ -49,7 +45,7 @@ tether:
                            # SK75: 109 to 132 GPa according to datasheet
 
 winch:
-    winch_model: "AsyncMachine" # or TorqueControlledMachine
+    winch_model: "TorqueControlledMachine" # or AsynchMachine
     max_force: 4000        # maximal (nominal) tether force; short overload allowed [N]
     v_ro_max:  8.0         # maximal reel-out speed                          [m/s]
     v_ro_min: -8.0         # minimal reel-out speed (=max reel-in speed)     [m/s]
diff --git a/examples/input_output_function.jl b/examples/input_output_function.jl
deleted file mode 100644
index 7a46875d..00000000
--- a/examples/input_output_function.jl
+++ /dev/null
@@ -1,211 +0,0 @@
-# Copyright (c) 2025 Bart van de Lint
-# SPDX-License-Identifier: MPL-2.0
-
-#=
-This example analyzes the input-output behavior of a ram air kite model by:
-
-1. Creating a ram air kite model and initializing it at 60° elevation
-2. Stabilizing the system
-3. Testing the system response to different control inputs:
-   - Left steering line torque
-   - Right steering line torque
-
-The script plots relationships between steering inputs and resulting angular velocities,
-providing insight into the kite's steering behavior and control characteristics.
-=#
-
-using Timers
-tic()
-@info "Loading packages "
-
-using KiteModels, LinearAlgebra, Statistics, OrdinaryDiffEqCore, OrdinaryDiffEqNonlinearSolve, OrdinaryDiffEqBDF
-using ModelingToolkit: setu, getu, setp, getp
-using ModelingToolkit
-
-PLOT = true
-if PLOT
-    using Pkg
-    if ! ("LaTeXStrings" ∈ keys(Pkg.project().dependencies))
-        using TestEnv; TestEnv.activate()
-    end
-    using ControlPlots, LaTeXStrings
-    import ControlPlots: plot
-end
-toc()
-
-
-include(joinpath(@__DIR__, "plotting.jl"))
-
-# Simulation parameters
-dt = 0.05
-total_time = 10  # Longer simulation to see oscillations
-vsm_interval = 3
-steps = Int(round(total_time / dt))
-
-# Steering parameters
-steering_freq = 1/2  # Hz - full left-right cycle frequency
-steering_magnitude = 10.0      # Magnitude of steering input [Nm]
-
-# Initialize model
-set = load_settings("system_ram.yaml")
-set.segments = 2
-set_values = [-50, 0.0, 0.0]  # Set values of the torques of the three winches. [Nm]
-set.quasi_static = true
-set.physical_model = "ram"
-
-@info "Creating wing, aero, vsm_solver, point_system and s:"
-s = RamAirKite(set)
-s.set.abs_tol = 1e-2
-s.set.rel_tol = 1e-2
-toc()
-
-# init_Q_b_w, R_b_w = KiteModels.measure_to_q(measure)
-# init_kite_pos = init!(s.point_system, s.set, R_b_w)
-# plot(s.point_system, 0.0; zoom=false, front=true)
-
-measure = Measurement()
-
-# Initialize at elevation
-measure.sphere_pos .= deg2rad.([70.0 70.0; 1.0 -1.0])
-KiteModels.init_sim!(s, measure; 
-    adaptive=false, remake=false, reload=true, 
-    solver=FBDF(nlsolve=OrdinaryDiffEqNonlinearSolve.NLNewton(relax=0.8, max_iter=1000))
-)
-OrdinaryDiffEqCore.set_proposed_dt!(s.integrator, dt)
-sys = s.sys
-
-@info "System initialized at:"
-toc()
-
-sys = s.sys
-# # Stabilize system
-s.integrator.ps[sys.stabilize] = true
-for i in 1:10÷dt
-    next_step!(s; dt, vsm_interval=1)
-end
-s.integrator.ps[sys.stabilize] = false
-plot(s, 0.0; zoom=true, front=false)
-
-# Function to step simulation with input u
-function step_with_input_integ(x, u, _, p)
-    (s, set_x, set_ix, set_sx, sx, set_u, get_x, dt) = p
-    set_x(s.integrator, x)
-    set_sx(s.integrator, sx)
-    set_u(s.integrator, u)
-    OrdinaryDiffEqCore.set_t!(s.integrator, 0.0)
-    OrdinaryDiffEqCore.reinit!(s.integrator, s.integrator.u; reinit_dae=false)
-    OrdinaryDiffEqCore.set_proposed_dt!(s.integrator, dt)
-    OrdinaryDiffEqCore.step!(s.integrator, dt)
-    return get_x(s.integrator)
-end
-
-# Get initial state
-x_vec = KiteModels.get_nonstiff_unknowns(s)
-sx_vec = KiteModels.get_stiff_unknowns(s)
-set_ix = setu(s.integrator, Initial.(x_vec)) # set_ix might not be needed anymore if prob is removed
-set_x = setu(s.integrator, x_vec)
-set_sx = setu(s.integrator, sx_vec)
-set_u = setu(s.integrator, collect(sys.set_values))
-get_x = getu(s.integrator, x_vec)
-get_sx = getu(s.integrator, sx_vec)
-sx = get_sx(s.integrator)
-x0 = get_x(s.integrator)
-
-function test_response(s, input_range, input_idx, step_fn, u0, x_idxs=nothing; steps=1)
-    # If no x_idxs specified, default to angular velocities
-    if isnothing(x_idxs)
-        x_idxs = (length(x0)-8):(length(x0)-6) # Assuming ω_b are the last 3 non-stiff states before stiff ones
-        output_size = 3
-    else
-        output_size = length(x_idxs)
-    end
-
-    output = zeros(output_size, length(input_range))
-    total_time = 0.0
-    iter = 0
-
-    for (i, input_val) in enumerate(input_range)
-        u = copy(u0)
-        u[input_idx] += input_val
-        x = copy(x0)
-        set_ix(s.integrator, x)
-        OrdinaryDiffEqCore.reinit!(s.integrator)
-        sx_ = get_sx(s.integrator)
-        for _ in 1:steps # Use _ if i is not used inside the loop
-            p = (s, set_x, set_ix, set_sx, sx_, set_u, get_x, dt) # Pass set_ix even if unused by integ function
-            total_time += @elapsed x = step_fn(x, u, nothing, p)
-            iter += 1
-        end
-        output[:, i] = x[x_idxs]
-    end
-
-    times_rt = dt*iter/total_time
-    @info "Number of steps: $iter, Times realtime: $times_rt, Total time: $total_time"
-    return input_range, output, times_rt
-end
-
-# Add helper function to find state indices
-function find_state_index(x_vec, symbol)
-    # Compare the variables using isequal for symbolic equality
-    idx = findfirst(x -> isequal(x, symbol), x_vec)
-    isnothing(idx) && error("Symbol $symbol not found in state vector")
-    return idx
-end
-
-function plot_input_output_relations(step_fn, suffix)
-    # Find relevant state indices
-    ω_idxs = [find_state_index(x_vec, sys.ω_b[i]) for i in 1:3]
-    twist_idx = find_state_index(sx_vec, sys.free_twist_angle[1])
-
-    # Test ranges
-    steer_range = range(-1, 1, length=100)
-    twist_range = range(-1, 1, length=100)
-
-    # Test steering input vs omega
-    @info "Testing steering input response for $suffix..."
-    _, ω_steer_left, _ = test_response(s, steer_range, 2, step_fn, measure.set_values, ω_idxs)
-    _, ω_steer_right, _ = test_response(s, steer_range, 3, step_fn, measure.set_values, ω_idxs)
-
-    # Test twist angle vs omega
-    @info "Testing twist angle response for $suffix..."
-    function step_with_twist(x, twist_val, _, p)
-        p[5][twist_idx] = twist_val[1]  # Set twist angle directly in stiff state vector copy
-        return step_fn(x, measure.set_values, nothing, p)
-    end
-    _, ω_twist, _ = test_response(s, twist_range, 1, step_with_twist, zeros(3), ω_idxs) # u0 is zeros(3) here, seems ok
-
-    return ω_steer_left, ω_steer_right, ω_twist, steer_range, twist_range
-end
-
-ω_steer_left_integ, ω_steer_right_integ, ω_twist_integ, steer_range, twist_range =
-    plot_input_output_relations(step_with_input_integ, "integ")
-
-# --- Plot combined results ---
-# Simplified data collection
-steering_data = [ω_steer_left_integ[1,:], ω_steer_right_integ[1,:]]
-steering_labels = ["Left Steering (integ)", "Right Steering (integ)"]
-twist_data = [ω_twist_integ[1,:]]
-twist_labels = ["Twist Input (integ)"]
-
-# Steering Plot
-steering_plot = plotx(steer_range,
-    steering_data, # ω_b[1] data
-    [ω_steer_left_integ[2,:], ω_steer_right_integ[2,:]], # ω_b[2] data
-    [ω_steer_left_integ[3,:], ω_steer_right_integ[3,:]]; # ω_b[3] data
-    ylabels=["ω_b[1]", "ω_b[2]", "ω_b[3]"],
-    labels=[steering_labels, steering_labels, steering_labels], # Repeat labels for each subplot
-    fig="Steering Input vs Angular Velocity (Integrator)",
-    xlabel="Steering Input [Nm]")
-
-# Twist Plot
-twist_plot = plotx(rad2deg.(twist_range),
-    twist_data, # ω_b[1] data
-    [ω_twist_integ[2,:]], # ω_b[2] data
-    [ω_twist_integ[3,:]]; # ω_b[3] data
-    ylabels=["ω_b[1]", "ω_b[2]", "ω_b[3]"],
-    labels=[twist_labels, twist_labels, twist_labels], # Repeat labels for each subplot
-    fig="Twist Angle vs Angular Velocity (Integrator)",
-    xlabel="Twist Angle [deg]")
-
-display(steering_plot)
-display(twist_plot)
diff --git a/examples/lin_ram_model.jl b/examples/lin_ram_model.jl
index e5bb3891..4bae048a 100644
--- a/examples/lin_ram_model.jl
+++ b/examples/lin_ram_model.jl
@@ -18,6 +18,7 @@ tic()
 
 using KiteModels, LinearAlgebra, Statistics, OrdinaryDiffEqCore 
 using ModelingToolkit
+using ModelingToolkit: t_nounits as t
 
 PLOT = true
 if PLOT
@@ -25,7 +26,7 @@ if PLOT
     if ! ("LaTeXStrings" ∈ keys(Pkg.project().dependencies))
         using TestEnv; TestEnv.activate()
     end
-    using ControlPlots, LaTeXStrings
+    using ControlPlots, LaTeXStrings, ControlSystemsBase
     import ControlPlots: plot
 end
 toc()
@@ -57,19 +58,16 @@ s.set.abs_tol = 1e-2
 s.set.rel_tol = 1e-2
 toc()
 
-measure = Measurement()
-
-# Define outputs for linearization - angular velocities
-@variables ω_b(ModelingToolkit.t_nounits)[1:3]
+# Define outputs for linearization - heading
+lin_outputs = @variables heading(t)
 
 # Initialize at elevation with linearization outputs
 s.point_system.winches[2].tether_length += 0.2
 s.point_system.winches[3].tether_length += 0.2
-measure.sphere_pos .= deg2rad.([65.0 65.0; 1.0 -1.0])
-KiteModels.init_sim!(s, measure; 
+KiteModels.init_sim!(s; 
     remake=false,
     reload=true,
-    lin_outputs=[ω_b...]  # Specify which outputs to track in linear model
+    lin_outputs  # Specify which outputs to track in linear model
 )
 sys = s.sys
 
@@ -95,207 +93,7 @@ s.integrator.ps[sys.stabilize] = false
 @show norm(A)
 @info "System linearized with matrix dimensions:" A=size(A) B=size(B) C=size(C) D=size(D)
 
-@show rad2deg(s.lin_prob[sys.elevation])
-
-# --- Get operating point values ---
-# Extract state and input at operating point
-u_op = copy(s.integrator[sys.set_values])
-# Create a SysState to capture state values at operating point
-sys_state_op = KiteModels.SysState(s)
-# Also get the direct state vector for linear model
-x_op = copy(s.integrator.u)
-
-# --- Create discretized system matrices ---
-# Function to compute discretized system matrices using matrix exponential method
-function discretize_linear_system(A, B, dt)
-    n = size(A, 1)
-    m = size(B, 2)
-    
-    # Create augmented matrix for simultaneous computation
-    M = [A B; zeros(m, n) zeros(m, m)]
-    
-    # Compute matrix exponential
-    expM = exp(M * dt)
-    
-    # Extract discretized matrices
-    Ad = expM[1:n, 1:n]
-    Bd = expM[1:n, n+1:n+m]
-    
-    return Ad, Bd
-end
-
-# Discretize the continuous-time system for more accurate integration
-@info "Discretizing system matrices..."
-Ad, Bd = discretize_linear_system(A, B, dt)
-
-# Verify operating point stability
-derivatives_at_op = A * zeros(size(A, 1)) + B * zeros(size(B, 2))
-@info "Derivatives at operating point (should be near zero):" norm(derivatives_at_op)
-
-# Create loggers
-logger_nonlinear = Logger(length(s.point_system.points), steps)
-logger_linear = deepcopy(logger_nonlinear)  # Same structure for comparison
-
-# Initialize system state trackers
-sys_state_nonlinear = KiteModels.SysState(s)
-# For the linear model, we'll use a deviation from operating point
-sys_state_linear = deepcopy(sys_state_op)
-
-# --- Prepare the simulation ---
-simulation_time_points = zeros(Float64, steps)
-# Pre-allocate arrays with fixed size instead of growing them dynamically
-input_history = Vector{Vector{Float64}}(undef, steps)
-perturbation_history = Vector{Vector{Float64}}(undef, steps)
-
-# --- Set up linear state tracking ---
-# Create a linear model state vector that starts at zero (deviations from op point)
-x_linear = zeros(length(x_op))
-# Create output vector for linear model
-y_linear = zeros(size(C, 1))
-
-@info "Starting side-by-side simulation..."
-# Begin simulation with fixed number of steps
-try
-    for i in 1:steps
-        global x_linear, sim_time
-        # Calculate current simulation time
-        sim_time = (i-1) * dt
-        simulation_time_points[i] = sim_time
-        
-        # --- Calculate time-varying steering inputs ---
-        # Use sinusoidal input for more realistic testing
-        steering = steering_magnitude * sin(2π * steering_freq * sim_time)
-        
-        # --- Nonlinear system simulation ---
-        # Compute control inputs: base winch force + steering
-        set_values_nonlinear = copy(u_op)
-        set_values_nonlinear .+= [0.0, steering, -steering]
-        input_history[i] = copy(set_values_nonlinear)
-        
-        # Compute perturbation from operating point
-        perturbation = set_values_nonlinear - u_op
-        perturbation_history[i] = copy(perturbation)
-        
-        # Step nonlinear simulation
-        (t_new, _) = next_step!(s, set_values_nonlinear; dt, vsm_interval=vsm_interval)
-        
-        # Log nonlinear state
-        KiteModels.update_sys_state!(sys_state_nonlinear, s)
-        log!(logger_nonlinear, sys_state_nonlinear)
-        
-        # --- Linear system simulation ---
-        # Use discretized state-space model
-        # x[k+1] = Ad*x[k] + Bd*u[k]
-        x_linear = Ad * x_linear + Bd * perturbation
-        
-        # Compute linear system output
-        y_linear .= C * x_linear + D * perturbation
-        sys_state_linear.turn_rates = sys_state_op.turn_rates .+ y_linear
-        
-        # Log linear state
-        log!(logger_linear, sys_state_linear)
-    end
-catch e
-    if isa(e, AssertionError)
-        @show i
-        println(e)
-    else
-        rethrow(e)
-    end
-end
-
-@info "Simulation completed:"
-toc()
-
-# --- Save logs ---
-save_log(logger_nonlinear, "nonlinear_model")
-save_log(logger_linear, "linear_model")
-lg_nonlinear = load_log("nonlinear_model")
-lg_linear = load_log("linear_model")
-sl_nonlinear = lg_nonlinear.syslog
-sl_linear = lg_linear.syslog
-
-# --- Plot comparison results ---
-# Extract necessary data
-turn_rates_nonlinear_deg = rad2deg.(hcat(sl_nonlinear.turn_rates...))
-turn_rates_linear_deg = rad2deg.(hcat(sl_linear.turn_rates...))
-
-# Format input data for plotting
-steering_inputs = [inputs[2] - u_op[2] for inputs in input_history]
-
-# Create comparison plots using plotx
-t_plot = sl_nonlinear.time
-
-# Prepare the data in the format expected by plotx
-ω_x_comparison = [turn_rates_nonlinear_deg[1,:], turn_rates_linear_deg[1,:]]
-ω_y_comparison = [turn_rates_nonlinear_deg[2,:], turn_rates_linear_deg[2,:]]
-ω_z_comparison = [turn_rates_nonlinear_deg[3,:], turn_rates_linear_deg[3,:]]
-steering_series = [steering_inputs]
-
-# Create the comparison plot
-p_comparison = plotx(t_plot,
-    ω_x_comparison,
-    ω_y_comparison,
-    ω_z_comparison,
-    steering_series;
-    ylabels=["ω_x [°/s]", "ω_y [°/s]", "ω_z [°/s]", "Steering [Nm]"],
-    labels=[
-        ["Nonlinear", "Linear"],
-        ["Nonlinear", "Linear"],
-        ["Nonlinear", "Linear"],
-        ["Input"]
-    ],
-    fig="Linear vs Nonlinear Model Comparison")
-display(p_comparison)
-
-# --- Calculate error metrics ---
-# Function to calculate normalized RMSE
-function calculate_nrmse(actual, predicted)
-    # Trim to same length if needed
-    min_length = min(length(actual), length(predicted))
-    # Calculate RMSE
-    rmse = sqrt(sum((actual[1:min_length] - predicted[1:min_length]).^2) / min_length)
-    # Normalize by range of actual values
-    range_actual = maximum(actual[1:min_length]) - minimum(actual[1:min_length])
-    # Avoid division by zero
-    if abs(range_actual) < 1e-10
-        return rmse  # Return unnormalized if range is too small
-    else
-        return rmse / range_actual
-    end
-end
-
-# Calculate both average L2 norm and NRMSE
-len = turn_rates_nonlinear_deg[1,:]
-error_ω_x = norm(turn_rates_nonlinear_deg[1,:] - turn_rates_linear_deg[1,:]) / len
-error_ω_y = norm(turn_rates_nonlinear_deg[2,:] - turn_rates_linear_deg[2,:]) / len
-error_ω_z = norm(turn_rates_nonlinear_deg[3,:] - turn_rates_linear_deg[3,:]) / len
-
-nrmse_ω_x = calculate_nrmse(turn_rates_nonlinear_deg[1,:], turn_rates_linear_deg[1,:])
-nrmse_ω_y = calculate_nrmse(turn_rates_nonlinear_deg[2,:], turn_rates_linear_deg[2,:])
-nrmse_ω_z = calculate_nrmse(turn_rates_nonlinear_deg[3,:], turn_rates_linear_deg[3,:])
-
-@info "Error metrics (average L2 norm):" error_ω_x error_ω_y error_ω_z
-@info "Normalized RMSE (lower is better):" nrmse_ω_x nrmse_ω_y nrmse_ω_z
-
-# --- Plot error over time ---
-# Calculate difference between linear and nonlinear models
-diff_ω_x = turn_rates_linear_deg[1,:] - turn_rates_nonlinear_deg[1,:]
-diff_ω_y = turn_rates_linear_deg[2,:] - turn_rates_nonlinear_deg[2,:]
-diff_ω_z = turn_rates_linear_deg[3,:] - turn_rates_nonlinear_deg[3,:]
-
-# Plot the differences
-p_error = plotx(t_plot,
-    [diff_ω_x],
-    [diff_ω_y],
-    [diff_ω_z],
-    [steering_inputs];
-    ylabels=["ω_x error [°/s]", "ω_y error [°/s]", "ω_z error [°/s]", "Steering [Nm]"],
-    labels=[
-        ["Error"],
-        ["Error"],
-        ["Error"],
-        ["Input"]
-    ],
-    fig="Linear vs Nonlinear Model Error")
-display(p_error)
\ No newline at end of file
+sys = ss(A,B,C,D)
+bode_plot(sys[1,1]; from=1e-4)
+bode_plot(sys[1,2]; from=1e-4)
+bode_plot(sys[1,3]; from=1e-4)
diff --git a/examples/ram_air_kite.jl b/examples/ram_air_kite.jl
index 99dc7e38..34042871 100644
--- a/examples/ram_air_kite.jl
+++ b/examples/ram_air_kite.jl
@@ -47,17 +47,14 @@ s.set.abs_tol = 1e-2
 s.set.rel_tol = 1e-2
 toc()
 
-measure = Measurement()
-measure.sphere_pos .= deg2rad.([70.0 70.0; 1.0 -1.0])
-
-# init_Q_b_w, R_b_w = KiteModels.measure_to_q(measure)
+# init_Q_b_w, R_b_w = KiteModels.initial_orient(s.set)
 # init_kite_pos = init!(s.point_system, s.set, R_b_w, init_Q_b_w)
 # plot(s.point_system, 0.0; zoom=false, front=false)
 
 # Initialize at elevation
 s.point_system.winches[2].tether_length += 0.2
 s.point_system.winches[3].tether_length += 0.2
-KiteModels.init_sim!(s, measure; remake=false, reload=true)
+KiteModels.init_sim!(s; remake=false, reload=true)
 sys = s.sys
 
 @info "System initialized at:"
diff --git a/examples/steering_test_1p.jl b/examples/steering_test_1p.jl
index ec9d5aac..080f48b3 100644
--- a/examples/steering_test_1p.jl
+++ b/examples/steering_test_1p.jl
@@ -12,7 +12,7 @@ set.rel_tol=0.000001
 set.cs_4p = 1.0
 set.v_wind = 12
 set.elevation = 69.4
-set.l_tether = 200
+set.l_tethers[1] = 200
 set.depower = set.depower_offset # fully powered kite
 # set.kcu_mass = 8.4
 # set.v_steering = 0.2*6
diff --git a/examples/steering_test_4p.jl b/examples/steering_test_4p.jl
index b93382df..7c0ccfe6 100644
--- a/examples/steering_test_4p.jl
+++ b/examples/steering_test_4p.jl
@@ -15,7 +15,7 @@ set.abs_tol=0.00006
 set.rel_tol=0.000001
 set.v_wind = 10
 set.elevation = 69.4
-set.l_tether = 200
+set.l_tethers[1] = 200
 set.depower = set.depower_offset # fully powered kite
 # set.kcu_mass = 8.4
 # set.v_steering = 0.2*6
diff --git a/src/KiteModels.jl b/src/KiteModels.jl
index ee0a988a..e841fdee 100644
--- a/src/KiteModels.jl
+++ b/src/KiteModels.jl
@@ -36,7 +36,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D
 using ADTypes: AutoFiniteDiff
 import ModelingToolkit.SciMLBase: successful_retcode
 
-export KPS3, KPS4, RamAirKite, KVec3, SimFloat, Measurement, PointMassSystem, ProfileLaw, EXP, LOG, EXPLOG     # constants and types
+export KPS3, KPS4, RamAirKite, KVec3, SimFloat, PointMassSystem, ProfileLaw, EXP, LOG, EXPLOG     # constants and types
 export calc_set_cl_cd!, copy_examples, copy_bin, update_sys_state!                            # helper functions
 export clear!, find_steady_state!, residual!                                                  # low level workers
 export init_sim!, init!, reinit!, next_step!, init_pos_vel                                    # high level workers
diff --git a/src/mtk_model.jl b/src/mtk_model.jl
index df070926..c3c34e66 100644
--- a/src/mtk_model.jl
+++ b/src/mtk_model.jl
@@ -650,14 +650,13 @@ function calc_R_t_w(elevation, azimuth)
 end
 
 """
-    scalar_eqs!(s, eqs, measure; R_b_w, wind_vec_gnd, va_kite_b, kite_pos, kite_vel, kite_acc, twist_angle, twist_ω)
+    scalar_eqs!(s, eqs; R_b_w, wind_vec_gnd, va_kite_b, kite_pos, kite_vel, kite_acc, twist_angle, twist_ω)
 
 Generate equations for scalar quantities like elevation, azimuth, heading and course angles.
     
     # Arguments
     - `s::RamAirKite`: The kite system state
     - `eqs`: Current system equations
-    - `measure::Measurement`: Current measurement data
     - `R_b_w`: Body to world rotation matrix
     - `wind_vec_gnd`: Ground wind vector
     - `va_kite_b`: Apparent wind velocity in body frame
@@ -673,9 +672,9 @@ Generate equations for scalar quantities like elevation, azimuth, heading and co
     - Course angle
     - Angular velocities and accelerations
     """
-function scalar_eqs!(s, eqs, measure; R_b_w, wind_vec_gnd, va_kite_b, kite_pos, kite_vel, kite_acc, twist_angle, twist_ω, ω_b, α_b)
+function scalar_eqs!(s, eqs; R_b_w, wind_vec_gnd, va_kite_b, kite_pos, kite_vel, kite_acc, twist_angle, twist_ω, ω_b, α_b)
     @parameters wind_scale_gnd = s.set.v_wind
-    @parameters measured_wind_dir_gnd = measure.wind_dir_gnd
+    @parameters upwind_dir = deg2rad(s.set.upwind_dir)
     @variables begin
         e_x(t)[1:3]
         e_y(t)[1:3]
@@ -688,7 +687,7 @@ function scalar_eqs!(s, eqs, measure; R_b_w, wind_vec_gnd, va_kite_b, kite_pos,
         e_x     ~ R_b_w[:,1]
         e_y     ~ R_b_w[:,2]
         e_z     ~ R_b_w[:,3]
-        wind_vec_gnd ~ wind_scale_gnd * rotate_around_z([1, 0, 0], measured_wind_dir_gnd)
+        wind_vec_gnd ~ wind_scale_gnd * rotate_around_z([0, -1, 0], -upwind_dir)
         wind_vel_kite ~ AtmosphericModels.calc_wind_factor(s.am, kite_pos[3], s.set.profile_law) * wind_vec_gnd
         va_kite ~ wind_vel_kite - kite_vel
         va_kite_b ~ R_b_w' * va_kite
@@ -817,19 +816,18 @@ function linear_vsm_eqs!(s, eqs, guesses; aero_force_b, aero_moment_b, group_aer
     return eqs, guesses
 end
 
-function create_sys!(s::RamAirKite, system::PointMassSystem, measure::Measurement; init_va_b)
+function create_sys!(s::RamAirKite, system::PointMassSystem; init_va_b)
     eqs = []
     defaults = Pair{Num, Real}[]
     guesses = Pair{Num, Real}[]
 
-    init_set_values = measure.set_values
     @parameters begin
         stabilize = false
         fix_nonstiff = false
     end
     @variables begin
         # potential differential variables
-        set_values(t)[eachindex(system.winches)] = init_set_values
+        set_values(t)[eachindex(system.winches)] = zeros(length(system.winches))
         kite_pos(t)[1:3] # xyz pos of kite in world frame
         kite_vel(t)[1:3]
         kite_acc(t)[1:3]
@@ -855,7 +853,7 @@ function create_sys!(s::RamAirKite, system::PointMassSystem, measure::Measuremen
     eqs, guesses = linear_vsm_eqs!(s, eqs, guesses; aero_force_b, aero_moment_b, group_aero_moment, init_va_b, twist_angle, va_kite_b, ω_b)
     eqs, defaults = diff_eqs!(s, eqs, defaults; tether_kite_force, tether_kite_moment, aero_force_b, aero_moment_b, 
         ω_b, α_b, R_b_w, kite_pos, kite_vel, kite_acc, stabilize, fix_nonstiff)
-    eqs = scalar_eqs!(s, eqs, measure; R_b_w, wind_vec_gnd, va_kite_b, kite_pos, kite_vel, kite_acc, twist_angle, twist_ω, ω_b, α_b)
+    eqs = scalar_eqs!(s, eqs; R_b_w, wind_vec_gnd, va_kite_b, kite_pos, kite_vel, kite_acc, twist_angle, twist_ω, ω_b, α_b)
     
     # te_I = (1/3 * (s.set.mass/8) * te_length^2)
     # # -damping / I * ω = α_damping
diff --git a/src/point_mass_system.jl b/src/point_mass_system.jl
index 9dd82b7d..e391f29a 100644
--- a/src/point_mass_system.jl
+++ b/src/point_mass_system.jl
@@ -493,58 +493,14 @@ function init!(system::PointMassSystem, set::Settings, R_b_w, Q_b_w)
     return nothing
 end
 
-# function update!(system::PointMassSystem, measure::Measurement)
-#     @unpack points, groups, segments, pulleys, tethers, winches, kite = system
-
-#     kite.pos .= R_b_w * [0.0, 0.0, -min_z]
-#     kite.orient .= Q_b_w
-#     kite.vel .= 0.0
-#     kite.angular_vel .= 0.0
-    
-#     for (winch, tether_length, tether_vel) in zip(winches, measure.tether_length, measure.tether_vel)
-#         winch.tether_length = tether_length
-#         winch.tether_vel = tether_vel
-#     end
-
-#     return nothing
-# end
-
-const MeasureFloat = Float32
-
-@with_kw mutable struct Measurement
-    set_values::MVector{3, MeasureFloat}    = [-50., -1., -1.]
-    tether_length::MVector{3, MeasureFloat} = [51., 51., 51.]
-    tether_vel::MVector{3, MeasureFloat}    = zeros(MeasureFloat, 3)
-    tether_acc::MVector{3, MeasureFloat}    = zeros(MeasureFloat, 3)
-    "elevation and azimuth in spherical coordinate system with columns (left, right) and rows (elevation, azimuth)"
-    sphere_pos::Matrix{MeasureFloat}            = deg2rad.([80.0 80.0; 1.0 -1.0])
-    sphere_vel::Matrix{MeasureFloat}            = zeros(MeasureFloat, 2, 2)
-    sphere_acc::Matrix{MeasureFloat}            = zeros(MeasureFloat, 2, 2)
-    "positive azimuth wind direction in right-handed ENU frame relative to east / x-axis"
-    wind_dir_gnd::MeasureFloat                  = zero(MeasureFloat)
-end
-
-function Base.getproperty(m::Measurement, val::Symbol)
-    if val === :elevation
-        sphere_pos = getfield(m, :sphere_pos)
-        return 0.5(sphere_pos[1, 1] + sphere_pos[1, 2])
-    elseif val === :azimuth
-        sphere_pos = getfield(m, :sphere_pos)
-        return 0.5(sphere_pos[2, 1] + sphere_pos[2, 2])
-    else
-        return getfield(m, val)
-    end
-end
-
-function measure_to_q(measure::Measurement, R_cad_body=I(3))
+function initial_orient(set::Settings, R_cad_body=I(3))
     x = [0, 0, -1] # laying flat along x axis
     z = [1, 0, 0] # laying flat along x axis
-    x = rotate_around_y(x, -measure.elevation)
-    z = rotate_around_y(z, -measure.elevation)
-    x = rotate_around_z(x, measure.azimuth)
-    z = rotate_around_z(z, measure.azimuth)
+    x = rotate_around_y(x, -deg2rad(set.elevation))
+    z = rotate_around_y(z, -deg2rad(set.elevation))
+    x = rotate_around_z(x, deg2rad(set.azimuth))
+    z = rotate_around_z(z, deg2rad(set.azimuth))
     R_b_w = R_cad_body' * hcat(x, z × x, z)
     Q_b_w = rotation_matrix_to_quaternion(R_b_w)
     return Q_b_w, R_b_w
 end
-
diff --git a/src/precompile.jl b/src/precompile.jl
index aaf928b7..5d92a3f4 100644
--- a/src/precompile.jl
+++ b/src/precompile.jl
@@ -87,11 +87,9 @@ end
             set.quasi_static = false
             set.physical_model = "ram"
             s = RamAirKite(set)
-            measure = Measurement()
 
             # Initialize at elevation
-            measure.sphere_pos .= deg2rad.([60.0 60.0; 1.0 -1.0])
-            KiteModels.init_sim!(s, measure; prn=false, precompile=true)
+            KiteModels.init_sim!(s; prn=false, precompile=true)
             @info "Copying $output_file to $prob_file !"
             cp(output_file, prob_file; force=true)
         end
diff --git a/src/ram_air_kite.jl b/src/ram_air_kite.jl
index 1050da4f..9857cc57 100644
--- a/src/ram_air_kite.jl
+++ b/src/ram_air_kite.jl
@@ -45,7 +45,7 @@ $(TYPEDFIELDS)
     guesses::Vector{Pair{Num, Real}} = Pair{Num, Real}[]
 
     set_set_values::Function       = () -> nothing
-    set_measure::Function          = () -> nothing
+    set_wind_dir::Function          = () -> nothing
     set_vsm::Function              = () -> nothing
     set_unknowns::Function         = () -> nothing
     set_nonstiff::Function         = () -> nothing
@@ -176,7 +176,7 @@ function SysState(s::RamAirKite, zoom=1.0)
 end
 
 """
-    init_sim!(s::RamAirKite, measure::Measurement; prn=true, precompile=false) -> Nothing
+    init_sim!(s::RamAirKite; prn=true, precompile=false) -> Nothing
 
 Initialize a kite power system model. 
 
@@ -186,7 +186,7 @@ and only update the state variables. Otherwise, it will create a new model from
 # Fast path (serialized model exists):
 1. Loads existing ODEProblem from disk
 2. Calls `reinit!` to update state variables
-3. Sets up integrator with current measurements
+3. Sets up integrator with initial settings
 
 # Slow path (no serialized model):
 1. Creates symbolic MTK system with all equations
@@ -196,24 +196,23 @@ and only update the state variables. Otherwise, it will create a new model from
 
 # Arguments
 - `s::RamAirKite`: The kite system state object  
-- `measure::Measurement`: Initial state measurements
 - `prn::Bool=true`: Whether to print progress information
 - `precompile::Bool=false`: Whether to build problem for precompilation
 
 # Returns
 `Nothing`
 """
-function init_sim!(s::RamAirKite, measure::Measurement;
+function init_sim!(s::RamAirKite;
     solver=ifelse(s.set.quasi_static, FBDF(nlsolve=OrdinaryDiffEqNonlinearSolve.NLNewton(relax=0.6)), FBDF()), 
     adaptive=true, prn=true, precompile=false, remake=false, reload=false, lin_outputs=Num[]
 )
-    function init(s, measure)
-        init_Q_b_w, R_b_w = measure_to_q(measure, s.wing.R_cad_body)
+    function init(s)
+        init_Q_b_w, R_b_w = initial_orient(s.set, s.wing.R_cad_body)
         init!(s.point_system, s.set, R_b_w, init_Q_b_w)
 
         init_va_b = R_b_w' * [s.set.v_wind, 0., 0.]
         
-        inputs = create_sys!(s, s.point_system, measure; init_va_b)
+        inputs = create_sys!(s, s.point_system; init_va_b)
         prn && @info "Simplifying the system"
         prn ? (@time (sys, _) = structural_simplify(s.full_sys, (inputs, []))) :
             ((sys, _) = structural_simplify(sys, (inputs, [])))
@@ -235,23 +234,18 @@ function init_sim!(s::RamAirKite, measure::Measurement;
     end
     prob_path = joinpath(KiteUtils.get_data_path(), get_prob_name(s.set; precompile))
     if !ispath(prob_path) || remake
-        init(s, measure)
+        init(s)
     end
-    success = reinit!(s, measure; solver, adaptive, precompile, reload, lin_outputs)
+    success = reinit!(s; solver, adaptive, precompile, reload, lin_outputs)
     if !success
         rm(prob_path)
         @info "Rebuilding the system. This can take some minutes..."
-        init(s, measure)
-        reinit!(s, measure; precompile, prn)
+        init(s)
+        reinit!(s; precompile, prn)
     end
     return nothing
 end
 
-
-function init_sim!(::RamAirKite; prn=true)
-    throw(ArgumentError("Use the function init_sim!(s::RamAirKite, measure::Measurement) instead."))
-end
-
 function linearize(s::RamAirKite; set_values=s.get_set_values(s.integrator))
     isnothing(s.lin_prob) && throw(ArgumentError("Run init_sim! with remake=true and lin_outputs=..."))
     s.set_lin_vsm(s.lin_prob, s.get_vsm(s.integrator))
@@ -270,14 +264,14 @@ This function performs the following operations:
    - Loads a serialized ODEProblem from disk
    - Initializes a new ODE integrator 
    - Generates getter/setter functions for the system
-2. Converts measurement data to quaternion orientation
+2. Converts initial settings to quaternion orientation
 3. Initializes the point mass system with new positions
 4. Sets initial values for all state variables
 5. Reinitializes the ODE integrator
 6. Updates the linearized aerodynamic model
 
 This is more efficient than `init!` as it reuses the existing model structure
-and only updates the state variables to match the current `measure`.
+and only updates the state variables to match the current initial settings.
 
 # Arguments
 - `s::RamAirKite`: The kite power system state object
@@ -290,7 +284,7 @@ and only updates the state variables to match the current `measure`.
 - `ArgumentError`: If no serialized problem exists (run `init_sim!` first)
 """
 function reinit!(
-    s::RamAirKite, measure::Measurement; 
+    s::RamAirKite; 
     solver=ifelse(s.set.quasi_static, FBDF(nlsolve=OrdinaryDiffEqNonlinearSolve.NLNewton(relax=0.4, max_iter=1000)), FBDF()),
     adaptive=true,
     prn=true, 
@@ -300,7 +294,7 @@ function reinit!(
 )
     isnothing(s.point_system) && throw(ArgumentError("PointMassSystem not defined"))
 
-    init_Q_b_w, R_b_w = measure_to_q(measure, s.wing.R_cad_body)
+    init_Q_b_w, R_b_w = initial_orient(s.set, s.wing.R_cad_body)
     init!(s.point_system, s.set, R_b_w, init_Q_b_w)
     
     if isnothing(s.prob) || reload
@@ -344,7 +338,7 @@ function generate_getters!(s, sym_vec)
     ])
 
     set_set_values = setp(sys, sys.set_values)
-    set_measure = setp(sys, sys.measured_wind_dir_gnd)
+    set_wind_dir = setp(sys, sys.upwind_dir)
     set_vsm = setp(sys, vsm_sym)
     set_unknowns = setu(sys, sym_vec)
     set_nonstiff = setu(sys, get_nonstiff_unknowns(s))
@@ -385,7 +379,7 @@ function generate_getters!(s, sym_vec)
     get_pos = getu(sys, sys.pos)
 
     s.set_set_values = (integ, val) -> set_set_values(integ, val)
-    s.set_measure = (integ, val) -> set_measure(integ, val)
+    s.set_wind_dir = (integ, val) -> set_wind_dir(integ, val)
     s.set_vsm = (integ, val) -> set_vsm(integ, val)
     s.set_unknowns = (integ, val) -> set_unknowns(integ, val)
     s.set_nonstiff = (integ, val) -> set_nonstiff(integ, val)
@@ -430,12 +424,12 @@ function linearize_vsm!(s::RamAirKite, integ=s.integrator)
     nothing
 end
 
-function next_step!(s::RamAirKite, set_values=nothing; measure::Union{Measurement, Nothing}=nothing, dt=1/s.set.sample_freq, vsm_interval=1)
+function next_step!(s::RamAirKite, set_values=nothing; upwind_dir=nothing, dt=1/s.set.sample_freq, vsm_interval=1)
     if (!isnothing(set_values)) 
         s.set_set_values(s.integrator, set_values)
     end
-    if (!isnothing(measure))
-        s.set_measure(s.integrator, measure.wind_dir_gnd)
+    if (!isnothing(upwind_dir))
+        s.set_wind_dir(s.integrator, upwind_dir)
     end
     if vsm_interval != 0 && s.iter % vsm_interval == 0
         linearize_vsm!(s)
diff --git a/test/bench3.jl b/test/bench3.jl
index 199b99a6..9dce165a 100644
--- a/test/bench3.jl
+++ b/test/bench3.jl
@@ -23,7 +23,7 @@ msg=""
 
 function set_defaults()
     KiteModels.clear!(kps)
-    kps.set.l_tether = 150.0
+    kps.set.l_tethers[1] = 150.0
     kps.set.elevation = 60.0
     kps.set.area = 20.0
     kps.set.rel_side_area = 50.0
@@ -36,7 +36,7 @@ end
 
 function init_392()
     KiteModels.clear!(kps)
-    kps.set.l_tether = 392.0
+    kps.set.l_tethers[1] = 392.0
     kps.set.elevation = 70.0
     kps.set.area = 10.0
     kps.set.rel_side_area = 50.0
diff --git a/test/bench4.jl b/test/bench4.jl
index 217421bc..97d9dd22 100644
--- a/test/bench4.jl
+++ b/test/bench4.jl
@@ -19,7 +19,7 @@ msg = String[]
 
     function set_defaults()
         KiteModels.clear!(kps4)
-        kps4.set.l_tether = 150.0
+        kps4.set.l_tethers[1] = 150.0
         kps4.set.elevation = 60.0
         kps4.set.area = 20.0
         kps4.set.rel_side_area = 50.0
@@ -32,7 +32,7 @@ msg = String[]
 
     function init_392()
         KiteModels.clear!(kps4)
-        kps4.set.l_tether = 392.0
+        kps4.set.l_tethers[1] = 392.0
         kps4.set.elevation = 70.0
         kps4.set.area = 10.0
         kps4.set.rel_side_area = 50.0
@@ -43,7 +43,7 @@ msg = String[]
 
     function init_150()
         KiteModels.clear!(kps4)
-        kps4.set.l_tether = 150.0
+        kps4.set.l_tethers[1] = 150.0
         kps4.set.elevation = 70.0
         kps4.set.area = 10.0
         kps4.set.rel_side_area = 50.0
diff --git a/test/create_xz_file.jl b/test/create_xz_file.jl
index 4be7ad91..c6407c85 100644
--- a/test/create_xz_file.jl
+++ b/test/create_xz_file.jl
@@ -32,17 +32,10 @@ s.set.abs_tol = 1e-2
 s.set.rel_tol = 1e-2
 toc()
 
-# init_Q_b_w, R_b_w = KiteModels.measure_to_q(measure)
-# init_kite_pos = init!(s.point_system, s.set, R_b_w, init_Q_b_w)
-# plot(s.point_system, 0.0; zoom=false, front=true)
-
-measure = Measurement()
-
 # Initialize at elevation
 s.point_system.winches[2].tether_length += 0.2
 s.point_system.winches[3].tether_length += 0.2
-measure.sphere_pos .= deg2rad.([70.0 70.0; 1.0 -1.0])
-KiteModels.init_sim!(s, measure; remake=false, reload=true)
+KiteModels.init_sim!(s; remake=false, reload=true)
 sys = s.sys
 
 @info "System initialized at:"
diff --git a/test/test_for_precompile.jl b/test/test_for_precompile.jl
index ebdfe367..6d119e17 100644
--- a/test/test_for_precompile.jl
+++ b/test/test_for_precompile.jl
@@ -111,8 +111,7 @@ if ! haskey(ENV, "NO_MTK")
     set.segments = 2
     set_values = [-50, -1.1, -1.1]
     mtk_kite = RamAirKite(set)
-    measure = Measurement()
-    KiteModels.init_sim!(mtk_kite, measure)
+    KiteModels.init_sim!(mtk_kite)
     logger = Logger(length(mtk_kite.point_system.points), 5)
 
     for i in 1:5
diff --git a/test/test_kps3.jl b/test/test_kps3.jl
index 873c9a27..0013c4e6 100644
--- a/test/test_kps3.jl
+++ b/test/test_kps3.jl
@@ -18,7 +18,7 @@ end
 
 function set_defaults()
     KiteModels.clear!(kps)
-    kps.set.l_tether = 150.0
+    kps.set.l_tethers[1] = 150.0
     kps.set.elevation = 60.0
     kps.set.area = 20.0
     kps.set.rel_side_area = 50.0
@@ -31,7 +31,7 @@ end
 
 function init_392()
     KiteModels.clear!(kps)
-    kps.set.l_tether = 392.0
+    kps.set.l_tethers[1] = 392.0
     kps.set.elevation = 70.0
     kps.set.area = 10.0
     kps.set.rel_side_area = 50.0
@@ -243,7 +243,7 @@ z= nothing
     res = test_initial_condition(initial_x)
 
     my_state = kps
-    kps.set.l_tether = 392.0
+    kps.set.l_tethers[1] = 392.0
     kps.set.elevation = 70.0
     kps.set.area = 10.0
     kps.set.v_wind = 9.1
diff --git a/test/test_kps3_fails.jl b/test/test_kps3_fails.jl
index edb3ef85..8da92236 100644
--- a/test/test_kps3_fails.jl
+++ b/test/test_kps3_fails.jl
@@ -25,7 +25,7 @@ end
 
 function set_defaults()
     KiteModels.clear!(kps)
-    kps.set.l_tether = 150.0
+    kps.set.l_tethers[1] = 150.0
     kps.set.elevation = 60.0
     kps.set.area = 20.0
     kps.set.rel_side_area = 50.0
diff --git a/test/test_kps4.jl b/test/test_kps4.jl
index 958d14a1..9864ce55 100644
--- a/test/test_kps4.jl
+++ b/test/test_kps4.jl
@@ -21,7 +21,7 @@ poss, vels = nothing, nothing
 
 function set_defaults()
     KiteModels.clear!(kps4)
-    kps4.set.l_tether = 150.0
+    kps4.set.l_tethers[1] = 150.0
     kps4.set.elevation = 60.0
     kps4.set.area = 20.0
     kps4.set.rel_side_area = 50.0
@@ -35,7 +35,7 @@ end
 
 function init_392()
     KiteModels.clear!(kps4)
-    kps4.set.l_tether = 392.0
+    kps4.set.l_tethers[1] = 392.0
     kps4.set.elevation = 70.0
     kps4.set.area = 10.0
     kps4.set.rel_side_area = 50.0
@@ -46,7 +46,7 @@ end
 
 function init_150()
     KiteModels.clear!(kps4)
-    kps4.set.l_tether = 150.0
+    kps4.set.l_tethers[1] = 150.0
     kps4.set.elevation = 70.0
     kps4.set.area = 10.18
     kps4.set.rel_side_area = 30.6
@@ -61,7 +61,7 @@ end
 function init3()
     kps4.set.alpha =  0.08163
     KiteModels.clear!(kps4)
-    kps4.set.l_tether = 150.0 # - kps4.set.height_k - kps4.set.h_bridle
+    kps4.set.l_tethers[1] = 150.0 # - kps4.set.height_k - kps4.set.h_bridle
     kps4.set.area = 10.18
     kps4.set.rel_side_area = 30.6
     kps4.set.mass = 6.21
diff --git a/test/test_ram_air_kite.jl b/test/test_ram_air_kite.jl
index ae1e2d5e..7e766537 100644
--- a/test/test_ram_air_kite.jl
+++ b/test/test_ram_air_kite.jl
@@ -26,12 +26,11 @@ const BUILD_SYS = true
     @info "Creating s:"
     @time s = RamAirKite(set)
 
-    measure = Measurement()
     s.set.abs_tol = 1e-2
     s.set.rel_tol = 1e-2
 
     # Initialize at elevation
-    measure.sphere_pos .= deg2rad.([80.0 80.0; 1.0 -1.0])
+    set.elevation = 80.0
 
     @testset "Model Initialization Chain" begin
         if BUILD_SYS
@@ -45,7 +44,7 @@ const BUILD_SYS = true
 
             # 1. First time initialization - should create new model
             @info "Testing initial init! (should create new model)..."
-            @time KiteModels.init_sim!(s, measure; prn=true)
+            @time KiteModels.init_sim!(s; prn=true)
 
             # Check that serialization worked
             @test isfile(prob_path)
@@ -64,7 +63,7 @@ const BUILD_SYS = true
 
         # 2. First init_sim! - should load from serialized file
         @info "Testing first init_sim! (should load serialized file)..."
-        @time KiteModels.init_sim!(s, measure; prn=true, reload=false)
+        @time KiteModels.init_sim!(s; prn=true, reload=false)
         next_step!(s)
 
         # Check that it's a new integrator
@@ -75,7 +74,7 @@ const BUILD_SYS = true
 
         # 3. Second init_sim! - should reuse existing integrator
         @info "Testing second init_sim! (should reuse integrator)..."
-        @time KiteModels.init_sim!(s, measure; prn=true, reload=false)
+        @time KiteModels.init_sim!(s; prn=true, reload=false)
 
         # This should create a new point system but reuse the existing integrator
         third_integrator_ptr = objectid(s.integrator)
@@ -108,19 +107,19 @@ const BUILD_SYS = true
     end
 
     @testset "State Consistency" begin
-        KiteModels.init_sim!(s, measure, prn=true, reload=false)
+        KiteModels.init_sim!(s, prn=true, reload=false)
         sys_state_before = KiteModels.SysState(s)
 
         # Check quaternion normalization
         @test isapprox(norm(s.integrator[s.sys.Q_b_w]), 1.0, atol=TOL)
 
         # Check elevation matches measurement
-        @test isapprox(sys_state_before.elevation, measure.elevation, atol=1e-2)
+        @test isapprox(sys_state_before.elevation, deg2rad(set.elevation), atol=1e-2)
 
         # Change measurement and reinitialize
-        old_elevation = measure.elevation
-        measure.sphere_pos[1,:] .= deg2rad(85.0)
-        KiteModels.init_sim!(s, measure, prn=true, reload=false)
+        old_elevation = set.elevation
+        set.elevation = 85.0
+        KiteModels.init_sim!(s, prn=true, reload=false)
 
         # Get new state using SysState
         sys_state_after = KiteModels.SysState(s)
@@ -149,7 +148,7 @@ const BUILD_SYS = true
 
     @testset "Simulation Step with SysState" begin
         # Basic step and time advancement test
-        KiteModels.init_sim!(s, measure; prn=true, reload=false)
+        KiteModels.init_sim!(s; prn=true, reload=false)
         sys_state_before = KiteModels.SysState(s)
 
         # Run a simulation step with zero set values
@@ -173,16 +172,14 @@ const BUILD_SYS = true
 
         @testset "Course Direction at 60 Degrees Elevation" begin # Corrected description
             # Initialize at 60 degrees elevation
-            measure.sphere_pos[1,:] .= deg2rad(60.0)
-            @test measure.elevation ≈ deg2rad(60.0) atol=1e-6
-            @test measure.azimuth ≈ 0.0 atol=1e-6
+            set.elevation = 60.0
 
-            KiteModels.init_sim!(s, measure; prn=true)
+            KiteModels.init_sim!(s; prn=true)
 
             # Verify initial conditions using SysState
             sys_state_init = KiteModels.SysState(s)
-            @test sys_state_init.elevation ≈ measure.elevation atol=1e-2 # Use relaxed tolerance consistent with other tests
-            @test sys_state_init.azimuth ≈ measure.azimuth atol=1e-2
+            @test sys_state_init.elevation ≈ deg2rad(set.elevation) atol=1e-2 # Use relaxed tolerance consistent with other tests
+            @test sys_state_init.azimuth ≈ deg2rad(set.azimuth) atol=1e-2
 
             # Run simulation steps
             test_step(s)
@@ -207,18 +204,18 @@ const BUILD_SYS = true
 
         @testset "Steering Response Using SysState" begin
             # Initialize model at moderate elevation
-            measure.sphere_pos .= deg2rad.([70.0 70.0; 1.0 -1.0])
-            KiteModels.init_sim!(s, measure; prn=true, reload=false)
+            set.elevation = 70.0
+            KiteModels.init_sim!(s; prn=true, reload=false)
             test_step(s)
             sys_state_initial = KiteModels.SysState(s)
 
             # steering right
-            KiteModels.init_sim!(s, measure; prn=true, reload=false)
+            KiteModels.init_sim!(s; prn=true, reload=false)
             test_step(s, [0, 10, -10]; steps=20)
             sys_state_right = KiteModels.SysState(s)
 
             # steering left
-            KiteModels.init_sim!(s, measure; prn=true, reload=false)
+            KiteModels.init_sim!(s; prn=true, reload=false)
             test_step(s, [0, -10, 10]; steps=20)
             sys_state_left = KiteModels.SysState(s)