OpenSourceAWE · 1-Bart-1 · Jun 18, 2025 · Jun 14, 2025 · Jun 14, 2025 · Jun 16, 2025
diff --git a/.gitignore b/.gitignore
@@ -1,3 +1,4 @@
+LocalPreferences.toml
 *.arrow
 *.bin
 *.xopp

diff --git a/AUTHORS.md b/AUTHORS.md
@@ -7,7 +7,7 @@ SPDX-License-Identifier: MIT
 This project was started by Uwe Fechner in 2022 with the goal to make the results 
 of his PhD thesis on the simulation and control of kite power systems, re-coded in a fast and modern programming language, available to the public.
 
-Bart van de Lint joined the project in 2024. His major contribution is the first modelling toolkit (MTK) based model, a ram-air kite model.
+Bart van de Lint joined the project in 2024. His major contribution is the first ModelingToolkit (MTK) based model, a ram-air kite model.
 
 ## Developers
 - Uwe Fechner, Delft/ Den Haag, The Netherlands -  original author of the `KPS3` and `KPS4` kite models

diff --git a/Manifest-v1.10.toml.default b/Manifest-v1.10.toml.default
@@ -2,7 +2,7 @@
 
 julia_version = "1.10.9"
 manifest_format = "2.0"
-project_hash = "e5116e83aee3763e3ebd89b0de5522b585f68fec"
+project_hash = "098d81951e8c706049a0f77aa39c5c7b3d03fe19"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "e2478490447631aedba0823d4d7a80b2cc8cdb32"
@@ -611,9 +611,9 @@ uuid = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
 version = "1.0.5"
 
 [[deps.EnzymeCore]]
-git-tree-sha1 = "7d7822a643c33bbff4eab9c87ca8459d7c688db0"
+git-tree-sha1 = "8272a687bca7b5c601c0c24fc0c71bff10aafdfd"
 uuid = "f151be2c-9106-41f4-ab19-57ee4f262869"
-version = "0.8.11"
+version = "0.8.12"
 weakdeps = ["Adapt"]
 
     [deps.EnzymeCore.extensions]
@@ -786,9 +786,9 @@ uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
 version = "0.1.5"
 
 [[deps.InlineStrings]]
-git-tree-sha1 = "6a9fde685a7ac1eb3495f8e812c5a7c3711c2d5e"
+git-tree-sha1 = "8594fac023c5ce1ef78260f24d1ad18b4327b420"
 uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48"
-version = "1.4.3"
+version = "1.4.4"
 weakdeps = ["ArrowTypes", "Parsers"]
 
     [deps.InlineStrings.extensions]
@@ -889,9 +889,9 @@ version = "0.3.8"
 
 [[deps.KiteUtils]]
 deps = ["Arrow", "CSV", "DocStringExtensions", "LinearAlgebra", "OrderedCollections", "Parameters", "Parsers", "Pkg", "PrecompileTools", "RecursiveArrayTools", "ReferenceFrameRotations", "Rotations", "StaticArrays", "StructArrays", "StructTypes", "YAML"]
-git-tree-sha1 = "ef1ab34bdb8c1bcf06f99b4364cae5df1be6936a"
+git-tree-sha1 = "93659cbb96514bafd053ae0e81d08d15248040e9"
 uuid = "90980105-b163-44e5-ba9f-8b1c83bb0533"
-version = "0.10.10"
+version = "0.10.11"
 
 [[deps.Krylov]]
 deps = ["LinearAlgebra", "Printf", "SparseArrays"]

diff --git a/Manifest-v1.11.toml.default b/Manifest-v1.11.toml.default
@@ -2,7 +2,7 @@
 
 julia_version = "1.11.5"
 manifest_format = "2.0"
-project_hash = "f9f7424be55e7c582f860e66512280fdde05365c"
+project_hash = "d42173c595cda511f92896765dacbdae834092a9"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "e2478490447631aedba0823d4d7a80b2cc8cdb32"

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
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
diff --git a/docs/src/ram_air_kite.md b/docs/src/ram_air_kite.md
@@ -18,6 +18,8 @@ Segment
 Tether
 Winch
 Group
+Wing
+Transform
 SystemStructure
 ```
 

diff --git a/docs/src/tether_during_sim.png b/docs/src/tether_during_sim.png
diff --git a/docs/src/tether_sys_struct.png b/docs/src/tether_sys_struct.png
diff --git a/docs/src/tutorial_system_structure.md b/docs/src/tutorial_system_structure.md
@@ -20,13 +20,16 @@ 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, [0.0, 0.0, set.l_tether], STATIC; wing_idx=0))
+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 rigid 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.
@@ -35,28 +38,35 @@ segment_idxs = Int[]
 for i in 1:set.segments
     global points, segments
     point_idx = i+1
-    pos = [0.0, 0.0, set.l_tether] - set.l_tether / set.segments * [0.0, 0.0, i]
+    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
 ```
-From these arrays of points and segments we can create a [`SystemStructure`](@ref), which can be plotted in 2d to quickly investigate if the model is correct.
+
+In order to describe the initial orientation of the structure, we define a [`Transform`](@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, [0.0, 0.0, 50.0], points[1].idx; 
+                            rot_point_idx=points[end].idx)]
+```
+
+From the points, segments and transform we create a [`SystemStructure`](@ref), which can be plotted in 2d to quickly investigate if the model is correct.
 ```julia
-sys_struct = SystemStructure("tether"; points, segments)
+sys_struct = SystemStructure("tether", set; points, segments, transforms)
 plot(sys_struct, 0.0)
-plt.gcf()
 ```
+![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
-model = SymbolicAWEModel(set, sys_struct)
+sam = SymbolicAWEModel(set, sys_struct)
 
-init_sim!(model; remake=false)
-for i in 1:100
-    plot(model, i/set.sample_freq)
-    next_step!(model; dt=1/set.sample_freq)
+init_sim!(sam; remake=false)
+for i in 1:80
+    plot(sam, i/set.sample_freq)
+    next_step!(sam)
 end
-plt.gcf()
 ```
+![Tether during simulation](tether_during_sim.png)
diff --git a/examples/ram_air_kite.jl b/examples/ram_air_kite.jl
@@ -5,7 +5,7 @@ using Timers
 tic()
 @info "Loading packages "
 
-PLOT = false
+PLOT = true
 using Pkg
 if ! ("LaTeXStrings" ∈ keys(Pkg.project().dependencies))
     using TestEnv; TestEnv.activate()
@@ -47,8 +47,8 @@ toc()
 # plot(sam.sys_struct, 0.0; zoom=false, front=false)
 
 # Initialize at elevation
-sam.sys_struct.winches[2].tether_length += 0.2
-sam.sys_struct.winches[3].tether_length += 0.2
+set.l_tethers[2] += 0.2
+set.l_tethers[3] += 0.2
 init_sim!(sam; remake=false, reload=false)
 sys = sam.sys
 
@@ -63,7 +63,7 @@ sys_state = SysState(sam)
 t = 0.0
 runtime = 0.0
 integ_runtime = 0.0
-bias = set.quasi_static ? 0.45 : 0.35
+bias = set.quasi_static ? 0.45 : 0.40
 t0 = sam.integrator.t
 
 try
@@ -147,3 +147,5 @@ p = plotx(sl.time,
 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)
diff --git a/examples/tether.jl b/examples/tether.jl
@@ -11,26 +11,27 @@ dynamics_type = DYNAMIC
 points = Point[]
 segments = Segment[]
 
-points = push!(points, Point(1, [0.0, 0.0, set.l_tether], STATIC; wing_idx=0))
+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, set.l_tether] - set.l_tether / set.segments * [0.0, 0.0, i]
+    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
 
-sys_struct = SystemStructure("tether"; points, segments)
+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)
 
-model = SymbolicAWEModel(set, sys_struct)
+sam = SymbolicAWEModel(set, sys_struct)
 
-init_sim!(model; remake=false)
-for i in 1:100
-    plot(model, i/set.sample_freq)
-    next_step!(model)
-end
+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
@@ -7,7 +7,7 @@ using ControlPlots, KiteModels
 
 export plot
 
-function ControlPlots.plot(sys::SystemStructure, reltime; l_tether=50.0, wing_pos=nothing, e_z=zeros(3), zoom=false, front=false)
+function ControlPlots.plot(sys::SystemStructure, reltime; l_tether=50.0, wing_pos=nothing, e_z=zeros(3), zoom=false, front=false, xy=nothing)
     pos = [sys.points[i].pos_w for i in eachindex(sys.points)]
     !isnothing(wing_pos) && (pos = [pos..., wing_pos...])
     seg = [[sys.segments[i].point_idxs[1], sys.segments[i].point_idxs[2]] for i in eachindex(sys.segments)]
@@ -24,7 +24,7 @@ function ControlPlots.plot(sys::SystemStructure, reltime; l_tether=50.0, wing_po
         xlim = (-12.5 - 0.5l_tether, 12.5 + 0.5l_tether)
         ylim = (-5, l_tether+10)
     end
-    ControlPlots.plot2d(pos, seg, reltime; zoom, front, xlim, ylim, dz_zoom=0.6)
+    ControlPlots.plot2d(pos, seg, reltime; zoom, front, xlim, ylim, dz_zoom=0.6, xy)
 end
 
 function ControlPlots.plot(s::SymbolicAWEModel, reltime; kwargs...)

diff --git a/src/KiteModels.jl b/src/KiteModels.jl
@@ -47,7 +47,7 @@ export calculate_rotational_inertia!
 export kite_ref_frame, orient_euler, spring_forces, upwind_dir, copy_model_settings, menu2
 export create_ram_sys_struct, create_simple_ram_sys_struct
 import LinearAlgebra: norm
-export SystemStructure, Point, Group, Segment, Pulley, Tether, Winch, Wing
+export SystemStructure, Point, Group, Segment, Pulley, Tether, Winch, Wing, Transform
 export DynamicsType, DYNAMIC, QUASI_STATIC, WING, STATIC
 export SegmentType, POWER, STEERING, BRIDLE
 
@@ -107,7 +107,7 @@ end
 
 include("KPS4.jl") # include code, specific for the four point kite model
 include("system_structure.jl")
-include("symbolic_awe_system.jl") # include code, specific for the ram air kite model
+include("symbolic_awe_model.jl") # include code, specific for the ram air kite model
 include("mtk_model.jl")
 include("KPS3.jl") # include code, specific for the one point kite model
 include("init.jl") # functions to calculate the initial state vector, the initial masses and initial springs

diff --git a/src/init.jl b/src/init.jl
@@ -144,6 +144,15 @@ function init(s::KPS4, X=zeros(2 * (s.set.segments+KITE_PARTICLES-1)+1); old=fal
     MVector{6*(s.set.segments+KITE_PARTICLES)+2, SimFloat}(res1), MVector{6*(s.set.segments+KITE_PARTICLES)+2, SimFloat}(res2)
 end
 
+# rotate a 3d vector around the x axis in the yz plane - following the right hand rule
+function rotate_around_x(vec, angle::T) where T
+    result = zeros(T, 3)
+    result[1] = vec[1]
+    result[2] = cos(angle) * vec[2] - sin(angle) * vec[3]
+    result[3] = sin(angle) * vec[2] + cos(angle) * vec[3]
+    result
+end
+
 # rotate a 3d vector around the y axis in the xz plane - following the right hand rule
 function rotate_around_y(vec, angle::T) where T
     result = zeros(T, 3)
@@ -161,9 +170,3 @@ function rotate_around_z(vec, angle::T) where T
     result[3] = vec[3]
     result
 end
-
-function rotate_v_around_k(v, k, θ)
-    k = normalize(k)
-    v_rot = v * cos(θ) + (k × v) * sin(θ)  + k * (k ⋅ v) * (1 - cos(θ))
-    return v_rot
-end
diff --git a/src/mtk_model.jl b/src/mtk_model.jl
@@ -65,6 +65,14 @@ 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
+end
+@register_array_symbolic get_pos_w(sys::SystemStructure, point_idx::Int16) begin
+    size=size(KVec3)
+    eltype=SimFloat
+end
+
 """
     force_eqs!(s, system, eqs, defaults, guesses; kwargs...)
 
@@ -103,6 +111,7 @@ 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)]
@@ -136,6 +145,8 @@ function force_eqs!(s, system, eqs, defaults, guesses;
                 end
             end
         end
+        @assert !iszero(mass)
+
         eqs = [
             eqs
             point_mass[point.idx] ~ mass
@@ -170,7 +181,7 @@ function force_eqs!(s, system, eqs, defaults, guesses;
                 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
-            else
+            elseif found == 0
                 pos_b = point.pos_b
                 eqs = [
                     eqs
@@ -189,7 +200,7 @@ function force_eqs!(s, system, eqs, defaults, guesses;
         elseif point.type == STATIC
             eqs = [
                 eqs
-                pos[:, point.idx]    ~ point.pos_w
+                pos[:, point.idx]    ~ get_pos_w(psys, point.idx)
                 vel[:, point.idx]    ~ zeros(3)
                 acc[:, point.idx]    ~ zeros(3)
             ]
@@ -214,7 +225,7 @@ function force_eqs!(s, system, eqs, defaults, guesses;
             ]
             defaults = [
                 defaults
-                [pos[j, point.idx] => point.pos_w[j] for j in 1:3]
+                [pos[j, point.idx] => get_pos_w(psys, point.idx)[j] for j in 1:3]
                 [vel[j, point.idx] => 0 for j in 1:3]
             ]
         elseif point.type == QUASI_STATIC
@@ -227,7 +238,7 @@ function force_eqs!(s, system, eqs, defaults, guesses;
             guesses = [
                 guesses
                 [acc[j, point.idx] => 0 for j in 1:3]
-                [pos[j, point.idx] => point.pos_w[j] for j in 1:3]
+                [pos[j, point.idx] => get_pos_w(psys, point.idx)[j] for j in 1:3]
                 [point_force[j, point.idx] => 0 for j in 1:3]
             ]
         else
@@ -653,14 +664,20 @@ function wing_eqs!(s, eqs, defaults; tether_wing_force, tether_wing_moment, aero
             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[i] for i in 1:3]
-            [wing_vel[wing.idx, i] => wing.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]
         ]
     end
 
     return eqs, defaults
 end
 
+function rotate_v_around_k(v, k, θ)
+    k = sym_normalize(k)
+    v_rot = v * cos(θ) + (k × v) * sin(θ)  + k * (k ⋅ v) * (1 - cos(θ))
+    return v_rot
+end
+
 function calc_R_v_w(wing_pos, e_x)
     z = sym_normalize(wing_pos)
     y = sym_normalize(z × e_x)
@@ -767,12 +784,12 @@ function scalar_eqs!(s, eqs; R_b_w, wind_vec_gnd, va_wing_b, wing_pos, wing_vel,
             distance_vel[wing.idx]    ~ wing_vel[wing.idx, :] ⋅ R_v_w[wing.idx, :, 3]
             distance_acc[wing.idx]    ~ wing_acc[wing.idx, :] ⋅ R_v_w[wing.idx, :, 3]
 
-            elevation[wing.idx]           ~ atan(z / x)
+            elevation[wing.idx]           ~ KiteUtils.calc_elevation(wing_pos[wing.idx, :])
             # elevation_vel = d/dt(atan(z/x)) = (x*ż' - z*ẋ')/(x^2 + z^2) according to wolframalpha
             elevation_vel[wing.idx]       ~ (x*z´ - z*x´) / 
                                     (x^2 + z^2)
             elevation_acc[wing.idx]       ~ ((x^2 + z^2)*(x*z´´ - z*x´´) + 2(z*x´ - x*z´)*(x*x´ + z*z´))/(x^2 + z^2)^2
-            azimuth[wing.idx]             ~ atan(y / x)
+            azimuth[wing.idx]             ~ -KiteUtils.azimuth_east(wing_pos[wing.idx, :])
             # azimuth_vel = d/dt(atan(y/x)) = (-y*x´ + x*y´)/(x^2 + y^2) # TODO: check if correct
             azimuth_vel[wing.idx]         ~ (-y*x´ + x*y´) / 
                                     (x^2 + y^2)
@@ -861,8 +878,8 @@ end
 
 function create_sys!(s::SymbolicAWEModel, system::SystemStructure; init_va_b)
     eqs = []
-    defaults = Pair{Num, Real}[]
-    guesses = Pair{Num, Real}[]
+    defaults = Pair{Num, Any}[]
+    guesses = Pair{Num, Any}[]
 
     @unpack wings, groups, winches = system
-Original file line number
+Diff line change
@@ -1,3 +1,4 @@
+    LocalPreferences.toml
     *.arrow
     *.bin
     *.xopp
@@ Expand Down @@