Simulating rnorm() using runif() I am trying to 'simulate' rnorm() using only runif().
I don't know if I should do:
sqrt(-2*log(U1))*cos(U2)

or
sqrt(-2*log(U1))*sin(U2)

Where U1 is a runif(0,1) and  U2 runif(0,6.28)
I do not know if I should do it using cos or sin, or is it that I need to sample from  one and the other consecutively? What is the mathematical logic behind it?
 A: It looks like you are trying to use the Box-Muller transform. 
The method is to use $U_1, U_2$ which are independently distribution Uniform(0,1). Then 
$Z_0 = \sqrt{-2 \ln(U_1)}\cos(2 \pi U_2) $
and
$Z_1 = \sqrt{-2 \ln(U_1)} \sin(2 \pi U_2) $
are a pair of independent $N(0,1)$. 
The derivation comes from using polar coordinates, but to be honest the exact steps escape me at the moment (see link if you are interested).
A: Let $N_1,N_2$ be independent $N(0,1)$ random variates. $(N_1,N_2)$ defines a point in Cartesian coordinates. We transform to polar by
$$ N_1 = R \cos \theta$$
$$ N_2 = R \sin \theta$$
This transformation is one to one and has continuous derivatives. So we can derive the joint distribution $f_{R,\theta}(r,\theta)$ by using 
$$f_{R,\theta}(r,\theta) = f_{N_1,N_2}(n_1,n_2)|J^{-1}|$$
where $J$ is the Jacobian of the transformation and
$$ J^{-1}(r,\theta )=
{\begin{bmatrix}{\dfrac {\partial n_1}{\partial r}}{\dfrac {\partial n_1}{\partial \theta }}\\{\dfrac {\partial n_2}{\partial r}}{\dfrac {\partial n_2}{\partial \theta }}\end{bmatrix}} 
= \begin{bmatrix}\cos \theta & -r\sin \theta \\\sin \theta & r\cos \theta \end{bmatrix}$$
Hence $$ |J^{-1}| = r \cos^2 \theta + r \sin^2 \theta = r   $$
and 
$$ f_{R,\theta}(r,\theta) = r \frac{1}{2\pi} e^{-\frac{1}{2}(n_1^2+n_2^2)} = \boxed{\frac{r}{2 \pi}e^{-\frac{1}{2}(r^2)}} \qquad 0 \leq \theta \leq 2 \pi, 0 \leq r \leq \infty \ $$
And we see $\theta \sim unif[0,2\pi]$ (since $f_{\Theta}(\theta) = \frac{1}{2\pi}$ ) and $R^2 \sim \exp(1/2)$
Hence, to simulate $\Theta$ we simply take $2 \pi U_2$ where $$U_2 \sim unif[0,1]$$ and to simulate $R$ we can take $$-2\ln U_1$$ where $U_1 \sim unif[0,1]$. ( to see why find the density of $X=-\ln U, U \sim  unif[0,1]$). 
So, Box and Muller simply inverted $N_1= R \cos \theta$, $N_2= R \sin \theta$ and moved from $(R,\Theta)$ to $(N_1,N_2)$ by simulating $\Theta$ from  $2 \pi U_2$, and an independent $R$ from $ \sqrt{- 2 \ln U_1}$ 
explicitly 
$$Z_0 = \sqrt{-2 \ln(U_1)}\cos(2 \pi U_2)$$
$$Z_1 = \sqrt{-2 \ln(U_1)} \sin(2 \pi U_2)$$
And that is that is the mathematical logic behind it. 

p.s. As I don't know whether it is clear through the mathematical justification (I believe it is, but you may have been lost in details),  note you need a pair of drawings from the uniform to get a pair of normals 
