Same program, different function calls:

Here's the code:

#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");

}