ความย่อ

การอนุญาตใช้สิทธิ

Source

//GPL #include <stdio.h> #include <math.h> #include <vector> using std::vector; const char *theader = "//Picture *** Use flashiness=1 !!! ***\n//\n// +w1024 +h1024 +a0.3 +am2\n// +w512 +h512 +a0.3 +am2\n//\n//Movie *** Use flashiness=0.25 !!! ***\n//\n// +kc +kff120 +w256 +h256 +a0.3 +am2\n// +kc +kff60 +w256 +h256 +a0.3 +am2\n//\"Fast\" preview\n// +w128 +h128\n#declare notwireframe=1;\n#declare withreflection=0;\n#declare flashiness=1; //Still pictures use 1, animated should probably be about 0.25.\n\n#declare rotation=seed(%d);\n\n#declare rot1=rand(rotation)*pi*2;\n#declare rot2=acos(1-2*rand(rotation));\n#declare rot3=(rand(rotation)+clock)*pi*2;\n#macro dorot()\n rotate rot1*180/pi*y\n rotate rot2*180/pi*x\n rotate rot3*180/pi*y\n#end\n\n"; const char *tline = "object {\n cylinder { <%lf,%lf,%lf>, <%lf,%lf,%lf>, .01 dorot() }\n pigment { colour <.3,.3,.3> }\n finish { ambient 0 diffuse 1 phong 1 }\n}\n\n"; const char *tvertex = "object {\n sphere { <%lf,%lf,%lf>, .01 dorot() }\n pigment { colour <.3,.3,.3> }\n finish { ambient 0 diffuse 1 phong 1 }\n}\n\n"; const char *tstartmesh = "object {\n mesh {\n"; const char *ttriangle = " triangle {\n <%lf,%lf,%lf>, <%lf,%lf,%lf>, <%lf,%lf,%lf>\n }\n"; const char *tendmesh = " //sphere { <0,0,0>, 1 }\n //sphere { <0,0,0>, ld+.01 inverse }\n dorot()\n }\n pigment { colour rgbt <.8,.8,.8,.4> }\n finish { ambient 0 diffuse 1 phong flashiness #if(withreflection) reflection { .2 } #end }\n //interior { ior 1.5 }\n photons {\n target on\n refraction on\n reflection on\n collect on\n }\n}\n\n"; const char *tfooter = "// CCC Y Y PP\n// C Y Y P P\n// C Y PP\n// C Y P\n// CCC Y P\n\n#local a=0;\n#while(a<11.0001)\n light_source { <4*sin(a*pi*2/11), 5*cos(a*pi*6/11), -4*cos(a*pi*2/11)> colour (1+<sin(a*pi*2/11),sin(a*pi*2/11+pi*2/3),sin(a*pi*2/11+pi*4/3)>)*2/11 }\n #local a=a+1;\n#end\n\nbackground { color <1,1,1> }\n\ncamera {\n perspective\n location <0,0,0>\n direction <0,0,1>\n right x/2\n up y/2\n sky <0,1,0>\n location <0,0,-4.8>\n look_at <0,0,0>\n}\n\nglobal_settings {\n max_trace_level 40\n photons {\n count 200000\n autostop 0\n }\n}\n"; #define PHI ((1+sqrt(5))/2) #define PI (3.14159265358979323846264338327) #define SQ2 (sqrt(2)) #define SQ3 (sqrt(3)) bool eq(double a, double b) { return a+0.00001>=b&&b+0.00001>=a; } bool eqt(double a1, double a2, double a3, double b1, double b2, double b3) { //printf("Tri: {%lf, %lf, %lf}, {%lf, %lf, %lf}\n", a1, a2, a3, b1, b2, b3); return eq(a1, b1)? eq(a2, b2)? eq(a3, b3):eq(a2, b3)&&eq(a3, b2):eq(a1, b2)? eq(a2, b3)? eq(a3, b1):eq(a2, b1)&&eq(a3, b3):eq(a1, b3)&&(eq(a2, b1)? eq(a3, b2):eq(a2, b3)&&eq(a3, b2)); } class vec { public: double x, y, z; vec() : x(0), y(0), z(0) {} vec(double nx, double ny, double nz) : x(nx), y(ny), z(nz) {} vec operator + (vec o) { return vec(x+o.x, y+o.y, z+o.z); } vec operator - (vec o) { return vec(x-o.x, y-o.y, z-o.z); } double operator * (vec o) { return x*o.x+y*o.y+z*o.z; } vec operator * (double o) { return vec(x*o, y*o, z*o); } vec operator ^ (vec o) { return vec(y*o.z-z*o.y, z*o.x-x*o.z, x*o.y-y*o.x); } double norm() { return sqrt(x*x+y*y+z*z); } }; class vec2 { public: double x, y; vec2() {} vec2(double nx, double ny) : x(nx), y(ny) {} vec2 operator + (vec2 o) { return vec2(x+o.x, y+o.y); } vec2 operator - (vec2 o) { return vec2(x-o.x, y-o.y); } double operator * (vec2 o) { return x*o.x+y*o.y; } vec2 operator * (double o) { return vec2(x*o, y*o); } vec2 operator ~ () { return vec2(y, -x); } double norm() { return sqrt(x*x+y*y); } }; vector<vec> cyclicperm(vector<vec> v) { vector<vec> r; vector<vec>::iterator i; for(i = v.begin(); i!=v.end(); ++i) { r.push_back(*i); r.push_back(vec(i->y, i->z, i->x)); r.push_back(vec(i->z, i->x, i->y)); } return r; } vector<vec> altperm(vector<vec> v) { vector<vec> r; vector<vec>::iterator i; for(i = v.begin(); i!=v.end(); ++i) { r.push_back(*i); r.push_back(vec(i->x, i->z, i->y)); } return r; } vector<vec> signperm(vector<vec> v) { vector<vec> r; vector<vec>::iterator i; for( i = v.begin(); i!=v.end(); ++i ) { int j; for(j = 0; j<8; ++j) if(((j&1)||i->x)&&((j&2)||i->y)&&((j&4)||i->z)) r.push_back(vec(j&1? i->x:-i->x, j&2? i->y:-i->y, j&4? i->z:-i->z)); } return r; } vector<vec> mvvec(double x, double y, double z) { vector<vec> v; v.push_back(vec(x, y, z)); return v; } vector<vec> mvvec(vec q) { vector<vec> v; v.push_back(q); return v; } vector<vec> concat(const vector<vec> a, const vector<vec> b) { vector<vec> r; r = a; r.insert(r.end(), b.begin(), b.end()); return r; } void printvvec(FILE *f, vector<vec> v) { vector<vec>::iterator i; for(i = v.begin(); i!=v.end(); ++i) fprintf(f, tvertex, i->x, i->y, i->z); } void printvveclines(FILE *f, vector<vec> v, double len) { vector<vec>::iterator i, j; len *= len; for(i = v.begin(); i!=v.end(); ++i) for(j = i+1; j!=v.end(); ++j) if(eq((*i-*j)*(*i-*j), len)) fprintf(f, tline, i->x, i->y, i->z, j->x, j->y, j->z); } void printvveclines(FILE *f, vector<vec> v) { vector<vec>::iterator i; for(i = v.begin(); i!=v.end(); i += 2) fprintf(f, tline, i->x, i->y, i->z, (i+1)->x, (i+1)->y, (i+1)->z); } void printvvecdottedlines(FILE *f, vector<vec> v) { vector<vec>::iterator i; int n, m; double s; for(i = v.begin(); i!=v.end(); i += 2) // for(i = v.begin(); i!=v.begin()+12; i += 2) { s = (*i-*(i+1)).norm(); m = (int)(s/0.04+.5); s = 1./(double)m; for(n = 1; n<m; ++n) { vec c = *i+(*(i+1)-*i)*(s*n); fprintf(f, tvertex, c.x, c.y, c.z); } } } void printvvectriangles(FILE *f, vector<vec> v, double len1, double len2, double len3) { vector<vec>::iterator i, j, k; len1 *= len1; len2 *= len2; len3 *= len3; for(i = v.begin(); i!=v.end(); ++i) for(j = i+1; j!=v.end(); ++j) for(k = j+1; k!=v.end(); ++k) if(eqt((*i-*j)*(*i-*j), (*j-*k)*(*j-*k), (*k-*i)*(*k-*i), len1, len2, len3)) fprintf(f, ttriangle, i->x, i->y, i->z, j->x, j->y, j->z, k->x, k->y, k->z); } void printvvectriangles(FILE *f, vector<vec> v) { vector<vec>::iterator i; for(i = v.begin(); i!=v.end(); i += 3) //i = v.begin(); fprintf(f, ttriangle, i->x, i->y, i->z, (i+1)->x, (i+1)->y, (i+1)->z, (i+2)->x, (i+2)->y, (i+2)->z); /*i += 3; fprintf(f, ttriangle, i->x, i->y, i->z, (i+1)->x, (i+1)->y, (i+1)->z, (i+2)->x, (i+2)->y, (i+2)->z); i += 3; fprintf(f, ttriangle, i->x, i->y, i->z, (i+1)->x, (i+1)->y, (i+1)->z, (i+2)->x, (i+2)->y, (i+2)->z); i += 3; fprintf(f, ttriangle, i->x, i->y, i->z, (i+1)->x, (i+1)->y, (i+1)->z, (i+2)->x, (i+2)->y, (i+2)->z); i += 3; fprintf(f, ttriangle, i->x, i->y, i->z, (i+1)->x, (i+1)->y, (i+1)->z, (i+2)->x, (i+2)->y, (i+2)->z); i += 3; fprintf(f, ttriangle, i->x, i->y, i->z, (i+1)->x, (i+1)->y, (i+1)->z, (i+2)->x, (i+2)->y, (i+2)->z); */} void SmallStellatedDodecahedron() { vector<vec> v; v = cyclicperm(signperm(mvvec(vec(0, PHI, 1)*(1/sqrt(PHI+2))))); FILE *f; f = fopen("SmallStellatedDodecahedron.pov", "wb"); fprintf(f, theader, 22491); printvvec(f, v); printvveclines(f, v, 2*PHI*(1/sqrt(PHI+2))); fprintf(f, tstartmesh); v = concat(v, cyclicperm(signperm(mvvec(vec(0, 2-PHI, 1)*(1/sqrt(PHI+2)))))); v = concat(v, signperm(mvvec(vec(PHI-1, PHI-1, PHI-1)*(1/sqrt(PHI+2))))); printvvectriangles(f, v, (2*PHI-2)*(1/sqrt(PHI+2)), (2*PHI-2)*(1/sqrt(PHI+2)), (4-2*PHI)*(1/sqrt(PHI+2))); fprintf(f, tendmesh); fprintf(f, tfooter); fclose(f); } void GreatStellatedDodecahedron() { vector<vec> v; v = concat(signperm(mvvec(vec(1, 1, 1)*(1/SQ3))), cyclicperm(signperm(mvvec(vec(0, PHI, 1/PHI)*(1/SQ3))))); FILE *f; f = fopen("GreatStellatedDodecahedron.pov", "wb"); fprintf(f, theader, 7409);//7412); printvvec(f, v); printvveclines(f, v, 2*PHI*(1/SQ3)); fprintf(f, tstartmesh); v = concat(v, cyclicperm(signperm(mvvec(vec(0, 2-PHI, PHI-1)*(1/SQ3))))); printvvectriangles(f, v, (2*PHI-2)*(1/SQ3), (2*PHI-2)*(1/SQ3), (4-2*PHI)*(1/SQ3)); fprintf(f, tendmesh); fprintf(f, tfooter); fclose(f); } void GreatDodecahedron() { vector<vec> v; v = cyclicperm(signperm(mvvec(vec(0, PHI, 1)*(1/sqrt(PHI+2))))); FILE *f; f = fopen("GreatDodecahedron.pov", "wb"); fprintf(f, theader, 11404); printvveclines(f, v, 2*(1/sqrt(PHI+2))); v = concat(v, concat(signperm(mvvec(vec(PHI-1, PHI-1, PHI-1)*(1/sqrt(PHI+2)))), cyclicperm(signperm(mvvec(vec(0, 2-PHI, 1)*(1/sqrt(PHI+2))))))); printvvec(f, v); fprintf(f, tstartmesh); printvvectriangles(f, v, (2*PHI-2)*(1/sqrt(PHI+2)), (2*PHI-2)*(1/sqrt(PHI+2)), (2)*(1/sqrt(PHI+2))); fprintf(f, tendmesh); fprintf(f, tfooter); fclose(f); } vector<vec> IcosaParse(const char *vs) { vector<vec> v, p; v = cyclicperm(signperm(mvvec(vec(0, PHI, 1)))); vec av; vector<vec>::iterator i, j, k; int q; static const vec2 rats[9] = {vec2(1, 0), vec2(PHI-1, 2-PHI), vec2(2-PHI, PHI-1), vec2(0, 1), vec2(0, PHI-1), vec2(0, 2-PHI), vec2(0, 0), vec2(2-PHI, 0), vec2(PHI-1, 0)}; for(i = v.begin(); i!=v.end(); ++i) for(j = v.begin(); j!=v.end(); ++j) for(k = v.begin(); k!=v.end(); ++k) if(eqt((*i-*j).norm(), (*j-*k).norm(), (*k-*i).norm(), 2, 2, 2)&&(*i^*j)**k>0) { vec t3 = *i*PHI*PHI+*j*PHI*PHI-*k*PHI*PHI*PHI, t1 = *j*PHI*PHI+*k*PHI*PHI-*i*PHI*PHI*PHI, t2 = *k*PHI*PHI+*i*PHI*PHI-*j*PHI*PHI*PHI; for(q = 0; vs[q]; ) { if(vs[q]<48) break; if(vs[q+1]<48) { p = concat(p, mvvec(t3+(t1-t3)*rats[vs[q]-'0'].x+(t2-t3)*rats[vs[q]-'0'].y)); q += 2; continue; } if(vs[q+4]<48) { vec2 a = rats[vs[q]-'0'], b = rats[vs[q+1]-'0'], c = rats[vs[q+2]-'0'], d = rats[vs[q+3]-'0']; double idet = 1/((a-b).x*(d-c).y-(a-b).y*(d-c).x); //fprintf(stderr, "%lf, %lf %lf, %lf %lf\n", (a-b).x, (d-c).x, (a-b).y, (d-c).y, idet); vec2 e = vec2(vec2((d-c).y, (d-c).x*-1)*(d-b), vec2((a-b).y*-1, (a-b).x)*(d-b))*idet; vec2 r = (a-b)*e.x+b; //fprintf(stderr, "%lf, %lf %lf, %lf %lf\n", r.x, r.y, t1.x, t1.y, idet); //fprintf(stderr, "(a-b)={%lf, %lf}, x=%lf, b={%lf, %lf}, e={%lf, %lf}\n(c-d)={%lf, %lf}, y=%lf, d={%lf, %lf}, e={%lf, %lf}\n", //(a-b).x, (a-b).y, e.x, b.x, b.y, ((a-b)*e.x+b).x, ((a-b)*e.x+b).y, //(c-d).x, (c-d).y, e.y, d.x, d.y, ((c-d)*e.y+d).x, ((c-d)*e.y+d).y //); //fprintf(stderr, "%lf %lf\n", r.x, r.y); p = concat(p, mvvec(t3+(t1-t3)*r.x+(t2-t3)*r.y)); av = av+(t3+(t1-t3)*r.x+(t2-t3)*r.y); //p = concat(p, mvvec(vec())); q += 5; continue; } break; } } //printf("%lf %lf %lf\n", av.x, av.y, av.z); double r = 0; for(i = p.begin(); i!=p.end(); ++i) //i = p.begin(); if(r<i->norm()) r = i->norm(); for(i = p.begin(); i!=p.end(); ++i) *i = *i*(1/r); return p; } void StellatedIcosahedron(const char *fn, int rs, const char *vs, const char *ls, const char *dls, const char *ts) { vector<vec> v; FILE *f; f = fopen(fn, "wb"); fprintf(f, theader, rs); printvvec(f, IcosaParse(vs)); printvvecdottedlines(f, IcosaParse(dls)); printvveclines(f, IcosaParse(ls)); fprintf(f, tstartmesh); printvvectriangles(f, IcosaParse(ts)); fprintf(f, tendmesh); fprintf(f, tfooter); fclose(f); } int main() { SmallStellatedDodecahedron(); GreatStellatedDodecahedron(); GreatDodecahedron(); StellatedIcosahedron("GreatIcosahedron.pov", 31234, "0 1 2 0417 1428 2538 ", "0 3 ", "0 0417 0417 1 1 1428 1428 2 2 2538 2538 3 ", "0 1 0417 1 2 1428 2 3 2538 "); StellatedIcosahedron("CompoundOfFiveTetrahedra.pov", 22113, "2 2514 1427 2715 1528 ", "2 5 ", "2 2 2514 1427 1427 2715 2715 1528 ", "2 2514 1427 2 2715 1528 "); return 0; }