#include <stdio.h>
#include <stdlib.h>
#include <strings.h>
#include <math.h>
#define L 400
int cnt[L][L];
#define SCALE 0.7
#define CMAX (9)
void dump()
{
FILE *fp=fopen("a.pgm","w");
fprintf(fp, "P5\n%d %d\n255\n",L,L);
int ix,iy;
for(ix=0;ix<L;ix++){
for(iy=0;iy<L;iy++){
fputc( (cnt[ix][iy]&255),fp);
}}
fclose(fp);
}
#define ISQ3 0.577350269189626
double n4vec[4][3]={
{ ISQ3, ISQ3, ISQ3},
{-ISQ3, -ISQ3, ISQ3},
{-ISQ3, ISQ3, -ISQ3},
{ ISQ3,-ISQ3, -ISQ3}
};
double nvec[6][3]={
{1,0,0},
{0,1,0},
{0,0,1},
{-1,0,0},
{0,-1,0},
{0,0,-1}
};
double n3vec[3][3]={
{1,0,0},{0,1,0},{0,0,1}};
void conv(double r[3], int ax)
{
double *v = nvec[ax];
double rv = r[0]*v[0]+r[1]*v[1]+r[2]*v[2];
double rr = r[0]*r[0] + r[1]*r[1] + r[2]*r[2];
r[0] = 2*rv * r[0] - rr *v[0] ;
r[1] = 2*rv * r[1] - rr *v[1] ;
r[2] = 2*rv * r[2] - rr *v[2] ;
}
double rfunc(double c0[3], int c)
{
int cc;
double rr[3]={c0[0],c0[1],c0[2]};
for(cc=0;cc<c;cc++)
{
conv(rr, 0);
conv(rr, 1);
conv(rr, 2);
rr[0]+=c0[0]; rr[1]+=c0[1]; rr[2]+=c0[2];
}
return rr[0]*rr[0]+rr[1]*rr[1]+rr[2]*rr[2] - 4.0;
}
double bisec_z(int lvl, double p0[3], double dz, double eps, double zmax)
{
while(p0[2]> -zmax)
{
double r2 = rfunc(p0,lvl);
if (r2<0.0) break;
p0[2] -= dz;
}
if (p0[2]< -zmax)return -999.0;
p0[2] += dz;
dz/=2;
while(dz>eps)
{
p0[2] -= dz;
double r2 = rfunc(p0,lvl);
if (r2<0.0) p0[2]+=dz;
dz/=2;
}
return p0[2];
}
int main(int argc, char **argv)
{
int ix,iy;
double c0[3],rr[3];
for(ix=0;ix<L;ix++){
printf("%d/%d\n",ix,L);
for(iy=0;iy<L;iy++){
c0[0]=SCALE*(2.0*ix/L -1.0);
c0[1]=SCALE*(2.0*iy/L -1.0);
double rr=c0[0]*c0[0]+c0[1]*c0[1];
if(rr>=4.0){ cnt[ix][iy]=0;continue;}
double zmax=sqrt(4-c0[0]*c0[0]-c0[1]*c0[1]);
double dz=0.1;
double zz=zmax+0.01;
c0[2]=zmax;
int c;
int skip=0;
for(c=0;c<=CMAX;c++)
{
double z0=bisec_z(c, c0, dz, 1e-7, zmax);
if (z0<-888.0)
{
skip=1;
break;
}
}
if(skip) { cnt[ix][iy]=0;continue;}
#define EPS 1e-9
double r2=rfunc(c0,CMAX);
double cz[3]={c0[0],c0[1],c0[2]+EPS};
double r2z=rfunc(cz,CMAX)-r2;
double cx[3]={c0[0]+EPS,c0[1],c0[2]};
double r2x=rfunc(cx,CMAX)-r2;
double cy[3]={c0[0],c0[1]+EPS,c0[2]};
double r2y=rfunc(cy,CMAX)-r2;
double nn=sqrt(r2z*r2z+r2x*r2x+r2y*r2y);
cnt[ix][iy]=(int)((0.2+0.8*fabs( r2z/nn))*255);
}}
dump();
}