Simulation results

src/FlightMechanicsSimulator.jl

@@ -10,7 +10,9 @@ using FlightMechanicsUtils
 import Base: eltype
 
 
-export ConstantInput, StepInput, DoubletInput, RampInput, SinusoidalInput, get_value
+export Input
+export get_value
+export ConstantInput, StepInput, DoubletInput, RampInput, SinusoidalInput
 include("models/inputs.jl")
 
 export ISA1976, F16StevensAtmosphere

src/models/atmosphere.jl

@@ -86,6 +86,11 @@ Calculate atmosphere for given height, `h` (m).
 # Returns
 - `F16StevensAtmosphere`
 """
+function F16StevensAtmosphere{T}(h) where T<:Number
+    T_, ρ, a, p = F16Stevens.atmosphere_f16(h)
+    return F16StevensAtmosphere{T}(T_, p, ρ, a, h)
+end
+
 function F16StevensAtmosphere(h)
     T, ρ, a, p = F16Stevens.atmosphere_f16(h)
     return F16StevensAtmosphere(T, p, ρ, a, h)
@@ -128,4 +133,5 @@ Calculate atmosphere for given height, `h` (m).
 # Returns
 - `ISA1976`
 """
+ISA1976{T}(h) where T<:Number = ISA1976{T}(atmosphere_isa(h)..., h)
 ISA1976(h) = ISA1976(atmosphere_isa(h)..., h)

src/simulation.jl

@@ -1,45 +1,84 @@
+mutable struct SimDEData{
+        T,
+        N,
+        DSSD<:DSStateDot,
+        C<:Input,
+        AC<:Aircraft,
+        ATM<:Atmosphere,
+        G<:Gravity
+    } <: DEDataVector{T}
+
+    x::SVector{N, T}
+    t::T
+    dssd::DSSD
+    controls::Array{C, 1}
+    aircraft::AC
+    atmosphere::ATM
+    gravity::G
+end
+
+
 """
     simulate(
         tini, tfin, dss, controls, aircraft, atmosphere, gravity;
         solver=TSit5(), solve_args=Dict()
     )
 
 """
-function simulate(tini, tfin, dss, controls, aircraft, atmosphere, gravity;
+function simulate(tini, tfin, dssd, controls, aircraft, atmosphere, gravity;
     solver=Tsit5(), solve_args=Dict()
     )
 
     tspan = (tini, tfin)
-    p = [dss, controls, aircraft, atmosphere, gravity]
 
-    x0 = get_x(dss)
-    prob = ODEProblem{false}(f, x0, tspan, p)
+    x0 = SimDEData(
+        get_x(dssd),
+        tini,
+        dssd,
+        controls,
+        aircraft,
+        atmosphere(get_height(dssd)),
+        gravity
+    )
+
+    prob = ODEProblem{false}(f, x0, tspan)
     sol = solve(prob, solver; solve_args...)
 
     df = DataFrame(sol')
-    rename!(df, get_x_names(dss))
+    rename!(df, get_x_names(dssd))
     df[!, :time] = sol.t
 
-    return df
+    return df, sol
 end
 
 
 function f(x, p, t)
-    dss = p[1]
-    controls = p[2]
-    aircraft = p[3]
-    atmosphere = p[4]
-    gravity = p[5]
+    dssd = x.dssd
+    controls = x.controls
+    aircraft = x.aircraft
+    atmosphere = x.atmosphere
+    gravity = x.gravity
 
-    controls_arr = get_value.(controls, t)
+    controls_arr = get_value.(x.controls, t)
 
+    dss = get_ds_state(dssd)
     dss = typeof(dss)(x)
 
-    atmosphere = atmosphere(get_height(dss))
+    atmosphere = typeof(x.atmosphere)(get_height(x.dssd))
+
+    dssd, outputs = f(time, dss, controls_arr, x.aircraft, x.atmosphere, x.gravity)
 
-    dssd, outputs = f(time, dss, controls_arr, aircraft, atmosphere, gravity)
+    x_dot = SimDEData(
+        get_xdot(dssd),
+        t,
+        dssd,
+        x.controls,
+        x.aircraft,
+        x.atmosphere,
+        x.gravity
+    )
 
-    return get_xdot(dssd)
+    return x_dot
 end
 
 

test/dynamics/sixdof_body_euler.jl

@@ -156,18 +156,26 @@ end
 
     kwargs_sim = Dict(:solver=>RK4(), :solve_args=>Dict(:reltol=>1e-10, :saveat=>dt))
 
-    r_base = simulate(
+    r_base, sol_base = simulate(
         t0,
         t1,
-        get_ds_state(dssd_base_trim),
+        dssd_base_trim,
         args_sim...;
         kwargs_sim...
     )
 
-    r = simulate(
+    r, sol = simulate(
         t0,
         t1,
-        SixDOFBodyEuler(get_ds_state(dssd_base_trim)),
+        DSStateDot(
+            SixDOFBodyEuler(get_ds_state(dssd_base_trim)),
+            [get_accel_body(dssd_base_trim)...,
+             get_euler_angles_rates(dssd_base_trim)[3:-1:1]...,
+             get_ang_accel_body(dssd_base_trim)...,
+             get_horizon_velocity(dssd_base_trim)...,
+             0.0
+             ]
+        ),
         args_sim...;
         kwargs_sim...
     )

test/f16/simulation.jl

@@ -77,7 +77,7 @@ cost =
 @test isapprox(x_dot[12], 0.0, atol = 1e-4)
 
 # RETRIM to refine flying condition
-dssd, controls_trim, outputs_trim, cost = trim(
+dssd_trim, controls_trim, outputs_trim, cost = trim(
     SixDOFAeroEuler(x_stev),
     controls_stev,
     F16(F16Stevens.MASS, F16Stevens.INERTIA, xcg),
@@ -87,8 +87,8 @@ dssd, controls_trim, outputs_trim, cost = trim(
     0.3,
 )
 
-x_trim = get_x(dssd)
-x_dot_trim = get_xdot(dssd)
+x_trim = get_x(dssd_trim)
+x_dot_trim = get_xdot(dssd_trim)
 
 # Linear acceleration
 @test isapprox(x_dot_trim[1:3], zeros(3), atol=1e-12)
@@ -107,16 +107,16 @@ x = x_trim
 # Transform to Input for simulate
 controls = ConstantInput.(controls_trim)
 
-results = simulate(
+results, sol = simulate(
     t0,
     t1,
-    SixDOFAeroEuler(x),
+    dssd_trim,
     controls,
     F16(F16Stevens.MASS, F16Stevens.INERTIA, xcg),
     F16StevensAtmosphere,
     LHDownGravity(FlightMechanicsSimulator.F16Stevens.GD*FT2M);
     solver=RK4(),
-    solve_args=Dict(:reltol=>1e-10, :saveat=>dt),
+    solve_args=Dict(:reltol=>1e-10, :saveat=>collect(t0:dt:t1)),
     )
 
 # Check X, Y against Stevens