# How can I simplify this code

How can I simplify the creation of a system of equations in a loop:

``````using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra

G_EARTH     = Float64[0.0, 0.0, -9.81]          # gravitational acceleration     [m/s²]
L0::Float64 = 5.0                               # initial segment length            [m]
V0::Float64 = 2                                 # initial velocity of lowest mass [m/s]
M0::Float64 = 0.5                               # mass per particle                [kg]
C_SPRING::Float64 = 50                          # spring constant
segments::Int64 = 5                             # number of tether segments         [-]
α0 = π/10                                       # initial tether angle            [rad]
duration = 30.0                                 # duration of the simulation        [s]
POS0 = zeros(3, segments+1)
VEL0 = zeros(3, segments+1)
ACC0 = zeros(3, segments+1)
SEGMENTS0 = zeros(3, segments)
UNIT_VECTORS0 = zeros(3, segments)
for i in 1:segments+1
local l0
l0 = -(i-1)*L0
v0 = (i-1)*V0/segments
POS0[:, i] .= [sin(α0) * l0, 0, cos(α0) * l0]
VEL0[:, i] .= [sin(α0) * v0, 0, cos(α0) * v0]
end
for i in 2:segments+1
ACC0[:, i] .= G_EARTH
end
for i in 1:segments
UNIT_VECTORS0[:, i] .= [0, 0, 1.0]
SEGMENTS0[:, i] .= POS0[:, i+1] - POS0[:, i]
end

# defining the model, Z component upwards
@parameters mass=M0 c_spring0=C_SPRING damping=0.5 l_seg=L0
@variables t
@variables pos(t)[1:3, 1:segments+1]  = POS0
@variables vel(t)[1:3, 1:segments+1]  = VEL0
@variables acc(t)[1:3, 1:segments+1]  = ACC0
@variables segment(t)[1:3, 1:segments]  = SEGMENTS0
@variables unit_vector(t)[1:3, 1:segments]  = UNIT_VECTORS0
@variables norm1(t)[1:segments] = l_seg * ones(segments)
@variables rel_vel(t)[1:3, 1:segments]  = zeros(3, segments)
@variables spring_vel(t)[1:segments] = zeros(segments)
@variables c_spring(t)[1:segments] = c_spring0 * ones(segments)
@variables spring_force(t)[1:3, 1:segments] = zeros(3, segments)
@variables total_force(t)[1:3, 1:segments] = zeros(3, segments)
D = Differential(t)

eqs1 = vcat(D.(pos) ~ vel,
D.(vel) ~ acc)
eqs2 = []
for i in segments:-1:1
global eqs2
eqs2 = vcat(eqs2, segment[:, i] ~ pos[:, i+1] - pos[:, i])
eqs2 = vcat(eqs2, norm1[i] ~ norm(segment[:, i]))
eqs2 = vcat(eqs2, unit_vector[:, i] ~ -segment[:, i]/norm1[i])
eqs2 = vcat(eqs2, rel_vel[:, i] ~ vel[:, i+1] - vel[:, i])
eqs2 = vcat(eqs2, spring_vel[i] ~ -unit_vector[:, i] ⋅ rel_vel[:, i])
eqs2 = vcat(eqs2, c_spring[i] ~ c_spring0 * (norm1[i] > l_seg))
eqs2 = vcat(eqs2, spring_force[:, i] ~ (c_spring[i] * (norm1[i] - l_seg) + damping * spring_vel[i]) * unit_vector[:, i])
if i == segments
eqs2 = vcat(eqs2, total_force[:, i] ~ spring_force[:, i])
else
eqs2 = vcat(eqs2, total_force[:, i] ~ spring_force[:, i]- spring_force[:, i+1])
end
eqs2 = vcat(eqs2, acc[:, i+1] .~ G_EARTH + total_force[:, i] / mass)
end
eqs2 = vcat(eqs2, acc[:, 1] .~ zeros(3))
eqs = vcat(eqs1..., eqs2)

@named sys = ODESystem(eqs, t)
simple_sys = structural_simplify(sys)
``````

For example:
Can I use only one vcat call for these two lines?
Or can I use some other function than vcat?

``````    eqs2 = vcat(eqs2, unit_vector[:, i] ~ -segment[:, i]/norm1[i])
eqs2 = vcat(eqs2, rel_vel[:, i] ~ vel[:, i+1] - vel[:, i])
``````

I tried `push!` and `append!` with no success…

I found out, some of the equations in the array eqs2 have the type `Equation` and some the type `BasicSymbol`. Perhaps that causes the issues?

Why is this not of type Equation:

``````julia> eqs2[1]
``````julia> eqs2[end]