Here's my quick&dirty solution. It's a bit similar to @azdahak's answer, but it uses an approximate mapping instead of cylindrical coordinates. On the other hand, there are no manually adjusted control parameters - the mapping coefficients are all determined automatically:
The label is brigh...
OK. I'd prefer that he do the summary, if he doesn't, I can clean up that SO answer...
It doesn't feel correct to just summarize someone else's answer for a day before, and get rep (though I don't care about rep, really). My motivation for answering was to popularize the site a bit.
It looks alright I think but obviously lacks something without the proper labels for the high-symmetry points. Do you have a good way of getting those, or just work them out for each case as required and do the Ticks manually?
@OleksandrR I work them out manually, but I have a couple of systems set up for use with analytical tight binding that given a sequence of points (i.e. a path in k-space), I can plot it. I used some of the ideas in generating the ticks for my plot. I'll send you them, if you want.
Anyone have any idea why this does not evaluate:NDSolve[{(\[Rho][s]^(5/3) + 5/3 s \[Rho][s]^(2/3) Derivative[1][\[Rho]][s])/ s + \[Rho][s] s == 0, \[Rho][0] == \[Rho]0, \[Rho]'[0] == 0}, \[Rho][s], {s, 0, 10}]