Convert between spiral coordinates and axial coordinates in hexagon grid

I'm making a hexagon tile based game, which generates a honeycomb arrangement of hexagons in a spiral shape. My goal is to convert a world position (e.g. a mouse click) to a spiral coordinate (index / layer / position around layer). I can't think of how to do this, so instead I've been looking to simplify the problem by converting to/from axial coordinates first. How can this be done?

My configuration is pointy-topped as follows:

And here are the spiral indexes:

Axial Coordinates for reference:

I already have these equations for spiral coordinates:

const Layer = (index: number): number => {
    return Math.floor((3 + Math.sqrt(12 * index - 3)) / 6);
}
const Position = (index: number, layer: number): number => {
    return index - 3 * layer * (layer - 1) - 1;
}
const Index = (layer: number, position: number): number => {
    return 3 * layer * (layer - 1) + 1 + position;
}

I was looking for an answer to this question which did not require conditionals. The solution I've come up with is based on the spiral algorithm from Amit Patel's Hex Grid Guide.

It relies on the fact that you can sum and scale axial coordinates. A hex has 6 neighbors which have the relative axial coordinates of [(1, 0), (1, -1), (0, -1), (-1, 0), (-1, 1), (0, 1)]. You can define any hex as a sum of two of these relative axial coordinates. I don't know quite how to explain this in words so diagram:

My diagram is for "flat-topped" hexagons, but the logic is the same for "point-topped". The steps to calculate the axial position are.

1. Calculate the layer number, and the position.
2. Add one to the position and modulus by 6*layer. This isn't strictly necessary, but makes it follow the same index structure from before.
3. divide the position by the layer and floor. This gives you the direction "around" the layer. In the diagram this is the vector [0, -1].
4. add 4 to the direction and modulus by 6. Depending on if you do flat-topped or point-topped this may change so experiment. Ensure that your list of axial neighbor vectors are in a "loop". The list I shared above works perfectly. This gives you the direction "out to" the layer. In the diagram this is the vector [-1, 0].
5. Multiply the "out to" vector by the layer number, and the "around" vector by the position modulus the layer number.
6. Sum the two final vectors together, and ta da axial coordinate!

A Python implementation of this could be:

neighbor_vec = ((1, 0), (1, -1), (0, -1), (-1, 0), (-1, 1), (0, 1))

layer = int(3 + Math.sqrt(12 * index - 3)) // 6
position = (index - 3 * layer * (layer - 1)) % (6 * layer)

direction = position // layer
along = neighbor_vec[direction]
out_to = neighbor_vec[(direction + 4) % 6]

along_scale = position % layer

q = along[0]*along_scale + out_to[0]*layer
r = along[1]*along_scale + out_to[1]*layer