function model = get_disc_model(r1, r2, q_target1, q_target2)
    import casadi.*
    m1 = 2;
    m2 = 1;
    cv = 2;
    eps = 1e-1;
    M = diag([m1;m1;m2;m2]); % inertia/mass matrix;
    a_n = 10;

    %% SOLUTION STARTS HEREE 
    %% Symbolic variables and bounds
    q = SX.sym('q',4);
    v = SX.sym('v',4);
    t = SX.sym('t',1);
    u = SX.sym('u',2);
    alpha = SX.sym('alpha', 2);
    theta = SX.sym('theta', 2);

    q1 = q(1:2);
    q2 = q(3:4);
    v1 = v(1:2);
    v2 = v(3:4);
    x = [q;v;t];

    %% Hand crafted time freezing
    theta_expr = [(alpha(1)+alpha(2))-(alpha(1)*alpha(2)); (1-alpha(1))*(1-alpha(2))];
    g_z = theta - theta_expr;

    f_c = norm(q1-q2)^2-(r1+r2)^2;
    J_normal = f_c.jacobian(q)';
    v_normal = J_normal'*v;

    c = [f_c;v_normal];

    f_drag = cv*[v1/norm(v1+eps);v2/norm(v2+eps)];
    f_v = [u;...
        zeros(2,1)]-f_drag;

    f_ode = [v;f_v;1];
    f_aux_n = [0;0;0;0;inv(M)*J_normal*a_n;0];

    f_1 = f_ode;
    f_2 = f_aux_n;
    
    f_x = theta(1)*f_1 + theta(2)*f_2;

    %% Initial values and LSQ
    q10 = [-1; 0];
    q20 = [1;0];
    v10 = [0;0];
    v20 = [0;0];

    x0 = [q10;q20;v10;v20;0];
    ubx = [10; 10;10; 10; 5; 5; 5; 5;inf]; 
    lbx = -ubx;
    ubu = [20;20];
    lbu = -ubu;

    x_ref = [q_target1;q_target2;zeros(4,1)];
    u_ref = [0;0];

    Q = diag([10;10;20;20;0*ones(4,1)]);
    R = diag([0.1 0.1]);
    Q_terminal = 100*diag([10;10;20;20;10;10;10;10]);

    %% Add components to modek
    model = NosnocModel();
    model.x = x;
    model.u = u;
    model.alpha = alpha;
    model.z = theta;
    model.g_z = g_z;
    model.x0 = x0;
    model.dims.n_dim_contact = 2;

    model.f_x = f_x;
    model.c = c;

    % box constraints on controls and states
    model.lbu = lbu;
    model.ubu = ubu;
    model.lbx = lbx;
    model.ubx = ubx;
    %% Objective
    model.f_q = ([q;v]-x_ref)'*Q*([q;v]-x_ref)+ u'*R*u;
    model.f_q_T = ([q;v]-x_ref)'*Q_terminal*([q;v]-x_ref);

    %% Extra options for custom
    model.custom_time_freezing = 1;
    model.general_inclusion = 1;
end
