#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
double pi = 3.1415926535;
struct point {
double x;
double y;
double z;
};
struct quaternion {
double a;
double b;
double c;
double d;
};
quaternion div (quaternion q1, double r){
quaternion s =
{q1.a/r,q1.b/r,q1.c/r,q1.d/r};
return(s);
}
quaternion add(quaternion q1, quaternion q2){
quaternion s =
{q1.a+q2.a,
q1.b+q2.b,
q1.c+q2.c,
q1.d+q2.d};
return(s);
}
quaternion mult(quaternion q1, quaternion q2){
quaternion s =
{q1.a*q2.a-q1.b*q2.b-q1.c*q2.c-q1.d*q2.d,
q1.a*q2.b+q2.a*q1.b+q1.c*q2.d-q1.d*q2.c,
q1.a*q2.c-q2.d*q1.b+q1.c*q2.a+q1.d*q2.b,
q1.a*q2.d+q2.c*q1.b-q1.c*q2.b+q1.d*q2.a};
return(s);
}
quaternion inv(quaternion q){
double r = q.a*q.a+q.b*q.b+q.c*q.c+q.d*q.d;
// if (r<0.001) {r=0.001;}
quaternion s = {q.a/r,-q.b/r,-q.c/r,-q.d/r};
return(s);
}
double norm(quaternion q){
double r = q.a*q.a+q.b*q.b+q.c*q.c+q.d*q.d;
if (r > 0.0001)
return(sqrt(r));
else
return(1);
}
int iterate(point p,
point shift,
double r, double scale, quaternion q,
point shift0,
FILE *f1, int depth, int ctr) {
int DEPTH_MAX = 1000;
double minrad = 0.0025;//25;
point p1;
quaternion q1,q2,q3,q4;
double c5 = cos(((double) depth) * pi/8);
double s5 = sin(((double) depth) * pi/8);
q1 = {0.0,shift0.x,shift0.y,shift0.z}; // original branch
q2 = {0.0,shift.x,shift.y,shift.z}; // next branch
double cosine =
(shift.x*shift0.x+
shift.y*shift0.y+
shift.z*shift0.z)/(norm(q1)*norm(q2));
// cross-product
q1.a = 0.0;
q1.b = shift0.y*shift.z-shift0.z*shift.y;
q1.c = shift0.z*shift.x-shift0.x*shift.z;
q1.d = shift0.x*shift.y-shift0.y*shift.x;
q1 = div(q1,norm(q1));
double theta = acos(cosine);
//printf("theta = %f\n",theta);
if (fabs(theta)>-0.0001) {
double s = sin(theta*0.5);
q2 = {cos(theta*0.5),s*q1.b,s*q1.c,s*q1.d};
// rotate shift to reference frame
q1 = {0,shift.x,shift.y,shift.z};
q4 = mult(q,mult(q1,inv(q)));
} else {
q4 = q1;
}
point shift1 = {q4.b,q4.c,q4.d};
p1.x = p.x + shift1.x*scale;
p1.y = p.y + shift1.y*scale;
p1.z = p.z + shift1.z*scale;
double r1 = r*scale;
// rotate reference frame
q = mult(q,q2);
fprintf(f1,"%s %f %f %f %s %f %s%d\n",
"sphere center",
p1.x,p1.y,p1.z,
"rad",
r1, "t",1+(ctr % 7));
if ((depth < DEPTH_MAX)&&(r1>minrad)) {
double y1,z1;
quaternion q1 = q;
q3 = q;
shift = {0,1,0};
double factor = 0.3;
point s1;
q2 = {c5,0,s5,0};
s1 = {0.15,0.3,0.025};
double t = iterate(p1,s1,r,scale*0.9,q3,shift,f1,depth+1,ctr+1);
s1 = {-0.7,0,0};
q1 = {0.0,s1.x,s1.y,s1.z};
q4 = mult(q2,mult(q1,inv(q2)));
s1 = {q4.b,q4.c,q4.d};
double t1 = iterate(p1,s1,r,scale*factor,q3,shift,f1,depth+1,ctr);
s1 = {+0.7,0.0,0.0};
q1 = {0.0,s1.x,s1.y,s1.z};
q4 = mult(q2,mult(q1,inv(q2)));
s1 = {q4.b,q4.c,q4.d};
double t5 = iterate(p1,s1,r,scale*factor,q3,shift,f1,depth+1,ctr);
s1 = {0.0,0.0,-0.7};
q1 = {0.0,s1.x,s1.y,s1.z};
q4 = mult(q2,mult(q1,inv(q2)));
s1 = {q4.b,q4.c,q4.d};
double t6 = iterate(p1,s1,r,scale*factor,q3,shift,f1,depth+1,ctr);
s1 = {0.0,0,+0.7};
q1 = {0.0,s1.x,s1.y,s1.z};
q4 = mult(q2,mult(q1,inv(q2)));
s1 = {q4.b,q4.c,q4.d};
double t7 = iterate(p1,s1,r,scale*factor,q3,shift,f1, depth+1,ctr);
}
}
int main() {
FILE *f1;
f1 = fopen("temp.dat","w");
point P = {0,-0.35,0};
point S = {0,0.1,0};
double c = cos(pi/15);
double s = sin(pi/15);
quaternion q = {c,0,s,0};
quaternion q2 = {c,s,0,0};
// q = mult(q,q2);
iterate(P,S,0.2,1.0,q,{0,1,0},f1,0,3);
fclose(f1);
system("cat header.dat temp.dat > temp2.dat");
system("./tachyon temp2.dat");
system("gimp outfile.tga");
}