Let's say that your surface is defined by $\vec{r}(\Theta) = [r_0 + s f(\Theta) ] \hat{r}$, where $\Theta$ represents the spherical angles. I'm going to conjecture that the normal vector in radial coordinates is proportional to $\vec{n} = (-s\nabla f,1)$, where the gradient is a 2-vector in the angle space, and the last component is the radial direction.
We can verify this by making sure the normal vector is orthogonal to the surface. Note that a tangent vector to the surface is $(1,s\nabla f)$. This tangent points along the gradient direction. A normal vector orthogonal to this one points along an isocontour, which by construction is $(\hat{\Theta}^\perp,0)$, where $\hat{\Theta}^\perp$ is 90 degree rotation of the unit angular vectors in the tangent plane to the sphere passing through the point. The dot product with both of these is zero; the first is obvious, while for the second, $\nabla f \cdot \hat{\Theta}^\perp = 0$.
TL;DR: In radial coordinates: $(-s\nabla f,1)$, then just normalize and convert to cartesian if needed.
Elaboration: Based on your comment, it seems that your $f(\vec{x})$ is a scalar function of a 3-vector. Then $\nabla f$ is a 3-vector. The notation I use lumps the spherical coordinates $(\phi,\theta)$ into a single abstract 2-vector-like quantity $\Theta$, so that when you say $f(\Theta)$, that means $f$ is a function of only the spherical coordinate angles (and not a function of radius). Then $\nabla f(\Theta)=\nabla_\Theta f(\Theta)$ is a gradient in this 2-dimensional angle space. Your $\nabla f$ needs to be projected onto the surface of a sphere first, since you only ever sample $f(\vec{x})$ for $\lVert \vec{x} \rVert = 1$.
So, to calculate your normals, given a point $\vec{P}(\vec{x}) = [R + s f(\vec{x})]\vec{x}$ with $\vec{x}$ such that $\lVert \vec{x} \rVert = 1$,
- Let $\vec{g}(\vec{x}) = \nabla f(\vec{x})$. Here, $\vec{g}$ is the true gradient of $f$ in 3-space, and you can calculate that in cartesian coordinates, so you get $\vec{g} = (g_x,g_y,g_z)$.
- Project out the radial component of $\vec{g}$ to get $\vec{h}(\vec{x})$. To do this, $\vec{h} = \vec{g} - \frac{\vec{g}\cdot\vec{x}}{\vec{x}\cdot\vec{x}}\vec{x}$. Note that the denominator should be 1. Here, $\vec{h}$ represents the component of $\vec{g}$ that should be tangential to a sphere centered at the origin and passing through $\vec{P}$.
- An outward normal vector to the surface is $\vec{n} = \vec{x} - s\cdot \vec{h}$. Normalize it to get a unit normal vector.
All these computations can be done in Cartesian coordinates, but notice that I never had to resort to referring to the Cartesian components of any vector; the basic operations are vector arithmetic and dot products.
I figured it out, since the sphere is centered at the origin, the cross product between $\vec{N}^{\,}$ and $\vec{Z}^{\,}$ will result in $\vec{R}^{\,}$.
Best Answer
If you know the normal vector $\vec N=(n_1,n_2,n_3)$ of a surface at a point $P=(x_0,y_0,z_0)$ then the plane tangent to the surface at $P$ is given by $n_1(x-x_0)+n_2(y-y_0)+n_3(z-z_0)=0$. This is just the equation of a plane given a point and it's normal vector (you should look into that). Any vector in this plane is tangent to the surface. This holds for any surface as long as you know the normal vector at $P$. If you don't know what the normal vector is then you have to use partial derivatives.