changing from global to local coordinate structure

I am going to put this question on math exchange and stack exchange as well (since it doesn't really fit in any of the specific fields, I don't know which it should go in), but here is the question:

I am seeking to evaluate spherical harmonics at a given point, but according to a local coordinate system. Let the actual point be called $(x_g, y_g, z_g)$ (the 'global' point). Now, I have an atom centered at, say, $(x_{at}, y_{at}, z_{at})$ and it has a local coordinate system defined in the following way:

(1) The z axis is defined by the vector from the given atom center to a chosen atom (so, for example, atom1's z axis is defined as the vector from atom1 to atom5)

(2) The y axis is defined by the orthogonal vector lieing in the plane defined by a second atom (so, for example, atom1's y axis is defined to lie in the plane defined by the first vector and the vector from atom1 to atom3). This is done via a simple gram schmidt orthogonalization.

(3) The x axis is defined by a simple cross product of the previous two vectors.

(note that all vectors are normalized).

Now, finding this information will leave you with a 3x3 matrix defined the points which define the 'end points' of the given axes. So the Local Coordinate matrix (for a given atom) may be defined as:

$$\left( \begin{array}{ccc} x_x & x_y & x_z \\ y_x & y_y & y_z \\ z_x & z_y & z_z \end{array} \right)$$

where it is assumed that the 'zero' point (or the center of the local coordinate system) is $(x_{at}, y_{at}, z_{at})$.

Now, given this information, what is the best way to evaluate spherical harmonics (or any function really) in the local coordinate system of a point in the global coordinate system. So, for example, I want to evaluate a point $(x_g, y_g, z_g)$ but with the 'x' direction defined by $(x_x, x_y, x_z)$ and so on.

$$\pmatrix{x_g \\ y_g \\ z_g} = \mathbf{R}_3(180^\circ - \lambda)\mathbf{R}_2(90^\circ - \phi)\mathbf{P}_2\pmatrix{x_x \\ x_y \\ x_z}$$
where $\lambda$ is the longitude of $(x_{at}, y_{at}, z_{at})$ and $\theta$ is the latitude. The angles 180 and 90 you will need to adapt to your problem, but the logic is the same. Matrix $\mathbf{P}_2$ is used to reflect the $y$ axis because the global system is right-handed and the Local Astronomical system is left-handed. If that is not your case, then just leave the matrix out.
You will have to transform each of your global points to the local coordinate system. From your description you really just have a translation and rotation. $$x_{local} = U(x_{global} -T)$$ where $U$ is unitary and $T$ is a translation vector. Now you have $Y_{m,l}^{local}(\theta(x_{local}),\eta(x_{local}))$.