import matplotlib.pyplot as plt
import numpy as np
import math
from mpl_toolkits.mplot3d import Axes3D

gold_ratio = (1+math.sqrt(5))/2
numpoints = 3000
turnfraction = gold_ratio
defaultcolour = 'blue'
highlightcolour = 'orange'
highlightoffset = 0
highlight = 3
power = -0.4794824

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
for i in range(numpoints):
    t = i/round(numpoints, 1)
    inclination = math.acos(1-2*t)
    azimuth = 2 * math.pi * turnfraction * i

    x = math.sin(inclination) * math.cos(azimuth)
    y = math.sin(inclination) * math.sin(azimuth)
    z = math.cos(inclination)

    colour=defaultcolour
    if ((i+highlightoffset)%highlight==0):
        colour=highlightcolour
    ax.scatter(x,y,zs=z,zdir='z',color=colour,s=1)
plt.show()