//GPL (as if anyone would want this anyway...)
#include <stdio.h>
#include <math.h>

#define PI 3.14159265358979323846264338327950288

//m>=0
double cfac(int l, int m, double z) {
  int x;
  int i;
  double r, s, t;
  r=log((2*l+1)*(.25/PI))*.5;
  r-=log(2)*l;
  i=(l+m+1)>>1;
  for(x=2;x<=l-i;++x) r-=   log(x);
 if(l-m<i) {
/*for(   ;x<=l-m;++x) r+= 0*log(x);*/ x=l-m+1;
  for(   ;x<=i  ;++x) r-=.5*log(x);
 } else {
/*for(   ;x<=i  ;++x) r-= 0*log(x);*/ x=i+1;
  for(   ;x<=l-m;++x) r+=   log(x);
 }
  for(   ;x<=l+m;++x) r+=.5*log(x);
  for(   ;x<=2*i;++x) r+=   log(x);

  if(m) r+=.5*log(2);

  r=exp(r)*((l^i)&1?-1:1)*(2*i-l-m?z:1);
  s=-z*z;
  t=r;
  for(i=((l+m+1)>>1)+1;i<=l;++i) t+=r*=s*(2*i-1)*(2*i)*(l-i+1)/((double)(2*i-l-m-1)*(2*i-l-m)*i);
  return(t);
}

double cval(int m, double x, double y) {
  int i;
  _Complex double d=x+1i*y, r;
  for(r=1;m;m>>=1,d*=d) if(m&1) r*=d;
  return(__real__ r);
}

void rend(unsigned char *c, int sx, int sy, int pit, int l, int m, double ang) {
  int a, b;
  double x, y, z, sa, ca, f, t, r;
  sa=sin(ang); ca=cos(ang);
  for(b=0;b<sy;++b) {
    for(a=0;a<sx;++a) {
    z=1-b*2/(sy-1.);
      x=a*2/(sx-1.)-1;
      t=1-x*x-z*z;
      if(t<0) t=1e99; else {
        y=-sqrt(t);
        t=cos(.2)*z-sin(.2)*y;
        y=sin(.2)*z+cos(.2)*y;
        z=t;
        t=ca*x+sa*y;
        y=-sa*x+ca*y;
        x=t;
    f=cfac(l, m, z);
        t=cval(m, x, y)*f;//printf("a%lf, %lf, %lf, %lf\n", t, x, y, z);
      }
      r=128+t*512; if(r<0) r=0; else if(r>255) r=255; c[a*3+b*pit  ]=(int)r|31;
      r=128+t*256; if(r<0) r=0; else if(r>255) r=255; c[a*3+b*pit+1]=(int)r|15;
      r=128+t*128; if(r<0) r=0; else if(r>255) r=255; c[a*3+b*pit+2]=(int)r|7;
    }
  }
}

unsigned char pic[256*256*3];

int main(void) {
  FILE *f;
  char fn[100];
  int a, l, m;
  for(a=0;a<256*256*3;++a) pic[a]=255;
  for(a=0;a<100;++a) {
    for(l=0;l<=4;++l) for(m=0;m<=l;++m) rend(pic+m*27*3+l*27*134*3, 26, 26, 134*3, l, m, PI*2/100*a);
//    for(l=0;l<=4;++l) for(m=0;m<=l;++m) rend(pic+m*51*3+l*51*256*3, 50, 50, 256*3, l, m, PI*2/100*a);
//    rend(pic+4*51*3+0*51*256*3, 50, 50, 256*3, 47, 1, PI*2/100*a);
//    rend(pic+4*51*3+0*51*256*3, 50, 50, 256*3, 47, 11, PI*2/100*a);
//    rend(pic, 256, 256, 256*3, 5, 2, PI*2/100*a);
    sprintf(fn, "sph%04d.ppm", a);
    f=fopen(fn, "wb");
    fprintf(f, "P6\n%d %d\n255\n", 134, 134);
    fwrite(pic, 134*134*3, 1, f);
    fclose(f);
  }
/*  int l, m, i;
  double r;
  for(l=0;l<=7;++l) for(m=0;m<=l;++m) {
    printf("%d,%d: %lf", l, m, r=cfac(l, m));
    for(i=((l+m+1)>>1)+1;i<=l;++i) printf(", %lf", r*=(double)(2*i-1)*(2*i)*(l-i+1)/((double)(2*i-l-m-1)*(2*i-l-m)*i));
    printf("\n");
  }*/
  return(0);
}