function [val, grads] = get_R(w)
    % input: w \in C^3
    % output: R \in C^(3x3)
    % grads = {3x3, 3x3, 3x3}

    beta = sqrt(w.'*w);
    w_hat = w/beta;

    f0 = get_fn(beta, 0);
    f2 = get_fn(beta, 2);

    val = 1/8*((f0+f2)*eye(3) + (f0-3*f2)*(w_hat*w_hat.'));     % Compute R(w) as a 3x3 matrix. 

    f1 = get_fn(beta, 1);
    f3 = get_fn(beta, 3);
    
    I = eye(3);
    
    grads = cell(3, 1);

    for idx = 1:3
        % Compute the derivative of R(w) w.r.t. w=(w1, w2, w3). 
        d_w_hat = (I(:,idx) - w_hat(idx)*w_hat)/beta;
        grads{idx}  = 1/8*(1i*(f1+f3)*w_hat(idx)*eye(3) + 1i*(f1-3*f3)*w_hat(idx)*(w_hat*w_hat.') ...
                + (f0-3*f2)*(d_w_hat*w_hat.' + w_hat*d_w_hat.'));   
    end
end

