/*********************/ /* Graphics routines */ /*********************/ #include "colors_waves.c" #define TIFF_FREE_PERIOD 1 int writetiff_new(char *filename, char *description, int x, int y, int width, int height, int compression) { TIFF *file; GLubyte *image, *p; int i; file = TIFFOpen(filename, "w"); if (file == NULL) { return 1; } image = (GLubyte *) malloc(width * height * sizeof(GLubyte) * 3); /* OpenGL's default 4 byte pack alignment would leave extra bytes at the end of each image row so that each full row contained a number of bytes divisible by 4. Ie, an RGB row with 3 pixels and 8-bit componets would be laid out like "RGBRGBRGBxxx" where the last three "xxx" bytes exist just to pad the row out to 12 bytes (12 is divisible by 4). To make sure the rows are packed as tight as possible (no row padding), set the pack alignment to 1. */ glPixelStorei(GL_PACK_ALIGNMENT, 1); glReadPixels(x, y, width, height, GL_RGB, GL_UNSIGNED_BYTE, image); TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32_t) width); TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32_t) height); TIFFSetField(file, TIFFTAG_BITSPERSAMPLE, 8); TIFFSetField(file, TIFFTAG_COMPRESSION, compression); TIFFSetField(file, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB); TIFFSetField(file, TIFFTAG_ORIENTATION, ORIENTATION_BOTLEFT); TIFFSetField(file, TIFFTAG_SAMPLESPERPIXEL, 3); TIFFSetField(file, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG); TIFFSetField(file, TIFFTAG_ROWSPERSTRIP, 1); TIFFSetField(file, TIFFTAG_IMAGEDESCRIPTION, description); p = image; for (i = height - 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); /* prenvents RAM consumption*/ TIFFClose(file); return 0; } int writetiff(char *filename, char *description, int x, int y, int width, int height, int compression) { TIFF *file; GLubyte *image, *p; int i; static int counter = 0; file = TIFFOpen(filename, "w"); if (file == NULL) { return 1; } image = (GLubyte *) malloc(width * height * sizeof(GLubyte) * 3); /* OpenGL's default 4 byte pack alignment would leave extra bytes at the end of each image row so that each full row contained a number of bytes divisible by 4. Ie, an RGB row with 3 pixels and 8-bit componets would be laid out like "RGBRGBRGBxxx" where the last three "xxx" bytes exist just to pad the row out to 12 bytes (12 is divisible by 4). To make sure the rows are packed as tight as possible (no row padding), set the pack alignment to 1. */ glPixelStorei(GL_PACK_ALIGNMENT, 1); glReadPixels(x, y, width, height, GL_RGB, GL_UNSIGNED_BYTE, image); TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32_t) width); TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32_t) height); TIFFSetField(file, TIFFTAG_BITSPERSAMPLE, 8); TIFFSetField(file, TIFFTAG_COMPRESSION, compression); TIFFSetField(file, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB); TIFFSetField(file, TIFFTAG_ORIENTATION, ORIENTATION_BOTLEFT); TIFFSetField(file, TIFFTAG_SAMPLESPERPIXEL, 3); TIFFSetField(file, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG); TIFFSetField(file, TIFFTAG_ROWSPERSTRIP, 1); TIFFSetField(file, TIFFTAG_IMAGEDESCRIPTION, description); p = image; for (i = height - 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; } /* added 9/9/22 and removed again, since it produces an unwanted "band" on the right */ /* readded 5/11/22 */ if (SAVE_MEMORY) free(image); /* prevents RAM consumption*/ // { // counter++; // if (counter%TIFF_FREE_PERIOD == 0) // { // free(image); /* prevents RAM consumption*/ // counter = 0; // } // } TIFFClose(file); 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); glOrtho(0.0, NX, 0.0, NY, -1.0, 1.0); } void blank() { 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 test_save_frame() /* some tests with various resolutions */ { static int counter = 0; char *name="wave.", 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, "Wave equation in a planar domain", 0, 0, WINWIDTH, WINHEIGHT, COMPRESSION_LZW); // works for 1080p -> "-50px" // choose one of the following according to the comment beside. // writetiff(n2, "Wave equation in a planar domain", 0, 0, WINWIDTH, WINHEIGHT-40, COMPRESSION_LZW); /* to use with 1080p in drop_billiard.c- probably the best because it's // generating 1080p image, lighter, and then cropping those 40 pixels to // avoid the strange band*/ // writetiff(n2, "Wave equation in a planar domain", 0, 0, WINWIDTH, WINHEIGHT-50, COMPRESSION_LZW); // works for 1080p -> "-50px" band!!! // writetiff(n2, "Wave equation in a planar domain", 0, 0, 1920, 1080-40, COMPRESSION_LZW); //another perfect 1080p from 1440p in setup // writetiff(n2, "Wave equation in a planar domain", -WINWIDTH/8+320, -WINHEIGHT/8+180, WINWIDTH-640, WINHEIGHT-400, COMPRESSION_LZW); // perfect 1040p from 1440p in setup } void test_save_frame_counter(int counter) /* some tests with various resolutions */ /* same as save_frame, but with imposed image number (for option DOUBLE_MOVIE) */ { char *name="wave.", n2[100]; char format[6]=".%05i"; strcpy(n2, name); sprintf(strstr(n2,"."), format, counter); strcat(n2, ".tif"); printf(" saving frame %s \n",n2); writetiff(n2, "Wave equation in a planar domain", 0, 0, WINWIDTH, WINHEIGHT, COMPRESSION_LZW); // works for 1080p -> "-50px" // choose one of the following according to the comment beside. // writetiff(n2, "Wave equation in a planar domain", 0, 0, WINWIDTH, WINHEIGHT-40, COMPRESSION_LZW); /* to use with 1080p in drop_billiard.c- probably the best because it's // generating 1080p image, lighter, and then cropping those 40 pixels to // avoid the strange band*/ // writetiff(n2, "Wave equation in a planar domain", 0, 0, WINWIDTH, WINHEIGHT-50, COMPRESSION_LZW); // works for 1080p -> "-50px" band!!! // writetiff(n2, "Wave equation in a planar domain", 0, 0, 1920, 1080-40, COMPRESSION_LZW); //another perfect 1080p from 1440p in setup // writetiff(n2, "BWave equation in a planar domain", -WINWIDTH/8+320, -WINHEIGHT/8+180, WINWIDTH-640, WINHEIGHT-400, COMPRESSION_LZW); // perfect 1040p from 1440p in setup } void save_frame() { static int counter = 0; char *name="wave.", 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, "Wave equation in a planar domain", 0, 0, WINWIDTH, WINHEIGHT, COMPRESSION_LZW); } void save_frame_counter(int counter) /* same as save_frame, but with imposed image number (for option DOUBLE_MOVIE) */ { char *name="wave.", n2[100]; char format[6]=".%05i"; strcpy(n2, name); sprintf(strstr(n2,"."), format, counter); strcat(n2, ".tif"); printf(" saving frame %s \n",n2); writetiff(n2, "Wave equation in a planar domain", 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 } } /*********************/ /* 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; } return(alph); } double iabs(int i) /* absolute value */ { int res; if (i<0) res = -i; else res = i; return(i); } double gaussian() /* returns standard normal random variable, using Box-Mueller algorithm */ { static double V1, V2, S; static int phase = 0; double X; if (phase == 0) { do { double U1 = (double)rand() / RAND_MAX; double U2 = (double)rand() / RAND_MAX; V1 = 2 * U1 - 1; V2 = 2 * U2 - 1; S = V1 * V1 + V2 * V2; } while(S >= 1 || S == 0); X = V1 * sqrt(-2 * log(S) / S); } else X = V2 * sqrt(-2 * log(S) / S); phase = 1 - phase; return X; } // 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= 1.0) return(0); omega = DPI/((double)polygon.nsides); cw = cos(omega*0.5); for (k=0; k NX-1) ij[0] = NX-1; if (ij[1] < 0) ij[1] = 0; if (ij[1] > NY-1) ij[1] = NY-1; } void xy_to_pos(double x, double y, double pos[2]) /* convert (x,y) position to double-valued position in table representing wave */ { double x1, y1; x1 = (x - XMIN)/(XMAX - XMIN); y1 = (y - YMIN)/(YMAX - YMIN); pos[0] = x1 * (double)NX; pos[1] = y1 * (double)NY; } void ij_to_xy(int i, int j, double xy[2]) /* convert (i,j) position in table representing wave to (x,y) */ { double x1, y1; xy[0] = XMIN + ((double)i)*(XMAX-XMIN)/((double)NX); xy[1] = YMIN + ((double)j)*(YMAX-YMIN)/((double)NY); } void erase_area(double x, double y, double dx, double dy) { double pos[2], rgb[3]; hsl_to_rgb(220.0, 0.8, 0.7, rgb); glColor3f(rgb[0], rgb[1], rgb[2]); glBegin(GL_QUADS); xy_to_pos(x - dx, y - dy, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(x + dx, y - dy, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(x + dx, y + dy, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(x - dx, y + dy, pos); glVertex2d(pos[0], pos[1]); glEnd(); } void erase_area_ij(int imin, int jmin, int imax, int jmax) { glColor3f(0.0, 0.0, 0.0); glBegin(GL_QUADS); glVertex2i(imin, jmin); glVertex2i(imax, jmin); glVertex2i(imax, jmax); glVertex2i(imin, jmax); glEnd(); } void erase_area_rgb(double x, double y, double dx, double dy, double rgb[3]) { double pos[2]; glColor3f(rgb[0], rgb[1], rgb[2]); glBegin(GL_QUADS); xy_to_pos(x - dx, y - dy, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(x + dx, y - dy, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(x + dx, y + dy, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(x - dx, y + dy, pos); glVertex2d(pos[0], pos[1]); glEnd(); } void erase_area_hsl(double x, double y, double dx, double dy, double h, double s, double l) { double pos[2], rgb[3]; hsl_to_rgb(h, s, l, rgb); erase_area_rgb(x, y, dx, dy, rgb); } void draw_line(double x1, double y1, double x2, double y2) { double pos[2]; glBegin(GL_LINE_STRIP); xy_to_pos(x1, y1, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(x2, y2, pos); glVertex2d(pos[0], pos[1]); glEnd(); } void draw_rectangle(double x1, double y1, double x2, double y2) { double pos[2]; glBegin(GL_LINE_LOOP); xy_to_pos(x1, y1, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(x2, y1, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(x2, y2, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(x1, y2, pos); glVertex2d(pos[0], pos[1]); glEnd(); } void draw_filled_rectangle(double x1, double y1, double x2, double y2) { double pos[2]; glBegin(GL_TRIANGLE_FAN); xy_to_pos(x1, y1, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(x2, y1, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(x2, y2, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(x1, y2, pos); glVertex2d(pos[0], pos[1]); glEnd(); } void draw_rotated_rectangle(double x1, double y1, double x2, double y2) { double pos[2]; double xa, ya, xb, yb, xc, yc; xa = 0.5*(x1 - y2); xb = 0.5*(x2 - y1); xc = 0.5*(x1 - y1); ya = 0.5*(x1 + y1); yb = 0.5*(x2 + y2); yc = 0.5*(x2 + y1); glBegin(GL_LINE_LOOP); xy_to_pos(xc, ya, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(xb, yc, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(xc, yb, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(xa, yc, pos); glVertex2d(pos[0], pos[1]); glEnd(); } void draw_circle(double x, double y, double r, int nseg) { int i; double pos[2], alpha, dalpha; dalpha = DPI/(double)nseg; glBegin(GL_LINE_LOOP); for (i=0; i<=nseg; i++) { alpha = (double)i*dalpha; xy_to_pos(x + r*cos(alpha), y + r*sin(alpha), pos); glVertex2d(pos[0], pos[1]); } glEnd(); } void draw_circle_arc(double x, double y, double r, double angle1, double dangle, int nseg) { int i; double pos[2], alpha, dalpha; dalpha = dangle/(double)nseg; glBegin(GL_LINE_STRIP); for (i=0; i<=nseg; i++) { alpha = angle1 + (double)i*dalpha; xy_to_pos(x + r*cos(alpha), y + r*sin(alpha), pos); glVertex2d(pos[0], pos[1]); } glEnd(); } void draw_colored_circle(double x, double y, double r, int nseg, double rgb[3]) { int i; double pos[2], alpha, dalpha; dalpha = DPI/(double)nseg; glColor3f(rgb[0], rgb[1], rgb[2]); glBegin(GL_TRIANGLE_FAN); xy_to_pos(x, y, pos); glVertex2d(pos[0], pos[1]); for (i=0; i<=nseg; i++) { alpha = (double)i*dalpha; xy_to_pos(x + r*cos(alpha), y + r*sin(alpha), pos); glVertex2d(pos[0], pos[1]); } glEnd(); } void draw_tpolygon(t_polygon polygon) { int i; double pos[2], alpha, dalpha; dalpha = DPI/(double)polygon.nsides; glBegin(GL_LINE_LOOP); for (i=0; i<=polygon.nsides; i++) { alpha = PID*polygon.angle + (double)i*dalpha; xy_to_pos(polygon.xc + polygon.radius*cos(alpha), polygon.yc + polygon.radius*sin(alpha), pos); glVertex2d(pos[0], pos[1]); } glEnd(); } double michelson_schedule(int i) { double t; t = (double)(i - INITIAL_TIME)/(double)NSTEPS; return(4.0*t*WALL_WIDTH); // return(2.0*t*WALL_WIDTH); } int init_circle_config_pattern(t_circle circles[NMAXCIRCLES], int circle_pattern) /* initialise the arrays circlex, circley, circlerad and circleactive */ /* for billiard shape D_CIRCLES */ { int i, j, k, n, ncirc0, n_p_active, ncandidates=5000, naccepted; double dx, dy, p, phi, r, r0, ra[5], sa[5], height, x, y = 0.0, gamma, dpoisson = 3.25*MU, xx[4], yy[4], dr, dphi; short int active_poisson[NMAXCIRCLES], far; switch (circle_pattern) { case (C_SQUARE): { ncircles = NGRIDX*NGRIDY; dy = (YMAX - YMIN)/((double)NGRIDY); for (i = 0; i < NGRIDX; i++) for (j = 0; j < NGRIDY; j++) { n = NGRIDY*i + j; circles[n].xc = ((double)(i-NGRIDX/2) + 0.5)*dy; circles[n].yc = YMIN + ((double)j + 0.5)*dy; circles[n].radius = MU; circles[n].active = 1; } break; } case (C_HEX): { ncircles = NGRIDX*(NGRIDY+1); dy = (YMAX - YMIN)/((double)NGRIDY); dx = dy*0.5*sqrt(3.0); for (i = 0; i < NGRIDX; i++) for (j = 0; j < NGRIDY+1; j++) { n = (NGRIDY+1)*i + j; circles[n].xc = ((double)(i-NGRIDX/2) + 0.5)*dy; /* is +0.5 needed? */ circles[n].yc = YMIN + ((double)j - 0.5)*dy; if ((i+NGRIDX)%2 == 1) circles[n].yc += 0.5*dy; circles[n].radius = MU; /* activate only circles that intersect the domain */ if ((circles[n].yc < YMAX + MU)&&(circles[n].yc > YMIN - MU)) circles[n].active = 1; else circles[n].active = 0; } break; } case (C_RAND_DISPLACED): { ncircles = NGRIDX*NGRIDY; dy = (YMAX - YMIN)/((double)NGRIDY); for (i = 0; i < NGRIDX; i++) for (j = 0; j < NGRIDY; j++) { n = NGRIDY*i + j; circles[n].xc = ((double)(i-NGRIDX/2) + 0.5 + 0.5*((double)rand()/RAND_MAX - 0.5))*dy; circles[n].yc = YMIN + 0.5 + ((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_PERCOL): { ncircles = NGRIDX*NGRIDY; dy = (YMAX - YMIN)/((double)NGRIDY); for (i = 0; i < NGRIDX; i++) for (j = 0; j < NGRIDY; j++) { n = NGRIDY*i + j; circles[n].xc = ((double)(i-NGRIDX/2) + 0.5)*dy; circles[n].yc = YMIN + ((double)j + 0.5)*dy; circles[n].radius = MU; p = (double)rand()/RAND_MAX; if (p < P_PERCOL) circles[n].active = 1; else circles[n].active = 0; } 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_CLOAK): { ncircles = 200; for (i = 0; i < 40; i++) for (j = 0; j < 5; j++) { n = 5*i + j; phi = (double)i*DPI/40.0; r = LAMBDA*0.5*(1.0 + (double)j/5.0); circles[n].xc = r*cos(phi); circles[n].yc = r*sin(phi); circles[n].radius = MU; circles[n].active = 1; } break; } case (C_CLOAK_A): /* optimized model A1 by C. Jo et al */ { ncircles = 200; ra[0] = 0.0731; sa[0] = 1.115; ra[1] = 0.0768; sa[1] = 1.292; ra[2] = 0.0652; sa[2] = 1.464; ra[3] = 0.056; sa[3] = 1.633; ra[4] = 0.0375; sa[4] = 1.794; for (i = 0; i < 40; i++) for (j = 0; j < 5; j++) { n = 5*i + j; phi = (double)i*DPI/40.0; r = LAMBDA*sa[j]; circles[n].xc = r*cos(phi); circles[n].yc = r*sin(phi); circles[n].radius = LAMBDA*ra[j]; circles[n].active = 1; } break; } case (C_LASER): { ncircles = 17; xx[0] = 0.5*(X_SHOOTER + X_TARGET); xx[1] = LAMBDA - 0.5*(X_TARGET - X_SHOOTER); 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); 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; i 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, circles[i].xc, circles[i].yc); /* 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].yc = 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_GOLDEN_MEAN): { ncircles = 300; gamma = (sqrt(5.0) - 1.0)*0.5; /* golden mean */ height = YMAX - YMIN; dx = 2.0*LAMBDA/((double)ncircles); for (n = 0; n < ncircles; n++) { circles[n].xc = -LAMBDA + n*dx; circles[n].yc = y; y += height*gamma; if (y > YMAX) y -= height; circles[n].radius = MU; circles[n].active = 1; } /* test for circles that overlap top or bottom boundary */ ncirc0 = ncircles; for (n=0; n < ncirc0; n++) { if (circles[n].yc + circles[n].radius > YMAX) { circles[ncircles].xc = circles[n].xc; circles[ncircles].yc = circles[n].yc - height; circles[ncircles].radius = MU; circles[ncircles].active = 1; ncircles ++; } else if (circles[n].yc - circles[n].radius < YMIN) { circles[ncircles].xc = circles[n].xc; circles[ncircles].yc = circles[n].yc + height; circles[ncircles].radius = MU; circles[ncircles].active = 1; ncircles ++; } } 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<1000; i++) { x = r*cos(phi); y = r*sin(phi); phi += gamma; r += MU*r0/r; if ((vabs(x) < LAMBDA)&&(vabs(y) < YMAX + MU)) { circles[ncircles].xc = x; circles[ncircles].yc = y; ncircles++; } } for (i=0; i YMIN - MU)) circles[i].active = 1; // printf("i = %i, circlex = %.3lg, circley = %.3lg\n", i, circles[i].xc, circles[i].yc); } break; } case (C_SQUARE_HEX): { ncircles = NGRIDX*(NGRIDY+1); dy = (YMAX - YMIN)/((double)NGRIDY); dx = dy*0.5*sqrt(3.0); for (i = 0; i < NGRIDX; i++) for (j = 0; j < NGRIDY+1; j++) { n = (NGRIDY+1)*i + j; circles[n].xc = ((double)(i-NGRIDX/2) + 0.5)*dy; /* is +0.5 needed? */ circles[n].yc = YMIN + ((double)j - 0.5)*dy; if (((i+NGRIDX)%4 == 2)||((i+NGRIDX)%4 == 3)) circles[n].yc += 0.5*dy; circles[n].radius = MU; /* activate only circles that intersect the domain */ if ((circles[n].yc < YMAX + MU)&&(circles[n].yc > YMIN - MU)) circles[n].active = 1; else circles[n].active = 0; } break; } case (C_RINGS): { ncircles = NGRIDX*NGRIDY; dphi = DPI/((double)NGRIDX); dr = 0.5*LAMBDA/(double)NGRIDY; for (i = 0; i < NGRIDX; i++) for (j = 0; j < NGRIDY; j++) { n = NGRIDY*i + j; phi = (double)i*dphi; r = 0.5*LAMBDA + (double)j*dr; circles[n].xc = r*cos(phi); circles[n].yc = r*sin(phi); circles[n].radius = MU; /* activate only circles that intersect the domain */ if ((circles[n].yc < YMAX + MU)&&(circles[n].yc > YMIN - MU)) circles[n].active = 1; else circles[n].active = 0; } break; } case (C_RINGS_T): { ncircles = NGRIDX*NGRIDY; dphi = DPI/((double)NGRIDX); dr = 0.5*LAMBDA/(double)NGRIDY; for (i = 0; i < NGRIDX; i++) for (j = 0; j < NGRIDY; j++) { n = NGRIDY*i + j; phi = (double)i*dphi; phi += 0.5*(double)j*dphi; r = 0.5*LAMBDA + (double)j*dr; circles[n].xc = r*cos(phi); circles[n].yc = r*sin(phi); circles[n].radius = MU; /* activate only circles that intersect the domain */ if ((circles[n].yc < YMAX + MU)&&(circles[n].yc > YMIN - MU)) circles[n].active = 1; else circles[n].active = 0; } break; } case (C_RINGS_SPIRAL): { ncircles = 0; // circles[0].xc = 0.5*LAMBDA; // circles[0].yc = 0.0; gamma = (sqrt(5.0) - 1.0)*PI; /* golden mean times 2Pi */ phi = 0.0; r0 = 0.5*LAMBDA; r = r0 + MU; for (i=0; i<1000; i++) { x = r*cos(phi); y = r*sin(phi); phi += gamma; r += 0.1*MU*r0/r; if (x*x + y*y < LAMBDA) { circles[ncircles].xc = x; circles[ncircles].yc = y; ncircles++; } } for (i=0; i YMIN - MU)) circles[i].active = 1; } break; } case (C_HEX_BOTTOM): { ncircles = NGRIDX*(NGRIDY+1); dy = (- YMIN)/((double)NGRIDY); // dx = dy*0.5*sqrt(3.0); dx = (XMAX - XMIN)/((double)NGRIDX); for (i = 0; i < NGRIDX; i++) for (j = 0; j < NGRIDY+1; j++) { n = (NGRIDY+1)*i + j; circles[n].xc = ((double)(i-NGRIDX/2) + 0.5)*dx; /* is +0.5 needed? */ circles[n].yc = YMIN + ((double)j - 0.5)*dy; if ((i+NGRIDX)%2 == 1) circles[n].yc += 0.5*dy; circles[n].radius = MU; /* activate only circles that intersect the domain */ if ((circles[n].yc < YMAX + MU)&&(circles[n].yc > YMIN - MU)) circles[n].active = 1; else circles[n].active = 0; } break; } case (C_HEX_BOTTOM2): { ncircles = NGRIDX*(NGRIDY+1); dy = (- YMIN)/((double)NGRIDY); dx = 2.0*LAMBDA/((double)NGRIDX); for (i = 0; i < NGRIDX; i++) for (j = 0; j < NGRIDY+1; j++) { n = (NGRIDY+1)*i + j; circles[n].xc = ((double)(i-NGRIDX/2) + 0.5)*dx; /* is +0.5 needed? */ circles[n].yc = YMIN + ((double)j - 0.5)*dy; if ((i+NGRIDX)%2 == 1) circles[n].yc += 0.5*dy; circles[n].radius = MU; /* activate only circles that intersect the domain */ if ((circles[n].yc < YMAX + MU)&&(circles[n].yc > YMIN - MU)) circles[n].active = 1; else circles[n].active = 0; } break; } case (C_SQUARE_BOTTOM): { ncircles = NGRIDX*(NGRIDY+1); dy = (- YMIN)/((double)NGRIDY); dx = 2.0*LAMBDA/((double)NGRIDX); for (i = 0; i < NGRIDX; i++) for (j = 0; j < NGRIDY+1; j++) { n = (NGRIDY+1)*i + j; circles[n].xc = ((double)(i-NGRIDX/2) + 0.5)*dx; circles[n].yc = YMIN + ((double)j - 0.5)*dy; circles[n].radius = MU; /* activate only circles that intersect the domain */ if ((circles[n].yc < YMAX + MU)&&(circles[n].yc > YMIN - MU)) circles[n].active = 1; else circles[n].active = 0; } break; } case (C_ONE): { circles[ncircles].xc = 0.0; circles[ncircles].yc = 0.0; circles[ncircles].radius = MU; circles[ncircles].active = 1; ncircles += 1; break; } case (C_TWO): /* used for comparison with cloak */ { circles[ncircles].xc = 0.0; circles[ncircles].yc = 0.0; circles[ncircles].radius = MU; circles[ncircles].active = 2; ncircles += 1; circles[ncircles].xc = 0.0; circles[ncircles].yc = 0.0; circles[ncircles].radius = 2.0*MU; circles[ncircles].active = 1; ncircles += 1; break; } case (C_NOTHING): { ncircles += 0; break; } default: { printf("Function init_circle_config not defined for this pattern \n"); } } return(ncircles); } int init_circle_config(t_circle circles[NMAXCIRCLES]) /* for backward compatibility */ { return (init_circle_config_pattern(circles, CIRCLE_PATTERN)); } int init_polygon_config_pattern(t_polygon polygons[NMAXCIRCLES], int circle_pattern) /* initialise the polygon configuration, for billiard shape D_CIRCLES */ /* uses init_circle_config, this is where C++ would be more elegant */ { int i, ncircles; t_circle circle[NMAXCIRCLES]; ncircles = init_circle_config_pattern(circle, circle_pattern); for (i=0; ix = z.x - 2.0*zperp.x; zprime->y = z.y - 2.0*zperp.y; return(1); } int compute_tokarsky_coordinates(double xshift, double yshift, double scaling, t_vertex polyline[NMAXPOLY]) /* compute positions of vertices of tokarsky room */ { int i; double pos[2]; polyline[0].x = 0.0; polyline[0].y = 2.0; polyline[1].x = 1.0; polyline[1].y = 3.0; polyline[2].x = 1.0; polyline[2].y = 4.0; polyline[3].x = 2.0; polyline[3].y = 4.0; polyline[4].x = 2.0; polyline[4].y = 3.0; polyline[5].x = 3.0; polyline[5].y = 3.0; polyline[6].x = 3.0; polyline[6].y = 2.0; polyline[7].x = 5.0; polyline[7].y = 2.0; polyline[8].x = 5.0; polyline[8].y = 3.0; polyline[9].x = 6.0; polyline[9].y = 3.0; polyline[10].x = 6.0; polyline[10].y = 4.0; polyline[11].x = 7.0; polyline[11].y = 3.0; polyline[12].x = 8.0; polyline[12].y = 3.0; polyline[13].x = 8.0; polyline[13].y = 2.0; polyline[14].x = 7.0; polyline[14].y = 2.0; polyline[15].x = 7.0; polyline[15].y = 1.0; polyline[16].x = 6.0; polyline[16].y = 0.0; polyline[17].x = 6.0; polyline[17].y = 1.0; polyline[18].x = 5.0; polyline[18].y = 1.0; polyline[19].x = 4.0; polyline[19].y = 0.0; polyline[20].x = 4.0; polyline[20].y = 1.0; polyline[21].x = 3.0; polyline[21].y = 1.0; polyline[22].x = 2.0; polyline[22].y = 0.0; polyline[23].x = 2.0; polyline[23].y = 1.0; polyline[24].x = 1.0; polyline[24].y = 1.0; polyline[25].x = 1.0; polyline[25].y = 2.0; for (i=0; i<26; i++) { polyline[i].x = (polyline[i].x + xshift)*scaling; polyline[i].y = (polyline[i].y + yshift)*scaling; xy_to_pos(polyline[i].x, polyline[i].y, pos); polyline[i].posi = pos[0]; polyline[i].posj = pos[1]; } return(26); } void compute_isospectral_coordinates(int type, int ishift, double xshift, double yshift, double scaling, t_vertex polyline[NMAXPOLY]) /* compute positions of vertices of isospectral billiards */ /* central triangle has coordinates (0,0), (1,0) and (LAMBDA,MU) fed into affine transformation */ /* defined by (xshift - 0.5), (yshift - 0.25) and scaling*/ { int i; double pos[2]; polyline[ishift].x = (xshift - 0.5)*scaling; polyline[ishift].y = (yshift - 0.25)*scaling; polyline[ishift+1].x = (0.5+xshift)*scaling; polyline[ishift+1].y = (yshift - 0.25)*scaling; polyline[ishift+2].x = (LAMBDA+xshift - 0.5)*scaling; polyline[ishift+2].y = (MU+yshift - 0.25)*scaling; axial_symmetry_tvertex(polyline[ishift], polyline[ishift+2], polyline[ishift+1], &polyline[ishift+3]); axial_symmetry_tvertex(polyline[ishift], polyline[ishift+1], polyline[ishift+2], &polyline[ishift+4]); axial_symmetry_tvertex(polyline[ishift+1], polyline[ishift+2], polyline[ishift], &polyline[ishift+5]); if (type == 0) { axial_symmetry_tvertex(polyline[ishift], polyline[ishift+3], polyline[ishift+2], &polyline[ishift+6]); axial_symmetry_tvertex(polyline[ishift+1], polyline[ishift+4], polyline[ishift], &polyline[ishift+7]); axial_symmetry_tvertex(polyline[ishift+2], polyline[ishift+5], polyline[ishift+1], &polyline[ishift+8]); } else { axial_symmetry_tvertex(polyline[ishift+2], polyline[ishift+3], polyline[ishift], &polyline[ishift+6]); axial_symmetry_tvertex(polyline[ishift], polyline[ishift+4], polyline[ishift+1], &polyline[ishift+7]); axial_symmetry_tvertex(polyline[ishift+1], polyline[ishift+5], polyline[ishift+2], &polyline[ishift+8]); } for (i=ishift; i NMAXPOLY) { printf("NMAXPOLY needs to be increased to %i\n", nsides); nsides = NMAXPOLY; } for (i=0; i DPI) angle -= DPI; while (angle < 0.0) angle += DPI; xy_to_pos(x*MU, y*MU, pos); polyline[i].posi = pos[0]; polyline[i].posj = pos[1]; } return(nsides); } int compute_star_coordinates(t_vertex polyline[NMAXPOLY]) /* compute positions of vertices of star-shaped domain */ { int i; double alpha, r, x, y, pos[2]; alpha = DPI/(double)NPOLY; for (i=0; i 0.0) x = -MU; else x = MU; polyline[0].x = x; polyline[0].y = -ymax; xy_to_pos(x, -ymax, pos); polyline[0].posi = pos[0]; polyline[0].posj = pos[1]; for (i=1; i 0.0) x = -MU; else x = MU; polyline[NSEG].x = x; polyline[NSEG].y = ymax; xy_to_pos(x, ymax, pos); polyline[NSEG].posi = pos[0]; polyline[NSEG].posj = pos[1]; return(NSEG+1); } int compute_double_fresnel_coordinates(t_vertex polyline[NMAXPOLY], double xshift) /* compute positions of vertices approximating two facing Fresnel lenses */ { int i; double pos[2]; compute_fresnel_coordinates(polyline); for (i=0; i<=NSEG; i++) { polyline[i].x -= xshift; xy_to_pos(polyline[i].x, polyline[i].y, pos); polyline[i].posi = pos[0]; polyline[i].posj = pos[1]; polyline[NSEG + 1 + i].x = -polyline[i].x; polyline[NSEG + 1 + i].y = polyline[i].y; xy_to_pos(polyline[NSEG + 1 + i].x, polyline[NSEG + 1 + i].y, pos); polyline[NSEG + 1 + i].posi = pos[0]; polyline[NSEG + 1 + i].posj = pos[1]; } return(2*NSEG+2); } int compute_noisepanel_coordinates(t_vertex polyline[NMAXPOLY]) /* compute positions of vertices of noise panel */ { int i, n, even; double ymax, dy, x, y, x1, pos[2]; /* find the leftmost point */ x = 0.0; n = 0; while (x > XMIN) { x -= LAMBDA; n++; } if (n%2 == 0) even = 1; else even = 0; i = 0; while (x <= XMAX + LAMBDA) { if (even) y = YMIN + 0.1; else y = YMIN + 0.1 + MU; x1 = x; if (x1 > XMAX) x1 = XMAX; else if (x1 < XMIN) x1 = XMIN; polyline[i].x = x1; polyline[i].y = y; xy_to_pos(x1, y, pos); polyline[i].posi = pos[0]; polyline[i].posj = pos[1]; x += LAMBDA; even = 1 - even; i++; } n = i; for (i=0; i XMIN) { x -= LAMBDA; n++; } if (n%2 == 0) even = 1; else even = 0; i = 0; while (x <= 0.0) { if (even) y = YMIN + 0.1; else y = YMIN + 0.1 + MU; x1 = x; if (x1 > XMAX) x1 = XMAX; else if (x1 < XMIN) x1 = XMIN; polyline[i].x = x1; polyline[i].y = y; xy_to_pos(x1, y, pos); polyline[i].posi = pos[0]; polyline[i].posj = pos[1]; x += LAMBDA; even = 1 - even; i++; } n = i; for (i=0; i0)||(j!=ropening))&&(maze[n].west)) { polyrect[nsides].x1 = x1 - width; polyrect[nsides].y1 = y1 - width; polyrect[nsides].x2 = x1 + width; polyrect[nsides].y2 = y1 + width + dy; nsides++; } if (maze[n].south) { polyrect[nsides].x1 = x1 - width; polyrect[nsides].y1 = y1 - width; polyrect[nsides].x2 = x1 + width + dx; polyrect[nsides].y2 = y1 + width; nsides++; } } /* top side of maze */ polyrect[nsides].x1 = YMIN + padding + MAZE_XSHIFT; polyrect[nsides].y1 = YMAX - padding - width; polyrect[nsides].x2 = YMAX - padding + MAZE_XSHIFT; polyrect[nsides].y2 = YMAX - padding + width; nsides++; /* right side of maze */ y1 = YMIN + padding + dy*((double)ropening); x1 = YMAX - padding + MAZE_XSHIFT; polyrect[nsides].x1 = x1 - width; polyrect[nsides].y1 = YMIN - 1.0; polyrect[nsides].x2 = x1 + width; polyrect[nsides].y2 = y1 - dy; nsides++; polyrect[nsides].x1 = x1 - width; polyrect[nsides].y1 = y1; polyrect[nsides].x2 = x1 + width; polyrect[nsides].y2 = YMAX + 1.0; nsides++; /* left side of maze */ x1 = YMIN + padding + MAZE_XSHIFT; polyrect[nsides].x1 = x1 - width; polyrect[nsides].y1 = YMIN - 1.0; polyrect[nsides].x2 = x1 + width; polyrect[nsides].y2 = YMIN + padding; nsides++; polyrect[nsides].x1 = x1 - width; polyrect[nsides].y1 = YMAX - padding; polyrect[nsides].x2 = x1 + width; polyrect[nsides].y2 = YMAX + 1.0; nsides++; if (type == 1) /* maze with closed sides */ { polyrect[nsides].x1 = XMIN - 0.5*width; polyrect[nsides].y1 = YMIN - 0.5*width; polyrect[nsides].x2 = XMIN + 0.5*width; polyrect[nsides].y2 = YMAX + 0.5*width; nsides++; polyrect[nsides].x1 = XMIN - 0.5*width; polyrect[nsides].y1 = YMIN - 0.5*width; polyrect[nsides].x2 = x1 + 0.5*width; polyrect[nsides].y2 = YMIN + 0.5*width; nsides++; polyrect[nsides].x1 = XMIN - 0.5*width; polyrect[nsides].y1 = YMAX - 0.5*width; polyrect[nsides].x2 = x1 + 0.5*width; polyrect[nsides].y2 = YMAX + 0.5*width; nsides++; } else if (type == 2) /* maze with channels */ { /* right channel */ y1 = YMIN + padding + dy*((double)ropening); x1 = YMAX - padding + MAZE_XSHIFT; polyrect[nsides].x1 = x1 - 0.5*width; polyrect[nsides].y1 = YMIN - padding; polyrect[nsides].x2 = XMAX + padding; polyrect[nsides].y2 = y1 - dy + width; nsides++; polyrect[nsides].x1 = x1 - 0.5*width; polyrect[nsides].y1 = y1 - width; polyrect[nsides].x2 = XMAX + padding; polyrect[nsides].y2 = YMAX + padding; nsides++; /* left channel */ x1 = YMIN + padding + MAZE_XSHIFT; polyrect[nsides].x1 = XMIN - padding; polyrect[nsides].y1 = YMIN - padding; polyrect[nsides].x2 = x1 + 0.5*width; polyrect[nsides].y2 = y1 - dy + width; nsides++; polyrect[nsides].x1 = XMIN - padding; polyrect[nsides].y1 = y1 - width; polyrect[nsides].x2 = x1 + 0.5*width; polyrect[nsides].y2 = YMAX + padding; nsides++; } for (i=0; i 0)&&(maze[n].west)) { polyarc[na].xc = MAZE_XSHIFT; polyarc[na].yc = 0.0; polyarc[na].r = r; polyarc[na].angle1 = phi; polyarc[na].dangle = dphi; polyarc[na].width = width; na++; } /* second circle */ r = rmin + dr; dphi *= 0.5; for (q=0; q<2; q++) { n = nmaze(1, block*NXMAZE + q); phi = (double)(block)*angle + (double)q*dphi; if (maze[n].west) { polyarc[na].xc = MAZE_XSHIFT; polyarc[na].yc = 0.0; polyarc[na].r = r; polyarc[na].angle1 = phi; polyarc[na].dangle = dphi; polyarc[na].width = width; na++; } } /* other circles */ ww = 2; i = 2; while (ww < NXMAZE) { dphi *= 0.5; for (p = 0; p < ww; p++) { r = rmin + (double)i*dr; printf("Circle, i = %i, dphi = %.2lg, r = %.2lg\n", i, dphi, r); for (q = 0; q < 2*ww; q++) { j = block*NXMAZE + q; n = nmaze(i,j); phi = (double)(block)*angle + (double)q*dphi; if (maze[n].west) { polyarc[na].xc = MAZE_XSHIFT; polyarc[na].yc = 0.0; polyarc[na].r = r; polyarc[na].angle1 = phi; polyarc[na].dangle = dphi; polyarc[na].width = width; na++; } } i++; } ww *= 2; } } /* outer boundary of maze */ polyarc[na].xc = MAZE_XSHIFT; polyarc[na].yc = 0.0; polyarc[na].r = rmax; polyarc[na].angle1 = dphi; polyarc[na].dangle = DPI - dphi; polyarc[na].width = width; na++; *npolyrect_rot = np; *npolyarc = na; free(maze); } int compute_interior_maze_coordinates(t_rectangle polyrect[NMAXPOLY], int type) /* compute positions of complement of maze */ { t_maze* maze; int i, j, n, nsides = 0, ropening; double dx, dy, x1, y1, x0, padding = 0.02, pos[2], width = MAZE_WIDTH; /* TODO */ maze = (t_maze *)malloc(NXMAZE*NYMAZE*sizeof(t_maze)); ropening = (NYMAZE+1)/2; init_maze(maze); /* move the entrance for maze type with two channels */ if (type == 2) { n = nmaze(0, ropening-1); maze[n].west = 0; n = nmaze(0, ropening); maze[n].west = 1; } /* build walls of maze */ // x0 = LAMBDA - 1.0; dx = (YMAX - YMIN - 2.0*padding)/(double)(NXMAZE); dy = (YMAX - YMIN - 2.0*padding)/(double)(NYMAZE); for (i=0; i= 0.0)&&(v2 >= 0.0)&&(v3 >= 0.0)) return(1); else return(0); } int xy_in_triangle_tvertex(double x, double y, t_vertex z1, t_vertex z2, t_vertex z3) /* returns 1 iff (x,y) is inside the triangle with vertices z1, z2, z3 */ { double v1, v2, v3; /* compute wedge products */ v1 = (z2.x - z1.x)*(y - z1.y) - (z2.y - z1.y)*(x - z1.x); v2 = (z3.x - z2.x)*(y - z2.y) - (z3.y - z2.y)*(x - z2.x); v3 = (z1.x - z3.x)*(y - z3.y) - (z1.y - z3.y)*(x - z3.x); if ((v1 >= 0.0)&&(v2 >= 0.0)&&(v3 >= 0.0)) return(1); else return(0); } int xy_in_polyrect(double x, double y, t_rectangle rectangle) /* returns 1 if (x,y) is in rectangle */ { double x1, y1, x2, y2; if (rectangle.x1 < rectangle.x2) { x1 = rectangle.x1; x2 = rectangle.x2; } else { x1 = rectangle.x2; x2 = rectangle.x1; } if (rectangle.y1 < rectangle.y2) { y1 = rectangle.y1; y2 = rectangle.y2; } else { y1 = rectangle.y2; y2 = rectangle.y1; } if (x < x1) return(0); if (x > x2) return(0); if (y < y1) return(0); if (y > y2) return(0); return(1); } int ij_in_polyrect(double i, double j, t_rectangle rectangle) /* returns 1 if (x,y) is in rectangle */ { int i1, i2, j1, j2; if (rectangle.posi1 < rectangle.posi2) { i1 = rectangle.posi1; i2 = rectangle.posi2; } else { i1 = rectangle.posi2; i2 = rectangle.posi1; } if (rectangle.posj1 < rectangle.posj2) { j1 = rectangle.posj1; j2 = rectangle.posj2; } else { j1 = rectangle.posj2; j2 = rectangle.posj1; } if (i < i1) return(0); if (i > i2) return(0); if (j < j1) return(0); if (j > j2) return(0); return(1); } int xy_in_rectrotated(double x, double y, t_rect_rotated rectrot) /* returns 1 if (x,y) is in rectangle */ { double l, u1, u2, v1, v2, pscal, h2; l = module2(rectrot.x2 - rectrot.x1, rectrot.y2 - rectrot.y1); if (l == 0.0) return(0); /* unit vector along axis */ u1 = (rectrot.x2 - rectrot.x1)/l; u2 = (rectrot.y2 - rectrot.y1)/l; /* vector from one extremity to (x,y) */ v1 = x - rectrot.x1; v2 = y - rectrot.y1; /* inner product */ pscal = u1*v1 + u2*v2; if (pscal < 0.0) return(0); if (pscal > l) return(0); h2 = v1*v1 + v2*v2 - pscal*pscal; return(4.0*h2 <= rectrot.width*rectrot.width); } int xy_in_arc(double x, double y, t_arc arc) /* returns 1 if (x,y) is in arc */ { double rho, phi, alpha; rho = module2(x - arc.xc, y - arc.yc); if (vabs(rho - arc.r) > 0.5*arc.width) return(0); phi = argument(x - arc.xc, y - arc.yc); alpha = phi - arc.angle1; while (alpha < 0.0) alpha += DPI; while (alpha > DPI) alpha -= DPI; return(alpha <= arc.dangle); } int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_circle *circles) /* returns 1 if (x,y) represents a point in the billiard */ { double l2, r1, r2, r2mu, omega, b, c, d, angle, z, x1, y1, x2, y2, y3, u, v, u1, v1, dx, dy, width, alpha, s, a, r, height, ca, sa, l, ht, xshift, zz[9][2], x5, x6, f, fp, h1; int i, j, k, k1, k2, condition = 0, m; static int first = 1, nsides; static double h, hh, ra, rb, ll, salpha; switch (b_domain) { case (D_NOTHING): { return(1); break; } case (D_RECTANGLE): { if ((vabs(x) 1.0) return(1); else return(0); break; } case (D_EXT_ELLIPSE_CURVED): { y1 = y + 0.4*x*x; if (x*x/(LAMBDA*LAMBDA) + y1*y1/(MU*MU) > 1.0) return(1); else return(0); break; } case (D_EXT_ELLIPSE_CURVED_BDRY): { if (y > YMAX - 0.05) return(0); if (y < YMIN + 0.05) return(0); y1 = y + 0.4*x*x; if (x*x/(LAMBDA*LAMBDA) + y1*y1/(MU*MU) > 1.0) return(1); else return(0); break; } case (D_STADIUM): { if ((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_FLAT): { if (y > -LAMBDA) return(1); else return(0); break; } case (D_ANNULUS): { l2 = LAMBDA*LAMBDA; r2 = x*x + y*y; if ((r2 > l2)&&(r2 < 1.0)) return(1); else return(0); } case (D_POLYGON): { condition = 1; omega = DPI/((double)NPOLY); c = cos(omega*0.5); for (k=0; k MU)) return(1); else if ((vabs(y-LAMBDA) < MU)||(vabs(y+LAMBDA) < MU)) return (1); else return(0); } case (D_GRATING): { k1 = -(int)((-YMIN)/LAMBDA); k2 = (int)(YMAX/LAMBDA); condition = 1; for (i=k1; i<= k2; i++) { z = (double)i*LAMBDA; condition = condition*(x*x + (y-z)*(y-z) > MU*MU); } // printf("x = %.3lg, y = %.3lg, k1 = %i, k2 = %i, condition = %i\n", x, y, k1, k2, condition); return(condition); } case (D_EHRENFEST): { return(((x-1.0)*(x-1.0) + y*y < LAMBDA*LAMBDA)||((x+1.0)*(x+1.0) + y*y < LAMBDA*LAMBDA)||((vabs(x) < 1.0)&&(vabs(y) < MU))); } case (D_DISK_GRID): { dy = (YMAX - YMIN)/((double)NGRIDY); for (i = -NGRIDX/2; i < NGRIDX/2; i++) for (j = 0; j < NGRIDY; j++) { x1 = ((double)i + 0.5)*dy; y1 = YMIN + ((double)j + 0.5)*dy; if ((x-x1)*(x-x1) + (y-y1)*(y-y1) < MU*MU) return(0); } return(1); } case (D_DISK_HEX): { dy = (YMAX - YMIN)/((double)NGRIDY); dx = dy*0.5*sqrt(3.0); for (i = -NGRIDX/2; i < NGRIDX/2; i++) for (j = -1; j < NGRIDY; j++) { x1 = ((double)i + 0.5)*dy; y1 = YMIN + ((double)j + 0.5)*dy; if ((i+NGRIDX)%2 == 1) y1 += 0.5*dy; if ((x-x1)*(x-x1) + (y-y1)*(y-y1) < MU*MU) return(0); } return(1); } case (D_PARABOLA): { return(x > 0.25*y*y/LAMBDA - LAMBDA); } case (D_TWO_PARABOLAS): { x1 = 0.25*y*y/MU - MU - LAMBDA; x2 = -x1; width = 0.25*MU; if (width > 0.2) width = 0.2; if (vabs(y) > 1.5*MU) return(1); else if ((x < x1 - width)||(x > x2 + width)) return(1); else if ((x > x1)&&(x < x2)) return(1); else return(0); } case (D_FOUR_PARABOLAS): { x1 = MU + LAMBDA - 0.25*y*y/MU; y1 = MU + LAMBDA - 0.25*x*x/MU; return((vabs(x) < x1)&&(vabs(y) < y1)); } case (D_POLY_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) + 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_HYPERBOLA): { b = MU*sqrt(1.0 + x*x/(LAMBDA*LAMBDA - MU*MU)); if (y > 1.02*b) return(1); else if (y < 0.98*b) return (1); else return(0); } case (D_TOKARSKY): { x1 = 4.0 + x/(XMAX - XMIN)*8.4; y1 = 2.0 + y/(XMAX - XMIN)*8.4; if ((x1 <= 0.0)||(x1 >= 8.0)) return(0); else if (x1 < 1.0) { if (y1 <= 2.0) return(0); else if (y1 >= x1 + 2.0) return(0); else return(1); } else if (x1 < 2.0) { if (y1 <= 1.0) return(0); else if (y1 >= 4.0) return(0); else return(1); } else if (x1 < 3.0) { if (y1 <= x1 - 2.0) return(0); else if (y1 >= 3.0) return(0); else return(1); } else if (x1 < 4.0) { if (y1 <= 1.0) return(0); else if (y1 >= 2.0) return(0); else return(1); } else if (x1 < 5.0) { if (y1 <= x1 - 4.0) return(0); else if (y1 >= 2.0) return(0); else return(1); } else if (x1 < 6.0) { if (y1 <= 1.0) return(0); else if (y1 >= 3.0) return(0); else return(1); } else if (x1 < 7.0) { if (y1 <= x1 - 6.0) return(0); else if (y1 >= 10.0 - x1) return(0); else return(1); } else { if (y1 <= 2.0) return(0); else if (y1 >= 3.0) return(0); else return(1); } } case (D_TOKA_PRIME): { // x1 = vabs(x); if (x + MU > 0.0) x1 = x; else x1 = -2.0*MU - x; condition = xy_in_triangle_tvertex(x1, y, polyline[0], polyline[1], polyline[2]); condition += xy_in_triangle_tvertex(x1, y, polyline[0], polyline[2], polyline[3]); condition += xy_in_triangle_tvertex(x1, y, polyline[i], polyline[3], polyline[4]); for (i=3; i<42; i++) condition += xy_in_triangle_tvertex(x1, y, polyline[i], polyline[43], polyline[i+1]); condition += xy_in_triangle_tvertex(x1, y, polyline[42], polyline[43], polyline[3]); return(condition >= 1); } case (D_ISOSPECTRAL): { /* 1st triangle */ condition = xy_in_triangle_tvertex(x, y, polyline[0], polyline[1], polyline[2]); condition += xy_in_triangle_tvertex(x, y, polyline[0], polyline[4], polyline[1]); condition += xy_in_triangle_tvertex(x, y, polyline[1], polyline[5], polyline[2]); condition += xy_in_triangle_tvertex(x, y, polyline[0], polyline[2], polyline[3]); condition += xy_in_triangle_tvertex(x, y, polyline[1], polyline[4], polyline[7]); condition += xy_in_triangle_tvertex(x, y, polyline[2], polyline[5], polyline[8]); condition += xy_in_triangle_tvertex(x, y, polyline[0], polyline[3], polyline[6]); /* 2nd triangle */ condition += xy_in_triangle_tvertex(x, y, polyline[9], polyline[10], polyline[11]); condition += xy_in_triangle_tvertex(x, y, polyline[9], polyline[13], polyline[10]); condition += xy_in_triangle_tvertex(x, y, polyline[10], polyline[14], polyline[11]); condition += xy_in_triangle_tvertex(x, y, polyline[9], polyline[11], polyline[12]); condition += xy_in_triangle_tvertex(x, y, polyline[9], polyline[16], polyline[13]); condition += xy_in_triangle_tvertex(x, y, polyline[10], polyline[17], polyline[14]); condition += xy_in_triangle_tvertex(x, y, polyline[11], polyline[15], polyline[12]); return(condition >= 1); } case (D_HOMOPHONIC): { /* conditions could be summarised in larger triangles, but this is to keep */ /* the option of using triangles with other angles than 30-60-90 */ /* 1st triangle */ condition = xy_in_triangle_tvertex(x, y, polyline[2], polyline[0], polyline[1]); condition += xy_in_triangle_tvertex(x, y, polyline[2], polyline[1], polyline[3]); condition += xy_in_triangle_tvertex(x, y, polyline[0], polyline[21], polyline[1]); condition += xy_in_triangle_tvertex(x, y, polyline[0], polyline[10], polyline[21]); condition += xy_in_triangle_tvertex(x, y, polyline[10], polyline[11], polyline[21]); condition += xy_in_triangle_tvertex(x, y, polyline[11], polyline[13], polyline[21]); condition += xy_in_triangle_tvertex(x, y, polyline[11], polyline[12], polyline[13]); condition += xy_in_triangle_tvertex(x, y, polyline[13], polyline[14], polyline[21]); condition += xy_in_triangle_tvertex(x, y, polyline[14], polyline[20], polyline[21]); condition += xy_in_triangle_tvertex(x, y, polyline[14], polyline[15], polyline[20]); condition += xy_in_triangle_tvertex(x, y, polyline[15], polyline[19], polyline[20]); condition += xy_in_triangle_tvertex(x, y, polyline[2], polyline[4], polyline[5]); condition += xy_in_triangle_tvertex(x, y, polyline[2], polyline[5], polyline[7]); condition += xy_in_triangle_tvertex(x, y, polyline[5], polyline[6], polyline[7]); condition += xy_in_triangle_tvertex(x, y, polyline[2], polyline[7], polyline[8]); condition += xy_in_triangle_tvertex(x, y, polyline[2], polyline[8], polyline[0]); condition += xy_in_triangle_tvertex(x, y, polyline[0], polyline[8], polyline[9]); condition += xy_in_triangle_tvertex(x, y, polyline[0], polyline[9], polyline[10]); condition += xy_in_triangle_tvertex(x, y, polyline[15], polyline[16], polyline[19]); condition += xy_in_triangle_tvertex(x, y, polyline[16], polyline[17], polyline[18]); condition += xy_in_triangle_tvertex(x, y, polyline[16], polyline[18], polyline[19]); /* 2nd triangle */ condition += xy_in_triangle_tvertex(x, y, polyline[22+2], polyline[22+0], polyline[22+1]); condition += xy_in_triangle_tvertex(x, y, polyline[22+2], polyline[22+1], polyline[22+3]); condition += xy_in_triangle_tvertex(x, y, polyline[22+0], polyline[22+21], polyline[22+1]); condition += xy_in_triangle_tvertex(x, y, polyline[22+0], polyline[22+10], polyline[22+21]); condition += xy_in_triangle_tvertex(x, y, polyline[22+10], polyline[22+11], polyline[22+21]); condition += xy_in_triangle_tvertex(x, y, polyline[22+11], polyline[22+13], polyline[22+21]); condition += xy_in_triangle_tvertex(x, y, polyline[22+11], polyline[22+12], polyline[22+13]); condition += xy_in_triangle_tvertex(x, y, polyline[22+13], polyline[22+14], polyline[22+21]); condition += xy_in_triangle_tvertex(x, y, polyline[22+14], polyline[22+20], polyline[22+21]); condition += xy_in_triangle_tvertex(x, y, polyline[22+14], polyline[22+15], polyline[22+20]); condition += xy_in_triangle_tvertex(x, y, polyline[22+15], polyline[22+19], polyline[22+20]); condition += xy_in_triangle_tvertex(x, y, polyline[22+2], polyline[22+3], polyline[22+5]); condition += xy_in_triangle_tvertex(x, y, polyline[22+3], polyline[22+4], polyline[22+5]); condition += xy_in_triangle_tvertex(x, y, polyline[22+2], polyline[22+5], polyline[22+6]); condition += xy_in_triangle_tvertex(x, y, polyline[22+2], polyline[22+6], polyline[22+8]); condition += xy_in_triangle_tvertex(x, y, polyline[22+2], polyline[22+8], polyline[22+9]); condition += xy_in_triangle_tvertex(x, y, polyline[22+6], polyline[22+7], polyline[22+8]); condition += xy_in_triangle_tvertex(x, y, polyline[22+11], polyline[22+10], polyline[22+16]); condition += xy_in_triangle_tvertex(x, y, polyline[22+11], polyline[22+16], polyline[22+18]); condition += xy_in_triangle_tvertex(x, y, polyline[22+11], polyline[22+18], polyline[22+12]); condition += xy_in_triangle_tvertex(x, y, polyline[22+16], polyline[22+17], polyline[22+18]); return(condition >= 1); } case (D_CIRCLES): { for (i = 0; i < ncirc; i++) if (circles[i].active) { x1 = circles[i].xc; y1 = circles[i].yc; r2 = circles[i].radius*circles[i].radius; if ((x-x1)*(x-x1) + (y-y1)*(y-y1) < r2) return(0); } return(1); } case (D_CIRCLES_IN_RECT): /* returns 2 inside circles, 0 outside rectangle */ { for (i = 0; i < ncirc; i++) if (circles[i].active) { x1 = circles[i].xc; y1 = circles[i].yc; r2 = circles[i].radius*circles[i].radius; if ((x-x1)*(x-x1) + (y-y1)*(y-y1) < r2) return(2); } if ((vabs(x) >= LAMBDA)||(vabs(y) >= 1.0)) return(0); else return(1); } case (D_POLYGONS): { for (i = 0; i < ncirc; i++) if ((polygons[i].active)&&(in_tpolygon(x, y, polygons[i]))) return(0); return(1); } case (D_VONKOCH): { condition = xy_in_triangle_tvertex(x, y, polyline[0], polyline[npolyline/3], polyline[2*npolyline/3]); m = 1; k = 1; for (i = 0; i < MDEPTH; i++) { m = m*4; for (j = 0; j < npolyline/m; j++) condition += xy_in_triangle_tvertex(x, y, polyline[j*m + k], polyline[j*m + 2*k], polyline[j*m + 3*k]); k = k*4; } return(condition >= 1); } case (D_STAR): { condition = xy_in_triangle_tvertex(x, y, polyline[NPOLY], polyline[NPOLY-1], polyline[0]); for (i = 0; i < NPOLY-1; i++) condition += xy_in_triangle_tvertex(x, y, polyline[NPOLY], polyline[i], polyline[i+1]); return(condition >= 1); } case (D_FRESNEL): { if (vabs(y) > 0.9*vabs(LAMBDA)) return(1); if (vabs(x) > MU) return(1); x1 = sqrt(LAMBDA*LAMBDA - y*y) - vabs(LAMBDA); while (x1 <= 0.0) x1 += MU; if (LAMBDA > 0.0) { if (x < x1) return(0); else return(1); } else { x1 = -x1; if (x > x1) return(0); else return(1); } } case (D_DOUBLE_FRESNEL): { if (vabs(y) > 0.9*vabs(LAMBDA)) return(1); if (LAMBDA > 0.0) { if (vabs(x) > LAMBDA + MU) return(1); x1 = sqrt(LAMBDA*LAMBDA - y*y) - LAMBDA; while (x1 <= 0.0) x1 += MU; x1 -= LAMBDA; if (vabs(x) > -x1) return(0); else return(1); } else { if (vabs(x) < -LAMBDA - MU) return(1); x1 = sqrt(LAMBDA*LAMBDA - y*y) + LAMBDA; while (x1 <= 0.0) x1 += MU; x1 -= LAMBDA; if (vabs(x) > x1) return(1); else return(0); } } case (D_NOISEPANEL): { x1 = vabs(x); while (x1 > 2.0*LAMBDA) x1 -= 2.0*LAMBDA; if (x1 <= LAMBDA) y1 = 0.1 + MU*x1/LAMBDA; else y1 = 0.1 + 2.0*MU - MU*x1/LAMBDA; return((y > YMIN + y1)&&(y < YMAX - y1)); } case (D_NOISEPANEL_RECT): { x1 = -x; if (x1 > NPWIDTH) { while (x1 > 2.0*LAMBDA) x1 -= 2.0*LAMBDA; if (x1 <= LAMBDA) y1 = 0.1 + MU*x1/LAMBDA; else y1 = 0.1 + 2.0*MU - MU*x1/LAMBDA; return((y > YMIN + y1)&&(y < YMAX - y1)&&(x > XMIN + 0.1)); } else if (x > NPWIDTH) { return((vabs(y) < YMAX - 0.1)&&(x < XMAX - 0.1)); } else return(0); } case (D_QRD): { x1 = vabs(x)/LAMBDA; k = (int)(x1 + 0.5); k1 = (k*k) % 13; y1 = (MU/13.0)*(14.0 - (double)k1); return ((y > YMIN + y1)&&(y < YMAX - y1)); } case (D_QRD_ASYM): { if (y > 0.0) { x1 = vabs(x)/LAMBDA; k = (int)(x1 + 0.5); k1 = (k*k) % 13; y1 = (MU/13.0)*(14.0 - (double)k1); return (y < YMAX - y1); } else { x1 = vabs(x + 1.0)/LAMBDA; k = (int)(x1 + 0.5); k1 = (k*k) % 17; y1 = (MU/17.0)*(18.0 - (double)k1); return (y > YMIN + y1); } } case (D_CIRCLE_SEGMENT): { if (vabs(y) > 0.9*vabs(LAMBDA)) return(1); y1 = 0.9*LAMBDA; x1 = sqrt(LAMBDA*LAMBDA - y1*y1) - vabs(LAMBDA) + MU; if ((LAMBDA > 0.0)&&(x < x1)) return(1); else if ((LAMBDA < 0.0)&&(x > -x1)) return(1); x1 = sqrt(LAMBDA*LAMBDA - y*y) - vabs(LAMBDA) + MU; if (LAMBDA > 0.0) { if (x < x1) return(0); else return(1); } else { if (x > -x1) return(0); else return(1); } } case (D_GROOVE): { s = 0.85*LAMBDA; a = 0.5*LAMBDA; x1 = x - XMIN - (double)((int)((x - XMIN)/LAMBDA))*LAMBDA; if (x1 < a) return (y > YMIN + LAMBDA); else return (y > YMIN + LAMBDA + s); } case (D_FABRY_PEROT): { return(vabs(x - y*LAMBDA/YMAX) > 0.5*MU); } case (D_LSHAPE): { if (vabs(x) > LAMBDA) return(0); else if (vabs(y) > 1.0) return(0); else if ((x > 0.0)&&(y > 0.0)) return(0); else return(1); } case (D_WAVEGUIDE): { x1 = XMIN + MU; x2 = XMAX - 2.0*MU - 1.5*LAMBDA; y1 = 0.5*LAMBDA; y2 = 1.5*LAMBDA; if (x < x1) return(0); if (x > x2 + 1.5*LAMBDA) return(0); if (vabs(y) > y2) return(0); if (x < x2) return(vabs(y) >= y1); r = module2(x-x2, y); if (r < 0.5*LAMBDA) return(0); if (r > 1.5*LAMBDA) return(0); return(1); } case (D_WAVEGUIDE_W): { x1 = vabs(x); width = LAMBDA - 2.0*MU; height = 0.5*MU; if (x1 > 2.0*LAMBDA - MU) return(0); if (vabs(y) > MU + width + height) return(0); if (y >= height) { r = module2(x1, y-height); if ((r > MU)&&(r < MU + width)) return(1); if (x1 > LAMBDA + MU) return(1); return(0); } if (y <= -height) { r = module2(x1-LAMBDA, y+height); if ((r > MU)&&(r < MU + width)) return(1); return(0); } if (x1 > LAMBDA + MU) return(1); if ((x1 > MU)&&(x1 < MU + width)) return(1); return(0); } case (D_WAVEGUIDES_W): { x1 = vabs(x); /* upper fiber */ width = LAMBDA - 2.0*MU; height = MU + width; if (y > MU + width + height) return(0); if (y >= height) { r = module2(x1, y-height); if ((r > MU)&&(r < MU + width)) return(1); if ((y < height + MU)&&(x1 > LAMBDA + MU)&&(x1 < LAMBDA + MU + width)) return(1); return(0); } else { r = module2(x1-LAMBDA, y-height); if ((r > MU)&&(r < MU + width)) return(1); } /* lower fiber */ height = -MU - width; width = LAMBDA - 2.0*MU_B; if (y < -MU_B - width + height) return(0); if (y >= +height) { r = module2(x1, y-height); if ((r > MU_B)&&(r < MU_B + width)) return(1); if ((y < height + MU_B)&&(x1 > LAMBDA + MU_B)&&(x1 < LAMBDA + MU_B + width)) return(1); return(0); } else { r = module2(x1-LAMBDA, y-height); if ((r > MU_B)&&(r < MU_B + width)) return(1); } return(0); } case (D_WAVEGUIDES_COUPLED): { x1 = vabs(x)*2.0/XMAX; y1 = vabs(y); width = 0.03; b = 1.0*(1.0 - 2.0*MU - width); a = MU + 0.5*width - b; x6 = x1*x1*x1; x5 = x6*x1*x1; x6 = x5*x1; f = a + b*(1.0 + 2.0*x6)/(1.0 + x6); fp = 6.0*b*x5/((1.0 + x6)*(1.0 + x6)); return(vabs(f - y1) < MU*sqrt(1.0 + fp*fp)); } case (D_WAVEGUIDES_COUPLED_N): { x1 = vabs(x); y1 = vabs(y); width = 0.01; b = 1.0*(1.0 - 2.0*MU - width); a = MU + 0.5*width - b; c = 15.0; x6 = x1*x1*x1; x5 = x6*x1*x1; x6 = x5*x1; f = a + b*(1.0 + 2.0*c*x6)/(1.0 + c*x6); fp = 6.0*b*c*x5/((1.0 + c*x6)*(1.0 + c*x6)); if (vabs(f - y1) < 0.9*MU*sqrt(1.0 + fp*fp)) return(1); if (vabs(f - y1) < MU*sqrt(1.0 + fp*fp)) return(2); return(0); } case (D_WAVEGUIDE_S): { a = 0.5*(1.0 - MU); if (x < 0.0) { x = -x; y = -y; } if (x < LAMBDA) { if (vabs(y - 2.0*a) < MU) return(1); if (vabs(y) < MU) return(1); if (vabs(y + 2.0*a) < MU) return(1); return(0); } r = module2(x-LAMBDA, y-a); return(vabs(r - a) < MU); } case (D_WAVEGUIDE_S_SHORT): { a = -1.0 + MU; if (y < 0.0) { x = -x; y = -y; } if (x < -LAMBDA) { if (vabs(y - 1.0 + MU) > MU) return(0); if (vabs(y - 1.0 + MU) > 0.8*MU) return(2); return(1); } r = module2(x-a, y); if (vabs(r - 1.0 + MU) < 0.8*MU) return(1); if (vabs(r - 1.0 + MU) < MU) return(2); return(0); } case (D_WAVEGUIDE_BADSPLICE): { if (x < 0.0) { y1 = vabs(vabs(y) - 0.5); } else { if (y > 0.0) y1 = vabs(y - 0.5 - MU); else y1 = vabs(y + 0.5 - MU_B); } if (y1 < 0.9*LAMBDA) return(1); if (y1 < LAMBDA) return(2); return(0); } case (D_MAZE): { for (i=0; i polyrect[i].x1)&&(x < polyrect[i].x2)&&(y > polyrect[i].y1)&&(y < polyrect[i].y2)) return(0); return(1); } case (D_MAZE_CLOSED): { for (i=0; i polyrect[i].x1)&&(x < polyrect[i].x2)&&(y > polyrect[i].y1)&&(y < polyrect[i].y2)) return(0); return(1); } case (D_MAZE_CHANNELS): { for (i=0; i polyrect[i].x1)&&(x < polyrect[i].x2)&&(y > polyrect[i].y1)&&(y < polyrect[i].y2)) return(0); return(1); } case (D_MAZE_CIRCULAR): { for (i=0; i 2.0)&&(y1 < 1.0)) return(1); if ((x1 < 1.0)&&(y1 > 2.0)) return(1); if (x1 + y1 < 1.0) return(1); if (x1 + y1 > 5.0) return(1); return(0); } case (D_FUNNELS): { y1 = y; if (y > 0.5*YMAX) y1 -= YMAX; if (y < -0.5*YMAX) y1 += YMAX; y1 = vabs(y1 - MU*x); y1 = vabs(y1 - 0.5*YMAX)*2.0/YMAX; if (y1 > 0.25*(1.0 + LAMBDA + x*x)) return(0); return(1); } case (D_ONE_FUNNEL): { y1 = vabs(y); if (y1 > MU + LAMBDA*x*x) return(0); return(1); } case (D_LENSES_RING): { if (first) { salpha = DPI/(double)NPOLY; h = LAMBDA*tan(PI/(double)NPOLY); if (h < MU) ll = sqrt(MU*MU - h*h); else ll = 0.0; first = 0; } for (i=0; i MU) return(1); a = LAMBDA - 0.25; if (x > 0) return ((x+a)*(x+a) + y*y > LAMBDA*LAMBDA); return ((x-a)*(x-a) + y*y > LAMBDA*LAMBDA); } case (D_LENS_WALL): { if (vabs(y) > MU) return(1); a = LAMBDA - 0.25; if (x > 0) return ((x+a)*(x+a) + y*y > LAMBDA*LAMBDA); return ((x-a)*(x-a) + y*y > LAMBDA*LAMBDA); } case (D_TWO_LENSES_WALL): { width = 0.2; a = 1.9; b = 1.9 - 2.0*LAMBDA + 2.0*width; x1 = vabs(x); if ((x1 < WALL_WIDTH)&&(vabs(y) > MU)) return(2); if (module2(x1 - a, y) > LAMBDA) return(1); if (module2(x1 - b, y) > LAMBDA) return(1); return(0); } case (D_TWO_LENSES_OBSTACLE): { width = 0.2; a = 1.9; b = 1.9 - 2.0*LAMBDA + 2.0*width; x1 = vabs(x); if (module2(x1 - a, y) > LAMBDA) return(1); if (module2(x1 - b, y) > LAMBDA) return(1); return(0); } case (D_FRESNEL_ZONE_PLATE): { a = 6.25; b = 0.125; if (vabs(x + LAMBDA) >= MU) return(1); return(cos(DPI*(a*y*y + b)) < 0.0); } case (D_FRESNEL_ZONE_PLATE_INV): { a = 6.25; if (vabs(x + LAMBDA) >= MU) return(1); return(cos(DPI*a*y*y) > 0.0); } case (D_LENS_ROTATED): { ca = cos(APOLY*PID); sa = sin(APOLY*PID); x1 = x*ca + y*sa; y1 = -x*sa + y*ca; if (vabs(y1) > MU) return(1); a = LAMBDA - 0.25; if (x1 > 0) return ((x1+a)*(x1+a) + y1*y1 > LAMBDA*LAMBDA); return ((x1-a)*(x1-a) + y1*y1 > LAMBDA*LAMBDA); } case (D_LENS_CONCAVE): { width = 0.1; x1 = vabs(x); if (vabs(y) > MU) return(1); return((x1 - LAMBDA - width)*(x1 - LAMBDA - width) + y*y < LAMBDA*LAMBDA); } case (D_LENS_CONVEX_CONCAVE): { if (vabs(y) > MU) return(1); xshift = 0.6; if (x > 0.0) { x1 = x - xshift; a = LAMBDA - 0.25; if (x1 > 0) return ((x1+a)*(x1+a) + y*y > LAMBDA*LAMBDA); return ((x1-a)*(x1-a) + y*y > LAMBDA*LAMBDA); } else { width = 0.1; x1 = vabs(x + xshift); return((x1 - LAMBDA - width)*(x1 - LAMBDA - width) + y*y < LAMBDA*LAMBDA); } } case (D_TREE): { /* upper triangle */ h = 1.5*LAMBDA; r = 1.5*LAMBDA; x1 = vabs(x); y1 = y - h; zz[0][0] = r*cos(PI/6.0); zz[0][1] = -0.5*r; zz[1][0] = 0.0; zz[1][1] = r; zz[2][0] = 0.0; zz[2][1] = -0.5*r; if (y1 > 0.95*zz[1][1]) return(1); /* middle triangle */ h = 0.25*LAMBDA; r = 2.0*LAMBDA; y2 = y - h; zz[3][0] = r*cos(PI/6.0); zz[3][1] = -0.5*r; zz[4][0] = 0.0; zz[4][1] = r; zz[5][0] = 0.0; zz[5][1] = -0.5*r; /* bottom triangle */ h = -1.5*LAMBDA; r = 3.0*LAMBDA; y3 = y - h; zz[6][0] = r*cos(PI/6.0); zz[6][1] = -0.5*r; zz[7][0] = 0.0; zz[7][1] = r; zz[8][0] = 0.0; zz[8][1] = -0.5*r; if (xy_in_triangle(x1, y1, zz[0], zz[1], zz[2])) return(0); if (xy_in_triangle(x1, y2, zz[3], zz[4], zz[5])) return(0); if (xy_in_triangle(x1, y3, zz[6], zz[7], zz[8])) return(0); if (xy_in_triangle(0.9*x1, 0.9*y1, zz[0], zz[1], zz[2])) return(1); if (xy_in_triangle(0.9*x1, 0.9*y2, zz[3], zz[4], zz[5])) return(1); if (xy_in_triangle(0.9*x1, 0.9*y3, zz[6], zz[7], zz[8])) return(1); if (xy_in_triangle(0.8*x1, 0.8*y1, zz[0], zz[1], zz[2])) return(2); if (xy_in_triangle(0.8*x1, 0.8*y2, zz[3], zz[4], zz[5])) return(2); if (xy_in_triangle(0.8*x1, 0.8*y3, zz[6], zz[7], zz[8])) return(2); return(1); } case (D_MICHELSON): { if ((vabs(x) > LAMBDA)&&(vabs(y) > LAMBDA)) return(2); r = 0.5*sqrt(2.0); h = YMAX - 0.5*WALL_WIDTH; x1 = (x + y)*r; y1 = (-x + y)*r; if ((vabs(y1) < 0.5*WALL_WIDTH)&&(vabs(x1) < LAMBDA*sqrt(2.0))) return(1); if ((vabs(x - h) < 0.5*WALL_WIDTH)&&(vabs(y) < LAMBDA)) return(2); if ((vabs(x) < LAMBDA)&&(vabs(y - h) < 0.5*WALL_WIDTH)) return(2); return(0); } case (D_MICHELSON_MOVING): { if ((vabs(x) > LAMBDA)&&(vabs(y) > LAMBDA)) return(2); r = 0.5*sqrt(2.0); h = YMAX - MU - 0.5*WALL_WIDTH; h1 = YMAX - MU - 0.5*WALL_WIDTH + michelson_position; x1 = (x + y)*r; y1 = (-x + y)*r; if ((vabs(y1) < 0.5*WALL_WIDTH)&&(vabs(x1) < LAMBDA*sqrt(2.0))) return(1); if ((vabs(x - h1) < 0.5*WALL_WIDTH)&&(vabs(y) < LAMBDA)) return(2); if ((vabs(x) < LAMBDA)&&(vabs(y - h) < 0.5*WALL_WIDTH)) return(2); return(0); } case (D_RITCHEY_CHRETIEN): { /* LAMBDA is magnification M */ d = 2.5; b = 3.5; r1 = 2.0*(d + b/LAMBDA); r2 = 6.0/(LAMBDA - 1.0); x1 = 0.5*d - r1; x2 = -0.5*d - r2; y1 = vabs(y); if (x < 0.0) { if (y1 > 0.3) return(1); if (x < -0.5*d - WALL_WIDTH) return(1); return(module2(x-x2, y1) > r2); } else { if (y1 < MU) return(1); if (x > 0.5*d + WALL_WIDTH) return(1); return(module2(x-x1, y1) < r1); } } case (D_MENGER): { x1 = 0.5*(x+1.0); y1 = 0.5*(y+1.0); for (k=0; k l2)&&(r2 < 1.0)) return(1); else if (r2mu <= l2) return(2); else return (0); } case (D_MENGER_HEATED): { if ((vabs(x) >= 1.0)||(vabs(y) >= 1.0)) return(0); else { x1 = 0.5*(x+1.0); y1 = 0.5*(y+1.0); for (k=0; k 1.2) return(2); else return(1); } case (D_MANDELBROT_CIRCLE): { u = 0.0; v = 0.0; i = 0; while ((i 1.2) return(2); // else if ((vabs(x) > XMAX - 0.01)||(vabs(y) > YMAX - 0.01)) return(2); else return(1); } case (D_VONKOCH_HEATED): { if (x*x + y*y > LAMBDA*LAMBDA) return(2); x1 = x; y1 = y; condition = xy_in_triangle_tvertex(x1, y1, polyline[0], polyline[npolyline/3], polyline[2*npolyline/3]); m = 1; k = 1; for (i = 0; i < MDEPTH; i++) { m = m*4; for (j = 0; j < npolyline/m; j++) condition += xy_in_triangle_tvertex(x1, y1, polyline[j*m + k], polyline[j*m + 2*k], polyline[j*m + 3*k]); k = k*4; } if (condition > 0) return(0); else return(1); } default: { printf("Function ij_in_billiard not defined for this billiard \n"); return(0); } } } int xy_in_billiard(double x, double y) /* returns 1 if (x,y) represents a point in the billiard */ { if (COMPARISON) { if (y > 0.0) return (xy_in_billiard_single_domain(x, y, B_DOMAIN, ncircles, circles)); else return (xy_in_billiard_single_domain(x, y, B_DOMAIN_B, ncircles_b, circles_b)); } else return (xy_in_billiard_single_domain(x, y, B_DOMAIN, ncircles, circles)); } int ij_in_billiard(int i, int j) /* returns 1 if (i,j) represents a point in the billiard */ { double xy[2]; ij_to_xy(i, j, xy); return(xy_in_billiard(xy[0], xy[1])); } void tvertex_lineto(t_vertex z) /* draws boundary segments of isospectral billiard */ { glVertex2d(z.posi, z.posj); } void hex_transfo(double u, double v, double *x, double *y) /* linear transformation of plane used for hex tiles */ { static double ra, rb; static int first = 1; if (first) { ra = 0.5; rb = 0.5*sqrt(3.0); first = 0; } *x = u + ra*v; *y = rb*v; } void draw_billiard(int fade, double fade_value) /* draws the billiard boundary */ { double x0, y0, x, y, x1, y1, x2, y2, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z, l, width, a, b, c, d, r1, r2, ymax, height, xmin, xmax, ca, sa, xshift, x5, x6, f, fp, xratio, w; int i, j, k, k1, k2, mr2, ntiles; static int first = 1, nsides; static double h, hh, sqr3, ll, salpha, arcangle; char message[100]; if (fade) { if (BLACK) glColor3f(fade_value, fade_value, fade_value); else glColor3f(1.0 - fade_value, 1.0 - fade_value, 1.0 - fade_value); } else { if (BLACK) glColor3f(1.0, 1.0, 1.0); else glColor3f(0.0, 0.0, 0.0); } glLineWidth(BOUNDARY_WIDTH); glEnable(GL_LINE_SMOOTH); switch (B_DOMAIN) { case (D_RECTANGLE): { glBegin(GL_LINE_LOOP); xy_to_pos(LAMBDA, -1.0, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(LAMBDA, 1.0, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(-LAMBDA, 1.0, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(-LAMBDA, -1.0, pos); glVertex2d(pos[0], pos[1]); glEnd(); break; } case (D_ELLIPSE): { glBegin(GL_LINE_LOOP); for (i=0; i<=NSEG; i++) { phi = (double)i*DPI/(double)NSEG; x = LAMBDA*cos(phi); y = sin(phi); xy_to_pos(x, y, pos); glVertex2d(pos[0], pos[1]); } glEnd (); /* draw foci */ if (FOCI) { if (fade) glColor3f(0.3*fade_value, 0.3*fade_value, 0.3*fade_value); else glColor3f(0.3, 0.3, 0.3); x0 = sqrt(LAMBDA*LAMBDA-1.0); glLineWidth(2); glEnable(GL_LINE_SMOOTH); draw_circle(x0, 0.0, r, NSEG); draw_circle(-x0, 0.0, r, NSEG); } break; } case (D_EXT_ELLIPSE): { glBegin(GL_LINE_LOOP); for (i=0; i<=NSEG; i++) { phi = (double)i*DPI/(double)NSEG; x = LAMBDA*cos(phi); y = MU*sin(phi); xy_to_pos(x, y, pos); glVertex2d(pos[0], pos[1]); } glEnd (); /* draw foci */ if (FOCI) { if (fade) glColor3f(0.3*fade_value, 0.3*fade_value, 0.3*fade_value); else glColor3f(0.3, 0.3, 0.3); x0 = sqrt(LAMBDA*LAMBDA-MU*MU); glLineWidth(2); glEnable(GL_LINE_SMOOTH); draw_circle(x0, 0.0, r, NSEG); draw_circle(-x0, 0.0, r, NSEG); } break; } case (D_EXT_ELLIPSE_CURVED): { glBegin(GL_LINE_LOOP); for (i=0; i<=NSEG; i++) { phi = (double)i*DPI/(double)NSEG; x = LAMBDA*cos(phi); y = MU*sin(phi) - 0.4*x*x; xy_to_pos(x, y, pos); glVertex2d(pos[0], pos[1]); } glEnd (); break; } case (D_EXT_ELLIPSE_CURVED_BDRY): { glBegin(GL_LINE_LOOP); for (i=0; i<=NSEG; i++) { phi = (double)i*DPI/(double)NSEG; x = LAMBDA*cos(phi); y = MU*sin(phi) - 0.4*x*x; xy_to_pos(x, y, pos); glVertex2d(pos[0], pos[1]); } glEnd (); draw_line(XMIN, YMAX - 0.05, XMAX, YMAX - 0.05); draw_line(XMIN, YMIN + 0.05, XMAX, YMIN + 0.05); break; } case (D_STADIUM): { glBegin(GL_LINE_LOOP); for (i=0; i<=NSEG; i++) { phi = -PID + (double)i*PI/(double)NSEG; x = 0.5*LAMBDA + cos(phi); y = sin(phi); xy_to_pos(x, y, pos); glVertex2d(pos[0], pos[1]); } for (i=0; i<=NSEG; i++) { phi = PID + (double)i*PI/(double)NSEG; x = -0.5*LAMBDA + cos(phi); y = sin(phi); xy_to_pos(x, y, pos); glVertex2d(pos[0], pos[1]); } glEnd(); break; } case (D_SINAI): { draw_circle(0.0, 0.0, LAMBDA, NSEG); break; } case (D_DIAMOND): { alpha = atan(1.0 - 1.0/LAMBDA); dphi = (PID - 2.0*alpha)/(double)NSEG; r = sqrt(LAMBDA*LAMBDA + (LAMBDA-1.0)*(LAMBDA-1.0)); glBegin(GL_LINE_LOOP); for (i=0; i<=NSEG; i++) { phi = alpha + (double)i*dphi; x = -LAMBDA + r*cos(phi); y = -LAMBDA + r*sin(phi); xy_to_pos(x, y, pos); glVertex2d(pos[0], pos[1]); } for (i=0; i<=NSEG; i++) { phi = alpha - PID + (double)i*dphi; x = -LAMBDA + r*cos(phi); y = LAMBDA + r*sin(phi); xy_to_pos(x, y, pos); glVertex2d(pos[0], pos[1]); } for (i=0; i<=NSEG; i++) { phi = alpha + PI + (double)i*dphi; x = LAMBDA + r*cos(phi); y = LAMBDA + r*sin(phi); xy_to_pos(x, y, pos); glVertex2d(pos[0], pos[1]); } for (i=0; i<=NSEG; i++) { phi = alpha + PID + (double)i*dphi; x = LAMBDA + r*cos(phi); y = -LAMBDA + r*sin(phi); xy_to_pos(x, y, pos); glVertex2d(pos[0], pos[1]); } glEnd(); break; } case (D_TRIANGLE): { glBegin(GL_LINE_LOOP); xy_to_pos(-LAMBDA, -1.0, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(LAMBDA, -1.0, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(-LAMBDA, 1.0, pos); glVertex2d(pos[0], pos[1]); glEnd(); break; } case (D_FLAT): { glBegin(GL_LINE_LOOP); xy_to_pos(XMIN, -LAMBDA, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(XMAX, -LAMBDA, pos); glVertex2d(pos[0], pos[1]); glEnd(); break; } case (D_ANNULUS): { draw_circle(0.0, 0.0, LAMBDA, NSEG); draw_circle(0.0, 0.0, 1.0, NSEG); break; } case (D_POLYGON): { 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); xy_to_pos(x, y, pos); glVertex2d(pos[0], pos[1]); } glEnd (); break; } case (D_YOUNG): { glBegin(GL_LINE_STRIP); xy_to_pos(-MU, YMIN, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(-MU, -LAMBDA-MU, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(MU, -LAMBDA-MU, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(MU, YMIN, pos); glVertex2d(pos[0], pos[1]); glEnd(); glBegin(GL_LINE_STRIP); xy_to_pos(-MU, YMAX, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(-MU, LAMBDA+MU, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(MU, LAMBDA+MU, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(MU, YMAX, pos); glVertex2d(pos[0], pos[1]); glEnd(); glBegin(GL_LINE_LOOP); xy_to_pos(-MU, -LAMBDA+MU, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(-MU, LAMBDA-MU, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(MU, LAMBDA-MU, pos); glVertex2d(pos[0], pos[1]); xy_to_pos(MU, -LAMBDA+MU, pos); glVertex2d(pos[0], pos[1]); glEnd(); break; } case (D_GRATING): { k1 = -(int)(-YMIN/LAMBDA); k2 = (int)(YMAX/LAMBDA); for (i=k1; i<= k2; i++) { z = (double)i*LAMBDA; draw_circle(0.0, z, MU, NSEG); } break; } case (D_EHRENFEST): { alpha = asin(MU/LAMBDA); x0 = 1.0 - sqrt(LAMBDA*LAMBDA - MU*MU); dphi = 2.0*(PI-alpha)/((double)NSEG); glBegin(GL_LINE_LOOP); for (i=0; i<=NSEG; i++) { phi = -PI + alpha + (double)i*dphi; x = 1.0 + LAMBDA*cos(phi); y = LAMBDA*sin(phi); xy_to_pos(x, y, pos); glVertex2d(pos[0], pos[1]); } for (i=0; i<=NSEG; i++) { phi = alpha + (double)i*dphi; x = -1.0 + LAMBDA*cos(phi); y = LAMBDA*sin(phi); xy_to_pos(x, y, pos); glVertex2d(pos[0], pos[1]); } glEnd (); break; } case (D_DISK_GRID): { glLineWidth(2); for (i = -NGRIDX/2; i < NGRIDX/2; i++) for (j = 0; j < NGRIDY; j++) { dy = (YMAX - YMIN)/((double)NGRIDY); dx = dy*0.5*sqrt(3.0); x1 = ((double)i + 0.5)*dy; y1 = YMIN + ((double)j + 0.5)*dy; draw_circle(x1, y1, MU, NSEG); } break; } case (D_DISK_HEX): { glLineWidth(2); for (i = -NGRIDX/2; i < NGRIDX/2; i++) for (j = -1; j < NGRIDY; j++) { dy = (YMAX - YMIN)/((double)NGRIDY); x1 = ((double)i + 0.5)*dy; y1 = YMIN + ((double)j + 0.5)*dy; if ((i+NGRIDX)%2 == 1) y1 += 0.5*dy; draw_circle(x1, y1, MU, NSEG); } break; } case (D_PARABOLA): { dy = (YMAX - YMIN)/(double)NSEG; glBegin(GL_LINE_STRIP); for (i = 0; i < NSEG+1; i++) { y = YMIN + dy*(double)i; x = 0.25*y*y/LAMBDA - LAMBDA; xy_to_pos(x, y, pos); glVertex2d(pos[0], pos[1]); } glEnd (); if (FOCI) { glColor3f(0.3, 0.3, 0.3); draw_circle(0.0, 0.0, r, NSEG); } break; } case (D_TWO_PARABOLAS): { dy = 3.0*MU/(double)NSEG; width = 0.25*MU; if (width > 0.2) width = 0.2; glBegin(GL_LINE_LOOP); for (i = 0; i < NSEG+1; i++) { y = -1.5*MU + dy*(double)i; x = 0.25*y*y/MU - MU - LAMBDA; xy_to_pos(x, y, pos); glVertex2d(pos[0], pos[1]); } for (i = 0; i < NSEG+1; i++) { y = 1.5*MU - dy*(double)i; x = 0.25*y*y/MU - (MU + width) - LAMBDA; xy_to_pos(x, y, pos); glVertex2d(pos[0], pos[1]); } glEnd (); glBegin(GL_LINE_LOOP); for (i = 0; i < NSEG+1; i++) { y = -1.5*MU + dy*(double)i; x = LAMBDA + MU - 0.25*y*y/MU; xy_to_pos(x, y, pos); glVertex2d(pos[0], pos[1]); } for (i = 0; i < NSEG+1; i++) { y = 1.5*MU - dy*(double)i; x = LAMBDA + (MU + width) - 0.25*y*y/MU; xy_to_pos(x, y, pos); glVertex2d(pos[0], pos[1]); } glEnd (); if (FOCI) { glColor3f(0.3, 0.3, 0.3); draw_circle(-LAMBDA, 0.0, r, NSEG); draw_circle(LAMBDA, 0.0, r, NSEG); } break; } case (D_FOUR_PARABOLAS): { x1 = 2.0*(sqrt(MU*(2.0*MU + LAMBDA)) - MU); dy = 2.0*x1/(double)NSEG; glBegin(GL_LINE_LOOP); for (i = 0; i < NSEG+1; i++) { y = -x1 + dy*(double)i; x = MU + LAMBDA - 0.25*y*y/MU; xy_to_pos(x, y, pos); glVertex2d(pos[0], pos[1]); } for (i = 0; i < NSEG+1; i++) { x = x1 - dy*(double)i; y = MU + LAMBDA - 0.25*x*x/MU; xy_to_pos(x, y, pos); glVertex2d(pos[0], pos[1]); } for (i = 0; i < NSEG+1; i++) { y = x1 - dy*(double)i; x = -MU - LAMBDA + 0.25*y*y/MU; xy_to_pos(x, y, pos); glVertex2d(pos[0], pos[1]); } for (i = 0; i < NSEG+1; i++) { x = -x1 + dy*(double)i; y = -MU - LAMBDA + 0.25*x*x/MU; xy_to_pos(x, y, pos); glVertex2d(pos[0], pos[1]); } glEnd (); if (FOCI) { glColor3f(0.3, 0.3, 0.3); draw_circle(-LAMBDA, 0.0, r, NSEG); draw_circle(LAMBDA, 0.0, r, NSEG); draw_circle(0.0, -LAMBDA, r, NSEG); draw_circle(0.0, LAMBDA, r, NSEG); } break; } case (D_POLY_PARABOLAS): { omega = PI/((double)NPOLY); a = 0.25/MU; b = 1.0/tan(omega); c = LAMBDA + MU; ymax = (-b + sqrt(b*b + 4.0*a*c))/(2.0*a); dy = 2.0*ymax/(double)NSEG; // printf("a = %.3lg, b = %.3lg, ymax = %.3lg\n", a, b,ymax); glBegin(GL_LINE_LOOP); for (k=0; k45; i--) tvertex_lineto(polyline[i]); glEnd(); /* inner lines */ // glLineWidth(BOUNDARY_WIDTH/2); glLineWidth(1); glColor3f(0.75, 0.75, 0.75); glBegin(GL_LINE_STRIP); tvertex_lineto(polyline[0]); tvertex_lineto(polyline[1]); tvertex_lineto(polyline[2]); tvertex_lineto(polyline[0]); tvertex_lineto(polyline[3]); tvertex_lineto(polyline[4]); glEnd(); glBegin(GL_LINE_STRIP); tvertex_lineto(polyline[0]); tvertex_lineto(polyline[44]); tvertex_lineto(polyline[45]); tvertex_lineto(polyline[0]); tvertex_lineto(polyline[46]); tvertex_lineto(polyline[45]); glEnd(); for (i=3; i<43; i++) { glBegin(GL_LINE_STRIP); tvertex_lineto(polyline[i]); tvertex_lineto(polyline[43]); glEnd(); glBegin(GL_LINE_STRIP); tvertex_lineto(polyline[i+42]); tvertex_lineto(polyline[85]); glEnd(); } break; } case (D_ISOSPECTRAL): { /* 1st triangle */ glBegin(GL_LINE_LOOP); tvertex_lineto(polyline[0]); tvertex_lineto(polyline[4]); tvertex_lineto(polyline[7]); tvertex_lineto(polyline[1]); tvertex_lineto(polyline[5]); tvertex_lineto(polyline[8]); tvertex_lineto(polyline[2]); tvertex_lineto(polyline[3]); tvertex_lineto(polyline[6]); glEnd(); /* inner lines */ glBegin(GL_LINE_LOOP); tvertex_lineto(polyline[0]); tvertex_lineto(polyline[1]); tvertex_lineto(polyline[2]); tvertex_lineto(polyline[0]); tvertex_lineto(polyline[3]); tvertex_lineto(polyline[2]); tvertex_lineto(polyline[5]); tvertex_lineto(polyline[1]); tvertex_lineto(polyline[4]); glEnd(); /* 2nd triangle */ glBegin(GL_LINE_LOOP); tvertex_lineto( polyline[9]); tvertex_lineto(polyline[16]); tvertex_lineto(polyline[13]); tvertex_lineto(polyline[10]); tvertex_lineto(polyline[17]); tvertex_lineto(polyline[14]); tvertex_lineto(polyline[11]); tvertex_lineto(polyline[15]); tvertex_lineto(polyline[12]); glEnd(); /* inner lines */ glBegin(GL_LINE_LOOP); tvertex_lineto( polyline[9]); tvertex_lineto(polyline[10]); tvertex_lineto(polyline[11]); tvertex_lineto( polyline[9]); tvertex_lineto(polyline[13]); tvertex_lineto(polyline[10]); tvertex_lineto(polyline[14]); tvertex_lineto(polyline[11]); tvertex_lineto(polyline[12]); glEnd(); break; } case (D_HOMOPHONIC): { /* 1st triangle */ glBegin(GL_LINE_LOOP); tvertex_lineto(polyline[1]); tvertex_lineto(polyline[3]); tvertex_lineto(polyline[4]); tvertex_lineto(polyline[5]); tvertex_lineto(polyline[6]); tvertex_lineto(polyline[8]); tvertex_lineto(polyline[9]); tvertex_lineto(polyline[10]); tvertex_lineto(polyline[12]); tvertex_lineto(polyline[13]); tvertex_lineto(polyline[15]); tvertex_lineto(polyline[16]); tvertex_lineto(polyline[17]); tvertex_lineto(polyline[18]); tvertex_lineto(polyline[20]); glEnd(); /* inner lines */ glLineWidth(BOUNDARY_WIDTH/2); glBegin(GL_LINE_STRIP); tvertex_lineto(polyline[9]); tvertex_lineto(polyline[1]); tvertex_lineto(polyline[2]); tvertex_lineto(polyline[5]); tvertex_lineto(polyline[7]); tvertex_lineto(polyline[2]); tvertex_lineto(polyline[8]); tvertex_lineto(polyline[21]); tvertex_lineto(polyline[10]); tvertex_lineto(polyline[2]); tvertex_lineto(polyline[21]); tvertex_lineto(polyline[11]); tvertex_lineto(polyline[13]); tvertex_lineto(polyline[21]); tvertex_lineto(polyline[14]); tvertex_lineto(polyline[20]); tvertex_lineto(polyline[15]); tvertex_lineto(polyline[19]); tvertex_lineto(polyline[16]); tvertex_lineto(polyline[18]); glEnd(); /* 2nd triangle */ glLineWidth(BOUNDARY_WIDTH); glBegin(GL_LINE_LOOP); tvertex_lineto(polyline[22+10]); tvertex_lineto(polyline[22+16]); tvertex_lineto(polyline[22+17]); tvertex_lineto(polyline[22+18]); tvertex_lineto(polyline[22+12]); tvertex_lineto(polyline[22+13]); tvertex_lineto(polyline[22+15]); tvertex_lineto(polyline[22+19]); tvertex_lineto(polyline[22+20]); tvertex_lineto(polyline[22+1]); tvertex_lineto(polyline[22+4]); tvertex_lineto(polyline[22+5]); tvertex_lineto(polyline[22+7]); tvertex_lineto(polyline[22+8]); tvertex_lineto(polyline[22+9]); glEnd(); /* inner lines */ glLineWidth(BOUNDARY_WIDTH/2); glBegin(GL_LINE_STRIP); tvertex_lineto(polyline[22+2]); tvertex_lineto(polyline[22+6]); tvertex_lineto(polyline[22+8]); tvertex_lineto(polyline[22+2]); tvertex_lineto(polyline[22+5]); tvertex_lineto(polyline[22+3]); tvertex_lineto(polyline[22+2]); tvertex_lineto(polyline[22+1]); tvertex_lineto(polyline[22+0]); tvertex_lineto(polyline[22+21]); tvertex_lineto(polyline[22+18]); tvertex_lineto(polyline[22+16]); tvertex_lineto(polyline[22+13]); tvertex_lineto(polyline[22+21]); tvertex_lineto(polyline[22+10]); tvertex_lineto(polyline[22+12]); tvertex_lineto(polyline[22+21]); tvertex_lineto(polyline[22+14]); tvertex_lineto(polyline[22+20]); tvertex_lineto(polyline[22+15]); glEnd(); break; } case (D_VONKOCH): { glLineWidth(BOUNDARY_WIDTH/2); glBegin(GL_LINE_LOOP); for (i=0; i 0.0) x = sqrt(LAMBDA*LAMBDA - y*y) - LAMBDA + MU; else x = -sqrt(LAMBDA*LAMBDA - y*y) - LAMBDA - MU; xy_to_pos(x, y, pos); glVertex2d(pos[0], pos[1]); } y = 0.9*LAMBDA; if (LAMBDA > 0.0) x = sqrt(LAMBDA*LAMBDA - y*y) - LAMBDA + MU; else x = -sqrt(LAMBDA*LAMBDA - y*y) - LAMBDA - MU; xy_to_pos(x, y, pos); glVertex2d(pos[0], pos[1]); glEnd(); break; } case (D_CIRCLES): { glLineWidth(BOUNDARY_WIDTH); for (i = 0; i < ncircles; i++) if (circles[i].active) draw_circle(circles[i].xc, circles[i].yc, circles[i].radius, NSEG); break; } case (D_CIRCLES_IN_RECT): { glLineWidth(BOUNDARY_WIDTH); for (i = 0; i < ncircles; i++) if (circles[i].active) draw_circle(circles[i].xc, circles[i].yc, circles[i].radius, NSEG); draw_rectangle(-LAMBDA, -1.0, LAMBDA, 1.0); if ((FOCI)&&(CIRCLE_PATTERN == C_LASER)) { glColor3f(0.3, 0.3, 0.3); draw_circle(X_SHOOTER, Y_SHOOTER, r, NSEG); } break; } case (D_POLYGONS): { glLineWidth(BOUNDARY_WIDTH); for (i = 0; i < ncircles; i++) if (polygons[i].active) draw_tpolygon(polygons[i]); break; } case (D_MAZE): { glLineWidth(BOUNDARY_WIDTH); if (fade) glColor3f(0.15*fade_value, 0.15*fade_value, 0.15*fade_value); else glColor3f(0.15, 0.15, 0.15); for (i=0; i XMIN) { draw_line(x, YMIN, x, YMAX); x -= LAMBDA; } y = 0.5*LAMBDA; while (y < YMAX) { draw_line(XMIN, y, XMAX, y); y += LAMBDA; } y = -0.5*LAMBDA; while (y > YMIN) { draw_line(XMIN, y, XMAX, y); y -= LAMBDA; } break; } case (D_TRIANGLE_TILES): { if (first) { h = LAMBDA/(2.0*sqrt(3.0)); hh = h*3.0; sqr3 = sqrt(3.0); first = 0; } glLineWidth(BOUNDARY_WIDTH); y = -h; while (y < YMAX) { draw_line(XMIN, y, XMAX, y); y += hh; } y = -h; while (y > YMIN) { draw_line(XMIN, y, XMAX, y); y -= hh; } x = -0.5*LAMBDA; y = -h; while (x < 1.5*XMAX) { draw_line(x - 10.0, y - 10.0*sqr3, x + 10.0, y + 10.0*sqr3); draw_line(x - 10.0, y + 10.0*sqr3, x + 10.0, y - 10.0*sqr3); x += LAMBDA; } x = -0.5*LAMBDA; while (x > 1.5*XMIN) { draw_line(x - 10.0, y - 10.0*sqr3, x + 10.0, y + 10.0*sqr3); draw_line(x - 10.0, y + 10.0*sqr3, x + 10.0, y - 10.0*sqr3); x -= LAMBDA; } break; } case (D_HEX_TILES): { ntiles = (int)(XMAX/LAMBDA) + 1; for (i=-ntiles; i 0) draw_line(x1, y1, x, y); x1 = x; y1 = y; } } } break; } case (D_ONE_FUNNEL): { for (k=-1; k<2; k+=2) { x1 = XMIN; y1 = (double)k*(MU + LAMBDA*x1*x1); for (i=0; i<=NSEG; i++) { x = XMIN + (XMAX - XMIN)*(double)i/(double)NSEG; y = (double)k*(MU + LAMBDA*x*x); draw_line(x1, y1, x, y); x1 = x; y1 = y; } } break; } case (D_LENSES_RING): { if (first) { salpha = DPI/(double)NPOLY; h = LAMBDA*tan(PI/(double)NPOLY); if (h < MU) ll = sqrt(MU*MU - h*h); else ll = 0.0; arcangle = atan(h/ll); first = 0; } for (i=0; i 0) { glLineWidth(2); x = 1.0/((double)MRATIO); draw_rectangle(x, x, -x, -x); } /* level 2 */ if (MDEPTH > 1) { glLineWidth(1); mr2 = MRATIO*MRATIO; l = 2.0/((double)mr2); for (i=0; i 2) { glLineWidth(1); l = 2.0/((double)(mr2*MRATIO)); for (i=0; i 0) { glLineWidth(2); x = 1.0/((double)MRATIO); draw_rotated_rectangle(x, x, -x, -x); } /* level 2 */ if (MDEPTH > 1) { glLineWidth(1); mr2 = MRATIO*MRATIO; l = 2.0/((double)mr2); for (i=0; i 2) { glLineWidth(1); l = 2.0/((double)(mr2*MRATIO)); for (i=0; i 0) { glLineWidth(2); x = 1.0/((double)MRATIO); draw_rectangle(x, x, -x, -x); } /* level 2 */ if (MDEPTH > 1) { glLineWidth(1); mr2 = MRATIO*MRATIO; l = 2.0/((double)mr2); for (i=0; i 2) { glLineWidth(1); l = 2.0/((double)(mr2*MRATIO)); for (i=0; i 0) { glLineWidth(2); x = 1.0/((double)MRATIO); draw_rectangle(x, x, -x, -x); } /* level 2 */ if (MDEPTH > 1) { glLineWidth(1); mr2 = MRATIO*MRATIO; l = 2.0/((double)mr2); for (i=0; i 2) { glLineWidth(1); l = 2.0/((double)(mr2*MRATIO)); for (i=0; i= COL_TURBO) color_scheme_asym(COLOR_SCHEME, value, 1.0, 1, rgb); else color_scheme(COLOR_SCHEME, value, 1.0, 1, rgb); break; } case (P_MEAN_ENERGY): { value = dy_e*(double)(j - jmin)*100.0/E_SCALE; if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym(COLOR_SCHEME, value, 1.0, 1, rgb); else color_scheme(COLOR_SCHEME, value, 1.0, 1, rgb); break; } case (P_LOG_ENERGY): { value = LOG_SHIFT + LOG_SCALE*log(dy_e*(double)(j - jmin)*100.0/E_SCALE); // if (value <= 0.0) value = 0.0; color_scheme(COLOR_SCHEME, value, 1.0, 1, rgb); break; } case (P_LOG_MEAN_ENERGY): { value = LOG_SHIFT + LOG_SCALE*log(dy_e*(double)(j - jmin)*100.0/E_SCALE); // if (value <= 0.0) value = 0.0; color_scheme(COLOR_SCHEME, value, 1.0, 1, rgb); break; } case (P_ENERGY_FLUX): { value = dy_e*(double)(j - jmin)*100.0/E_SCALE; if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, COLOR_PALETTE, value, 1.0, 1, rgb); else color_scheme_palette(COLOR_SCHEME, COLOR_PALETTE, value, 1.0, 1, rgb); // value = min + 1.0*dy*(double)(j - jmin); // amp = 0.7*color_amplitude_linear(value, 1.0, 1); // while (amp > 1.0) amp -= 2.0; // while (amp < -1.0) amp += 2.0; // amp_to_rgb(0.5*(1.0 + amp), rgb); break; } case (P_TOTAL_ENERGY_FLUX): { // value = min + 1.0*dy*(double)(j - jmin); // amp = 0.7*color_amplitude_linear(value, 1.0, 1); // while (amp > 1.0) amp -= 2.0; // while (amp < -1.0) amp += 2.0; // amp_to_rgb(0.5*(1.0 + amp), rgb); value = dy_phase*(double)(j - jmin); color_scheme_palette(C_ONEDIM_LINEAR, COLOR_PALETTE, value, 1.0, 1, rgb); break; } case (P_PHASE): { value = min + 1.0*dy*(double)(j - jmin); // lum = (color_amplitude(value, 1.0, 1))*0.5; // if (lum < 0.0) lum = 0.0; // hsl_to_rgb(value*360.0, 0.9, 0.5, rgb); // color_scheme(COLOR_SCHEME, value, 1.0, 1, rgb); // amp = color_amplitude_linear(value, 1.0, 1); amp = 0.5*color_amplitude_linear(value, 1.0, 1); while (amp > 1.0) amp -= 2.0; while (amp < -1.0) amp += 2.0; amp_to_rgb(0.5*(1.0 + amp), rgb); break; } } glColor3f(rgb[0], rgb[1], rgb[2]); if (ROTATE_COLOR_SCHEME) { glVertex2i(j, imin); glVertex2i(j, imax); glVertex2i(j+1, imax); glVertex2i(j+1, imin); } else { glVertex2i(imin, j); glVertex2i(imax, j); glVertex2i(imax, j+1); glVertex2i(imin, j+1); } } glEnd (); glColor3f(1.0, 1.0, 1.0); glLineWidth(BOUNDARY_WIDTH); draw_rectangle(x1, y1, x2, y2); } void draw_color_scheme_palette(double x1, double y1, double x2, double y2, int plot, double min, double max, int palette) { int j, k, ij_botleft[2], ij_topright[2], imin, imax, jmin, jmax; double y, dy, dy_e, rgb[3], value, lum, amp, dy_phase; xy_to_ij(x1, y1, ij_botleft); xy_to_ij(x2, y2, ij_topright); rgb[0] = 0.0; rgb[1] = 0.0; rgb[2] = 0.0; // erase_area_rgb(0.5*(x1 + x2), x2 - x1, 0.5*(y1 + y2), y2 - y1, rgb); if (ROTATE_COLOR_SCHEME) { jmin = ij_botleft[0]; jmax = ij_topright[0]; imin = ij_botleft[1]; imax = ij_topright[1]; } else { imin = ij_botleft[0]; imax = ij_topright[0]; jmin = ij_botleft[1]; jmax = ij_topright[1]; } glBegin(GL_QUADS); dy = (max - min)/((double)(jmax - jmin)); dy_e = max/((double)(jmax - jmin)); dy_phase = 1.0/((double)(jmax - jmin)); for (j = jmin; j < jmax; j++) { switch (plot) { case (P_AMPLITUDE): { value = min + 1.0*dy*(double)(j - jmin); color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); break; } case (P_ENERGY): { value = dy_e*(double)(j - jmin)*100.0/E_SCALE; if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); break; } case (P_MEAN_ENERGY): { value = dy_e*(double)(j - jmin)*100.0/E_SCALE; if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); break; } case (P_LOG_ENERGY): { value = LOG_SCALE*log(dy_e*(double)(j - jmin)*100.0/E_SCALE); // if (value <= 0.0) value = 0.0; color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); break; } case (P_LOG_MEAN_ENERGY): { value = LOG_SHIFT + LOG_SCALE*log(dy_e*(double)(j - jmin)*100.0/E_SCALE); // if (value <= 0.0) value = 0.0; color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); break; } case (P_ENERGY_FLUX): { value = dy_e*(double)(j - jmin)*100.0/E_SCALE; if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); // value = min + 1.0*dy*(double)(j - jmin); // amp = 0.7*color_amplitude_linear(value, 1.0, 1); // while (amp > 1.0) amp -= 2.0; // while (amp < -1.0) amp += 2.0; // amp_to_rgb(0.5*(1.0 + amp), rgb); break; } case (P_TOTAL_ENERGY_FLUX): { // value = min + 1.0*dy*(double)(j - jmin); // amp = 0.7*color_amplitude_linear(value, 1.0, 1); // while (amp > 1.0) amp -= 2.0; // while (amp < -1.0) amp += 2.0; // amp_to_rgb(0.5*(1.0 + amp), rgb); value = dy_phase*(double)(j - jmin); color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb); break; } case (P_PHASE): { value = min + 1.0*dy*(double)(j - jmin); // lum = (color_amplitude(value, 1.0, 1))*0.5; // if (lum < 0.0) lum = 0.0; // hsl_to_rgb(value*360.0, 0.9, 0.5, rgb); // color_scheme(COLOR_SCHEME, value, 1.0, 1, rgb); // amp = color_amplitude_linear(value, 1.0, 1); amp = 0.5*color_amplitude_linear(value, 1.0, 1); while (amp > 1.0) amp -= 2.0; while (amp < -1.0) amp += 2.0; amp_to_rgb(0.5*(1.0 + amp), rgb); break; } } glColor3f(rgb[0], rgb[1], rgb[2]); if (ROTATE_COLOR_SCHEME) { glVertex2i(j, imin); glVertex2i(j, imax); glVertex2i(j+1, imax); glVertex2i(j+1, imin); } else { glVertex2i(imin, j); glVertex2i(imax, j); glVertex2i(imax, j+1); glVertex2i(imin, j+1); } } glEnd (); glColor3f(1.0, 1.0, 1.0); glLineWidth(BOUNDARY_WIDTH); draw_rectangle(x1, y1, x2, y2); } void draw_color_scheme_palette_fade(double x1, double y1, double x2, double y2, int plot, double min, double max, int palette, int fade, double fade_value) { int j, k, ij_botleft[2], ij_topright[2], imin, imax, jmin, jmax; double y, dy, dy_e, rgb[3], value, lum, amp, dy_phase; xy_to_ij(x1, y1, ij_botleft); xy_to_ij(x2, y2, ij_topright); rgb[0] = 0.0; rgb[1] = 0.0; rgb[2] = 0.0; // erase_area_rgb(0.5*(x1 + x2), x2 - x1, 0.5*(y1 + y2), y2 - y1, rgb); if (ROTATE_COLOR_SCHEME) { jmin = ij_botleft[0]; jmax = ij_topright[0]; imin = ij_botleft[1]; imax = ij_topright[1]; } else { imin = ij_botleft[0]; imax = ij_topright[0]; jmin = ij_botleft[1]; jmax = ij_topright[1]; } glBegin(GL_QUADS); dy = (max - min)/((double)(jmax - jmin)); dy_e = max/((double)(jmax - jmin)); dy_phase = 1.0/((double)(jmax - jmin)); for (j = jmin; j < jmax; j++) { switch (plot) { case (P_AMPLITUDE): { value = min + 1.0*dy*(double)(j - jmin); color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); break; } case (P_ENERGY): { value = dy_e*(double)(j - jmin)*100.0/E_SCALE; if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); break; } case (P_MEAN_ENERGY): { value = dy_e*(double)(j - jmin)*100.0/E_SCALE; if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); break; } case (P_LOG_ENERGY): { value = LOG_SCALE*log(dy_e*(double)(j - jmin)*100.0/E_SCALE); // printf("value = %.2lg\n", value); // if (value <= 0.0) value = 0.0; color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); break; } case (P_LOG_MEAN_ENERGY): { value = LOG_SHIFT + LOG_SCALE*log(dy_e*(double)(j - jmin)*100.0/E_SCALE); // if (value <= 0.0) value = 0.0; color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); break; } case (P_ENERGY_FLUX): { value = dy_e*(double)(j - jmin); if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); // value = min + 1.0*dy*(double)(j - jmin); // amp = 0.7*color_amplitude_linear(value, 1.0, 1); // while (amp > 1.0) amp -= 2.0; // while (amp < -1.0) amp += 2.0; // amp_to_rgb(0.5*(1.0 + amp), rgb); break; } case (P_TOTAL_ENERGY_FLUX): { // value = min + 1.0*dy*(double)(j - jmin); // amp = 0.7*color_amplitude_linear(value, 1.0, 1); // while (amp > 1.0) amp -= 2.0; // while (amp < -1.0) amp += 2.0; // amp_to_rgb(0.5*(1.0 + amp), rgb); value = dy_phase*(double)(j - jmin); color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb); break; } case (P_AVERAGE_ENERGY): { value = dy_e*(double)(j - jmin)*100.0/E_SCALE; if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); break; } case (P_LOG_AVERAGE_ENERGY): { value = LOG_SCALE*log(dy_e*(double)(j - jmin)*100.0/E_SCALE); color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); break; } case (P_PHASE): { value = min + 1.0*dy*(double)(j - jmin); // lum = (color_amplitude(value, 1.0, 1))*0.5; // if (lum < 0.0) lum = 0.0; // hsl_to_rgb(value*360.0, 0.9, 0.5, rgb); // color_scheme(COLOR_SCHEME, value, 1.0, 1, rgb); // amp = color_amplitude_linear(value, 1.0, 1); amp = 0.5*color_amplitude_linear(value, 1.0, 1); while (amp > 1.0) amp -= 2.0; while (amp < -1.0) amp += 2.0; amp_to_rgb(0.5*(1.0 + amp), rgb); break; } case (Z_EULER_VORTICITY): { value = min + 1.0*dy*(double)(j - jmin); color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); break; } case (Z_EULER_LOG_VORTICITY): { value = min + 1.0*dy*(double)(j - jmin); color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); break; } case (Z_EULER_VORTICITY_ASYM): { value = min + 1.0*dy*(double)(j - jmin); color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); break; } case (Z_EULER_LPRESSURE): { value = min + 1.0*dy*(double)(j - jmin); color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); break; } case (Z_EULER_PRESSURE): { value = min + 1.0*dy*(double)(j - jmin); color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); break; } case (Z_EULER_DENSITY): { value = min + 1.0*dy*(double)(j - jmin); color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); break; } case (Z_EULER_SPEED): { value = min + 1.0*dy*(double)(j - jmin); color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); break; } case (Z_EULERC_VORTICITY): { value = min + 1.0*dy*(double)(j - jmin); color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); break; } case (P_3D_LOG_MEAN_ENERGY): { value = LOG_SCALE*dy_e*(double)(j - jmin)*100.0/E_SCALE; color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); break; } default: { value = min + 1.0*dy*(double)(j - jmin); color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); break; } } if (fade) for (k=0; k<3; k++) rgb[k] *= fade_value; glColor3f(rgb[0], rgb[1], rgb[2]); if (ROTATE_COLOR_SCHEME) { glVertex2i(j, imin); glVertex2i(j, imax); glVertex2i(j+1, imax); glVertex2i(j+1, imin); } else { glVertex2i(imin, j); glVertex2i(imax, j); glVertex2i(imax, j+1); glVertex2i(imin, j+1); } } glEnd (); if (fade) glColor3f(fade_value, fade_value, fade_value); else glColor3f(1.0, 1.0, 1.0); glLineWidth(BOUNDARY_WIDTH); draw_rectangle(x1, y1, x2, y2); } void draw_circular_color_scheme_palette_fade(double x1, double y1, double radius, int plot, double min, double max, int palette, int fade, double fade_value) { int j, k; double x, y, dy, dy_e, dy_phase, rgb[3], value, lum, amp, dphi, pos[2], phi, xy[2], zscale = 0.9; // glBegin(GL_TRIANGLE_FAN); xy_to_pos(x1, y1, xy); glVertex2d(xy[0], xy[1]); dy = (max - min)/360.0; dy_e = max/360.0; dy_phase = 1.0/360.0; dphi = DPI/360.0; for (j = 0; j <= 360; j++) { switch (plot) { case (P_AMPLITUDE): { value = min + 1.0*dy*(double)(j); color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); break; } case (P_ENERGY): { value = dy_e*(double)(j)*100.0/E_SCALE; if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); break; } case (P_MEAN_ENERGY): { value = dy_e*(double)(j)*100.0/E_SCALE; if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); break; } case (P_LOG_ENERGY): { value = LOG_SCALE*log(dy_e*(double)(j)*100.0/E_SCALE); color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); break; } case (P_LOG_MEAN_ENERGY): { value = LOG_SHIFT + LOG_SCALE*log(dy_e*(double)(j)*100.0/E_SCALE); color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); break; } case (P_ENERGY_FLUX): { value = dy_e*(double)(j); if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); break; } case (P_TOTAL_ENERGY_FLUX): { value = dy_phase*(double)(j); color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb); break; } case (P_AVERAGE_ENERGY): { value = dy_e*(double)(j)*100.0/E_SCALE; if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); break; } case (P_LOG_AVERAGE_ENERGY): { value = LOG_SCALE*log(dy_e*(double)(j)*100.0/E_SCALE); color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); break; } case (P_PHASE): { value = min + 1.0*dy*(double)(j); amp = 0.5*color_amplitude_linear(value, 1.0, 1); while (amp > 1.0) amp -= 2.0; while (amp < -1.0) amp += 2.0; amp_to_rgb(0.5*(1.0 + amp), rgb); break; } case (Z_EULER_VORTICITY): { value = min + 1.0*dy*(double)(j); color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); break; } case (Z_EULER_LOG_VORTICITY): { value = min + 1.0*dy*(double)(j); color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); break; } case (Z_EULER_VORTICITY_ASYM): { value = min + 1.0*dy*(double)(j); color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); break; } case (Z_EULER_LPRESSURE): { value = min + 1.0*dy*(double)(j); color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); break; } case (Z_EULER_PRESSURE): { value = min + 1.0*dy*(double)(j); color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); break; } case (Z_EULER_DENSITY): { value = min + 1.0*dy*(double)(j); color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); break; } case (Z_EULER_SPEED): { value = min + 1.0*dy*(double)(j); color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); break; } case (Z_EULERC_VORTICITY): { value = min + 1.0*dy*(double)(j); color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); break; } default: { value = min + 1.0*dy*(double)(j); color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); break; } } if (fade) for (k=0; k<3; k++) rgb[k] *= fade_value; glColor3f(rgb[0], rgb[1], rgb[2]); glBegin(GL_TRIANGLE_FAN); xy_to_pos(x1, y1, xy); glVertex2d(xy[0], xy[1]); xy_to_pos(x1 + radius*cos(dphi*(double)(j-1)), y1 + zscale*radius*sin(dphi*(double)(j-1)), xy); glVertex2d(xy[0], xy[1]); xy_to_pos(x1 + radius*cos(dphi*(double)j), y1 + zscale*radius*sin(dphi*(double)j), xy); glVertex2d(xy[0], xy[1]); glEnd (); } xy_to_pos(x1 + radius*cos(dphi), y1 + radius*sin(dphi), xy); // glVertex2d(xy[0], xy[1]); // glEnd (); if (fade) glColor3f(fade_value, fade_value, fade_value); else glColor3f(1.0, 1.0, 1.0); glLineWidth(BOUNDARY_WIDTH*3/2); glEnable(GL_LINE_SMOOTH); dphi = DPI/(double)NSEG; glBegin(GL_LINE_LOOP); for (j = 0; j < NSEG; j++) { xy_to_pos(x1 + radius*cos(dphi*(double)j), y1 + zscale*radius*sin(dphi*(double)j), xy); glVertex2d(xy[0], xy[1]); glVertex2d(xy[0], xy[1]); } glEnd (); } void print_speed(double speed, int fade, double fade_value) { char message[100]; double y = YMAX - 0.1, pos[2]; static double xleftbox, xlefttext; static int first = 1; if (first) { xleftbox = XMIN + 0.3; xlefttext = xleftbox - 0.45; first = 0; } erase_area_hsl(xleftbox, y + 0.025, 0.22, 0.05, 0.0, 0.9, 0.0); if (fade) glColor3f(fade_value, fade_value, fade_value); else glColor3f(1.0, 1.0, 1.0); xy_to_pos(xlefttext + 0.28, y, pos); sprintf(message, "Mach %.3lg", speed); write_text(pos[0], pos[1], message); } void print_frequency(double phase_shift, int fade, double fade_value) { char message[100]; double y = YMAX - 0.1, pos[2]; static double xleftbox, xlefttext; static int first = 1; if (first) { xleftbox = XMIN + 0.3; xlefttext = xleftbox - 0.45; first = 0; } erase_area_hsl(xleftbox, y + 0.025, 0.22, 0.05, 0.0, 0.9, 0.0); if (fade) glColor3f(fade_value, fade_value, fade_value); else glColor3f(1.0, 1.0, 1.0); xy_to_pos(xlefttext + 0.28, y, pos); sprintf(message, "Frequency %.2f", 25.0*phase_shift); write_text(pos[0], pos[1], message); } void init_laplacian_coords(t_laplacian laplace[NX*NY], double phi[NX*NY]) /* compute coordinates of neighbours to compute Laplacian */ { int i, j, iplus, iminus, i1, i2, i3, j1, j2, j3, ij[2];; printf("Initialising Laplacian table\n"); /* Laplacian in the bulk */ #pragma omp parallel for private(i,j) for (i=1; i 0.0) t = (double)time*OMEGA - (double)j*OSCIL_LEFT_YSHIFT/(double)NY; else t = (double)time*OMEGA + (double)(NY-j)*OSCIL_LEFT_YSHIFT/(double)NY; if (t < 0.0) return(0.0); else return(AMPLITUDE*cos(t)*exp(-(double)t*DAMPING)); } case (OSC_SLOWING): { a = 0.0025; t = (double)time*OMEGA; phase = t - a*t*t; // if (time%1000 == 0) printf("time = %i, phase = %.4lg\n", time, phase); return(AMPLITUDE*cos(phase)*exp(-phase*DAMPING)); } case (OSC_WAVE_PACKET): { t = (double)time*OMEGA; // a = 10.0; // a = 0.02/OMEGA; a = sqrt(INITIAL_VARIANCE)/OMEGA; phase = AMPLITUDE*cos(t); envelope = exp(-(t-0.2)*(t-0.2)/(a*a))*sqrt(DPI/a); return(phase*envelope); } case (OSC_CHIRP): { // a = 0.25; t = (double)time*OMEGA; phase = t + ACHIRP*t*t; return(AMPLITUDE*sin(phase)*exp(-phase*DAMPING)); } case (OSC_BEAM): { xy_to_ij(0.0, OSCIL_YMAX, ij); jmax = ij[1]; jmin = NY - jmax; if (j < jmin) dist2 = (double)(jmin-1-j)*(YMAX-YMIN)/(double)NY; else if (j > jmax) dist2 = (double)(j - jmax)*(YMAX-YMIN)/(double)NY; else dist2 = 0.0; dist2 = 500.0*dist2*dist2; t = (double)time*OMEGA; amp = AMPLITUDE*cos(t)*exp(-(double)t*DAMPING); amp *= exp(-dist2/INITIAL_VARIANCE); return(amp); } case (OSC_BEAM_GAUSSIAN): { dist2 = (double)(j - NY/2)*(YMAX-YMIN)/(double)NY; dist2 = dist2*dist2; t = (double)time*OMEGA; amp = AMPLITUDE*cos(t)*exp(-(double)t*DAMPING); amp *= exp(-dist2/INITIAL_VARIANCE); return(amp); } case (OSC_BEAM_SINE): { xy_to_ij(0.0, OSCIL_YMAX, ij); jmax = ij[1]; jmin = NY - jmax; dist2 = (double)(j - NY/2)/(double)(jmax - NY/2); if (j > jmax) return(0.0); if (j < jmin) return(0.0); t = (double)time*OMEGA; amp = AMPLITUDE*cos(t)*exp(-(double)t*DAMPING); amp *= cos(PID*dist2); return(amp); } case (OSC_BEAM_TWOPERIODS): { xy_to_ij(0.0, OSCIL_YMAX, ij); jmax = ij[1]; jmin = NY - jmax; if (j < jmin) dist2 = (double)(jmin-1-j)*(YMAX-YMIN)/(double)NY; else if (j > jmax) dist2 = (double)(j - jmax)*(YMAX-YMIN)/(double)NY; else dist2 = 0.0; dist2 = 500.0*dist2*dist2; t = (double)time*OMEGA; amp = AMPLITUDE*(0.5*cos(t) + 0.5*cos(sqrt(7.0)*t))*exp(-(double)t*DAMPING); amp *= exp(-dist2/INITIAL_VARIANCE); return(amp); } } } void init_ior_2d(short int *xy_in[NX], double *tcc_table[NX], double *tgamma_table[NX], double ior_angle) /* compute variable index of refraction */ /* should be at some point merged with 3D version in suv_wave_3d.c */ { int i, j, k, n, inlens, ncircles; double courant2 = COURANT*COURANT, courantb2 = COURANTB*COURANTB, lambda1, mu1, courant_med, courant_med2; double a, u, v, u1, x, y, xy[2], norm2, speed, r2, c, salpha, h, ll, ca, sa, x1, y1, dx, dy, sum, sigma, x0, y0, rgb[3]; double xc[NGRIDX*NGRIDY], yc[NGRIDX*NGRIDY], height[NGRIDX*NGRIDY]; static double xc_stat[NGRIDX*NGRIDY], yc_stat[NGRIDX*NGRIDY], sigma_stat; static int first = 1; t_circle circles[NMAXCIRCLES]; rgb[0] = 1.0; rgb[1] = 1.0; rgb[2] = 1.0; if (VARIABLE_IOR) { switch (IOR) { case (IOR_MANDELBROT): { #pragma omp parallel for private(i,j) for (i=0; i= MANDELLEVEL) { // tc[i*NY+j] = COURANT; tcc_table[i][j] = courant2; tgamma_table[i][j] = GAMMA; } else { speed = (double)k/(double)MANDELLEVEL; if (speed < 1.0e-10) speed = 1.0e-10; else if (speed > 10.0) speed = 10.0; tcc_table[i][j] = courant2*speed; // tc[i*NY+j] = COURANT*sqrt(speed); tgamma_table[i][j] = GAMMA; } } } break; } case (IOR_EARTH): { for (i=0; i 1.0) c = 0.0; else if (r2 < 0.25*0.25) c = 0.8*COURANT; else if (r2 < 0.58*0.58) c = COURANT*(0.68 - 0.55*r2); else c = COURANT*(1.3 - 0.9*r2); // tc[i*NY+j] = c; tcc_table[i][j] = c; tgamma_table[i][j] = GAMMA; } } break; } case (IOR_EXPLO_LENSING): { salpha = DPI/(double)NPOLY; // lambda1 = LAMBDA; // mu1 = LAMBDA; lambda1 = 0.5*LAMBDA; mu1 = 0.5*LAMBDA; h = lambda1*tan(PI/(double)NPOLY); if (h < mu1) ll = sqrt(mu1*mu1 - h*h); else ll = 0.0; // #pragma omp parallel for private(i,j) for (i=0; i 1.0) height[n] = 1.0; if (height[n] < 0.0) height[n] = 0.0; } // #pragma omp parallel for private(i,j) for (i=0; i 1.0) height[n] = 1.0; if (height[n] < 0.0) height[n] = 0.0; } for (n = 0; n 1.0) height[n] = 1.0; if (height[n] < 0.0) height[n] = 0.0; } for (n = 0; n MU)) { tcc_table[i][j] = 0.0; tgamma_table[i][j] = GAMMA; } else { if (xy_in[i][j] == 1) { tcc_table[i][j] = courant2; tgamma_table[i][j] = GAMMA; } else if (xy_in[i][j] == 0) { tcc_table[i][j] = courantb2; tgamma_table[i][j] = GAMMAB; } else { tcc_table[i][j] = 0.0; tgamma_table[i][j] = 1.0; } } } } break; } case (IOR_LENS_OBSTACLE): { printf("Initializing IOR_LENS_OBSTACLE\n"); for (i=0; i MU)) { tcc_table[i][j] = 0.0; tgamma_table[i][j] = GAMMA; } else { if (xy_in[i][j] != 0) { tcc_table[i][j] = courant2; tgamma_table[i][j] = GAMMA; } else { tcc_table[i][j] = courantb2; tgamma_table[i][j] = GAMMAB; } } } } break; } case (IOR_LENS_CONVEX_CONCAVE): { a = LAMBDA + 0.1 - sqrt(LAMBDA*LAMBDA - MU*MU) + 0.6; for (i=0; i MU)) { tcc_table[i][j] = 0.0; tgamma_table[i][j] = GAMMA; } else { if (xy_in[i][j] != 0) { tcc_table[i][j] = courant2; tgamma_table[i][j] = GAMMA; } else { tcc_table[i][j] = courantb2; tgamma_table[i][j] = GAMMAB; } } } } break; } case (IOR_TREE): { for (i=0; i LAMBDA) { tcc_table[i][j] = 0.0; tgamma_table[i][j] = 1.0; } else { tcc_table[i][j] = courantb2; tgamma_table[i][j] = GAMMAB; } } } break; } case (IOR_WAVE_GUIDES_COUPLED_B): { for (i=0; i LAMBDA)&&(vabs(y) > 0.7)) { tcc_table[i][j] = 0.0; tgamma_table[i][j] = 1.0; } else { tcc_table[i][j] = courantb2; tgamma_table[i][j] = GAMMAB; } } } break; } case (IOR_WAVE_GUIDE_COATED): { courant_med = 0.8*COURANT + 0.2*COURANTB; courant_med2 = courant_med*courant_med; for (i=0; i rmin2)) { r2 = (double)(irad2); wave_height1 = wave_height*exp(-r2/variance); phi[packet[i].ix + j][packet[i].iy + k] = wave_height1; phi[packet[i].ix - j][packet[i].iy + k] = wave_height1; phi[packet[i].ix + j][packet[i].iy - k] = wave_height1; phi[packet[i].ix - j][packet[i].iy - k] = wave_height1; } else if (irad2 <= rmin2) { phi[packet[i].ix + j][packet[i].iy + k] = 0.0; phi[packet[i].ix - j][packet[i].iy + k] = 0.0; phi[packet[i].ix + j][packet[i].iy - k] = 0.0; phi[packet[i].ix - j][packet[i].iy - k] = 0.0; } } // printf("Adding wave packet %i with value %.3lg\n", i, wave_height); // if (add==0) { // t -= 1.0/(double)NVID; t -= 0.1/(double)NVID; wave_height = wave_packet_height(t, packet[i], type_envelope); for (j=0; j rmin2)) { r2 = (double)(irad2); wave_height1 = wave_height*exp(-r2/variance); psi[packet[i].ix + j][packet[i].iy + k] = wave_height1; psi[packet[i].ix - j][packet[i].iy + k] = wave_height1; psi[packet[i].ix + j][packet[i].iy - k] = wave_height1; psi[packet[i].ix - j][packet[i].iy - k] = wave_height1; } else if (irad2 <= rmin2) { psi[packet[i].ix + j][packet[i].iy + k] = 0.0; psi[packet[i].ix - j][packet[i].iy + k] = 0.0; psi[packet[i].ix + j][packet[i].iy - k] = 0.0; psi[packet[i].ix - j][packet[i].iy - k] = 0.0; } } } } } void add_circular_wave_loc(double factor, double x, double y, double *phi[NX], double *psi[NX], short int * xy_in[NX], int xmin, int xmax, int ymin, int ymax) /* add drop at (x,y) to the field with given prefactor */ { int i, j; double xy[2], dist2; // for (i=0; i (double)packet[i].var_envelope) amp *= exp(-0.1*(t - (double)packet[i].var_envelope)); // if (t < (double)packet[i].var_envelope + 100.0) // add_circular_wave_loc(amp, packet[i].xc, packet[i].yc, phi, psi, xy_in, xmin, xmax, ymin, ymax); } } void add_wave_packets(double *phi[NX], double *psi[NX], short int * xy_in[NX], t_wave_packet *packet, int time, int radius, int local, int add_period, int type_envelope) /* add some wave packet sources */ { if (local) add_wave_packets_locally(phi, psi, packet, time, radius, add_period, type_envelope); else add_wave_packets_globally(phi, psi, xy_in, packet, time, add_period, type_envelope); } int old_source_schedule(int i) { int mod; if (i < 200) return(0); mod = i%(10*OSCILLATING_SOURCE_PERIOD); if ((mod < 2*OSCILLATING_SOURCE_PERIOD)&&(rand()%3 < 2)) return(1); return(0); } void old_init_input_signal() { int i; for (i=0; i