Charcoal, 162 146 bytes

Ｎθ⊞υ⟦⁰¦¹⊘⊕₂⁵⟧≔⟦⟧ηＦυＦＥ²⊞ＯΦιν×⊖⊗κ§ι⁰Ｆ¬№υκＦΦ⊞Ｏυκ⁼⁴ΣＸＥμ⁻ξ§κπ²⊞η⟦κμ⟧ＦηＦ…¹θ⊞υＥ§ι⁰⁺×κλ×⁻θꧧι¹μＦηＦΦυ⬤ι№η⟦μκ⟧Ｆ…²θＦ…¹λ⊞υＥ§ι⁰⁺⁺×μν×⁻λ짧ι¹ξ×⁻θλ§κξＩＥυ∕ι₂ΣＸι² 

Try it online! Link is to verbose version of code. Explanation:

Ｎθ 

Input b.

⊞υ⟦⁰¦¹⊘⊕₂⁵⟧≔⟦⟧η 

Start with one vertex and no edges of an icosahedron of side 2.

Ｆυ 

Process all vertices as they are discovered.

ＦＥ²⊞ＯΦιν×⊖⊗κ§ι⁰ 

Generate more vertices by rotating about the line x=y=z and reflecting in the xy plane. (This is the same procedure as for the regular dodecahedron in the linked question.)

Ｆ¬№υκ 

Ignore previously discovered vertices.

ＦΦ⊞Ｏυκ⁼⁴ΣＸＥμ⁻ξ§κπ²⊞η⟦κμ⟧ 

Save the new vertex and also discover any edges by checking for previously found vertices that are 2 away. (Very fortunately, ϕ²+1²+(ϕ-1)²=2² under floating-point arithmetic.) Since only previously found vertices are checked, each edge is only detected once, and has a specific direction.

Ｆη 

Loop over the edges.

Ｆ…¹θ 

Loop over the intermediate points.

⊞υＥ§ι⁰⁺×κλ×⁻θꧧι¹μ 

Extrapolate a point along the edge. (Sadly Charcoal doesn't have a way of adding two vectors which would really simplify this. Maybe if it had octonions...)

Ｆη 

Loop over the edges again.

ＦΦυ⬤ι№η⟦μκ⟧ 

Loop over the vertices that have edges to both vertices of that edge. (Note that because the edges have a specific direction, each face can only be detected once.)

Ｆ…²θＦ…¹λ 

Loop over the interior points of the face.

⊞υＥ§ι⁰⁺⁺×μν×⁻λ짧ι¹ξ×⁻θλ§κξ 

Extrapolate a point inside the face. (Due to the normalisation, it's not necessary to interpolate.)

ＩＥυ∕ι₂ΣＸι² 

Normalise all the points.

Pyth, 99 bytes

mcR.add+Ksm.>L+0d3*_B1_B.n3s+mm+V*Rkhd*R-QkedStQJf&SITq.aT2^K2mmsMC.b*RlNYkdfq3lT./UQf.A}RJ>3PyT^K3 

Try it online!

We begin with the vertices of the regular icosahedron made by all cyclic permutations and sign flips of \$(0,1,\phi)\$. Then we find all unordered pairs of edges that have distance 2, interpolate along the edge they specify. Then we find all unordered triplets where all the pairs are in the list of edges, and using partitions of the input we interpolate along the face they specify. Finally we normalize all points.

Ruby, 161 153 148 146 bytes

->n{(n*n*10+2).times{|k|b=k/10%n e=b-a=k/10/n r=n*c=1i**0.4 r+=e<0?b/c-c*a:a/c-b*c+e p 2*r*c**k/s=(r.abs*2i+z=n-e*=e<0?5**0.5-1:2).abs,z*~0**k/s}} 

Try it online!

Golfed version of the original 175-byte concept. Now competing with the golfing languages!

• Eliminated variable j=k/10. Introduced new variable e=a-b.

• a & b are now integers (not normalised by dividing by m=n*1.0, with other values scaled up by n accordingly.)

• Conditionals used for selection instead of arrays. d was revised true/false instead of 0/1, then eliminated completely. e<0? (used twice) is overall shorter than d ? because it doesn't require whitespace, or prior definition of d=0>e.

• Instead of using pythagoras (squaring then square root) to calculate scale factor s, multiply r.abs by 1i, add z and take abs of the result. This allowed z to be calculated on the same line as the scale factor calculation (but means the output is now x+yi before z.)

• Moved ~0**k to allow deletion of parentheses. Similarly redistributed calculation of c and r (Multiplying r.abs by 2i instead of 1i as previously allows movement of 2*r*c**k and deletion of parentheses.)

• Removed [] around output. This alters readability because it prints each point over 2 lines instead of 1, but saves 2 bytes.

Ruby, 175 bytes

first working version

->n{(n*n*10+2).times{|k|j=k/10;a=j/n/m=n*1.0;b=j%n/m z=(1+(a-b)*[5**0.5-1,2][d=a>b ?0:1])*~0**k r=((c=1i**0.4)+[b/c-a*c,b-b*c+a/c-a][d])*c**k*2 p [z/s=(r.abs2+z*z)**0.5,r/s]}} 

Try it online!

A function that takes an argument n and prints the vertices in the format [z,x+yi]. One pair of opposite vertices is aligned with the z axis; these vertices are printed last. The code generates each vertex exactly once with no duplicates, per OP's comment in the post.

Explanation

According to https://mathworld.wolfram.com/RegularIcosahedron.html

A construction for a icosahedron... places the end vertices at (0,0,+/-1) and the central vertices around two staggered circles of radii 2/sqrt(5) and heights +/- 1/sqrt(5)

This code produces the zcoordinates of the icosahedron at a scaled up size of 0,0,sqrt(5) for the end vertices, which puts the other vertices at z=+/-1. For calculation of the x and y coordinates, the radius of the circle is conveniently set to 1, but the coordinates are doubled for compatability with the z coordinates prior to normalization.

The icosahedron is split into 10 identical diamonds as below. We include only the points marked with an X because the points marked with O are included as part of the adjacent diamond. We do however have to include the poles (marked as ?) once.

The diamond is folded at the join of the two triangles. Looking down the z axis at the x+yi plane, the polar apex is at 0+0i, and the other apex is at 1+0i. The right vertex is at a 1/10th rotation, c=i**0.4 = 0.809+0.588i. The left vertex is at the conjugate, 0.809-0.588i. As the absolute value of c is 1, the conjugate is conveniently generated by 1/c.

 North Pole ? x+yi=0 z=sqrt(5) / \ O X / \ / \ O X X x+yi / \ / \ / \ x+yi =1/c O---X---X---X =c=i**0.4 z=1 =0.809 \ / \ / \ / =0.809 -0.588i O X X +0.588i Note: x,y coordinates at half \ / \ / the scale of z coordinates O X (they are doubled before normalization) \ / O x+yi=1 z=-1 

We work through the n*n points in the diamond with increments in the value of j. we start at [r=x+yi=c,z=1], with increasing value of b=j%n moving diagonally towards the bottom corner and a=j/n moving diagonally towards the top corner. different vectors (selected by d=a>b) must be added to the starting point for the upper and lower triangle.

 x+yi z Upper triangle a>b a*(0-c)+b*(1/c-0) (a-b)*(5**0.5-1) Lower triangle a<=b b*(1-c)+a*(1/c-1) (a-b)*2 

For each value of j we generate 10 images by 1/10 rotation of the value r=x+yi by multiplying by c**k and flipping the z value by multiplying by -1**k (due to priority of operators in Ruby, this is represented as ~0**k)

After n*n*10 iterations, we have covered all points up to a=b=n-1. It remains to add the polar vertices. This just requires another 2 iterations where the value of j is n*n, a=n and b=0.

Finally, the last line normalizes the points to a sphere and outputs them.

Python3, 421 bytes:

def f(l,r=[]): if[]==l:yield r;return for i in l[0]:yield from f(l[1:],r+[i]) C=lambda n=0,b=[0,1,(1+5**.5)/2]:[]if n==3 else[*f([[[i,-i],[0]][i==0]for i in b[n:]+b[:n]])]+C(n+1) def G(w): c,R=C(),[] for x in c: for y in c: if sum((a-b)**2 for a,b in zip(x,y))**.5==2: R+=[x,y] for v in range(1,w):R+=[[((w-v)*a+v*b)/w for a,b in zip(x,y)]] return{tuple(i/sum(j**2 for j in p)**.5 for i in p)for p in R} 

Try it online!