clear
close all
clc


syms amin amax Vmax
syms dav dvt dxt t
syms aav avt axt
syms vel


% Równanie hamowania
decel_firstEquation = struct('x', 0, 'y', 0)
decel_secondEquation = struct('x', Vmax/2, 'y', amax)
decel_thirdEquation = struct('x', Vmax, 'y', amin)
[dav, dvt, dxt] = parabolicAccelerationTrajectory(decel_firstEquation, decel_secondEquation, decel_thirdEquation);

syms decel_v(t) decel_a(vel) decel_x(t)
decel_v(t) = dvt
decel_a(vel) = dav
decel_x(t) = dxt


% Równanie przyspieszenia
accel_firstEquation = struct('x', 0, 'y', amin)
accel_secondEquation = struct('x', Vmax/2, 'y', amax)
accel_thirdEquation = struct('x', Vmax, 'y', 0)
[aav, avt, axt] = parabolicAccelerationTrajectory(accel_firstEquation, accel_secondEquation, accel_thirdEquation);


syms accel_v(t) accel_a(vel) accel_x(t)
accel_v(t) = avt
accel_a(vel) = aav
accel_x(t) = axt






function [av, vt, xt] = parabolicAccelerationTrajectory(firstPointEquation, secondPointEquation, thirdPointEquation)
    % Wyprowadzenie funkcji zależności przyspieszenia od prędkości
    % Deklaracja zmiennych symbolicznych
    syms amin amax Vmax
    syms a2 a1 a0

    syms a(v)
    syms v(t)
    syms x(t)
    syms v

    % Zależność w formie funkcji kwadratowej
    a(v) = sym(a2*v^2 + a1*v + a0);

    % Określamy trzy punkty przez które ma przechodzić parabola:

    eq1 = a(firstPointEquation.x) == firstPointEquation.y;
    eq2 = a(secondPointEquation.x) == secondPointEquation.y;
    eq3 = a(thirdPointEquation.x) == thirdPointEquation.y;

    eqs = [eq1 eq2 eq3];

    % Znalezienie współczynników a0 a1 a2
    vars = [a2 a1 a0];
    sol = solve(eqs, vars);

    a0 = sol.a0;
    a1 = sol.a1;
    a2 = sol.a2;

    % Rozwiązanie
    a(v) = a2*v^2 + a1*v + a0;


    % Wyznaczanie wzoru funkcji v(t), gdzie v(t) = integ(a(v))
    syms v(t) v0

    % Równanie różniczkowe powstałe z całkowania równania v(t) = integ(a(v))
    ode = diff(v, t) == a(v);
    cond = v(0) == v0;

    v(t) = dsolve(ode,cond);


    % Wyznaczanie wzoru funkcji x(t), gdzie x(t) = integ(v(t))
    syms x(t) x0
    ode = diff(x, t) == v(t);
    cond = x(0) == x0;

    x(t) = dsolve(ode, cond);

    syms vel

    av = subs(a2*vel^2 + a1*vel + a0);
    vt = v(t);
    xt = x(t);

    parabolicAccelerationTrajectory = [av, vt, xt];
end