Using the question posted by @einbandi and answer posted by @acl here the general solution seems to be:
Module[{vars}, mat = {(*Some code that constructs a matrix out of expressions of vars*)} fmat = Function @@{{vars}, mat}; ] Applying to my specific case, the following works:
Module[{\[Theta], \[Phi]}, row1 = {}; row2 = {}; For[ii = 1, ii <= Length[modes], ii++, {l, m} = modes[[ii]]; G1lm = sglm[-2, l, m, \[Theta]] + (-1)^l sglm[-2, l, -m, \[Theta]]; G2lm = sglm[-2, l, m, \[Theta]] - (-1)^l sglm[-2, l, -m, \[Theta]]; AppendTo[row1, G1lm Cos[m \[Phi]]]; AppendTo[row1, G1lm Sin[m \[Phi]]]; AppendTo[row2, G2lm Sin[m \[Phi]]]; AppendTo[row2, -G2lm Cos[m \[Phi]]]; ]; proj = Function @@ {{\[Theta], \[Phi]}, {row1, row2}}; ] Enclosing the construction of the matrix in a Module ensures that \[Theta] and \[Phi] are protected, and the last line turns proj into a function of the two variables. This method is also fast: it's about twice as fast as the replacement method I was originally using.
I'm sure there's much more concise ways to construct the matrix than the For loop I'm using here, but since this part of the code is only called once, any computational time savings are negligible.
