function model = get_particle_model(friction)
    import casadi.*
    g = 10;
    a_n = 20;
    mu = 0.2;
    % Symbolic variables and bounds
    q = SX.sym('q',2); 
    v = SX.sym('v',2);
    t = SX.sym('t',1);

    if friction
        % multipliers
        alpha = SX.sym('alpha', 3);
        theta = SX.sym('theta', 3);

        theta_expr = [(alpha(1)+alpha(2))-(alpha(1)*alpha(2));...
            (1-alpha(1))*(1-alpha(2))*alpha(3);...
            (1-alpha(1))*(1-alpha(2))*(1-alpha(3))];

        g_z = theta - theta_expr;

        % Switching function
        c = [q(2); v(2); v(1)];

        f_ode = [v;0;-g;1];
        f_aux_n = [0;0;0;a_n;0];
        f_aux_t = [0;0;a_n*mu;0;0];
        
        f_1 = f_ode;
        f_2 = f_aux_n - f_aux_t;
        f_3 = f_aux_n + f_aux_t;
        % dynamics
        f_x = theta(1)*f_1 + theta(2)*f_2 + theta(3)*f_3;
    else
        % multipliers
        alpha = SX.sym('alpha', 2);
        theta = SX.sym('theta', 2);

        theta_expr = [(alpha(1)+alpha(2))-(alpha(1)*alpha(2));...
            (1-alpha(1))*(1-alpha(2))];

        g_z = theta - theta_expr;

        % Switching function
        c = [q(2); v(2)];

        f_ode = [v;0;-g;1];
        f_aux_n = [0;0;0;a_n;0];
        
        f_1 = f_ode;
        f_2 = f_aux_n;
        % dynamics
        f_x = theta(1)*f_1 + theta(2)*f_2;
    end
    model = NosnocModel();
    model.x = [q;v;t];
    model.alpha = alpha;
    model.z = theta;
    model.g_z = g_z;
    model.c = c;
    model.f_x = f_x;

    model.x0 = [0;1;3;0;0];
    model.custom_time_freezing = 1;
    model.general_inclusion = 1;
end
