Using the spherical coordinates, create the icosahedral tessellation and project it with a suitable projection for your work, then apply any GIS to solve this point in polygon problem.

If you want to write your own code, note that each pair of vertices determines a reflection plane (containing the origin) and that points can be classified according to which side of each of the 15 distinct reflection planes they lie on. This classification is actually a finer one than requested, because each of the 20 faces of the icosahedron is further divided into six congruent triangles, giving 120 regions (equal to a fundamental domain for the action of the reflections together with its 119 images).

For the record, a set of normal vectors to these planes (in geocentric Cartesian coordinates) is

X   Y   Z
1-f -f  1
f   -1  f-1
-2  0   0
f-1 -f  1
-f  -1  f-1
-1  f-1 f
1-f f   1
0   -2  0
-1  1-f f
-f  1   f-1
1   1-f f
0   0   2
f-1 f   1
f   1   f-1
1   f-1 f

where f = (1 + sqrt(5))/2 = 1.61803398874989... is the Golden Ratio. (These all point towards the Equator or northwards.) To use them, convert any point (on the sphere or off) into geocentric Cartesian coordinates (x,y,z) and form the 15 dot products with the 15 normals (which essentially multiplies the 15 by 3 matrix by the 3 by 1 vector (x,y,z)'). Classify those products as negative or non-negative, yielding a 15-element Boolean vector. Considering that as a 15-bit binary integer, there are only 120 possible values, one for each of the 120 little triangles. These values encode the desired information about the face (and which part of the face) where the point (x,y,z) is located.

To illustrate, I associated 120 random colors with these codes to display 10,000 random points on the unit sphere:

Both the icosahedron and its dual, the dodecahedron, are evident.

Using the spherical coordinates, create the icosahedral tessellation and project it with a suitable projection for your work, then apply any GIS to solve this point in polygon problem.

If you want to write your own code, note that each pair of vertices determines a reflection plane (containing the origin) and that points can be classified according to which side of each of the 15 distinct reflection planes they lie on. This classification is actually a finer one than requested, because each of the 20 faces of the icosahedron is further divided into six congruent triangles, giving 120 regions (equal to a fundamental domain for the action of the reflections together with its 119 images).

For the record, a set of normal vectors to these planes (in geocentric Cartesian coordinates) is

X   Y   Z
1-f -f  1
f   -1  f-1
-2  0   0
f-1 -f  1
-f  -1  f-1
-1  f-1 f
1-f f   1
0   -2  0
-1  1-f f
-f  1   f-1
1   1-f f
0   0   2
f-1 f   1
f   1   f-1
1   f-1 f

where f = (1 + sqrt(5))/2 = 1.61803398874989... is the Golden Ratio. (These all point towards the Equator or northwards.) To use them, convert any point (on the sphere or off) into geocentric Cartesian coordinates (x,y,z) and form the 15 dot products with the 15 normals (which essentially multiplies the 15 by 3 matrix by the 3 by 1 vector (x,y,z)'). Classify those products as negative or non-negative, yielding a 15-element Boolean vector. Considering that as a 15-bit binary integer, there are only 120 possible values, one for each of the 120 little triangles. These values encode the desired information about the face (and which part of the face) where the point (x,y,z) is located.

To illustrate, I associated 120 random colors with these codes to display 10,000 random points on the unit sphere:

Both the icosahedron and its dual, the dodecahedron, are evident.