/*========================================================================= //========================================================================= ##### #### # # ##### ###### ##### ## #### ##### # # # # # # # # # # # # # # ##### # # ##### # # # # # # # # # # # ##### ###### # # # # # # # # # # # # # # # # #### ##### ##### # # # # # #### # Stephan LeBohec University of Utah Department of Physics 115 South 1400 East Salt-Lake-City, UT 84112-0830 Tel:801-587-9923 http://www.physics.utah.edu/~lebohec/ //========================================================================= //=========================================================================*/ #include #include #include #include "quaternions.c" #define PI 3.141592653589793 #define THREEDOUTPUT 0 /*Set to one to generate x,y,z output. Otherwise generates angle output from view point defined in function output()*/ struct generator { int nbseg; //Number of segments int nbsub; //Number of subdivision double *a; struct quaternion *q; }; int output(double vd, int i,double x,double y, double z); int read_generator(struct generator *g); int orthonormal(struct quaternion *dirx, struct quaternion *diry, struct quaternion *dirz); //========================================================================= //========================================================================= int main(int nbarg, char* args[]) { //Deal with the command line argument which specifies the number of iterations int nbiter; if(nbarg==2) nbiter=atoi(args[1]); else { fprintf(stderr,"You must provide the number of iterations as an argument\n"); exit(0); } //Starting point as a quaternion struct quaternion pt; pt=qset(0,0,-1.2,-1); /*Basis of vectors attached to the point. The fractal will start along direction dirx. */ struct quaternion dirx; dirx=qset(0,0,1,0); //Start along the y direction struct quaternion diry; diry=qset(0,-1,0,0); // struct quaternion dirz; //Construct an orthonormal and direct basis of vectors (dirx,diry,dirz) orthonormal(&dirx,&diry,&dirz); //Sets the distance of the viewing point along the x axis. double vdx=5; //Reads the generator struct generator g; read_generator(&g); double resol=pow(1.0/g.nbsub,nbiter); //Length of segments int nbrseg=(int)floor(pow(g.nbseg,nbiter)); //Number of segments /*Initialise the array of integers used to keep track of the curvilinear abcissa in base nbseg*/ int *sint; sint=new int[nbiter]; for(int i=0;i=0;k--) if(sint[k]!=0) { rotin=g.q[sint[k]%g.nbseg]; angle=g.a[sint[k]%g.nbseg]; }; struct quaternion rot; /*Computes the vector axis of rotation in the global coordinate system*/ rot=qadd(qscale(rotin.b,dirx),qscale(rotin.c,diry)); rot=qadd(rot,qscale(rotin.d,dirz)); rot=qrotset(angle,rot.b,rot.c,rot.d); //Apply the rotation to the vector of the local basis of vectors dirx=qrot(rot,dirx); diry=qrot(rot,diry); dirz=qrot(rot,dirz); } return 0; } //========================================================================= //========================================================================= int read_generator(struct generator *g) { /*The generator consists of a sequence of directions and angle of rotation around these directions. In the construction of the fractal curve the displacements are done along the x direction of the coordinate system attached to the point traveling along the curve. For fractal curves involving 90 degree angles only, the generator can be constructed from the following correspondance: Left turn: 90 0 0 1 Right turn: -90 0 0 1 Up turn: -90 0 1 0 Down turn: 90 0 1 0 The first rotation entry should be something that does nothing, like: 0 1 0 0 The first two entries should be the number of subdivision made in each segment when going from one iteration to the next. The second entry should be the number of segments in the generator. */ fscanf(stdin,"%d",&(g->nbsub)); fscanf(stdin,"%d",&(g->nbseg)); g->q=new struct quaternion [g->nbseg]; g->a=new double [g->nbseg]; for(int i=0;inbseg;i++) { double a,x,y,z; fscanf(stdin,"%lf%lf%lf%lf",&a,&x,&y,&z); g->q[i]=qset(0,x,y,z); g->a[i]=a; } return 0; } //========================================================================= //========================================================================= int orthonormal(struct quaternion *dirx, struct quaternion *diry, struct quaternion *dirz) { /*This function construct an orthonormal basis of vectors (dirx,diry,dirz). dirx keeps the same direction. dirz is taken perpendicular to dirx and diry. diry is then set perpendicular to dirx and dirz*/ //Ensures the three input quaternions are pure imaginary dirx->a=0; diry->a=0; dirz->a=0; //Normalize dirx *dirx=qunit(*dirx); //Construct dirz as the cross product dirx X diry *dirz=qcrossp(*dirx,*diry); //Normalize dirz *dirz=qunit(*dirz); //Construct diry (no need for normalization since |dirx|=|dirz|=1) *diry=qcrossp(*dirz,*dirx); return 0; } //========================================================================= //========================================================================= int output(double vd,int i,double x,double y, double z) { /*This function prints out the point x,y,z or its spherical coordinates from a viewing point on the x axis in -vd */ fprintf(stdout,"%d ",i); if(THREEDOUTPUT) fprintf(stdout,"%lf %lf %lf\n",x,y,z); else { /*The view point is on negative x axis at distance vd from origin. We calculate the angle theta between the x axis and the direction of the point and the angle phi corresponding to the position of the point around the x axis. View point in x=-vd,y=0,z=0*/ double phi=atan2(z,-y); //phi=0 for points with z=0 and y<0 double dx=vd+x; double theta; if(dx!=0) theta=atan(sqrt(y*y+z*z)/dx); else theta=0.5*PI; double r=sqrt(dx*dx+y*y+z*z); fprintf(stdout,"%lf %lf %lf\n",r,theta,phi); } return 0; } //========================================================================= //=========================================================================