/* Calculate periodic orbits */ #include #include #include /* Macros for doing vector arithmatic */ #define setv(v, x1, y1) {v.x = x1; v.y = y1;} #define copyv(v, w) {v.x = w.x; v.y = w.y;} #define dotv(v1, v2) (v1.x*v2.x + v1.y*v2.y) #define hypotv(v) hypot(v.x, v.y) #define smulv(v1, a, v2) {v1.x=a*v2.x; v1.y=a*v2.y;} #define subv(v1, v2, v3) {v1.x=v2.x-v3.x; v1.y=v2.y-v3.y;} #define addv(v1, v2, v3) {v1.x=v2.x+v3.x; v1.y=v2.y+v3.y;} #define addvx(v1, v2, x1, y1) {v1.x=v2.x+x1; v1.y=v2.y+y1;} #define addtov(v, x1, y1) {v.x=v.x+x1; v.y=v.y+y1;} #define linear(v, a, va, b, vb) {v.x=a*va.x+b*vb.x; v.y=a*va.y+b*vb.y;} #define polar(v, r, theta) {v.x=r*cos(theta); v.y=r*sin(theta);} #define printv(v) printf("(%f, %f)", v.x, v.y) #define SMALL 1.0e-15 #define ND 3 #define PI 3.14159265358979323846 typedef struct _Vector { double x; double y; } Vector; char *progname, *orbit_str; main(argc, argv) /* driver */ int argc; char *argv[]; { int N, *orbit, old, i; void periodic(); double a=1.0/6.0; progname = argv[0]; while (argc>1 && argv[1][0]=='-') { switch (argv[1][1]) { case 'a': a = atof(argv[2]); break; case 'R': a = 1.0/atof(argv[2]); break; default: fprintf(stderr, "usage: %s [-a x.xx|-R x.xx] orbit_string\n", progname); fprintf(stderr, " e.g. > orbit -R 6 123\n"); exit(0); break; } argc-=2; argv+=2; } if (argc!=2) { fprintf(stderr, "usage: %s [-a x.xx|-R x.xx] orbit_string\n", progname); fprintf(stderr, " e.g. > orbit -R 6 123\n"); exit(0); } orbit_str = argv[1]; /* Convert orbit_string to list of disks and check for illegal orbit */ N = strlen(orbit_str); if (N<2) { fprintf(stderr, "error: %s is an illegal orbit\n", argv[1]); exit(0); } orbit = (int *) calloc(N, sizeof(int)); old = orbit_str[N-1] - '1'; for(i=0; i2 || orbit[i]==old) { fprintf(stderr, "error: %s is an illegal orbit\n", orbit_str); exit(0); } old = orbit[i]; } periodic(orbit, N, a); } /* calculate the periodic orbit */ void periodic(orbit, N, a) int *orbit; /* list of disk bounces */ int N; /* length of orbit */ double a; /* size of disk */ { Vector tmp_v, tmp_v1, A, B, C, disk[3], start[3], *Pos; int i, j, old, next; double len, old_len, tmp, tau, r, J[4], lam; double tmp1, tmp2, tmp3, tmp4, old_Bx; double l=0.5/sqrt(3.0), ia, asq; /* Set things up */ Pos = (Vector *) calloc(N+1, sizeof(Vector)); ia = 1.0/a; asq = a*a; setv(disk[0], -0.5, -l); setv(disk[1], 0.5, -l); setv(disk[2], 0, 2.0*l); addvx(start[0], disk[0], a, 0); addvx(start[1], disk[1], -a, 0); addvx(start[2], disk[2], 0, -a); /* starting positions */ for(i=0; iSMALL); addv(Pos[i], B, disk[orbit[i]]); subv(tmp_v, A, B); len += hypotv(tmp_v); subv(tmp_v, C, B); len += hypotv(tmp_v); } } while(fabs(old_len-len) > SMALL); printf("Orbit = %s, a = %f\n", orbit_str, a); for(i=0; i0.0) lam = 0.5*(tmp+sqrt(tmp*tmp-4.0)); else lam = 0.5*(tmp-sqrt(tmp*tmp-4.0)); /* Note that we have not remove the symmetry factor */ printf("Time = %f Stability = %f\n", ia*len, lam); }