Fix ram model

.gitignore

@@ -1,3 +1,4 @@
+*.arrow
 *.bin
 *.xopp
 *.so

REUSE.toml

@@ -1,49 +1,31 @@
 version = 1
 
 [[annotations]]
-path = "data/prob*"
-
-SPDX-FileCopyrightText = "2025 Bart van de Lint"
-SPDX-License-Identifier = "MPL-2.0"
-
-[[annotations]]
-path = "data/ram_air_kite*"
-
+path = ["data/prob*", "data/ram*", "examples/discrete_jacobian.jl"]
 SPDX-FileCopyrightText = "2025 Bart van de Lint"
 SPDX-License-Identifier = "MPL-2.0"
 
 [[annotations]]
-path = "Manifest*"
-
-SPDX-FileCopyrightText = "2025 Uwe Fechner, Bart van de Lint"
-SPDX-License-Identifier = "MIT"
-
-[[annotations]]
-path = "Project.toml"
-
-SPDX-FileCopyrightText = "2025 Uwe Fechner, Bart van de Lint"
-SPDX-License-Identifier = "MIT"
+path = "data/*ram.yaml"
+SPDX-FileCopyrightText = "2025 Bart van de Lint, Uwe Fechner"
+SPDX-License-Identifier = "CC-BY-4.0"
 
 [[annotations]]
-path = "CITATION.cff"
-
+path = ["data/settings.yaml", "data/system.yaml", "data/vsm_settings.yaml"]
 SPDX-FileCopyrightText = "2025 Uwe Fechner"
-SPDX-License-Identifier = "MIT"
+SPDX-License-Identifier = "CC-BY-4.0"
 
 [[annotations]]
-path = ".gitignore"
-
+path = ["Manifest*", "Project.toml", ".gitignore"]
 SPDX-FileCopyrightText = "2025 Uwe Fechner, Bart van de Lint"
 SPDX-License-Identifier = "MIT"
 
 [[annotations]]
 path = ["docs/src/*.png", "docs/src/*.md", "docs/src/kite_power_tools.drawio"]
-
 SPDX-FileCopyrightText = "2025 Uwe Fechner"
 SPDX-License-Identifier = "CC-BY-4.0"
 
 [[annotations]]
-path = ["docs/.gitignore", "docs/Project.toml"]
-
+path = ["docs/.gitignore", "docs/Project.toml", "CITATION.cff"]
 SPDX-FileCopyrightText = "2025 Uwe Fechner"
 SPDX-License-Identifier = "MIT"

docs/src/examples_ram_air.md

@@ -73,4 +73,15 @@ SIMPLE=true; include("examples/ram_air_kite.jl")
 The simple model has a very simple bridle system without pulleys and with less attachment points on the wing. 
 While the default model has a [speed system](https://kiteboarding.com/proddetail.asp?prod=ozone-r1v4-pro-tune-speedsystem-complete) with pulleys and more attachment points on the wing.
 
-![Oscillating steering input response, simple system](oscillating_steering_simple.png)
+![Oscillating steering input response, simple system](oscillating_steering_simple.png)
+
+## How to create a RamAirKite
+The following code is a minimal example that shows how to create a ram air kite struct:
+```
+using KiteModels
+
+# Initialize model
+set = load_settings("system_ram.yaml")
+
+rak = RamAirKite(set)
+```

docs/src/functions.md

@@ -68,7 +68,6 @@ calculate_rotational_inertia!
 calc_set_cl_cd!
 calc_aero_forces!
 calc_particle_forces!
-find_z_axis_point
 inner_loop!
 loop!
 ```

examples/ram_air_kite.jl

@@ -26,7 +26,7 @@ include(joinpath(@__DIR__, "plotting.jl"))
 
 # Simulation parameters
 dt = 0.05
-total_time = 10  # Longer simulation to see oscillations
+total_time = 10.0  # Longer simulation to see oscillations
 vsm_interval = 3
 steps = Int(round(total_time / dt))
 
@@ -47,25 +47,25 @@ 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_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()
+# 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
-measure.sphere_pos .= deg2rad.([70.0 70.0; 1.0 -1.0])
 KiteModels.init_sim!(s, measure; remake=false, reload=true)
 sys = s.sys
 
 @info "System initialized at:"
 toc()
 
-# # Stabilize system
+# Stabilize system
 s.integrator.ps[sys.stabilize] = true
-for i in 1:10÷dt
+for i in 1:1÷dt
     next_step!(s; dt, vsm_interval=1)
 end
 s.integrator.ps[sys.stabilize] = false
@@ -75,7 +75,8 @@ sys_state = KiteModels.SysState(s)
 t = 0.0
 runtime = 0.0
 integ_runtime = 0.0
-bias = 0.5
+bias = 0.35
+t0 = s.integrator.t
 
 try
     while t < total_time
@@ -84,7 +85,7 @@ try
         PLOT && plot(s, t; zoom=false, front=false)
 
         # Calculate steering inputs based on cosine wave
-        steering = steering_magnitude * cos(2π * steering_freq * t+bias)
+        steering = steering_magnitude * cos(2π * steering_freq * t + bias)
         set_values = -s.set.drum_radius .* s.integrator[sys.winch_force]
         _vsm_interval = 1
         if t > 1.0
@@ -93,8 +94,8 @@ try
         end
 
         # Step simulation
-        steptime = @elapsed (t_new, integ_steptime) = next_step!(s, set_values; dt, vsm_interval=_vsm_interval)
-        t = t_new - 10.0  # Adjust for initial stabilization time
+        steptime = @elapsed (t_new, integ_steptime) = next_step!(s, set_values; dt, vsm_interval=vsm_interval)
+        t = t_new - t0  # Adjust for initial stabilization time
 
         # Track performance after initial transient
         if (t > total_time/2)
@@ -105,6 +106,7 @@ try
 
         # Log state variables
         KiteModels.update_sys_state!(sys_state, s)
+        sys_state.time = t
         log!(logger, sys_state)
     end
 catch e
@@ -129,15 +131,15 @@ sl = lg.syslog
 turn_rates_deg = rad2deg.(hcat(sl.turn_rates...))
 v_reelout_23 = [sl.v_reelout[i][2] for i in eachindex(sl.v_reelout)], [sl.v_reelout[i][3] for i in eachindex(sl.v_reelout)] # Winch 2 and 3
 aero_force_z = [sl.aero_force_b[i][3] for i in eachindex(sl.aero_force_b)]
-aero_moment_y = [sl.aero_moment_b[i][2] for i in eachindex(sl.aero_moment_b)]
+aero_moment_z = [sl.aero_moment_b[i][3] for i in eachindex(sl.aero_moment_b)]
 twist_angles_deg = rad2deg.(hcat(sl.twist_angles...))
 AoA_deg = rad2deg.(sl.AoA)
 heading_deg = rad2deg.(sl.heading)
 
-p = plotx(sl.time .- 10,
+p = plotx(sl.time,
     [turn_rates_deg[1,:], turn_rates_deg[2,:], turn_rates_deg[3,:]],
     v_reelout_23,
-    [aero_force_z, aero_moment_y],
+    [aero_force_z, aero_moment_z],
     [twist_angles_deg[1,:], twist_angles_deg[2,:], twist_angles_deg[3,:], twist_angles_deg[4,:]],
     [AoA_deg],
     [heading_deg];
@@ -146,7 +148,7 @@ p = plotx(sl.time .- 10,
     labels=[
         [L"\omega_x", L"\omega_y", L"\omega_z"],
         ["v_ro[2]", "v_ro[3]"],
-        [L"F_{aero,z}", L"M_{aero,y}"],
+        [L"F_{aero,z}", L"M_{aero,z}"],
         ["twist[1]", "twist[2]", "twist[3]", "twist[4]"],
         ["AoA"],
         ["heading"]

src/mtk_model.jl

@@ -271,12 +271,13 @@ function force_eqs!(s, system, eqs, defaults, guesses;
         inertia = 1/3 * (s.set.mass/length(groups)) * (norm(group.chord))^2 # plate inertia around leading edge
         @assert !(inertia ≈ 0.0)
         @parameters twist_damp = 50
+        @parameters max_twist = deg2rad(90)
 
         eqs = [
             eqs
             group_tether_moment[group.idx] ~ sum(tether_moment[group.idx, :])
             twist_α[group.idx] ~ (group_aero_moment[group.idx] + group_tether_moment[group.idx]) / inertia
-            twist_angle[group.idx] ~ clamp(free_twist_angle[group.idx], -π/2, π/2)
+            twist_angle[group.idx] ~ clamp(free_twist_angle[group.idx], -max_twist, max_twist)
         ]
         if group.type == DYNAMIC
             eqs = [
@@ -618,6 +619,20 @@ function diff_eqs!(s, eqs, defaults; tether_kite_force, tether_kite_moment, aero
     return eqs, defaults
 end
 
+function calc_R_v_w(kite_pos, e_x)
+    z = sym_normalize(kite_pos)
+    y = sym_normalize(z × e_x)
+    x = y × z
+    return [x y z]
+end
+
+function calc_R_t_w(elevation, azimuth)
+    x = rotate_around_z(rotate_around_y([0, 0, -1], -elevation), azimuth)
+    z = rotate_around_z(rotate_around_y([1, 0, 0], -elevation), azimuth)
+    y = z × x
+    return [x y z]
+end
+
 """
     scalar_eqs!(s, eqs, measure; R_b_w, wind_vec_gnd, va_kite_b, kite_pos, kite_vel, kite_acc, twist_angle, twist_ω)
 
@@ -663,7 +678,7 @@ function scalar_eqs!(s, eqs, measure; R_b_w, wind_vec_gnd, va_kite_b, kite_pos,
         va_kite_b ~ R_b_w' * va_kite
     ]
     @variables begin
-        heading_x(t)
+        heading(t)
         turn_rate(t)[1:3]
         turn_acc(t)[1:3]
         azimuth(t)
@@ -682,9 +697,7 @@ function scalar_eqs!(s, eqs, measure; R_b_w, wind_vec_gnd, va_kite_b, kite_pos,
         simple_twist_angle(t)[1:2]
         simple_twist_ω(t)[1:2]
         R_v_w(t)[1:3, 1:3]
-        view_z(t)[1:3]
-        view_y(t)[1:3]
-        view_x(t)[1:3]
+        R_t_w(t)[1:3, 1:3]
         distance(t)
         distance_vel(t)
         distance_acc(t)
@@ -695,21 +708,18 @@ function scalar_eqs!(s, eqs, measure; R_b_w, wind_vec_gnd, va_kite_b, kite_pos,
     x´´, y´´, z´´ = kite_acc
 
     half_len = length(twist_angle)÷2
+    heading_vec = R_t_w' * R_v_w[:,1]
 
     eqs = [
         eqs
-        view_z          ~ sym_normalize(kite_pos)
-        view_y          ~ view_z × e_y
-        view_x          ~ view_y × view_z
-        R_v_w[:,1]      ~ view_x
-        R_v_w[:,2]      ~ view_y
-        R_v_w[:,3]      ~ view_z
-        heading_x       ~ calc_heading(e_x)
+        vec(R_v_w)     .~ vec(calc_R_v_w(kite_pos, e_x))
+        vec(R_t_w)     .~ vec(calc_R_t_w(elevation, azimuth))
+        heading         ~ atan(heading_vec[2], heading_vec[1])
         turn_rate       ~ R_v_w' * (R_b_w * ω_b) # Project angular velocity onto view frame
         turn_acc        ~ R_v_w' * (R_b_w * α_b)
         distance        ~ norm(kite_pos)
-        distance_vel    ~ kite_vel ⋅ view_z
-        distance_acc    ~ kite_acc ⋅ view_z
+        distance_vel    ~ kite_vel ⋅ R_v_w[:,3]
+        distance_acc    ~ kite_acc ⋅ R_v_w[:,3]
 
         elevation           ~ atan(z / x)
         # elevation_vel = d/dt(atan(z/x)) = (x*ż' - z*ẋ')/(x^2 + z^2) according to wolframalpha
@@ -721,9 +731,9 @@ function scalar_eqs!(s, eqs, measure; R_b_w, wind_vec_gnd, va_kite_b, kite_pos,
         azimuth_vel         ~ (-y*x´ + x*y´) / 
                                 (x^2 + y^2)
         azimuth_acc         ~ ((x^2 + y^2)*(-y*x´´ + x*y´´) + 2(y*x´ - x*y´)*(x*x´ + y*y´))/(x^2 + y^2)^2
+        course              ~ atan(-azimuth_vel, elevation_vel)
         x_acc               ~ kite_acc ⋅ e_x
         y_acc               ~ kite_acc ⋅ e_y
-        course              ~ atan(-azimuth_vel, elevation_vel)
 
         angle_of_attack     ~ calc_angle_of_attack(va_kite_b) + 0.5twist_angle[half_len] + 0.5twist_angle[half_len+1]
         simple_twist_angle[1] ~ sum(twist_angle[1:half_len]) / half_len
@@ -763,14 +773,8 @@ function linear_vsm_eqs!(s, eqs, guesses; aero_force_b, aero_moment_b, group_aer
     @assert length(s.point_system.groups) == length(sol.group_moment_dist)
 
     y_ = [init_va_b; zeros(length(s.point_system.groups)); zeros(3)]
-    jac_, x_ = VortexStepMethod.linearize(
-        s.vsm_solver, 
-        s.aero, 
-        y_;
-        va_idxs=1:3, 
-        omega_idxs=4:6,
-        theta_idxs=7:6+length(s.point_system.groups), 
-        moment_frac=s.point_system.groups[1].moment_frac)
+    x_ = zeros(3+3+length(s.point_system.groups))
+    jac_ = zeros(length(x_), length(y_))
 
     @parameters begin
         last_y[eachindex(y_)] = y_