+vec3 rand33(vec3 p3)
+{
+ p3 = fract(p3 * vec3(.1031, .1030, .0973));
+ p3 += dot(p3, p3.yxz+33.33);
+ return fract((p3.xxy + p3.yxx)*p3.zyx);
+}
+
+float stars( vec3 rd, float rr, float size ){
+ vec3 co = rd * rr;
+
+ float a = atan(co.y, length(co.xz)) + 4.0 * PI;
+
+ float spaces = 1.0 / rr;
+ size = (rr * 0.0015) * fwidth(a) * 1000.0 * size;
+ a -= mod(a, spaces) - spaces * 0.5;
+
+ float count = floor(sqrt(pow(rr, 2.0) * (1.0 - pow(sin(a), 2.0))) * 3.0);
+
+ float plane = atan(co.z, co.x) + 4.0 * PI;
+ plane = plane - mod(plane, PI / count);
+
+ vec2 delta = rand33(vec3(plane, a, 0.0)).xy;
+
+ float level = sin(a + spaces * (delta.y - 0.5) * (1.0 - size)) * rr;
+ float ydist = sqrt(rr * rr - level * level);
+ float angle = plane + (PI * (delta.x * (1.0-size) + size * 0.5) / count);
+ vec3 center = vec3(cos(angle) * ydist, level, sin(angle) * ydist);
+ float star = smoothstep(size, 0.0, distance(center, co));
+ return star;
+}
+