diff --git a/sub_part_billiard_phasespace.c b/sub_part_billiard_phasespace.c new file mode 100644 index 0000000..bda0376 --- /dev/null +++ b/sub_part_billiard_phasespace.c @@ -0,0 +1,6187 @@ +/* variant of sub_part_billiard.c with lambda as global variable */ + +#include "colormaps.c" + +#define DUMMY_ABSORBING -1000.0 /* dummy value of config[0] for absorbing circles */ +#define BOUNDARY_SHIFT 100000.0 /* shift of boundary parametrisation for circles in domain */ +#define DUMMY_SIDE_ABS -10000 /* dummy value of returned side for absorbing circles */ + +long int global_time = 0; /* counter to keep track of global time of simulation */ +int nparticles=NPART; + +double lambda, mu, penrose_ratio, scaling_factor; + +/*********************/ +/* some basic math */ +/*********************/ + + double vabs(double x) /* absolute value */ + { + double res; + + if (x<0.0) res = -x; + else res = x; + return(res); + } + + double module2(double x, double y) /* Euclidean norm */ + { + double m; + + m = sqrt(x*x + y*y); + return(m); + } + + double argument(double x, double y) + { + double alph; + + if (x!=0.0) + { + alph = atan(y/x); + if (x<0.0) + alph += PI; + } + else + { + alph = PID; + if (y<0.0) + alph = PI*1.5; + } +// if (alph < 0.0) alph += DPI; + return(alph); + } + + int polynome(double a, double b, double c, double r[2]) + { + double delta, rdelta; + int im = 1; + + delta = b*b - 4*a*c; + if (delta<0.0) + { + /* printf("ca deconne!");*/ + rdelta = 0.0; + im = 0; + } + else rdelta = sqrt(delta); + + r[0] = (-b + rdelta)/(2.0*a); + r[1] = (-b - rdelta)/(2.0*a); + + return(im); + } + + double ipow(double x, int n) + { + double y; + int i; + + y = x; + for (i=1; i= 0; i--) + { +// if (TIFFWriteScanline(file, p, height - i - 1, 0) < 0) + if (TIFFWriteScanline(file, p, i, 0) < 0) + { + free(image); + TIFFClose(file); + return 1; + } + p += width * sizeof(GLubyte) * 3; + } +// free(image); + TIFFClose(file); + if (SAVE_MEMORY) free(image); + return 0; +} + + +void init() /* initialisation of window */ +{ + glLineWidth(3); + + glClearColor(0.0, 0.0, 0.0, 1.0); + glClear(GL_COLOR_BUFFER_BIT); + + glOrtho(XMIN, XMAX, YMIN, YMAX , -1.0, 1.0); +} + + + +void rgb_color_scheme(int i, double rgb[3]) /* color scheme */ +{ + double hue, y, r; + + hue = (double)(COLORSHIFT + i*360/NCOLORS); + r = 0.9; + + while (hue < 0.0) hue += 360.0; + while (hue >= 360.0) hue -= 360.0; + + hsl_to_rgb(hue, r, 0.5, rgb); + /* saturation = r, luminosity = 0.5 */ +} + + +void rgb_color_scheme_lum(int i, double lum, double rgb[3]) /* color scheme */ +{ + double hue, y, r; + + hue = (double)(COLORSHIFT + i*360/NCOLORS); + r = 0.9; + + while (hue < 0.0) hue += 360.0; + while (hue >= 360.0) hue -= 360.0; + + hsl_to_rgb(hue, r, lum, rgb); + /* saturation = r */ +} + +void blank() +{ + double rgb[3]; + + if (COLOR_OUTSIDE) + { + hsl_to_rgb(OUTER_COLOR, 0.9, 0.15, rgb); + glClearColor(rgb[0], rgb[1], rgb[2], 1.0); + } + else if (BLACK) glClearColor(0.0, 0.0, 0.0, 1.0); + else glClearColor(1.0, 1.0, 1.0, 1.0); + glClear(GL_COLOR_BUFFER_BIT); +} + + +void save_frame() +{ + static int counter = 0; + char *name="part.", n2[100]; + char format[6]=".%05i"; + + counter++; +// printf (" p2 counter = %d \n",counter); + strcpy(n2, name); + sprintf(strstr(n2,"."), format, counter); + strcat(n2, ".tif"); + printf(" saving frame %s \n",n2); + writetiff(n2, "Billiard in an ellipse", 0, 0, + WINWIDTH, WINHEIGHT, COMPRESSION_LZW); + +} + +void save_frame_counter(int i) +{ + char *name="part.", n2[100]; + char format[6]=".%05i"; + + strcpy(n2, name); + sprintf(strstr(n2,"."), format, i); + strcat(n2, ".tif"); + printf(" saving frame %s \n",n2); + writetiff(n2, "Billiard in an ellipse", 0, 0, + WINWIDTH, WINHEIGHT, COMPRESSION_LZW); + +} + + +void write_text_fixedwidth( double x, double y, char *st) +{ + int l, i; + + l=strlen( st ); // see how many characters are in text string. + glRasterPos2d( x, y); // location to start printing text + for( i=0; i < l; i++) // loop until i is greater then l + { +// glutBitmapCharacter(GLUT_BITMAP_TIMES_ROMAN_24, st[i]); // Print a character on the screen +// glutBitmapCharacter(GLUT_BITMAP_8_BY_13, st[i]); // Print a character on the screen + glutBitmapCharacter(GLUT_BITMAP_9_BY_15, st[i]); // Print a character on the screen + } +} + +void write_text( double x, double y, char *st) +{ + int l,i; + + l=strlen( st ); // see how many characters are in text string. + glRasterPos2d( x, y); // location to start printing text + for( i=0; i < l; i++) // loop until i is greater then l + { + glutBitmapCharacter(GLUT_BITMAP_TIMES_ROMAN_24, st[i]); // Print a character on the screen +// glutBitmapCharacter(GLUT_BITMAP_8_BY_13, st[i]); // Print a character on the screen + } +} + +void erase_area(double x, double y, double dx, double dy, double rgb[3]) +{ + glColor3f(rgb[0], rgb[1], rgb[2]); + glBegin(GL_QUADS); + glVertex2d(x - dx, y - dy); + glVertex2d(x + dx, y - dy); + glVertex2d(x + dx, y + dy); + glVertex2d(x - dx, y + dy); + glEnd(); +} + +void erase_rectangle(double x1, double y1, double x2, double y2, double rgb[3]) +{ + glColor3f(rgb[0], rgb[1], rgb[2]); + glBegin(GL_QUADS); + glVertex2d(x1, y1); + glVertex2d(x2, y1); + glVertex2d(x2, y2); + glVertex2d(x1, y2); + glEnd(); +} + +void erase_rectangle_outside(double h, double s, double l) +{ + double rgb[3], dx; + int k; + + dx = 0.5*(XMAX - lambda); + hsl_to_rgb(h, s, l, rgb); + erase_rectangle(XMIN, YMIN, XMAX, -1.0, rgb); + erase_rectangle(XMIN, 1.0, XMAX, YMAX, rgb); + erase_rectangle(XMIN, YMIN, -lambda, YMAX, rgb); + erase_rectangle(lambda, YMIN, XMAX, YMAX, rgb); +// erase_area(0.0, 1.1, 2.0, 0.1, rgb); +// erase_area(0.0, -1.1, 2.0, 0.1, rgb); +// erase_area(lambda + dx, 0.0, dx, 2.0, rgb); +// erase_area(-lambda - dx, 0.0, dx, 2.0, rgb); + + glColor3f(1.0, 1.0, 1.0); + glLineWidth(BILLIARD_WIDTH); + glBegin(GL_LINE_LOOP); + glVertex2d(lambda, -1.0); + glVertex2d(lambda, 1.0); + glVertex2d(-lambda, 1.0); + glVertex2d(-lambda, -1.0); + glEnd(); +} + +void draw_line(double x1, double y1, double x2, double y2) +{ + glBegin(GL_LINE_STRIP); + glVertex2d(x1, y1); + glVertex2d(x2, y2); + glEnd(); +} + +void draw_circle(double x, double y, double r, int nseg) +{ + int i; + double phi, dphi, x1, y1; + + dphi = DPI/(double)nseg; + + glEnable(GL_LINE_SMOOTH); + glBegin(GL_LINE_LOOP); + for (i=0; i<=nseg; i++) + { + phi = (double)i*dphi; + x1 = x + r*cos(phi); + y1 = y + r*sin(phi); + glVertex2d(x1, y1); + } + glEnd (); +} + +void draw_colored_circle(double x, double y, double r, int nseg, double rgb[3]) +{ + int i, ij[2]; + double pos[2], alpha, dalpha; + + dalpha = DPI/(double)nseg; + +// glLineWidth(2); + + glColor3f(rgb[0], rgb[1], rgb[2]); + glBegin(GL_TRIANGLE_FAN); + glVertex2d(x,y); + for (i=0; i<=nseg; i++) + { + alpha = (double)i*dalpha; + glVertex2d(x + r*cos(alpha), y + r*sin(alpha)); + } + + glEnd(); +} + + +void draw_initial_condition_circle(double x, double y, double r, int color) +/* draws a colored circle to mark initial condition */ +{ + double rgb[3]; + + rgb_color_scheme(color, rgb); + draw_colored_circle(x, y, r, NSEG, rgb); + rgb[0] = 0.0; rgb[1] = 0.0; rgb[2] = 0.0; + glColor3f(rgb[0], rgb[1], rgb[2]); + glLineWidth(4); + draw_circle(x, y, r, NSEG); +} + +int in_circle(double x, double y, double r) +/* test whether (x,y) is in circle of radius r */ +{ + return(x*x + y*y < r*r); +} + +int out_circle(double x, double y, double r) +/* test whether (x,y) is in circle of radius r */ +{ + return(x*x + y*y > r*r); +} + +int in_polygon(double x, double y, double r, int npoly, double apoly) +/* test whether (x,y) is in regular polygon of npoly sides inscribed in circle of radious r, turned by apoly Pi/2 */ +{ + int condition = 1, k; + double omega, cw, angle; + + omega = DPI/((double)npoly); + cw = cos(omega*0.5); + for (k=0; k 0.0) x2 = cos(omega*0.5) + sqrt(lambda*lambda - sin(omega*0.5)*sin(omega*0.5)); + else x2 = cos(omega*0.5) - sqrt(lambda*lambda - sin(omega*0.5)*sin(omega*0.5)); + + if (PAINT_INT) + { + if (BLACK) glColor3f(1.0, 1.0, 1.0); + else glColor3f(0.0, 0.0, 0.0); + + glBegin(GL_TRIANGLE_FAN); + glVertex2d(0.0, 0.0); + for (i=0; i<=NPOLY; i++) + { + for (j=0; j 0.0) x2 = cos(omega*0.5) + sqrt(lambda*lambda - sin(omega*0.5)*sin(omega*0.5)); + else x2 = cos(omega*0.5) - sqrt(lambda*lambda - sin(omega*0.5)*sin(omega*0.5)); + glBegin(GL_LINE_STRIP); + for (i=0; i<=NPOLY; i++) + { + for (j=0; j= 1) + { + rgb_color_scheme_lum(circles[k].color, 0.85, rgb1); + circles[k].new--; + } + else rgb_color_scheme(circles[k].color, rgb1); + draw_colored_circle(circles[k].xc, circles[k].yc, circles[k].radius, NSEG, rgb1); + } + } + init_billiard_color(); + for (k=0; k= 1) + { + rgb_color_scheme_lum(circles[k].color, 0.85, rgb1); + circles[k].new--; + } + else rgb_color_scheme(circles[k].color, rgb1); + draw_colored_circle(circles[k].xc, circles[k].yc, circles[k].radius, NSEG, rgb1); + } + } + init_billiard_color(); + for (k=0; k= 1) + { + if (circles[k].active == 2) + { +// hsl_to_rgb(150.0, 0.9, 0.4, rgb); +// glColor3f(rgb[0], rgb[1], rgb[2]); + glColor3f(0.0, 1.0, 0.0); + rgb[0] = 0.0; rgb[1] = 0.9; rgb[2] = 0.0; + draw_colored_circle(circles[k].xc, circles[k].yc, circles[k].radius, NSEG, rgb); + } + else draw_circle(circles[k].xc, circles[k].yc, circles[k].radius, NSEG); + init_billiard_color(); + } + + /* draw shooter position for laser pattern */ + if (CIRCLE_PATTERN == C_LASER) + { + hsl_to_rgb(0.0, 0.9, 0.5, rgb); + glColor3f(rgb[0], rgb[1], rgb[2]); + + draw_circle(x_shooter, y_shooter, circles[ncircles-1].radius, NSEG); + } + + init_billiard_color(); + + glBegin(GL_LINE_LOOP); + glVertex2d(lambda, -1.0); + glVertex2d(lambda, 1.0); + glVertex2d(-lambda, 1.0); + glVertex2d(-lambda, -1.0); + glEnd(); + break; + } + case (D_CIRCLES_IN_GENUSN): + { +// for (k=0; k= 1) + for (k=0; k= 1)&&(in_polygon(circles[k].xc, circles[k].yc, 1.0, NPOLY, APOLY))) + { + if (circles[k].active == 2) + { + hsl_to_rgb(150.0, 0.9, 0.4, rgb); + glColor3f(rgb[0], rgb[1], rgb[2]); + } + draw_circle(circles[k].xc, circles[k].yc, circles[k].radius, NSEG); + init_billiard_color(); + } + + /* draw shooter position for laser pattern */ + if ((CIRCLE_PATTERN == C_LASER)||(CIRCLE_PATTERN == C_LASER_GENUSN)) + { + hsl_to_rgb(0.0, 0.9, 0.5, rgb); + glColor3f(rgb[0], rgb[1], rgb[2]); + + draw_circle(x_shooter, y_shooter, circles[ncircles-1].radius, NSEG); + } + + /* draw polygon */ + init_billiard_color(); + + omega = DPI/((double)NPOLY); + glBegin(GL_LINE_LOOP); + for (i=0; i<=NPOLY; i++) + { + x = cos(i*omega + APOLY*PID); + y = sin(i*omega + APOLY*PID); + glVertex2d(scaling_factor*x, scaling_factor*y); + } + glEnd (); + break; + } + case (D_CIRCLES_IN_TORUS): + { + rgb[0] = 0.0; rgb[1] = 0.0; rgb[2] = 0.0; + for (k=0; k= 1) + { + rgb_color_scheme_lum(circles[k].color, 0.85, rgb1); + circles[k].new--; + } + else rgb_color_scheme(circles[k].color, rgb1); + draw_colored_circle(circles[k].xc, circles[k].yc, circles[k].radius, NSEG, rgb1); + } + } + init_billiard_color(); + for (k=0; k= 1) + { + if (circles[k].active == 2) + { + hsl_to_rgb(150.0, 0.9, 0.4, rgb); + glColor3f(rgb[0], rgb[1], rgb[2]); + } + draw_circle(circles[k].xc, circles[k].yc, circles[k].radius, NSEG); + init_billiard_color(); + } + + /* draw shooter position for laser pattern */ + if (CIRCLE_PATTERN == C_LASER) + { + hsl_to_rgb(0.0, 0.9, 0.5, rgb); + glColor3f(rgb[0], rgb[1], rgb[2]); + + draw_circle(x_shooter, y_shooter, circles[ncircles-1].radius, NSEG); + } + + init_billiard_color(); + + glBegin(GL_LINE_LOOP); + glVertex2d(lambda, -1.0); + glVertex2d(lambda, 1.0); + glVertex2d(-lambda, 1.0); + glVertex2d(-lambda, -1.0); + glEnd(); + break; + } + case (D_POLYLINE): + { + for (k=0; k4.0*lambda + 4.0) s = 4.0*lambda + 4.0; + + theta = conf[1]; + + /* we treat the cases of starting in one corner separately */ + /* to avoid numerical problems due to hitting a corner */ + + /* bottom left corner */ + if ((s==0)||(s==4.0*lambda + 4.0)) + { + pos[0] = -lambda; + pos[1] = -1.0; + if (theta > PID) theta = PID; + *alpha = theta; + return(0); + } + /* bottom right corner */ + else if (s==2.0*lambda) + { + pos[0] = lambda; + pos[1] = -1.0; + if (theta > PID) theta = PID; + *alpha = theta + PID; + return(0); + } + /* top right corner */ + else if (s==2.0*lambda + 2.0) + { + pos[0] = lambda; + pos[1] = -1.0; + if (theta > PID) theta = PID; + *alpha = theta + PI; + return(0); + } + /* top left corner */ + else if (s==4.0*lambda + 2.0) + { + pos[0] = lambda; + pos[1] = -1.0; + if (theta > PID) theta = PID; + *alpha = theta + 3.0*PID; + return(0); + } + /* bottom side */ + else if ((s>0)&&(s<2.0*lambda)) + { + pos[0] = s - lambda; + pos[1] = -1.0; + *alpha = theta; + return(1); + } + /* right side */ + else if (s<2.0*lambda + 2.0) + { + pos[0] = lambda; + pos[1] = s - 2.0*lambda - 1.0; + *alpha = theta + PID; + return(2); + } + /* top side */ + else if (s<4.0*lambda + 2.0) + { + pos[0] = 3.0*lambda + 2.0 - s; + pos[1] = 1.0; + *alpha = theta + PI; + return(3); + } + /* left side */ + else + { + pos[0] = -lambda; + pos[1] = 4.0*lambda + 3.0 - s; + *alpha = theta + 3.0*PID; + return(4); + } + + } + + + int vrectangle_xy(double config[8], double alpha, double pos[2]) + /* determine initial configuration for start at point pos = (x,y) */ + { + double l, s0, c0, x1, y1, margin = 1e-12; + int c, intb=1; + + /* initial position and velocity */ + + s0 = sin(alpha); + c0 = cos(alpha); + + /* intersection with lower part of boundary */ + if (s0<0.0) + { + x1 = pos[0] - c0*(1.0 + pos[1])/s0; + y1 = -1.0; + if ((x1>=-lambda)&&(x1<=lambda)) + { + config[0] = x1 + lambda; + if ((x1 <= -lambda + margin)||(x1 >= lambda -margin)) config[1] = alpha + PI; /* corners */ +// if ((x1 == -lambda)||(x1 == lambda)) config[1] = alpha + PI; /* corners */ + else config[1] = -alpha; + intb = 0; + } + } + /* intersection with right-hand part of boundary */ + if (intb&&(c0>0.0)) + { + x1 = lambda; + y1 = pos[1] + s0*(lambda - pos[0])/c0; + if ((y1>=-1.0)&&(y1<=1.0)) + { + config[0] = 2.0*lambda + 1.0 + y1; + if ((y1 <= -1.0 + margin)||(y1 >= 1.0 -margin)) config[1] = alpha + PI; /* corners */ +// if ((y1 == -1.0)||(y1 == 1.0)) config[1] = alpha + PI; /* corners */ + else config[1] = PID-alpha; + intb = 0; + } + } + /* intersection with upper part of boundary */ + if (intb&&(s0>0.0)) + { + x1 = pos[0] + c0*(1.0 - pos[1])/s0; + y1 = 1.0; + if ((x1>=-lambda)&&(x1<=lambda)) + { + config[0] = 3.0*lambda + 2.0 - x1; + if ((x1 <= -lambda + margin)||(x1 >= lambda -margin)) config[1] = alpha + PI; /* corners */ +// if ((x1 == -lambda)||(x1 == lambda)) config[1] = alpha + PI; /* corners */ + else config[1] = PI-alpha; + intb = 0; + } + } + /* intersection with left-hand part of boundary */ + if (intb&&(c0<0.0)) + { + x1 = -lambda; + y1 = pos[1] + s0*(-lambda - pos[0])/c0; + if ((y1>=-1.0)&&(y1<=1.0)) + { + config[0] = 4.0*lambda + 3.0 - y1; + if ((y1 <= -1.0 + margin)||(y1 >= 1.0 -margin)) config[1] = alpha + PI; /* corners */ +// if ((y1 == -1.0)||(y1 == 1.0)) config[1] = alpha + PI; /* corners */ + else config[1] = 3.0*PID-alpha; + intb = 0; + } + } + + if (config[1] < 0.0) config[1] += DPI; + + config[2] = 0.0; /* running time */ + config[3] = module2(x1-pos[0], y1-pos[1]); + config[4] = pos[0]; + config[5] = pos[1]; + config[6] = x1; + config[7] = y1; + + return(c); + } + + int vrectangle(double config[8]) + { + double pos[2], alpha; + int c; + + /* position et vitesse de depart */ + + c = pos_rectangle(config, pos, &alpha); + + vrectangle_xy(config, alpha, pos); + + return(c); + } + + +/****************************************************************************************/ +/* elliptic billiard */ +/****************************************************************************************/ + + int pos_ellipse(double conf[2], double pos[2], double *alpha) + /* determine position on boundary of ellipse */ + { + double theta; + + pos[0] = lambda*cos(conf[0]); + pos[1] = sin(conf[0]); + + theta = argument(-lambda*pos[1],pos[0]/lambda); + *alpha = theta + conf[1]; + + return(1); + } + + + int vellipse_xy(double config[8], double alpha, double pos[2]) + /* determine initial configuration for start at point pos = (x,y) */ + { + double c0, s0, lam2, a, b, c, x1, y1, t, theta; + int i; + + c0 = cos(alpha); + s0 = sin(alpha); + lam2 = lambda*lambda; + + /* intersection with ellipse, using parametric equation of line */ + a = c0*c0 + lam2*s0*s0; + b = pos[0]*c0 + lam2*pos[1]*s0; + c = pos[0]*pos[0] + lam2*pos[1]*pos[1] - lam2; + + t = (-b+sqrt(b*b - a*c))/a; + x1 = pos[0] + t*c0; + y1 = pos[1] + t*s0; + + /* parameter of intersection with boundary ellipse */ + config[0] = argument(x1/lambda, y1); + while (config[0] < 0.0) config[0] += DPI; + while (config[0] > DPI) config[0] -= DPI; + + /* computation of outgoing angle after collision with boundary */ + theta = argument(-lambda*y1,x1/lambda); + config[1] = theta - alpha; + while (config[1] < 0.0) config[1] += DPI; + while (config[1] > DPI) config[1] -= DPI; + + config[2] = 0.0; /* running time */ + config[3] = module2(x1-pos[0], y1-pos[1]); /* distance to collision */ + config[4] = pos[0]; /* start position */ + config[5] = pos[1]; + config[6] = x1; /* position of collision */ + config[7] = y1; + + return(1); + } + + int vellipse(double config[8]) + /* determine initial configuration when starting from boundary */ + { + double pos[2], theta, alpha; + int i; + + pos[0] = lambda*cos(config[0]); + pos[1] = sin(config[0]); + + theta = argument(-lambda*pos[1],pos[0]/lambda); + alpha = theta + config[1]; + + vellipse_xy(config, alpha, pos); + + return(1); + } + +/****************************************************************************************/ +/* stadium billiard */ +/****************************************************************************************/ + + int pos_stade(double conf[2], double pos[2], double *alpha) + /* determine position on boundary of stadium */ + { + double s, theta, l, psi0, psi; + + s = conf[0]; + theta = conf[1]; + l = lambda/2.0; + + if (l >= 0.0) + { + if ((s>=0)&&(s<=lambda)) + { + pos[0] = s - l; + pos[1] = -1.0; + *alpha = theta; + return(0); + } + else if (s<=lambda+PI) + { + pos[0] = l + sin(s - lambda); + pos[1] = -cos(s - lambda); + *alpha = theta + s - lambda; + return(1); + } + else if (s<=2.0*lambda+PI) + { + pos[0] = 3.0*l + PI - s; + pos[1] = 1.0; + *alpha = theta + PI; + return(2); + } + else + { + pos[0] = -l - sin(s - 2.0*lambda - PI); + pos[1] = cos(s - 2.0*lambda - PI); + *alpha = theta + s - 2.0*lambda; + return(3); + } + } + else /* for lens-shaped billiard, to be checked */ + { + psi0 = asin(-l); + if ((s>=0)&&(s<=PI-2.0*psi0)) + { + psi = s + psi0; + pos[0] = sin(psi) + l; + pos[1] = -cos(psi); + *alpha = theta + psi; + return(0); + } + else + { + psi = s + 3.0*psi0 - PI; + pos[0] = - sin(psi) - l; + pos[1] = cos(psi); + *alpha = theta + psi + PI; + return(2); + } + } + } + + + int vstade_xy(double config[8], double alpha, double pos[2]) + /* determine initial configuration for start at point pos = (x,y) */ + { + double l, s0, c0, t, x, y, x1, y1, a, b, res[2]; + double smin, psi, margin = 1e-12; + int c, intb=1, intc, i; + + /* initial position and velocity */ + + l = lambda/2.0; + if (l>=0.0) smin = 0.0; else smin = -l; + s0 = sin(alpha); + c0 = cos(alpha); + + /* intersection with lower straight part of boundary */ + if ((s0<0.0)&&(l>0)) + { + x1 = pos[0] + c0*(-1.0 - pos[1])/s0; + y1 = -1.0; + if ((x1>=-l)&&(x1<=l)) + { + config[0] = x1 + l; + config[1] = -alpha; + intb = 0; + } + } + /* intersection with upper straight part of boundary */ + if (intb&&(s0>0.0)&&(l>0)) + { + x1 = pos[0] + c0*(1.0 - pos[1])/s0; + y1 = 1.0; + if ((x1>=-l)&&(x1<=l)) + { + config[0] = 3.0*l + PI - x1; + config[1] = PI-alpha; + intb = 0; + } + } + /* intersection with right-hand arc of boundary */ + if (intb) + { + a = 2.0*pos[0]*c0 + 2.0*pos[1]*s0 - lambda*c0; + b = pos[0]*pos[0] + pos[1]*pos[1] + l*l - lambda*pos[0] - 1.0; + intc = polynome(1.0, a, b, res); + if (intc) for(i=0; i<2; i++) + { + x = pos[0] + c0*res[i]; + y = pos[1] + s0*res[i]; + psi = argument(-y, x-l); + if (intb&&(sin(psi) >= smin)&&(res[i]>margin)) + { + if (l>0.0) config[0] = lambda + psi; + else config[0] = psi - asin(-l); + config[1] = -alpha + psi; + intb = 0; + x1 = x; y1 = y; + } + } + } + /* intersection with left-hand arc of boundary */ + if (intb) + { + a = 2.0*pos[0]*c0 + 2.0*pos[1]*s0 + lambda*c0; + b = pos[0]*pos[0] + pos[1]*pos[1] + l*l + lambda*pos[0] - 1.0; + intc = polynome(1.0, a, b, res); + if (intc) for(i=0; i<2; i++) + { + x = pos[0] + c0*res[i]; + y = pos[1] + s0*res[i]; + psi = argument(y, -l-x); + if (intb&&(sin(psi) >= smin)&&(res[i]>margin)) + { + if (l>0.0) config[0] = 2.0*lambda + PI + psi; + else config[0] = psi - 3.0*asin(-l) + PI; + config[1] = -alpha + psi + PI; + intb = 0; + x1 = x; y1 = y; + } + } + } + + config[2] = 0.0; /* running time */ + config[3] = module2(x1-pos[0], y1-pos[1]); + config[4] = pos[0]; + config[5] = pos[1]; + config[6] = x1; + config[7] = y1; + + return(c); + } + + int vstade(double config[8]) + { + double alpha, pos[2]; + int c; + + c = pos_stade(config, pos, &alpha); + + vstade_xy(config, alpha, pos); + + return(c); + } + +/****************************************************************************************/ +/* Sinai billiard */ +/****************************************************************************************/ + + int pos_sinai(double conf[2], double pos[2], double *alpha) + /* determine position on boundary of Sinai billiard */ + /* s in [0,2 Pi) is on the circle, other s are on boundary of window */ + { + double s, theta, psi0, psi, s1, s2, s3, s4; + + s = conf[0]; + theta = conf[1]; + if (conf[1] < 0.0) conf[1] += DPI; + + s1 = DPI + XMAX - XMIN; + s2 = s1 + YMAX - YMIN; + s3 = s2 + XMAX - XMIN; + s4 = s3 + YMAX - YMIN; + + if (s < DPI) /* circle */ + { + pos[0] = lambda*cos(s); + pos[1] = lambda*sin(s); + theta = PID + s; + *alpha = theta - conf[1]; + return(0); + } + else if (s < s1) /* boundary of window */ + { + pos[0] = XMIN + s - DPI; + pos[1] = YMIN; + *alpha = conf[1]; + return(-1); + } + else if (s < s2) + { + pos[0] = XMAX; + pos[1] = YMIN + s - s1; + *alpha = conf[1]; + return(-2); + } + else if (s < s3) + { + pos[0] = XMAX - s + s2; + pos[1] = YMAX; + *alpha = conf[1]; + return(-3); + } + else + { + pos[0] = XMIN; + pos[1] = YMAX - s + s3; + *alpha = conf[1]; + return(-4); + } + } + + + int vsinai_xy(double config[8], double alpha, double pos[2]) + /* determine initial configuration for start at point pos = (x,y) */ + { + double l, s0, c0, t, t1, x, y, x1, y1, a, b, delta, res[2], s1, s2, s3, s4, s, r; + double psi, lam2, margin = 1e-12; + int c, intb=1, intc, i; + + /* initial position and velocity */ + + c0 = cos(alpha); + s0 = sin(alpha); + s1 = DPI + XMAX - XMIN; + s2 = s1 + YMAX - YMIN; + s3 = s2 + XMAX - XMIN; + s4 = s3 + YMAX - YMIN; + + /* intersection with circle, using parametric equation of line */ + b = pos[0]*c0 + pos[1]*s0; + a = pos[0]*pos[0] + pos[1]*pos[1] - lambda*lambda; + delta = b*b - a; + + if ((delta > margin)&&(a > margin)) + { + t = - b - sqrt(delta); + x1 = pos[0] + t*c0; + y1 = pos[1] + t*s0; + s = argument(x1,y1); + if (s<0.0) s += DPI; + if (s>=DPI) s -= DPI; + config[0] = s; + config[1] = 3.0*PID - s + alpha; + c = 0; + } + else if (c0 > 0.0) /* intersection with boundary of window */ + { + y1 = pos[1] + (XMAX - pos[0])*s0/c0; + if ((y1 >= YMIN)&&(y1 <= YMAX)) /* hitting right boundary */ + { + x1 = XMAX; + config[0] = s3 + YMAX - y1; + config[1] = alpha; + c = 2; + } + else if (s0 > 0.0) /* hitting upper boundary */ + { + x1 = pos[0] + (YMAX - pos[1])*c0/s0; + y1 = YMAX; + config[0] = DPI + x1 - XMIN; + config[1] = alpha; + c = 3; + } + else /* hitting lower boundary */ + { + x1 = pos[0] + (YMIN - pos[1])*c0/s0; + y1 = YMIN; + config[0] = s2 + XMAX - x1; + config[1] = alpha; + c = 1; + } + } + else if (c0 < 0.0) + { + y1 = pos[1] + (XMIN - pos[0])*s0/c0; + if ((y1 >= YMIN)&&(y1 <= YMAX)) /* hitting left boundary */ + { + x1 = XMIN; + config[0] = s1 + y1 - YMIN; + config[1] = alpha; + c = 4; + } + else if (s0 > 0.0) /* hitting upper boundary */ + { + x1 = pos[0] + (YMAX - pos[1])*c0/s0; + y1 = YMAX; + config[0] = DPI + x1 - XMIN; + config[1] = alpha; + c = 3; + } + else /* hitting lower boundary */ + { + x1 = pos[0] + (YMIN - pos[1])*c0/s0; + y1 = YMIN; + config[0] = s2 + XMAX - x1; + config[1] = alpha; + c = 1; + } + } + else /* vertical motion */ + { + if (s0 > 0.0) + { + x1 = pos[0]; + y1 = YMAX; + config[0] = DPI + x1 - XMIN; + config[1] = alpha; + c = 3; + } + else + { + x1 = pos[0]; + y1 = YMIN; + config[0] = s2 + XMAX - x1; + config[1] = alpha; + c = 1; + } + } + + if (config[1] < 0.0) config[1] += DPI; + + config[2] = 0.0; /* running time */ + config[3] = module2(x1-pos[0], y1-pos[1]); + config[4] = pos[0]; + config[5] = pos[1]; + config[6] = x1; + config[7] = y1; + + return(-c); + /* return a negative value if the disc is not hit, for color scheme */ + } + + int vsinai(double config[8]) + { + double alpha, pos[2]; + int c; + + /* position et vitesse de depart */ + + c = pos_sinai(config, pos, &alpha); + + vsinai_xy(config, alpha, pos); + + return(c); + } + + +/****************************************************************************************/ +/* triangle billiard */ +/****************************************************************************************/ + + + int pos_triangle(double conf[2], double pos[2], double *alpha) + /* determine position on boundary of triangle */ + /* the corners of the triangle are (-lambda,-1), (lambda,-1), (-lambda,1) */ + /* we use arclength for horizontal and vertical side, x for diagonal */ + { + double s, theta; + + s = conf[0]; + theta = conf[1]; + + if ((s>=0)&&(s<=2.0*lambda)) + { + pos[0] = s - lambda; + pos[1] = -1.0; + *alpha = theta; + return(0); + } + else if (s<=4.0*lambda) + { + pos[0] = 3.0*lambda - s; + pos[1] = -3.0 + s/lambda; + *alpha = theta + PI - argument(lambda, 1.0); + return(1); + } + else + { + pos[0] = -lambda; + pos[1] = 4.0*lambda + 1.0 - s; + *alpha = theta + 3.0*PID; + return(2); + } + + } + + + int vtriangle_xy(double config[8], double alpha, double pos[2]) + /* determine initial configuration for start at point pos = (x,y) */ + /* Warning: reflection in the corners is not yet implemented correctly */ + { + double s0, c0, t, x, y, x1, y1, psi; + int c, intb=1, intc, i; + + /* initial position and velocity */ + + s0 = sin(alpha); + c0 = cos(alpha); + + /* intersection with lower part of boundary */ +// if ((s0<0.0)&&(pos[1]>0.0)) + if (s0<0.0) + { + x1 = pos[0] - c0*(1.0 + pos[1])/s0; + y1 = -1.0; + if ((x1>=-lambda)&&(x1<=lambda)) + { + config[0] = x1 + lambda; + config[1] = -alpha; + intb = 0; + } + } + /* intersection with left-hand part of boundary */ + if (intb&&(c0<0.0)) + { + x1 = -lambda; + y1 = pos[1] + s0*(-lambda - pos[0])/c0; + if ((y1>=-1.0)&&(y1<=1.0)) + { + config[0] = 4.0*lambda + 1.0 - y1; + config[1] = 3.0*PID-alpha; + intb = 0; + } + } + /* intersection with diagonal part of boundary */ + if (intb) + { + t = -(pos[0] + lambda*pos[1])/(c0 + lambda*s0); + x1 = pos[0] + t*c0; + y1 = pos[1] + t*s0; + if ((x1>=-lambda)&&(x1<=lambda)) + { + psi = argument(lambda, 1.0); + config[0] = 3.0*lambda - x1; + config[1] = PI - alpha - psi; +// config[1] = PI - alpha - atan(1.0/lambda); + intb = 0; + } + } + + config[2] = 0.0; /* running time */ + config[3] = module2(x1-pos[0], y1-pos[1]); + config[4] = pos[0]; + config[5] = pos[1]; + config[6] = x1; + config[7] = y1; + + return(c); + } + + int vtriangle(double config[8]) + { + double alpha, pos[2]; + int c; + + /* position et vitesse de depart */ + + c = pos_triangle(config, pos, &alpha); + + vtriangle_xy(config, alpha, pos); + + return(c); + } + + +/****************************************************************************************/ + /* annulus billiard */ +/****************************************************************************************/ + + int pos_annulus(double conf[2], double pos[2], double *alpha) + /* determine position on boundary of annulus */ + { + double s, theta, psi0, psi, s1, s2, s3, s4; + + s = conf[0]; + theta = conf[1]; + if (conf[1] < 0.0) conf[1] += DPI; + + if (conf[0] < DPI) /* inner circle */ + { + pos[0] = lambda*cos(conf[0]) + mu; + pos[1] = lambda*sin(conf[0]); + + theta = PID + conf[0]; + *alpha = theta - conf[1]; + return(0); + } + else /* outer circle */ + { + pos[0] = cos(conf[0]); + pos[1] = sin(conf[0]); + + theta = argument(-pos[1],pos[0]); + *alpha = theta + conf[1]; + return(1); + } + } + + int vannulus_xy(double config[8], double alpha, double pos[2]) + /* determine initial configuration for start at point pos = (x,y) */ + { + double l, s0, c0, t, t1, x, y, x1, y1, a, b, delta, res[2], s, r; + double psi, lam2, margin = 1.0e-14, theta; +// double psi, lam2, margin = 1.0e-14, theta; + int c, intb=1, intc, i; + + /* initial position and velocity */ + + c0 = cos(alpha); + s0 = sin(alpha); + + /* intersection with inner circle, using parametric equation of line */ + b = (pos[0]-mu)*c0 + pos[1]*s0; + a = (pos[0]-mu)*(pos[0]-mu) + pos[1]*pos[1] - lambda*lambda; + delta = b*b - a; + + if ((delta > margin)&&(a > margin)) + { + t = - b - sqrt(delta); + x1 = pos[0] + t*c0; + y1 = pos[1] + t*s0; + s = argument(x1-mu,y1); + while (s<0.0) s += DPI; + while (s>=DPI) s -= DPI; + config[0] = s; + config[1] = 3.0*PID - s + alpha; + c = 0; + } + else /* intersection with outer circle, using parametric equation of line */ + { + b = pos[0]*c0 + pos[1]*s0; + a = pos[0]*pos[0] + pos[1]*pos[1] - 1.0; + + t = (-b+sqrt(b*b - a)); + x1 = pos[0] + t*c0; + y1 = pos[1] + t*s0; + + /* parameter of intersection with outer circle */ + config[0] = argument(x1, y1); + while (config[0] < DPI) config[0] += DPI; + while (config[0] >= 2.0*DPI) config[0] -= DPI; + + /* computation of outgoing angle after collision with outer circle */ + theta = argument(-y1,x1); + config[1] = theta - alpha; +// while (config[1] < 0.0) config[1] += DPI; +// while (config[1] > DPI) config[1] -= DPI; + c = 1; + } + + if (config[1] < 0.0) config[1] += DPI; + +// config[2] = 1.0e-12; /* running time */ + config[2] = 0.0; /* running time */ + config[3] = module2(x1-pos[0], y1-pos[1]); /* distance to collision */ + config[4] = pos[0]; /* start position */ + config[5] = pos[1]; + config[6] = x1; /* position of collision */ + config[7] = y1; + + return(c); + } + + int vannulus(double config[8]) + /* determine initial configuration when starting from boundary */ + { + double pos[2], alpha; + int c; + + c = pos_annulus(config, pos, &alpha); + + vannulus_xy(config, alpha, pos); + + return(c); + } + +/****************************************************************************************/ +/* polygonal billiard */ +/****************************************************************************************/ + + int pos_polygon(double conf[2], double pos[2], double *alpha) + /* determine position on boundary of polygon */ + /* conf[0] is arclength on boundary */ + { + double s, theta, omega, length, s1, angle, x, y; + int c; + + s = conf[0]; + theta = conf[1]; + + omega = DPI/((double)NPOLY); + length = 2.0*sin(0.5*omega); + + c = (int)(s/length); /* side of polygon */ + + s1 = s - ((double)c)*length; + + x = 1.0 + (cos(omega) - 1.0)*s1/length; + y = sin(omega)*s1/length; + + angle = (double)c*omega + PID*APOLY; + + pos[0] = x*cos(angle) - y*sin(angle); + pos[1] = x*sin(angle) + y*cos(angle); + + *alpha = (0.5 + (double)c)*omega + theta + PID*(1.0 + APOLY); + + return(c); + } + + int vpolygon_xy(double config[8], double alpha, double pos[2]) + /* determine initial configuration for start at point pos = (x,y) */ + { + double s, theta, omega, length, rlength, s1, rangle, x, y, xp, yp, x1, y1, ca, sa; + int k, c, intb=1, intc, i; + + /* dimensions/angles of polygon */ + omega = DPI/((double)NPOLY); + length = 2.0*sin(0.5*omega); + rlength = cos(0.5*omega); + + for (k=0; k 0.0)&&(intb)) + { + ca = cos(rangle); + sa = sin(rangle); + + x = pos[0]*ca + pos[1]*sa; + y = -pos[0]*sa + pos[1]*ca; + + xp = rlength; + yp = y + (xp-x)*tan(theta); + + if (vabs(yp) < 0.5*length) + { + /* rotate back */ + x1 = xp*ca - yp*sa; + y1 = xp*sa + yp*ca; + + intb = 0; + c = k; + config[0] = ((double)k + 0.5)*length + yp; + config[1] = PID - theta; + } + } + } + + if (config[1] < 0.0) config[1] += DPI; + + config[2] = 0.0; /* running time */ + config[3] = module2(x1-pos[0], y1-pos[1]); /* distance to collision */ + config[4] = pos[0]; /* start position */ + config[5] = pos[1]; + config[6] = x1; /* position of collision */ + config[7] = y1; + + return(c); + } + + int vpolygon(double config[8]) + /* determine initial configuration when starting from boundary */ + { + double pos[2], alpha; + int c; + + c = pos_polygon(config, pos, &alpha); + + vpolygon_xy(config, alpha, pos); + + return(c); + } + +/****************************************************************************************/ +/* Reuleaux-type and star-shaped billiard */ +/****************************************************************************************/ + + int pos_reuleaux(double conf[2], double pos[2], double *alpha) + /* determine position on boundary of polygon */ + /* conf[0] is arclength on boundary */ + { + double s, theta, omega2, beta2, beta, s1, angle, x2, x, y; + int c; + + s = conf[0]; + theta = conf[1]; + + omega2 = PI/((double)NPOLY); + beta2 = asin(sin(omega2)/vabs(lambda)); + beta = beta2*2.0; + + c = (int)(s/beta); /* side of shape */ + + s1 = s - ((double)c)*beta; + + if (lambda > 0.0) x2 = cos(omega2) + sqrt(lambda*lambda - sin(omega2)*sin(omega2)); + else x2 = cos(omega2) - sqrt(lambda*lambda - sin(omega2)*sin(omega2)); + + x = x2 - lambda*cos(s1 - beta2); + y = lambda*sin(s1 - beta2); + + angle = 2.0*((double)c)*omega2 + PID*APOLY; + + pos[0] = x*cos(angle) - y*sin(angle); + pos[1] = x*sin(angle) + y*cos(angle); + + *alpha = PID - s1 + beta2 + theta + 2.0*(double)c*omega2 + APOLY*PID; + +// printf("alpha = %.5lg\t", *alpha); + + return(c); + } + + int vreuleaux_xy(double config[8], double alpha, double pos[2]) + /* determine initial configuration for start at point pos = (x,y) */ + { + double s, theta, omega2, beta, s1, rangle, x, y, x1[NPOLY], y1[NPOLY], xi, yi, t, x2; + double ca, sa, a, b, margin = 1.0e-14, tmin, tval[NPOLY], tempconf[NPOLY][2]; + int k, c, intb=1, intc, i, nt = 0, cval[NPOLY], ntmin; + + /* dimensions/angles of polygon */ + omega2 = PI/((double)NPOLY); + beta = 2.0*asin(sin(omega2)/vabs(lambda)); +// printf("beta = %.5lg\n", beta); + + if (lambda > 0.0) x2 = cos(omega2) + sqrt(lambda*lambda - sin(omega2)*sin(omega2)); + else x2 = cos(omega2) - sqrt(lambda*lambda - sin(omega2)*sin(omega2)); +// printf("x2 = %.5lg\n", x2); + + for (k=0; k margin) + { + if (lambda > 0.0) t = -a - sqrt(a*a - b); + else t = -a + sqrt(a*a - b); + + xi = x + t*cos(theta); + yi = y + t*sin(theta); + + if ((t > margin)&&(vabs(yi) <= sin(omega2))) + { + cval[nt] = k; + tval[nt] = t; + + /* rotate back */ + x1[nt] = xi*ca - yi*sa; + y1[nt] = xi*sa + yi*ca; + + tempconf[nt][0] = ((double)k + 0.5)*beta + asin(yi/lambda); + tempconf[nt][1] = PID - asin(yi/lambda) - theta; + nt++; + } + } + } + + /* find earliest intersection */ + tmin = tval[0]; + ntmin = 0; + for (i=1; i margin) + if (bb*bb - aa*cc >= 0.0) + { + t = (-bb + sqrt(bb*bb - aa*cc))/aa; + + xi = x + t*cos(theta); + yi = y + t*sin(theta); + + if (yi >= 0.0) phi = argument((xi - co)/axis1, yi/axis2); + else phi = -argument((xi - co)/axis1, -yi/axis2); + + if ((t > margin)&&((vabs(phi) <= phimax)||(vabs(phi-DPI) <= phimax))) + { + intb = 0; + c = k; + + /* rotate back */ + x1 = xi*ca - yi*sa; + y1 = xi*sa + yi*ca; + + config[0] = (double)(2*k + 1)*phimax + phi; + config[1] = argument(-axis1*sin(phi), axis2*cos(phi)) - theta; + } + } + } + +// if (nt == 0) printf("nt = %i\t ntmin = %i \tcmin = %i\n", nt, ntmin, c); + + if (config[1] < 0.0) config[1] += DPI; + + config[2] = 0.0; /* running time */ + config[3] = module2(x1-pos[0], y1-pos[1]); /* distance to collision */ + config[4] = pos[0]; /* start position */ + config[5] = pos[1]; + config[6] = x1; /* position of collision */ + config[7] = y1; + + return(c); + } + + + int vflower(double config[8]) + /* determine initial configuration when starting from boundary */ + { + double pos[2], alpha; + int c; + + c = pos_flower(config, pos, &alpha); + + vflower_xy(config, alpha, pos); + + return(c); + } + + /****************************************************************************************/ +/* Alternating between Reuleaux-type and star-shaped billiard */ +/****************************************************************************************/ + + int pos_alt_reuleaux(double conf[2], double pos[2], double *alpha) + /* determine position on boundary of polygon */ + /* conf[0] is arclength on boundary */ + { + double s, theta, omega2, beta2, beta, s1, angle, x2plus, x2minus, x, y; + int c; + + s = conf[0]; + theta = conf[1]; + + omega2 = PI/((double)NPOLY); + beta2 = asin(sin(omega2)/vabs(lambda)); + beta = beta2*2.0; + + c = (int)(s/beta); /* side of shape */ + + s1 = s - ((double)c)*beta; + + x2plus = cos(omega2) + sqrt(lambda*lambda - sin(omega2)*sin(omega2)); + x2minus = cos(omega2) - sqrt(lambda*lambda - sin(omega2)*sin(omega2)); + + if (c%2 == 0) x = x2plus - lambda*cos(s1 - beta2); + else x = x2minus + lambda*cos(s1 - beta2); + + if (c%2 == 0) y = lambda*sin(s1 - beta2); + else y = -lambda*sin(s1 - beta2); + + /* test, to be removed */ +// x = x2plus - lambda*cos(s1 - beta2); +// y = lambda*sin(s1 - beta2); + + angle = 2.0*((double)c)*omega2 + PID*APOLY; + + pos[0] = x*cos(angle) - y*sin(angle); + pos[1] = x*sin(angle) + y*cos(angle); + + *alpha = PID - s1 + beta2 + theta + 2.0*(double)c*omega2 + APOLY*PID; + +// printf("alpha = %.5lg\t", *alpha); + + return(c); + } + + int valt_reuleaux_xy(double config[8], double alpha, double pos[2]) + /* determine initial configuration for start at point pos = (x,y) */ + { + double s, theta, omega2, beta, s1, rangle, x, y, x1[NPOLY], y1[NPOLY], xi, yi, t, x2plus, x2minus, arcsine; + double ca, sa, a, b, margin = 1.0e-14, tmin, tval[NPOLY], tempconf[NPOLY][2]; + int k, c, intb=1, intc, i, nt = 0, cval[NPOLY], ntmin; + + /* dimensions/angles of polygon */ + omega2 = PI/((double)NPOLY); + beta = 2.0*asin(sin(omega2)/vabs(lambda)); +// printf("beta = %.5lg\n", beta); + + x2plus = cos(omega2) + sqrt(lambda*lambda - sin(omega2)*sin(omega2)); + x2minus = cos(omega2) - sqrt(lambda*lambda - sin(omega2)*sin(omega2)); +// printf("x2 = %.5lg\n", x2); + + for (k=0; k margin) + { + if (k%2==0) t = -a - sqrt(a*a - b); + else t = -a + sqrt(a*a - b); + + xi = x + t*cos(theta); + yi = y + t*sin(theta); + + if ((t > margin)&&(vabs(yi) <= sin(omega2))) + { + cval[nt] = k; + tval[nt] = t; + + /* rotate back */ + x1[nt] = xi*ca - yi*sa; + y1[nt] = xi*sa + yi*ca; + + if (k%2==0) arcsine = asin(yi/lambda); + else arcsine = -asin(yi/lambda); + + tempconf[nt][0] = ((double)k + 0.5)*beta + arcsine; + tempconf[nt][1] = PID - arcsine - theta; + nt++; + } + } + } + + /* find earliest intersection */ + tmin = tval[0]; + ntmin = 0; + for (i=1; i0.0)) + { + phi = lambda*PI; + t = -(pos[0]*sin(phi) + pos[1]*cos(phi))/sin(alpha+phi); + x1 = pos[0] + t*c0; + y1 = pos[1] + t*s0; + if (y1 > 0.0) + { + config[0] = y1; + config[1] = PI - alpha - phi; + intb = 0; + } + } + /* other cases, assign a dummy coordinate */ + if (intb) + { + if (s0 > 0.0) + { + t = (1000.0 - pos[0])/s0; + x1 = pos[0] + t*c0; + y1 = 1000.0; + } + else + { + t = -(1000.0 + pos[0])/c0; + x1 = -1000.0; + y1 = pos[1] + t*s0; + } + config[0] = -1000.0; + config[1] = PID; + } + + config[2] = 0.0; /* running time */ + config[3] = module2(x1-pos[0], y1-pos[1]); + config[4] = pos[0]; + config[5] = pos[1]; + config[6] = x1; + config[7] = y1; + + return(c); + } + + int vangle(double config[8]) + { + double alpha, pos[2]; + int c; + + /* position et vitesse de depart */ + + c = pos_angle(config, pos, &alpha); + + vangle_xy(config, alpha, pos); + + return(c); + } + +/****************************************************************************************/ +/* L-shaped billiard (conical singularity) */ +/****************************************************************************************/ + + + int pos_lshape(double conf[2], double pos[2], double *alpha) + /* determine position on boundary of L-shaped billiard */ + /* the corners of the L shape are (-1,-1), (1,-1), (1,0), (0,0), (0,1), (-1,1) */ + /* returns the number of the side hit, or 0 if hitting a corner */ + { + double s, theta; + + s = conf[0]; + if (s<0.0) s = 0.0; + if (s>8.0) s = 8.0; + + theta = conf[1]; + + + /* bottom side */ + if ((s>0)&&(s<2.0)) + { + pos[0] = s -1.0; + pos[1] = -1.0; + *alpha = theta; + return(1); + } + /* lower right side */ + else if (s<3.0) + { + pos[0] = 1.0; + pos[1] = s - 3.0; + *alpha = theta + PID; + return(2); + } + /* lower top side */ + else if (s<4.0) + { + pos[0] = 4.0 - s; + pos[1] = 0.0; + *alpha = theta + PI; + return(3); + } + /* upper right side */ + else if (s<5.0) + { + pos[0] = 0.0; + pos[1] = s - 4.0; + *alpha = theta + PID; + return(4); + } + /* upper top side */ + else if (s<6.0) + { + pos[0] = 5.0 - s; + pos[1] = 1.0; + *alpha = theta + PI; + return(5); + } + /* left side */ + else + { + pos[0] = -1.0; + pos[1] = 7.0 - s; + *alpha = theta - PID; + return(6); + } + + } + + + int vlshape_xy(double config[8], double alpha, double pos[2]) + /* determine initial configuration for start at point pos = (x,y) */ + { + double t, l, s0, c0, margin = 1e-12, tmin; + double x1[6], y1[6], tval[6], tempconf[6][2]; + int i, c, intb=1, cval[6], nt = 0, ntmin; + + /* initial position and velocity */ + + s0 = sin(alpha); + c0 = cos(alpha); + +// print_config(config); + + /* intersection with lower part of boundary */ + if (s0<0.0) + { + t = - (1.0 + pos[1])/s0; + x1[nt] = pos[0] + c0*t; + y1[nt] = -1.0; + if ((x1[nt]>=-1.0)&&(x1[nt]<=0.0)) + { + tempconf[nt][0] = 5.0 - x1[nt]; + tempconf[nt][1] = alpha + PI; + cval[nt] = 5; + tval[nt] = t; + nt++; + } + else if ((x1[nt]>0.0)&&(x1[nt]<=1.0)) + { + tempconf[nt][0] = 4.0 - x1[nt]; + tempconf[nt][1] = alpha + PI; + cval[nt] = 3; + tval[nt] = t; + nt++; + } + } + /* intersection with lower part of right-hand boundary */ + if (c0>0.0) + { + t = (1.0 - pos[0])/c0; + x1[nt] = 1.0; + y1[nt] = pos[1] + s0*t; + if ((y1[nt]>=-1.0)&&(y1[nt]<=0.0)) + { + tempconf[nt][0] = 7.0 - y1[nt]; + tempconf[nt][1] = PID + alpha; + cval[nt] = 6; + tval[nt] = t; + nt++; + } + } + /* intersection with upper part of right-hand boundary */ + if ((pos[0] < 0.0)&&(c0>0.0)) + { + t = - pos[0]/c0; + x1[nt] = 0.0; + y1[nt] = pos[1] + s0*t; + if ((y1[nt]>=0.0)&&(y1[nt]<=1.0)) + { + tempconf[nt][0] = 7.0 - y1[nt]; + tempconf[nt][1] = PID + alpha; + cval[nt] = 4; + tval[nt] = t; + nt++; + } + } + /* intersection with right-hand part of upper boundary */ + if ((pos[1] < 0.0)&&(s0>0.0)) + { + t = - pos[1]/s0; + x1[nt] = pos[0] + c0*t; + y1[nt] = 0.0; + if ((x1[nt]>=-0.0)&&(x1[nt]<=1.0)) + { + tempconf[nt][0] = 1.0 + x1[nt]; + tempconf[nt][1] = alpha; + cval[nt] = 1; + tval[nt] = t; + nt++; + } + } + /* intersection with left-hand part of upper boundary */ + if (s0>0.0) + { + t = (1.0 - pos[1])/s0; + x1[nt] = pos[0] + c0*t; + y1[nt] = 1.0; + if ((x1[nt]>=-1.0)&&(x1[nt]<=0.0)) + { + tempconf[nt][0] = 1.0 + x1[nt]; + tempconf[nt][1] = alpha; + cval[nt] = 1; + tval[nt] = t; + nt++; + } + } + /* intersection with left-hand part of boundary */ + if (c0<0.0) + { + t = (-1.0 - pos[0])/c0; + x1[nt] = -1.0; + y1[nt] = pos[1] + s0*t; + if ((y1[nt]>=0.0)&&(y1[nt]<=1.0)) + { + tempconf[nt][0] = 4.0 + y1[nt]; + tempconf[nt][1] = alpha - PID; + cval[nt] = 4; + tval[nt] = t; + nt++; + } + else if ((y1[nt]>=-1.0)&&(y1[nt]<0.0)) + { + tempconf[nt][0] = 3.0 + y1[nt]; + tempconf[nt][1] = alpha - PID; + cval[nt] = 2; + tval[nt] = t; + nt++; + } + } + + /* find earliest intersection */ + tmin = tval[0]; + ntmin = 0; + for (i=1; i 0.0)&&(intb)) + { + ca = cos(rangle); + sa = sin(rangle); + + x = pos[0]*ca + pos[1]*sa; + y = -pos[0]*sa + pos[1]*ca; + +// xp = -rlength; + xp = rlength; + yp = y + (xp-x)*tan(theta); + + if (vabs(yp) < 0.5*length - margin) + { + /* rotate back */ + x1 = xp*ca - yp*sa; + y1 = xp*sa + yp*ca; + + intb = 0; + c = k + NPOLY/2; + if (c > NPOLY) c -= NPOLY; + config[0] = ((double)c + 0.5)*length - yp; +// if (config[0] > (double)NPOLY*length) config[0] -= (double)NPOLY*length; + config[1] = PID + theta; +// config[0] = ((double)k + 0.5)*length + yp; +// config[1] = PID - theta; + } + } + } + if (intb) x1 = 1.0e12; /* to make inactive particles too close to a corner */ + + if (config[1] < 0.0) config[1] += DPI; + + config[2] = 0.0; /* running time */ + config[3] = module2(x1-pos[0], y1-pos[1]); /* distance to collision */ + config[4] = pos[0]; /* start position */ + config[5] = pos[1]; + config[6] = x1; /* position of collision */ + config[7] = y1; + + return(c); + } + + int vgenusn(double config[8]) + /* determine initial configuration when starting from boundary */ + { + double pos[2], alpha; + int c; + + c = pos_genusn(config, pos, &alpha); + + vgenusn_xy(config, alpha, pos); + + return(c); + } + +/****************************************************************************************/ +/* Polygonal billiard with parabolic sides */ +/****************************************************************************************/ + + int pos_parabolas(double conf[2], double pos[2], double *alpha) + /* determine position on boundary of polygon */ + /* conf[0] is y + ymax on first side */ + { + double s, theta, omega2, aa, bb, cc, ymax, angle, x2, x, y; + int c; + + s = conf[0]; + theta = conf[1]; + + omega2 = PI/((double)NPOLY); + aa = 0.25/mu; + bb = 1.0/tan(omega2); + cc = lambda + mu; + ymax = ( - bb + sqrt(bb*bb + 4.0*aa*cc))/(2.0*aa); + + c = (int)(s/(2.0*ymax)); /* side of shape */ + + y = s - (2.0*(double)c + 1.0)*ymax; + x = lambda + mu - 0.25*y*y/mu; + +// printf("y = %.3lg\n", y); + + angle = 2.0*((double)c)*omega2 + PID*APOLY; + + pos[0] = x*cos(angle) - y*sin(angle); + pos[1] = x*sin(angle) + y*cos(angle); + + *alpha = PID + theta + atan(0.5*y/mu) + angle; + +// *alpha = PID + theta + atan(0.5*y/mu) + 2.0*(double)c*omega2 + APOLY*PID; + +// printf("alpha = %.5lg\t", *alpha); + + return(c); + } + + int vparabolas_xy(double config[8], double alpha, double pos[2]) + /* determine initial configuration for start at point pos = (x,y) */ + { + double s, theta, omega2, beta, s1, rangle, x, y, x1, y1, xi, yi, t, x2; + double ca, sa, a, b, margin = 1.0e-14, tmin, aa, bb, cc, ymax; + int k, c, intb=1, intc, i, nt = 0, ntmin; + + /* dimensions/angles of polygon */ + omega2 = PI/((double)NPOLY); + aa = 0.25/mu; + bb = 1.0/tan(omega2); + cc = lambda + mu; + ymax = ( - bb + sqrt(bb*bb + 4.0*aa*cc))/(2.0*aa); +// printf("ymax = %.3lg\n", ymax); + +// print_config(config); + + for (k=0; k= 1.0 - margin) + { + +// t = -cc/(2.0*bb); + xi = lambda + mu - 0.25*y*y/mu; + yi = y; + +// if ((t > margin)&&(vabs(yi) <= ymax)) + if (vabs(yi) <= ymax) + { + intb = 0; + + /* rotate back */ + x1 = xi*ca - yi*sa; + y1 = xi*sa + yi*ca; + + config[0] = yi + (1.0 + 2.0*(double)k)*ymax; + config[1] = PID - theta + atan(0.5*yi/mu); + } + +// printf("s = %.3lg\n", config[0]); + } + else if (bb*bb - aa*cc > margin) + { + t = (-bb + sqrt(bb*bb - aa*cc))/aa; + + xi = x + t*cos(theta); + yi = y + t*sin(theta); + + if ((t > margin)&&(vabs(yi) <= ymax)) +// if (vabs(yi) <= ymax) + { + intb = 0; + +// printf("t = %.3lg\n", t); + + /* rotate back */ + x1 = xi*ca - yi*sa; + y1 = xi*sa + yi*ca; + + config[0] = yi + (1.0 + 2.0*(double)k)*ymax; + config[1] = PID - theta + atan(0.5*yi/mu); + } + } +// if (!intb) printf("Intersection found with side %i\n", k); + } + +// if (intb) printf("No intersection found!\n"); + + if (config[1] < 0.0) config[1] += DPI; + + config[2] = 0.0; /* running time */ + config[3] = module2(x1-pos[0], y1-pos[1]); /* distance to collision */ + config[4] = pos[0]; /* start position */ + config[5] = pos[1]; + config[6] = x1; /* position of collision */ + config[7] = y1; + +// print_config(config); + + return(k); + } + + int vparabolas(double config[8]) + /* determine initial configuration when starting from boundary */ + { + double pos[2], alpha; + int c; + + c = pos_parabolas(config, pos, &alpha); + + vparabolas_xy(config, alpha, pos); + + return(c); + } + +/****************************************************************************************/ +/* Penrose solution to illumination problem */ +/****************************************************************************************/ + + int pos_penrose(double conf[2], double pos[2], double *alpha) + /* determine position on boundary of domain */ + /* conf[0] parametrization of boundary by arclength or angle */ + { + double s, s1, theta, cc, l1, l2, width, x, y, phi; + int c, i; + static int first = 1; + static double sval[17]; + + s = conf[0]; + theta = conf[1]; + + cc = sqrt(lambda*lambda - (1.0-mu)*(1.0-mu)); + width = 0.1*mu; + l1 = mu - width; + l2 = lambda - cc; + + /* s values of different boundary parts */ + if (first) + { + sval[0] = 0.0; sval[1] = PI; sval[2] = sval[1] + l1; + sval[3] = sval[2] + l2; sval[4] = sval[3] + l1; + sval[5] = sval[4] + PI; sval[6] = sval[5] + l1; + sval[7] = sval[6] + l2; sval[8] = sval[7] + l1; + for (i=1; i<=8; i++) sval[8+i] = sval[8] + sval[i]; +// for (i=0; i<16; i++) printf("sval[%i] = %.3lg\n", i, sval[i]); + first = 0; + } + + if (s < sval[1]) /* upper ellipse */ + { + pos[0] = lambda*cos(s); + pos[1] = mu + (1.0-mu)*sin(s); + phi = argument(-lambda*sin(s),(1.0-mu)*cos(s)); + *alpha = phi + theta; + return(0); + } + else if (s < sval[2]) /* upper left straight parts */ + { + pos[0] = -lambda; + pos[1] = mu - s + PI; + *alpha = theta - PID; + return(1); + } + else if (s < sval[3]) + { + pos[0] = -lambda + s - sval[2]; + pos[1] = width; + *alpha = theta; + return(2); + } + else if (s < sval[4]) + { + pos[0] = -cc; + pos[1] = width + s -sval[3]; + *alpha = theta + PID; + return(3); + } + else if (s < sval[5]) /* left ellipse/mushroom head */ + { + s1 = s - sval[4]; + pos[0] = -cc + 0.5*penrose_ratio*mu*sin(s1); + pos[1] = mu*cos(s1); + phi = argument(0.5*penrose_ratio*mu*cos(s1), -mu*sin(s1)); + *alpha = phi + theta; + return(4); + } + else if (s < sval[6]) /* lower left straight parts */ + { + s1 = s - sval[5]; + pos[0] = -cc; + pos[1] = -mu + s1; + *alpha = theta + PID; + return(5); + } + else if (s < sval[7]) + { + s1 = s - sval[6]; + pos[0] = -cc - s1; + pos[1] = -width; + *alpha = theta + PI; + return(6); + } + else if (s < sval[8]) + { + s1 = s - sval[7]; + pos[0] = -lambda; + pos[1] = -width - s1; + *alpha = theta - PID; + return(7); + } + + else if (s < sval[9]) /* lower ellipse */ + { + s1 = s - sval[8]; + pos[0] = -lambda*cos(s1); + pos[1] = -mu - (1.0-mu)*sin(s1); + phi = argument(lambda*sin(s1),-(1.0-mu)*cos(s1)); + *alpha = phi + theta; + return(8); + } + else if (s < sval[10]) /* lower right straight parts */ + { + s1 = s - sval[9]; + pos[0] = lambda; + pos[1] = -mu + s1; + *alpha = theta + PID; + return(9); + } + else if (s < sval[11]) + { + s1 = s - sval[10]; + pos[0] = lambda - s1; + pos[1] = -width; + *alpha = theta + PI; + return(10); + } + else if (s < sval[12]) + { + s1 = s - sval[11]; + pos[0] = cc; + pos[1] = -width -s1; + *alpha = theta - PID; + return(11); + } + else if (s < sval[13]) /* right ellipse/mushroom head */ + { + s1 = s - sval[12]; + pos[0] = cc - 0.5*penrose_ratio*mu*sin(s1); + pos[1] = -mu*cos(s1); + phi = argument(-0.5*penrose_ratio*mu*cos(s1), mu*sin(s1)); + *alpha = phi + theta; + return(12); + } + else if (s < sval[14]) /* upper right straight parts */ + { + s1 = s - sval[13]; + pos[0] = cc; + pos[1] = mu - s1; + *alpha = theta - PID; + return(13); + } + else if (s < sval[15]) + { + s1 = s - sval[14]; + pos[0] = cc + s1; + pos[1] = width; + *alpha = theta; + return(14); + } + else + { + s1 = s - sval[15]; + pos[0] = lambda; + pos[1] = width + s1; + *alpha = theta + PID; + return(15); + } + } + + int vpenrose_xy(double config[8], double alpha, double pos[2]) + /* determine initial configuration for start at point pos = (x,y) */ + { + double s, theta, cc, width, l1, l2, s1, rangle, x, y, x1[30], y1[30], xi, yi, t, x2; + double ca, sa, a, b, d, margin = 1.0e-14, tmin, tval[30], tempconf[30][2], lam2, mu2; + int k, c, intb=1, intc, i, nt = 0, cval[30], ntmin; + static int first = 1; + static double sval[17]; + + /* dimensions of domain */ + cc = sqrt(lambda*lambda - (1.0-mu)*(1.0-mu)); + width = 0.1*mu; + l1 = mu - width; + l2 = lambda - cc; + lam2 = lambda*lambda; + mu2 = (1.0-mu)*(1.0-mu); + + /* s values of different boundary parts */ + /* USE STATIC DOUBLES ? */ + if (first) + { + sval[0] = 0.0; sval[1] = PI; sval[2] = sval[1] + l1; + sval[3] = sval[2] + l2; sval[4] = sval[3] + l1; + sval[5] = sval[4] + PI; sval[6] = sval[5] + l1; + sval[7] = sval[6] + l2; sval[8] = sval[7] + l1; + for (i=1; i<=8; i++) sval[8+i] = sval[8] + sval[i]; +// for (i=0; i<16; i++) printf("sval[%i] = %.3lg\n", i, sval[i]); + first = 0; + } + + ca = cos(alpha); + sa = sin(alpha); + + /* intersection with upper ellipse */ + a = mu2*ca*ca + lam2*sa*sa; + b = mu2*pos[0]*ca + lam2*(pos[1]-mu)*sa; + d = mu2*pos[0]*pos[0] + lam2*(pos[1]-mu)*(pos[1]-mu) - lam2*mu2; + + if (b*b > a*d) + { + t = (-b + sqrt(b*b - a*d))/a; + x = pos[0] + t*ca; + y = pos[1] + t*sa; + +// printf("x = %.3lg, y = %.3lg\n", x, y); + + if ((t > margin)&&(y >= mu)) + { + cval[nt] = 0; + tval[nt] = t; + x1[nt] = x; + y1[nt] = y; + tempconf[nt][0] = argument(x/lambda, (y-mu)/(1.0-mu)); + tempconf[nt][1] = argument(-lambda*(y-mu)/(1.0-mu), (1.0-mu)*x/lambda) - alpha; +// tempconf[nt][1] = argument(-lambda*(y-mu)/(1.0-mu), x/mu) - alpha; + nt++; + } + } + + /* intersection with lower ellipse */ + b = mu2*pos[0]*ca + lam2*(pos[1]+mu)*sa; + d = mu2*pos[0]*pos[0] + lam2*(pos[1]+mu)*(pos[1]+mu) - lam2*mu2; + + if (b*b > a*d) + { + t = (-b + sqrt(b*b - a*d))/a; + x = pos[0] + t*ca; + y = pos[1] + t*sa; + + if ((t > margin)&&(y <= -mu)) + { + cval[nt] = 8; + tval[nt] = t; + x1[nt] = x; + y1[nt] = y; + tempconf[nt][0] = sval[8] + argument(-x/lambda, -(y+mu)/(1.0-mu)); + tempconf[nt][1] = argument(-lambda*(y+mu)/(1.0-mu), (1.0-mu)*x/lambda) - alpha; + nt++; + } + } + + /* intersection with right ellipse */ + a = 4.0*ca*ca + sa*sa*penrose_ratio*penrose_ratio; + b = 4.0*(pos[0] - cc)*ca + pos[1]*sa*penrose_ratio*penrose_ratio; + d = 4.0*(pos[0] - cc)*(pos[0] - cc) + (pos[1]*pos[1] - mu*mu)*penrose_ratio*penrose_ratio; + + if (b*b > a*d) + { + t = (-b + sqrt(b*b - a*d))/a; + x = pos[0] + t*ca; + y = pos[1] + t*sa; + + if ((t > margin)&&(x <= cc)) + { + cval[nt] = 12; + tval[nt] = t; + x1[nt] = x; + y1[nt] = y; + tempconf[nt][0] = sval[12] + argument(-y/mu, 2.0*(cc-x)/(mu*penrose_ratio)); + tempconf[nt][1] = argument(0.5*penrose_ratio*y, 2.0*(cc-x)/penrose_ratio) - alpha; + nt++; + } + + t = (-b - sqrt(b*b - a*d))/a; + x = pos[0] + t*ca; + y = pos[1] + t*sa; + + if ((t > margin)&&(x <= cc)) + { + cval[nt] = 12; + tval[nt] = t; + x1[nt] = x; + y1[nt] = y; + tempconf[nt][0] = sval[12] + argument(-y/mu, 2.0*(cc-x)/(mu*penrose_ratio)); + tempconf[nt][1] = argument(0.5*penrose_ratio*y, 2.0*(cc-x)/penrose_ratio) - alpha; + nt++; + } + + } + + /* intersection with left ellipse */ + b = 4.0*(pos[0] + cc)*ca + pos[1]*sa*penrose_ratio*penrose_ratio; + d = 4.0*(pos[0] + cc)*(pos[0] + cc) + (pos[1]*pos[1] - mu*mu)*penrose_ratio*penrose_ratio; + + if (b*b > a*d) + { + t = (-b + sqrt(b*b - a*d))/a; + x = pos[0] + t*ca; + y = pos[1] + t*sa; + + if ((t > margin)&&(x >= -cc)) + { + cval[nt] = 4; + tval[nt] = t; + x1[nt] = x; + y1[nt] = y; + tempconf[nt][0] = sval[4] + argument(y/mu, 2.0*(cc+x)/(mu*penrose_ratio)); + tempconf[nt][1] = argument(0.5*penrose_ratio*y, -2.0*(cc+x)/penrose_ratio) - alpha; + nt++; + } + + t = (-b - sqrt(b*b - a*d))/a; + x = pos[0] + t*ca; + y = pos[1] + t*sa; + + if ((t > margin)&&(x >= -cc)) + { + cval[nt] = 4; + tval[nt] = t; + x1[nt] = x; + y1[nt] = y; + tempconf[nt][0] = sval[4] + argument(y/mu, 2.0*(cc+x)/(mu*penrose_ratio)); + tempconf[nt][1] = argument(0.5*penrose_ratio*y, -2.0*(cc+x)/penrose_ratio) - alpha; + nt++; + } + + } + + /* rightmost vertical segments */ + if (ca > 0.0) + { + t = (lambda - pos[0])/ca; + x = lambda; + y = pos[1] + t*sa; + + if ((y >= width)&&(y < mu)) + { + cval[nt] = 15; + tval[nt] = t; + x1[nt] = x; + y1[nt] = y; + tempconf[nt][0] = sval[15] + y - width; + tempconf[nt][1] = PID - alpha; + nt++; + } + + else if ((y <= -width)&&(y > -mu)) + { + cval[nt] = 9; + tval[nt] = t; + x1[nt] = x; + y1[nt] = y; + tempconf[nt][0] = sval[9] + y + mu; + tempconf[nt][1] = PID - alpha; + nt++; + } + } + + /* leftmost vertical segments */ + if (ca < 0.0) + { + t = -(lambda + pos[0])/ca; + x = -lambda; + y = pos[1] + t*sa; + + if ((y >= width)&&(y < mu)) + { + cval[nt] = 1; + tval[nt] = t; + x1[nt] = x; + y1[nt] = y; + tempconf[nt][0] = sval[1] + mu - y; + tempconf[nt][1] = 3.0*PID - alpha; + nt++; + } + + else if ((y <= -width)&&(y > -mu)) + { + cval[nt] = 7; + tval[nt] = t; + x1[nt] = x; + y1[nt] = y; + tempconf[nt][0] = sval[7] - width - y; + tempconf[nt][1] = 3.0*PID - alpha; + nt++; + } + } + + /* vertical segments of right mushroom head */ + if (ca < 0.0) + { + t = (cc - pos[0])/ca; + x = cc; + y = pos[1] + t*sa; + + if ((t > 0.0)&&(y >= width)&&(y < mu)) + { + cval[nt] = 13; + tval[nt] = t; + x1[nt] = x; + y1[nt] = y; + tempconf[nt][0] = sval[13] - y + mu; + tempconf[nt][1] = 3.0*PID - alpha; + nt++; + } + + else if ((t > 0.0)&&(y <= -width)&&(y > -mu)) + { + cval[nt] = 11; + tval[nt] = t; + x1[nt] = x; + y1[nt] = y; + tempconf[nt][0] = sval[11] - y - width; + tempconf[nt][1] = 3.0*PID - alpha; + nt++; + } + } + + /* vertical segments of left mushroom head */ + if (ca > 0.0) + { + t = (-cc - pos[0])/ca; + x = -cc; + y = pos[1] + t*sa; + + if ((t > 0.0)&&(y >= width)&&(y < mu)) + { + cval[nt] = 3; + tval[nt] = t; + x1[nt] = x; + y1[nt] = y; + tempconf[nt][0] = sval[3] + y - width; + tempconf[nt][1] = PID - alpha; + nt++; + } + + else if ((t > 0.0)&&(y <= -width)&&(y > -mu)) + { + cval[nt] = 5; + tval[nt] = t; + x1[nt] = x; + y1[nt] = y; + tempconf[nt][0] = sval[5] + y + mu; + tempconf[nt][1] = PID - alpha; + nt++; + } + } + + /* upper horizontal segments */ + if (sa < 0.0) + { + t = (width - pos[1])/sa; + x = pos[0] + t*ca; + y = width; + + if ((t > 0.0)&&(x >= cc)&&(x < lambda)) + { + cval[nt] = 14; + tval[nt] = t; + x1[nt] = x; + y1[nt] = y; + tempconf[nt][0] = sval[14] + x - cc; + tempconf[nt][1] = - alpha; + nt++; + } + + else if ((t > 0.0)&&(x <= -cc)&&(x > -lambda)) + { + cval[nt] = 2; + tval[nt] = t; + x1[nt] = x; + y1[nt] = y; + tempconf[nt][0] = sval[2] + x + lambda; + tempconf[nt][1] = - alpha; + nt++; + } + } + + /* lower horizontal segments - TO BE CORRECTED */ + if (sa > 0.0) + { + t = (-width - pos[1])/sa; + x = pos[0] + t*ca; + y = -width; + + if ((t > 0.0)&&(x >= cc)&&(x < lambda)) + { + cval[nt] = 10; + tval[nt] = t; + x1[nt] = x; + y1[nt] = y; + tempconf[nt][0] = sval[10] - x + lambda; + tempconf[nt][1] = PI - alpha; + nt++; + } + + else if ((t > 0.0)&&(x <= -cc)&&(x > -lambda)) + { + cval[nt] = 6; + tval[nt] = t; + x1[nt] = x; + y1[nt] = y; + tempconf[nt][0] = sval[6] - x - cc; + tempconf[nt][1] = PI - alpha; + nt++; + } + } + /* find earliest intersection */ + tmin = tval[0]; + ntmin = 0; + for (i=1; i= ncircles) ncirc = ncircles - 1; + + angle = conf[0] - (double)ncirc*DPI; + + pos[0] = circles[ncirc].xc + circles[ncirc].radius*cos(angle); + pos[1] = circles[ncirc].yc + circles[ncirc].radius*sin(angle); + + *alpha = angle + PID + conf[1]; + + return(ncirc); + } + + + int vcircles_xy(double config[8], double alpha, double pos[2]) + /* determine initial configuration for start at point pos = (x,y) */ + { + double c0, s0, b, c, t, theta, delta, margin = 1.0e-12, tmin, rlarge = 1000.0; + double tval[ncircles], xint[ncircles], yint[ncircles], phiint[ncircles]; + int i, nt = 0, nscat[ncircles], ntmin; + + c0 = cos(alpha); + s0 = sin(alpha); + + for (i=0; i margin) /* there is an intersection with circle i */ + { + t = -b - sqrt(delta); + if (t > margin) + { + nscat[nt] = i; + + tval[nt] = t; + xint[nt] = pos[0] + t*c0; + yint[nt] = pos[1] + t*s0; + phiint[nt] = argument(xint[nt] - circles[i].xc, yint[nt] - circles[i].yc); + + nt++; + } + } + } + + if (nt > 0) /* there is at least one intersection */ + { + /* find earliest intersection */ + tmin = tval[0]; + ntmin = 0; + for (i=1; i= DPI) phiint[ntmin] -= DPI; + + config[0] = (double)nscat[ntmin]*DPI + phiint[ntmin]; + config[1] = PID - alpha + phiint[ntmin]; /* CHECK */ + if (config[1] < 0.0) config[1] += DPI; + if (config[1] >= PI) config[1] -= DPI; + + config[2] = 0.0; /* running time */ + config[3] = module2(xint[ntmin]-pos[0], yint[ntmin]-pos[1]); /* distance to collision */ + config[4] = pos[0]; /* start position */ + config[5] = pos[1]; + config[6] = xint[ntmin]; /* position of collision */ + config[7] = yint[ntmin]; + + /* set dummy coordinates if circles are absorbing */ + if (ABSORBING_CIRCLES) + { + config[0] = DUMMY_ABSORBING; + config[1] = PI; + } + + return(nscat[ntmin]); + } + else /* there is no intersection - set dummy values */ + { + config[0] = DUMMY_ABSORBING; + config[1] = PI; + config[2] = 0.0; + config[3] = rlarge; + config[4] = pos[0]; /* start position */ + config[5] = pos[1]; + config[6] = rlarge*cos(alpha); + config[7] = rlarge*sin(alpha); + + return(-1); + } + } + + int vcircles(double config[8]) + /* determine initial configuration when starting from boundary */ + { + double pos[2], alpha; + int c; + + c = pos_circles(config, pos, &alpha); + + c = vcircles_xy(config, alpha, pos); + + return(c); + } + +/****************************************************************************************/ +/* billiard with circular scatterers in a rectangle */ +/****************************************************************************************/ + + int pos_circles_in_rect(double conf[2], double pos[2], double *alpha) + /* determine position on boundary of circle */ + /* position varies between 0 and ncircles*2Pi for circles and between -BOUNDARY_SHIFT and 0 for boundary*/ + /* returns number of hit circle */ + { + double angle; + int ncirc, c; + + if (conf[0] >= 0) + { + ncirc = (int)(conf[0]/DPI); + if (ncirc >= ncircles) ncirc = ncircles - 1; + + angle = conf[0] - (double)ncirc*DPI; + + pos[0] = circles[ncirc].xc + circles[ncirc].radius*cos(angle); + pos[1] = circles[ncirc].yc + circles[ncirc].radius*sin(angle); + + *alpha = angle + PID + conf[1]; + + return(ncirc); + } + else /* particle starts on boundary */ + { +// conf[0] += 4.0*(lambda + 1.0); + conf[0] += BOUNDARY_SHIFT; + c = pos_rectangle(conf, pos, alpha); + +// conf[0] -= 4.0*(lambda + 1.0); + conf[0] -= BOUNDARY_SHIFT; + + return(-c-1); + } + } + + + int vcircles_in_rect_xy(double config[8], double alpha, double pos[2]) + /* determine initial configuration for start at point pos = (x,y) */ +// double config[8], alpha, pos[2]; + + { + double c0, s0, b, c, t, theta, delta, margin = 1.0e-12, tmin, rlarge = 1.0e10; + double tval[ncircles], xint[ncircles], yint[ncircles], phiint[ncircles]; + int i, nt = 0, nscat[ncircles], ntmin, side; + + c0 = cos(alpha); + s0 = sin(alpha); + + for (i=0; i margin) /* there is an intersection with circle i */ + { + t = -b - sqrt(delta); + if (t > margin) + { + nscat[nt] = i; + + tval[nt] = t; + xint[nt] = pos[0] + t*c0; + yint[nt] = pos[1] + t*s0; + phiint[nt] = argument(xint[nt] - circles[i].xc, yint[nt] - circles[i].yc); + + /* test wether intersection is in rectangle */ + if ((vabs(xint[nt]) < lambda)&&(vabs(yint[nt]) < 1.0)) nt++; + } + } + } + + if (nt > 0) /* there is at least one intersection */ + { + /* find earliest intersection */ + tmin = tval[0]; + ntmin = 0; + for (i=1; i= DPI) phiint[ntmin] -= DPI; + + config[0] = (double)nscat[ntmin]*DPI + phiint[ntmin]; + config[1] = PID - alpha + phiint[ntmin]; /* CHECK */ + if (config[1] < 0.0) config[1] += DPI; + if (config[1] >= PI) config[1] -= DPI; + + config[2] = 0.0; /* running time */ + config[3] = module2(xint[ntmin]-pos[0], yint[ntmin]-pos[1]); /* distance to collision */ + config[4] = pos[0]; /* start position */ + config[5] = pos[1]; + config[6] = xint[ntmin]; /* position of collision */ + config[7] = yint[ntmin]; + + + /* set dummy coordinates if circles are absorbing */ + if (ABSORBING_CIRCLES) + { + config[0] = DUMMY_ABSORBING; + config[1] = PI; +// return(DUMMY_SIDE_ABS); + } + + if (ABSORBING_CIRCLES) return(DUMMY_SIDE_ABS); + else return(nscat[ntmin]); + } + else /* there is no intersection with the circles - compute intersection with boundary */ + { + + side = vrectangle_xy(config, alpha, pos); + config[0] -= BOUNDARY_SHIFT; +// config[0] -= 4.0*(lambda+1.0); + +// printf("Hit side %i\n", side); +// print_config(config); + return(side - 5); + } + } + + int vcircles_in_rect(double config[8]) + /* determine initial configuration when starting from boundary */ + { + double pos[2], alpha; + int c; + + c = pos_circles_in_rect(config, pos, &alpha); + +// vcircles_in_rect_xy(config, alpha, pos); + c = vcircles_in_rect_xy(config, alpha, pos); + + return(c); + } + + +/****************************************************************************************/ +/* billiard with circular scatterers in a genus n surface (polygon with identifies opposite sides) */ +/****************************************************************************************/ + + int pos_circles_in_genusn(double conf[2], double pos[2], double *alpha) + /* determine position on boundary of circle */ + /* position varies between 0 and ncircles*2Pi for circles and between */ + /* BOUNDARY_SHIFT and 0 for boundary*/ + /* returns number of hit circle */ + { + double s, theta, omega, length, s1, angle, x, y; + int ncirc, c; + + if (conf[0] >= 0) + { + ncirc = (int)(conf[0]/DPI); + if (ncirc >= ncircles) ncirc = ncircles - 1; + + angle = conf[0] - (double)ncirc*DPI; + + pos[0] = circles[ncirc].xc + circles[ncirc].radius*cos(angle); + pos[1] = circles[ncirc].yc + circles[ncirc].radius*sin(angle); + + *alpha = angle + PID + conf[1]; + + return(ncirc); + } + else /* particle starts on boundary */ + { + omega = DPI/((double)NPOLY); + length = 2.0*sin(0.5*omega); + +// conf[0] += 2.0*length*(double)NPOLY; + conf[0] += BOUNDARY_SHIFT; + c = pos_genusn(conf, pos, alpha); + +// conf[0] -= 2.0*length*(double)NPOLY; + conf[0] -= BOUNDARY_SHIFT; + + return(-c); + } + } + + + int vcircles_in_genusn_xy(double config[8], double alpha, double pos[2]) + /* determine initial configuration for start at point pos = (x,y) */ + { + double c0, s0, b, c, t, theta, delta, margin = 1.0e-12, tmin, rlarge = 1.0e10, omega, length, angle, cw; + double tval[ncircles], xint[ncircles], yint[ncircles], phiint[ncircles]; + int i, k, nt = 0, nscat[ncircles], ntmin, side, condition; + + c0 = cos(alpha); + s0 = sin(alpha); + omega = DPI/((double)NPOLY); + length = 2.0*sin(0.5*omega); + cw = cos(omega*0.5); + + for (i=0; i margin) /* there is an intersection with circle i */ + { + t = -b - sqrt(delta); + if (t > margin) + { + nscat[nt] = i; + + tval[nt] = t; + xint[nt] = pos[0] + t*c0; + yint[nt] = pos[1] + t*s0; + phiint[nt] = argument(xint[nt] - circles[i].xc, yint[nt] - circles[i].yc); + + /* test wether intersection is in polygon */ + if (in_polygon(xint[nt], yint[nt], 1.0, NPOLY, APOLY)) nt++; + } + } + } + + if (nt > 0) /* there is at least one intersection */ + { + /* find earliest intersection */ + tmin = tval[0]; + ntmin = 0; + for (i=1; i= DPI) phiint[ntmin] -= DPI; + + config[0] = (double)nscat[ntmin]*DPI + phiint[ntmin]; + config[1] = PID - alpha + phiint[ntmin]; /* CHECK */ + if (config[1] < 0.0) config[1] += DPI; + if (config[1] >= PI) config[1] -= DPI; + + config[2] = 0.0; /* running time */ + config[3] = module2(xint[ntmin]-pos[0], yint[ntmin]-pos[1]); /* distance to collision */ + config[4] = pos[0]; /* start position */ + config[5] = pos[1]; + config[6] = xint[ntmin]; /* position of collision */ + config[7] = yint[ntmin]; + + + /* set dummy coordinates if circles are absorbing */ + if (ABSORBING_CIRCLES) + { + config[0] = DUMMY_ABSORBING; + config[1] = PI; + } + + return(nscat[ntmin]); + } + else /* there is no intersection with the circles - compute intersection with boundary */ + { + side = vgenusn_xy(config, alpha, pos); + +// config[0] -= 2.0*length*(double)NPOLY; + config[0] -= BOUNDARY_SHIFT; + + return(side); + } + } + + int vcircles_in_genusn(double config[8]) + /* determine initial configuration when starting from boundary */ + { + double pos[2], alpha; + int c; + + c = pos_circles_in_genusn(config, pos, &alpha); + + vcircles_in_genusn_xy(config, alpha, pos); + + return(c); + } + + +/****************************************************************************************/ +/* billiard with circular scatterers in a torus (rectangle with periodic boundary conditions) */ +/****************************************************************************************/ + + int pos_circles_in_torus(double conf[2], double pos[2], double *alpha) + /* determine position on boundary of circle */ + /* position varies between 0 and ncircles*2Pi for circles and between -BOUNDARY_SHIFT and 0 for boundary*/ + /* returns number of hit circle */ + { + double angle; + int ncirc, c; + + if (conf[0] >= 0) + { + ncirc = (int)(conf[0]/DPI); + if (ncirc >= ncircles) ncirc = ncircles - 1; + + angle = conf[0] - (double)ncirc*DPI; + + pos[0] = circles[ncirc].xc + circles[ncirc].radius*cos(angle); + pos[1] = circles[ncirc].yc + circles[ncirc].radius*sin(angle); + + *alpha = angle + PID + conf[1]; + + return(ncirc); + } + else /* particle starts on boundary */ + { +// conf[0] += 4.0*(lambda + 1.0); + conf[0] += BOUNDARY_SHIFT; + c = pos_rectangle(conf, pos, alpha); + +// conf[0] -= 4.0*(lambda + 1.0); + conf[0] -= BOUNDARY_SHIFT; + + return(-c-1); + } + } + + + int vcircles_in_torus_xy(double config[8], double alpha, double pos[2]) + /* determine initial configuration for start at point pos = (x,y) */ +// double config[8], alpha, pos[2]; + + { + double c0, s0, b, c, t, theta, delta, margin = 1.0e-12, tmin, rlarge = 1.0e10; + double tval[ncircles], xint[ncircles], yint[ncircles], phiint[ncircles]; + int i, nt = 0, nscat[ncircles], ntmin, side; + + c0 = cos(alpha); + s0 = sin(alpha); + + for (i=0; i margin) /* there is an intersection with circle i */ + { + t = -b - sqrt(delta); + if (t > margin) + { + nscat[nt] = i; + + tval[nt] = t; + xint[nt] = pos[0] + t*c0; + yint[nt] = pos[1] + t*s0; + phiint[nt] = argument(xint[nt] - circles[i].xc, yint[nt] - circles[i].yc); + + /* test wether intersection is in rectangle */ + if ((vabs(xint[nt]) < lambda)&&(vabs(yint[nt]) < 1.0)) nt++; + } + } + } + + if (nt > 0) /* there is at least one intersection */ + { + /* find earliest intersection */ + tmin = tval[0]; + ntmin = 0; + for (i=1; i= DPI) phiint[ntmin] -= DPI; + + config[0] = (double)nscat[ntmin]*DPI + phiint[ntmin]; + config[1] = PID - alpha + phiint[ntmin]; /* CHECK */ + if (config[1] < 0.0) config[1] += DPI; + if (config[1] >= PI) config[1] -= DPI; + + config[2] = 0.0; /* running time */ + config[3] = module2(xint[ntmin]-pos[0], yint[ntmin]-pos[1]); /* distance to collision */ + config[4] = pos[0]; /* start position */ + config[5] = pos[1]; + config[6] = xint[ntmin]; /* position of collision */ + config[7] = yint[ntmin]; + + + /* set dummy coordinates if circles are absorbing */ + if (ABSORBING_CIRCLES) + { + config[0] = DUMMY_ABSORBING; + config[1] = PI; + } + + return(nscat[ntmin]); + } + else /* there is no intersection with the circles - compute intersection with boundary */ + { + + side = vrectangle_xy(config, alpha, pos); + + if (config[0] < 2.0*lambda) config[0] = 4.0*lambda + 2.0 - config[0]; + else if (config[0] < 2.0*lambda + 2.0) config[0] = 6.0*lambda + 4.0 - config[0]; + else if (config[0] < 4.0*lambda + 2.0) config[0] = 4.0*lambda + 2.0 - config[0]; + else config[0] = 6.0*lambda + 4.0 - config[0]; + + config[0] -= BOUNDARY_SHIFT; + config[1] = PI - config[1]; +// config[0] -= 4.0*(lambda+1.0); + +// printf("Hit side %i\n", side); +// print_config(config); + return(side - 5); + } + } + + int vcircles_in_torus(double config[8]) + /* determine initial configuration when starting from boundary */ + { + double pos[2], alpha; + int c; + + c = pos_circles_in_torus(config, pos, &alpha); + + vcircles_in_torus_xy(config, alpha, pos); + + return(c); + } + + +/****************************************************************************************/ +/* billiard in poylgonal line */ +/****************************************************************************************/ + + int pos_polyline(double conf[2], double pos[2], double *alpha) + /* determine position on boundary of domain */ + /* position varies between 0 and nsides */ + /* returns number of hit setment */ + { + double len, rlarge = 1000.0; + int nside; + + nside = (int)conf[0]; + if (nside >= nsides) nside = nsides - 1; + + if (nside >= 0) /* position is on a side of the polygonal line */ + { + len = conf[0] - (double)nside; + + pos[0] = polyline[nside].x1 + len*(polyline[nside].x2 - polyline[nside].x1); + pos[1] = polyline[nside].y1 + len*(polyline[nside].y2 - polyline[nside].y1); + + *alpha = polyline[nside].angle + conf[1]; + return(nside); + } + else /* position is on an absorbing circle */ + { + pos[0] = rlarge; + pos[1] = 0.0; + + *alpha = 0.0; + return(-1); + } + } + + + int vpolyline_xy(double config[8], double alpha, double pos[2]) + /* determine initial configuration for start at point pos = (x,y) */ + { + double c0, s0, a, b, c, t, dx, delta, s, xi, yi, margin = 1.0e-12, tmin, rlarge = 1000.0; + double tval[nsides + ncircles], xint[nsides + ncircles], yint[nsides + ncircles], sint[nsides + ncircles]; + int i, nt = 0, nsegment[nsides + ncircles], ntmin; + + c0 = cos(alpha); + s0 = sin(alpha); + + for (i=0; i margin) /* there is an intersection with the line containing segment i */ + { + t = -(a*pos[0] + b*pos[1] + c)/delta; + if (t > margin) + { + xi = pos[0] + t*c0; + yi = pos[1] + t*s0; + dx = polyline[i].x2 - polyline[i].x1; + + if (vabs(dx) > margin) s = (xi - polyline[i].x1)/dx; + else s = (yi - polyline[i].y1)/(polyline[i].y2 - polyline[i].y1); +// printf("s = %.2f\n", s); + + if ((s >= 0.0)&&(s <= 1.0)) + /* the intersection is on the segment */ + { + nsegment[nt] = i; + tval[nt] = t; + sint[nt] = s; + xint[nt] = pos[0] + t*c0; + yint[nt] = pos[1] + t*s0; +// printf("s = %.2f, x = %.2f, y = %.2f\n", s, xint[nt], yint[nt]); + nt++; + } + } + } + } + + if (ABSORBING_CIRCLES) for (i=0; i margin) /* there is an intersection with circle i */ + { + t = -b - sqrt(delta); + if (t > margin) + { + nsegment[nt] = -1-i; + + tval[nt] = t; + xint[nt] = pos[0] + t*c0; + yint[nt] = pos[1] + t*s0; + sint[nt] = argument(xint[nt] - circles[i].xc, yint[nt] - circles[i].yc); + + nt++; + } + } + } + + if (nt > 0) /* there is at least one intersection */ + { + /* find earliest intersection */ + tmin = tval[0]; + ntmin = 0; + for (i=1; i= 0) + { + config[0] = (double)nsegment[ntmin] + sint[ntmin]; + config[1] = polyline[nsegment[ntmin]].angle - alpha; + if (config[1] < 0.0) config[1] += DPI; + if (config[1] >= PI) config[1] -= DPI; + } + /* set dummy coordinates if circles are absorbing */ + else if ((ABSORBING_CIRCLES)&&(nsegment[ntmin] < 0)) + { + config[0] = DUMMY_ABSORBING; + config[1] = PI; + } + config[2] = 0.0; /* running time */ + config[3] = module2(xint[ntmin]-pos[0], yint[ntmin]-pos[1]); /* distance to collision */ + config[4] = pos[0]; /* start position */ + config[5] = pos[1]; + config[6] = xint[ntmin]; /* position of collision */ + config[7] = yint[ntmin]; + + +// print_config(config); + + return(nsegment[ntmin]); + } + else /* there is no intersection - set dummy values */ + { + config[0] = DUMMY_ABSORBING; + config[1] = PI; + config[2] = 0.0; + config[3] = rlarge; + config[4] = pos[0]; /* start position */ + config[5] = pos[1]; + config[6] = rlarge*cos(alpha); + config[7] = rlarge*sin(alpha); + + return(-1); + } + } + + int vpolyline(double config[8]) + /* determine initial configuration when starting from boundary */ + { + double pos[2], alpha; + int c; + + c = pos_polyline(config, pos, &alpha); + + vpolyline_xy(config, alpha, pos); +// c = vpolyline_xy(config, alpha, pos); + + return(c); + } + +/****************************************************************************************/ +/* general billiard */ +/****************************************************************************************/ + + int pos_billiard(double conf[8], double pos[2], double *alpha) + /* determine initial configuration for start at point pos = (x,y) */ + { + switch (B_DOMAIN) { + case (D_RECTANGLE): + { + return(pos_rectangle(conf, pos, alpha)); + break; + } + case (D_ELLIPSE): + { + return(pos_ellipse(conf, pos, alpha)); + break; + } + case (D_STADIUM): + { + return(pos_stade(conf, pos, alpha)); + break; + } + case (D_SINAI): + { + return(pos_sinai(conf, pos, alpha)); + break; + } + case (D_TRIANGLE): + { + return(pos_triangle(conf, pos, alpha)); + break; + } + case (D_ANNULUS): + { + return(pos_annulus(conf, pos, alpha)); + break; + } + case (D_POLYGON): + { + return(pos_polygon(conf, pos, alpha)); + break; + } + case (D_REULEAUX): + { + return(pos_reuleaux(conf, pos, alpha)); + break; + } + case (D_FLOWER): + { + return(pos_flower(conf, pos, alpha)); + break; + } + case (D_ALT_REU): + { + return(pos_alt_reuleaux(conf, pos, alpha)); + break; + } + case (D_ANGLE): + { + return(pos_angle(conf, pos, alpha)); + break; + } + case (D_LSHAPE): + { + return(pos_lshape(conf, pos, alpha)); + break; + } + case (D_GENUSN): + { + return(pos_genusn(conf, pos, alpha)); + break; + } + case (D_PARABOLAS): + { + return(pos_parabolas(conf, pos, alpha)); + break; + } + case (D_PENROSE): + { + return(pos_penrose(conf, pos, alpha)); + break; + } + case (D_CIRCLES): + { + return(pos_circles(conf, pos, alpha)); + break; + } + case (D_CIRCLES_IN_RECT): + { + return(pos_circles_in_rect(conf, pos, alpha)); + break; + } + case (D_CIRCLES_IN_GENUSN): + { + return(pos_circles_in_genusn(conf, pos, alpha)); + break; + } + case (D_CIRCLES_IN_TORUS): + { + return(pos_circles_in_torus(conf, pos, alpha)); + break; + } + case (D_POLYLINE): + { + return(pos_polyline(conf, pos, alpha)); + break; + } + default: + { + printf("Function pos_billiard not defined for this billiard \n"); + } + } + } + + + + int vbilliard_xy(double config[8], double alpha, double pos[2]) + /* determine initial configuration for start at point pos = (x,y) */ + { + switch (B_DOMAIN) { + case (D_RECTANGLE): + { + return(vrectangle_xy(config, alpha, pos)); + break; + } + case (D_ELLIPSE): + { + return(vellipse_xy(config, alpha, pos)); + break; + } + case (D_STADIUM): + { + return(vstade_xy(config, alpha, pos)); + break; + } + case (D_SINAI): + { + return(vsinai_xy(config, alpha, pos)); + break; + } + case (D_TRIANGLE): + { + return(vtriangle_xy(config, alpha, pos)); + break; + } + case (D_ANNULUS): + { + return(vannulus_xy(config, alpha, pos)); + break; + } + case (D_POLYGON): + { + return(vpolygon_xy(config, alpha, pos)); + break; + } + case (D_REULEAUX): + { + return(vreuleaux_xy(config, alpha, pos)); + break; + } + case (D_FLOWER): + { + return(vflower_xy(config, alpha, pos)); + break; + } + case (D_ALT_REU): + { + return(valt_reuleaux_xy(config, alpha, pos)); + break; + } + case (D_ANGLE): + { + return(vangle_xy(config, alpha, pos)); + break; + } + case (D_LSHAPE): + { + return(vlshape_xy(config, alpha, pos)); + break; + } + case (D_GENUSN): + { + return(vgenusn_xy(config, alpha, pos)); + break; + } + case (D_PARABOLAS): + { + return(vparabolas_xy(config, alpha, pos)); + break; + } + case (D_PENROSE): + { + return(vpenrose_xy(config, alpha, pos)); + break; + } + case (D_CIRCLES): + { + return(vcircles_xy(config, alpha, pos)); + break; + } + case (D_CIRCLES_IN_RECT): + { + return(vcircles_in_rect_xy(config, alpha, pos)); + break; + } + case (D_CIRCLES_IN_GENUSN): + { + return(vcircles_in_genusn_xy(config, alpha, pos)); + break; + } + case (D_CIRCLES_IN_TORUS): + { + return(vcircles_in_torus_xy(config, alpha, pos)); + break; + } + case (D_POLYLINE): + { + return(vpolyline_xy(config, alpha, pos)); + break; + } + default: + { + printf("Function vbilliard_xy not defined for this billiard \n"); + } + } + } + + /* TO DO: fix returned value */ + + int vbilliard(double config[8]) + /* determine initial configuration when starting from boundary */ + { + double pos[2], theta, alpha; + int c; + + switch (B_DOMAIN) { + case (D_RECTANGLE): + { + c = pos_rectangle(config, pos, &alpha); + + return(vrectangle(config)); + break; + } + case (D_ELLIPSE): + { + c = pos_ellipse(config, pos, &alpha); + + return(vellipse(config)); + break; + } + case (D_STADIUM): + { + c = pos_stade(config, pos, &alpha); + + return(vstade(config)); + break; + } + case (D_SINAI): + { + c = pos_sinai(config, pos, &alpha); + + return(vsinai(config)); + break; + } + case (D_TRIANGLE): + { + c = pos_triangle(config, pos, &alpha); + + return(vtriangle(config)); + break; + } + case (D_ANNULUS): + { + c = pos_annulus(config, pos, &alpha); + + return(vannulus(config)); + break; + } + case (D_POLYGON): + { + c = pos_polygon(config, pos, &alpha); + + return(vpolygon(config)); + break; + } + case (D_REULEAUX): + { + c = pos_reuleaux(config, pos, &alpha); + + return(vreuleaux(config)); + break; + } + case (D_FLOWER): + { + c = pos_flower(config, pos, &alpha); + + return(vflower(config)); + break; + } + case (D_ALT_REU): + { + c = pos_alt_reuleaux(config, pos, &alpha); + + return(valt_reuleaux(config)); + break; + } + case (D_ANGLE): + { + c = pos_angle(config, pos, &alpha); + + return(vangle(config)); + break; + } + case (D_LSHAPE): + { + c = pos_lshape(config, pos, &alpha); + + return(vlshape(config)); + break; + } + case (D_GENUSN): + { + c = pos_genusn(config, pos, &alpha); + + return(vgenusn(config)); + break; + } + case (D_PARABOLAS): + { + c = pos_parabolas(config, pos, &alpha); + + return(vparabolas(config)); + break; + } + case (D_PENROSE): + { + c = pos_penrose(config, pos, &alpha); + + return(vpenrose(config)); + break; + } + case (D_CIRCLES): + { + c = pos_circles(config, pos, &alpha); + + return(vcircles(config)); + break; + } + case (D_CIRCLES_IN_RECT): + { + c = pos_circles_in_rect(config, pos, &alpha); + + return(vcircles_in_rect(config)); + break; + } + case (D_CIRCLES_IN_GENUSN): + { + c = pos_circles_in_genusn(config, pos, &alpha); + + return(vcircles_in_genusn(config)); + break; + } + case (D_CIRCLES_IN_TORUS): + { + c = pos_circles_in_torus(config, pos, &alpha); + + return(vcircles_in_torus(config)); + break; + } + case (D_POLYLINE): + { + c = pos_polyline(config, pos, &alpha); + + return(vpolyline(config)); + break; + } + default: + { + printf("Function vbilliard not defined for this billiard \n"); + } + } + } + + int xy_in_billiard(double x, double y) + /* returns 1 if (x,y) represents a point in the billiard */ + { + double l2, r1, r2, omega, omega2, c, angle, x1, y1, x2, co, so, x2plus, x2minus, width; + int condition, k; + + switch (B_DOMAIN) { + case D_RECTANGLE: + { + if ((vabs(x) -0.5*lambda)&&(x < 0.5*lambda)&&(y > -1.0)&&(y < 1.0)) return(1); + else if (module2(x+0.5*lambda, y) < 1.0) return(1); + else if (module2(x-0.5*lambda, y) < 1.0) return(1); + else return(0); + break; + } + case D_SINAI: + { + if (x*x + y*y > lambda*lambda) return(1); + else return(0); + break; + } + case D_DIAMOND: + { + l2 = lambda*lambda; + r2 = l2 + (lambda-1.0)*(lambda-1.0); + if ((x*x + y*y < 1.0)&&((x-lambda)*(x-lambda) + (y-lambda)*(y-lambda) > r2) + &&((x-lambda)*(x-lambda) + (y+lambda)*(y+lambda) > r2) + &&((x+lambda)*(x+lambda) + (y-lambda)*(y-lambda) > r2) + &&((x+lambda)*(x+lambda) + (y+lambda)*(y+lambda) > r2)) return(1); + else return(0); + break; + } + case D_TRIANGLE: + { + if ((x>-lambda)&&(y>-1.0)&&(lambda*y+x<0.0)) return(1); + else return(0); + break; + } + case D_ANNULUS: + { + l2 = lambda*lambda; + r1 = x*x + y*y; + r2 = (x-mu)*(x-mu) + y*y; + if ((r2 > l2)&&(r1 < 1.0)) return(1); + else return(0); + break; + } + case D_POLYGON: + { + return(in_polygon(x, y, 1.0, NPOLY, APOLY)); + break; + } + case D_REULEAUX: + { + condition = 1; + omega2 = PI/((double)NPOLY); + co = cos(omega2); + so = sin(omega2); + if (lambda > 0.0) x2 = co + sqrt(lambda*lambda - so*so); + else x2 = co - sqrt(lambda*lambda - so*so); + + for (k=0; k 0.0) condition = condition*((x1-x2)*(x1-x2) + y1*y1 > lambda*lambda); + else condition = condition*((x1-x2)*(x1-x2) + y1*y1 < lambda*lambda); + } + return(condition); + break; + } + /* D_REULEAUX : distance to all centers of arcs should be larger than lambda */ + case D_FLOWER: + { + /* TO DO */ + return(1); + break; + } + case D_ALT_REU: + { + condition = 1; + omega2 = PI/((double)NPOLY); + co = cos(omega2); + so = sin(omega2); + x2plus = co + sqrt(lambda*lambda - so*so); + x2minus = co - sqrt(lambda*lambda - so*so); + + for (k=0; k lambda*lambda); + else condition = condition*((x1-x2minus)*(x1-x2minus) + y1*y1 < lambda*lambda); + } + return(condition); + break; + } + case D_ANGLE: + { + if (y <= 0.0) return(0); + else if (x*sin(lambda*PI) < -y*cos(lambda*PI)) return(1); + else return(0); + break; + } + case D_LSHAPE: + { + if ((y <= 0.0)&&(x >= -1.0)&&(x <= 1.0)) return(1); + else if ((x >= -1.0)&&(x <= 0.0)) return(1); + else return(0); + break; + } + case D_GENUSN: /* same as polygon */ + { + return(in_polygon(x, y, 1.0, NPOLY, APOLY)); + break; + } + case D_PARABOLAS: + { + condition = 1; + omega = DPI/((double)NPOLY); + for (k=0; k= lambda) return(0); + /* upper and lower ellipse */ + else if ((vabs(y) >= mu)&&(x*x/(lambda*lambda) + (y1-mu)*(y1-mu)/((1.0-mu)*(1.0-mu)) >= 1.0)) return(0); + /* small ellipses */ + else if ((vabs(x) <= c)&&(4.0*(x1-c)*(x1-c)/(mu*mu*penrose_ratio*penrose_ratio) + y*y/(mu*mu) <= 1.0)) return(0); + /* straight parts */ + else if ((vabs(x) >= c)&&(vabs(y) <= width)) return(0); + else return(1); + } + case D_CIRCLES: + { + condition = 1; + for (k=0; k= lambda)||(vabs(y) >= 1.0)) return(0); + else + { + condition = 1; + for (k=0; k= lambda)||(vabs(y) >= 1.0)) return(0); + else + { + condition = 1; + for (k=0; k YMAX) y -= height; + circles[n].radius = mu; + circles[n].active = 1; + } + + break; + } + case (C_GOLDEN_SPIRAL): + { + ncircles = 1; + circles[0].xc = 0.0; + circles[0].yc = 0.0; + + gamma = (sqrt(5.0) - 1.0)*PI; /* golden mean times 2Pi */ + phi = 0.0; + r0 = 2.0*mu; + r = r0 + mu; + + for (i=0; i YMIN - mu)) circles[i].active = 1; + } + break; + } + case (C_RAND_DISPLACED): + { + ncircles = NCX*NCY; + dy = (YMAX - YMIN)/((double)NCY); + for (i = 0; i < NCX; i++) + for (j = 0; j < NCY; j++) + { + n = NCY*i + j; + circles[n].xc = ((double)(i-NCX/2) + 0.5*((double)rand()/RAND_MAX - 0.5))*dy; + circles[n].yc = YMIN + ((double)j + 0.5 + 0.5*((double)rand()/RAND_MAX - 0.5))*dy; + circles[n].radius = mu*sqrt(1.0 + 0.8*((double)rand()/RAND_MAX - 0.5)); + circles[n].active = 1; + } + break; + } + case (C_RAND_POISSON): + { + ncircles = NPOISSON; + for (n = 0; n < NPOISSON; n++) + { + circles[n].xc = lambda*(2.0*(double)rand()/RAND_MAX - 1.0); + circles[n].yc = (YMAX - YMIN)*(double)rand()/RAND_MAX + YMIN; + circles[n].radius = mu; + circles[n].active = 1; + } + break; + } + case (C_POISSON_DISC): + { + printf("Generating Poisson disc sample\n"); + /* generate first circle */ + circles[0].xc = lambda*(2.0*(double)rand()/RAND_MAX - 1.0); + circles[0].yc = (YMAX - YMIN)*(double)rand()/RAND_MAX + YMIN; + active_poisson[0] = 1; + n_p_active = 1; + ncircles = 1; + + while ((n_p_active > 0)&&(ncircles < NMAXCIRCLES)) + { + /* randomly select an active circle */ + i = rand()%(ncircles); + while (!active_poisson[i]) i = rand()%(ncircles); +// printf("Starting from circle %i at (%.3f,%.3f)\n", i, circlex[i], circley[i]); + /* generate new candidates */ + naccepted = 0; + for (j=0; j= dpoisson*dpoisson); + /* new circle is in domain */ + far = far*(vabs(x) < lambda)*(y < YMAX)*(y > YMIN); + } + if (far) /* accept new circle */ + { + printf("New circle at (%.3f,%.3f) accepted\n", x, y); + circles[ncircles].xc = x; + circles[ncircles].xc = y; + circles[ncircles].radius = mu; + circles[ncircles].active = 1; + active_poisson[ncircles] = 1; + ncircles++; + n_p_active++; + naccepted++; + } +// else printf("Rejected\n"); + } + if (naccepted == 0) /* inactivate circle i */ + { +// printf("No candidates work, inactivate circle %i\n", i); + active_poisson[i] = 0; + n_p_active--; + } + printf("%i active circles\n", n_p_active); + } + + printf("Generated %i circles\n", ncircles); + break; + } + + case (C_LASER): + { + ncircles = 17; + + xx[0] = 0.5*(x_shooter + x_target); + xx[1] = lambda - 0.5*(x_target - x_shooter); + if (xx[1] > lambda) xx[1] = 2.0*lambda - xx[1]; + xx[2] = -xx[0]; + xx[3] = -xx[1]; + + yy[0] = 0.5*(y_shooter + y_target); + yy[1] = 1.0 - 0.5*(y_target - y_shooter); + if (yy[1] > 1.0) yy[1] = 2.0 - yy[1]; + yy[2] = -yy[0]; + yy[3] = -yy[1]; + + for (i = 0; i < 4; i++) + for (j = 0; j < 4; j++) + { + circles[4*i + j].xc = xx[i]; + circles[4*i + j].yc = yy[j]; + + } + + circles[ncircles - 1].xc = x_target; + circles[ncircles - 1].yc = y_target; + + for (i=0; ix1 = z.x1 - 2.0*zperp.x1; + zprime->y1 = z.y1 - 2.0*zperp.y1; + + return(1); +} + +void init_polyline(t_segment polyline[NMAXPOLY], t_circle circles[NMAXCIRCLES]) +{ + int i, j, k, l, n, z, ii, jj, terni[SDEPTH], ternj[SDEPTH], quater[SDEPTH], cond; + short int vkoch[NMAXCIRCLES], turnright; + double ratio, omega, angle, sw, length, dist, x, y, ta, tb, a, b; + + switch (POLYLINE_PATTERN) { + case (P_RECTANGLE): + { + add_rectangle_to_polyline(0.0, 0.0, 2.0*lambda, 2.0, polyline, circles); + break; + } + case (P_TOKARSKY): + { + nsides = 26; + ncircles = 26; + + polyline[0].x1 = 0.0; polyline[0].y1 = 2.0; polyline[0].angle = 0.0; + polyline[1].x1 = 1.0; polyline[1].y1 = 2.0; polyline[1].angle = -PID; + polyline[2].x1 = 1.0; polyline[2].y1 = 1.0; polyline[2].angle = 0.0; + polyline[3].x1 = 2.0; polyline[3].y1 = 1.0; polyline[3].angle = -PID; + polyline[4].x1 = 2.0; polyline[4].y1 = 0.0; polyline[4].angle = 0.5*PID; + polyline[5].x1 = 3.0; polyline[5].y1 = 1.0; polyline[5].angle = 0.0; + polyline[6].x1 = 4.0; polyline[6].y1 = 1.0; polyline[6].angle = -PID; + polyline[7].x1 = 4.0; polyline[7].y1 = 0.0; polyline[7].angle = 0.5*PID; + polyline[8].x1 = 5.0; polyline[8].y1 = 1.0; polyline[8].angle = 0.0; + polyline[9].x1 = 6.0; polyline[9].y1 = 1.0; polyline[9].angle = -PID; + polyline[10].x1 = 6.0; polyline[10].y1 = 0.0; polyline[10].angle = 0.5*PID; + polyline[11].x1 = 7.0; polyline[11].y1 = 1.0; polyline[11].angle = PID; + polyline[12].x1 = 7.0; polyline[12].y1 = 2.0; polyline[12].angle = 0.0; + polyline[13].x1 = 8.0; polyline[13].y1 = 2.0; polyline[13].angle = PID; + polyline[14].x1 = 8.0; polyline[14].y1 = 3.0; polyline[14].angle = PI; + polyline[15].x1 = 7.0; polyline[15].y1 = 3.0; polyline[15].angle = 1.5*PID; + polyline[16].x1 = 6.0; polyline[16].y1 = 4.0; polyline[16].angle = -PID; + polyline[17].x1 = 6.0; polyline[17].y1 = 3.0; polyline[17].angle = PI; + polyline[18].x1 = 5.0; polyline[18].y1 = 3.0; polyline[18].angle = -PID; + polyline[19].x1 = 5.0; polyline[19].y1 = 2.0; polyline[19].angle = PI; + polyline[20].x1 = 3.0; polyline[20].y1 = 2.0; polyline[20].angle = PID; + polyline[21].x1 = 3.0; polyline[21].y1 = 3.0; polyline[21].angle = PI; + polyline[22].x1 = 2.0; polyline[22].y1 = 3.0; polyline[22].angle = PID; + polyline[23].x1 = 2.0; polyline[23].y1 = 4.0; polyline[23].angle = PI; + polyline[24].x1 = 1.0; polyline[24].y1 = 4.0; polyline[24].angle = -PID; + polyline[25].x1 = 1.0; polyline[25].y1 = 3.0; polyline[25].angle = 2.5*PID; + + ratio = (XMAX - XMIN)/8.4; + for (i=0; i NMAXPOLY) + { + printf("NMAXPOLY has to be increased to at least %i\n", nsides); + nsides = NMAXPOLY; + } + + for (i=0; i DPI) angle -= DPI; + while (angle < 0.0) angle += DPI; + } + + for (i=0; i NMAXPOLY) + { + printf("NMAXPOLY has to be increased to at least %i\n", nsides); + nsides = NMAXPOLY; + } + + for (i=0; i DPI) angle -= DPI; + while (angle < 0.0) angle += DPI; + } + + for (i=0; i