Why is this allowed? ("Fourier's Trick"; finding the coefficients in a Fourier Series)

In my textbook (Introduction to Electrodynamics, D. Griffiths), we derive the equation for some strange potential function. Eventually, we get to this (for $n \in \mathbb{Z}^+$):

$$ V_0(y) = \sum_{n=0}^{\infty} C_n\sin{\frac{n\pi}{a}y} \tag{3.31}$$

Here's where things go awry for me.

... how do we actually determine the coefficients $C_n$, buried as they are in that infinite sum? The device for accomplishing this is so lovely it deserves a name—I call it Fourier's trick, though it seems Euler had used essentially the same idea somewhat earlier. Here's how it goes: Multiply Eq. 3.31 by $\sin{n'\pi y/a}$ (where $n'$ is a positive integer), and integrate from 0 to a:

$$ \displaystyle \sum_{n=0}^{\infty} C_n \int_0^a\sin{\frac{n\pi}{a}y} \sin{\frac{n'\pi}{a}y} dy ~~~=~~~ \int_0^a V_0(y)\sin{\frac{n'\pi}{a}y} dy$$

The answer understandably comes out to something very nice and convenient. But... why is this something you can do? There's no obvious reason for why that doesn't intrinsically change the problem (in the same way that I can say "Multiply both sides by $0$. You've successfully reduced the problem to zero. Well done!)

(While typing out the above, I suspect that it has something to do with the inner product of a function and an orthonormal basis? The infinite $\sin$ functions create an orthonormal basis, and taking that integral over all possible values effectively extracts the coefficients for each basis function. When it is suggested that we multiply by $\sin{\frac{n'\pi}{a}y}$ and integrate, this isn't changing the basis at all, it's just (sneakily) extracting the coefficients, which only exist when $n = n'$ (because the $\sin$ functions are all orthogonal). It's like taking the coefficients of a basis with itself... right?

I think this may be one of those cases where, in the process of asking the question, I figure out the answer—but this is all fairly new to me, and I'd like to ask it anyway for confirmation and, possibly, a clearer explanation).

Your suspicion about inner product is entirely correct. The trigonometic polynomials $\{\sin(nx), \cos(mx)\}$ (possibly translated and scaled) are known to form an orthogonal system with respect to the scalar product given by $\langle f, g\rangle:=\int fg$, and if you choose the function space correctly (usually one uses a space called $L^2$) and use properly chosen norming factors then one can show they actally form a complete orthonormal system $e_k$, which simply means you can express any function in that space as $$ f = \sum_k\langle f, e_k \rangle e_k$$ (where convergence is to be understood in that space wrt the norm derived from the scalar product). The coefficients in that sum are what you are looking at. You may know that kind of representation from the finite dimensional case.

I deliberately did not specify the constants which make the orthogonal system orthonormal, nor an interval as domain of definition -- by translation and scaling you can do something like that on any bounded interval in $\mathbb{R}$ There is also a complex version of this, in which case one would use the scalar product $\int f\bar{g}$, and $\{e^{ikx}\}$ as orthogonal system.

If you want to look up the details, then most introductions to real analysis will have a section on that topic. Rudin's books, for example do explain this.

It might help to break this up into smaller steps.

\begin{align*} &V_0(y) = \sum_{n=1}^{\infty} C_n \sin \left(\frac{n \pi y}{a} \right) \\ \implies & V_0(y) \sin \left(\frac{n' \pi y}{a} \right) = \sum_{n=1}^{\infty} C_n \sin \left(\frac{n \pi y}{a} \right)\sin \left(\frac{n' \pi y}{a} \right) \\ \implies & \int_0^a V_0(y) \sin \left(\frac{n' \pi y}{a} \right) \, dy = \sum_{n=1}^{\infty} C_n \int_0^a \sin \left(\frac{n \pi y}{a} \right)\sin \left(\frac{n' \pi y}{a} \right) \, dy = \frac{a}{2} C_{n'}. \end{align*}

We now solve for $C_{n'}$ to obtain \begin{equation*} C_{n'} = \frac{2}{a} \int_0^a V_0(y) \sin \left(\frac{n' \pi y}{a} \right) \, dy. \end{equation*}