Zerocrossing parsing

Author: BenLauwens

src/odes/macros.jl

@@ -14,19 +14,46 @@ macro model(expr::Expr)
   expr.head != :function && error("Expression is not a function definition!")
   args = getArguments(expr)
   func_name = shift!(args)
+  c_vec = Vector{Expr}()
+  a_vec = Vector{Expr}()
   f_vec = Vector{Expr}()
+  zc_vec = Vector{Expr}()
+  neg_vec = Vector{Expr}()
+  pos_vec = Vector{Expr}()
   dx = expr.args[1].args[5]
+  p = 0
   for ex in expr.args[2].args
-    ex isa Expr && ex.head == Symbol("=") && ex.args[1] isa Expr && ex.args[1].head == :ref && ex.args[1].args[1] == dx && push!(f_vec, ex)
+    if ex isa Expr
+      if ex.head == Symbol("=")
+        if ex.args[1] isa Expr && ex.args[1].head == :ref && ex.args[1].args[1] == dx
+          push!(f_vec, ex)
+        elseif ex.args[1] isa Symbol
+          ex.args[2] isa Expr ? push!(a_vec, ex) : push!(c_vec, ex)
+        end
+      elseif ex isa Expr && ex.head == :if && ex.args[1].head == :call && ex.args[1].args[1] == :<
+        push!(zc_vec, ex.args[1].args[2])
+        push!(neg_vec, ex.args[2])
+        length(ex.args) == 3 ? push!(pos_vec, ex.args[3]) : Expr()
+      end
+    end
   end
+  println(zc_vec)
+  println(neg_vec)
+  println(pos_vec)
   n = length(f_vec)
   deps = zeros(Bool, n, n)
+  param_deps = zeros(Bool, n, p)
+  zc_deps = zeros(Bool, p, n)
   for ex in expr.args[2].args
     ex isa Expr && ex.head == Symbol("=") && ex.args[1] isa Expr && ex.args[1].head == :ref && ex.args[1].args[1] == dx && check_dependencies(ex.args[2], deps, ex.args[1].args[2], expr.args[1].args[3])
   end
   esc(:(function $func_name()
       f = Array{Function}($n)
-      $((:(f[$(f_vec[i].args[1].args[2])] = ($(expr.args[1].args[2])::TaylorSeries.Taylor1, $(expr.args[1].args[3])::Vector{TaylorSeries.Taylor1}, $(expr.args[1].args[4])::Vector{Float64})->$(f_vec[i].args[2])) for i in 1:length(:($f_vec)))...)
+      $((:(f[$(f_vec[i].args[1].args[2])] = ($(expr.args[1].args[2])::TaylorSeries.Taylor1, $(expr.args[1].args[3])::Vector{TaylorSeries.Taylor1}, $(expr.args[1].args[4])::Vector{Float64})->begin
+      $((:($(c)) for c in :($c_vec))...)
+      $((:($(a)) for a in :($a_vec))...)
+      $(f_vec[i].args[2])
+      end) for i in 1:length(:($f_vec)))...)
       Model(f, $deps)
     end))
 end

test/continuous.jl

@@ -18,29 +18,21 @@ cont = @continuous diffeq(sim, [0.0, 20.0], [2020.0]; stiff=false, order=4)
 run(sim, 71)
 
 @model function bouncing_ball(t, x, p, dx)
-  dx[1] = x[3]
-  dx[3] = zero(t)
-  dx[2] = x[4]
-  dx[4] = -9.81 + p[1]*(2.0e4/0.1*(x[2]-(3.0-ceil(evaluate(x[1])/0.2)*0.2)))
-end
-
-@zerocrossing function contact(t, x, p, res)
-  yf = 3.0 - ceil(x[1]/0.2)*0.2
-  res = x[4] < 0.0 && x[2] - yf < 0.05
-  if res
+  g = 9.81
+  m = 1.0
+  b = 30.0
+  k = 1.0e6
+  if x[1] < 0.0
     p[1] = 1.0
-  end
-end
-
-@zerocrossing function up(t, x, p, res)
-  yf = 3.0 - ceil(x[1]/0.2)*0.2
-  res = x[4] > 0.0 && x[2] - yf > 0.05
-  if res
+  else
     p[1] = 0.0
   end
+  f = k*x[1] + b*x[2]
+  dx[1] = x[2]
+  dx[2] = -g - p[1]*f/m
 end
 
 sim = Simulation()
-cont = @continuous bouncing_ball(sim, [0.0, 3.05, 1.0, 0.0], [0.0])
+cont = @continuous bouncing_ball(sim, [2.0, 0.0], [0.0])
 @process report(sim, cont)
 run(sim, 5)