Help find the bug in my SPLL code
I've written articles like that, and it's incredibly hard -- and sometimes not the best exposition -- to get the graphs to match exactly with the text.  From the looks of things, the author was intentionally showing you idealized graphics rather than the actual output of code.


  My phase-difference signale_D[n] is flipped vertically from what the tutorial shows...


Eh, don't worry.  It's a control loop, there's always a stray minus sign in there someplace.  As long as you match it up with another stray minus sign so that it's stable, you're OK.


  I have output that is close to what the article shows, except that my cos_out[n] seems to be "leaning" to the left: cos_out


That's because your K_p is so high.  You can get a simple loop, fast settling, or a nice sinusoidal response -- not all three.  If you're using the author's numbers, or if you've selected gains to make your loop settle as fast as his idealized example (which was probably chosen more to make things fit in the page rather than because it's the best tuning) then what's happening is that strong signal at the sum frequency is bleeding into the NCO phase, and making the "lean".

Typically a PLL will be set up to take hundreds or thousands of cycles to settle, which implies a much smaller K_p (certainly less than 0.01, and almost arbitrarily small depending on what you're doing).  How fast a PLL settles depends on what you're doing, and what you're doing it with, so if your next question is "how fast should it settle" -- I can't answer.  But I would say try it for settling in around 100 cycles, 1000, and 10,000.  Keep in mind that the capture range (the frequency error between signal and NCO initial frequency over which it can actually acquire a lock) and its tracking performance (its ability to follow a varying input frequency) both go down with slower settling.