diff --git a/REUSE.toml b/REUSE.toml
index 04307e59..c51e9c18 100644
--- a/REUSE.toml
+++ b/REUSE.toml
@@ -21,7 +21,8 @@ SPDX-FileCopyrightText = "2025 Uwe Fechner, Bart van de Lint"
 SPDX-License-Identifier = "MIT"
 
 [[annotations]]
-path = ["docs/src/*.png", "docs/src/*.svg", "docs/src/*.md", "docs/src/kite_power_tools.drawio"]
+path = ["docs/src/*.png", "docs/src/*.svg", "docs/src/*.md", 
+        "docs/src/symbolic_awe_model/*", "docs/src/kite_power_tools.drawio"]
 SPDX-FileCopyrightText = "2025 Uwe Fechner, Bart van de Lint"
 SPDX-License-Identifier = "CC-BY-4.0"
 
diff --git a/data/model_1.10_ram_dynamic_3_seg.bin.default.xz b/data/model_1.10_ram_dynamic_3_seg.bin.default.xz
index 3d138739..b968c62c 100644
Binary files a/data/model_1.10_ram_dynamic_3_seg.bin.default.xz and b/data/model_1.10_ram_dynamic_3_seg.bin.default.xz differ
diff --git a/data/model_1.11_ram_dynamic_3_seg.bin.default.xz b/data/model_1.11_ram_dynamic_3_seg.bin.default.xz
index f330084c..b7f62751 100644
Binary files a/data/model_1.11_ram_dynamic_3_seg.bin.default.xz and b/data/model_1.11_ram_dynamic_3_seg.bin.default.xz differ
diff --git a/docs/make.jl b/docs/make.jl
index d39485ae..a31a3485 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -28,12 +28,14 @@ makedocs(;
         "Home" => "index.md",
         "Types" => "types.md",
         "Functions" => "functions.md",
-        "SymbolicAWEModel" => "ram_air_kite.md",
+        "SymbolicAWEModel" => [
+		"Overview" => "symbolic_awe_model/ram_air_kite.md",
+		"Examples" => "symbolic_awe_model/examples_ram_air.md",
+		"Tutorial" => "symbolic_awe_model/tutorial_system_structure.md"
+	],
         "Parameters" => "parameters.md",
         "Examples 1p" => "examples.md",
         "Examples 4p" => "examples_4p.md",
-        "Examples SymbolicAWEModel" => "examples_ram_air.md",
-        "SystemStructure for custom models" => "tutorial_system_structure.md",
         "Quickstart" => "quickstart.md",
         "Advanced usage" => "advanced.md",
     ],
diff --git a/docs/src/examples_ram_air.md b/docs/src/symbolic_awe_model/examples_ram_air.md
similarity index 100%
rename from docs/src/examples_ram_air.md
rename to docs/src/symbolic_awe_model/examples_ram_air.md
diff --git a/docs/src/group-slice.svg b/docs/src/symbolic_awe_model/group_slice.svg
similarity index 99%
rename from docs/src/group-slice.svg
rename to docs/src/symbolic_awe_model/group_slice.svg
index 4a5fb540..a25c8e71 100644
--- a/docs/src/group-slice.svg
+++ b/docs/src/symbolic_awe_model/group_slice.svg
@@ -8,7 +8,7 @@
    version="1.1"
    id="svg1"
    inkscape:version="1.4.2 (1:1.4.2+202505120737+ebf0e940d0)"
-   sodipodi:docname="group-slice.svg"
+   sodipodi:docname="group_slice.svg"
    xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape"
    xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd"
    xmlns="http://www.w3.org/2000/svg"
diff --git a/docs/src/oscillating_steering.png b/docs/src/symbolic_awe_model/oscillating_steering.png
similarity index 100%
rename from docs/src/oscillating_steering.png
rename to docs/src/symbolic_awe_model/oscillating_steering.png
diff --git a/docs/src/oscillating_steering_simple.png b/docs/src/symbolic_awe_model/oscillating_steering_simple.png
similarity index 100%
rename from docs/src/oscillating_steering_simple.png
rename to docs/src/symbolic_awe_model/oscillating_steering_simple.png
diff --git a/docs/src/protune-speed-system.jpg b/docs/src/symbolic_awe_model/protune-speed-system.jpg
similarity index 100%
rename from docs/src/protune-speed-system.jpg
rename to docs/src/symbolic_awe_model/protune-speed-system.jpg
diff --git a/docs/src/protune-speed-system.jpg.license b/docs/src/symbolic_awe_model/protune-speed-system.jpg.license
similarity index 100%
rename from docs/src/protune-speed-system.jpg.license
rename to docs/src/symbolic_awe_model/protune-speed-system.jpg.license
diff --git a/docs/src/ram_air_kite.md b/docs/src/symbolic_awe_model/ram_air_kite.md
similarity index 96%
rename from docs/src/ram_air_kite.md
rename to docs/src/symbolic_awe_model/ram_air_kite.md
index de813e6f..8fbb1226 100644
--- a/docs/src/ram_air_kite.md
+++ b/docs/src/symbolic_awe_model/ram_air_kite.md
@@ -35,7 +35,7 @@ Segment
 Pulley(idx, segment_idxs, type)
 Pulley
 Tether
-Winch(idx, model, tether_idxs, tether_length; tether_vel=0.0)
+Winch(idx, model, tether_idxs; tether_length=0.0, tether_vel=0.0)
 Winch
 Wing(idx, group_idxs, R_b_c, pos_cad)
 Wing
@@ -50,4 +50,4 @@ reinit!
 scalar_eqs!
 linear_vsm_eqs!
 force_eqs!
-```
\ No newline at end of file
+```
diff --git a/docs/src/tether_during_sim.png b/docs/src/symbolic_awe_model/tether_during_sim.png
similarity index 100%
rename from docs/src/tether_during_sim.png
rename to docs/src/symbolic_awe_model/tether_during_sim.png
diff --git a/docs/src/tether_sys_struct.png b/docs/src/symbolic_awe_model/tether_sys_struct.png
similarity index 100%
rename from docs/src/tether_sys_struct.png
rename to docs/src/symbolic_awe_model/tether_sys_struct.png
diff --git a/docs/src/symbolic_awe_model/tutorial_system_structure.md b/docs/src/symbolic_awe_model/tutorial_system_structure.md
new file mode 100644
index 00000000..3ee117ae
--- /dev/null
+++ b/docs/src/symbolic_awe_model/tutorial_system_structure.md
@@ -0,0 +1,175 @@
+```@meta
+CurrentModule = KiteModels
+```
+# Custom SystemStructure and SymbolicAWESystem
+
+A custom `SystemStructure` can be used to create models of kite power systems of almost any configuration.
+- custom amount of tethers
+- custom bridle configurations
+- quasi-static or dynamic point masses
+- different amounts of stiffness, damping and diameter on different tether segments
+
+## Precondition
+First, following the [Quickstart](@ref) section up to the installation of the examples. Make sure that
+at least `KiteModels` version 0.8 is installed by typing `using Pkg; Pkg.status()`. To start Julia,
+either use `julia --project`, or `./bin/run_julia`.
+
+## Creating a simple tether
+
+We start by loading the necessary packages and defining settings and parameters.
+
+```julia
+using KiteModels, VortexStepMethod, ControlPlots
+
+set = Settings("system_ram.yaml")
+set.segments = 20
+set.l_tether = 50.0
+dynamics_type = DYNAMIC
+```
+
+Then, we define vectors of the system structure types we are going to use. For this simple example we only need points and segments.
+
+```julia
+points = Point[]
+segments = Segment[]
+
+points = push!(points, Point(1, zeros(3), STATIC; wing_idx=0))
+```
+
+The first point we add is a static point. There are four different [`DynamicsType`](@ref)s to choose from: `STATIC`, `QUASI_STATIC`, `DYNAMIC` and `WING`. `STATIC` just means that the point doesn't move. `DYNAMIC` is a point modeled with acceleration, while `QUASI_STATIC` constrains this acceleration to be zero at all times. A `WING` point is connected to a wing body.
+
+Now we can add `DYNAMIC` points and connect them to each other with segments. `BRIDLE` segments don't need to have a tether, because they have a constant unstretched length.
+```julia
+segment_idxs = Int[]
+for i in 1:set.segments
+    global points, segments
+    point_idx = i+1
+    pos = [0.0, 0.0, i * set.l_tether / set.segments]
+    if i < set.segments
+        push!(points, Point(point_idx, pos, dynamics_type; wing_idx=0))
+    else
+        push!(points, Point(point_idx, pos, dynamics_type; mass=1.0, wing_idx=0))
+    end
+    segment_idx = i
+    push!(segments, Segment(segment_idx, (point_idx-1, point_idx), BRIDLE))
+    push!(segment_idxs, segment_idx)
+end
+```
+
+In order to describe the initial orientation of the structure, we define a [`Transform(idx, elevation, azimuth, heading)`](@ref) with an elevation (-80 degrees), azimuth and heading, and a base position `[0.0, 0.0, 50.0]`.
+```julia
+transforms = [Transform(1, deg2rad(-80), 0.0, 0.0; 
+              base_pos = [0.0, 0.0, 50.0], base_point_idx=points[1].idx,
+              rot_point_idx=points[end].idx)]
+```
+
+From the points, segments and transform we create a [`SystemStructure(name, set)`](@ref), which can be plotted in 2d to quickly investigate if the model is correct.
+```julia
+sys_struct = SystemStructure("tether", set; points, segments, transforms)
+plot(sys_struct, 0.0)
+```
+![SystemStructure visualization](tether_sys_struct.png)
+
+If the system looks good, we can easily model it, by first creating a [`SymbolicAWEModel`](@ref), initializing it and stepping through time.
+```julia
+sam = SymbolicAWEModel(set, sys_struct)
+
+init_sim!(sam; remake=false)
+for i in 1:80
+    plot(sam, i/set.sample_freq)
+    next_step!(sam)
+end
+```
+![Tether during simulation](tether_during_sim.png)
+
+# Using a winch and a tether
+
+Let's try to adjust the length of the tether in the last example. To do this we first need to create a set of segments with a common changing `l0`, called a [`Tether`](@ref).
+```julia
+set.v_wind = 0.0
+tethers = [Tether(1,[segment.idx for segment in segments])]
+```
+As you can see, we just add all of the segments from the simple tether to our [`Tether`](@ref) struct.
+The next step is to create a winch. Each winch can be connected to one or more tethers, so it is possible to connect multiple tethers to the same winch. We have to specify which kind of winch we want to use. For now, only the `TorqueControlledMachine` from `WinchModels.jl` is supported.
+
+```julia
+using WinchModels
+wm = TorqueControlledMachine(set)
+winches = [Winch(1, wm, [1])]
+```
+
+The 2d plot of the [`SystemStructure`](@ref) should still look the same, so we don't have to plot that. We can just create the system, and simulate it. We just need to be sure that we call plot with `t=0.0` to reset the plot.
+
+```julia
+sys_struct = SystemStructure("winch", set; points, segments, tethers, winches, transforms)
+@show set.v_wind
+sam = SymbolicAWEModel(set, sys_struct)
+init_sim!(sam; remake=false)
+ss = SysState(sam)
+
+for i in 1:80
+    plot(sam, (i-1)/set.sample_freq)
+    next_step!(sam; set_values=[-20.0])
+    update_sys_state!(ss, sam)
+end
+@show ss.l_tether[1]
+```
+
+# Using a pulley
+
+First, we need to update some settings. `l_tether` is specified such that the plot window is zoomed in correctly.
+
+```julia
+using KiteModels, VortexStepMethod, ControlPlots
+
+set = se("system_ram.yaml")
+set.v_wind = 10.0
+set.l_tether = 5.0
+set.abs_tol = 1e-4
+set.rel_tol = 1e-4
+dynamics_type = DYNAMIC
+```
+
+Now we create points and segments in a similar manner as in the last example. The `mass` keyword can be used to specify the mass of the point itself. When `mass=0.0` the mass of the point just consists of the tether/segment mass.
+
+```julia
+points = Point[]
+segments = Segment[]
+pulleys = Pulley[]
+
+push!(points, Point(1, [0.0, 0.0, 2.0], STATIC))
+push!(points, Point(2, [2.0, 0.0, 2.0], STATIC))
+push!(points, Point(3, [0.1, 0.0, 1.0], DYNAMIC))
+push!(points, Point(4, [0.1, 0.0, 0.0], DYNAMIC; mass=0.1))
+
+push!(segments, Segment(1, (3,1), BRIDLE))
+push!(segments, Segment(2, (3,2), BRIDLE))
+push!(segments, Segment(3, (3,4), BRIDLE))
+```
+
+Pulleys can be modeled when three or more [`Segment`](@ref)s are connected to a common [`Point`](@ref). When creating a pulley, only two segments are specified: these are the segments of the tether moving through the pulley.
+
+```julia
+push!(pulleys, Pulley(1, (1,2), DYNAMIC))
+```
+
+We can then use a [`Transform`](@ref) to describe the orientation of the initial system. 
+
+```julia
+transforms = [Transform(1, -deg2rad(0.0), 0.0, 0.0; base_pos=[1.0, 0.0, 4.0], base_point_idx=1, rot_point_idx=2)]
+sys_struct = KiteModels.SystemStructure("pulley", set; points, segments, pulleys, transforms)
+plot(sys_struct, 0.0; zoom=false, l_tether=set.l_tether)
+```
+
+If the plot of the [`SystemStructure`](@ref) looks good, we can continue by creating a [`SymbolicAWEModel`](@ref) and simulating through time.
+
+```julia
+sam = SymbolicAWEModel(set, sys_struct)
+init_sim!(sam; remake=false)
+for i in 1:100
+    plot(sam, i/set.sample_freq; zoom=false)
+    next_step!(sam)
+end
+```
+
+
diff --git a/docs/src/tutorial_system_structure.md b/docs/src/tutorial_system_structure.md
deleted file mode 100644
index 396b661d..00000000
--- a/docs/src/tutorial_system_structure.md
+++ /dev/null
@@ -1,78 +0,0 @@
-```@meta
-CurrentModule = KiteModels
-```
-# Custom SystemStructure and SymbolicAWESystem
-
-A custom `SystemStructure` can be used to create models of kite power systems of almost any configuration.
-- custom amount of tethers
-- custom bridle configurations
-- quasi-static or dynamic point masses
-- different amounts of stiffness, damping and diameter on different tether segments
-
-## Precondition
-First, following the [Quickstart](@ref) section up to the installation of the examples. Make sure that
-at least `KiteModels` version 0.8 is installed by typing `using Pkg; Pkg.status()`. To start Julia,
-either use `julia --project`, or `./bin/run_julia`.
-
-## Creating a simple tether
-
-We start by loading the necessary packages and defining settings and parameters.
-
-```julia
-using KiteModels, VortexStepMethod, ControlPlots
-
-set = se("system_ram.yaml")
-set.segments = 20
-dynamics_type = DYNAMIC
-```
-
-Then, we define vectors of the system structure types we are going to use. For this simple example we only need points and segments.
-
-```julia
-points = Point[]
-segments = Segment[]
-
-points = push!(points, Point(1, zeros(3), STATIC; wing_idx=0))
-```
-
-The first point we add is a static point. There are four different [`DynamicsType`](@ref)s to choose from: `STATIC`, `QUASI_STATIC`, `DYNAMIC` and `WING`. `STATIC` just means that the point doesn't move. `DYNAMIC` is a point modeled with acceleration, while `QUASI_STATIC` constrains this acceleration to be zero at all times. A `WING` point is connected to a wing body.
-
-Now we can add `DYNAMIC` points and connect them to each other with segments. `BRIDLE` segments don't need to have a tether, because they have a constant unstretched length.
-```julia
-segment_idxs = Int[]
-for i in 1:set.segments
-    global points, segments
-    point_idx = i+1
-    pos = [0.0, 0.0, i * set.l_tether / set.segments]
-    push!(points, Point(point_idx, pos, dynamics_type; wing_idx=0))
-    segment_idx = i
-    push!(segments, Segment(segment_idx, (point_idx-1, point_idx), BRIDLE))
-    push!(segment_idxs, segment_idx)
-end
-```
-
-In order to describe the initial orientation of the structure, we define a [`Transform(idx, elevation, azimuth, heading)`](@ref) with an elevation (-80 degrees), azimuth and heading, and a base position `[0.0, 0.0, 50.0]`.
-```julia
-transforms = [Transform(1, deg2rad(-80), 0.0, 0.0; 
-              base_pos = [0.0, 0.0, 50.0], base_point_idx=points[1].idx,
-              rot_point_idx=points[end].idx)]
-```
-
-From the points, segments and transform we create a [`SystemStructure(name, set)`](@ref), which can be plotted in 2d to quickly investigate if the model is correct.
-```julia
-sys_struct = SystemStructure("tether", set; points, segments, transforms)
-plot(sys_struct, 0.0)
-```
-![SystemStructure visualization](tether_sys_struct.png)
-
-If the system looks good, we can easily model it, by first creating a [`SymbolicAWEModel`](@ref), initializing it and stepping through time.
-```julia
-sam = SymbolicAWEModel(set, sys_struct)
-
-init_sim!(sam; remake=false)
-for i in 1:80
-    plot(sam, i/set.sample_freq)
-    next_step!(sam)
-end
-```
-![Tether during simulation](tether_during_sim.png)
diff --git a/examples/pulley.jl b/examples/pulley.jl
new file mode 100644
index 00000000..0229738d
--- /dev/null
+++ b/examples/pulley.jl
@@ -0,0 +1,39 @@
+# SPDX-FileCopyrightText: 2025 Bart van de Lint
+#
+# SPDX-License-Identifier: MPL-2.0
+
+using KiteModels, VortexStepMethod, ControlPlots
+
+set = se("system_ram.yaml")
+set.v_wind = 10.0
+set.l_tether = 5.0
+set.abs_tol = 1e-4
+set.rel_tol = 1e-4
+dynamics_type = DYNAMIC
+
+points = Point[]
+segments = Segment[]
+pulleys = Pulley[]
+
+push!(points, Point(1, [0.0, 0.0, 2.0], STATIC))
+push!(points, Point(2, [2.0, 0.0, 2.0], STATIC))
+push!(points, Point(3, [0.1, 0.0, 1.0], DYNAMIC))
+push!(points, Point(4, [0.1, 0.0, 0.0], DYNAMIC; mass=0.1))
+
+push!(segments, Segment(1, (3,1), BRIDLE))
+push!(segments, Segment(2, (3,2), BRIDLE))
+push!(segments, Segment(3, (3,4), BRIDLE))
+
+push!(pulleys, Pulley(1, (1,2), DYNAMIC))
+
+transforms = [Transform(1, -deg2rad(0.0), 0.0, 0.0; base_pos=[1.0, 0.0, 4.0], base_point_idx=1, rot_point_idx=2)]
+sys_struct = KiteModels.SystemStructure("pulley", set; points, segments, pulleys, transforms)
+plot(sys_struct, 0.0; zoom=false, l_tether=set.l_tether)
+
+sam = SymbolicAWEModel(set, sys_struct)
+
+init_sim!(sam; remake=false)
+for i in 1:100
+    plot(sam, i/set.sample_freq; zoom=false)
+    next_step!(sam)
+end
diff --git a/examples/ram_air_kite.jl b/examples/ram_air_kite.jl
index 5ec2bd98..24998b80 100644
--- a/examples/ram_air_kite.jl
+++ b/examples/ram_air_kite.jl
@@ -5,7 +5,7 @@ using Timers
 tic()
 @info "Loading packages "
 
-PLOT = true
+PLOT = false
 using Pkg
 if ! ("LaTeXStrings" ∈ keys(Pkg.project().dependencies))
     using TestEnv; TestEnv.activate()
@@ -30,7 +30,7 @@ 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 = Settings("system_ram.yaml")
 set.segments = 3
 set_values = [-50, 0.0, 0.0]  # Set values of the torques of the three winches. [Nm]
 set.quasi_static = false
@@ -42,13 +42,9 @@ sam.set.abs_tol = 1e-2
 sam.set.rel_tol = 1e-2
 toc()
 
-# init_Q_b_w, R_b_w = KiteModelsam.initial_orient(sam.set)
-# init_kite_pos = init!(sam.sys_struct, sam.set, R_b_w, init_Q_b_w)
-# plot(sam.sys_struct, 0.0; zoom=false, front=false)
-
 # Initialize at elevation
-set.l_tethers[2] += 0.2
-set.l_tethers[3] += 0.2
+set.l_tethers[2] += 0.4
+set.l_tethers[3] += 0.4
 init_sim!(sam; remake=false, reload=false)
 sys = sam.sys
 
@@ -63,7 +59,7 @@ sys_state = SysState(sam)
 t = 0.0
 runtime = 0.0
 integ_runtime = 0.0
-bias = set.quasi_static ? 0.45 : 0.40
+bias = set.quasi_static ? 0.45 : 0.35
 t0 = sam.integrator.t
 
 try
@@ -149,3 +145,5 @@ display(p)
 @info "Performance:" times_realtime=(total_time/2)/runtime integrator_times_realtime=(total_time/2)/integ_runtime
 
 # 55x realtime (PLOT=false, CPU: Intel i9-9980HK (16) @ 5.000GHz)
+# 40-65x realtime (PLOT=false, CPU: Intel i9-9980HK (16) @ 5.000GHz) - commit 6620ed5d0a38e96930615aad9a66e4cd666955f2
+# 40x realtime (PLOT=false, CPU: Intel i9-9980HK (16) @ 5.000GHz) - commit 88a78894038d3cbd50fbff83dfbe5c26266b0637
diff --git a/examples/sam_tutorial.jl b/examples/sam_tutorial.jl
new file mode 100644
index 00000000..620c9f9b
--- /dev/null
+++ b/examples/sam_tutorial.jl
@@ -0,0 +1,86 @@
+# SPDX-FileCopyrightText: 2025 Bart van de Lint
+#
+# SPDX-License-Identifier: MPL-2.0
+using KiteModels, VortexStepMethod, ControlPlots
+
+set = Settings("system_ram.yaml")
+set.segments = 20
+set.l_tether = 50.0
+dynamics_type = DYNAMIC
+
+points = Point[]
+segments = Segment[]
+
+points = push!(points, Point(1, zeros(3), STATIC; wing_idx=0))
+
+segment_idxs = Int[]
+for i in 1:set.segments
+    global points, segments
+    point_idx = i+1
+    pos = [0.0, 0.0, -i * set.l_tether / set.segments]
+    if i < set.segments
+        push!(points, Point(point_idx, pos, dynamics_type; wing_idx=0))
+    else
+        push!(points, Point(point_idx, pos, dynamics_type; mass=1.0, wing_idx=0))
+    end
+    segment_idx = i
+    push!(segments, Segment(segment_idx, (point_idx-1, point_idx), BRIDLE))
+    push!(segment_idxs, segment_idx)
+end
+
+transforms = [Transform(1, deg2rad(-80), 0.0, 0.0; 
+              base_pos = [0.0, 0.0, 50.0], base_point_idx=points[1].idx,
+              rot_point_idx=points[end].idx)]
+
+sys_struct = SystemStructure("tether", set; points, segments, transforms)
+plot(sys_struct, 0.0)
+
+sam = SymbolicAWEModel(set, sys_struct)
+
+init_sim!(sam; remake=false)
+for i in 1:80
+    plot(sam, i/set.sample_freq)
+    next_step!(sam)
+end
+
+# ADDING A WINCH
+set.v_wind = 0.0
+tethers = [Tether(1,[segment.idx for segment in segments])]
+
+using WinchModels
+wm = TorqueControlledMachine(set)
+winches = [Winch(1, wm, [1])]
+
+sys_struct = SystemStructure("winch", set; points, segments, tethers, winches, transforms)
+@show set.v_wind
+sam = SymbolicAWEModel(set, sys_struct)
+init_sim!(sam; remake=false)
+ss = SysState(sam)
+
+for i in 1:80
+    plot(sam, (i-1)/set.sample_freq)
+    next_step!(sam; set_values=[-20.0])
+    update_sys_state!(ss, sam)
+end
+@show ss.l_tether[1]
+
+# ADDING A PULLEY
+last_pos = points[end].pos_cad
+@show last_pos
+push!(points, Point(22, last_pos + [0, 0, -5], DYNAMIC))
+push!(points, Point(23, last_pos + [0, 0, -5], STATIC))
+push!(segments, Segment(21, (21,22), BRIDLE))
+push!(segments, Segment(22, (21,23), BRIDLE))
+pulleys = [Pulley(1, (21,22), DYNAMIC)]
+transforms[1].elevation = deg2rad(-85.0)
+sys_struct = SystemStructure("pulley", set; points, segments, tethers, winches, pulleys, transforms)
+plot(sys_struct, 0.0)
+
+sam = SymbolicAWEModel(set, sys_struct)
+
+init_sim!(sam; remake=false)
+for i in 1:80
+    plot(sam, i/set.sample_freq)
+    next_step!(sam; set_values=[-10.0])
+end
+
diff --git a/examples/sparse.jl b/examples/sparse.jl
new file mode 100644
index 00000000..bfa6cfc2
--- /dev/null
+++ b/examples/sparse.jl
@@ -0,0 +1,41 @@
+# SPDX-FileCopyrightText: 2025 Bart van de Lint
+#
+# SPDX-License-Identifier: MPL-2.0
+
+using DifferentiationInterface
+using SparseArrays
+using ADTypes
+
+# 1. Define your function
+function f(x)
+    y = similar(x, length(x)-1)
+    for i in 1:length(y)
+        y[i] = x[i+1]^2 - x[i]^2
+    end
+    return y
+end
+
+# 2. Choose a finite difference backend
+backend = AutoFiniteDiff()
+
+# 3. Provide the known sparsity pattern as a SparseMatrixCSC{Bool}
+#    For f: ℝⁿ → ℝⁿ⁻¹, each output depends only on x[i] and x[i+1]
+n = 5
+rows = Int[]
+cols = Int[]
+for i in 1:n-1
+    push!(rows, i); push!(cols, i)     # y[i] depends on x[i]
+    push!(rows, i); push!(cols, i+1)   # y[i] depends on x[i+1]
+end
+S = sparse(rows, cols, trues(length(rows)), n-1, n) # Bool-valued sparsity pattern
+
+# 4. Prepare the Jacobian using the known sparsity pattern
+jac_prep = prepare_jacobian(f, backend, zeros(n); sparsity_pattern=S)
+
+# 5. Compute the sparse Jacobian at a point
+x = randn(n)
+J = jacobian(f, jac_prep, backend, x)
+
+# 6. Show the result
+println("Jacobian (sparse):")
+display(J)
diff --git a/examples/tether.jl b/examples/tether.jl
deleted file mode 100644
index e1983659..00000000
--- a/examples/tether.jl
+++ /dev/null
@@ -1,37 +0,0 @@
-# SPDX-FileCopyrightText: 2025 Bart van de Lint
-#
-# SPDX-License-Identifier: MPL-2.0
-
-using KiteModels, VortexStepMethod, ControlPlots
-
-set = se("system_ram.yaml")
-set.segments = 20
-dynamics_type = DYNAMIC
-
-points = Point[]
-segments = Segment[]
-
-points = push!(points, Point(1, zeros(3), STATIC; wing_idx=0))
-
-segment_idxs = Int[]
-for i in 1:set.segments
-    global points, segments
-    point_idx = i+1
-    pos = [0.0, 0.0, i * set.l_tether / set.segments]
-    push!(points, Point(point_idx, pos, dynamics_type; wing_idx=0))
-    segment_idx = i
-    push!(segments, Segment(segment_idx, (point_idx-1, point_idx), BRIDLE))
-    push!(segment_idxs, segment_idx)
-end
-
-transforms = [Transform(1, deg2rad(-80), 0.0, 0.0, [0.0, 0.0, 50.0], points[1].idx; rot_point_idx=points[end].idx)]
-sys_struct = SystemStructure("tether", set; points, segments, transforms)
-plot(sys_struct, 0.0)
-
-sam = SymbolicAWEModel(set, sys_struct)
-
-init_sim!(sam; remake=false)
-for i in 1:80
-    plot(sam, i/set.sample_freq)
-    next_step!(sam)
-end
diff --git a/ext/KiteModelsControlPlotsExt.jl b/ext/KiteModelsControlPlotsExt.jl
index aa91c2e0..281c8020 100644
--- a/ext/KiteModelsControlPlotsExt.jl
+++ b/ext/KiteModelsControlPlotsExt.jl
@@ -18,11 +18,11 @@ function ControlPlots.plot(sys::SystemStructure, reltime; l_tether=50.0, wing_po
         xlim = (pos[end][2] - 6, pos[end][2]+6)
         ylim = (pos[end][3] - 10, pos[end][3]+2)
     elseif !zoom && !front
-        xlim = (-5, l_tether+10)
-        ylim = (-5, l_tether+10)
+        xlim = (-0.1l_tether, 1.2l_tether)
+        ylim = (-0.1l_tether, 1.2l_tether)
     elseif !zoom && front
-        xlim = (-12.5 - 0.5l_tether, 12.5 + 0.5l_tether)
-        ylim = (-5, l_tether+10)
+        xlim = (-0.75l_tether, 0.75l_tether)
+        ylim = (-0.1l_tether, 1.2l_tether)
     end
     ControlPlots.plot2d(pos, seg, reltime; zoom, front, xlim, ylim, dz_zoom=0.6, xy)
 end
diff --git a/src/mtk_model.jl b/src/mtk_model.jl
index 6b34c5a9..61b6e8bc 100644
--- a/src/mtk_model.jl
+++ b/src/mtk_model.jl
@@ -65,13 +65,72 @@ function rotation_matrix_to_quaternion(R)
     return [w, x, y, z]
 end
 
-function get_pos_w(sys_struct::SystemStructure, idx::Int16)
-    return sys_struct.points[idx].pos_w
+get_pos_w(sys_struct::SystemStructure, idx::Int16) = sys_struct.points[idx].pos_w
+@register_array_symbolic get_pos_w(sys::SystemStructure, idx::Int16) begin
+    size=(3,)
+    eltype=SimFloat
+end
+get_pos_b(sys_struct::SystemStructure, idx::Int16) = sys_struct.points[idx].pos_b
+@register_array_symbolic get_pos_b(sys::SystemStructure, idx::Int16) begin
+    size=(3,)
+    eltype=SimFloat
+end
+get_wing_pos_w(sys_struct::SystemStructure, idx::Int16) = sys_struct.wings[idx].pos_w
+@register_array_symbolic get_wing_pos_w(sys::SystemStructure, idx::Int16) begin
+    size=(3,)
+    eltype=SimFloat
 end
-@register_array_symbolic get_pos_w(sys::SystemStructure, point_idx::Int16) begin
-    size=size(KVec3)
+get_wing_vel_w(sys_struct::SystemStructure, idx::Int16) = sys_struct.wings[idx].vel_w
+@register_array_symbolic get_wing_vel_w(sys::SystemStructure, idx::Int16) begin
+    size=(3,)
     eltype=SimFloat
 end
+get_orient(sys_struct::SystemStructure, idx::Int16) = sys_struct.wings[idx].orient
+@register_array_symbolic get_orient(sys::SystemStructure, idx::Int16) begin
+    size=(4,)
+    eltype=SimFloat
+end
+get_angular_vel(sys_struct::SystemStructure, idx::Int16) = sys_struct.wings[idx].angular_vel
+@register_array_symbolic get_angular_vel(sys::SystemStructure, idx::Int16) begin
+    size=(3,)
+    eltype=SimFloat
+end
+get_mass(sys_struct::SystemStructure, idx::Int16) = sys_struct.points[idx].mass
+@register_symbolic get_mass(sys::SystemStructure, idx::Int16)
+get_l0(sys_struct::SystemStructure, idx::Int16) = sys_struct.segments[idx].l0
+@register_symbolic get_l0(sys::SystemStructure, idx::Int16)
+get_diameter(sys_struct::SystemStructure, idx::Int16) = sys_struct.segments[idx].diameter
+@register_symbolic get_diameter(sys::SystemStructure, idx::Int16)
+get_compression_frac(sys_struct::SystemStructure, idx::Int16) = sys_struct.segments[idx].compression_frac
+@register_symbolic get_compression_frac(sys::SystemStructure, idx::Int16)
+get_moment_frac(sys_struct::SystemStructure, idx::Int16) = sys_struct.groups[idx].moment_frac
+@register_symbolic get_moment_frac(sys::SystemStructure, idx::Int16)
+get_sum_length(sys_struct::SystemStructure, idx::Int16) = sys_struct.pulleys[idx].sum_length
+@register_symbolic get_sum_length(sys::SystemStructure, idx::Int16)
+get_tether_length(sys_struct::SystemStructure, idx::Int16) = sys_struct.winches[idx].tether_length
+@register_symbolic get_tether_length(sys::SystemStructure, idx::Int16)
+get_tether_vel(sys_struct::SystemStructure, idx::Int16) = sys_struct.winches[idx].tether_vel
+@register_symbolic get_tether_vel(sys::SystemStructure, idx::Int16)
+
+get_set_mass(set) = set.mass
+@register_symbolic get_set_mass(set::Settings)
+get_rho_tether(set) = set.rho_tether
+@register_symbolic get_rho_tether(set::Settings)
+get_e_tether(set) = set.e_tether
+@register_symbolic get_e_tether(set::Settings)
+get_damping(set) = set.damping
+@register_symbolic get_damping(set::Settings)
+get_c_spring(set) = set.c_spring
+@register_symbolic get_c_spring(set::Settings)
+get_cd_tether(set) = set.cd_tether
+@register_symbolic get_cd_tether(set::Settings)
+
+
+get_v_wind(set::Settings) = set.v_wind
+@register_symbolic get_v_wind(set::Settings)
+get_upwind_dir(set::Settings) = set.upwind_dir
+@register_symbolic get_upwind_dir(set::Settings)
+
 
 """
     force_eqs!(s, system, eqs, defaults, guesses; kwargs...)
@@ -101,7 +160,7 @@ Tuple containing:
 - Tether forces on wing
 - Tether moments on wing
 """
-function force_eqs!(s, system, eqs, defaults, guesses; 
+function force_eqs!(s, system, psys, pset, eqs, defaults, guesses; 
         R_b_w, wing_pos, wing_vel, wind_vec_gnd, group_aero_moment, twist_angle, twist_ω, stabilize, set_values, fix_nonstiff)
 
     @parameters acc_multiplier = 1
@@ -111,7 +170,6 @@ function force_eqs!(s, system, eqs, defaults, guesses;
     # ==================== POINTS ==================== #
     tether_wing_force = zeros(Num, length(wings), 3)
     tether_wing_moment = zeros(Num, length(wings), 3)
-    @parameters psys::SystemStructure = system
     @variables begin
         pos(t)[1:3, eachindex(points)]
         vel(t)[1:3, eachindex(points)]
@@ -119,6 +177,9 @@ function force_eqs!(s, system, eqs, defaults, guesses;
         point_force(t)[1:3, eachindex(points)]
         tether_r(t)[1:3, eachindex(points)]
         point_mass(t)[eachindex(points)]
+        chord_b(t)[1:3, eachindex(points)]
+        normal(t)[1:3, eachindex(points)]
+        pos_b(t)[1:3, eachindex(points)]
 
         spring_force_vec(t)[1:3, eachindex(segments)]
         drag_force(t)[1:3, eachindex(segments)]
@@ -126,11 +187,11 @@ function force_eqs!(s, system, eqs, defaults, guesses;
     end
     for point in points
         F::Vector{Num} = zeros(Num, 3)
-        mass = 0.0
+        mass = get_mass(psys, point.idx)
         in_bridle = false
         for segment in segments
             if point.idx in segment.point_idxs
-                mass_per_meter = s.set.rho_tether * π * (segment.diameter/2)^2
+                mass_per_meter = get_rho_tether(pset) * π * (get_diameter(psys, segment.idx)/2)^2
                 inverted = segment.point_idxs[2] == point.idx
                 if inverted
                     F .-= spring_force_vec[:, segment.idx]
@@ -178,14 +239,17 @@ function force_eqs!(s, system, eqs, defaults, guesses;
                     and should be part of exactly 1 wing.")
 
                 fixed_pos = group.le_pos
-                chord_b = point.pos_b - fixed_pos
-                normal = chord_b × group.y_airf
-                pos_b = fixed_pos + cos(twist_angle[group.idx]) * chord_b - sin(twist_angle[group.idx]) * normal
+                eqs = [
+                    eqs
+                    chord_b[:, point.idx]   ~ get_pos_b(psys, point.idx) .- fixed_pos
+                    normal[:, point.idx]    ~ chord_b[:, point.idx] × group.y_airf
+                    pos_b[:, point.idx]     ~ fixed_pos .+ cos(twist_angle[group.idx]) * chord_b[:, point.idx] - sin(twist_angle[group.idx]) * normal[:, point.idx]
+                ]
             elseif found == 0
-                pos_b = point.pos_b
                 eqs = [
                     eqs
-                    tether_r[:, point.idx] ~ pos[:, point.idx] - wing_pos[point.wing_idx, :]
+                    pos_b[:, point.idx]     ~ get_pos_b(psys, point.idx)
+                    tether_r[:, point.idx]  ~ pos[:, point.idx] - wing_pos[point.wing_idx, :]
                 ]
                 tether_wing_moment[point.wing_idx, :] .+= tether_r[:, point.idx] × point_force[:, point.idx]
             end
@@ -193,7 +257,7 @@ function force_eqs!(s, system, eqs, defaults, guesses;
             
             eqs = [
                 eqs
-                pos[:, point.idx]    ~ wing_pos[point.wing_idx, :] + R_b_w[point.wing_idx, :, :] * pos_b
+                pos[:, point.idx]    ~ wing_pos[point.wing_idx, :] + R_b_w[point.wing_idx, :, :] * pos_b[:, point.idx]
                 vel[:, point.idx]    ~ zeros(3)
                 acc[:, point.idx]    ~ zeros(3)
             ]
@@ -257,6 +321,8 @@ function force_eqs!(s, system, eqs, defaults, guesses;
             group_tether_moment(t)[eachindex(groups)]
             tether_force(t)[eachindex(groups), eachindex(groups[1].point_idxs)]
             tether_moment(t)[eachindex(groups), eachindex(groups[1].point_idxs)]
+            r_group(t)[eachindex(groups), eachindex(groups[1].point_idxs)]
+            r_vec(t)[eachindex(groups), eachindex(groups[1].point_idxs), 1:3]
         end
     end
     
@@ -279,16 +345,16 @@ function force_eqs!(s, system, eqs, defaults, guesses;
         init_z_airf = x_airf × group.y_airf
         z_airf = x_airf * sin(twist_angle[group.idx]) + init_z_airf * cos(twist_angle[group.idx])
         for (i, point_idx) in enumerate(group.point_idxs)
-            r = (points[point_idx].pos_b - (group.le_pos + group.moment_frac*group.chord)) ⋅ normalize(group.chord)
             eqs = [
                 eqs
-                tether_force[group.idx, i] ~ (point_force[:, point_idx] ⋅ (R_b_w[wing.idx, :, :] * -z_airf))
-                tether_moment[group.idx, i] ~ r * tether_force[group.idx, i]
+                r_vec[group.idx, i, :]      ~ (get_pos_b(psys, point_idx) .- (group.le_pos + get_moment_frac(psys, group.idx)*group.chord))
+                r_group[group.idx, i]       ~ r_vec[group.idx, i, :] ⋅ normalize(group.chord)
+                tether_force[group.idx, i]  ~ (point_force[:, point_idx] ⋅ (R_b_w[wing.idx, :, :] * -z_airf))
+                tether_moment[group.idx, i] ~ r_group[group.idx, i] * tether_force[group.idx, i]
             ]
         end
         
-        inertia = 1/3 * (s.set.mass/length(groups)) * (norm(group.chord))^2 # plate inertia around leading edge
-        @assert !(inertia ≈ 0.0)
+        inertia = 1/3 * (get_set_mass(pset)/length(groups)) * (norm(group.chord))^2 # plate inertia around leading edge
         @parameters twist_damp = 50
         @parameters max_twist = deg2rad(90)
 
@@ -354,37 +420,32 @@ function force_eqs!(s, system, eqs, defaults, guesses;
         # if s.set.quasi_static
             guesses = [
                 guesses
-                [segment_vec[i, segment.idx] => points[p2].pos_w[i] - points[p1].pos_w[i] for i in 1:3]
+                [segment_vec[i, segment.idx] => get_pos_w(psys, p2)[i] - get_pos_w(psys, p1)[i] for i in 1:3]
             ]
         # end
 
-        if segment.type == BRIDLE
-            in_pulley = 0
-            for pulley in pulleys
-                if segment.idx == pulley.segment_idxs[1] # each bridle segment has to be part of no pulley or one pulley
-                    eqs = [
-                        eqs
-                        l0[segment.idx] ~ pulley_l0[pulley.idx]
-                    ]
-                    in_pulley += 1
-                end
-                if segment.idx == pulley.segment_idxs[2]
-                    eqs = [
-                        eqs
-                        l0[segment.idx] ~ pulley.sum_length - pulley_l0[pulley.idx]
-                    ]
-                    in_pulley += 1
-                end
+        in_pulley = 0
+        for pulley in pulleys
+            if segment.idx == pulley.segment_idxs[1] # each bridle segment has to be part of no pulley or one pulley
+                eqs = [
+                    eqs
+                    l0[segment.idx] ~ pulley_l0[pulley.idx]
+                ]
+                in_pulley += 1
             end
-            if in_pulley == 0
+            if segment.idx == pulley.segment_idxs[2]
                 eqs = [
                     eqs
-                    l0[segment.idx] ~ segment.l0
+                    l0[segment.idx] ~ get_sum_length(psys, pulley.idx) - pulley_l0[pulley.idx]
                 ]
+                in_pulley += 1
             end
-            (in_pulley > 1) && error("Bridle segment number $(segment.idx) is part of
-                $in_pulley pulleys, and should be part of either 0 or 1 pulleys.")
-        elseif segment.type == POWER_LINE || segment.type == STEERING_LINE
+        end
+        (in_pulley > 1) && error("Bridle segment number $(segment.idx) is part of
+            $in_pulley pulleys, and should be part of either 0 or 1 pulleys.")
+
+        #TODO: Segments cannot be part of a tether if they are part of a pulley.
+        if in_pulley == 0
             in_tether = 0
             for tether in tethers
                 if segment.idx in tether.segment_idxs # each tether segment has to be part of exactly one tether
@@ -396,42 +457,32 @@ function force_eqs!(s, system, eqs, defaults, guesses;
                             in_winch += 1
                         end
                     end
-                    !(in_winch in [0,1]) && error("Tether number $(tether.idx) is part of
-                        $(in_winch) winches, and should be part of exactly 0 or 1 winch.")
-
-                    if in_winch == 1
-                        eqs = [
-                            eqs
-                            l0[segment.idx] ~ tether_length[winch_idx] / length(tether.segment_idxs)
-                        ]
-                    elseif in_winch == 0
-                        eqs = [
-                            eqs
-                            l0[segment.idx] ~ segment.l0
-                        ]
-                    end
+                    (in_winch != 1) && error("Tether number $(tether.idx) is connected to
+                        $(in_winch) winches, and should have 1 winch connected.")
+
+                    eqs = [
+                        eqs
+                        l0[segment.idx] ~ tether_length[winch_idx] / length(tether.segment_idxs)
+                    ]
                     in_tether += 1
                 end
             end
-            (in_tether != 1) && error("Segment number $(segment.idx) is part of 
-                $in_tether tethers, and should be part of exactly 1 tether.")
-        else
-            error("Unknown segment type: $(segment.type)")
+            !(in_tether in [0,1]) && error("Segment number $(segment.idx) is part of 
+                $in_tether tethers, and should be part of exactly 0 or 1 tether.")
+            if in_tether == 0
+                eqs = [
+                    eqs
+                    l0[segment.idx] ~ get_l0(psys, segment.idx)
+                ]
+            end
         end
 
-        stiffness_m = s.set.e_tether * (segment.diameter/2)^2 * pi
+        stiffness_m = get_e_tether(pset) * (get_diameter(psys, segment.idx)/2)^2 * pi
         @parameters stiffness_frac = 0.01
         (segment.type == BRIDLE) && (stiffness_m = stiffness_frac * stiffness_m)
 
-        damping_m = (s.set.damping / s.set.c_spring) * stiffness_m
+        damping_m = (get_damping(pset) / get_c_spring(pset)) * stiffness_m
         
-        if segment.compression_frac ≈ 1.0
-            eqs = [eqs; stiffness[segment.idx]       ~ stiffness_m / len[segment.idx]]
-        else
-            eqs = [eqs; stiffness[segment.idx]       ~ ifelse(len[segment.idx] > l0[segment.idx],
-                                        stiffness_m / len[segment.idx],
-                                        segment.compression_frac * stiffness_m / len[segment.idx])]
-        end
         eqs = [
             eqs
             # spring force equations
@@ -441,6 +492,9 @@ function force_eqs!(s, system, eqs, defaults, guesses;
             rel_vel[:, segment.idx]      ~ vel[:, p1] - vel[:, p2]
             spring_vel[segment.idx]      ~ rel_vel[:, segment.idx] ⋅ unit_vec[:, segment.idx]
             damping[segment.idx]         ~ damping_m / len[segment.idx]
+            stiffness[segment.idx]       ~ ifelse(len[segment.idx] > l0[segment.idx],
+                                        stiffness_m / len[segment.idx],
+                                        get_compression_frac(psys, segment.idx) * stiffness_m / len[segment.idx])
             spring_force[segment.idx] ~  (stiffness[segment.idx] * (len[segment.idx] - l0[segment.idx]) - 
                             damping[segment.idx] * spring_vel[segment.idx])
             spring_force_vec[:, segment.idx]  ~ spring_force[segment.idx] * unit_vec[:, segment.idx]
@@ -449,12 +503,12 @@ function force_eqs!(s, system, eqs, defaults, guesses;
             height[segment.idx]          ~ max(0.0, 0.5(pos[:, p1][3] + pos[:, p2][3]))
             segment_vel[:, segment.idx]  ~ 0.5(vel[:, p1] + vel[:, p2])
             segment_rho[segment.idx]     ~ calc_rho(s.am, height[segment.idx])
-            wind_vel[:, segment.idx]     ~ AtmosphericModels.calc_wind_factor(s.am, height[segment.idx], s.set.profile_law) * wind_vec_gnd
+            wind_vel[:, segment.idx]     ~ AtmosphericModels.calc_wind_factor(s.am, max(height[segment.idx], 1e-3), s.set.profile_law) * wind_vec_gnd
             va[:, segment.idx]           ~ wind_vel[:, segment.idx] - segment_vel[:, segment.idx]
-            area[segment.idx]            ~ len[segment.idx] * segment.diameter
+            area[segment.idx]            ~ len[segment.idx] * get_diameter(psys, segment.idx)
             app_perp_vel[:, segment.idx] ~ va[:, segment.idx] - 
                                         (va[:, segment.idx] ⋅ unit_vec[:, segment.idx]) * unit_vec[:, segment.idx]
-            drag_force[:, segment.idx]   ~ (0.5 * segment_rho[segment.idx] * s.set.cd_tether * norm(va[:, segment.idx]) * 
+            drag_force[:, segment.idx]   ~ (0.5 * segment_rho[segment.idx] * get_cd_tether(pset) * norm(va[:, segment.idx]) * 
                                         area[segment.idx]) * app_perp_vel[:, segment.idx]
         ]
     end
@@ -466,12 +520,11 @@ function force_eqs!(s, system, eqs, defaults, guesses;
         pulley_force(t)[eachindex(pulleys)]
         pulley_acc(t)[eachindex(pulleys)]
     end
-    @parameters pulley_damp = 1.0
+    @parameters pulley_damp = 5.0
     for pulley in pulleys
         segment = segments[pulley.segment_idxs[1]]
-        mass_per_meter = s.set.rho_tether * π * (segment.diameter/2)^2
-        mass = pulley.sum_length * mass_per_meter
-        @assert !(mass ≈ 0.0)
+        mass_per_meter = get_rho_tether(pset) * π * (get_diameter(psys, segment.idx)/2)^2
+        mass = get_sum_length(psys, pulley.idx) * mass_per_meter
         eqs = [
             eqs
             pulley_force[pulley.idx]    ~ spring_force[pulley.segment_idxs[1]] - spring_force[pulley.segment_idxs[2]]
@@ -485,7 +538,7 @@ function force_eqs!(s, system, eqs, defaults, guesses;
             ]
             defaults = [
                 defaults
-                pulley_l0[pulley.idx] => segments[pulley.segment_idxs[1]].l0
+                pulley_l0[pulley.idx] => get_l0(psys, pulley.segment_idxs[1])
                 pulley_vel[pulley.idx] => 0
             ]
         elseif pulley.type == QUASI_STATIC
@@ -496,7 +549,7 @@ function force_eqs!(s, system, eqs, defaults, guesses;
             ]
             guesses = [
                 guesses
-                pulley_l0[pulley.idx] => segments[pulley.segment_idxs[1]].l0
+                pulley_l0[pulley.idx] => get_l0(psys, pulley.segment_idxs[1])
             ]
         else
             error("Wrong pulley type")
@@ -540,8 +593,8 @@ function force_eqs!(s, system, eqs, defaults, guesses;
         ]
         defaults = [
             defaults
-            tether_length[winch.idx] => winch.tether_length
-            tether_vel[winch.idx] => 0
+            tether_length[winch.idx] => get_tether_length(psys, winch.idx)
+            tether_vel[winch.idx] => get_tether_vel(psys, winch.idx)
         ]
     end
 
@@ -564,7 +617,7 @@ function force_eqs!(s, system, eqs, defaults, guesses;
 end
 
 """
-    wing_eqs!(s, eqs, defaults; kwargs...)
+    wing_eqs!(s, eqs, pset, defaults; kwargs...)
 
 Generate the differential equations for wing dynamics including quaternion kinematics,
 angular velocities and accelerations, and forces/moments.
@@ -587,7 +640,7 @@ angular velocities and accelerations, and forces/moments.
 # Returns
 Tuple of updated equations and defaults
 """
-function wing_eqs!(s, eqs, defaults; tether_wing_force, tether_wing_moment, aero_force_b, 
+function wing_eqs!(s, eqs, psys, pset, defaults; tether_wing_force, tether_wing_moment, aero_force_b, 
     aero_moment_b, ω_b, α_b, R_b_w, wing_pos, wing_vel, wing_acc, stabilize, fix_nonstiff
 )
     wings = s.sys_struct.wings
@@ -603,6 +656,7 @@ function wing_eqs!(s, eqs, defaults; tether_wing_force, tether_wing_moment, aero
 
         # rest: forces, moments, vectors and scalar values
         moment_b(t)[eachindex(wings), 1:3] # moment in principal frame
+        wing_mass(t)[eachindex(wings)]
     end
     @parameters ω_damp = 150
 
@@ -614,7 +668,6 @@ function wing_eqs!(s, eqs, defaults; tether_wing_force, tether_wing_moment, aero
     for wing in wings
         vsm_wing = s.vsm_wings[wing.idx]
         I_b = [vsm_wing.inertia_tensor[1,1], vsm_wing.inertia_tensor[2,2], vsm_wing.inertia_tensor[3,3]]
-        @assert !(s.set.mass ≈ 0)
         axis = sym_normalize(wing_pos[wing.idx, :])
         axis_b = R_b_w[wing.idx, :, :]' * axis
         eqs = [
@@ -658,14 +711,15 @@ function wing_eqs!(s, eqs, defaults; tether_wing_force, tether_wing_moment, aero
                     wing_acc[wing.idx, :]
                 )
             )
-            wing_acc[wing.idx, :] ~ (tether_wing_force[wing.idx, :] + R_b_w[wing.idx, :, :] * aero_force_b[wing.idx, :]) / s.set.mass
+            wing_mass[wing.idx] ~ get_set_mass(pset)
+            wing_acc[wing.idx, :] ~ (tether_wing_force[wing.idx, :] + R_b_w[wing.idx, :, :] * aero_force_b[wing.idx, :]) / wing_mass[wing.idx]
         ]
         defaults = [
             defaults
-            [Q_b_w[wing.idx, i] => wing.orient[i] for i in 1:4]
-            [ω_b[wing.idx, i] => wing.angular_vel[i] for i in 1:3]
-            [wing_pos[wing.idx, i] => wing.pos_w[i] for i in 1:3]
-            [wing_vel[wing.idx, i] => wing.vel_w[i] for i in 1:3]
+            [Q_b_w[wing.idx, i] => get_orient(psys, wing.idx)[i] for i in 1:4]
+            [ω_b[wing.idx, i] => get_angular_vel(psys, wing.idx)[i] for i in 1:3]
+            [wing_pos[wing.idx, i] => get_wing_pos_w(psys, wing.idx)[i] for i in 1:3]
+            [wing_vel[wing.idx, i] => get_wing_vel_w(psys, wing.idx)[i] for i in 1:3]
         ]
     end
     
@@ -715,20 +769,21 @@ Generate equations for scalar quantities like elevation, azimuth, heading and co
     - Course angle
     - Angular velocities and accelerations
     """
-function scalar_eqs!(s, eqs; R_b_w, wind_vec_gnd, va_wing_b, wing_pos, wing_vel, wing_acc, twist_angle, twist_ω, ω_b, α_b)
+function scalar_eqs!(s, eqs, pset; R_b_w, wind_vec_gnd, va_wing_b, wing_pos, wing_vel, wing_acc, twist_angle, twist_ω, ω_b, α_b)
     @unpack wings = s.sys_struct
-    @parameters wind_scale_gnd = s.set.v_wind
-    @parameters upwind_dir = deg2rad(s.set.upwind_dir)
+    wind_scale_gnd = get_v_wind(pset)
     @variables begin
         e_x(t)[eachindex(wings), 1:3]
         e_y(t)[eachindex(wings), 1:3]
         e_z(t)[eachindex(wings), 1:3]
         wind_vel_wing(t)[eachindex(wings), 1:3]
         va_wing(t)[eachindex(wings), 1:3]
+        upwind_dir(t)
     end
     eqs = [
         eqs
-        wind_vec_gnd ~ wind_scale_gnd * rotate_around_z([0, -1, 0], -upwind_dir)
+        upwind_dir ~ deg2rad(get_upwind_dir(pset))
+        wind_vec_gnd ~ max(wind_scale_gnd, 1e-6) * rotate_around_z([0, -1, 0], -upwind_dir)
     ]
     for wing in wings
         eqs = [
@@ -884,6 +939,8 @@ function create_sys!(s::SymbolicAWEModel, system::SystemStructure; init_va_b)
     @unpack wings, groups, winches = system
 
     @parameters begin
+        psys::SystemStructure = system
+        pset::Settings = s.set
         stabilize = false
         fix_nonstiff = false
     end
@@ -910,14 +967,14 @@ function create_sys!(s::SymbolicAWEModel, system::SystemStructure; init_va_b)
     end
 
     eqs, defaults, guesses, tether_wing_force, tether_wing_moment = 
-        force_eqs!(s, system, eqs, defaults, guesses; 
+        force_eqs!(s, system, psys, pset, eqs, defaults, guesses; 
             R_b_w, wing_pos, wing_vel, wind_vec_gnd, group_aero_moment, twist_angle, twist_ω, stabilize, set_values, fix_nonstiff)
     eqs, guesses = linear_vsm_eqs!(s, eqs, guesses; aero_force_b, aero_moment_b, group_aero_moment, init_va_b, twist_angle, va_wing_b, ω_b)
-    eqs, defaults = wing_eqs!(s, eqs, defaults; tether_wing_force, tether_wing_moment, aero_force_b, aero_moment_b, 
+    eqs, defaults = wing_eqs!(s, eqs, psys, pset, defaults; tether_wing_force, tether_wing_moment, aero_force_b, aero_moment_b, 
         ω_b, α_b, R_b_w, wing_pos, wing_vel, wing_acc, stabilize, fix_nonstiff)
-    eqs = scalar_eqs!(s, eqs; R_b_w, wind_vec_gnd, va_wing_b, wing_pos, wing_vel, wing_acc, twist_angle, twist_ω, ω_b, α_b)
+    eqs = scalar_eqs!(s, eqs, pset; R_b_w, wind_vec_gnd, va_wing_b, wing_pos, wing_vel, wing_acc, twist_angle, twist_ω, ω_b, α_b)
     
-    # te_I = (1/3 * (s.set.mass/8) * te_length^2)
+    # te_I = (1/3 * (get_set_mass(pset)/8) * te_length^2)
     # # -damping / I * ω = α_damping
     # # solve for c: (c * (k*m/s^2) / (k*m^2)) * (m/s)=m/s^2 in wolframalpha
     # # damping should be in N*m*s
diff --git a/src/symbolic_awe_model.jl b/src/symbolic_awe_model.jl
index afea5f06..2e46a863 100644
--- a/src/symbolic_awe_model.jl
+++ b/src/symbolic_awe_model.jl
@@ -25,31 +25,31 @@
     defaults::Vector{Pair} = Pair[]
     guesses::Vector{Pair} = Pair[]
 
-    set_psys::Function             = (_, _) -> nothing
-    set_set_values::Function       = (_, _) -> nothing
-    set_wind::Function         = (_, _) -> nothing
-    set_vsm::Function              = (_, _) -> nothing
-    set_unknowns::Function         = (_, _) -> nothing
-    set_nonstiff::Function         = (_, _) -> nothing
-    set_lin_vsm::Function          = (_, _) -> nothing
-    set_lin_set_values::Function   = (_, _) -> nothing
-    set_lin_unknowns::Function     = (_, _) -> nothing
-    set_stabilize::Function        = (_, _) -> nothing
+    set_psys::Union{Function, Nothing}             = nothing
+    set_set_values::Union{Function, Nothing}       = nothing
+    set_set::Union{Function, Nothing}             = nothing
+    set_vsm::Union{Function, Nothing}              = nothing
+    set_unknowns::Union{Function, Nothing}         = nothing
+    set_nonstiff::Union{Function, Nothing}         = nothing
+    set_lin_vsm::Union{Function, Nothing}          = nothing
+    set_lin_set_values::Union{Function, Nothing}   = nothing
+    set_lin_unknowns::Union{Function, Nothing}     = nothing
+    set_stabilize::Union{Function, Nothing}        = nothing
     
-    get_vsm::Function              = (_) -> nothing
-    get_set_values::Function       = (_) -> nothing
-    get_unknowns::Function         = (_) -> nothing
-    get_wing_state::Function       = (_) -> nothing
-    get_winch_state::Function      = (_) -> nothing
-    get_point_state::Function      = (_) -> nothing
-    get_y::Function                = (_) -> nothing
-    get_unstretched_length::Function = (_) -> nothing
-    get_tether_length::Function    = (_) -> nothing
-    get_wing_pos::Function         = (_) -> nothing
-    get_winch_force::Function      = (_) -> nothing
-    get_spring_force::Function     = (_) -> nothing
-    get_stabilize::Function        = (_) -> nothing
-    get_pos::Function              = (_) -> nothing
+    get_vsm::Union{Function, Nothing}              = nothing
+    get_set_values::Union{Function, Nothing}       = nothing
+    get_unknowns::Union{Function, Nothing}         = nothing
+    get_wing_state::Union{Function, Nothing}       = nothing
+    get_winch_state::Union{Function, Nothing}      = nothing
+    get_point_state::Union{Function, Nothing}      = nothing
+    get_y::Union{Function, Nothing}                = nothing
+    get_unstretched_length::Union{Function, Nothing} = nothing
+    get_tether_length::Union{Function, Nothing}    = nothing
+    get_wing_pos::Union{Function, Nothing}         = nothing
+    get_winch_force::Union{Function, Nothing}      = nothing
+    get_spring_force::Union{Function, Nothing}     = nothing
+    get_stabilize::Union{Function, Nothing}        = nothing
+    get_pos::Union{Function, Nothing}              = nothing
 end
 
 """
@@ -185,73 +185,72 @@ function update_sys_state!(ss::SysState, s::SymbolicAWEModel, zoom=1.0)
     ss.time = s.integrator.t # Use integrator time
 
     # Get the state vectors from the integrator
-    set_values, tether_length, tether_vel, winch_force = s.get_winch_state(s.integrator)
-    Q_b_w, elevation, azimuth, course, heading, twist_angle, wing_vel, aero_force_b, 
-        aero_moment_b, turn_rate, va_wing_b, wind_vel_wing = s.get_wing_state(s.integrator)
-    pos, acc, wind_vec_gnd = s.get_point_state(s.integrator)
+    if !isnothing(s.get_winch_state)
+        nw = length(s.sys_struct.winches)
+        set_values, tether_length, tether_vel, winch_force = s.get_winch_state(s.integrator)
+        ss.l_tether[1:nw] .= tether_length
+        ss.v_reelout[1:nw] .= tether_vel
+        ss.force[1:nw] .= winch_force
+        ss.set_torque[1:nw] .= set_values
+    end
+    if !isnothing(s.get_wing_state)
+        Q_b_w, elevation, azimuth, course, heading, twist_angle, wing_vel, aero_force_b, 
+            aero_moment_b, turn_rate, va_wing_b, wind_vel_wing = s.get_wing_state(s.integrator)
+        ss.orient .= Q_b_w[1, :]
+        ss.turn_rates .= turn_rate[1, :]
+        ss.elevation = elevation[1]
+        ss.azimuth = azimuth[1]
+        # Depower and Steering from twist angles
+        num_groups = length(s.sys_struct.wings[1].group_idxs)
+        ss.twist_angles[1:num_groups] .= twist_angle[1:num_groups]
+        ss.depower = rad2deg(mean(ss.twist_angles)) # Average twist for depower
+        ss.steering = rad2deg(ss.twist_angles[num_groups] - ss.twist_angles[1])
+        ss.heading = heading[1] # Use heading from MTK model
+        ss.course = course[1]
+        # Apparent Wind and Aerodynamics
+        ss.v_app = norm(va_wing_b[1, :])
+        ss.v_wind_kite .= wind_vel_wing[1, :]
+        # Calculate AoA and Side Slip from apparent wind in body frame
+        # AoA: angle between v_app projected onto xz-plane and x-axis
+        # Side Slip: angle between v_app and the xz-plane
+        if ss.v_app > 1e-6 # Avoid division by zero
+            ss.AoA = atan(va_wing_b[1, 3], va_wing_b[1, 1])
+            ss.side_slip = asin(va_wing_b[1, 2] / ss.v_app)
+        else
+            ss.AoA = 0.0
+            ss.side_slip = 0.0
+        end
+        ss.aero_force_b .= aero_force_b[1, :]
+        ss.aero_moment_b .= aero_moment_b[1, :]
+        ss.vel_kite .= wing_vel[1, :]
+        # Calculate Roll, Pitch, Yaw from Quaternion
+        q = Q_b_w[1, :]
+        # roll (x-axis rotation)
+        sinr_cosp = 2 * (q[1] * q[2] + q[3] * q[4])
+        cosr_cosp = 1 - 2 * (q[2] * q[2] + q[3] * q[3])
+        ss.roll = atan(sinr_cosp, cosr_cosp)
+        # pitch (y-axis rotation)
+        sinp = 2 * (q[1] * q[3] - q[4] * q[2])
+        if abs(sinp) >= 1
+            ss.pitch = copysign(pi / 2, sinp) # use 90 degrees if out of range
+        else
+            ss.pitch = asin(sinp)
+        end
+        # yaw (z-axis rotation)
+        siny_cosp = 2 * (q[1] * q[4] + q[2] * q[3])
+        cosy_cosp = 1 - 2 * (q[3] * q[3] + q[4] * q[4])
+        ss.yaw = atan(siny_cosp, cosy_cosp)
+    end
 
+    pos, acc, wind_vec_gnd = s.get_point_state(s.integrator)
     P = length(s.sys_struct.points)
     for i in 1:P
         ss.X[i] = pos[1, i] * zoom
         ss.Y[i] = pos[2, i] * zoom
         ss.Z[i] = pos[3, i] * zoom
     end
-
-    # --- Populate SysState fields ---
     ss.acc = norm(acc) # Use the norm of the wing's acceleration vector
-    ss.orient .= Q_b_w[1, :]
-    ss.turn_rates .= turn_rate[1, :]
-    ss.elevation = elevation[1]
-    ss.azimuth = azimuth[1]
-
-    # Handle potential size mismatch for tether/winch related arrays
-    num_winches = length(s.sys_struct.winches)
-    ss.l_tether[1:num_winches] .= tether_length
-    ss.v_reelout[1:num_winches] .= tether_vel
-    ss.force[1:num_winches] .= winch_force
-
-    # Depower and Steering from twist angles
-    num_groups = length(s.sys_struct.wings[1].group_idxs)
-    ss.twist_angles[1:num_groups] .= twist_angle[1:num_groups]
-    ss.depower = rad2deg(mean(ss.twist_angles)) # Average twist for depower
-    ss.steering = rad2deg(ss.twist_angles[num_groups] - ss.twist_angles[1])
-    ss.heading = heading[1] # Use heading from MTK model
-    ss.course = course[1]
-    # Apparent Wind and Aerodynamics
-    ss.v_app = norm(va_wing_b[1, :])
     ss.v_wind_gnd .= wind_vec_gnd
-    ss.v_wind_kite .= wind_vel_wing[1, :]
-    # Calculate AoA and Side Slip from apparent wind in body frame
-    # AoA: angle between v_app projected onto xz-plane and x-axis
-    # Side Slip: angle between v_app and the xz-plane
-    if ss.v_app > 1e-6 # Avoid division by zero
-        ss.AoA = atan(va_wing_b[1, 3], va_wing_b[1, 1])
-        ss.side_slip = asin(va_wing_b[1, 2] / ss.v_app)
-    else
-        ss.AoA = 0.0
-        ss.side_slip = 0.0
-    end
-    ss.aero_force_b .= aero_force_b[1, :]
-    ss.aero_moment_b .= aero_moment_b[1, :]
-    ss.vel_kite .= wing_vel[1, :]
-    # Calculate Roll, Pitch, Yaw from Quaternion
-    q = Q_b_w[1, :]
-    # roll (x-axis rotation)
-    sinr_cosp = 2 * (q[1] * q[2] + q[3] * q[4])
-    cosr_cosp = 1 - 2 * (q[2] * q[2] + q[3] * q[3])
-    ss.roll = atan(sinr_cosp, cosr_cosp)
-    # pitch (y-axis rotation)
-    sinp = 2 * (q[1] * q[3] - q[4] * q[2])
-    if abs(sinp) >= 1
-        ss.pitch = copysign(pi / 2, sinp) # use 90 degrees if out of range
-    else
-        ss.pitch = asin(sinp)
-    end
-    # yaw (z-axis rotation)
-    siny_cosp = 2 * (q[1] * q[4] + q[2] * q[3])
-    cosy_cosp = 1 - 2 * (q[3] * q[3] + q[4] * q[4])
-    ss.yaw = atan(siny_cosp, cosy_cosp)
-    ss.set_torque[1:3] .= set_values
     nothing
 end
 
@@ -411,7 +410,6 @@ function reinit!(
     isnothing(s.sys_struct) && error("SystemStructure not defined")
 
     # init_Q_b_w, R_b_w, init_va_b = initial_orient(s)
-    init!(s.sys_struct, s.set)
     
     if isnothing(s.prob) || reload
         model_path = joinpath(KiteUtils.get_data_path(), get_model_name(s.set; precompile))
@@ -446,6 +444,7 @@ function reinit!(
         prn && @info "Initialized integrator in $t seconds"
     end
 
+    init!(s.sys_struct, s.set)
     init_unknowns_vec!(s, s.sys_struct, s.unknowns_vec)
     s.set_unknowns(s.integrator, s.unknowns_vec)
     s.set_psys(s.integrator, s.sys_struct)
@@ -533,8 +532,8 @@ function generate_getters!(s, sym_vec)
     
     set_psys = setp(sys, sys.psys)
     s.set_psys = (integ, val) -> set_psys(integ, val)
-    set_wind = setp(sys, [sys.wind_scale_gnd, sys.upwind_dir])
-    s.set_wind = (integ, val) -> set_wind(integ, val)
+    set_set = setp(sys, sys.pset)
+    s.set_set = (integ, val) -> set_set(integ, val)
     set_unknowns = setu(sys, sym_vec)
     s.set_unknowns = (integ, val) -> set_unknowns(integ, val)
     set_nonstiff = setu(sys, get_nonstiff_unknowns(s.sys_struct, s.sys))
@@ -835,7 +834,9 @@ Set ground wind speed (m/s) and upwind direction (radians). Direction: 0=north,
 π=zouth, -π/2=west (default).
 """
 function set_v_wind_ground!(s::SymbolicAWEModel, v_wind_gnd=s.set.v_wind, upwind_dir=-pi/2)
-    s.set_wind(s.integrator, [v_wind_gnd, upwind_dir])
+    s.set.v_wind = v_wind_gnd
+    s.set.upwind_dir = rad2deg(upwind_dir)
+    s.set_set(s.integrator, s.set)
     return nothing
 end
 
diff --git a/src/system_structure.jl b/src/system_structure.jl
index df855531..2f7166bf 100644
--- a/src/system_structure.jl
+++ b/src/system_structure.jl
@@ -52,19 +52,20 @@ A point mass.
 
 $(TYPEDFIELDS)
 """
-struct Point
-    idx::Int16
-    transform_idx::Int16 # idx of wing used for initial orientation
-    wing_idx::Int16
-    pos_cad::KVec3
-    pos_b::KVec3 # pos relative to wing COM in body frame
-    pos_w::KVec3 # pos in world frame
-    vel_w::KVec3 # vel in world frame
-    type::DynamicsType
+mutable struct Point
+    const idx::Int16
+    const transform_idx::Int16 # idx of wing used for initial orientation
+    const wing_idx::Int16
+    const pos_cad::KVec3
+    const pos_b::KVec3 # pos relative to wing COM in body frame
+    const pos_w::KVec3 # pos in world frame
+    const vel_w::KVec3 # vel in world frame
+    const type::DynamicsType
+    mass::SimFloat
 end
 
 """
-    Point(idx, pos_cad, type; wing_idx=1, vel_w=zeros(KVec3), transform_idx=1)
+    Point(idx, pos_cad, type; wing_idx=1, vel_w=zeros(KVec3), transform_idx=1, mass=0.0)
 
 Constructs a Point object. A point can be of four different [`DynamicsType`](@ref)s:
 - `STATIC`: the point doesn't move. ``\\ddot{\\mathbf{r}} = \\mathbf{0}``
@@ -100,8 +101,8 @@ To create a Point:
     point = Point(1, [1.0, 2.0, 3.0], DYNAMIC; wing_idx=1)
 ```
 """
-function Point(idx, pos_cad, type; wing_idx=1, vel_w=zeros(KVec3), transform_idx=1)
-    Point(idx, transform_idx, wing_idx, pos_cad, zeros(KVec3), zeros(KVec3), vel_w, type)
+function Point(idx, pos_cad, type; wing_idx=1, vel_w=zeros(KVec3), transform_idx=1, mass=0.0)
+    Point(idx, transform_idx, wing_idx, pos_cad, zeros(KVec3), zeros(KVec3), vel_w, type, mass)
 end
 
 """
@@ -139,7 +140,7 @@ The governing equation is:
 \\end{aligned}
 ```
 
-![System Overview](group-slice.svg)
+![System Overview](group_slice.svg)
 
 where:
 - ``\\tau`` is the total torque about the twist axis
@@ -374,12 +375,12 @@ mutable struct Winch
     const idx::Int16
     const model::AbstractWinchModel
     const tether_idxs::Vector{Int16}
-    tether_length::SimFloat
+    tether_length::Union{SimFloat, Nothing}
     tether_vel::SimFloat
 end
 
 """
-    Winch(idx, model, tether_idxs, tether_length; tether_vel=0.0)
+    Winch(idx, model, tether_idxs; tether_length=nothing, tether_vel=0.0)
 
 Constructs a Winch object that controls tether length through torque or speed regulation.
 
@@ -405,10 +406,10 @@ see the [WinchModels.jl documentation](https://github.com/aenarete/WinchModels.j
 - `idx::Int16`: Unique identifier for the winch.
 - `model::AbstractWinchModel`: The winch model (TorqueControlledMachine, AsyncMachine, etc.).
 - `tether_idxs::Vector{Int16}`: Vector containing the indices of the tethers connected to this winch.
-- `tether_length::SimFloat`: Initial tether length.
 
 # Keyword Arguments
 - `tether_vel::SimFloat=0.0`: Initial tether velocity (reel-out rate).
+- `tether_length::SimFloat`: Initial tether length.
 
 # Returns
 - `Winch`: A new Winch object.
@@ -419,7 +420,7 @@ To create a Winch:
     winch = Winch(1, TorqueControlledMachine(set), [1, 2], 100.0)
 ```
 """
-function Winch(idx, model, tether_idxs, tether_length; tether_vel=0.0)
+function Winch(idx, model, tether_idxs; tether_length=0.0, tether_vel=0.0)
     return Winch(idx, model, tether_idxs, tether_length, tether_vel)
 end
 
@@ -454,6 +455,7 @@ struct Wing
     group_idxs::Vector{Int16}
     transform_idx::Int16
     R_b_c::Matrix{SimFloat}
+    R_b_w::Matrix{SimFloat}
     angular_vel::KVec3
     pos_w::KVec3
     pos_cad::KVec3
@@ -461,7 +463,7 @@ struct Wing
 end
 function Base.getproperty(wing::Wing, sym::Symbol)
     if sym == :orient
-        return rotation_matrix_to_quaternion(wing.R_b_c)
+        return rotation_matrix_to_quaternion(wing.R_b_w)
     else
         return getfield(wing, sym)
     end
@@ -537,7 +539,7 @@ Create a wing with identity orientation and two attached groups:
 """
 function Wing(idx, group_idxs, R_b_c, pos_cad; transform_idx=1, angular_vel=zeros(KVec3), 
         pos_w=zeros(KVec3), vel_w=zeros(KVec3))
-    return Wing(idx, group_idxs, transform_idx, R_b_c, angular_vel, pos_w, pos_cad, vel_w)
+    return Wing(idx, group_idxs, transform_idx, R_b_c, zeros(3,3), angular_vel, pos_w, pos_cad, vel_w)
 end
 
 """
@@ -729,6 +731,7 @@ function SystemStructure(name, set;
     end
     for (i, segment) in enumerate(segments)
         @assert segment.idx == i
+        (segment.l0 ≈ 0) && (segment.l0 = norm(points[segment.point_idxs[1]].pos_cad - points[segment.point_idxs[2]].pos_cad))
     end
     for (i, pulley) in enumerate(pulleys)
         @assert pulley.idx == i
@@ -738,6 +741,11 @@ function SystemStructure(name, set;
     end
     for (i, winch) in enumerate(winches)
         @assert winch.idx == i
+        if iszero(winch.tether_length)
+            for segment_idx in tethers[winch.tether_idxs[1]].segment_idxs
+                winch.tether_length += segments[segment_idx].l0
+            end
+        end
         set.l_tethers[i]   = winch.tether_length
         set.v_reel_outs[i] = winch.tether_vel
     end
@@ -778,6 +786,11 @@ function SystemStructure(set::Settings, wing::RamAirWing)
     end
 end
 
+function apply_heading(vec, R_t_w, curr_R_t_w, heading)
+    vec_along_z = rotate_around_z(curr_R_t_w' * vec, heading)
+    return R_t_w * vec_along_z
+end
+
 function init!(transforms::Vector{Transform}, sys_struct::SystemStructure)
     @unpack points, wings = sys_struct
     for transform in transforms
@@ -805,8 +818,7 @@ function init!(transforms::Vector{Transform}, sys_struct::SystemStructure)
         for point in points
             if point.transform_idx == transform.idx
                 vec = point.pos_w - base_pos
-                vec_along_z = rotate_around_z(curr_R_t_w' * vec, transform.heading)
-                point.pos_w .= base_pos + R_t_w * vec_along_z
+                point.pos_w .= base_pos + apply_heading(vec, R_t_w, curr_R_t_w, transform.heading)
             end
             if point.type == WING
                 wing = wings[point.wing_idx]
@@ -816,10 +828,9 @@ function init!(transforms::Vector{Transform}, sys_struct::SystemStructure)
         for wing in wings
             if wing.transform_idx == transform.idx
                 vec = wing.pos_w - base_pos
-                vec_along_z = rotate_around_x(curr_R_t_w' * vec, transform.heading)
-                wing.pos_w .= base_pos + R_t_w * vec_along_z
+                wing.pos_w .= base_pos + apply_heading(vec, R_t_w, curr_R_t_w, transform.heading)
                 for i in 1:3
-                    wing.R_b_c[:, i] .= R_t_w * rotate_around_x(curr_R_t_w' * wing.R_b_c[:, i], transform.heading)
+                    wing.R_b_w[:, i] .= apply_heading(wing.R_b_c[:, i], R_t_w, curr_R_t_w, transform.heading)
                 end
             end
         end
@@ -974,9 +985,9 @@ function create_ram_sys_struct(set::Settings, vsm_wing::RamAirWing)
     points, segments, tethers, left_steering_idx = create_tether(3, set, points, segments, tethers, attach_points[2], STEERING_LINE, dynamics_type, z)
     points, segments, tethers, right_steering_idx = create_tether(4, set, points, segments, tethers, attach_points[4], STEERING_LINE, dynamics_type, z)
 
-    winches = [winches; Winch(1, TorqueControlledMachine(set), [left_power_idx, right_power_idx], set.l_tether)]
-    winches = [winches; Winch(2, TorqueControlledMachine(set), [left_steering_idx], set.l_tether)]
-    winches = [winches; Winch(3, TorqueControlledMachine(set), [right_steering_idx], set.l_tether)]
+    winches = [winches; Winch(1, TorqueControlledMachine(set), [left_power_idx, right_power_idx])]
+    winches = [winches; Winch(2, TorqueControlledMachine(set), [left_steering_idx])]
+    winches = [winches; Winch(3, TorqueControlledMachine(set), [right_steering_idx])]
 
     wings = [Wing(1, [1,2,3,4], I(3), zeros(3))]
     transforms = [Transform(1, deg2rad(set.elevation), deg2rad(set.azimuth), deg2rad(set.heading);
diff --git a/test/test_ram_air_kite.jl b/test/test_ram_air_kite.jl
index 6d84958c..012eca36 100644
--- a/test/test_ram_air_kite.jl
+++ b/test/test_ram_air_kite.jl
@@ -155,12 +155,11 @@ const BUILD_SYS = true
 
             # Check if wind speed and direction have been updated correctly
             @test norm(s.integrator[s.sys.wind_vec_gnd]) ≈ new_wind_speed
-            @test s.integrator.ps[s.sys.upwind_dir] ≈ -pi/4 # Default upwind_dir
             @test s.integrator[s.sys.wind_vec_gnd[1]] ≈ -s.integrator[s.sys.wind_vec_gnd[2]]
 
             KiteModels.set_v_wind_ground!(s, initial_wind_speed, initial_upwind_dir)
             @test s.integrator[s.sys.wind_vec_gnd[1]] ≈ initial_wind_speed
-            @test s.integrator.ps[s.sys.upwind_dir] ≈ -pi/2 # Default upwind_dir
+            @test norm(s.integrator[s.sys.wind_vec_gnd]) ≈ initial_wind_speed
         end
     end
 
@@ -325,4 +324,4 @@ end
 # Restore original data path
 set_data_path(old_path)
 
-nothing
\ No newline at end of file
+nothing