T = 4*radius    # VIV period

# initial condition FSI
p0=radius/3; v0=0.1*L; a0=0; t0=0

# motion function uses global var to adjust
posx(t) = p0 + (t-t0)*v0

function sdf_1(x,t)
    √sum(abs2,(2*L >= x[1] .> L) ? ((1.0*x[2] - L)) : x ) - thk/2
end

function sdf_2(x,t)
    √sum(abs2,(2*L <= x[1] .< 3*L) ? (x[1] .- (2.0*x[2])) : x ) - thk/2
end

sdf_list = [sdf_1,sdf_2]

function map(x,t)
    α = amp*cos(t*U/L); R = SA[cos(α) sin(α); -sin(α) cos(α)]
    R * (x - SA[3L-L*sin(0*t*U/L),4L])
end

body = AutoBody(map) do x,t  # signed distance function
        minimum(sdf_list) do SDF
            SDF(x,t)
        end
    end

# generate sim
sim = Simulation((10L,10L), (U,0), radius; ν=U*radius/Re, body=body)

# get start time
duration=10; step=0.1; t₀=round(sim_time(sim))

@time @gif for tᵢ in range(t₀,t₀+duration;step)

    # update until time tᵢ in the background
    t = sum(sim.flow.Δt[1:end-1])
    while t < tᵢ*sim.L/sim.U

        # measure body
        measure!(sim,t)

        # update flow
        mom_step!(sim.flow,sim.pois)
        
        Δt = sim.flow.Δt[end]
        p0 += Δt*v0 
        t0 = t; t += Δt
    end

    # plot vorticity
    @inside sim.flow.σ[I] = WaterLily.curl(3,I,sim.flow.u)*sim.L/sim.U
    flood(sim.flow.σ; shift=(-0.5,-0.5),clims=(-5,5))
    body_plot!(sim); plot!(title="tU/L $tᵢ")
    
    # print time step
    println("tU/L=",round(tᵢ,digits=4),", Δt=",round(sim.flow.Δt[end],digits=3))
end

using WaterLily,StaticArrays
function two_circles(n,m;Re=250,U=1)
    body1 = AutoBody((x,t)->√(x'*x) - m/8, map=(x,t)->x .- m/3 - SA_F32[t,0])
    body2 = AutoBody((x,t)->√(x'*x) - m/6, map=(x,t)->x .- 2m/3 + SA_F32[t,0])
    Simulation((n,m), (U,0), radius; ν=U*radius/Re, body=body1+body2)
end