How can one establish that the inverse stereographic projection $f : \mathbb R^2 \to \mathbb R^3$ given by $$ f(x,y):= \frac{(2x,2y,x^2+y^2-1)}{x^2+y^2+1} $$ is a conformal (angle-preserving) map?
I wanted to show by a direct calculation that this map preserves angles. I considered two vectors $(x_1,y_1),(x_2,y_2) \in \mathbb R^2$ as well as two vectors $f(x_1,y_1),f(x_2,y_2) \in \mathbb R^3$.
If I denote $\theta$ to be the angle between $(x_1,y_1),(x_2,y_2)$ and $f(\theta)$ to be the angle between $f(x_1,y_1),f(x_2,y_2)$, then I want to establish $f(\theta)=\theta$.
First, I observe that $f(x,y)$ has unit magnitude: \begin{align*} \|f(x,y)\| &= \frac{\sqrt{(2x)^2+(2y)^2+(x^2+y^2-1)^2}}{x^2+y^2+1} \\ &= \frac{\sqrt{4(x^2+y^2)+((x^2+y^2)^2-2(x^2+y^2)+1)}}{x^2+y^2+1} \\ &= \frac{\sqrt{(x^2+y^2)^2+2(x^2+y^2)+1}}{x^2+y^2+1} \\ &= \frac{\sqrt{(x^2+y^2+1)^2}}{x^2+y^2+1} \\ &= 1, \end{align*} confirming that the inverse stereographic projection maps the plane back to the unit sphere except the north and south poles.
Well, now I have \begin{align*} \cos f(\theta) &= \frac{f(x_1,y_1)\cdot f(x_2,y_2)}{\|f(x_1,y_1)\|\|f(x_2,y_2)\|} \\ &= \frac{(2x_1,2y_1,x_1^2+y_1^2-1) \cdot(2x_2,2y_2,x_2^2+y_2^2-1)}{1 \cdot 1} \\ &= \frac{(2x_1)(2x_2)+(2y_1)(2y_2)+(x_1^2+y_1^2-1)(x_2^2+y_2^2-1)}{(x_1^2+y_1^2+1)(x_2^2+y_2^2+1)} \end{align*} and \begin{align*} \cos\theta &= \frac{(x_1,y_1)\cdot (x_2,y_2)}{\|(x_1,y_1)\|\|(x_2,y_2)\|} \\ &= \frac{x_1x_2+y_1y_2}{\sqrt{x_1^2+y_1^2}\sqrt{x_2^2+y_2^2}}. \end{align*} However, I am stuck on how to get my expressions of $\cos \theta$ and of $\cos f(\theta)$ to equal.
Parametrizing curves and taking derivatives will work, but I prefer to work directly with the tangent spaces. In this case, the tangent space of your domain, which is $\mathbb{R}^2$, it actually $\mathbb{R}^2$ itself. If you want to show that your map $f$ is conformal, you need to show that for vectors $u$ and $v$ in the tangent space of $(x,y)\in\mathbb{R}$ there is a number $\lambda(x,y)$ so that $u\cdot v=\lambda(x,y) df_{(x,y)}(u)\cdot df_{(x,y)}(v)$. The notation here is that $df_{(x,y)}$ is the Jacobian of $f$ evaluated at $(x,y)$. Computing the Jacobian is not so hard: $$df_{(x,y)}=\frac{1}{(x^2+y^2+1)^2}\begin{bmatrix} 2(y^2-x^2+1)& -4xy\\ -4xy & 2(x^2-y^2+1) \\ -4x& -4y\end{bmatrix}$$ Disregard the coefficient in front of this matrix for now and label the entries $f_{i,x}$ and $f_{i,y}$. If $u=(u_1,u_2)$ and $v=(v_1,v_2)$, then evaluating $df_{(x,y)}$ on these vectors in $\mathbb{R}^2$ gives two vectors in $\mathbb{R}^3$ which we can write compactly as $$df_{(x,y)}(u)=\frac{1}{(x^2+y^2+1)^2}(f_{i,x}u_1+f_{i,y}u_2)_{i=1}^3\quad \quad df_{(x,y)}(v)=\frac{1}{(x^2+y^2+1)^2}(f_{i,x}v_1+f_{i,y}v_2)_{i=1}^3$$ The dot product of these vectors is therefore the sum of the componentwise products \begin{align}df_{(x,y)}(u)\cdot df_{(x,y)}(v)&=\left[\sum_{i=1}^3 f_{i,x}^2u_1v_1+f_{i,x}f_{i,y}(u_1v_2+u_2v_1)+f_{i,y}^2u_2v_2\right]\frac{1}{(x^2+y^2+1)^4} \end{align} But for clarity we will rearrange as $$\left[u_1v_1\left(\sum_{i=1}^3 f_{i,x}^2\right)+u_2v_2\left(\sum_{i=1}^3f_{i,y}^2\right)+(u_1v_2+u_2v_1)\sum_{i=1}^3 f_{i,x}f_{i,y}\right]\frac{1}{(x^2+y^2+1)^4}$$ It takes some elbow grease, but you can show that $\sum_{i=1}^3 f_{i,x}^2=\sum_{i=1}^3 f_{i,y}^2=4(x^2+y^2+1)^2$ and also that $\sum_{i=1}^3 f_{i,x}f_{i,y}=0$. I will leave the former to you, for the latter: \begin{align} \sum_{i=1}^3 f_{i,x}f_{i,y}&=2(y^2-x^2+1)(-4xy)+(-4xy)2(x^2-y^2+1)+(-4x)(-4y)\\ &=(-8xy)\left[y^2-x^2+1+x^2-y^2+1\right]+16xy\\ &=-16xy+16xy=0 \end{align} In conclusion, we see that $df_{(x,y)}(u)\cdot df_{(x,y)}(v)=\frac{4(u\cdot v)}{(x^2+y^2+1)^2}$. In particular, letting $\lambda(x,y)=\frac{4}{(x^2+y^2+1)^2}$, we see that $f$ is a conformal map. Note that the function $\lambda(x,y)$ is the stretching factor. Near the origin, thinking about the projection itself, we expect little to no stretching. Indeed, we see that near the origin $f$ is nearly the identity, only scaling by $4$. On the other hand, when $(x,y)$ is very far away from the origin, making $x^2+y^2$ large, $\lambda(x,y)$ is very small. This is because the projection is trying to squeeze all of these far away points in $\mathbb{R}^2$ near the north pole of $S^2$.