I am trying to understand the following code in this answer regarding Möbius transformation:

(* Projection from the sphere to the plane *)
stereo = Compile[{{xyz, _Real, 1}, {XYZ, _Real, 1}}, Module[{
     r = Sqrt[(xyz[[1]] - XYZ[[1]])^2 + (xyz[[2]] - XYZ[[2]])^2],
     theta = ArcTan[(xyz[[1]] - XYZ[[1]]), (xyz[[2]] - XYZ[[2]])]},
    {(r (1 + xyz[[3]]))/(1 - XYZ[[3]] + xyz[[3]]) Cos[theta + Pi] + xyz[[1]],
     (r (1 + xyz[[3]]))/(1 - XYZ[[3]] + xyz[[3]]) Sin[theta + Pi] + xyz[[2]], 0}]];

I have read the documents for the functions Compile and Module, but still cannot figure out how the code works.

Could anyone elaborate the following?

  • What does {xyz,_Real,1} do and what does xyz[[1]] mean?

    [Added:] In the document for Compile, it is said that

    Compile[{{x1,t1,n1},…},expr] assumes that xi is a rank ni array of objects, each of a type that matches ti.

    But what is 1 in {xyz,_Real,1}? If it means rank 1, then why later it is written that xyz[[1]], xyz[[2]], xyz[[3]], which suggests that xyz is an array?

  • How do Compile and Module work together to give the definition of a function?

  • What is the formula for the function that this code is really defining? What is the input and what it is the output?
  • $\begingroup$ Part and Compile will pretty quickly help to explain your first bulleted question. $\endgroup$
    – ktm
    Commented Feb 13, 2017 at 20:07
  • 1
    $\begingroup$ Be sure to always read under "Details and Options" on documentation pages. You can select [[ and press F1 (or Command-Shift-F) to see what it is (it's array indexing). The formula is exactly as written there, just substitute in r and theta from the first part of the Module. xyz in this case is a three-element real array representing a 3D vector with x, y, z components. $\endgroup$
    – Szabolcs
    Commented Feb 13, 2017 at 20:07
  • $\begingroup$ @Szabolcs: Thank you for your comment. I'm slowly learning things and I might come back later with my own answer. I would like to leave the question undeleted (if it is not a very stupid trivial one :-)). $\endgroup$
    – user664
    Commented Feb 13, 2017 at 21:13
  • 1
    $\begingroup$ Rank 1 is an array, you're thinking about rank 0 which is a point. Rank 2 is a matrix, etc. Module does the same thing inside Compile as it does outside Compile. $\endgroup$
    – C. E.
    Commented Feb 14, 2017 at 2:19
  • $\begingroup$ @C.E.: the rank is indeed confusing at first. Thanks for pointing that out. $\endgroup$
    – user664
    Commented Feb 14, 2017 at 2:29

1 Answer 1

  • As the document says,

    Compile[{{x1,t1,n1},…},expr] assumes that xi is a rank ni array of objects, each of a type that matches ti.

    Hence, {xyz,_Real,1} means xyz is a rank 1 array of "real numbers". Since xyz is a vector, xyz[[1]] means the first component of xyz. Therefore, if xyz={2,4,3}, then xyz[[1]]=2.

  • The expr in Compile[{{x1,t1,n1},…},expr] is given by the block of Module, a simple example of which for seeing how it works could be

  • Here is how the function really works. Given $x=(x_1,x_2,x_3)$ and $y=(y_1,y_2,y_3)$ in $\mathbb{R}^3$, the output of the function is $$ \bigg(\frac{r(1+x_3)}{1+x_3-y_3}\cos(\theta+\pi)+x_1, \frac{r(1+x_3)}{1+x_3-y_3}\sin(\theta+\pi)+x_2,0\bigg) $$ where $$ r:=\sqrt{(x_1-y_1)^2+(x_2-y_2)^2},\quad \theta:=\arctan\frac{x_2-y_2}{x_1-y_1}. $$

    The coordinate of the center of the sphere $S$ is given by $x=(x_1,x_2,x_3)$ and the radius of the sphere is $1$. This makes sense since the output of the function is of the form $(R\cos\psi+x_1,R\sin\psi+x_2,0)$, which implies that $x_1,x_2$ are the first two coordinates of the center of the sphere.

    On the other hand, $y=(y_1,y_2,y_3)$ is the coordinate of a point on the sphere. The output of the function gives the coordinate of the "sterographic projection" of $y$ (the intersection of the line (which goes through $y$ and the north pole of the sphere) with the complex plane). When $x=(0,0,0),$ the function coincides with the one defined in Wikipedia.

One needs to be careful with the definition $\theta$ since in Mathematica, ArcTan[x,y] is different from the function $\arctan\frac{y}{x}$. For instance one has in Mathematica ArcTan[1,1]=$\frac{\pi}{4}$ and ArcTan[-1,-1]=$\frac{-3\pi}{4}$. This is why $\theta+\pi$ instead of $\theta$ is used in the output of the function. For the sake of a function, one could define $\theta$ as the angle of the vector $(y_1-x_1,y_2-x_2,0)$ with the $X_1$-axis and replace $\theta+\pi$ with $\theta$ in the output of the function.

In the code, one could replace theta = ArcTan[(xyz[[1]] - XYZ[[1]]), (xyz[[2]] - XYZ[[2]])] with theta = ArcTan[(XYZ[[1]] - xyz[[1]]), (XYZ[[2]] - xyz[[2]])] and theta + Pi in the Cos and Sin functions with theta.

Having not checked it though, I think one can simply replace Cos[theta+Pi] with XYZ[[1]]-xyz[[1]] and Sin[theta+Pi] with XYZ[[2]]-xyz[[2]] so that one does not need to introduce the variable theta at all.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.