For the hexagonal grid maths I used this article: Coordinates in Hexagon-Based Tile Maps, which is an easy read and explains things quite well, although I did have to correct some of the code because the pseudo-code was incorrect.
This is what the demo looks like:
I won't be showing the code here because it is a direct implementation of the pseudo-code described in the article I've linked above. There are some optimisations like pre-calculating sin(30deg) and cos(30deg). There is still quite a lot of room for improvement in terms of performance, for example r*2 and h+s are calculated several times when they don't need to be.
So without further delay, here's the demo...