How to calculate the area covered by any spherical rectangle?

Is there any analytic or generalized formula to calculate area covered by any rectangle having length $l$ & width $b$ each as a great circle arc on a spherical surface with a radius $R$? i.e. How to find the area $A$ of rectangle in terms of length $l$, width $b$ and radius $R$ ($A=f(l, b,R)$)?

Note: Spherical rectangle is a quadrilateral having equal opposite sides but non-parallel & all the interior angles are equal in magnitude & each one is greater than $90^\circ$.

Assume we are working on a sphere of radius $1$, or consider the lengths in radians and the areas in steradians.

Extend the sides of length $l$ until they meet. This results in a triangle with sides $$ w,\quad\frac\pi2-\frac l2,\quad\frac\pi2-\frac l2 $$

The Spherical Law of Cosines says that $$ \begin{align} \cos(A) &=\frac{\cos\left(\frac\pi2-\frac l2\right)-\cos\left(\frac\pi2-\frac l2\right)\cos(w)}{\sin\left(\frac\pi2-\frac l2\right)\sin(w)}\\ &=\frac{\sin\left(\frac l2\right)}{\cos\left(\frac l2\right)}\frac{1-\cos(w)}{\sin(w)}\\[6pt] &=\tan\left(\frac l2\right)\tan\left(\vphantom{\frac l2}\frac w2\right) \end{align} $$ One quarter of the spherical excess of the rectangle is $D-\frac\pi2$ and $$ \sin\left(D-\frac\pi2\right)=\tan\left(\frac l2\right)\tan\left(\vphantom{\frac l2}\frac w2\right) $$ Therefore, the area of the rectangle is $$ \bbox[5px,border:2px solid #C0A000]{4\sin^{-1}\left(\tan\left(\frac l2\right)\tan\left(\vphantom{\frac l2}\frac w2\right)\right)} $$ Note that for small $l$ and $w$, this is approximately $lw$; and if $l+w=\pi$ (that is, the rectangle is a great circle), we get an area of $2\pi$ (one half the sphere).

Scaling for a sphere of radius $R$ gives $$ \bbox[5px,border:2px solid #C0A000]{4R^2\sin^{-1}\left(\tan\left(\frac l{2R}\right)\tan\left(\vphantom{\frac l2}\frac w{2R}\right)\right)} $$

Note the similarity to the formula for the area of a spherical right triangle with legs $a$ and $b$: $$ 2\tan^{-1}\left(\tan\left(\vphantom{\frac b2}\frac a2\right)\tan\left(\frac b2\right)\right) $$ or for a sphere of radius $R$, $$ 2R^2\tan^{-1}\left(\tan\left(\vphantom{\frac b2}\frac a{2R}\right)\tan\left(\frac b{2R}\right)\right) $$

On a sphere of radius $R > 0$, a geodesic triangle with interior angles $\theta_{1}$, $\theta_{2}$, and $\theta_{3}$ has area $R^{2}(\theta_{1} + \theta_{2} + \theta_{3} - \pi)$. One way you might proceed, therefore, it to triangulate your geodesic polygon (whatever its actual shape) and sum the areas of its triangular pieces. For a quadrilateral with interior angles $\theta_{1}$, $\theta_{2}$, $\theta_{3}$, and $\theta_{4}$, the area is $$ R^{2}(\theta_{1} + \theta_{2} + \theta_{3} + \theta_{4} - 2\pi). $$

Added (in light of OP's clarifications): If $\theta \leq \pi$ is the interior angle at each vertex of the quadrilateral, then $$ \theta = \pi - \cos^{-1}\left\lvert\tan\frac{\ell}{2R} \tan\frac{b}{2R}\right\rvert, $$ so $$ \text{Area} = R^{2}(4\theta - 2\pi) = R^{2} \left[2\pi - 4\cos^{-1}\left\lvert\tan\frac{\ell}{2R} \tan\frac{b}{2R}\right\rvert\right]. $$ (Particularly, the absolute value of the product of tangents inside the arccos does not exceed unity.)

To see this, it's convenient to work in Cartesian coordinates with the sphere centered at the origin. Denote the vertices of the quadrilateral by $$ v_{1} = (A, B, C),\qquad v_{2} = (A, -B, C),\qquad v_{3} = (A, -B, -C),\qquad v_{4} = (A, B, -C). $$ Of course, $A^{2} + B^{2} + C^{2} = R^{2}$.

The plane through the origin, $v_{1}$, and $v_{2}$ has equation $Cx - Az = 0$, and so has unit normal vector $$ n_{1} = \frac{(C, 0, -A)}{\sqrt{A^{2} + C^{2}}} = \frac{(C, 0, -A)}{\sqrt{R^{2} - B^{2}}}. $$ The plane through the origin, $v_{1}$, and $v_{4}$ has equation $Bx - Ay = 0$, and so has unit normal vector $$ n_{2} = \frac{(B, -A, 0)}{\sqrt{A^{2} + B^{2}}} = \frac{(B, -A, 0)}{\sqrt{R^{2} - C^{2}}}. $$

The "large" angle between these planes is the "large" angle between the great circles they determine (because each normal vector $n_{i}$ is tangent to the sphere at $v_{1}$), i.e., the interior angle $\theta$ of the quadrilateral. Taking the ordinary dot product of the unit normals, $$ \cos\theta = n_{1} \cdot n_{2} = \frac{B}{\sqrt{R^{2} - B^{2}}}\, \frac{C}{\sqrt{R^{2} - C^{2}}}. \tag{1} $$

Let $2\psi_{1}$ denote the angle subtended at the center of the sphere by the side from $v_{1}$ to $v_{2}$, and let $\ell_{1} = 2R\psi_{1}$ denote the corresponding side length, so that $\psi_{1} = \ell_{1}/(2R)$. Thinking of $v_{1}$ and $v_{2}$ as vectors in space, the ordinary dot product gives $$ \frac{v_{1} \cdot v_{2}}{R^{2}} = \cos(2\psi_{1}) = 1 - 2\sin^{2} \psi_{1}. $$ On the other hand, using the Cartesian components of these vectors, we have $$ \frac{v_{1} \cdot v_{2}}{R^{2}} = \frac{A^{2} - B^{2} + C^{2}}{A^{2} + B^{2} + C^{2}} = 1 - 2\frac{B^{2}}{R^{2}}. $$ Equating, it follows at once that $B/R = \pm\sin\psi_{1}$. (As a consistency check, the isoceles triangle with the origin, $v_{1}$, and $v_{2}$ as vertices has apex angle $2\psi_{1}$ and base $2B$.)

Similarly, letting $2\psi_{2}$ denote the angle subtended at the center of the sphere by the side from $v_{1}$ to $v_{4}$, and letting $\ell_{2} = 2R\psi_{2}$ denote the corresponding side length, we find $C/R = \pm\sin\psi_{2}$.

That is, $$ B = \pm R\sin\frac{\ell_{1}}{2R},\qquad C = \pm R\sin\frac{\ell_{2}}{2R}, \tag{2} $$ and consequently $$ \frac{B}{\sqrt{R^{2} - B^{2}}} = \pm\tan\frac{\ell_{1}}{2R},\qquad \frac{C}{\sqrt{R^{2} - C^{2}}} = \pm\tan\frac{\ell_{1}}{2R}. \tag{3} $$

Substituting (3) into (1), $$ \cos\theta = \frac{B}{\sqrt{R^{2} - B^{2}}}\, \frac{C}{\sqrt{R^{2} - C^{2}}} = \pm \tan\frac{\ell_{1}}{2R} \tan\frac{\ell_{2}}{2R}. $$ To get the large angle, either take the arccos of the negative value, or $\pi$ minus the arccos of the positive value. (The angle/area formulas above do the latter.)

Let the north-pole $N=(0,0,1)$ be the center of the rectangle, and the planes $x_1=0$, $\>x_2=0$ be its symmetry planes. Let the point $$P=(\cos\phi\sin\theta,\sin\phi\sin\theta,\cos\theta)$$ with $0<\phi<{\pi\over2}$, $\>0<\theta<\pi$ ($\theta=0$ at the north pole!) be one of its vertices. Then we see two rightangled spherical triangles, both with hypotenuse $NP$ of spherical length $\theta$, one of them with angle $\phi$ at $N$ and the other with angle ${\pi\over2}-\phi$ at $N$. The legs opposite these angles have spherical lengths ${a\over2}$ and ${b\over2}$ respectively, where $a$ and $b$ were given in advance.

Napier's rule for rightangled spherical triangles then give $$\sin{a\over2}=\sin\phi\sin\theta,\quad \sin{b\over2}=\cos\phi\sin\theta\ .$$ From these two equations one can compute $\theta$ and $\phi$, and other Napier's formulas then give the missing angles at $P$. Knowing all angles one obtains the area at once; see user86418's answer.