You can't do what you are hoping because you can't cover a sphere with squares. In fact, you can't even put one square on the sphere. If you look closely, the angles of your squares add up to more than $360^\circ$, by an amount that is determined by the solid angle taken up by the square.
A few choices that maintain some of what you want. 1) The black grid in your first figure comes from the classic lines of latitude and longitude. You get quadrilaterals except near the poles. You can choose what latitude to make the cells almost square, but they will be much larger near the equator than near the poles. 2) Use the latitude/longitude lines as in 1, but have fewer cells as you get near the poles. You can keep the cells close to square, but won't have a nice grid. 3) Use the square grid as you are, but recognize that you can't have one grid that covers the sphere. Maybe a grid that covers about half the sphere centered on a point of current interest meets your needs.
I'd ignore the radius $R$ (i.e. take the earth's radius as unit of length) and then convert to Cartesian coordinates:
\begin{align*}
x&=\cos\phi\cos\lambda\\
y&=\cos\phi\sin\lambda\\
z&=\sin\phi
\end{align*}
These vectors point from the center of the sphere to the vertices of your polygon. Two such vertices are joined by a greatcircle arc. That greatcircle is the intersection of your sphere with a plane through the origin. That plane through the origin is defined by its normal vector. The cross product of two vectors is perpendicular to both of them, so you can use the cross product to compute these normal vectors. Writing $v_i=(x_i,y_iz_i)$ for the vertices, you get the normals of the edges as
$$n_i = v_i\times v_{i+1}$$
The angle at a vertex is equal to the angle between the corresponding normal vectors. You can compute that using the dot product, dividing by the lengths to ensure normalization:
$$\alpha_i = \pm\arccos-
\frac{n_{i-1}\cdot n_{i}}{\lVert n_{i-1}\rVert\,\lVert n_i\rVert}$$
The minus sign in there is because for a normalized dot product of $1$, both normals point in the same direction, so you have zero change in direction but the inner angle there is $180$.
So which sign should you choose for the $\pm$ in there? For that, look at the determinant of the $3\times3$ matrix formed by $v_{i-1}$, $v_i$ and $v_{i+1}$. The sign of that matrix will tell you whether the triangle these three vectors form on the surface is oriented clockwise or counter-clockwise. Take the sign from this and the value from the $\arccos$ above and you should be almost done.
As usual with signs, you have a decision to make: which sign is which? Well, one choice of sign (e.g. exactly as described above) will lead you one polygon, the other (with the sign of the determinant flipped, or eqivalently the order of vectors inside the determinant reversed) leads to the complementary polygon. One will have all its vertices in clockwise order, the other in counter-clockwise order. Trying this out on a tiny example (e.g. three points forming three right angles) will tell you which one is which.
Best Answer
The planes carve a regular spherical polygon in the surface of the (unit) sphere, such that the dihedral angles $\delta$ between the planes match the exterior angles of the polygon. The semi-vertical angle $\alpha$ of the cone gives the (spherical) circumradius of the polygon.
If $P$ is the center of the polygon, $Q$ one of its vertices, and $A$ the midpoint of a side adjacent to $Q$, then $\triangle PQA$ is a spherical right triangle with hypotenuse $\alpha$ and acute angles $\pi/n$ and $(\pi-\delta)/2$.
By the Spherical Law of Cosines for Angles, we have
$$\begin{align} \cos A &= -\cos P \cos Q + \sin P \sin Q \cos\alpha \\[4pt] 0 &= -\cos\frac{\pi}{n}\sin\frac\delta2+\sin\frac{\pi}{n}\cos\frac\delta2\cos\alpha \\[4pt] \tan\frac{\delta}{2} &= \tan\frac{\pi}{n}\;\cos\alpha \end{align}$$ as desired. $\square$