Your formulas have already been confirmed. But as for the meaning behind the trig: "Ask Dr. Math" has a derivation of that bearing forumla at http://mathforum.org/library/drmath/view/55417.html in the "Date: 06/19/2002 at 08:45:23" message. After reading it a few times, the steps I see:

1. Given the points on the sphere NorthPole, Location1=(LON1, LAT1), and Location2=(LON2, LAT2)
2. In general, a (LON,LAT) pair turns into (x,y,z) = ( cos(LAT)*cos(LON), cos(LAT)*sin(LON), sin(LAT)*1 ) on a sphere
   → we're wrongly assuming earth is a sphere, for the sake of calculation: it's close enough that you won't accidentally end up "unemployed, in Greenland"</Vizzini>)
3. Rotate the earth (or, easier, your mathematical frame of reference) so that Location1=(0, LAT1), and Location2=(LON2-LON1, LAT2)=(DLON, LAT2)
   → this allows you to drop out a couple sines and cosines by knowing that cos(0)=1 and sin(0)=0
4. Your coordinates become

       N = <NORTH> = ( 0 , 0 , 1 )
       A = <LOC#1> = ( cos(LAT1), 0, sin(LAT1) )  → this is greatly simplified by LON1 = 0
       B = <LOC#2> = ( cos(LAT2)*cos(DLON), cos(LAT2)*sin(DLON), sin(LAT2) )
   

5. Create two normal vectors for two planes (one plane thru N and A; the other thru A and B) by taking the cross products p=NxA and q=BxA (the direction of q keeps the bearing oriented properly)

     p = NxA = (0, cos(LAT1), 0)
     q = BxA = (sin(LAT1)*cos(LAT2)*sin(DLON),
            cos(LAT1)*sin(LAT2)-sin(LAT1)*cos(LAT2)*cos(DLON),
            -cos(LAT1)*cos(LAT2)*sin(DLON))
   

6. The angle between these planes is your bearing. Normally that angle would be from cos(BEARING) = p DOT q / (|p|*|q|). But because NxA is all in the y-direction of your rotated earth, you can simplify to the arctangent function, already used in your code.