From 416366d8daa29458846749459c761eb65a840099 Mon Sep 17 00:00:00 2001 From: nilsberglund-orleans <83530463+nilsberglund-orleans@users.noreply.github.com> Date: Tue, 28 Dec 2021 15:42:56 +0100 Subject: [PATCH] Add files via upload --- drop_billiard.c | 36 +- global_ljones.c | 122 ++++++ global_particles.c | 5 +- global_pdes.c | 13 +- heat.c | 3 +- lennardjones.c | 897 ++++++++++++++++++++++++++++++++++++++++++++ mangrove.c | 1 + particle_billiard.c | 147 +++++--- schrodinger.c | 1 + sub_lj.c | 666 ++++++++++++++++++++++++++++++++ sub_part_billiard.c | 237 +++++++++--- sub_wave.c | 323 +++++++++++++++- sub_wave_comp.c | 150 ++++++-- wave_billiard.c | 90 ++--- wave_common.c | 14 + wave_comparison.c | 236 +++--------- wave_energy.c | 9 + 17 files changed, 2571 insertions(+), 379 deletions(-) create mode 100644 global_ljones.c create mode 100644 lennardjones.c create mode 100644 sub_lj.c diff --git a/drop_billiard.c b/drop_billiard.c index ae9f302..cf0a7cd 100644 --- a/drop_billiard.c +++ b/drop_billiard.c @@ -16,10 +16,10 @@ /* It may be possible to increase parameter PAUSE */ /* */ /* create movie using */ -/* ffmpeg -i tif_drop/part.%05d.tif -vcodec libx264 drop.mp4 */ +/* ffmpeg -i part.%05d.tif -vcodec libx264 drop.mp4 */ /* */ /*********************************************************************************/ - + #include #include #include @@ -28,16 +28,10 @@ #include #include /* Sam Leffler's libtiff library. */ -#define MOVIE 1 /* set to 1 to generate movie */ +#define MOVIE 0 /* set to 1 to generate movie */ -// #define WINWIDTH 1280 /* window width */ -// #define WINHEIGHT 720 /* window height */ - -#define WINWIDTH 1920 /* window width */ -#define WINHEIGHT 1080 /* window height */ - -// #define WINWIDTH 2560 /* window width */ -// #define WINHEIGHT 1440 /* window height */ +#define WINWIDTH 1280 /* window width */ +#define WINHEIGHT 720 /* window height */ #define XMIN -2.0 #define XMAX 2.0 /* x interval */ @@ -65,7 +59,7 @@ #define NGOLDENSPIRAL 2000 /* max number of points for C_GOLDEN_SPIRAL arrandement */ #define SDEPTH 1 /* Sierpinski gastket depth */ -#define LAMBDA 0.3 /* parameter controlling the dimensions of domain - with 0.5 @1440p it cut out top-bottom borders. 0.4 is @limit */ +#define LAMBDA 0.3 /* parameter controlling the dimensions of domain */ // #define LAMBDA 1.124950941 /* sin(36°)/sin(31.5°) for 5-star shape with 45° angles */ // #define LAMBDA 1.445124904 /* sin(36°)/sin(24°) for 5-star shape with 60° angles */ // #define LAMBDA 3.75738973 /* sin(36°)/sin(9°) for 5-star shape with 90° angles */ @@ -85,7 +79,7 @@ #define NPARTMAX 100000 /* maximal number of particles after resampling */ #define NSTEPS 5000 /* number of frames of movie */ -#define TIME 150 /* time between movie frames, for fluidity of real-time simulation */ +#define TIME 150 /* time between movie frames, for fluidity of real-time simulation */ #define DPHI 0.0001 /* integration step */ #define NVID 75 /* number of iterations between images displayed on screen */ @@ -107,29 +101,31 @@ #define COLOR_PALETTE 1 /* Color palette, see list in global_pdes.c */ -#define NCOLORS 14 /* number of colors */ +#define NCOLORS 14 /* number of colors */ #define COLORSHIFT 3 /* hue of initial color */ #define RAINBOW_COLOR 0 /* set to 1 to use different colors for all particles */ -#define NSEG 150 /* number of segments of boundary */ +#define NSEG 100 /* number of segments of boundary */ #define BILLIARD_WIDTH 4 /* width of billiard */ #define FRONT_WIDTH 4 /* width of wave front */ -#define BLACK 1 /* set to 1 for black background */ -#define COLOR_OUTSIDE 1 /* set to 1 for colored outside */ +#define BLACK 0 /* set to 1 for black background */ +#define COLOR_OUTSIDE 0 /* set to 1 for colored outside */ #define OUTER_COLOR 300.0 /* color outside billiard */ #define PAINT_INT 0 /* set to 1 to paint interior in other color (for polygon) */ #define PAINT_EXT 1 /* set to 1 to paint exterior of billiard */ + #define PAUSE 1000 /* number of frames after which to pause */ #define PSLEEP 1 /* sleep time during pause */ #define SLEEP1 1 /* initial sleeping time */ #define SLEEP2 100 /* final sleeping time */ #define END_FRAMES 0 /* number of frames at end of movie */ -#define PI 3.141592654 +#define PI 3.141592654 #define DPI 6.283185307 #define PID 1.570796327 + #include "global_particles.c" #include "sub_part_billiard.c" @@ -200,6 +196,8 @@ void init_partial_drop_config(int i1, int i2, double x0, double y0, double angle } } + + int resample(int color[NPARTMAX], double *configs[NPARTMAX]) /* add particles where the front is stretched too thin */ { @@ -412,7 +410,7 @@ void graph_movie(int time, int color[NPARTMAX], double *configs[NPARTMAX]) { // print_config(configs[i]); - if (configs[i][2]<0.0) + if (configs[i][2]<0.0) { vbilliard(configs[i]); if (!RAINBOW_COLOR) diff --git a/global_ljones.c b/global_ljones.c new file mode 100644 index 0000000..67d3dbf --- /dev/null +++ b/global_ljones.c @@ -0,0 +1,122 @@ +/* Global variables and parameters for lennardjones */ + +/* Basic math */ + +#define PI 3.141592654 +#define DPI 6.283185307 +#define PID 1.570796327 + +/* shape of domain */ + +#define D_CIRCLES 20 /* several circles */ +#define D_CIRCLES_IN_RECT 201 /* several circles in a rectangle */ + +#define NMAXCIRCLES 20000 /* total number of circles/polygons (must be at least NCX*NCY for square grid) */ +#define MAXNEIGH 20 /* max number of neighbours kept in memory */ + +#define C_SQUARE 0 /* square grid of circles */ +#define C_HEX 1 /* hexagonal/triangular grid of circles */ +#define C_RAND_DISPLACED 2 /* randomly displaced square grid */ +#define C_RAND_PERCOL 3 /* random percolation arrangement */ +#define C_RAND_POISSON 4 /* random Poisson point process */ +#define C_CLOAK 5 /* invisibility cloak */ +#define C_CLOAK_A 6 /* first optimized invisibility cloak */ +#define C_LASER 7 /* laser fight in a room of mirrors */ +#define C_POISSON_DISC 8 /* Poisson disc sampling */ + +#define C_GOLDEN_MEAN 10 /* pattern based on vertical shifts by golden mean */ +#define C_GOLDEN_SPIRAL 11 /* spiral pattern based on golden mean */ +#define C_SQUARE_HEX 12 /* alternating between square and hexagonal/triangular */ + +#define C_ONE 97 /* one single circle, as for Sinai */ +#define C_TWO 98 /* two concentric circles of different type */ +#define C_NOTHING 99 /* no circle at all, for comparisons */ + +/* particle interaction */ + +#define I_COULOMB 0 /* Coulomb force */ +#define I_LENNARD_JONES 1 /* Lennard-Jones force */ +#define I_LJ_DIRECTIONAL 2 /* Lennard-Jones with direction dependence */ +#define I_LJ_PENTA 3 /* Lennard-Jones with pentagonal symmetry */ +#define I_GOLDENRATIO 4 /* Lennard-Jones type with equilibria related by golden ratio */ + +/* Boundary conditions */ + +#define BC_SCREEN 0 /* harmonic boundary conditions outside screen area */ +#define BC_RECTANGLE 1 /* harmonic boundary conditions on a resizeable rectangle */ +#define BC_CIRCLE 2 /* harmonic boundary conditions outside a moving circle */ + +/* Plot types */ + +#define P_KINETIC 0 /* colors represent kinetic energy of particles */ +#define P_NEIGHBOURS 1 /* colors represent number of neighbours */ +#define P_HEALTH 2 /* colors represent health (for SIR model) */ +#define P_BONDS 3 /* draw lattice based on neighbours */ + + +/* Color schemes */ + +#define C_LUM 0 /* color scheme modifies luminosity (with slow drift of hue) */ +#define C_HUE 1 /* color scheme modifies hue */ +#define C_PHASE 2 /* color scheme shows phase */ +#define C_ONEDIM 3 /* use preset 1d color scheme (for Turbo, Viridis, Magma, Inferno, Plasma) */ + +/* Color palettes */ + +#define COL_JET 0 /* JET color palette */ +#define COL_HSLUV 1 /* HSLUV color palette (perceptually uniform) */ + +#define COL_TURBO 10 /* TURBO color palette (by Anton Mikhailov) */ +#define COL_VIRIDIS 11 /* Viridis color palette */ +#define COL_MAGMA 12 /* Magma color palette */ +#define COL_INFERNO 13 /* Inferno color palette */ +#define COL_PLASMA 14 /* Plasma color palette */ +#define COL_CIVIDIS 15 /* Cividis color palette */ +#define COL_PARULA 16 /* Parula color palette */ + + +typedef struct +{ + double xc, yc, radius; /* center and radius of circle */ + short int active; /* circle is active */ + double energy; /* dissipated energy */ + double vx; /* x velocity of particle */ + double vy; /* y velocity of particle */ + double mass_inv; /* inverse of particle mass */ + double fx; /* x component of force on particle */ + double fy; /* y component of force on particle */ + int hashx; /* hash grid positions of particles */ + int hashy; /* hash grid positions of particles */ + int neighb; /* number of neighbours */ + int neighbours[MAXNEIGH]; /* coordinates of neighbours */ + double nghangle[MAXNEIGH]; /* angles of neighbours */ +} t_particle; + +typedef struct +{ + int number; /* total number of particles in cell */ + int particles[HASHMAX]; /* numbers of particles in cell */ +} t_hashgrid; + +typedef struct +{ + double xc, yc, radius; /* center and radius of circle */ + short int active; /* circle is active */ + double energy; /* dissipated energy */ + double vx; /* x velocity of particle */ + double vy; /* y velocity of particle */ + double mass_inv; /* inverse of particle mass */ + double fx; /* x component of force on particle */ + double fy; /* y component of force on particle */ + int hashx; /* hash grid positions of particles */ + int hashy; /* hash grid positions of particles */ + int neighb; /* number of neighbours */ + int health; /* 0 = sane, 1 = infected, 2 = recovered */ + double infected_time; /* time since infected */ + int protected; /* 0 = not protected, 1 = protected */ +} t_person; + + + +int ncircles, counter = 0; + diff --git a/global_particles.c b/global_particles.c index add4c29..dfba2a8 100644 --- a/global_particles.c +++ b/global_particles.c @@ -72,8 +72,11 @@ double x_shooter = -0.2, y_shooter = -0.6, x_target = 0.4, y_target = 0.7; #define P_RECTANGLE 0 /* rectangle (for test purposes) */ #define P_TOKARSKY 1 /* Tokarsky unilluminable room */ #define P_POLYRING 2 /* polygonal ring */ -#define P_SIERPINSKI 3 /* polygonal ring */ +#define P_SIERPINSKI 3 /* sierpinski carpet */ #define P_VONKOCH 4 /* von Koch curve */ +#define P_POLYGON 5 /* regular polygon, alternative for D_POLYGON */ +#define P_TOKA_PRIME 6 /* Tokarsky room made of 86 triangles */ +#define P_TREE 7 /* pine tree */ /* Color palettes */ diff --git a/global_pdes.c b/global_pdes.c index becb62e..8025bea 100644 --- a/global_pdes.c +++ b/global_pdes.c @@ -43,15 +43,16 @@ #define D_PENROSE 33 /* Penrose unilluminable room */ #define D_HYPERBOLA 34 /* one branch of hyperbola */ #define D_TOKARSKY 35 /* Tokarsky unilluminable room */ - +#define D_TOKA_PRIME 36 /* Tokarsky room made of 86 triangles */ #define D_ISOSPECTRAL 37 /* isospectral billiards */ #define D_HOMOPHONIC 38 /* homophonic billiards */ #define D_POLYGONS 40 /* several polygons */ #define D_VONKOCH 41 /* von Koch snowflake fractal */ +#define D_STAR 42 /* star shape */ #define NMAXCIRCLES 10000 /* total number of circles/polygons (must be at least NCX*NCY for square grid) */ -#define NMAXPOLY 10000 /* maximal number of vertices of polygonal lines (for von Koch et al) */ +#define NMAXPOLY 50000 /* maximal number of vertices of polygonal lines (for von Koch et al) */ // #define NMAXCIRCLES 10000 /* total number of circles/polygons (must be at least NCX*NCY for square grid) */ #define C_SQUARE 0 /* square grid of circles */ @@ -68,6 +69,8 @@ #define C_GOLDEN_SPIRAL 11 /* spiral pattern based on golden mean */ #define C_SQUARE_HEX 12 /* alternating between square and hexagonal/triangular */ +#define C_RINGS 20 /* obstacles arranged in concentruc rings */ + #define C_ONE 97 /* one single circle, as for Sinai */ #define C_TWO 98 /* two concentric circles of different type */ #define C_NOTHING 99 /* no circle at all, for comparisons */ @@ -103,6 +106,8 @@ #define P_ENERGY 1 /* plot energy of wave */ #define P_MIXED 2 /* plot amplitude in upper half, energy in lower half */ #define P_MEAN_ENERGY 3 /* energy averaged over time */ +#define P_LOG_ENERGY 4 /* log of energy averaged over time */ +#define P_LOG_MEAN_ENERGY 5 /* log of energy averaged over time */ /* For Schrodinger equation */ #define P_MODULE 10 /* plot module of wave function squared */ @@ -134,14 +139,14 @@ typedef struct { double xc, yc, radius; /* center and radius of circle */ - short int active; /* circle is active */ + short int active, top; /* circle is active, circle is in top half */ } t_circle; typedef struct { double xc, yc, radius, angle; /* center, radius and angle of polygon */ int nsides; /* number of sides of polygon */ - short int active; /* polygon is active */ + short int active, top; /* polygon is active, polygon is in top half */ } t_polygon; typedef struct diff --git a/heat.c b/heat.c index 4e64234..118cea8 100644 --- a/heat.c +++ b/heat.c @@ -154,7 +154,8 @@ // #define SLOPE 0.1 /* sensitivity of color on wave amplitude */ #define SLOPE 0.2 /* sensitivity of color on wave amplitude */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ -#define E_SCALE 100.0 /* scaling factor for energy representation */ +#define E_SCALE 100.0 /* scaling factor for energy representation */ +#define LOG_SCALE 1.0 /* scaling factor for energy log representation */ #define COLORHUE 260 /* initial hue of water color for scheme C_LUM */ #define COLORDRIFT 0.0 /* how much the color hue drifts during the whole simulation */ diff --git a/lennardjones.c b/lennardjones.c new file mode 100644 index 0000000..225ac6d --- /dev/null +++ b/lennardjones.c @@ -0,0 +1,897 @@ +/*********************************************************************************/ +/* */ +/* Animation of interacting particles in a planar domain */ +/* */ +/* N. Berglund, november 2021 */ +/* */ +/* UPDATE 24/04: distinction between damping and "elasticity" parameters */ +/* UPDATE 27/04: new billiard shapes, bug in color scheme fixed */ +/* UPDATE 28/04: code made more efficient, with help of Marco Mancini */ +/* */ +/* 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 lennardjones lennardjones.c */ +/* -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lXmu -lglut -O3 -fopenmp */ +/* */ +/* OMP acceleration may be more effective after executing */ +/* export OMP_NUM_THREADS=2 in the shell before running the program */ +/* */ +/* To make a video, set MOVIE to 1 and create subfolder tif_ljones */ +/* It may be possible to increase parameter PAUSE */ +/* */ +/* create movie using */ +/* ffmpeg -i lj.%05d.tif -vcodec libx264 lj.mp4 */ +/* */ +/*********************************************************************************/ + +#include +#include +#include +#include +#include +#include +#include /* Sam Leffler's libtiff library. */ +#include + +#define MOVIE 1 /* set to 1 to generate movie */ + +/* General geometrical parameters */ + +#define WINWIDTH 1280 /* window width */ +#define WINHEIGHT 720 /* window height */ + +#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 INITXMIN -2.0 +#define INITXMAX 2.0 /* x interval for initial condition */ +#define INITYMIN -1.125 +#define INITYMAX 1.125 /* y interval for initial condition */ + +#define CIRCLE_PATTERN 8 /* pattern of circles, see list in global_ljones.c */ + +#define INTERACTION 3 /* particle interaction, see list in global_ljones.c */ + +#define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */ +#define NPOISSON 100 /* number of points for Poisson C_RAND_POISSON arrangement */ +#define PDISC_DISTANCE 2.5 /* minimal distance in Poisson disc process, controls density of particles */ +#define PDISC_CANDIDATES 100 /* number of candidates in construction of Poisson disc process */ +#define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */ + +#define LAMBDA 2.0 /* parameter controlling the dimensions of domain */ +#define MU 0.02 /* parameter controlling radius of particles */ +// #define MU 0.015 /* parameter controlling radius of particles */ +#define NPOLY 3 /* number of sides of polygon */ +#define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */ +#define MDEPTH 4 /* depth of computation of Menger gasket */ +#define MRATIO 3 /* ratio defining Menger gasket */ +#define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ +#define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */ +#define FOCI 1 /* set to 1 to draw focal points of ellipse */ +// #define NGRIDX 50 /* number of grid point for grid of disks */ +#define NGRIDX 32 /* number of grid point for grid of disks */ +#define NGRIDY 26 /* number of grid point for grid of disks */ + +#define X_SHOOTER -0.2 +#define Y_SHOOTER -0.6 +#define X_TARGET 0.4 +#define Y_TARGET 0.7 /* shooter and target positions in laser fight */ + +/* Boundary conditions, see list in global_ljones.c */ + +#define B_COND 3 + +/* Parameters for length and speed of simulation */ + +#define NSTEPS 3600 /* number of frames of movie */ +// #define NSTEPS 1300 /* number of frames of movie */ +#define NVID 150 /* number of iterations between images displayed on screen */ +#define NSEG 100 /* number of segments of boundary */ +#define INITIAL_TIME 0 /* time after which to start saving frames */ +#define BOUNDARY_WIDTH 2 /* width of billiard boundary */ + +#define PAUSE 1000 /* number of frames after which to pause */ +#define PSLEEP 1 /* sleep time during pause */ +#define SLEEP1 1 /* initial sleeping time */ +#define SLEEP2 1 /* final sleeping time */ +#define END_FRAMES 100 /* number of still frames at end of movie */ + +/* Parameters of initial condition */ + + +/* Plot type, see list in global_ljones.c */ + +#define PLOT 3 + +/* Color schemes */ + +#define COLOR_PALETTE 10 /* Color palette, see list in global_ljones.c */ + +#define BLACK 1 /* background */ + +#define COLOR_SCHEME 1 /* choice of color scheme, see list in global_ljones.c */ + +#define SCALE 0 /* set to 1 to adjust color scheme to variance of field */ +#define SLOPE 0.5 /* sensitivity of color on wave amplitude */ +#define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ + +#define COLORHUE 260 /* initial hue of water color for scheme C_LUM */ +#define COLORDRIFT 0.0 /* how much the color hue drifts during the whole simulation */ +#define LUMMEAN 0.5 /* amplitude of luminosity variation for scheme C_LUM */ +#define LUMAMP 0.3 /* amplitude of luminosity variation for scheme C_LUM */ +#define HUEMEAN 220.0 /* mean value of hue for color scheme C_HUE */ +#define HUEAMP -50.0 /* amplitude of variation of hue for color scheme C_HUE */ + +/* particle properties */ + +#define PARTICLE_HUE_MIN 330.0 /* color of original particle */ +#define PARTICLE_HUE_MAX 30.0 /* color of saturated particle */ +#define PARTICLE_EMAX 2.0e1 /* max energy for particle to survive */ + +#define RANDOM_RADIUS 0 /* set to 1 for random circle radius */ +#define MOVE_PARTICLES 1 /* set to 1 for mobile particles */ +#define INERTIA 1 /* set to 1 for taking inertia into account */ +#define DT_PARTICLE 2.0e-6 /* time step for particle displacement */ +#define KREPEL 20.0 /* constant in repelling force between particles */ +#define EQUILIBRIUM_DIST 5.0 /* Lennard-Jones equilibrium distance */ +// #define EQUILIBRIUM_DIST 15.0 /* Lennard-Jones equilibrium distance */ +#define REPEL_RADIUS 20.0 /* radius in which repelling force acts (in units of particle radius) */ +#define DAMPING 1.5e5 /* damping coefficient of particles */ +// #define DAMPING 1.0e-10 /* damping coefficient of particles */ +#define PARTICLE_MASS 1.0 /* mass of particle of radius MU */ +// #define V_INITIAL 0.0 /* initial velocity range */ +#define V_INITIAL 5.0 /* initial velocity range */ +#define SIGMA 5.0 /* noise intensity in thermostat */ +#define BETA 5.0e-3 /* initial inverse temperature */ +#define MU_XI 0.1 /* friction constant in thermostat */ +#define KSPRING_BOUNDARY 1.0e5 /* confining harmonic potential outside simulation region */ +#define NBH_DIST_FACTOR 4.5 /* radius in which to count neighbours */ + +#define INCREASE_BETA 1 /* set to 1 to increase BETA during simulation */ +#define BETA_FACTOR 5.0e2 /* factor by which to change BETA during simulation */ +#define N_TOSCILLATIONS 2.5 /* number of temperature oscillations in BETA schedule */ +// #define BETA_FACTOR 2.0e3 /* factor by which to change BETA during simulation */ + +#define INCREASE_KREPEL 0 /* set to 1 to increase KREPEL during simulation */ +#define KREPEL_FACTOR 1000.0 /* factor by which to change KREPEL during simulation */ + +#define ADD_PARTICLES 0 /* set to 1 to add particles */ +#define ADD_TIME 300 /* time at which to add first particle */ +#define ADD_PERIOD 2000 /* time interval between adding further particles */ + +#define FLOOR_FORCE 0 /* set to 1 to limit force on particle to FMAX */ +#define FMAX 2.0e10 /* maximal force */ + +#define HASHX 32 /* size of hashgrid in x direction */ +#define HASHY 18 /* size of hashgrid in y direction */ +#define HASHMAX 100 /* maximal number of particles per hashgrid cell */ +#define HASHGRID_PADDING 0.1 /* padding of hashgrid outside simulation window */ + +#define DRAW_COLOR_SCHEME 0 /* set to 1 to plot the color scheme */ +#define COLORBAR_RANGE 8.0 /* scale of color scheme bar */ +#define COLORBAR_RANGE_B 12.0 /* scale of color scheme bar for 2nd part */ +#define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ + +/* For debugging purposes only */ +#define FLOOR 1 /* set to 1 to limit wave amplitude to VMAX */ +#define VMAX 10.0 /* max value of wave amplitude */ + +#include "global_ljones.c" +#include "sub_lj.c" + + +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; +} + + +/*********************/ +/* animation part */ +/*********************/ + +void hash_xy_to_ij(double x, double y, int ij[2]) +{ + static int first = 1; + static double lx, ly; + int i, j; + + if (first) + { + lx = XMAX - XMIN + 2.0*HASHGRID_PADDING; + ly = YMAX - YMIN + 2.0*HASHGRID_PADDING; + first = 0; + } + + i = (int)((double)HASHX*(x - XMIN + HASHGRID_PADDING)/lx); + j = (int)((double)HASHY*(y - YMIN + HASHGRID_PADDING)/ly); + + if (i<0) i = 0; + if (i>=HASHX) i = HASHX-1; + if (j<0) j = 0; + if (j>=HASHY) j = HASHY-1; + + ij[0] = i; + ij[1] = j; + +// printf("Mapped (%.3f,%.3f) to (%i, %i)\n", x, y, ij[0], ij[1]); +} + +double lennard_jones_force_aniso(double r, double req) +{ + int i; + double rmin = 0.01, rplus, ratio = 1.0; + + if (r > REPEL_RADIUS*MU) return(0.0); + else + { + if (r > rmin) rplus = r; + else rplus = rmin; + + for (i=0; i<6; i++) ratio *= req*MU/rplus; + + return((ratio - 2.0*ratio*ratio)/rplus); + } +} + +double lennard_jones_force(double r) +{ + int i; + double rmin = 0.01, rplus, ratio = 1.0; + + if (r > REPEL_RADIUS*MU) return(0.0); + else + { + if (r > rmin) rplus = r; + else rplus = rmin; + +// ratio = pow(EQUILIBRIUM_DIST*MU/rplus, 6.0); + for (i=0; i<6; i++) ratio *= EQUILIBRIUM_DIST*MU/rplus; + + return((ratio - 2.0*ratio*ratio)/rplus); + } +} + +void aniso_lj_force(double r, double ca, double sa, double force[2]) +{ + int i; + double rmin = 0.01, rplus, ratio = 1.0, c2, s2, c4, s4, a, aprime, f1, f2; + + if (r > REPEL_RADIUS*MU) + { + force[0] = 0.0; + force[1] = 0.0; + } + else + { + if (r > rmin) rplus = r; + else rplus = rmin; + + for (i=0; i<6; i++) ratio *= EQUILIBRIUM_DIST*MU/rplus; + + /* cos(2phi) and sin(2phi) */ + c2 = ca*ca - sa*sa; + s2 = 2.0*ca*sa; + + /* cos(4phi) and sin(4phi) */ + c4 = c2*c2 - s2*s2; + s4 = 2.0*c2*s2; + + a = 0.5*(9.0 - 7.0*c4); + aprime = 14.0*s4; + + f1 = ratio*(a - ratio)/rplus; + f2 = ratio*aprime/rplus; + + force[0] = f1*ca - f2*sa; + force[1] = f1*sa + f2*ca; + } +} + +void penta_lj_force(double r, double ca, double sa, double force[2]) +{ + int i; + double rmin = 0.01, rplus, ratio = 1.0, c2, s2, c4, s4, c5, s5, a, aprime, f1, f2; + static double a0, b0; + static int first = 1; + + if (first) + { + a0 = cos(0.1*PI) + 0.5; + b0 = a0 - 1.0; + first = 0; + } + + if (r > REPEL_RADIUS*MU) + { + force[0] = 0.0; + force[1] = 0.0; + } + else + { + if (r > rmin) rplus = r; + else rplus = rmin; + + for (i=0; i<6; i++) ratio *= EQUILIBRIUM_DIST*MU/rplus; + + /* cos(2phi) and sin(2phi) */ + c2 = ca*ca - sa*sa; + s2 = 2.0*ca*sa; + + /* cos(4phi) and sin(4phi) */ + c4 = c2*c2 - s2*s2; + s4 = 2.0*c2*s2; + + /* cos(5phi) and sin(5phi) */ + c5 = ca*c4 - sa*s4; + s5 = sa*c4 + ca*s4; + + a = a0 - b0*c5; + aprime = 5.0*b0*s5; + + f1 = ratio*(a - ratio)/rplus; + f2 = ratio*aprime/rplus; + + force[0] = f1*ca - f2*sa; + force[1] = f1*sa + f2*ca; + } +} + + +int compute_repelling_force(int i, int j, double force[2], t_particle* particle, double krepel) +/* compute repelling force of particle j on particle i */ +/* returns 1 if distance between particles is smaller than NBH_DIST_FACTOR*MU */ +{ + double x1, y1, x2, y2, distance, r, f, angle, ca, sa, aniso, fx, fy, ff[2]; + + x1 = particle[i].xc; + y1 = particle[i].yc; + x2 = particle[j].xc; + y2 = particle[j].yc; + + distance = module2(x2 - x1, y2 - y1); + + if (distance == 0.0) + { + force[0] = 0.0; + force[1] = 0.0; + return(1); + } + else + { + ca = (x2 - x1)/distance; + sa = (y2 - y1)/distance; + + switch (INTERACTION) { + case (I_COULOMB): + { + f = krepel/(1.0e-8 + distance*distance); + force[0] = f*ca; + force[1] = f*sa; + break; + } + case (I_LENNARD_JONES): + { + f = krepel*lennard_jones_force(distance); + force[0] = f*ca; + force[1] = f*sa; + break; + } + case (I_LJ_DIRECTIONAL): + { + aniso_lj_force(distance, ca, sa, ff); + force[0] = krepel*ff[0]; + force[1] = krepel*ff[1]; + break; + } + case (I_LJ_PENTA): + { + penta_lj_force(distance, ca, sa, ff); + force[0] = krepel*ff[0]; + force[1] = krepel*ff[1]; + break; + } + } + } + + if ((distance < NBH_DIST_FACTOR*MU)&&(j != i)) return(1); + else return(0); +} + + +void update_hashgrid(t_particle* particle, int* hashgrid_number, int* hashgrid_particles) +{ + int i, j, k, n, m, ij[2], max = 0; + +// printf("Updating hashgrid_number\n"); + for (i=0; i max) max = n; +// printf("Placed particle %i at (%i,%i) in hashgrid\n", k, ij[0], ij[1]); +// printf("%i particles at (%i,%i)\n", hashgrid_number[ij[0]][ij[1]], ij[0], ij[1]); + } + + printf("Maximal number of particles per hash cell: %i\n", max); +} + +void add_particle(double x, double y, double vx, double vy, t_particle particle[NMAXCIRCLES]) +{ + int i = ncircles; + + particle[i].xc = x; + particle[i].yc = y; + particle[i].radius = MU; + particle[i].active = 1; + particle[i].neighb = 0; + + particle[i].energy = 0.0; + + if (RANDOM_RADIUS) particle[i].radius = particle[i].radius*(0.75 + 0.5*((double)rand()/RAND_MAX)); + + particle[i].mass_inv = 1.0; + + particle[i].vx = vx; + particle[i].vy = vy; + particle[i].energy = (particle[i].vx*particle[i].vx + particle[i].vy*particle[i].vy)*particle[i].mass_inv; + + ncircles++; +} + +double neighbour_color(int nnbg) +{ + if (nnbg > 7) nnbg = 7; + switch(nnbg){ + case (7): return(350.0); + case (6): return(320.0); + case (5): return(260.0); + case (4): return(200.0); + case (3): return(140.0); + case (2): return(100.0); + case (1): return(70.0); + default: return(30.0); + } +} + +void draw_particles(t_particle particle[NMAXCIRCLES]) +{ + int j, k, m, width, nnbg; + double ej, hue, rgb[3], radius, x1, y1, x2, y2, angle; + + for (j=0; j 0.0) + { + hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*ej/PARTICLE_EMAX; + if (hue > PARTICLE_HUE_MIN) hue = PARTICLE_HUE_MIN; + if (hue < PARTICLE_HUE_MAX) hue = PARTICLE_HUE_MAX; + } + radius = particle[j].radius; + width = BOUNDARY_WIDTH; + break; + } + case (P_NEIGHBOURS): + { + hue = neighbour_color(particle[j].neighb); + radius = particle[j].radius; + width = BOUNDARY_WIDTH; + break; + } + case (P_BONDS): + { + glLineWidth(BOUNDARY_WIDTH); + hue = neighbour_color(particle[j].neighb); + radius = 0.015; +// radius = particle[j].radius; +// radius = 0.8*particle[j].radius; + width = 1; + for (k = 0; k < particle[j].neighb; k++) + { + m = particle[j].neighbours[k]; + angle = particle[j].nghangle[k]; + x1 = particle[j].xc + radius*cos(angle); + y1 = particle[j].yc + radius*sin(angle); + x2 = particle[m].xc - radius*cos(angle); + y2 = particle[m].yc - radius*sin(angle); + draw_line(x1, y1, x2, y2); +// draw_line(particle[j].xc, particle[j].yc, particle[m].xc, particle[m].yc); + } + break; + } + } + + hsl_to_rgb(hue, 0.9, 0.5, rgb); + draw_colored_circle(particle[j].xc, particle[j].yc, radius, NSEG, rgb); + + glLineWidth(width); + glColor3f(1.0, 1.0, 1.0); + draw_circle(particle[j].xc, particle[j].yc, radius, NSEG); + } +} + + +double repel_schedule(int i) +{ + static double kexponent; + static int first = 1; + double krepel; + + if (first) + { + kexponent = log(KREPEL_FACTOR)/(double)(INITIAL_TIME + NSTEPS); + first = 0; + } + krepel = KREPEL*exp(kexponent*(double)i); + printf("krepel = %.3lg\n", krepel); + return(krepel); +} + + +double temperature_schedule(int i) +{ + static double bexponent, omega; + static int first = 1; + double beta; + + if (first) + { + bexponent = log(BETA_FACTOR)/(double)(INITIAL_TIME + NSTEPS); + omega = N_TOSCILLATIONS*DPI/(double)(INITIAL_TIME + NSTEPS); + first = 0; + } + beta = BETA*exp(bexponent*(double)i); + beta = beta*2.0/(1.0 + cos(omega*(double)i)); + printf("beta = %.3lg\n", beta); + return(beta); +} + +void animation() +{ + double time, scale, diss, rgb[3], dissip, gradient[2], x, y, dx, dy, dt, xleft, xright, a, b, + length, fx, fy, force[2], totalenergy = 0.0, krepel = KREPEL, pos[2], + beta = BETA, xi = 0.0; + double *qx, *qy, *px, *py; + int i, j, k, n, m, s, ij[2], i0, iplus, iminus, j0, jplus, jminus, p, q, total_neighbours = 0, min_nb, max_nb, close; + static int imin, imax; + static short int first = 1; + t_particle *particle; + int *hashgrid_number, *hashgrid_particles; + t_hashgrid *hashgrid; + char message[100]; + + + particle = (t_particle *)malloc(NMAXCIRCLES*sizeof(t_particle)); /* particles */ + + hashgrid = (t_hashgrid *)malloc(HASHX*HASHY*sizeof(t_hashgrid)); /* hashgrid */ + + hashgrid_number = (int *)malloc(HASHX*HASHY*sizeof(int)); /* total number of particles in each hash grid cell */ + hashgrid_particles = (int *)malloc(HASHX*HASHY*HASHMAX*sizeof(int)); /* numbers of particles in each hash grid cell */ + + qx = (double *)malloc(NMAXCIRCLES*sizeof(double)); + qy = (double *)malloc(NMAXCIRCLES*sizeof(double)); + px = (double *)malloc(NMAXCIRCLES*sizeof(double)); + py = (double *)malloc(NMAXCIRCLES*sizeof(double)); + + /* initialise positions and radii of circles */ + init_particle_config(particle); + + + /* initialise particles */ + for (i=0; i < NMAXCIRCLES; i++) + { + particle[i].neighb = 0; + + particle[i].energy = 0.0; + y = particle[i].yc; + if (y >= YMAX) y -= particle[i].radius; + if (y <= YMIN) y += particle[i].radius; + + if (RANDOM_RADIUS) particle[i].radius = particle[i].radius*(0.75 + 0.5*((double)rand()/RAND_MAX)); + + particle[i].mass_inv = 1.0; + + particle[i].vx = V_INITIAL*gaussian(); + particle[i].vy = V_INITIAL*gaussian(); + particle[i].energy = (particle[i].vx*particle[i].vx + particle[i].vy*particle[i].vy)*particle[i].mass_inv; + + px[i] = particle[i].vx; + py[i] = particle[i].vy; + } + for (i=ncircles; i < NMAXCIRCLES; i++) + { + particle[i].active = 0; + particle[i].neighb = 0; + particle[i].energy = 0.0; + particle[i].mass_inv = 1.0; + particle[i].vx = 0.0; + particle[i].vy = 0.0; + px[i] = 0.0; + py[i] = 0.0; + } + + xi = 0.0; + + for (i=0; i= HASHX) iplus = HASHX-1; + + j0 = particle[j].hashy; + jminus = j0 - 1; if (jminus < 0) jminus = 0; + jplus = j0 + 1; if (jplus >= HASHY) jplus = HASHY-1; + + fx = 0.0; + fy = 0.0; + for (p=iminus; p<= iplus; p++) + for (q=jminus; q<= jplus; q++) + for (k=0; k XMAX) fx -= KSPRING_BOUNDARY*(particle[j].xc - XMAX); + else if (particle[j].xc < XMIN) fx += KSPRING_BOUNDARY*(XMIN - particle[j].xc); + if (particle[j].yc > YMAX) fy -= KSPRING_BOUNDARY*(particle[j].yc - YMAX); + else if (particle[j].yc < YMIN) fy += KSPRING_BOUNDARY*(YMIN - particle[j].yc); + + if (FLOOR_FORCE) + { + if (fx > FMAX) fx = FMAX; + if (fx < -FMAX) fx = -FMAX; + if (fy > FMAX) fy = FMAX; + if (fy < -FMAX) fy = -FMAX; + } + + particle[j].fx = fx; + particle[j].fy = fy; + } + + /* timestep of thermostat algorithm */ + for (j=0; j max_nb) max_nb = particle[j].neighb; + if (particle[j].neighb < min_nb) min_nb = particle[j].neighb; + } + printf("Mean number of neighbours: %.3f\n", (double)total_neighbours/(double)ncircles); + printf("Min number of neighbours: %i\n", min_nb); + printf("Max number of neighbours: %i\n", max_nb); + + draw_particles(particle); + + /* add a particle */ + if ((ADD_PARTICLES)&&((i - ADD_TIME + 1)%ADD_PERIOD == 0)) + { + j = 0; + while (module2(particle[j].xc,particle[j].yc) > 0.7) j = rand()%ncircles; + x = particle[j].xc + 2.5*MU; + y = particle[j].yc; + add_particle(x, y, 0.0, 0.0, particle); + } + + update_hashgrid(particle, hashgrid_number, hashgrid_particles); + + if (INCREASE_BETA) /* print force constant */ + { + erase_area_hsl(XMAX - 0.39, YMAX - 0.1 + 0.025, 0.37, 0.05, 0.0, 0.9, 0.0); + glColor3f(1.0, 1.0, 1.0); + sprintf(message, "Temperature %.3f", 1.0/beta); + write_text(XMAX - 0.7, YMAX - 0.1, message); + } + else if (INCREASE_KREPEL) /* print force constant */ + { + erase_area_hsl(XMAX - 0.24, YMAX - 0.1 + 0.025, 0.22, 0.05, 0.0, 0.9, 0.0); + glColor3f(1.0, 1.0, 1.0); + sprintf(message, "Force %.0f", krepel); + write_text(XMAX - 0.42, YMAX - 0.1, message); + } + + glutSwapBuffers(); + + + if (MOVIE) + { + if (i >= INITIAL_TIME) save_frame_lj(); + else printf("Initial phase time %i of %i\n", i, INITIAL_TIME); + + /* it seems that saving too many files too fast can cause trouble with the file system */ + /* so this is to make a pause from time to time - parameter PAUSE may need adjusting */ + if (i % PAUSE == PAUSE - 1) + { + printf("Making a short pause\n"); + sleep(PSLEEP); + s = system("mv lj*.tif tif_ljones/"); + } + } + + } + + if (MOVIE) + { + for (i=0; i=0) color[i]++; if ((!RAINBOW_COLOR)&&(c>=0)) color[i]++; if (!RAINBOW_COLOR) @@ -606,13 +613,25 @@ void print_part_number(double *configs[NPARTMAX], int active[NPARTMAX], double x } - +void print_collision_number(int ncollisions, double x, double y) +{ + char message[50]; + double rgb[3]; + + hsl_to_rgb(0.0, 0.0, 0.0, rgb); + erase_area(x, y, 0.5, 0.1, rgb); + + glColor3f(1.0, 1.0, 1.0); + sprintf(message, "%i collisions", ncollisions); + write_text(x, y, message); + +} void animation() { - double time, dt, alpha, r, rgb[3]; + double time, dt, alpha, r, rgb[3], alphamax; double *configs[NPARTMAX]; - int i, j, resamp = 1, s, i1, i2; + int i, j, resamp = 1, s, i1, i2, c, lengthmax; int *color, *newcolor, *active; // t_circle *circles; /* experimental */ @@ -630,7 +649,15 @@ void animation() ||(B_DOMAIN == D_CIRCLES_IN_TORUS)) init_circles(circles); else if (B_DOMAIN == D_POLYLINE) init_polyline(polyline, circles); - + + if (POLYLINE_PATTERN == P_TOKA_PRIME) + { + x_shooter = -polyline[84].x1; + y_shooter = polyline[84].y1; + x_target = polyline[84].x1; + y_target = polyline[84].y1; + } + /* initialize system by putting particles in a given point with a range of velocities */ r = cos(PI/(double)NPOLY)/cos(DPI/(double)NPOLY); @@ -642,8 +669,35 @@ void animation() // init_drop_config(-1.0 + 0.3*sqrt(2.0), -1.0 + 0.5*sqrt(2.0), 0.0, DPI, configs); // init_line_config(-1.25, -0.5, -1.25, 0.5, 0.0, configs); -// init_drop_config(0.0, 0.0, -PI, PI, configs); - init_drop_config(-0.95, 0.0, PI, 3.0*PI, configs); +// init_drop_config(-0.95, 0.0, -0.103 + DPI/15.0, -0.1 + DPI/15.0, configs); + + /* find long trajectory */ +// alphamax = 0.0; +// lengthmax = 1; +// for (alpha = 0.0; alpha < DPI; alpha += 0.00001) +// { +// init_drop_config(x_shooter, y_shooter, alpha, alpha + DPI, configs); +// i = 0; +// c = 1; +// while ((c >= 0)&&(i<=1000)) +// { +// c = vbilliard(configs[0]); +// i++; +// } +// if (i > 100) printf("Angle %.6lg, length %i\n", alpha, i); +// if (i > lengthmax) +// { +// lengthmax = i; +// alphamax = alpha; +// } +// } +// printf("Angle %.6lg, max length %i\n", alphamax, lengthmax); + + alphamax = 2.50949; + init_drop_config(x_shooter, y_shooter, alphamax, alphamax + DPI, configs); + + +// init_drop_config(-0.95, 0.0, 1.0, 1.0 + DPI, configs); // init_drop_config(-1.3, -0.1, 0.0, DPI, configs); // init_drop_config(1.4, 0.1, 0.0, DPI, configs); // init_drop_config(0.5, 0.5, -1.0, 1.0, configs); @@ -663,7 +717,8 @@ void animation() glColor3f(0.0, 0.0, 0.0); if (DRAW_BILLIARD) draw_billiard(); if (PRINT_PARTICLE_NUMBER) print_part_number(configs, active, XMIN + 0.1, YMIN + 0.1); - + else if (PRINT_COLLISION_NUMBER) print_collision_number(ncollisions, XMIN + 0.1, YMIN + 0.1); + glutSwapBuffers(); @@ -716,6 +771,8 @@ void animation() // draw_config(newcolor, configs, active); if (DRAW_BILLIARD) draw_billiard(); if (PRINT_PARTICLE_NUMBER) print_part_number(configs, active, XMIN + 0.1, YMIN + 0.1); + else if (PRINT_COLLISION_NUMBER) print_collision_number(ncollisions, XMIN + 0.1, YMIN + 0.1); + for (j=0; j= 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; + } + 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 save_frame_lj() +{ + static int counter = 0; + char *name="lj.", 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, "Particles with Lennard-Jones interaction in a planar domain", 0, 0, + WINWIDTH, WINHEIGHT, COMPRESSION_LZW); + +} + +void save_frame_lj_counter(int counter) +{ + char *name="lj.", 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, "Particles with Lennard-Jones interaction 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); + } + + + +/*********************/ +/* drawing routines */ +/*********************/ + + + +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); + glVertex2d(x - dx, y - dy); + glVertex2d(x + dx, y - dy); + glVertex2d(x + dx, y + dy); + glVertex2d(x - dx, y + dy); + 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); + glVertex2d(x - dx, y - dy); + glVertex2d(x + dx, y - dy); + glVertex2d(x + dx, y + dy); + glVertex2d(x - dx, y + dy); + 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) +{ + glBegin(GL_LINE_STRIP); + glVertex2d(x1, y1); + glVertex2d(x2, y2); + glEnd(); +} + +void draw_rectangle(double x1, double y1, double x2, double y2) +{ + glBegin(GL_LINE_LOOP); + glVertex2d(x1, y1); + glVertex2d(x2, y1); + glVertex2d(x2, y2); + glVertex2d(x1, y2); + glEnd(); +} + +void draw_colored_triangle(double x1, double y1, double x2, double y2, double x3, double y3, double rgb[3]) +{ + glColor3f(rgb[0], rgb[1], rgb[2]); + glBegin(GL_TRIANGLE_FAN); + glVertex2d(x1, y1); + glVertex2d(x2, y2); + glVertex2d(x3, y3); + 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; + glVertex2d(x + r*cos(alpha), y + r*sin(alpha)); + } + 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); + glVertex2d(x, y); + for (i=0; i<=nseg; i++) + { + alpha = (double)i*dalpha; + glVertex2d(x + r*cos(alpha), y + r*sin(alpha)); + } + + glEnd(); +} + + + +void init_particle_config(t_particle particles[NMAXCIRCLES]) +/* initialise particle configuration */ +{ + int i, j, k, n, ncirc0, n_p_active, ncandidates = PDISC_CANDIDATES, naccepted; + double dx, dy, p, phi, r, r0, ra[5], sa[5], height, x, y = 0.0, gamma, dpoisson = PDISC_DISTANCE*MU, xx[4], yy[4]; + 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; + particles[n].xc = ((double)(i-NGRIDX/2) + 0.5)*dy; + particles[n].yc = YMIN + ((double)j + 0.5)*dy; + particles[n].radius = MU; + particles[n].active = 1; + } + break; + } + case (C_HEX): + { + ncircles = NGRIDX*(NGRIDY+1); + dx = (INITXMAX - INITXMIN)/((double)NGRIDX); + dy = (INITYMAX - INITYMIN)/((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; + particles[n].xc = ((double)(i-NGRIDX/2) + 0.5)*dx; /* is +0.5 needed? */ + particles[n].yc = INITYMIN + ((double)j - 0.5)*dy; + if ((i+NGRIDX)%2 == 1) particles[n].yc += 0.5*dy; + particles[n].radius = MU; + /* activate only circles that intersect the domain */ + if ((particles[n].yc < INITYMAX + MU)&&(particles[n].yc > INITYMIN - MU)) particles[n].active = 1; + else particles[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; + particles[n].xc = ((double)(i-NGRIDX/2) + 0.5*((double)rand()/RAND_MAX - 0.5))*dy; + particles[n].yc = YMIN + ((double)j + 0.5 + 0.5*((double)rand()/RAND_MAX - 0.5))*dy; + particles[n].radius = MU; +// particles[n].radius = MU*sqrt(1.0 + 0.8*((double)rand()/RAND_MAX - 0.5)); + particles[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; + particles[n].xc = ((double)(i-NGRIDX/2) + 0.5)*dy; + particles[n].yc = YMIN + ((double)j + 0.5)*dy; + particles[n].radius = MU; + p = (double)rand()/RAND_MAX; + if (p < P_PERCOL) particles[n].active = 1; + else particles[n].active = 0; + } + break; + } + case (C_RAND_POISSON): + { + ncircles = NPOISSON; + for (n = 0; n < NPOISSON; n++) + { +// particles[n].xc = LAMBDA*(2.0*(double)rand()/RAND_MAX - 1.0); + particles[n].xc = (XMAX - XMIN)*(double)rand()/RAND_MAX + XMIN; + particles[n].yc = (YMAX - YMIN)*(double)rand()/RAND_MAX + YMIN; + particles[n].radius = MU; + particles[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); + particles[n].xc = r*cos(phi); + particles[n].yc = r*sin(phi); + particles[n].radius = MU; + particles[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]; + particles[n].xc = r*cos(phi); + particles[n].yc = r*sin(phi); + particles[n].radius = LAMBDA*ra[j]; + particles[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++) + { + particles[4*i + j].xc = xx[i]; + particles[4*i + j].yc = yy[j]; + + } + + particles[ncircles - 1].xc = X_TARGET; + particles[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, particles[i].xc, particles[i].yc); + /* generate new candidates */ + naccepted = 0; + for (j=0; j= dpoisson*dpoisson); + /* new circle is in domain */ + far = far*(x < INITXMAX)*(x > INITXMIN)*(y < INITYMAX)*(y > INITYMIN); +// far = far*(vabs(x) < LAMBDA)*(y < INITYMAX)*(y > INITYMIN); + } + if (far) /* accept new circle */ + { + printf("New circle at (%.3f,%.3f) accepted\n", x, y); + particles[ncircles].xc = x; + particles[ncircles].yc = y; + particles[ncircles].radius = MU; + particles[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++) + { + particles[n].xc = -LAMBDA + n*dx; + particles[n].yc = y; + y += height*gamma; + if (y > YMAX) y -= height; + particles[n].radius = MU; + particles[n].active = 1; + } + + /* test for circles that overlap top or bottom boundary */ + ncirc0 = ncircles; + for (n=0; n < ncirc0; n++) + { + if (particles[n].yc + particles[n].radius > YMAX) + { + particles[ncircles].xc = particles[n].xc; + particles[ncircles].yc = particles[n].yc - height; + particles[ncircles].radius = MU; + particles[ncircles].active = 1; + ncircles ++; + } + else if (particles[n].yc - particles[n].radius < YMIN) + { + particles[ncircles].xc = particles[n].xc; + particles[ncircles].yc = particles[n].yc + height; + particles[ncircles].radius = MU; + particles[ncircles].active = 1; + ncircles ++; + } + } + break; + } + case (C_GOLDEN_SPIRAL): + { + ncircles = 1; + particles[0].xc = 0.0; + particles[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)) + { + particles[ncircles].xc = x; + particles[ncircles].yc = y; + ncircles++; + } + } + + for (i=0; i YMIN - MU)) particles[i].active = 1; +// printf("i = %i, circlex = %.3lg, circley = %.3lg\n", i, particles[i].xc, particles[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; + particles[n].xc = ((double)(i-NGRIDX/2) + 0.5)*dy; /* is +0.5 needed? */ + particles[n].yc = YMIN + ((double)j - 0.5)*dy; + if (((i+NGRIDX)%4 == 2)||((i+NGRIDX)%4 == 3)) particles[n].yc += 0.5*dy; + particles[n].radius = MU; + /* activate only circles that intersect the domain */ + if ((particles[n].yc < YMAX + MU)&&(particles[n].yc > YMIN - MU)) particles[n].active = 1; + else particles[n].active = 0; + } + break; + } + case (C_ONE): + { + particles[ncircles].xc = 0.0; + particles[ncircles].yc = 0.0; + particles[ncircles].radius = MU; + particles[ncircles].active = 1; + ncircles += 1; + break; + } + case (C_TWO): /* used for comparison with cloak */ + { + particles[ncircles].xc = 0.0; + particles[ncircles].yc = 0.0; + particles[ncircles].radius = MU; + particles[ncircles].active = 2; + ncircles += 1; + + particles[ncircles].xc = 0.0; + particles[ncircles].yc = 0.0; + particles[ncircles].radius = 2.0*MU; + particles[ncircles].active = 1; + ncircles += 1; + break; + } + case (C_NOTHING): + { + ncircles += 0; + break; + } + default: + { + printf("Function init_circle_config not defined for this pattern \n"); + } + } +} + +void init_people_config(t_person people[NMAXCIRCLES]) +/* initialise particle configuration */ +{ + t_particle particles[NMAXCIRCLES]; + int n; + + init_particle_config(particles); + + for (n=0; n= 0; i--) + for (i = height - 1; i >= 0; i--) { // if (TIFFWriteScanline(file, p, height - i - 1, 0) < 0) if (TIFFWriteScanline(file, p, i, 0) < 0) @@ -121,10 +120,10 @@ int writetiff(char *filename, char *description, int x, int y, int width, int he } p += width * sizeof(GLubyte) * 3; } - free(image); /* prenvents RAM consumption*/ TIFFClose(file); return 0; } + void init() /* initialisation of window */ { @@ -136,6 +135,8 @@ void init() /* initialisation of window */ glOrtho(XMIN, XMAX, YMIN, YMAX , -1.0, 1.0); } + + void rgb_color_scheme(int i, double rgb[3]) /* color scheme */ { double hue, y, r; @@ -171,7 +172,7 @@ void blank() if (COLOR_OUTSIDE) { - hsl_to_rgb(OUTER_COLOR, 0.9, 0.15, rgb); + hsl_to_rgb(OUTER_COLOR, 0.9, 0.15, rgb); glClearColor(rgb[0], rgb[1], rgb[2], 1.0); } else if (BLACK) glClearColor(0.0, 0.0, 0.0, 1.0); @@ -179,6 +180,7 @@ void blank() glClear(GL_COLOR_BUFFER_BIT); } + void save_frame() { static int counter = 0; @@ -191,16 +193,12 @@ void save_frame() sprintf(strstr(n2,"."), format, counter); strcat(n2, ".tif"); printf(" saving frame %s \n",n2); - - // choose one of the following according to the comment beside. - writetiff(n2, "Billiard in an ellipse", 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, "Billiard in an ellipse", 0, 0, WINWIDTH, WINHEIGHT-50, COMPRESSION_LZW); // works for 1080p -> "-50px" band!!! - // writetiff(n2, "Billiard in an ellipse", 0, 0, 1920, 1080-40, COMPRESSION_LZW); //another perfect 1080p from 1440p in setup - // writetiff(n2, "Billiard in an ellipse", -WINWIDTH/8+320, -WINHEIGHT/8+180, WINWIDTH-640, WINHEIGHT-400, COMPRESSION_LZW); // perfect 1040p from 1440p in setup + writetiff(n2, "Billiard in an ellipse", 0, 0, + WINWIDTH, WINHEIGHT, COMPRESSION_LZW); + } + void write_text_fixedwidth( double x, double y, char *st) { int l, i; @@ -1173,28 +1171,6 @@ void draw_billiard() /* draws the billiard boundary */ draw_circle(circles[k].xc, circles[k].yc, circles[k].radius, NSEG); break; } -// { -// rgb[0] = 0.0; rgb[1] = 0.0; rgb[2] = 0.0; -// for (k=0; k= 1) -// { -// rgb_color_scheme_lum(circlecolor[k], 0.85, rgb1); -// newcircle[k]--; -// } -// else rgb_color_scheme(circlecolor[k], rgb1); -// draw_colored_circle(circlex[k], circley[k], circlerad[k], NSEG, rgb1); -// } -// } -// init_billiard_color(); -// for (k=0; kx1 = z.x1 - 2.0*zperp.x1; + zprime->y1 = z.y1 - 2.0*zperp.y1; + + return(1); +} + void init_polyline(t_segment polyline[NMAXPOLY], t_circle circles[NMAXCIRCLES]) { int i, j, k, l, n, z, ii, jj, terni[SDEPTH], ternj[SDEPTH], quater[SDEPTH], cond; short int vkoch[NMAXCIRCLES], turnright; - double ratio, omega, angle, sw, length, dist, x, y; + double ratio, omega, angle, sw, length, dist, x, y, ta, tb, a, b; switch (POLYLINE_PATTERN) { case (P_RECTANGLE): @@ -5931,6 +5951,139 @@ void init_polyline(t_segment polyline[NMAXPOLY], t_circle circles[NMAXCIRCLES]) break; } + case (P_POLYGON): + { + nsides = NPOLY; + ncircles = 0; + omega = DPI/(double)NPOLY; + sw = sin(omega/2.0); + + for (i=0; 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) { @@ -68,12 +121,6 @@ void init() /* initialisation of window */ glOrtho(0.0, NX, 0.0, NY, -1.0, 1.0); } - - - - - - void blank() { if (BLACK) glClearColor(0.0, 0.0, 0.0, 1.0); @@ -82,6 +129,53 @@ void blank() } +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() { @@ -419,7 +513,7 @@ void init_circle_config(t_circle circles[NMAXCIRCLES]) /* 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]; + 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) { @@ -728,6 +822,26 @@ void init_circle_config(t_circle circles[NMAXCIRCLES]) } 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_ONE): { circles[ncircles].xc = 0.0; @@ -785,7 +899,11 @@ void init_polygon_config(t_polygon polygons[NMAXCIRCLES]) /* if (i < ncircles) printf("(x,y) = (%.2f, %.2f), r = %.2f, angle = %.2f, sides = %i\n", polygons[i].xc, polygons[i].yc, polygons[i].radius, polygons[i].angle, polygons[i].nsides);*/ } - + + /* adjust angles for C_RINGS configuration */ + if (CIRCLE_PATTERN == C_RINGS) + for (i=0; i 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 */ @@ -1522,6 +1740,13 @@ int xy_in_billiard(double x, double y) } 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_MENGER): { x1 = 0.5*(x+1.0); @@ -2247,6 +2472,56 @@ void draw_billiard() /* draws the billiard boundary */ } break; } + case (D_TOKA_PRIME): + { + glBegin(GL_LINE_LOOP); + tvertex_lineto(polyline[0]); + for (i=4; i<43; i++) tvertex_lineto(polyline[i]); + tvertex_lineto(polyline[3]); + tvertex_lineto(polyline[2]); + tvertex_lineto(polyline[1]); + + tvertex_lineto(polyline[44]); + tvertex_lineto(polyline[45]); + for (i=84; i>45; 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 */ @@ -2402,6 +2677,14 @@ void draw_billiard() /* draws the billiard boundary */ glEnd(); break; } + case (D_STAR): + { + glLineWidth(BOUNDARY_WIDTH); + glBegin(GL_LINE_LOOP); + for (i=0; i= COL_TURBO) color_scheme_asym(COLOR_SCHEME, energy, scale, time, rgb); - else color_scheme(COLOR_SCHEME, energy, scale, time, rgb); - } - else if (PLOT == P_MIXED) - { - if (j > NY/2) color_scheme(COLOR_SCHEME, phi[i][j], scale, time, rgb); - else color_scheme(COLOR_SCHEME, compute_energy(phi, psi, xy_in, i, j), scale, time, rgb); + switch (PLOT) { + case (P_AMPLITUDE): + { + color_scheme(COLOR_SCHEME, phi[i][j], scale, time, rgb); + break; + } + case (P_ENERGY): + { + energy = compute_energy(phi, psi, xy_in, i, j); + if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym(COLOR_SCHEME, energy, scale, time, rgb); + else color_scheme(COLOR_SCHEME, energy, scale, time, rgb); + break; + } + case (P_MIXED): + { + if (j > NY/2) color_scheme(COLOR_SCHEME, phi[i][j], scale, time, rgb); + else color_scheme(COLOR_SCHEME, compute_energy(phi, psi, xy_in, i, j), scale, time, rgb); + break; + } +// case (P_MEAN_ENERGY): +// { +// energy = compute_energy(phi, psi, xy_in, i, j); +// total_energy[i][j] += energy; +// if (COLOR_PALETTE >= COL_TURBO) +// color_scheme_asym(COLOR_SCHEME, total_energy[i][j]/(double)(time+1), scale, time, rgb); +// else color_scheme(COLOR_SCHEME, total_energy[i][j]/(double)(time+1), scale, time, rgb); +// break; +// } + case (P_LOG_ENERGY): + { + energy = compute_energy(phi, psi, xy_in, i, j); + color_scheme(COLOR_SCHEME, LOG_SHIFT + LOG_SCALE*log(energy), scale, time, rgb); + break; + } +// case (P_LOG_MEAN_ENERGY): +// { +// energy = compute_energy(phi, psi, xy_in, i, j); +// if (energy == 0.0) energy = 1.0e-20; +// total_energy[i][j] += energy; +// color_scheme(COLOR_SCHEME, LOG_SCALE*log(total_energy[i][j]/(double)(time+1)), scale, time, rgb); +// break; +// } } + +// if (PLOT == P_AMPLITUDE) +// color_scheme(COLOR_SCHEME, phi[i][j], scale, time, rgb); +// else if (PLOT == P_ENERGY) +// { +// energy = compute_energy(phi, psi, xy_in, i, j); +// if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym(COLOR_SCHEME, energy, scale, time, rgb); +// else color_scheme(COLOR_SCHEME, energy, scale, time, rgb); +// } +// else if (PLOT == P_MIXED) +// { +// if (j > NY/2) color_scheme(COLOR_SCHEME, phi[i][j], scale, time, rgb); +// else color_scheme(COLOR_SCHEME, compute_energy(phi, psi, xy_in, i, j), scale, time, rgb); +// } glColor3f(rgb[0], rgb[1], rgb[2]); glVertex2i(i, j); @@ -745,20 +805,29 @@ void compute_energy_tblr(double *phi[NX], double *psi[NX], short int *xy_in[NX], /* compute total energy in top/bottom left/right boxes */ { int i, j, ij[2]; - double energy = 0.0, rescale, pos, xleft = XMAX, xright = XMIN, emax = 1.2; - double energy_ij[NX][NY]; + double energy = 0.0, rescale, pos, xleft = XMAX, xright = XMIN, emax = 1.2, xc, r; + double *energy_ij[NX]; static short int first = 1; static int ileft, iright, jmid = NY/2; static double sqremax; + for (i=0; i xright)) xright = circles[i].xc + circles[i].radius; + if ((polygons[i].active)&&(polygons[i].xc + polygons[i].radius > xright)) xright = polygons[i].xc + polygons[i].radius; + } +// for (i=0; i xright)) xright = circles[i].xc + circles[i].radius; xy_to_ij(xleft, 0.0, ij); ileft = ij[0]; @@ -837,15 +906,26 @@ void compute_energy_tblr(double *phi[NX], double *psi[NX], short int *xy_in[NX], // printf("Energies: %.5lg, %.5lg, %.5lg\n %.5lg, %.5lg, %.5lg\n", energies[0], energies[1], energies[2], energies[3], energies[4], energies[5]); + for (i=0; i 1280) + { + boxheight = 0.035; +// centerx -= 0.2; + boxwidth = 0.08; + leftboxshift = 0.16; + centerboxshift = 0.06; + } ytop = YMAX - 0.1; -// ybot = -0.1; ybot = YMIN + 0.05; if (DRAW_COLOR_SCHEME) { @@ -853,42 +933,42 @@ void print_energies(double energies[6], double top_energy, double bottom_energy) textxright -= 0.35; } - erase_area(XMIN + 0.185, ytop + 0.025, 0.1, 0.05); + erase_area(XMIN + leftboxshift, ytop + 0.025, boxwidth, boxheight); glColor3f(1.0, 1.0, 1.0); sprintf(message, "%.3f", energies[0]/top_energy); xy_to_pos(XMIN + 0.1, ytop, pos); write_text(pos[0], pos[1], message); - erase_area(centerx + 0.085, ytop + 0.025, 0.1, 0.05); + erase_area(centerx + centerboxshift, ytop + 0.025, boxwidth, boxheight); glColor3f(1.0, 1.0, 1.0); sprintf(message, "%.3f", energies[1]/top_energy); xy_to_pos(centerx, ytop, pos); write_text(pos[0], pos[1], message); - erase_area(boxxright, ytop + 0.025, 0.15, 0.05); + erase_area(boxxright, ytop + 0.025, boxwidth + 0.05, boxheight); glColor3f(1.0, 1.0, 1.0); sprintf(message, "%.5f", energies[2]/top_energy); xy_to_pos(textxright, ytop, pos); write_text(pos[0], pos[1], message); - erase_area(XMIN + 0.185, ybot + 0.025, 0.1, 0.05); + erase_area(XMIN + leftboxshift, ybot + 0.025, boxwidth, boxheight); glColor3f(1.0, 1.0, 1.0); sprintf(message, "%.3f", energies[3]/bottom_energy); xy_to_pos(XMIN + 0.1, ybot, pos); write_text(pos[0], pos[1], message); - erase_area(centerx + 0.085, ybot + 0.025, 0.1, 0.05); + erase_area(centerx + centerboxshift, ybot + 0.025, boxwidth, boxheight); glColor3f(1.0, 1.0, 1.0); sprintf(message, "%.3f", energies[4]/bottom_energy); xy_to_pos(centerx, ybot, pos); write_text(pos[0], pos[1], message); - erase_area(boxxright, ybot + 0.025, 0.15, 0.05); + erase_area(boxxright, ybot + 0.025, boxwidth + 0.05, boxheight); glColor3f(1.0, 1.0, 1.0); sprintf(message, "%.5f", energies[5]/bottom_energy); diff --git a/wave_billiard.c b/wave_billiard.c index ad685c8..066c33f 100644 --- a/wave_billiard.c +++ b/wave_billiard.c @@ -42,27 +42,37 @@ #include /* Sam Leffler's libtiff library. */ #include -#define MOVIE 0 /* set to 1 to generate movie */ -#define DOUBLE_MOVIE 0 /* set to 1 to produce movies for wave height and energy simultaneously */ +#define MOVIE 1 /* set to 1 to generate movie */ +#define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */ /* General geometrical parameters */ -#define WINWIDTH 1280 /* window width */ -#define WINHEIGHT 720 /* window height */ - -#define NX 1280 /* number of grid points on x axis */ -#define NY 720 /* number of grid points on y axis */ +#define WINWIDTH 1920 /* window width */ +#define WINHEIGHT 1000 /* window height */ +#define NX 1920 /* number of grid points on x axis */ +#define NY 1000 /* number of grid points on y axis */ #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 YMIN -1.0329 +#define YMAX 1.1129 /* y interval for 9/16 aspect ratio */ + +// #define WINWIDTH 1280 /* window width */ +// #define WINHEIGHT 720 /* window height */ + +// #define NX 1280 /* number of grid points on x axis */ +// #define NY 720 /* number of grid points on y axis */ +// #define NX 2560 /* number of grid points on x axis */ +// #define NY 1440 /* number of grid points on y axis */ + +// #define YMIN -1.125 +// #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ #define JULIA_SCALE 1.0 /* scaling for Julia sets */ /* Choice of the billiard table */ -#define B_DOMAIN 38 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN 42 /* choice of domain shape, see list in global_pdes.c */ #define CIRCLE_PATTERN 2 /* pattern of circles or polygons, see list in global_pdes.c */ @@ -70,10 +80,10 @@ #define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */ #define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */ -#define LAMBDA 0.9 /* parameter controlling the dimensions of domain */ -#define MU 0.03 /* parameter controlling the dimensions of domain */ -#define NPOLY 6 /* number of sides of polygon */ -#define APOLY 2.0 /* angle by which to turn polygon, in units of Pi/2 */ +#define LAMBDA 1.1 /* parameter controlling the dimensions of domain */ +#define MU 0.4 /* parameter controlling the dimensions of domain */ +#define NPOLY 14 /* number of sides of polygon */ +#define APOLY -1.0 /* angle by which to turn polygon, in units of Pi/2 */ #define MDEPTH 5 /* depth of computation of Menger gasket */ #define MRATIO 3 /* ratio defining Menger gasket */ #define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ @@ -81,8 +91,6 @@ #define FOCI 1 /* set to 1 to draw focal points of ellipse */ #define NGRIDX 10 /* number of grid point for grid of disks */ #define NGRIDY 12 /* number of grid point for grid of disks */ -// #define NGRIDX 16 /* number of grid point for grid of disks */ -// #define NGRIDY 20 /* number of grid point for grid of disks */ #define X_SHOOTER -0.2 #define Y_SHOOTER -0.6 @@ -106,11 +114,9 @@ #define OMEGA 0.002 /* frequency of periodic excitation */ #define AMPLITUDE 1.0 /* amplitude of periodic excitation */ -#define COURANT 0.02 /* Courant number */ +#define COURANT 0.04 /* Courant number */ #define COURANTB 0.01 /* Courant number in medium B */ #define GAMMA 0.0 /* damping factor in wave equation */ -// #define GAMMAB 5.0e-3 /* damping factor in wave equation */ -// #define GAMMAB 1.0e-2 /* damping factor in wave equation */ #define GAMMAB 1.0e-6 /* damping factor in wave equation */ #define GAMMA_SIDES 1.0e-4 /* damping factor on boundary */ #define GAMMA_TOPBOT 1.0e-7 /* damping factor on boundary */ @@ -128,12 +134,11 @@ /* Parameters for length and speed of simulation */ -// #define NSTEPS 500 /* number of frames of movie */ -#define NSTEPS 1500 /* number of frames of movie */ -#define NVID 40 /* number of iterations between images displayed on screen */ +#define NSTEPS 2100 /* number of frames of movie */ +#define NVID 25 /* number of iterations between images displayed on screen */ #define NSEG 100 /* number of segments of boundary */ #define INITIAL_TIME 0 /* time after which to start saving frames */ -#define BOUNDARY_WIDTH 2 /* width of billiard boundary */ +#define BOUNDARY_WIDTH 3 /* width of billiard boundary */ #define PAUSE 1000 /* number of frames after which to pause */ #define PSLEEP 1 /* sleep time during pause */ @@ -145,8 +150,6 @@ /* Parameters of initial condition */ #define INITIAL_AMP 0.75 /* amplitude of initial condition */ -// #define INITIAL_VARIANCE 0.0003 /* variance of initial condition */ -// #define INITIAL_WAVELENGTH 0.015 /* wavelength of initial condition */ #define INITIAL_VARIANCE 0.0003 /* variance of initial condition */ #define INITIAL_WAVELENGTH 0.015 /* wavelength of initial condition */ @@ -158,32 +161,29 @@ /* Color schemes */ -#define COLOR_PALETTE 14 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE 12 /* Color palette, see list in global_pdes.c */ #define BLACK 1 /* background */ #define COLOR_SCHEME 3 /* choice of color scheme, see list in global_pdes.c */ #define SCALE 0 /* set to 1 to adjust color scheme to variance of field */ -// #define SLOPE 0.25 /* sensitivity of color on wave amplitude */ -#define SLOPE 1.0 /* sensitivity of color on wave amplitude */ +#define SLOPE 0.5 /* sensitivity of color on wave amplitude */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ -// #define E_SCALE 150.0 /* scaling factor for energy representation */ -#define E_SCALE 100.0 /* scaling factor for energy representation */ +#define E_SCALE 150.0 /* scaling factor for energy representation */ +#define LOG_SCALE 2.5 /* scaling factor for energy log representation */ #define COLORHUE 260 /* initial hue of water color for scheme C_LUM */ #define COLORDRIFT 0.0 /* how much the color hue drifts during the whole simulation */ #define LUMMEAN 0.5 /* amplitude of luminosity variation for scheme C_LUM */ #define LUMAMP 0.3 /* amplitude of luminosity variation for scheme C_LUM */ #define HUEMEAN 180.0 /* mean value of hue for color scheme C_HUE */ -// #define HUEMEAN 210.0 /* mean value of hue for color scheme C_HUE */ #define HUEAMP -180.0 /* amplitude of variation of hue for color scheme C_HUE */ -// #define HUEAMP -180.0 /* amplitude of variation of hue for color scheme C_HUE */ #define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */ -#define COLORBAR_RANGE 2.0 /* scale of color scheme bar */ -#define COLORBAR_RANGE_B 12.0 /* scale of color scheme bar for 2nd part */ -#define ROTATE_COLOR_SCHEME 1 /* set to 1 to draw color scheme horizontally */ +#define COLORBAR_RANGE 6.0 /* scale of color scheme bar */ +#define COLORBAR_RANGE_B 10.0 /* scale of color scheme bar for 2nd part */ +#define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ #define SAVE_TIME_SERIES 0 /* set to 1 to save wave time series at a point */ @@ -606,7 +606,7 @@ void draw_color_bar(int plot, double range) void animation() { - double time, scale, ratio, startleft[2], startright[2]; + double time, scale, ratio, startleft[2], startright[2], sign; double *phi[NX], *psi[NX], *phi_tmp[NX], *psi_tmp[NX], *total_energy[NX]; short int *xy_in[NX]; int i, j, s, sample_left[2], sample_right[2]; @@ -635,7 +635,7 @@ void animation() else if (B_DOMAIN == D_POLYGONS) init_polygon_config(polygons); printf("Polygons initialized\n"); - /* initialise polyline for von Koch and simular domains */ + /* initialise polyline for von Koch and similar domains */ npolyline = init_polyline(MDEPTH, polyline); for (i=0; i= NY) jplus -= jmid; - jminus = j-1; - if (jminus < jmid) jminus += jmid; - } - } - else if (B_COND == BC_VPER_HABS) - { - iplus = (i+1); if (iplus == NX) iplus = NX-1; - iminus = (i-1); if (iminus == -1) iminus = 0; - if (j < jmid) /* lower half */ - { - jplus = (j+1); - if (jplus >= jmid) jplus -= jmid; - jminus = (j-1); - if (jminus < 0) jminus += jmid; - } - else /* upper half */ - { - jplus = j+1; - if (jplus >= NY) jplus -= jmid; - jminus = j-1; - if (jminus < jmid) jminus += jmid; - } - } - - /* imposing linear wave on top and bottom by making Laplacian 1d */ - if (OSCILLATE_TOPBOT) - { - if (j == NY-1) jminus = NY-1; - else if (j == 0) jplus = 0; - } - - delta = phi_in[iplus][j] + phi_in[iminus][j] + phi_in[i][jplus] + phi_in[i][jminus] - 4.0*phi_in[i][j]; - - x = phi_in[i][j]; - y = psi_in[i][j]; - - /* evolve phi */ - if ((B_COND == BC_PERIODIC)||(B_COND == BC_DIRICHLET)) - phi_out[i][j] = -y + 2*x + cc*delta - KAPPA*x - gamma*(x-y); - else if ((B_COND == BC_ABSORBING)||(B_COND == BC_ABS_REFLECT)) - { - if ((i>0)&&(i0)&&(j0)&&(i VMAX) phi_out[i][j] = VMAX; - if (phi_out[i][j] < -VMAX) phi_out[i][j] = -VMAX; - if (psi_out[i][j] > VMAX) psi_out[i][j] = VMAX; - if (psi_out[i][j] < -VMAX) psi_out[i][j] = -VMAX; - } - } - } - } -// printf("phi(0,0) = %.3lg, psi(0,0) = %.3lg\n", phi[NX/2][NY/2], psi[NX/2][NY/2]); -} - void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX], double *psi_out[NX], short int *xy_in[NX]) /* time step of field evolution */ @@ -717,6 +588,13 @@ void evolve_wave(double *phi[NX], double *psi[NX], double *phi_tmp[NX], double * } +void draw_color_bar(int plot, double range) +{ + if (ROTATE_COLOR_SCHEME) draw_color_scheme(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range); + else draw_color_scheme(1.7, YMIN + 0.1, 1.9, YMAX - 0.1, plot, -range, range); +} + + void animation() { @@ -739,13 +617,14 @@ void animation() printf("initializing circle configuration\n"); if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN_B == D_CIRCLES)) init_circle_config_comp(circles); if ((B_DOMAIN == D_POLYGONS)|(B_DOMAIN_B == D_POLYGONS)) init_polygon_config_comp(polygons); +// for (i=0; i