function mesh = create_spherical_cap(R, phi, nR)

% Create a circle first
mesh = create_circle(R*sin(phi), nR);

[th, r] = cart2pol(mesh.Nodes(:,2), mesh.Nodes(:,3));

th0 = atan2(r, R*cos(phi));

% Rescale r for equivalent arc lengths
%r = cos(th0) / cos(phi) .* r; 
r = (cos(phi) ./ cos(th0)) .*r;

mesh.Nodes(:, 2) = r .* cos(th);
mesh.Nodes(:, 3) = r .* sin(th);

%mesh = translate_mesh(mesh, [0 0 R*cos(phi)]);

l = sqrt(r.^2 + (R*cos(phi))^2);
%th0 = acos(R*cos(phi) ./ l);
th0 = atan2(r, R*cos(phi));

z = R*cos(th0);
mesh.Nodes(:,2:3) = bsxfun(@times, mesh.Nodes(:,2:3), R./l);
mesh.Nodes(:,4) = z;

end