N_FLAT_BINS=40;dlong = pi*(3-sqrt(5.0));  % ~2.39996323dz    = 2.0/N_FLAT_BINS;longitude = 0;z    = 1 - dz/2;directions=zeros(3,N_FLAT_BINS);for  k = 1:N_FLAT_BINS    r = sqrt(1-z*z);    v=[1*cos(longitude)*r 1*sin(longitude)*r z]‘;    directions(:,k)=v;    z = z - dz;    longitude = longitude + dlong;endfigure(1)plot3(directions(1,:),directions(2,:),directions(3,:),‘r.‘,0,0,0,‘b.‘), hold on, axis equalfor  k = 1:N_FLAT_BINS-1    plot3([directions(1,k) 0],[directions(2,k) 0],[directions(3,k) 0],‘y-‘,‘LineWidth‘,1);    plot3([directions(1,k) directions(1,k+1)],[directions(2,k) directions(2,k+1)],[directions(3,k) directions(3,k+1)],‘LineWidth‘,1);end