function [bc_params, d1, d2] = opt_bezier(P, S, vp, vs, aR_max, err_max)

    vpN = norm(vp); vsN = norm(vs);
    theta1 = atan2(vp(2),vp(1)); theta2 = atan2(vs(2), vs(2));
    
    Rp = vpN^2 / aR_max;
    Rs = vsN^2 / aR_max;
    
    g(1) = S(1) - P(1);
    h(1) = S(2) - P(2);
    g(2) = 2*sin(theta2-theta1);
    h(2) = 2*(h(1)*cos(theta1)-g(1)*sin(theta1));
    g(3) = 2*sin(theta2-theta1);
    h(3) = 2*(h(1)*cos(theta2)-g(1)*sin(theta2));
    
    d1 = 1; d2 = 1;
    while 1
        rho_p = 3*d1^2 / (h(2) + d2*g(2));
        rho_s = 3*d2^2 / (h(3) + d1*g(3));
        
        err1 = Rp - rho_p;
        err2 = Rs - rho_s;
        error = max(abs(err1), abs(err2));
        
        fprintf('error = %.4e | d1 = %.4e | d2 = %.4e\n', error, d1, d2);
        pause;
        if error < err_max
            break;
        end
        
        d1 = d1 + err1 / Rp;
        d2 = d2 + err2 / Rs;
    end
    
    bc_params = [P (P+d1*[cos(theta1);sin(theta1)]) (S-d2*[cos(theta2);sin(theta2)]) S];
end
