From d79439c0087e9e575ec5fc1ffbd2c467f299e946 Mon Sep 17 00:00:00 2001 From: nilsberglund-orleans <83530463+nilsberglund-orleans@users.noreply.github.com> Date: Mon, 27 Sep 2021 18:15:19 +0200 Subject: [PATCH] Add files via upload --- particle_pinball.c | 886 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 886 insertions(+) create mode 100644 particle_pinball.c diff --git a/particle_pinball.c b/particle_pinball.c new file mode 100644 index 0000000..602ae1f --- /dev/null +++ b/particle_pinball.c @@ -0,0 +1,886 @@ +/*********************************************************************************/ +/* */ +/* Animation of particles in billiard */ +/* */ +/* N. Berglund, december 2012, april 2021 */ +/* UPDATE 14 April 21 : graphics files go to subfolder, */ +/* Switch MOVIE to decide whether to create a movie */ +/* UPDATE 3 May 21 : new domains */ +/* */ +/* Feel free to reuse, but if doing so it would be nice to drop a */ +/* line to nils.berglund@univ-orleans.fr - Thanks! */ +/* */ +/* compile with */ +/* gcc -o particle_billiard particle_billiard.c */ +/* -O3 -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lXmu -lglut */ +/* */ +/* To make a video, set MOVIE to 1 and create subfolder tif_part */ +/* It may be possible to increase parameter PAUSE */ +/* */ +/* create movie using */ +/* ffmpeg -i part.%05d.tif -vcodec libx264 part.mp4 */ +/* */ +/*********************************************************************************/ + +#include +#include +#include +#include +#include +#include +#include /* Sam Leffler's libtiff library. */ + +#define MOVIE 0 /* set to 1 to generate movie */ + +#define WINWIDTH 1280 /* window width */ +#define WINHEIGHT 720 /* window height */ + +#define XMIN -4.0 +#define XMAX 4.0 /* x interval */ +#define YMIN -1.25 +#define YMAX 3.25 /* y interval for 9/16 aspect ratio */ +// #define XMIN -2.0 +// #define XMAX 2.0 /* x interval */ +// #define YMIN -1.125 +// #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ + +#define BOXYMIN -1.0 +#define BOXYMAX 1.0 /* y dimensions of box (for circles in rectangle) */ + +#define SCALING_FACTOR 1.0 /* scaling factor of drawing, needed for flower billiards, otherwise set to 1.0 */ + +/* Choice of the billiard table, see global_particles.c */ + +#define B_DOMAIN 21 /* choice of domain shape */ + +#define CIRCLE_PATTERN 6 /* pattern of circles */ + +#define ABSORBING_CIRCLES 0 /* set to 1 for circular scatterers to be absorbing */ + +#define NMAXCIRCLES 5000 /* total number of circles (must be at least NCX*NCY for square grid) */ +#define NCX 40 /* number of circles in x direction */ +#define NCY 10 /* number of circles in y direction */ +#define NPOISSON 500 /* number of points for Poisson C_RAND_POISSON arrangement */ +#define NGOLDENSPIRAL 2000 /* max number of points for C_GOLDEN_SPIRAL arrandement */ + +#define LAMBDA 3.8 /* parameter controlling shape of domain */ +#define MU 0.05 /* second parameter controlling shape of billiard */ +#define FOCI 1 /* set to 1 to draw focal points of ellipse */ +#define NPOLY 4 /* number of sides of polygon */ +#define APOLY 0.5 /* angle by which to turn polygon, in units of Pi/2 */ +#define DRAW_BILLIARD 1 /* set to 1 to draw billiard */ +#define DRAW_CONSTRUCTION_LINES 0 /* set to 1 to draw additional construction lines for billiard */ +#define PERIODIC_BC 0 /* set to 1 to enforce periodic boundary conditions when drawing particles */ + +#define RESAMPLE 0 /* set to 1 if particles should be added when dispersion too large */ +#define DEBUG 0 /* draw trajectories, for debugging purposes */ + +/* Simulation parameters */ + +#define NPART 1 /* number of particles */ +#define NPARTMAX 100000 /* maximal number of particles after resampling */ +#define LMAX 0.01 /* minimal segment length triggering resampling */ +#define DMIN 0.02 /* minimal distance to boundary for triggering resampling */ +#define CYCLE 1 /* set to 1 for closed curve (start in all directions) */ +#define SHOWTRAILS 0 /* set to 1 to keep trails of the particles */ +#define TEST_ACTIVE 0 /* set to 1 to test whether particle is in billiard */ + +#define NSTEPS 7700 /* number of frames of movie */ +#define TIME 4000 /* time between movie frames, for fluidity of real-time simulation */ +#define DPHI 0.00007 /* integration step */ +#define NVID 150 /* number of iterations between images displayed on screen */ + +/* Decreasing TIME accelerates the animation and the movie */ +/* For constant speed of movie, TIME*DPHI should be kept constant */ +/* However, increasing DPHI too much deterioriates quality of simulation */ +/* NVID tells how often a picture is drawn in the animation, increase it for faster anim */ +/* For a good quality movie, take for instance TIME = 400, DPHI = 0.00005, NVID = 100 */ + +/* Colors and other graphical parameters */ + +#define NCOLORS 32 /* number of colors */ +#define COLORSHIFT 0 /* hue of initial color */ +#define RAINBOW_COLOR 1 /* set to 1 to use different colors for all particles */ +#define SINGLE_COLOR 1 /* set to 1 to make all particles a single color */ +#define FLOWER_COLOR 0 /* set to 1 to adapt initial colors to flower billiard (tracks vs core) */ +#define NSEG 100 /* number of segments of boundary */ +#define LENGTH 0.1 /* length of velocity vectors */ +#define BILLIARD_WIDTH 2 /* width of billiard */ +#define PARTICLE_WIDTH 4 /* width of particles */ +#define FRONT_WIDTH 3 /* width of wave front */ +#define COLOR_TRAJECTORY 8 /* hue for single color */ + +#define BLACK 1 /* set to 1 for black background */ +#define COLOR_OUTSIDE 0 /* set to 1 for colored outside */ +#define OUTER_COLOR 270.0 /* color outside billiard */ +#define PAINT_INT 0 /* set to 1 to paint interior in other color (for polygon/Reuleaux) */ +#define PAINT_EXT 0 /* set to 1 to paint exterior in other color */ +#define ERASE_OUTSIDE 1 /* set to 1 to erase outside of rectangular billiard (beta) */ + + +#define PAUSE 1000 /* number of frames after which to pause */ +#define PSLEEP 5 /* sleep time during pause */ +#define SLEEP1 1 /* initial sleeping time */ +#define SLEEP2 1000 /* final sleeping time */ +#define END_FRAMES 100 /* number of still frames at end of movie */ + +#define NPATHBINS 200 /* number of bins for path length histogramm */ +#define PATHLMAX 2.5 /* max free path on graph */ + +#include "global_particles.c" +#include "sub_part_billiard.c" +#include "sub_part_pinball.c" + +int ncol = 0, nobst = 0, nmaxpeg = 0; +int npath[NPATHBINS]; + +/*********************/ +/* animation part */ +/*********************/ + +void init_boundary_config(double smin, double smax, double anglemin, double anglemax, double *configs[NPARTMAX]) +/* initialize configuration: drop on the boundary, beta version */ +/* WORKS FOR ELLIPSE, HAS TO BE ADAPTED TO GENERAL BILLIARD */ +{ + int i; + double ds, da, s, angle, theta, alpha, pos[2]; + + if (anglemin <= 0.0) anglemin = PI/((double)NPART); + if (anglemax >= PI) anglemax = PI*(1.0 - 1.0/((double)NPART)); + ds = (smax - smin)/((double)NPART); + da = (anglemax - anglemin)/((double)NPART); + for (i=0; i 1) dalpha = (angle2 - angle1)/((double)(NPART-1)); + else dalpha = 0.0; + for (i=0; i= NCOLORS) color[i] -= NCOLORS; +// } +// } + + configs[i][2] += DPHI; + + cosphi = (configs[i][6] - configs[i][4])/configs[i][3]; + sinphi = (configs[i][7] - configs[i][5])/configs[i][3]; + len = configs[i][2] + LENGTH; + if (len > configs[i][3]) len = configs[i][3]; + + x0 = configs[i][4]; + y0 = configs[i][5]; + x1 = configs[i][4] + configs[i][2]*cosphi; + y1 = configs[i][5] + configs[i][2]*sinphi; + x2 = configs[i][4] + len*cosphi; + y2 = configs[i][5] + len*sinphi; + + /* test whether particle does not escape billiard */ + if ((TEST_ACTIVE)&&(active[i])) active[i] = xy_in_billiard(x1, y1); + + if (active[i]) + { + rgb_color_scheme(color[i], rgb); + glColor3f(rgb[0], rgb[1], rgb[2]); + + glBegin(GL_LINE_STRIP); + glVertex2d(SCALING_FACTOR*x0, SCALING_FACTOR*y0); + glVertex2d(SCALING_FACTOR*x2, SCALING_FACTOR*y2); + glEnd (); + } + +// if (configs[i][2] > configs[i][3] - DPHI) +// { +// glBegin(GL_LINE_STRIP); +// glVertex2d(SCALING_FACTOR*x0, SCALING_FACTOR*y0); +// glVertex2d(SCALING_FACTOR*configs[i][6], SCALING_FACTOR*configs[i][7]); +// glEnd (); +// } + } + if (DRAW_BILLIARD) draw_billiard(); +} + + + +void draw_config(int color[NPARTMAX], double *configs[NPARTMAX], int active[NPARTMAX]) +/* draw the particles */ +{ + int i; + double x0, y0, x1, y1, x2, y2, cosphi, sinphi, rgb[3]; + + glutSwapBuffers(); + if (!SHOWTRAILS) blank(); + if (PAINT_INT) paint_billiard_interior(); + + glLineWidth(PARTICLE_WIDTH); + + glEnable(GL_LINE_SMOOTH); + + for (i=0; i= NCOLORS) color[i] -= NCOLORS; +// } +// } + + configs[i][2] += DPHI; + + cosphi = (configs[i][6] - configs[i][4])/configs[i][3]; + sinphi = (configs[i][7] - configs[i][5])/configs[i][3]; + x1 = configs[i][4] + configs[i][2]*cosphi; + y1 = configs[i][5] + configs[i][2]*sinphi; + x2 = configs[i][4] + (configs[i][2] + LENGTH)*cosphi; + y2 = configs[i][5] + (configs[i][2] + LENGTH)*sinphi; + + /* test whether particle does not escape billiard */ + if ((TEST_ACTIVE)&&(active[i])) active[i] = xy_in_billiard(x1, y1); + + if (active[i]) + { + rgb_color_scheme(color[i], rgb); + glColor3f(rgb[0], rgb[1], rgb[2]); + + glBegin(GL_LINE_STRIP); + glVertex2d(SCALING_FACTOR*x1, SCALING_FACTOR*y1); + glVertex2d(SCALING_FACTOR*x2, SCALING_FACTOR*y2); + glEnd (); + + /* taking care of boundary conditions - only needed for periodic boundary conditions */ + if (PERIODIC_BC) + { + if (SCALING_FACTOR*x2 > XMAX) + { + glBegin(GL_LINE_STRIP); + glVertex2d(SCALING_FACTOR*(x1+XMIN-XMAX), SCALING_FACTOR*y1); + glVertex2d(SCALING_FACTOR*(x2+XMIN-XMAX), SCALING_FACTOR*y2); + glEnd (); + } + + if (SCALING_FACTOR*x2 < XMIN) + { + glBegin(GL_LINE_STRIP); + glVertex2d(SCALING_FACTOR*(x1-XMIN+XMAX), SCALING_FACTOR*y1); + glVertex2d(SCALING_FACTOR*(x2-XMIN+XMAX), SCALING_FACTOR*y2); + glEnd (); + } + + if (SCALING_FACTOR*y2 > YMAX) + { + glBegin(GL_LINE_STRIP); + glVertex2d(SCALING_FACTOR*x1, SCALING_FACTOR*(y1+YMIN-YMAX)); + glVertex2d(SCALING_FACTOR*x2, SCALING_FACTOR*(y2+YMIN-YMAX)); + glEnd (); + } + + if (SCALING_FACTOR*y2 < YMIN) + { + glBegin(GL_LINE_STRIP); + glVertex2d(SCALING_FACTOR*x1, SCALING_FACTOR*(y1+YMAX-YMIN)); + glVertex2d(SCALING_FACTOR*x2, SCALING_FACTOR*(y2+YMAX-YMIN)); + glEnd (); + } + } + } + + /* draw trajectories, for debugging purpose */ + if (DEBUG) + { + glLineWidth(1.0); + glBegin(GL_LINES); + glVertex2d(SCALING_FACTOR*configs[i][4], SCALING_FACTOR*configs[i][5]); + glVertex2d(SCALING_FACTOR*configs[i][6], SCALING_FACTOR*configs[i][7]); + glEnd (); + glLineWidth(3.0); + } + +// if (configs[i][2] > configs[i][3] - DPHI) configs[i][2] -= configs[i][3]; + } + if (DRAW_BILLIARD) draw_billiard(); +} + + +void graph_movie(int time, int color[NPARTMAX], double *configs[NPARTMAX], int active[NPARTMAX]) +/* compute next movie frame */ +{ + int i, j, k, c; + double rgb[3]; + + for (j=0; j= 0)&&(circlecolor[c] == 0)) nobst++; + circlecolor[c]++; + newcircle[c] = 10; + + /* update free path statistics */ + k = (int)((double)NPATHBINS*configs[i][3]/PATHLMAX); + if (k < NPATHBINS) npath[k]++; + +// if (circlecolor[c] > nmaxpeg) nmaxpeg = circlecolor[c]; + + +// if (circlecolor[c] > NCOLORS) circlecolor[c] = 1; +// if (c>=0) color[i]++; + + + if (!RAINBOW_COLOR) + { + color[i]++; + if (color[i] >= NCOLORS) color[i] -= NCOLORS; + } + } + + configs[i][2] += DPHI; + + if (configs[i][2] > configs[i][3] - DPHI) + { + configs[i][2] -= configs[i][3]; + } + } + } + +// draw_config(color, configs); +} + + +void print_particle_numbers(double *configs[NPARTMAX]) +{ + char message[50], message1[50]; + double cosphi, x1; + static double rgb[3], xleft, xright; + static short int first = 1; + int i, nleft = 0, nmid = 0, nright = 0; + + rgb[0] = 0.0; rgb[1] = 0.0; rgb[2] = 0.0; + + if (first) /* compute box limits */ + { + /* find leftmost and rightmost circle */ + for (i=0; i xright)) xright = circlex[i] + circlerad[i]; + + first = 0; + + printf("xleft = %.3lg, xright = %.3lg", xleft, xright); + } + + for (i=0; i xright) nright++; + else if (x1 < xleft) nleft++; + else nmid++; + +// if (i == nparticles-1) printf("x1 = %.3lg, nleft = %i, nright = %i\n", x1, nleft, nright); + } + + erase_area(XMIN + 0.31, YMIN + 0.07, 0.24, 0.05, rgb); + sprintf(message, "%4d particles", nleft); + glColor3f(1.0, 1.0, 1.0); + write_text_fixedwidth(XMIN + 0.1, YMIN + 0.04, message); + + erase_area(0.0, YMIN + 0.07, 0.24, 0.05, rgb); + sprintf(message, "%4d particles", nmid); + glColor3f(1.0, 1.0, 1.0); + write_text_fixedwidth(-0.21, YMIN + 0.04, message); + + erase_area(XMAX - 0.29, YMIN + 0.07, 0.24, 0.05, rgb); + sprintf(message, "%4d particles", nright); + glColor3f(1.0, 1.0, 1.0); + write_text_fixedwidth(XMAX - 0.5, YMIN + 0.04, message); +} + +void draw_statistics() +{ + int i, n, colmax = 35, pegcollisions[70], nypegs = 70, meanpegs = 0, total_coll = 0, ymax = 0, meanbins = 0, total_bin = 0; + double x, y, yscale = 120.0, y0, dx, rgb[3], xshift; + char message[50]; + + glLineWidth(1); + + y0 = 0.5*(YMAX + YMIN) + 0.2; + dx = (XMAX-0.6)/(double)colmax; + xshift = XMIN + 0.3; + rgb[0] = 0.0; rgb[1] = 0.0; rgb[2] = 1.0; + + /* histogramm of number of collisions per peg */ + for (i=0; i ymax) ymax = npath[i]; + total_bin += npath[i]; + meanbins += i*npath[i]; + } + + yscale = 0.9*(YMAX-y0)*(double)ymax; + dx = (XMAX-0.6)/(double)NPATHBINS; + xshift = 0.3; + rgb[0] = 1.0; rgb[1] = 0.0; rgb[2] = 0.0; + + for (i=1; i 0.99) circleactive[i] = 0; + if (vabs(circlex[i]) + circlerad[i] > 0.99*LAMBDA) circleactive[i] = 0; + } +// if (vabs(circley[i]) > 1.0) circleactive[i] = 0; + + /* initialize system by putting particles in a given point with a range of velocities */ + r = cos(PI/(double)NPOLY)/cos(DPI/(double)NPOLY); + +// init_line_config(-1.25, -0.5, -1.25, 0.5, 0.0, configs); + init_drop_config(0.5, 0.1, -0.5*PID, 0.5*PID, configs); +// init_drop_config(0.0, 0.0, -0.5*PID, 0.5*PID, configs); +// init_drop_config(-1.4, 0.0, -0.5*PID, 0.5*PID, configs); +// init_drop_config(0.5, 0.5, -1.0, 1.0, configs); +// init_sym_drop_config(-1.0, 0.5, -PID, PID, configs); +// init_drop_config(-0.999, 0.0, -alpha, alpha, configs); + +// other possible initial conditions : +// init_line_config(0.0, -0.5, 0.0, 0.5, 0.0, configs); +// init_line_config(-1.25, -0.5, -1.25, 0.5, 0.0*PID, configs); +// init_line_config(-1.0, -0.3, -1.0, 0.3, 0.0, configs); +// init_line_config(-0.7, -0.45, -0.7, 0.45, 0.0, configs); +// init_line_config(-1.5, 0.1, -0.1, 1.0, -0.5*PID, configs); + +// if (!SHOWTRAILS) blank(); + blank(); + glColor3f(0.0, 0.0, 0.0); + if (DRAW_BILLIARD) draw_billiard(); + if (ERASE_OUTSIDE) erase_rectangle_outside(270.0, 0.1, 0.15); +// print_particle_numbers(configs); + + glutSwapBuffers(); + +// if (MOVIE) +// { +// for (i=0; i<20; i++) save_frame(); +// s = system("mv part*.tif tif_part/"); +// } + + + for (i=0; i nmaxpeg) nmaxpeg = circlecolor[j]; + + sprintf(message, "max hits per peg: %d", nmaxpeg); + glColor3f(1.0, 1.0, 1.0); + write_text(-0.6, YMIN + 0.08, message); +// write_text(-0.3, YMIN + 0.04, message); + + draw_statistics(); + + for (j=0; j