+struct
+{
+ int n;
+ double darkness;
+ double power;
+ double density;
+ int g;
+ double gpower;
+ double gfactor;
+ double gdpower;
+}
+color_starfield_parms;
+typedef struct
+{
+ double x, y, z, e;
+ double R, G, B, A;
+} starfield_t;
+starfield_t *starfield = NULL;
+int starfield_cmp(const void *a_, const void *b_)
+{
+ const starfield_t *a = a_;
+ const starfield_t *b = b_;
+ if(a->z < b->z)
+ return -1;
+ if(a->z > b->z)
+ return +1;
+ return 0;
+}
+void color_starfield(double x, double y, double z, double *r, double *g, double *b)
+{
+ int i, j, k;
+ double f, mindot, d;
+
+ if(!starfield)
+ {
+ fprintf(stderr, "Initializing starfield...\n");
+ starfield = malloc(sizeof(*starfield) * color_starfield_parms.n);
+ for(i = 0; i < color_starfield_parms.n; ++i)
+ {
+ double r;
+ do
+ {
+ starfield[i].x = rnd() * 2 - 1;
+ starfield[i].y = rnd() * 2 - 1;
+ starfield[i].z = rnd() * 2 - 1;
+ r = starfield[i].x * starfield[i].x
+ + starfield[i].y * starfield[i].y
+ + starfield[i].z * starfield[i].z;
+ }
+ while(r > 1);
+
+ starfield[i].R = rnd();
+ starfield[i].G = rnd();
+ starfield[i].B = rnd();
+ f = starfield[i].R * 0.299 + starfield[i].G * 0.587 + starfield[i].B * 0.114;
+ starfield[i].R /= f;
+ starfield[i].G /= f;
+ starfield[i].B /= f;
+ starfield[i].A = rnd();
+ starfield[i].e = color_starfield_parms.density * pow(rnd(), -color_starfield_parms.power);
+ }
+ fprintf(stderr, "Gravitating starfield...\n");
+ for(k = 0; k < color_starfield_parms.g; ++k)
+ {
+ i = rand() % color_starfield_parms.n;
+ j = rand() % color_starfield_parms.n;
+ f = pow(rnd(), color_starfield_parms.gpower);
+ f = f * color_starfield_parms.gfactor;
+ d = pow(starfield[j].x - starfield[i].x, 2)
+ + pow(starfield[j].y - starfield[i].y, 2)
+ + pow(starfield[j].z - starfield[i].z, 2);
+ f *= pow(1 / (1 + d), color_starfield_parms.gdpower);
+ starfield[i].x += f * (starfield[j].x - starfield[i].x);
+ starfield[i].y += f * (starfield[j].y - starfield[i].y);
+ starfield[i].z += f * (starfield[j].z - starfield[i].z);
+ double r = starfield[i].x * starfield[i].x + starfield[i].y * starfield[i].y + starfield[i].z * starfield[i].z;
+ r = sqrt(r);
+ starfield[i].x /= r;
+ starfield[i].y /= r;
+ starfield[i].z /= r;
+ }
+ fprintf(stderr, "Sorting starfield...\n");
+ qsort(starfield, sizeof(*starfield), color_starfield_parms.n, starfield_cmp);
+ fprintf(stderr, "Done.\n");
+ }
+
+ *r = *g = *b = 0;
+ mindot = pow(1/256.0, 1.0/color_starfield_parms.density);
+ for(i = 0; i < color_starfield_parms.n; ++i)
+ {
+ double dot = x * starfield[i].x + y * starfield[i].y + z * starfield[i].z;
+ if(dot <= mindot)
+ continue;
+ double f = pow(dot, starfield[i].e) * starfield[i].A;
+ *r += starfield[i].R * f;
+ *g += starfield[i].G * f;
+ *b += starfield[i].B * f;
+ }
+ // make fit in 0..1
+ *r = *r / (color_starfield_parms.darkness + *r);
+ *g = *g / (color_starfield_parms.darkness + *g);
+ *b = *b / (color_starfield_parms.darkness + *b);
+}
+