# Modeling a Pendulum in Cartesion Coordinates

model Pendulum
    constant Real PI = 3.141596;
    parameter Real m = 1;
    parameter Real L = 0.5;
    constant Real g= 9.81;
    Real F;
    output Real x(start=0.5), y(start=0);
    output Real vx, vy;
initial equation
    vx = 0;
    vy = 0;
equation
  m*der(vx) = -(x/L)*F;
  m*der(vy) = -(y/L)*F - m *g;
  der(x) = vx;
  der(y) = vy;
  x^2 + y^2 = L^2;
end Pendulum;

using DifferentialEquations
using ModelingToolkit
using Plots

@variables t x(t) y(t) vx(t) vy(t) F(t)
D = Differential(t)

@parameters m L
@constants g = 9.81

eqs = [m*D(vx) ~ -(x/L)*F,
       m*D(vy) ~ -(y/L)*F - m*g,
       D(x) ~ vx,
       D(y) ~ vy,
       x^2 + y^2 ~ L^2]

@named sys = ODESystem(eqs, t)

sys = structural_simplify(sys)

u0 = [vx => 0, vy => 0, x => 0.5, y => 0, F=>1, D(x)=>1.0]
p = [m => 1.0, L => 0.5]

tspan = (0.0, 0.4)
prob = ODEProblem(sys, u0, tspan, p, jac=true)
@time sol = solve(prob)

display(plot(sol, idx = [sys.x, sys.y]))