daviddoji/YouTube-simulations
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Add files via upload

Browse Source
This commit is contained in:
nilsberglund-orleans 2021-10-24 15:20:56 +02:00 committed by GitHub
parent 660e3d15fd
commit dadfb985ed
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
18 changed files with 3207 additions and 502 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1183
colormaps.c Normal file
View File

File diff suppressed because it is too large Load Diff

130
colors_waves.c Normal file
View File

@ -0,0 +1,130 @@
#include "colormaps.c"
double color_amplitude(double value, double scale, int time)
/* transforms the wave amplitude into a double in [-1,1] to feed into color scheme */
{
return(tanh(SLOPE*value/scale)*exp(-((double)time*ATTENUATION)));
}
double color_amplitude_asym(double value, double scale, int time)
/* transforms the wave amplitude into a double in [-1,1] to feed into color scheme */
{
return(2.0*tanh(SLOPE*value/scale)*exp(-((double)time*ATTENUATION)) - 1.0);
}
void color_scheme(int scheme, double value, double scale, int time, double rgb[3]) /* color scheme */
{
double hue, y, r, amplitude;
int intpart;
/* saturation = r, luminosity = y */
switch (scheme) {
case (C_LUM):
{
hue = COLORHUE + (double)time*COLORDRIFT/(double)NSTEPS;
if (hue < 0.0) hue += 360.0;
if (hue >= 360.0) hue -= 360.0;
r = 0.9;
amplitude = color_amplitude(value, scale, time);
y = LUMMEAN + amplitude*LUMAMP;
intpart = (int)y;
y -= (double)intpart;
hsl_to_rgb(hue, r, y, rgb);
break;
}
case (C_HUE):
{
r = 0.9;
amplitude = color_amplitude(value, scale, time);
y = 0.5;
hue = HUEMEAN + amplitude*HUEAMP;
if (hue < 0.0) hue += 360.0;
if (hue >= 360.0) hue -= 360.0;
hsl_to_rgb(hue, r, y, rgb);
break;
}
case (C_ONEDIM):
{
amplitude = color_amplitude(value, scale, time);
amp_to_rgb(0.5*(1.0 + amplitude), rgb);
break;
}
}
}
void color_scheme_lum(int scheme, double value, double scale, int time, double lum, double rgb[3]) /* color scheme */
{
double hue, y, r, amplitude;
int intpart;
/* saturation = r, luminosity = y */
switch (scheme) {
case (C_LUM):
{
hue = COLORHUE + (double)time*COLORDRIFT/(double)NSTEPS;
if (hue < 0.0) hue += 360.0;
if (hue >= 360.0) hue -= 360.0;
r = 0.9;
amplitude = color_amplitude(value, scale, time);
y = LUMMEAN + amplitude*LUMAMP;
intpart = (int)y;
y -= (double)intpart;
hsl_to_rgb(hue, r, y, rgb);
break;
}
case (C_HUE):
{
r = 0.9;
amplitude = color_amplitude(value, scale, time);
y = lum;
hue = HUEMEAN + amplitude*HUEAMP;
if (hue < 0.0) hue += 360.0;
if (hue >= 360.0) hue -= 360.0;
hsl_to_rgb(hue, r, y, rgb);
break;
}
}
}
void color_scheme_asym(int scheme, double value, double scale, int time, double rgb[3]) /* color scheme */
{
double hue, y, r, amplitude;
int intpart;
/* saturation = r, luminosity = y */
switch (scheme) {
case (C_LUM):
{
hue = COLORHUE + (double)time*COLORDRIFT/(double)NSTEPS;
if (hue < 0.0) hue += 360.0;
if (hue >= 360.0) hue -= 360.0;
r = 0.9;
amplitude = color_amplitude(value, scale, time);
y = LUMMEAN + amplitude*LUMAMP;
intpart = (int)y;
y -= (double)intpart;
hsl_to_rgb(hue, r, y, rgb);
break;
}
case (C_HUE):
{
r = 0.9;
amplitude = color_amplitude_asym(value, scale, time);
y = 0.5;
hue = HUEMEAN + 0.8*amplitude*HUEAMP;
if (hue < 0.0) hue += 360.0;
if (hue >= 360.0) hue -= 360.0;
hsl_to_rgb(hue, r, y, rgb);
break;
}
case (C_ONEDIM):
{
amplitude = color_amplitude(value, scale, time);
amp_to_rgb(amplitude, rgb);
break;
}
}
}

30
drop_billiard.c
View File

@ -30,14 +30,14 @@
#define MOVIE 0 /* set to 1 to generate movie */
// #define WINWIDTH 720 /* window width */
#define WINWIDTH 1280 /* window width */
#define WINWIDTH 720 /* window width */
// #define WINWIDTH 1280 /* window width */
#define WINHEIGHT 720 /* window height */
#define XMIN -2.0
#define XMAX 2.0 /* x interval */
// #define XMIN -1.125
// #define XMAX 1.125 /* x interval */
// #define XMIN -2.0
// #define XMAX 2.0 /* x interval */
#define XMIN -1.125
#define XMAX 1.125 /* x interval */
#define YMIN -1.125
#define YMAX 1.125 /* y interval for 9/16 aspect ratio */
@ -65,9 +65,9 @@
// #define LAMBDA 3.75738973 /* sin(36°)/sin(9°) for 5-star shape with 90° angles */
// #define LAMBDA -1.73205080756888 /* -sqrt(3) for Reuleaux triangle */
// #define LAMBDA 1.73205080756888 /* sqrt(3) for triangle tiling plane */
#define MU 0.9 /* second parameter controlling shape of billiard */
#define MU 1.0 /* 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 NPOLY 6 /* number of sides of polygon */
#define APOLY 0.0 /* 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 1 /* set to 1 to draw additional construction lines for billiard */
@ -78,10 +78,10 @@
#define NPART 10000 /* number of particles */
#define NPARTMAX 100000 /* maximal number of particles after resampling */
#define NSTEPS 4200 /* number of frames of movie */
#define TIME 100 /* time between movie frames, for fluidity of real-time simulation */
#define NSTEPS 5000 /* number of frames of movie */
#define TIME 150 /* time between movie frames, for fluidity of real-time simulation */
#define DPHI 0.0001 /* integration step */
#define NVID 50 /* number of iterations between images displayed on screen */
#define NVID 75 /* 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 */
@ -99,12 +99,14 @@
/* color and other graphical parameters */
#define NCOLORS 12 /* number of colors */
#define COLORSHIFT 0 /* hue of initial color */
#define COLOR_PALETTE 1 /* Color palette, see list in global_pdes.c */
#define NCOLORS 6 /* 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 100 /* number of segments of boundary */
#define BILLIARD_WIDTH 4 /* width of billiard */
#define FRONT_WIDTH 3 /* width of wave front */
#define FRONT_WIDTH 4 /* width of wave front */
#define BLACK 1 /* set to 1 for black background */
#define COLOR_OUTSIDE 0 /* set to 1 for colored outside */

16
global_particles.c
View File

@ -36,6 +36,7 @@ double x_shooter = -0.2, y_shooter = -0.6, x_target = 0.4, y_target = 0.7;
#define D_CIRCLES 20 /* several circles */
#define D_CIRCLES_IN_RECT 21 /* several circles inside a rectangle */
#define D_CIRCLES_IN_GENUSN 22 /* several circles in polygon with identified opposite sides */
#define D_CIRCLES_IN_TORUS 23 /* several circles in a rectangle with periodic boundary conditions */
#define C_FOUR_CIRCLES 0 /* four circles almost touching each other */
#define C_SQUARE 1 /* square grid of circles */
@ -49,4 +50,17 @@ double x_shooter = -0.2, y_shooter = -0.6, x_target = 0.4, y_target = 0.7;
#define C_LASER 11 /* laser fight in a room of mirrors */
#define C_LASER_GENUSN 12 /* laser fight in a translation surface */
#define C_LASER_GENUSN 12 /* laser fight in a translation surface */
/* 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 */

25
global_pdes.c
View File

@ -1,5 +1,8 @@
/* Global variables and parameters for wave_billiard, heat and schrodinger */
// #include "hsluv.c"
/* Basic math */
#define PI 3.141592654
@ -33,13 +36,17 @@
#define D_TWO_PARABOLAS 19 /* two facing parabolic antennas */
#define D_CIRCLES 20 /* several circles */
#define D_CIRCLES_IN_RECT 201 /* several circles in a rectangle */
#define D_FOUR_PARABOLAS 31 /* four parabolas with axes in NSEW directions */
#define D_POLY_PARABOLAS 32 /* polygon with parabolic sides */
#define D_PENROSE 33 /* Penrose illumination problem */
#define D_PENROSE 33 /* Penrose unilluminable room */
#define D_HYPERBOLA 34 /* one branch of hyperbola */
#define D_TOKARSKY 35 /* Tokarsky unilluminable room */
#define NMAXCIRCLES 1000 /* total number of circles (must be at least NCX*NCY for square grid) */
#define D_ISOSPECTRAL 37 /* isospectral billiards */
#define NMAXCIRCLES 600 /* total number of circles (must be at least NCX*NCY for square grid) */
#define C_SQUARE 0 /* square grid of circles */
#define C_HEX 1 /* hexagonal/triangular grid of circles */
@ -48,7 +55,7 @@
#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 */
@ -102,6 +109,18 @@
#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 */
double circlex[NMAXCIRCLES], circley[NMAXCIRCLES], circlerad[NMAXCIRCLES]; /* position and radius of circular scatterers */

14
heat.c
View File

@ -80,6 +80,17 @@
#define NGRIDX 15 /* 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
#define X_TARGET 0.4
#define Y_TARGET 0.7 /* shooter and target positions in laser fight */
#define ISO_XSHIFT_LEFT -1.65
#define ISO_XSHIFT_RIGHT 0.4
#define ISO_YSHIFT_LEFT -0.05
#define ISO_YSHIFT_RIGHT -0.05
#define ISO_SCALE 0.85 /* coordinates for isospectral billiards */
/* You can add more billiard tables by adapting the functions */
/* xy_in_billiard and draw_billiard in sub_wave.c */
@ -153,11 +164,10 @@
// #define HUEAMP -130.0 /* amplitude of variation of hue for color scheme C_HUE */
#include "hsluv.c"
#include "global_pdes.c"
#include "sub_wave.c"
double courant2; /* Courant parameter squared */
double dx2; /* spatial step size squared */
double intstep; /* integration step */

467
mangrove.c
View File

@ -49,8 +49,10 @@
#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 640 /* number of grid points on x axis */
#define NY 360 /* number of grid points on y axis */
// #define NX 1280 /* number of grid points on x axis */
// #define NY 720 /* number of grid points on y axis */
#define XMIN -2.0
#define XMAX 2.0 /* x interval */
@ -63,6 +65,7 @@
#define B_DOMAIN 20 /* choice of domain shape, see list in global_pdes.c */
// #define CIRCLE_PATTERN 1 /* pattern of circles, see list in global_pdes.c */
#define CIRCLE_PATTERN 8 /* pattern of circles, see list in global_pdes.c */
#define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */
@ -80,6 +83,17 @@
#define NGRIDX 15 /* 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
#define X_TARGET 0.4
#define Y_TARGET 0.7 /* shooter and target positions in laser fight */
#define ISO_XSHIFT_LEFT -1.65
#define ISO_XSHIFT_RIGHT 0.4
#define ISO_YSHIFT_LEFT -0.05
#define ISO_YSHIFT_RIGHT -0.05
#define ISO_SCALE 0.85 /* coordinates for isospectral billiards */
/* You can add more billiard tables by adapting the functions */
/* xy_in_billiard and draw_billiard below */
@ -90,18 +104,24 @@
#define OSCILLATE_TOPBOT 1 /* set to 1 to enforce a planar wave on top and bottom boundary */
#define X_SHIFT -0.9 /* x range on which to apply OSCILLATE_TOPBOT */
#define OMEGA 0.002 /* frequency of periodic excitation */
#define K_BC 3.0 /* spatial period of periodic excitation in y direction */
#define KX_BC 30.0 /* spatial period of periodic excitation in x direction */
#define KY_BC 10.0 /* spatial period of periodic excitation in y direction */
#define OMEGA 0.00133333333 /* frequency of periodic excitation */
#define K_BC 3.0 /* spatial period of periodic excitation in y direction */
#define KX_BC 10.0 /* spatial period of periodic excitation in x direction */
#define KY_BC 3.3333 /* spatial period of periodic excitation in y direction */
// #define KX_BC 20.0 /* spatial period of periodic excitation in x direction */
// #define KY_BC 6.66666 /* spatial period of periodic excitation in y direction */
// #define OMEGA 0.002 /* frequency of periodic excitation */
// #define K_BC 3.0 /* spatial period of periodic excitation in y direction */
// #define KX_BC 30.0 /* spatial period of periodic excitation in x direction */
// #define KY_BC 10.0 /* spatial period of periodic excitation in y direction */
#define AMPLITUDE 1.0 /* amplitude of periodic excitation */
#define COURANT 0.02 /* Courant number */
#define COURANTB 0.01 /* Courant number in medium B */
#define COURANTB 0.015 /* Courant number in medium B */
// #define COURANTB 0.00666 /* Courant number in medium B */
#define GAMMA 2.0e-6 /* damping factor in wave equation */
#define GAMMAB 2.5e-4 /* damping factor in wave equation */
// #define GAMMAB 5.0e-4 /* damping factor in wave equation */
// #define GAMMAB 1.0e-4 /* damping factor in wave equation */
#define GAMMA 3.0e-6 /* damping factor in wave equation */
#define GAMMAB 5.0e-4 /* damping factor in wave equation */
// #define GAMMA 2.0e-6 /* damping factor in wave equation */
// #define GAMMAB 2.5e-4 /* damping factor in wave equation */
#define GAMMA_SIDES 1.0e-4 /* damping factor on boundary */
#define GAMMA_TOPBOT 1.0e-6 /* damping factor on boundary */
#define KAPPA 0.0 /* "elasticity" term enforcing oscillations */
@ -119,8 +139,8 @@
/* Parameters for length and speed of simulation */
// #define NSTEPS 1000 /* number of frames of movie */
#define NSTEPS 4500 /* number of frames of movie */
#define NSTEPS 1000 /* number of frames of movie */
// #define NSTEPS 5500 /* number of frames of movie */
#define NVID 60 /* number of iterations between images displayed on screen */
#define NSEG 100 /* number of segments of boundary */
#define INITIAL_TIME 100 /* time after which to start saving frames */
@ -167,42 +187,40 @@
#define MANGROVE_HUE_MIN 180.0 /* color of original mangrove */
#define MANGROVE_HUE_MAX -50.0 /* color of saturated mangrove */
// #define MANGROVE_EMAX 5.0e-3 /* max energy for mangrove to survive */
#define MANGROVE_EMAX 1.1e-3 /* max energy for mangrove to survive */
#define MANGROVE_EMAX 1.5e-3 /* max energy for mangrove to survive */
#define RANDOM_RADIUS 1 /* set to 1 for random circle radius */
#define ERODE_MANGROVES 0 /* set to 1 for mangroves to be eroded */
#define ERODE_MANGROVES 1 /* set to 1 for mangroves to be eroded */
#define RECOVER_MANGROVES 1 /* set to 1 to allow mangroves to recover */
#define MOVE_MANGROVES 1 /* set to 1 for mobile mangroves */
#define DETACH_MANGROVES 1 /* set to 1 for mangroves to be able to detach */
#define INERTIA 1 /* set to 1 for taking inertia into account */
#define REPELL_MANGROVES 1 /* set to 1 for mangroves to repell each other */
#define DT_MANGROVE 0.1 /* time step for mangrove displacement */
#define KSPRING 0.25 /* spring constant of mangroves */
#define KWAVE 2.0 /* constant in force due to wave gradient */
#define KSPRING 0.05 /* spring constant of mangroves */
#define KWAVE 4.0 /* constant in force due to wave gradient */
#define KREPEL 5.0 /* constant in repelling force between mangroves */
#define REPEL_RADIUS 1.1 /* radius in which repelling force acts (in units of mangrove radius) */
#define DXMAX 0.02 /* max displacement of mangrove in one time step */
#define L_DETACH 0.2 /* spring length beyond which mangroves detach */
#define DAMP_MANGROVE 0.2 /* damping coefficient of mangroves */
#define L_DETACH 0.25 /* spring length beyond which mangroves detach */
#define DAMP_MANGROVE 0.1 /* damping coefficient of mangroves */
#define MANGROVE_MASS 1.5 /* mass of mangrove of radius MU */
#define HASHX 25 /* size of hashgrid in x direction */
#define HASHY 15 /* size of hashgrid in y direction */
#define HASHMAX 10 /* maximal number of mangroves per hashgrid cell */
#define HASHGRID_PADDING 0.1 /* padding of hashgrid outside simulation window */
/* 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 "hsluv.c"
#include "global_pdes.c"
#include "sub_wave.c"
#include "wave_common.c"
double courant2, courantb2; /* Courant parameters squared */
double circle_energy[NMAXCIRCLES]; /* energy dissipated by the circles */
double circley_wrapped[NMAXCIRCLES]; /* position of circle centers wrapped vertically */
double anchor_x[NMAXCIRCLES]; /* points moving circles are attached to */
double anchor_y[NMAXCIRCLES]; /* points moving circles are attached to */
double vx[NMAXCIRCLES]; /* x velocity of circles */
double vy[NMAXCIRCLES]; /* y velocity of circles */
double circlerad_initial[NMAXCIRCLES]; /* initial circle radii */
double mass_inverse[NMAXCIRCLES]; /* inverse of mangrove mass */
short int circle_attached[NMAXCIRCLES]; /* has value 1 if the circle is attached to its anchor */
/*********************/
/* animation part */
@ -411,15 +429,109 @@ void evolve_wave(double *phi[NX], double *psi[NX], double *phi_tmp[NX], double *
evolve_wave_half(phi_tmp, psi_tmp, phi, psi, xy_in);
}
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]);
}
void compute_repelling_force(int i, int j, double force[2])
/* compute repelling force of mangrove j on mangrove i */
{
double x1, y1, x2, y2, distance, r, f;
x1 = circlex[i];
y1 = circley[i];
x2 = circlex[j];
y2 = circley[j];
distance = module2(x2 - x1, y2 - y1);
r = circlerad[i] + circlerad[j];
if (r <= 0.0) r = 0.001*MU;
f = KREPEL/(0.001 + distance*distance);
if ((distance > 0.0)&&(distance < REPEL_RADIUS*r))
{
force[0] = f*(x1 - x2)/distance;
force[1] = f*(y1 - y2)/distance;
}
else
{
force[0] = 0.0;
force[1] = 0.0;
}
}
void update_hashgrid(int* mangrove_hashx, int* mangrove_hashy, int* hashgrid_number, int* hashgrid_mangroves)
{
int i, j, k, n, m, ij[2], max = 0;
printf("Updating hashgrid_number\n");
for (i=0; i<HASHX*HASHY; i++) hashgrid_number[i] = 0;
printf("Updated hashgrid_number\n");
/* place each mangrove in hash grid */
for (k=1; k<ncircles; k++)
// if (circleactive[k])
{
// printf("placing circle %i\t", k);
hash_xy_to_ij(circlex[k], circley[k], ij);
i = ij[0]; j = ij[1];
// printf("ij = (%i, %i)\t", i, j);
n = hashgrid_number[i*HASHY + j];
m = i*HASHY*HASHMAX + j*HASHMAX + n;
// printf("n = %i, m = %i\n", n, m);
if (m < HASHX*HASHY*HASHMAX) hashgrid_mangroves[m] = k;
else printf("Too many mangroves in hash cell, try increasing HASHMAX\n");
hashgrid_number[i*HASHY + j]++;
mangrove_hashx[k] = i;
mangrove_hashy[k] = j;
if (n > max) max = n;
// printf("Placed mangrove %i at (%i,%i) in hashgrid\n", k, ij[0], ij[1]);
// printf("%i mangroves at (%i,%i)\n", hashgrid_number[ij[0]][ij[1]], ij[0], ij[1]);
}
printf("Maximal number of mangroves per hash cell: %i\n", max);
}
void animation()
{
double time, scale, diss, rgb[3], hue, y, dissip, ej, gradient[2], dx, dy, dt, xleft, xright, length;
double time, scale, diss, rgb[3], hue, y, dissip, ej, gradient[2], dx, dy, dt, xleft, xright,
length, fx, fy, force[2];
double *phi[NX], *psi[NX], *phi_tmp[NX], *psi_tmp[NX];
short int *xy_in[NX], redraw = 0;
int i, j, s, ij[2];
int i, j, k, n, s, ij[2], i0, iplus, iminus, j0, jplus, jminus, p, q;
static int imin, imax;
static short int first = 1;
double *circleenergy, *circley_wrapped, *anchor_x, *anchor_y, *vx, *vy, *circlerad_initial, *mass_inverse;
short int *circle_attached;
int *mangrove_hashx, *mangrove_hashy, *hashgrid_number, *hashgrid_mangroves;
/* Since NX and NY are big, it seemed wiser to use some memory allocation here */
for (i=0; i<NX; i++)
@ -431,6 +543,21 @@ void animation()
xy_in[i] = (short int *)malloc(NY*sizeof(short int));
}
circleenergy = (double *)malloc(NMAXCIRCLES*sizeof(double)); /* energy dissipated by the circles */
circley_wrapped = (double *)malloc(NMAXCIRCLES*sizeof(double)); /* position of circle centers wrapped vertically */
anchor_x = (double *)malloc(NMAXCIRCLES*sizeof(double)); /* points moving circles are attached to */
anchor_y = (double *)malloc(NMAXCIRCLES*sizeof(double)); /* points moving circles are attached to */
vx = (double *)malloc(NMAXCIRCLES*sizeof(double)); /* x velocity of circles */
vy = (double *)malloc(NMAXCIRCLES*sizeof(double)); /* y velocity of circles */
circlerad_initial = (double *)malloc(NMAXCIRCLES*sizeof(double)); /* initial circle radii */
mass_inverse = (double *)malloc(NMAXCIRCLES*sizeof(double)); /* inverse of mangrove mass */
circle_attached = (short int *)malloc(NMAXCIRCLES*sizeof(short int)); /* has value 1 if the circle is attached to its anchor */
mangrove_hashx = (int *)malloc(NMAXCIRCLES*sizeof(int)); /* hash grid positions of mangroves */
mangrove_hashy = (int *)malloc(NMAXCIRCLES*sizeof(int)); /* hash grid positions of mangroves */
hashgrid_number = (int *)malloc(HASHX*HASHY*sizeof(int)); /* total number of mangroves in each hash grid cell */
hashgrid_mangroves = (int *)malloc(HASHX*HASHY*HASHMAX*sizeof(int)); /* numbers of mangoves in each hash grid cell */
/* initialise positions and radii of circles */
if (B_DOMAIN == D_CIRCLES) init_circle_config();
@ -455,13 +582,14 @@ void animation()
/* initialise mangroves */
for (i=0; i < ncircles; i++)
{
circle_energy[i] = 0.0;
circleenergy[i] = 0.0;
y = circley[i];
if (y >= YMAX) y -= circlerad[i];
if (y <= YMIN) y += circlerad[i];
// if (y >= YMAX) y -= (YMAX - YMIN);
// if (y <= YMIN) y += (YMAX - YMIN);
circley_wrapped[i] = y;
// circleactive[i] = 1;
if (RANDOM_RADIUS) circlerad[i] = circlerad[i]*(0.75 + 0.5*((double)rand()/RAND_MAX));
@ -483,6 +611,9 @@ void animation()
}
}
/* initialise hash table for interacting mangroves */
if (REPELL_MANGROVES) update_hashgrid(mangrove_hashx, mangrove_hashy, hashgrid_number, hashgrid_mangroves);
if (first) /* compute box limits where circles are reset */
{
/* find leftmost and rightmost circle */
@ -515,7 +646,7 @@ void animation()
for (i=0; i<=INITIAL_TIME + NSTEPS; i++)
{
//printf("%d\n",i);
printf("Computing frame %d\n",i);
/* compute the variance of the field to adjust color scheme */
/* the color depends on the field divided by sqrt(1 + variance) */
if (SCALE)
@ -525,83 +656,24 @@ void animation()
}
else scale = 1.0;
printf("Drawing wave\n");
draw_wave(phi, psi, xy_in, scale, i, PLOT);
printf("Evolving wave\n");
for (j=0; j<NVID; j++)
{
// printf("%i ", j);
evolve_wave(phi, psi, phi_tmp, psi_tmp, xy_in);
// if (i % 10 == 9) oscillate_linear_wave(0.2*scale, 0.15*(double)(i*NVID + j), -1.5, YMIN, -1.5, YMAX, phi, psi);
}
/* compute energy dissipated in obstacles */
if (ERODE_MANGROVES) for (j=0; j<ncircles; j++)
{
dissip = compute_dissipation(phi, psi, xy_in, circlex[j], circley_wrapped[j]);
/* make sure the dissipation does not grow too fast because of round-off/blow-up */
if (dissip > 0.1*MANGROVE_EMAX)
{
dissip = 0.1*MANGROVE_EMAX;
printf("Flooring dissipation!\n");
}
if (circleactive[j])
{
circle_energy[j] += dissip;
ej = circle_energy[j];
if (ej <= MANGROVE_EMAX)
{
if (ej > 0.0)
{
hue = MANGROVE_HUE_MIN + (MANGROVE_HUE_MAX - MANGROVE_HUE_MIN)*ej/MANGROVE_EMAX;
if (hue < 0.0) hue += 360.0;
}
else hue = MANGROVE_HUE_MIN;
hsl_to_rgb(hue, 0.9, 0.5, rgb);
if (j%NGRIDY == 0) printf("Circle %i, energy %.5lg, hue %.5lg\n", j, ej, hue);
draw_colored_circle(circlex[j], circley[j], circlerad[j], NSEG, rgb);
/* shrink mangrove */
if (ej > 0.0)
{
// circlerad[j] -= MU*ej*ej/(MANGROVE_EMAX*MANGROVE_EMAX);
// if (circlerad[j] < 0.0) circlerad[j] = 0.0;
circlerad[j] = circlerad_initial[j]*(1.0 - ej*ej/(MANGROVE_EMAX*MANGROVE_EMAX));
redraw = 1;
}
else circlerad[j] = circlerad_initial[j];
}
else /* remove mangrove */
{
circleactive[j] = 0;
/* reinitialize table xy_in */
redraw = 1;
}
}
else /* allow disabled mangroves to recover */
{
circle_energy[j] -= 0.15*dissip;
// circlerad[j] += 0.005*MU;
// if (circlerad[j] > MU) circlerad[j] = MU;
// if ((circle_energy[j] < 0.0)&&(circlerad[j] > 0.0))
if (circle_energy[j] < 0.0)
{
circleactive[j] = 1;
// circlerad[j] = circlerad[j]*(0.75 + 0.5*((double)rand()/RAND_MAX));
circlerad[j] = circlerad_initial[j];
circle_energy[j] = -MANGROVE_EMAX;
/* reinitialize table xy_in */
redraw = 1;
}
}
// printf("Circle %i, energy %.5lg\n", j, circle_energy[j]);
}
/* move mangroves */
if (MOVE_MANGROVES) for (j=0; j<ncircles; j++) if (circleactive[j])
{
compute_gradient(phi, psi, circlex[j], circley_wrapped[j], gradient);
// printf("gradient = (%.3lg, %.3lg)\t", gradient[0], gradient[1]);
// if (j%NGRIDY == 0) printf("gradient (%.3lg, %.3lg)\n", gradient[0], gradient[1]);
// if (j%NGRIDY == 0) printf("circle %i (%.3lg, %.3lg) -> ", j, circlex[j], circley[j]);
@ -617,11 +689,41 @@ void animation()
dy += DT_MANGROVE*(-KSPRING*(circley_wrapped[j] - anchor_y[j]));
}
/* compute repelling force from other mangroves */
if (REPELL_MANGROVES)
{
/* determine neighboring grid points */
i0 = mangrove_hashx[j];
iminus = i0 - 1; if (iminus < 0) iminus = 0;
iplus = i0 + 1; if (iplus >= HASHX) iplus = HASHX-1;
j0 = mangrove_hashy[j];
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<hashgrid_number[p*HASHY+q]; k++)
if (circleactive[hashgrid_mangroves[p*HASHY*HASHMAX + q*HASHMAX + k]])
{
compute_repelling_force(j, hashgrid_mangroves[p*HASHY*HASHMAX + q*HASHMAX + k], force);
fx += force[0];
fy += force[1];
}
// if (fx*fx + fy*fy > 0.001) printf("Force on mangrove %i: (%.3f, %.3f)\n", j, fx, fy);
dx += DT_MANGROVE*fx;
dy += DT_MANGROVE*fy;
}
/* detach mangrove if spring is too long */
if (DETACH_MANGROVES)
{
length = module2(circlex[j] - anchor_x[j], circley_wrapped[j] - anchor_y[j]);
if (j%NGRIDY == 0) printf("spring length %.i: %.3lg\n", j, length);
// if (j%NGRIDY == 0) printf("spring length %.i: %.3lg\n", j, length);
// if (length > L_DETACH) circle_attached[j] = 0;
if (length*mass_inverse[j] > L_DETACH) circle_attached[j] = 0;
}
@ -638,9 +740,9 @@ void animation()
circlex[j] += vx[j]*DT_MANGROVE;
circley[j] += vy[j]*DT_MANGROVE;
circley_wrapped[j] += vy[j]*DT_MANGROVE;
if (j%NGRIDY == 0)
printf("circle %.i: (dx,dy) = (%.3lg,%.3lg), (vx,vy) = (%.3lg,%.3lg)\n",
j, circlex[j]-anchor_x[j], circley[j]-anchor_y[j], vx[j], vy[j]);
// if (j%NGRIDY == 0)
// printf("circle %.i: (dx,dy) = (%.3lg,%.3lg), (vx,vy) = (%.3lg,%.3lg)\n",
// j, circlex[j]-anchor_x[j], circley[j]-anchor_y[j], vx[j], vy[j]);
}
else
{
@ -655,10 +757,148 @@ void animation()
if (circley_wrapped[j] >= YMAX) circley_wrapped[j] = YMAX;
// if (j%NGRIDY == 0) printf("(%.3lg, %.3lg)\n", circlex[j], circley[j]);
redraw = 1;
}
/* test for debugging */
if (1) for (j=0; j<ncircles; j++)
{
dissip = compute_dissipation(phi, psi, xy_in, circlex[j], circley_wrapped[j]);
/* make sure the dissipation does not grow too fast because of round-off/blow-up */
if (dissip > 0.1*MANGROVE_EMAX)
{
dissip = 0.1*MANGROVE_EMAX;
printf("Flooring dissipation!\n");
}
if (circleactive[j])
{
circleenergy[j] += dissip;
ej = circleenergy[j];
// printf("ej = %.3f\n", ej);
if (ej <= MANGROVE_EMAX)
{
if (ej > 0.0)
{
hue = MANGROVE_HUE_MIN + (MANGROVE_HUE_MAX - MANGROVE_HUE_MIN)*ej/MANGROVE_EMAX;
if (hue < 0.0) hue += 360.0;
}
else hue = MANGROVE_HUE_MIN;
hsl_to_rgb(hue, 0.9, 0.5, rgb);
// if (j%NGRIDY == 0) printf("Circle %i, energy %.5lg, hue %.5lg\n", j, ej, hue);
draw_colored_circle(circlex[j], circley[j], circlerad[j], NSEG, rgb);
/* shrink mangrove */
if ((ERODE_MANGROVES)&&(ej > 0.0))
{
circlerad[j] = circlerad_initial[j]*(1.0 - ej*ej/(MANGROVE_EMAX*MANGROVE_EMAX));
redraw = 1;
}
else circlerad[j] = circlerad_initial[j];
}
else /* remove mangrove */
{
circleactive[j] = 0;
/* reinitialize table xy_in */
redraw = 1;
}
}
else if (RECOVER_MANGROVES) /* allow disabled mangroves to recover */
{
circleenergy[j] -= 0.15*dissip;
printf("Circle %i energy %.3lg\n", j, circleenergy[j]);
if (circleenergy[j] < 0.0)
{
printf("Reactivating circle %i?\n", j);
/* THE PROBLEM occurs when circleactive[0] is set to 1 again */
if (j>0) circleactive[j] = 1;
circlerad[j] = circlerad_initial[j];
circleenergy[j] = -MANGROVE_EMAX;
/* reinitialize table xy_in */
redraw = 1;
}
}
}
/* compute energy dissipated in obstacles */
/* if (ERODE_MANGROVES) for (j=0; j<ncircles; j++)
{
// printf("j = %i\t", j);
dissip = compute_dissipation(phi, psi, xy_in, circlex[j], circley_wrapped[j]);
printf("dissip = %.3f\t", dissip);
/* make sure the dissipation does not grow too fast because of round-off/blow-up */
// if (dissip > 0.1*MANGROVE_EMAX)
// {
// dissip = 0.1*MANGROVE_EMAX;
// printf("Flooring dissipation!\n");
// }
//
// if (circleactive[j])
// {
// circleenergy[j] += dissip;
// ej = circleenergy[j];
// printf("ej = %.3f\n", ej);
// if (ej <= MANGROVE_EMAX)
// {
// if (ej > 0.0)
// {
// hue = MANGROVE_HUE_MIN + (MANGROVE_HUE_MAX - MANGROVE_HUE_MIN)*ej/MANGROVE_EMAX;
// if (hue < 0.0) hue += 360.0;
// }
// else hue = MANGROVE_HUE_MIN;
// hsl_to_rgb(hue, 0.9, 0.5, rgb);
// if (j%NGRIDY == 0) printf("Circle %i, energy %.5lg, hue %.5lg\n", j, ej, hue);
// draw_colored_circle(circlex[j], circley[j], circlerad[j], NSEG, rgb);
//
// /* shrink mangrove */
// if (ej > 0.0)
// {
// circlerad[j] -= MU*ej*ej/(MANGROVE_EMAX*MANGROVE_EMAX);
// if (circlerad[j] < 0.0) circlerad[j] = 0.0;
// circlerad[j] = circlerad_initial[j]*(1.0 - ej*ej/(MANGROVE_EMAX*MANGROVE_EMAX));
// redraw = 1;
// }
// else circlerad[j] = circlerad_initial[j];
// }
// else /* remove mangrove */
// {
// circleactive[j] = 0;
/* reinitialize table xy_in */
// redraw = 1;
// }
// }
// else /* allow disabled mangroves to recover */
// {
// circleenergy[j] -= 0.15*dissip;
// printf("ej = %.3f\n", circleenergy[j]);
// circlerad[j] += 0.005*MU;
// if (circlerad[j] > MU) circlerad[j] = MU;
// if ((circleenergy[j] < 0.0)&&(circlerad[j] > 0.0))
// if (circleenergy[j] < 0.0)
// {
// circleactive[j] = 1;
// circlerad[j] = circlerad[j]*(0.75 + 0.5*((double)rand()/RAND_MAX));
// circlerad[j] = circlerad_initial[j];
// circleenergy[j] = -MANGROVE_EMAX;
/* reinitialize table xy_in */
// redraw = 1;
// }
// }
// printf("Circle %i, energy %.5lg\n", j, circleenergy[j]);
// }
printf("Updating hashgrid\n");
if (REPELL_MANGROVES) update_hashgrid(mangrove_hashx, mangrove_hashy, hashgrid_number, hashgrid_mangroves);
printf("Drawing billiard\n");
draw_billiard();
glutSwapBuffers();
@ -666,7 +906,7 @@ void animation()
if (redraw)
{
printf("Reinitializing xy_in\n");
init_xyin_xrange(xy_in, imin, NX-1);
init_xyin_xrange(xy_in, imin, NX);
// init_xyin_xrange(xy_in, imin, imax);
}
redraw = 0;
@ -701,6 +941,21 @@ void animation()
free(psi_tmp[i]);
free(xy_in[i]);
}
free(circleenergy);
free(circley_wrapped);
free(anchor_x);
free(anchor_y);
free(vx);
free(vy);
free(circlerad_initial);
free(mass_inverse);
free(circle_attached);
free(mangrove_hashx);
free(mangrove_hashy);
free(hashgrid_number);
free(hashgrid_mangroves);
}

145
particle_billiard.c
View File

@ -32,13 +32,18 @@
#define MOVIE 0 /* set to 1 to generate movie */
#define WINWIDTH 1280 /* window width */
// #define WINWIDTH 1280 /* window width */
#define WINWIDTH 720 /* 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 XMIN -1.4
#define XMAX 1.4 /* x interval */
#define YMIN -1.4
#define YMAX 1.4 /* 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 SCALING_FACTOR 1.0 /* scaling factor of drawing, needed for flower billiards, otherwise set to 1.0 */
@ -60,7 +65,7 @@
// #define LAMBDA 1.4 /* parameter controlling shape of domain */
// #define MU 0.2 /* second parameter controlling shape of billiard */
#define LAMBDA 1.5 /* parameter controlling shape of domain */
#define LAMBDA 1.3 /* parameter controlling shape of domain */
#define MU 0.3 /* 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 */
@ -74,15 +79,16 @@
/* Simulation parameters */
#define NPART 10000 /* number of particles */
#define NPART 100 /* 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 SHOWTRAILS 1 /* set to 1 to keep trails of the particles */
#define TEST_ACTIVE 1 /* set to 1 to test whether particle is in billiard */
#define NSTEPS 3500 /* number of frames of movie */
#define TIME 1200 /* time between movie frames, for fluidity of real-time simulation */
#define NSTEPS 1300 /* number of frames of movie */
#define TIME 2500 /* time between movie frames, for fluidity of real-time simulation */
#define DPHI 0.00001 /* integration step */
#define NVID 150 /* number of iterations between images displayed on screen */
@ -94,12 +100,14 @@
/* Colors and other graphical parameters */
#define NCOLORS 16 /* number of colors */
#define COLOR_PALETTE 0 /* Color palette, see list in global_pdes.c */
#define NCOLORS 64 /* number of colors */
#define COLORSHIFT 0 /* hue of initial color */
#define RAINBOW_COLOR 0 /* set to 1 to use different colors for all particles */
#define RAINBOW_COLOR 1 /* set to 1 to use different colors for all particles */
#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.01 /* length of velocity vectors */
#define LENGTH 0.03 /* length of velocity vectors */
#define BILLIARD_WIDTH 2 /* width of billiard */
#define PARTICLE_WIDTH 2 /* width of particles */
#define FRONT_WIDTH 3 /* width of wave front */
@ -170,6 +178,31 @@ void init_drop_config(double x0, double y0, double angle1, double angle2, double
}
}
void init_partial_drop_config(double x0, double y0, double angle1, double angle2, int particle1, int particle2,
int col, double *configs[NPARTMAX], int color[NPARTMAX], int newcolor[NPARTMAX])
/* initialize configuration: drop at (x0,y0) for a range of particles */
{
int i;
double dalpha, alpha;
double conf[2], pos[2];
while (angle2 < angle1) angle2 += DPI;
if (particle2 - particle1 > 1) dalpha = (angle2 - angle1)/((double)(particle2 - particle1-1));
else dalpha = 0.0;
for (i=particle1; i<particle2; i++)
{
alpha = angle1 + dalpha*((double)i);
// printf("alpha=%.5lg\n", alpha);
pos[0] = x0;
pos[1] = y0;
vbilliard_xy(configs[i], alpha, pos);
color[i] = col;
newcolor[i] = col;
}
}
void init_sym_drop_config(double x0, double y0, double angle1, double angle2, double *configs[NPARTMAX])
/* initialize configuration with two symmetric partial drops */
{
@ -217,6 +250,71 @@ void init_line_config(double x0, double y0, double x1, double y1, double angle,
}
void draw_config_showtrails(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], len;
glutSwapBuffers();
if (PAINT_INT) paint_billiard_interior();
glLineWidth(PARTICLE_WIDTH);
glEnable(GL_LINE_SMOOTH);
for (i=0; i<nparticles; i++)
{
// if (configs[i][2]<0.0)
// {
// vbilliard(configs[i]);
// if (!RAINBOW_COLOR)
// {
// color[i]++;
// if (color[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 */
{
@ -332,7 +430,8 @@ void graph_movie(int time, int color[NPARTMAX], double *configs[NPARTMAX], int a
{
// printf("reflecting particle %i\n", i);
c = vbilliard(configs[i]);
if (c>=0) color[i]++;
// if (c>=0) color[i]++;
if ((!RAINBOW_COLOR)&&(c>=0)) color[i]++;
if (!RAINBOW_COLOR)
{
color[i]++;
@ -354,7 +453,7 @@ void graph_movie(int time, int color[NPARTMAX], double *configs[NPARTMAX], int a
void animation()
{
double time, dt, alpha, r;
double time, dt, alpha, r, rgb[3];
double *configs[NPARTMAX];
int i, j, resamp = 1, s, i1, i2;
int *color, *newcolor, *active;
@ -430,15 +529,27 @@ void animation()
}
sleep(SLEEP1);
/* initialize drops in different colors */
init_partial_drop_config(0.0, 0.0, 0.0, DPI, 0, 2*NPART/5, 0, configs, color, newcolor);
init_partial_drop_config(0.0, 0.8, 0.0, DPI, 2*NPART/5, 4*NPART/5, 10, configs, color, newcolor);
init_partial_drop_config(1.2, 0.1, 0.0, DPI, 4*NPART/5, NPART, 36, configs, color, newcolor);
for (i=0; i<=NSTEPS; i++)
{
graph_movie(TIME, newcolor, configs, active);
draw_config(newcolor, configs, active);
if (SHOWTRAILS) draw_config_showtrails(newcolor, configs, active);
else draw_config(newcolor, configs, active);
// draw_config(newcolor, configs, active);
if (DRAW_BILLIARD) draw_billiard();
for (j=0; j<NPARTMAX; j++) color[j] = newcolor[j];
/* draw initial points */
draw_initial_condition_circle(0.0, 0.0, 0.02, 0);
draw_initial_condition_circle(0.0, 0.8, 0.02, 10);
draw_initial_condition_circle(1.2, 0.1, 0.02, 36);
if (MOVIE)
{

185
particle_pinball.c
View File

@ -38,6 +38,7 @@
#define XMIN -4.0
#define XMAX 4.0 /* x interval */
#define YMIN -1.25
// #define YMID 1.25
#define YMAX 3.25 /* y interval for 9/16 aspect ratio */
// #define XMIN -2.0
// #define XMAX 2.0 /* x interval */
@ -51,20 +52,30 @@
/* Choice of the billiard table, see global_particles.c */
#define B_DOMAIN 21 /* choice of domain shape */
#define B_DOMAIN 23 /* choice of domain shape */
#define CIRCLE_PATTERN 6 /* pattern of circles */
// #define CIRCLE_PATTERN 21 /* 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 NMAXCIRCLES 1000 /* total number of circles (must be at least NCX*NCY for square grid) */
#define NCX 44 /* 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 NCX 50 /* number of circles in x direction */
// #define NCY 20 /* number of circles in y direction */
// #define NCX 52 /* number of circles in x direction */
// #define NCY 30 /* number of circles in y direction */
#define NPOISSON 350 /* 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 LAMBDA 3.8105 /* parameter controlling shape of domain */
// #define MU 0.1 /* second parameter controlling shape of billiard */
// #define MU 0.085 /* second parameter controlling shape of billiard */
// #define MU 0.09 /* second parameter controlling shape of billiard */
#define MU 0.05 /* second parameter controlling shape of billiard */
// #define MU 0.034 /* 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 */
@ -85,9 +96,14 @@
#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 NSTEPS 11700 /* number of frames of movie */
// #define NSTEPS 1000 /* number of frames of movie */
// #define TIME 1500 /* time between movie frames, for fluidity of real-time simulation */
#define TIME 4000 /* time between movie frames, for fluidity of real-time simulation */
// #define DPHI 0.0001 /* integration step */
// #define DPHI 0.00002 /* integration step */
#define DPHI 0.00007 /* integration step */
// #define DPHI 0.000035 /* integration step */
#define NVID 150 /* number of iterations between images displayed on screen */
/* Decreasing TIME accelerates the animation and the movie */
@ -98,8 +114,11 @@
/* Colors and other graphical parameters */
#define NCOLORS 32 /* number of colors */
#define COLORSHIFT 0 /* hue of initial color */
#define COLOR_PALETTE 1 /* Color palette, see list in global_pdes.c */
// #define NCOLORS 256 /* number of colors */
#define NCOLORS 48 /* number of colors */
#define COLORSHIFT 2 /* 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) */
@ -108,6 +127,7 @@
#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 0 /* hue for single color */
#define COLOR_TRAJECTORY 8 /* hue for single color */
#define BLACK 1 /* set to 1 for black background */
@ -118,7 +138,7 @@
#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 PAUSE 500 /* 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 */
@ -133,6 +153,7 @@
int ncol = 0, nobst = 0, nmaxpeg = 0;
int npath[NPATHBINS];
double max_free_path = 0.0;
/*********************/
/* animation part */
@ -404,6 +425,7 @@ void graph_movie(int time, int color[NPARTMAX], double *configs[NPARTMAX], int a
{
int i, j, k, c;
double rgb[3];
static double total_pathlength = 0.0;
for (j=0; j<time; j++)
{
@ -424,22 +446,44 @@ void graph_movie(int time, int color[NPARTMAX], double *configs[NPARTMAX], int a
glEnd ();
}
total_pathlength += configs[i][3];
c = vbilliard(configs[i]);
ncol++;
/* update number of collisions, not counting boundary for periodic b.c. */
if (B_DOMAIN != D_CIRCLES_IN_TORUS) ncol++;
else if (c >= 0) ncol++;
if ((c >= 0)&&(circlecolor[c] == 0)) nobst++;
circlecolor[c]++;
/* take care of circles doubled because of periodic boundary conditions */
if ((circleactive[c])&&(B_DOMAIN == D_CIRCLES_IN_TORUS)&&(partner_circle[c] != c)) circlecolor[partner_circle[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 (ncol > 1) /* disregard very first collision */
{
if (B_DOMAIN != D_CIRCLES_IN_TORUS)
{
k = (int)((double)NPATHBINS*configs[i][3]/PATHLMAX);
if (k < NPATHBINS) npath[k]++;
if (total_pathlength > max_free_path) max_free_path = total_pathlength;
total_pathlength = 0.0;
}
else /* case with periodic boundary conditions */
{
if (c >= 0) /* a circle is hit, update histogram */
{
// printf("total path length %.3lg\n", total_pathlength);
k = (int)((double)NPATHBINS*total_pathlength/PATHLMAX);
if (k < NPATHBINS) npath[k]++;
if (total_pathlength > max_free_path) max_free_path = total_pathlength;
total_pathlength = 0.0;
}
}
}
if (!RAINBOW_COLOR)
{
@ -514,8 +558,10 @@ void print_particle_numbers(double *configs[NPARTMAX])
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;
int i, n, colmax = 55, pegcollisions[90], nypegs = 70, meanpegs = 0, meansquarepegs = 0, total_coll = 0, ymax = 0,
meanbins = 0, meansquarebins, total_bin = 0, stephits = 10, tickstephits = 1;
double x, y, yscale = 110.0, y0, dx, rgb[3], xshift, coll_mean, coll_stdv, path_mean, path_stdv, ebin, len_over_bin,
steppath = 0.5;
char message[50];
glLineWidth(1);
@ -525,10 +571,10 @@ void draw_statistics()
xshift = XMIN + 0.3;
rgb[0] = 0.0; rgb[1] = 0.0; rgb[2] = 1.0;
/* histogramm of number of collisions per peg */
/* histogram of number of collisions per peg */
for (i=0; i<colmax; i++) pegcollisions[i] = 0;
for (i=0; i<ncircles; i++)
for (i=0; i<ncircles; i++) if ((circleactive[i])&&(!double_circle[i]))
{
n = circlecolor[i];
if (n < colmax) pegcollisions[n]++;
@ -537,6 +583,7 @@ void draw_statistics()
{
total_coll += pegcollisions[i];
meanpegs += i*pegcollisions[i];
meansquarepegs += i*i*pegcollisions[i];
}
for (i=1; i<colmax; i++)
@ -550,8 +597,8 @@ void draw_statistics()
glColor3f(1.0, 1.0, 1.0);
glBegin(GL_LINE_STRIP);
/* histogramm */
for (i=1; i<colmax; i++)
/* histogram */
for (i=0; i<colmax; i++)
{
x = xshift + (double)i*dx;
y = y0 + (double)pegcollisions[i]*YMAX/yscale;
@ -574,7 +621,7 @@ void draw_statistics()
glEnd ();
}
for (i=1; i<colmax; i++)
for (i=tickstephits; i<colmax; i+=tickstephits)
{
glBegin(GL_LINE_STRIP);
glVertex2d(xshift + (double)i*dx, y0 - 0.025);
@ -582,10 +629,10 @@ void draw_statistics()
glEnd ();
}
for (i=10; i<colmax; i+=10)
for (i=stephits; i<nypegs - stephits; i+=stephits)
{
sprintf(message, "%4d", i);
write_text_fixedwidth(xshift + (double)i*dx - 0.12, y0 - 0.12, message);
write_text_fixedwidth(xshift + (double)i*dx - 0.15, y0 - 0.12, message);
}
for (i=10; i<nypegs; i+=10)
@ -594,20 +641,34 @@ void draw_statistics()
write_text_fixedwidth(xshift - 0.3, y0 - 0.025 + (double)i*YMAX/yscale, message);
}
coll_mean = (double)meanpegs/(double)total_coll;
coll_stdv = sqrt((double)meansquarepegs/(double)total_coll - coll_mean*coll_mean);
sprintf(message, "hits");
write_text_fixedwidth(xshift + (double)(colmax-3)*dx, y0 - 0.12, message);
sprintf(message, "pegs");
write_text_fixedwidth(xshift - 0.25, y0 - 0.025 + (double)(nypegs - 3)*YMAX/yscale, message);
sprintf(message, "Mean %.4lg", (double)meanpegs/(double)total_coll);
write_text(-1.0, YMAX - 0.3, message);
// sprintf(message, "Mean %.4lg", (double)meanpegs/(double)total_coll);
/* histogramm of path lengths */
sprintf(message, "Max hits/peg %d", nmaxpeg);
write_text(-1.5, YMAX - 0.3, message);
sprintf(message, "Mean hits/peg %.4lg", coll_mean);
write_text(-1.5, YMAX - 0.5, message);
sprintf(message, "Stdv hits/peg %.4lg", coll_stdv);
write_text(-1.5, YMAX - 0.7, message);
/* histogram of path lengths */
for (i=1; i<NPATHBINS; i++)
{
if (npath[i] > ymax) ymax = npath[i];
total_bin += npath[i];
meanbins += i*npath[i];
meansquarebins += i*i*npath[i];
}
yscale = 0.9*(YMAX-y0)*(double)ymax;
@ -620,13 +681,13 @@ void draw_statistics()
x = xshift + (double)i*dx;
y = y0 + (double)npath[i]*YMAX/yscale;
erase_rectangle(x, y0, x+dx, y, rgb);
if (y > y0) erase_rectangle(x, y0, x+dx, y, rgb);
}
glColor3f(1.0, 1.0, 1.0);
glBegin(GL_LINE_STRIP);
/* histogramm */
for (i=1; i<NPATHBINS; i++)
/* histogram */
for (i=1; i<NPATHBINS; i++) if (npath[i] > 0)
{
x = xshift + (double)i*dx;
y = y0 + (double)npath[i]*YMAX/yscale;
@ -636,19 +697,23 @@ void draw_statistics()
glVertex2d(x+dx, y0);
glVertex2d(x, y0);
}
glVertex2d(xshift, y0);
glEnd ();
glBegin(GL_LINE_STRIP);
glVertex2d(xshift, YMAX - 0.1);
glVertex2d(xshift, y0);
glVertex2d(xshift + (double)NPATHBINS*dx, y0);
glEnd ();
for (x = 0.5; x < PATHLMAX; x+=0.5)
for (x = steppath; x < PATHLMAX; x+=steppath)
{
i = (int)(x*(double)NPATHBINS/PATHLMAX);
sprintf(message, "%.2f", x);
write_text_fixedwidth(xshift + (double)i*dx - 0.1, y0 - 0.12, message);
}
for (x = 0.1; x < PATHLMAX; x+=0.1)
for (x = steppath; x < PATHLMAX; x+=steppath)
{
i = (int)(x*(double)NPATHBINS/PATHLMAX);
glBegin(GL_LINE_STRIP);
@ -656,12 +721,24 @@ void draw_statistics()
glVertex2d(xshift + (double)i*dx, y0 + 0.025);
glEnd ();
}
ebin = (double)meanbins/(double)total_bin; /* mean bin */
len_over_bin = PATHLMAX/(double)NPATHBINS; /* conversion from bin to path length */
path_mean = ebin*len_over_bin; /* mean free path */
path_stdv = sqrt((double)meansquarebins/(double)total_bin - ebin*ebin)*len_over_bin;
sprintf(message, "free path");
write_text_fixedwidth(XMAX - 0.6, y0 - 0.12, message);
sprintf(message, "Mean free path %.4lg", (double)meanbins*PATHLMAX/(double)(total_bin*NPATHBINS));
sprintf(message, "Max free path %.4lg", max_free_path);
write_text(2.2, YMAX - 0.3, message);
sprintf(message, "Mean free path %.4lg", path_mean);
write_text(2.2, YMAX - 0.5, message);
sprintf(message, "Stdv free path %.4lg", path_stdv);
write_text(2.2, YMAX - 0.7, message);
}
@ -681,24 +758,38 @@ void animation(int circle_config)
configs[i] = (double *)malloc(8*sizeof(double));
/* init circle configuration if the domain is D_CIRCLES */
if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN == D_CIRCLES_IN_RECT)||(B_DOMAIN == D_CIRCLES_IN_GENUSN))
if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN == D_CIRCLES_IN_RECT)||(B_DOMAIN == D_CIRCLES_IN_GENUSN)
||(B_DOMAIN == D_CIRCLES_IN_TORUS))
init_circle_config_pinball(circle_config);
/* remove discs that are not in domain */
if ((B_DOMAIN == D_CIRCLES_IN_RECT)||(B_DOMAIN == D_CIRCLES_IN_GENUSN))
if ((B_DOMAIN == D_CIRCLES_IN_RECT)||(B_DOMAIN == D_CIRCLES_IN_GENUSN))
// ||(B_DOMAIN == D_CIRCLES_IN_TORUS))
for (i=0; i<ncircles; i++)
{
if (vabs(circley[i]) + circlerad[i] > 0.99) circleactive[i] = 0;
if (vabs(circlex[i]) + circlerad[i] > 0.99*LAMBDA) circleactive[i] = 0;
}
else if (B_DOMAIN == D_CIRCLES_IN_TORUS)
for (i=0; i<ncircles; i++)
{
if (vabs(circley[i]) - circlerad[i] > 1.0) circleactive[i] = 0;
if (vabs(circlex[i]) - circlerad[i] > LAMBDA) circleactive[i] = 0;
}
// for (i=0; i<ncircles; i++)
// printf("Circle %i at (%.2f, %.2f), double %i, partner %i\n", i, circlex[i], circley[i], double_circle[i], partner_circle[i]);
// 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_drop_config(0.4, 0.0, PID, 0.5*PID, configs);
init_drop_config(0.0, 0.0, -0.45*PID, 0.45*PID, configs);
// 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(0.5, 0.1, -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);
@ -795,9 +886,9 @@ void animation(int circle_config)
for (j=0; j<ncircles; j++)
if (circlecolor[j] > 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);
// 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();
@ -849,8 +940,8 @@ void display(void)
}
animation(CIRCLE_PATTERN);
// animation(C_TRI);
// animation(C_GOLDEN_SPIRAL);
// animation(C_SQUARE);
// animation(C_HEX);

43
schrodinger.c
View File

@ -39,10 +39,12 @@
/* General geometrical parameters */
#define WINWIDTH 1280 /* window width */
// #define WINWIDTH 1280 /* window width */
#define WINWIDTH 720 /* window width */
#define WINHEIGHT 720 /* window height */
#define NX 1280 /* number of grid points on x axis */
// #define NX 1280 /* number of grid points on x axis */
#define NX 720 /* number of grid points on x axis */
#define NY 720 /* number of grid points on y axis */
// #define NX 640 /* number of grid points on x axis */
// #define NY 360 /* number of grid points on y axis */
@ -52,22 +54,24 @@
#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 -2.0
#define YMAX 2.0 /* y interval for 9/16 aspect ratio */
// #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, see list in global_pdes.c */
#define B_DOMAIN 9 /* choice of domain shape */
#define B_DOMAIN 19 /* choice of domain shape */
#define CIRCLE_PATTERN 0 /* pattern of circles, see list in global_pdes.c */
#define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */
#define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */
#define LAMBDA 0.3 /* parameter controlling the dimensions of domain */
#define MU 0.05 /* parameter controlling the dimensions of domain */
#define LAMBDA 0.0 /* parameter controlling the dimensions of domain */
#define MU 1.75 /* parameter controlling the dimensions of domain */
#define NPOLY 6 /* number of sides of polygon */
#define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */
#define MDEPTH 3 /* depth of computation of Menger gasket */
@ -78,13 +82,25 @@
#define NGRIDX 15 /* 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
#define X_TARGET 0.4
#define Y_TARGET 0.7 /* shooter and target positions in laser fight */
#define ISO_XSHIFT_LEFT -1.65
#define ISO_XSHIFT_RIGHT 0.4
#define ISO_YSHIFT_LEFT -0.05
#define ISO_YSHIFT_RIGHT -0.05
#define ISO_SCALE 0.85 /* coordinates for isospectral billiards */
/* You can add more billiard tables by adapting the functions */
/* xy_in_billiard and draw_billiard in sub_wave.c */
/* Physical patameters of wave equation */
// #define DT 0.00000005
#define DT 0.000000005
#define DT 0.00000001
// #define DT 0.000000005
// #define DT 0.000000005
#define HBAR 1.0
@ -94,8 +110,9 @@
/* Parameters for length and speed of simulation */
#define NSTEPS 200 /* number of frames of movie */
#define NVID 1200 /* number of iterations between images displayed on screen */
#define NSTEPS 1400 /* number of frames of movie */
#define NVID 2000 /* number of iterations between images displayed on screen */
// #define NVID 1200 /* number of iterations between images displayed on screen */
#define NSEG 100 /* number of segments of boundary */
#define BOUNDARY_WIDTH 2 /* width of billiard boundary */
@ -111,7 +128,7 @@
/* Plot type, see list in global_pdes.c */
#define PLOT 10
#define PLOT 11
/* Color schemes, see list in global_pdes.c */
@ -133,8 +150,6 @@
#define HUEMEAN 150.0 /* mean value of hue for color scheme C_HUE */
#define HUEAMP -150.0 /* amplitude of variation of hue for color scheme C_HUE */
#include "hsluv.c"
#include "global_pdes.c"
#include "sub_wave.c"
@ -401,7 +416,7 @@ void animation()
printf("Integration step %.3lg\n", intstep);
/* initialize wave wave function */
init_coherent_state(-1.2, 0.0, 20.0, 0.0, 0.25, phi, psi, xy_in);
init_coherent_state(0.5, 0.0, 40.0, 0.0, 0.25, phi, psi, xy_in);
// init_coherent_state(0.0, 0.0, 0.0, 5.0, 0.03, phi, psi, xy_in);
// init_coherent_state(-0.5, 0.0, 1.0, 1.0, 0.05, phi, psi, xy_in);

382
sub_part_billiard.c
View File

@ -1,5 +1,8 @@
#include "colormaps.c"
#define DUMMY_ABSORBING -1000.0 /* dummy value of config[0] for absorbing circles */
#define BOUNDARY_SHIFT 100000.0 /* shift of boundary parametrisation for circles in domain */
#define DUMMY_SIDE_ABS -10000 /* dummy value of returned side for absorbing circles */
long int global_time = 0; /* counter to keep track of global time of simulation */
int nparticles=NPART;
@ -133,44 +136,6 @@ void init() /* initialisation of window */
}
void hsl_to_rgb(double h, double s, double l, double rgb[3]) /* color conversion from HSL to RGB */
/* h = hue, s = saturation, l = luminosity */
{
double c = 0.0, m = 0.0, x = 0.0;
c = (1.0 - fabs(2.0 * l - 1.0)) * s;
m = 1.0 * (l - 0.5 * c);
x = c * (1.0 - fabs(fmod(h / 60.0, 2) - 1.0));
if (h >= 0.0 && h < 60.0)
{
rgb[0] = c+m; rgb[1] = x+m; rgb[2] = m;
}
else if (h < 120.0)
{
rgb[0] = x+m; rgb[1] = c+m; rgb[2] = m;
}
else if (h < 180.0)
{
rgb[0] = m; rgb[1] = c+m; rgb[2] = x+m;
}
else if (h < 240.0)
{
rgb[0] = m; rgb[1] = x+m; rgb[2] = c+m;
}
else if (h < 300.0)
{
rgb[0] = x+m; rgb[1] = m; rgb[2] = c+m;
}
else if (h < 360.0)
{
rgb[0] = c+m; rgb[1] = m; rgb[2] = x+m;
}
else
{
rgb[0] = m; rgb[1] = m; rgb[2] = m;
}
}
void rgb_color_scheme(int i, double rgb[3]) /* color scheme */
{
@ -186,6 +151,7 @@ void rgb_color_scheme(int i, double rgb[3]) /* color scheme */
/* saturation = r, luminosity = 0.5 */
}
void rgb_color_scheme_lum(int i, double lum, double rgb[3]) /* color scheme */
{
double hue, y, r;
@ -349,6 +315,18 @@ void draw_colored_circle(double x, double y, double r, int nseg, double rgb[3])
}
void draw_initial_condition_circle(double x, double y, double r, int color)
/* draws a colored circle to mark initial condition */
{
double rgb[3];
rgb_color_scheme(color, rgb);
draw_colored_circle(x, y, r, NSEG, rgb);
rgb[0] = 0.0; rgb[1] = 0.0; rgb[2] = 0.0;
glColor3f(rgb[0], rgb[1], rgb[2]);
glLineWidth(4);
draw_circle(x, y, r, NSEG);
}
int in_circle(double x, double y, double r)
/* test whether (x,y) is in circle of radius r */
@ -1163,10 +1141,13 @@ void draw_billiard() /* draws the billiard boundary */
{
if (circleactive[k] == 2)
{
hsl_to_rgb(150.0, 0.9, 0.4, rgb);
glColor3f(rgb[0], rgb[1], rgb[2]);
// hsl_to_rgb(150.0, 0.9, 0.4, rgb);
// glColor3f(rgb[0], rgb[1], rgb[2]);
glColor3f(0.0, 1.0, 0.0);
rgb[0] = 0.0; rgb[1] = 0.9; rgb[2] = 0.0;
draw_colored_circle(circlex[k], circley[k], circlerad[k], NSEG, rgb);
}
draw_circle(circlex[k], circley[k], circlerad[k], NSEG);
else draw_circle(circlex[k], circley[k], circlerad[k], NSEG);
init_billiard_color();
}
@ -1227,6 +1208,55 @@ void draw_billiard() /* draws the billiard boundary */
glEnd ();
break;
}
case (D_CIRCLES_IN_TORUS):
{
rgb[0] = 0.0; rgb[1] = 0.0; rgb[2] = 0.0;
for (k=0; k<ncircles; k++) if (circleactive[k])
{
// printf("k = %i, color = %i\n", k, circlecolor[k]);
if (circlecolor[k] == 0) draw_colored_circle(circlex[k], circley[k], circlerad[k], NSEG, rgb);
else
{
if (newcircle[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; k<ncircles; k++) if (circleactive[k] >= 1)
{
if (circleactive[k] == 2)
{
hsl_to_rgb(150.0, 0.9, 0.4, rgb);
glColor3f(rgb[0], rgb[1], rgb[2]);
}
draw_circle(circlex[k], circley[k], circlerad[k], NSEG);
init_billiard_color();
}
/* draw shooter position for laser pattern */
if (CIRCLE_PATTERN == C_LASER)
{
hsl_to_rgb(0.0, 0.9, 0.5, rgb);
glColor3f(rgb[0], rgb[1], rgb[2]);
draw_circle(x_shooter, y_shooter, circlerad[ncircles-1], NSEG);
}
init_billiard_color();
glBegin(GL_LINE_LOOP);
glVertex2d(LAMBDA, -1.0);
glVertex2d(LAMBDA, 1.0);
glVertex2d(-LAMBDA, 1.0);
glVertex2d(-LAMBDA, -1.0);
glEnd();
break;
}
default:
{
printf("Function draw_billiard not defined for this billiard \n");
@ -3924,7 +3954,7 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */
int vcircles_xy(double config[8], double alpha, double pos[2])
/* determine initial configuration for start at point pos = (x,y) */
{
double c0, s0, b, c, t, theta, delta, margin = 1.0e-12, tmin, rlarge = 1.0e10;
double c0, s0, b, c, t, theta, delta, margin = 1.0e-12, tmin, rlarge = 1000.0;
double tval[ncircles], xint[ncircles], yint[ncircles], phiint[ncircles];
int i, nt = 0, nscat[ncircles], ntmin;
@ -3991,8 +4021,8 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */
}
else /* there is no intersection - set dummy values */
{
config[0] = 0.0;
config[1] = 0.0;
config[0] = DUMMY_ABSORBING;
config[1] = PI;
config[2] = 0.0;
config[3] = rlarge;
config[4] = pos[0]; /* start position */
@ -4012,7 +4042,7 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */
c = pos_circles(config, pos, &alpha);
vcircles_xy(config, alpha, pos);
c = vcircles_xy(config, alpha, pos);
return(c);
}
@ -4125,9 +4155,11 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */
{
config[0] = DUMMY_ABSORBING;
config[1] = PI;
// return(DUMMY_SIDE_ABS);
}
return(nscat[ntmin]);
if (ABSORBING_CIRCLES) return(DUMMY_SIDE_ABS);
else return(nscat[ntmin]);
}
else /* there is no intersection with the circles - compute intersection with boundary */
{
@ -4150,13 +4182,14 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */
c = pos_circles_in_rect(config, pos, &alpha);
vcircles_in_rect_xy(config, alpha, pos);
// vcircles_in_rect_xy(config, alpha, pos);
c = vcircles_in_rect_xy(config, alpha, pos);
return(c);
}
/****************************************************************************************/
/****************************************************************************************/
/* billiard with circular scatterers in a genus n surface (polygon with identifies opposite sides) */
/****************************************************************************************/
@ -4298,6 +4331,151 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */
}
/****************************************************************************************/
/* billiard with circular scatterers in a torus (rectangle with periodic boundary conditions) */
/****************************************************************************************/
int pos_circles_in_torus(double conf[2], double pos[2], double *alpha)
/* determine position on boundary of circle */
/* position varies between 0 and ncircles*2Pi for circles and between -BOUNDARY_SHIFT and 0 for boundary*/
/* returns number of hit circle */
{
double angle;
int ncirc, c;
if (conf[0] >= 0)
{
ncirc = (int)(conf[0]/DPI);
if (ncirc >= ncircles) ncirc = ncircles - 1;
angle = conf[0] - (double)ncirc*DPI;
pos[0] = circlex[ncirc] + circlerad[ncirc]*cos(angle);
pos[1] = circley[ncirc] + circlerad[ncirc]*sin(angle);
*alpha = angle + PID + conf[1];
return(ncirc);
}
else /* particle starts on boundary */
{
// conf[0] += 4.0*(LAMBDA + 1.0);
conf[0] += BOUNDARY_SHIFT;
c = pos_rectangle(conf, pos, alpha);
// conf[0] -= 4.0*(LAMBDA + 1.0);
conf[0] -= BOUNDARY_SHIFT;
return(-c-1);
}
}
int vcircles_in_torus_xy(double config[8], double alpha, double pos[2])
/* determine initial configuration for start at point pos = (x,y) */
// double config[8], alpha, pos[2];
{
double c0, s0, b, c, t, theta, delta, margin = 1.0e-12, tmin, rlarge = 1.0e10;
double tval[ncircles], xint[ncircles], yint[ncircles], phiint[ncircles];
int i, nt = 0, nscat[ncircles], ntmin, side;
c0 = cos(alpha);
s0 = sin(alpha);
for (i=0; i<ncircles; i++) if (circleactive[i])
{
b = (pos[0]-circlex[i])*c0 + (pos[1]-circley[i])*s0;
c = (pos[0]-circlex[i])*(pos[0]-circlex[i]) + (pos[1]-circley[i])*(pos[1]-circley[i]) - circlerad[i]*circlerad[i];
delta = b*b - c;
if (delta > margin) /* there is an intersection with circle i */
{
t = -b - sqrt(delta);
if (t > margin)
{
nscat[nt] = i;
tval[nt] = t;
xint[nt] = pos[0] + t*c0;
yint[nt] = pos[1] + t*s0;
phiint[nt] = argument(xint[nt] - circlex[i], yint[nt] - circley[i]);
/* test wether intersection is in rectangle */
if ((vabs(xint[nt]) < LAMBDA)&&(vabs(yint[nt]) < 1.0)) nt++;
}
}
}
if (nt > 0) /* there is at least one intersection */
{
/* find earliest intersection */
tmin = tval[0];
ntmin = 0;
for (i=1; i<nt; i++)
if (tval[i] < tmin)
{
tmin = tval[i];
ntmin = i;
}
while (phiint[ntmin] < 0.0) phiint[ntmin] += DPI;
while (phiint[ntmin] >= DPI) phiint[ntmin] -= DPI;
config[0] = (double)nscat[ntmin]*DPI + phiint[ntmin];
config[1] = PID - alpha + phiint[ntmin]; /* CHECK */
if (config[1] < 0.0) config[1] += DPI;
if (config[1] >= PI) config[1] -= DPI;
config[2] = 0.0; /* running time */
config[3] = module2(xint[ntmin]-pos[0], yint[ntmin]-pos[1]); /* distance to collision */
config[4] = pos[0]; /* start position */
config[5] = pos[1];
config[6] = xint[ntmin]; /* position of collision */
config[7] = yint[ntmin];
/* set dummy coordinates if circles are absorbing */
if (ABSORBING_CIRCLES)
{
config[0] = DUMMY_ABSORBING;
config[1] = PI;
}
return(nscat[ntmin]);
}
else /* there is no intersection with the circles - compute intersection with boundary */
{
side = vrectangle_xy(config, alpha, pos);
if (config[0] < 2.0*LAMBDA) config[0] = 4.0*LAMBDA + 2.0 - config[0];
else if (config[0] < 2.0*LAMBDA + 2.0) config[0] = 6.0*LAMBDA + 4.0 - config[0];
else if (config[0] < 4.0*LAMBDA + 2.0) config[0] = 4.0*LAMBDA + 2.0 - config[0];
else config[0] = 6.0*LAMBDA + 4.0 - config[0];
config[0] -= BOUNDARY_SHIFT;
config[1] = PI - config[1];
// config[0] -= 4.0*(LAMBDA+1.0);
// printf("Hit side %i\n", side);
// print_config(config);
return(side - 5);
}
}
int vcircles_in_torus(double config[8])
/* determine initial configuration when starting from boundary */
{
double pos[2], alpha;
int c;
c = pos_circles_in_torus(config, pos, &alpha);
vcircles_in_torus_xy(config, alpha, pos);
return(c);
}
/****************************************************************************************/
/* general billiard */
@ -4397,6 +4575,11 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */
return(pos_circles_in_genusn(conf, pos, alpha));
break;
}
case (D_CIRCLES_IN_TORUS):
{
return(pos_circles_in_torus(conf, pos, alpha));
break;
}
default:
{
printf("Function pos_billiard not defined for this billiard \n");
@ -4500,6 +4683,11 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */
return(vcircles_in_genusn_xy(config, alpha, pos));
break;
}
case (D_CIRCLES_IN_TORUS):
{
return(vcircles_in_torus_xy(config, alpha, pos));
break;
}
default:
{
printf("Function vbilliard_xy not defined for this billiard \n");
@ -4642,6 +4830,13 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */
return(vcircles_in_genusn(config));
break;
}
case (D_CIRCLES_IN_TORUS):
{
c = pos_circles_in_torus(config, pos, &alpha);
return(vcircles_in_torus(config));
break;
}
default:
{
printf("Function vbilliard not defined for this billiard \n");
@ -4843,6 +5038,18 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */
}
break;
}
case D_CIRCLES_IN_TORUS: /* same as D_CIRCLES_IN_RECT */
{
if ((vabs(x) >= LAMBDA)||(vabs(y) >= 1.0)) return(0);
else
{
condition = 1;
for (k=0; k<ncircles; k++)
if (circleactive[k]) condition = condition*out_circle(x-circlex[k], y-circley[k], circlerad[k]);
return(condition);
}
break;
}
default:
{
printf("Function ij_in_billiard not defined for this billiard \n");
@ -4919,6 +5126,8 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */
ncircles = NCX*(NCY+1);
dy = (YMAX - YMIN)/((double)NCY);
dx = dy*0.5*sqrt(3.0);
// dx = (XMAX - XMIN)/((double)NCX);
// dy = dx/(0.5*sqrt(3.0));
for (i = 0; i < NCX; i++)
for (j = 0; j < NCY+1; j++)
{
@ -5071,42 +5280,54 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */
break;
}
// case (C_LASER):
// {
// ncircles = 17;
//
// xx[0] = 0.5*(X_SHOOTER + X_TARGET);
// xx[1] = LAMBDA - 0.5*(X_TARGET - X_SHOOTER);
case (C_LASER):
{
ncircles = 17;
xx[0] = 0.5*(x_shooter + x_target);
xx[1] = LAMBDA - 0.5*(x_target - x_shooter);
if (xx[1] > LAMBDA) xx[1] = 2.0*LAMBDA - xx[1];
xx[2] = -xx[0];
xx[3] = -xx[1];
yy[0] = 0.5*(y_shooter + y_target);
yy[1] = 1.0 - 0.5*(y_target - y_shooter);
if (yy[1] > 1.0) yy[1] = 2.0 - yy[1];
yy[2] = -yy[0];
yy[3] = -yy[1];
// 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[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++)
// {
// circlex[4*i + j] = xx[i];
// circley[4*i + j] = yy[j];
//
// }
//
// circlex[ncircles - 1] = X_TARGET;
// circley[ncircles - 1] = Y_TARGET;
//
// for (i=0; i<ncircles - 1; i++)
// {
// circlerad[i] = MU;
// circleactive[i] = 1;
// }
//
// circlerad[ncircles - 1] = 0.5*MU;
// circleactive[ncircles - 1] = 2;
//
// break;
// }
for (i = 0; i < 4; i++)
for (j = 0; j < 4; j++)
{
circlex[4*i + j] = xx[i];
circley[4*i + j] = yy[j];
}
circlex[ncircles - 1] = x_target;
circley[ncircles - 1] = y_target;
for (i=0; i<ncircles - 1; i++)
{
circlerad[i] = MU;
circleactive[i] = 1;
}
circlerad[ncircles - 1] = 0.5*MU;
circleactive[ncircles - 1] = 2;
break;
}
default:
{
printf("Function init_circle_config not defined for this pattern \n");
@ -5114,3 +5335,6 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */
}
}

360
sub_part_pinball.c Normal file
View File

@ -0,0 +1,360 @@
/* global variables needed for circle configuration with periodic boundary conditions */
short int double_circle[NMAXCIRCLES]; /* set to 1 if a circle is a translate of another one on the boundary */
int partner_circle[NMAXCIRCLES]; /* number of circle of which current circle is a copy */
void init_circle_config_pinball(int circle_pattern)
{
int i, j, k, n, ncirc0, n_p_active, ncandidates=5000, naccepted;
double dx, dy, xx[4], yy[4], x, y, gamma, height, phi, r0, r, dpoisson = 3.25*MU;
short int active_poisson[NMAXCIRCLES], far;
switch (circle_pattern) {
case (C_FOUR_CIRCLES):
{
ncircles = 4;
circlex[0] = 1.0;
circley[0] = 0.0;
circlerad[0] = 0.8;
circlex[1] = -1.0;
circley[1] = 0.0;
circlerad[1] = 0.8;
circlex[2] = 0.0;
circley[2] = 0.8;
circlerad[2] = 0.4;
circlex[3] = 0.0;
circley[3] = -0.8;
circlerad[3] = 0.4;
for (i=0; i<4; i++)
{
circleactive[i] = 1;
circlecolor[i] = 0;
newcircle[i] = 0;
}
break;
}
case (C_SQUARE):
{
ncircles = NCX*NCY;
dy = (BOXYMAX - BOXYMIN)/((double)NCY);
// dy = (YMAX - YMIN)/((double)NCY);
for (i = 0; i < NCX; i++)
for (j = 0; j < NCY; j++)
{
n = NCY*i + j;
circlex[n] = ((double)(i-NCX/2) + 0.5)*dy;
circley[n] = BOXYMIN + ((double)j + 0.5)*dy;
circlerad[n] = MU;
circleactive[n] = 1;
circlecolor[n] = 0;
newcircle[n] = 0;
}
break;
}
case (C_HEX):
{
ncircles = NCX*(NCY+1);
dy = (YMAX - YMIN)/((double)NCY);
dx = dy*0.5*sqrt(3.0);
for (i = 0; i < NCX; i++)
for (j = 0; j < NCY+1; j++)
{
n = (NCY+1)*i + j;
// circlex[n] = ((double)(i-NCX/2) + 0.5)*dy;
circlex[n] = ((double)(i-NCX/2))*dy;
if (NCX % 2 == 0) circlex[n] += 0.5*dy;
circley[n] = YMIN + ((double)j - 0.5)*dy;
if ((i+NCX)%2 == 1) circley[n] += 0.5*dy;
circlerad[n] = MU;
circleactive[n] = 1;
circlecolor[n] = 0;
newcircle[n] = 0;
}
break;
}
case (C_TRI):
{
ncircles = NCX*(NCY+1);
dy = (BOXYMAX - BOXYMIN)/((double)NCY);
dx = dy*0.5*sqrt(3.0);
for (i = 0; i < NCX; i++)
for (j = 0; j < NCY+1; j++)
{
n = (NCY+1)*i + j;
circlex[n] = ((double)(i-NCX/2) + 0.5)*dx;
circley[n] = BOXYMIN + ((double)j)*dy;
if ((i+NCX)%2 == 1) circley[n] += 0.5*dy;
circlerad[n] = MU;
circleactive[n] = 1;
circlecolor[n] = 0;
newcircle[n] = 0;
/* take care of periodic boundary conditions */
if (B_DOMAIN == D_CIRCLES_IN_TORUS)
{
if ((j == NCY)&&((i+NCX)%2 == 0))
{
double_circle[n] = 1;
partner_circle[n] = (NCY+1)*i;
}
else
{
double_circle[n] = 0;
if ((j == 0)&&((i+NCX)%2 == 0)) partner_circle[n] = (NCY+1)*i + NCY;
else partner_circle[n] = n;
}
}
else double_circle[n] = 0;
}
break;
}
case (C_GOLDEN_MEAN):
{
ncircles = NCX*NCY;
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++)
{
circlex[n] = -LAMBDA + n*dx;
circley[n] = y;
y += height*gamma;
if (y > YMAX) y -= height;
circlerad[n] = MU;
circleactive[n] = 1;
circlecolor[n] = 0;
newcircle[n] = 0;
}
break;
}
case (C_GOLDEN_SPIRAL):
{
ncircles = 1;
circlex[0] = 0.0;
circley[0] = 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<NGOLDENSPIRAL; i++)
{
x = r*cos(phi);
y = r*sin(phi);
phi += gamma;
r += MU*r0/r;
if ((vabs(x) < LAMBDA)&&(vabs(y) < YMAX + MU))
{
circlex[ncircles] = x;
circley[ncircles] = y;
ncircles++;
}
}
for (i=0; i<ncircles; i++)
{
circlerad[i] = MU;
circlecolor[i] = 0;
newcircle[i] = 0;
/* inactivate circles outside the domain */
if ((circley[i] < YMAX + MU)&&(circley[i] > YMIN - MU)) circleactive[i] = 1;
}
break;
}
case (C_RAND_DISPLACED):
{
ncircles = NCX*NCY;
dy = (YMAX - YMIN)/((double)NCY);
for (i = 0; i < NCX; i++)
for (j = 0; j < NCY; j++)
{
n = NCY*i + j;
circlex[n] = ((double)(i-NCX/2) + 0.5*((double)rand()/RAND_MAX - 0.5))*dy;
circley[n] = YMIN + ((double)j + 0.5 + 0.5*((double)rand()/RAND_MAX - 0.5))*dy;
circlerad[n] = MU*sqrt(1.0 + 0.8*((double)rand()/RAND_MAX - 0.5));
circleactive[n] = 1;
circlecolor[n] = 0;
newcircle[n] = 0;
}
break;
}
case (C_RAND_POISSON):
{
ncircles = NPOISSON;
for (n = 0; n < NPOISSON; n++)
{
circlex[n] = LAMBDA*(2.0*(double)rand()/RAND_MAX - 1.0);
circley[n] = (BOXYMAX - BOXYMIN)*(double)rand()/RAND_MAX + BOXYMIN;
circlerad[n] = MU;
circleactive[n] = 1;
circlecolor[n] = 0;
newcircle[n] = 0;
double_circle[n] = 0;
partner_circle[n] = n;
/* take care of periodic boundary conditions */
if (B_DOMAIN == D_CIRCLES_IN_TORUS)
{
/* inactivate circles in corners */
if ((vabs(circlex[n]) > LAMBDA - MU)&&((circley[n] < BOXYMIN + MU)||(circley[n] > BOXYMAX - MU)))
circleactive[n] = 0;
if (circlex[n] < - LAMBDA + MU)
{
circlex[ncircles] = circlex[n] + 2.0*LAMBDA;
circley[ncircles] = circley[n];
partner_circle[ncircles] = n;
partner_circle[n] = ncircles;
ncircles++;
}
else if (circlex[n] > LAMBDA - MU)
{
circlex[ncircles] = circlex[n] - 2.0*LAMBDA;
circley[ncircles] = circley[n];
partner_circle[ncircles] = n;
partner_circle[n] = ncircles;
ncircles++;
}
if (circley[n] < BOXYMIN + MU)
{
circlex[ncircles] = circlex[n];
circley[ncircles] = circley[n] + BOXYMAX - BOXYMIN;
partner_circle[ncircles] = n;
partner_circle[n] = ncircles;
ncircles++;
}
else if (circley[n] > BOXYMAX - MU)
{
circlex[ncircles] = circlex[n];
circley[ncircles] = circley[n] - BOXYMAX + BOXYMIN;
partner_circle[ncircles] = n;
partner_circle[n] = ncircles;
ncircles++;
}
}
}
printf("%i circles\n", ncircles);
if (B_DOMAIN == D_CIRCLES_IN_TORUS) for (n = NPOISSON; n < ncircles; n++)
{
// printf("circle %i at (%.3f, %.3f)\n", n, circlex[n], circley[n]);
circlerad[n] = MU;
if (circleactive[partner_circle[n]]) circleactive[n] = 1;
circlecolor[n] = 0;
newcircle[n] = 0;
double_circle[n] = 1;
}
break;
}
case (C_POISSON_DISC):
{
printf("Generating Poisson disc sample\n");
/* generate first circle */
circlex[0] = LAMBDA*(2.0*(double)rand()/RAND_MAX - 1.0);
circley[0] = (YMAX - YMIN)*(double)rand()/RAND_MAX + YMIN;
active_poisson[0] = 1;
n_p_active = 1;
ncircles = 1;
while ((n_p_active > 0)&&(ncircles < NMAXCIRCLES))
{
/* randomly select an active circle */
i = rand()%(ncircles);
while (!active_poisson[i]) i = rand()%(ncircles);
// printf("Starting from circle %i at (%.3f,%.3f)\n", i, circlex[i], circley[i]);
/* generate new candidates */
naccepted = 0;
for (j=0; j<ncandidates; j++)
{
r = dpoisson*(2.0*(double)rand()/RAND_MAX + 1.0);
phi = DPI*(double)rand()/RAND_MAX;
x = circlex[i] + r*cos(phi);
y = circley[i] + r*sin(phi);
// printf("Testing new circle at (%.3f,%.3f)\t", x, y);
far = 1;
for (k=0; k<ncircles; k++) if ((k!=i))
{
/* new circle is far away from circle k */
far = far*((x - circlex[k])*(x - circlex[k]) + (y - circley[k])*(y - circley[k]) >= 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);
circlex[ncircles] = x;
circley[ncircles] = y;
circlerad[ncircles] = MU;
circleactive[ncircles] = 1;
circlecolor[ncircles] = 0;
newcircle[ncircles] = 0;
active_poisson[ncircles] = 1;
ncircles++;
n_p_active++;
naccepted++;
}
// else printf("Rejected\n");
}
if (naccepted == 0) /* inactivate circle i */
{
// printf("No candidates work, inactivate circle %i\n", i);
active_poisson[i] = 0;
n_p_active--;
}
printf("%i active circles\n", n_p_active);
}
printf("Generated %i circles\n", ncircles);
break;
}
// case (C_LASER):
// {
// ncircles = 17;
//
// xx[0] = 0.5*(X_SHOOTER + X_TARGET);
// xx[1] = LAMBDA - 0.5*(X_TARGET - X_SHOOTER);
// 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++)
// {
// circlex[4*i + j] = xx[i];
// circley[4*i + j] = yy[j];
//
// }
//
// circlex[ncircles - 1] = X_TARGET;
// circley[ncircles - 1] = Y_TARGET;
//
// for (i=0; i<ncircles - 1; i++)
// {
// circlerad[i] = MU;
// circleactive[i] = 1;
// }
//
// circlerad[ncircles - 1] = 0.5*MU;
// circleactive[ncircles - 1] = 2;
//
// break;
// }
default:
{
printf("Function init_circle_config not defined for this pattern \n");
}
}
}

406
sub_wave.c
View File

@ -2,6 +2,9 @@
/* Graphics routines */
/*********************/
#include "colors_waves.c"
int writetiff(char *filename, char *description, int x, int y, int width, int height, int compression)
{
TIFF *file;
@ -67,149 +70,8 @@ void init() /* initialisation of window */
void hsl_to_rgb_jet(double h, double s, double l, double rgb[3]) /* color conversion from HSL to RGB */
/* h = hue, s = saturation, l = luminosity */
{
double c = 0.0, m = 0.0, x = 0.0;
c = (1.0 - fabs(2.0 * l - 1.0)) * s;
m = 1.0 * (l - 0.5 * c);
x = c * (1.0 - fabs(fmod(h / 60.0, 2) - 1.0));
if (h >= 0.0 && h < 60.0)
{
rgb[0] = c+m; rgb[1] = x+m; rgb[2] = m;
}
else if (h < 120.0)
{
rgb[0] = x+m; rgb[1] = c+m; rgb[2] = m;
}
else if (h < 180.0)
{
rgb[0] = m; rgb[1] = c+m; rgb[2] = x+m;
}
else if (h < 240.0)
{
rgb[0] = m; rgb[1] = x+m; rgb[2] = c+m;
}
else if (h < 300.0)
{
rgb[0] = x+m; rgb[1] = m; rgb[2] = c+m;
}
else if (h < 360.0)
{
rgb[0] = c+m; rgb[1] = m; rgb[2] = x+m;
}
else
{
rgb[0] = m; rgb[1] = m; rgb[2] = m;
}
}
void hsl_to_rgb(double h, double s, double l, double rgb[3]) /* color conversion from HSL to RGB */
{
double r, g, b;
switch (COLOR_PALETTE) {
case (COL_JET):
{
hsl_to_rgb_jet(h, s, l, rgb);
break;
}
case (COL_HSLUV):
{
hsluv2rgb(h, 100.0*s, l, &r, &g, &b);
rgb[0] = r*100.0;
rgb[1] = g*100.0;
rgb[2] = b*100.0;
break;
}
}
// if (HSLUV_COLORS)
// {
// hsluv2rgb(h, 100.0*s, l, &r, &g, &b);
// rgb[0] = r*100.0;
// rgb[1] = g*100.0;
// rgb[2] = b*100.0;
// }
// else hsl_to_rgb_jet(h, s, l, rgb);
}
double color_amplitude(double value, double scale, int time)
/* transforms the wave amplitude into a double in [-1,1] to feed into color scheme */
{
return(tanh(SLOPE*value/scale)*exp(-((double)time*ATTENUATION)));
}
void color_scheme(int scheme, double value, double scale, int time, double rgb[3]) /* color scheme */
{
double hue, y, r, amplitude;
int intpart;
/* saturation = r, luminosity = y */
switch (scheme) {
case (C_LUM):
{
hue = COLORHUE + (double)time*COLORDRIFT/(double)NSTEPS;
if (hue < 0.0) hue += 360.0;
if (hue >= 360.0) hue -= 360.0;
r = 0.9;
amplitude = color_amplitude(value, scale, time);
y = LUMMEAN + amplitude*LUMAMP;
intpart = (int)y;
y -= (double)intpart;
hsl_to_rgb(hue, r, y, rgb);
break;
}
case (C_HUE):
{
r = 0.9;
amplitude = color_amplitude(value, scale, time);
y = 0.5;
hue = HUEMEAN + amplitude*HUEAMP;
if (hue < 0.0) hue += 360.0;
if (hue >= 360.0) hue -= 360.0;
hsl_to_rgb(hue, r, y, rgb);
break;
}
}
}
void color_scheme_lum(int scheme, double value, double scale, int time, double lum, double rgb[3]) /* color scheme */
{
double hue, y, r, amplitude;
int intpart;
/* saturation = r, luminosity = y */
switch (scheme) {
case (C_LUM):
{
hue = COLORHUE + (double)time*COLORDRIFT/(double)NSTEPS;
if (hue < 0.0) hue += 360.0;
if (hue >= 360.0) hue -= 360.0;
r = 0.9;
amplitude = color_amplitude(value, scale, time);
y = LUMMEAN + amplitude*LUMAMP;
intpart = (int)y;
y -= (double)intpart;
hsl_to_rgb(hue, r, y, rgb);
break;
}
case (C_HUE):
{
r = 0.9;
amplitude = color_amplitude(value, scale, time);
y = lum;
hue = HUEMEAN + amplitude*HUEAMP;
if (hue < 0.0) hue += 360.0;
if (hue >= 360.0) hue -= 360.0;
hsl_to_rgb(hue, r, y, rgb);
break;
}
}
}
void blank()
@ -502,7 +364,7 @@ void init_circle_config()
/* 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;
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];
short int active_poisson[NMAXCIRCLES], far;
switch (CIRCLE_PATTERN) {
@ -621,6 +483,42 @@ void init_circle_config()
}
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++)
{
circlex[4*i + j] = xx[i];
circley[4*i + j] = yy[j];
}
circlex[ncircles - 1] = X_TARGET;
circley[ncircles - 1] = Y_TARGET;
for (i=0; i<ncircles - 1; i++)
{
circlerad[i] = MU;
circleactive[i] = 1;
}
circlerad[ncircles - 1] = 0.5*MU;
circleactive[ncircles - 1] = 2;
break;
}
case (C_POISSON_DISC):
{
printf("Generating Poisson disc sample\n");
@ -628,6 +526,7 @@ void init_circle_config()
circlex[0] = LAMBDA*(2.0*(double)rand()/RAND_MAX - 1.0);
circley[0] = (YMAX - YMIN)*(double)rand()/RAND_MAX + YMIN;
active_poisson[0] = 1;
// circleactive[0] = 1;
n_p_active = 1;
ncircles = 1;
@ -811,12 +710,104 @@ void init_circle_config()
}
int axial_symmetry(double z1[2], double z2[2], double z[2], double zprime[2])
/* compute reflection of point z wrt axis through z1 and z2 */
{
double u[2], r, zdotu, zparallel[2], zperp[2];
/* compute unit vector parallel to z1-z2 */
u[0] = z2[0] - z1[0];
u[1] = z2[1] - z1[1];
r = module2(u[0], u[1]);
if (r == 0) return(0); /* z1 and z2 are the same */
u[0] = u[0]/r;
u[1] = u[1]/r;
// printf("u = (%.2f, %.2f)\n", u[0], u[1]);
/* projection of z1z on z1z2 */
zdotu = (z[0] - z1[0])*u[0] + (z[1] - z1[1])*u[1];
zparallel[0] = zdotu*u[0];
zparallel[1] = zdotu*u[1];
// printf("zparallel = (%.2f, %.2f)\n", zparallel[0], zparallel[1]);
/* normal vector to z1z2 */
zperp[0] = z[0] - z1[0] - zparallel[0];
zperp[1] = z[1] - z1[1] - zparallel[1];
// printf("zperp = (%.2f, %.2f)\n", zperp[0], zperp[1]);
/* reflected point */
zprime[0] = z[0] - 2.0*zperp[0];
zprime[1] = z[1] - 2.0*zperp[1];
return(1);
}
void compute_isospectral_coordinates(int type, double xshift, double yshift, double scaling, double vertex[9][2])
/* 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*/
{
vertex[0][0] = (xshift - 0.5)*scaling;
vertex[0][1] = (yshift - 0.25)*scaling;
vertex[1][0] = (0.5 + xshift)*scaling;
vertex[1][1] = (yshift - 0.25)*scaling;
vertex[2][0] = (LAMBDA + xshift - 0.5)*scaling;
vertex[2][1] = (MU + yshift - 0.25)*scaling;
axial_symmetry(vertex[0], vertex[2], vertex[1], vertex[3]);
axial_symmetry(vertex[0], vertex[1], vertex[2], vertex[4]);
axial_symmetry(vertex[1], vertex[2], vertex[0], vertex[5]);
if (type == 0)
{
axial_symmetry(vertex[0], vertex[3], vertex[2], vertex[6]);
axial_symmetry(vertex[1], vertex[4], vertex[0], vertex[7]);
axial_symmetry(vertex[2], vertex[5], vertex[1], vertex[8]);
}
else
{
axial_symmetry(vertex[2], vertex[3], vertex[0], vertex[6]);
axial_symmetry(vertex[0], vertex[4], vertex[1], vertex[7]);
axial_symmetry(vertex[1], vertex[5], vertex[2], vertex[8]);
}
}
void isospectral_initial_point(double x, double y, double left[2], double right[2])
/* compute initial coordinates in isospectral billiards */
{
left[0] = (x + ISO_XSHIFT_LEFT)*ISO_SCALE;
left[1] = (y + ISO_YSHIFT_LEFT)*ISO_SCALE;
right[0] = (x + ISO_XSHIFT_RIGHT)*ISO_SCALE;
right[1] = (y + ISO_YSHIFT_RIGHT)*ISO_SCALE;
}
int xy_in_triangle(double x, double y, double z1[2], double z2[2], double z3[2])
/* returns 1 iff (x,y) is inside the triangle with vertices z1, z2, z3 */
{
double v1, v2, v3;
/* compute wedge products */
v1 = (z2[0] - z1[0])*(y - z1[1]) - (z2[1] - z1[1])*(x - z1[0]);
v2 = (z3[0] - z2[0])*(y - z2[1]) - (z3[1] - z2[1])*(x - z2[0]);
v3 = (z1[0] - z3[0])*(y - z3[1]) - (z1[1] - z3[1])*(x - z3[0]);
if ((v1 >= 0.0)&&(v2 >= 0.0)&&(v3 >= 0.0)) return(1);
else return(0);
}
int xy_in_billiard(double x, double y)
/* returns 1 if (x,y) represents a point in the billiard */
// double x, y;
{
double l2, r2, r2mu, omega, b, c, angle, z, x1, y1, x2, y2, u, v, u1, v1, dx, dy, width;
int i, j, k, k1, k2, condition;
static double vertex[9][2], wertex[9][2];
static int first = 1;
switch (B_DOMAIN) {
case (D_RECTANGLE):
@ -1048,6 +1039,33 @@ int xy_in_billiard(double x, double y)
else return(1);
}
}
case (D_ISOSPECTRAL):
{
if (first)
{
compute_isospectral_coordinates(0, ISO_XSHIFT_LEFT, ISO_YSHIFT_LEFT, ISO_SCALE, vertex);
compute_isospectral_coordinates(1, ISO_XSHIFT_RIGHT, ISO_YSHIFT_RIGHT, ISO_SCALE, wertex);
first = 0;
}
/* 1st triangle */
condition = xy_in_triangle(x, y, vertex[0], vertex[1], vertex[2]);
condition += xy_in_triangle(x, y, vertex[0], vertex[4], vertex[1]);
condition += xy_in_triangle(x, y, vertex[1], vertex[5], vertex[2]);
condition += xy_in_triangle(x, y, vertex[0], vertex[2], vertex[3]);
condition += xy_in_triangle(x, y, vertex[1], vertex[4], vertex[7]);
condition += xy_in_triangle(x, y, vertex[2], vertex[5], vertex[8]);
condition += xy_in_triangle(x, y, vertex[0], vertex[3], vertex[6]);
/* 2nd triangle */
condition += xy_in_triangle(x, y, wertex[0], wertex[1], wertex[2]);
condition += xy_in_triangle(x, y, wertex[0], wertex[4], wertex[1]);
condition += xy_in_triangle(x, y, wertex[1], wertex[5], wertex[2]);
condition += xy_in_triangle(x, y, wertex[0], wertex[2], wertex[3]);
condition += xy_in_triangle(x, y, wertex[0], wertex[7], wertex[4]);
condition += xy_in_triangle(x, y, wertex[1], wertex[8], wertex[5]);
condition += xy_in_triangle(x, y, wertex[2], wertex[6], wertex[3]);
return(condition >= 1);
}
case (D_CIRCLES):
{
for (i = 0; i < ncircles; i++)
@ -1060,6 +1078,19 @@ int xy_in_billiard(double x, double y)
}
return(1);
}
case (D_CIRCLES_IN_RECT): /* returns 2 inside circles, 0 outside rectangle */
{
for (i = 0; i < ncircles; i++)
if (circleactive[i])
{
x1 = circlex[i];
y1 = circley[i];
r2 = circlerad[i]*circlerad[i];
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_MENGER):
{
x1 = 0.5*(x+1.0);
@ -1225,10 +1256,22 @@ void toka_lineto(double x1, double y1)
glVertex2d(pos[0], pos[1]);
}
void iso_lineto(double z[2])
/* draws boundary segments of isospectral billiard */
{
double pos[2];
xy_to_pos(z[0], z[1], pos);
glVertex2d(pos[0], pos[1]);
}
void draw_billiard() /* draws the billiard boundary */
{
double x0, x, y, x1, y1, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z, l, width, a, b, c, ymax;
static double vertex[9][2], wertex[9][2];
int i, j, k, k1, k2, mr2;
static int first = 1;
if (BLACK) glColor3f(1.0, 1.0, 1.0);
else glColor3f(0.0, 0.0, 0.0);
@ -1635,7 +1678,7 @@ void draw_billiard() /* draws the billiard boundary */
glColor3f(0.3, 0.3, 0.3);
for (k=0; k<NPOLY; k++)
{
alpha = APOLY*PID + (k+0.5)*omega;
alpha = APOLY*PID + (2.0*(double)k+1.0)*omega;
draw_circle(LAMBDA*cos(alpha), LAMBDA*sin(alpha), r, NSEG);
}
}
@ -1795,6 +1838,70 @@ void draw_billiard() /* draws the billiard boundary */
}
break;
}
case (D_ISOSPECTRAL):
{
if (first)
{
compute_isospectral_coordinates(0, ISO_XSHIFT_LEFT, ISO_YSHIFT_LEFT, ISO_SCALE, vertex);
compute_isospectral_coordinates(1, ISO_XSHIFT_RIGHT, ISO_YSHIFT_RIGHT, ISO_SCALE, wertex);
// compute_isospectral_coordinates(0, -1.0, 0.05, 0.9, vertex);
// compute_isospectral_coordinates(1, 0.9, 0.05, 0.9, wertex);
// for (i=0; i<9; i++) printf("(x%i, y%i) = (%.2f, %.2f)\n", i, i, vertex[i][0], vertex[i][1]);
first = 0;
}
/* 1st triangle */
glBegin(GL_LINE_LOOP);
iso_lineto(vertex[0]);
iso_lineto(vertex[4]);
iso_lineto(vertex[7]);
iso_lineto(vertex[1]);
iso_lineto(vertex[5]);
iso_lineto(vertex[8]);
iso_lineto(vertex[2]);
iso_lineto(vertex[3]);
iso_lineto(vertex[6]);
glEnd();
/* inner lines */
glBegin(GL_LINE_LOOP);
iso_lineto(vertex[0]);
iso_lineto(vertex[1]);
iso_lineto(vertex[2]);
iso_lineto(vertex[0]);
iso_lineto(vertex[3]);
iso_lineto(vertex[2]);
iso_lineto(vertex[5]);
iso_lineto(vertex[1]);
iso_lineto(vertex[4]);
glEnd();
/* 2nd triangle */
glBegin(GL_LINE_LOOP);
iso_lineto(wertex[0]);
iso_lineto(wertex[7]);
iso_lineto(wertex[4]);
iso_lineto(wertex[1]);
iso_lineto(wertex[8]);
iso_lineto(wertex[5]);
iso_lineto(wertex[2]);
iso_lineto(wertex[6]);
iso_lineto(wertex[3]);
glEnd();
/* inner lines */
glBegin(GL_LINE_LOOP);
iso_lineto(wertex[0]);
iso_lineto(wertex[1]);
iso_lineto(wertex[2]);
iso_lineto(wertex[0]);
iso_lineto(wertex[4]);
iso_lineto(wertex[1]);
iso_lineto(wertex[5]);
iso_lineto(wertex[2]);
iso_lineto(wertex[3]);
glEnd();
break;
}
case (D_CIRCLES):
{
glLineWidth(BOUNDARY_WIDTH);
@ -1802,6 +1909,19 @@ void draw_billiard() /* draws the billiard boundary */
if (circleactive[i]) draw_circle(circlex[i], circley[i], circlerad[i], NSEG);
break;
}
case (D_CIRCLES_IN_RECT):
{
glLineWidth(BOUNDARY_WIDTH);
for (i = 0; i < ncircles; i++)
if (circleactive[i]) draw_circle(circlex[i], circley[i], circlerad[i], 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_MENGER):
{
glLineWidth(3);

10
turbo_colormap.c Normal file
View File

File diff suppressed because one or more lines are too long

163
wave_billiard.c
View File

@ -48,37 +48,31 @@
/* General geometrical parameters */
#define WINWIDTH 1280 /* window width */
// #define WINWIDTH 720 /* window width */
#define WINHEIGHT 720 /* window height */
#define NX 1280 /* number of grid points on x axis */
// #define NX 720 /* number of grid points on x axis */
#define NY 720 /* 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 XMIN -1.6
// #define XMAX 1.6 /* x interval */
// #define YMIN -1.6
// #define YMAX 1.6 /* y interval for 9/16 aspect ratio */
#define JULIA_SCALE 1.0 /* scaling for Julia sets */
/* Choice of the billiard table */
#define B_DOMAIN 34 /* choice of domain shape, see list in global_pdes.c */
#define B_DOMAIN 32 /* choice of domain shape, see list in global_pdes.c */
#define CIRCLE_PATTERN 8 /* pattern of circles, see list in global_pdes.c */
#define CIRCLE_PATTERN 7 /* pattern of circles, see list in global_pdes.c */
#define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */
#define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */
#define LAMBDA 0.6 /* parameter controlling the dimensions of domain */
#define MU 0.3 /* parameter controlling the dimensions of domain */
#define NPOLY 3 /* number of sides of polygon */
#define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */
#define LAMBDA 0.0 /* parameter controlling the dimensions of domain */
#define MU 1.0 /* parameter controlling the dimensions of domain */
#define NPOLY 7 /* 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 */
@ -87,6 +81,17 @@
#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
#define X_TARGET 0.4
#define Y_TARGET 0.7 /* shooter and target positions in laser fight */
#define ISO_XSHIFT_LEFT -1.65
#define ISO_XSHIFT_RIGHT 0.4
#define ISO_YSHIFT_LEFT -0.05
#define ISO_YSHIFT_RIGHT -0.05
#define ISO_SCALE 0.85 /* coordinates for isospectral billiards */
/* You can add more billiard tables by adapting the functions */
/* xy_in_billiard and draw_billiard below */
@ -99,9 +104,9 @@
#define OMEGA 0.002 /* frequency of periodic excitation */
#define AMPLITUDE 1.0 /* amplitude of periodic excitation */
#define COURANT 0.02 /* Courant number */
#define COURANTB 0.01 /* Courant number in medium B */
#define COURANTB 0.02 /* Courant number in medium B */
#define GAMMA 0.0 /* damping factor in wave equation */
#define GAMMAB 1.0e-6 /* damping factor in wave equation */
#define GAMMAB 5.0e-3 /* 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 */
#define KAPPA 0.0 /* "elasticity" term enforcing oscillations */
@ -118,8 +123,8 @@
/* Parameters for length and speed of simulation */
#define NSTEPS 4050 /* number of frames of movie */
#define NVID 32 /* number of iterations between images displayed on screen */
#define NSTEPS 5000 /* number of frames of movie */
#define NVID 50 /* 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 */
@ -133,48 +138,49 @@
/* Parameters of initial condition */
#define INITIAL_AMP 0.2 /* amplitude of initial condition */
#define INITIAL_VARIANCE 0.002 /* variance of initial condition */
#define INITIAL_WAVELENGTH 0.1 /* wavelength of initial condition */
#define INITIAL_AMP 0.75 /* amplitude of initial condition */
#define INITIAL_VARIANCE 0.0005 /* variance of initial condition */
#define INITIAL_WAVELENGTH 0.02 /* wavelength of initial condition */
/* Plot type, see list in global_pdes.c */
#define PLOT 1
#define PLOT_B 0 /* plot type for second movie */
#define PLOT_B 3 /* plot type for second movie */
/* Color schemes */
#define COLOR_PALETTE 0 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE 13 /* Color palette, see list in global_pdes.c */
#define BLACK 1 /* background */
#define COLOR_SCHEME 1 /* choice of color scheme, see list in global_pdes.c */
#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.08 /* sensitivity of color on wave amplitude */
// #define SLOPE 0.05 /* sensitivity of color on wave amplitude */
#define SLOPE 0.15 /* sensitivity of color on wave amplitude */
#define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */
#define E_SCALE 200.0 /* scaling factor for energy representation */
// #define E_SCALE 2500.0 /* scaling factor for energy representation */
#define E_SCALE 100.0 /* scaling factor for energy 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 220.0 /* mean value of hue for color scheme C_HUE */
#define HUEAMP -230.0 /* amplitude of variation of hue for color scheme C_HUE */
#define HUEMEAN 180.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 DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */
#define SAVE_TIME_SERIES 0 /* set to 1 to save wave time series at a point */
/* For debugging purposes only */
#define FLOOR 0 /* set to 1 to limit wave amplitude to VMAX */
#define VMAX 10.0 /* max value of wave amplitude */
#include "hsluv.c"
#include "global_pdes.c" /* constants and global variables */
#include "sub_wave.c" /* common functions for wave_billiard, heat and schrodinger */
#include "wave_common.c" /* common functions for wave_billiard, wave_comparison, etc */
FILE *time_series_left, *time_series_right;
double courant2, courantb2; /* Courant parameters squared */
@ -200,11 +206,18 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
#pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta,x,y,c,cc,gamma)
for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
if (xy_in[i][j])
// if (xy_in[i][j])
// {
// c = COURANT;
// cc = courant2;
// gamma = GAMMA;
// }
if (xy_in[i][j] != 0)
{
c = COURANT;
cc = courant2;
gamma = GAMMA;
if (xy_in[i][j] == 1) gamma = GAMMA;
else gamma = GAMMAB;
}
else if (TWOSPEEDS)
{
@ -213,7 +226,7 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
gamma = GAMMAB;
}
if ((TWOSPEEDS)||(xy_in[i][j])){
if ((TWOSPEEDS)||(xy_in[i][j] != 0)){
/* discretized Laplacian for various boundary conditions */
if ((B_COND == BC_DIRICHLET)||(B_COND == BC_ABSORBING))
{
@ -320,11 +333,18 @@ void evolve_wave(double *phi[NX], double *psi[NX], double *phi_tmp[NX], double *
void animation()
{
double time, scale;
double *phi[NX], *psi[NX], *phi_tmp[NX], *psi_tmp[NX];
double time, scale, ratio, startleft[2], startright[2];
double *phi[NX], *psi[NX], *phi_tmp[NX], *psi_tmp[NX], *total_energy[NX];
short int *xy_in[NX];
int i, j, s;
int i, j, s, sample_left[2], sample_right[2];
static int counter = 0;
long int wave_value;
if (SAVE_TIME_SERIES)
{
time_series_left = fopen("wave_left.dat", "w");
time_series_right = fopen("wave_right.dat", "w");
}
/* Since NX and NY are big, it seemed wiser to use some memory allocation here */
for (i=0; i<NX; i++)
@ -333,24 +353,46 @@ void animation()
psi[i] = (double *)malloc(NY*sizeof(double));
phi_tmp[i] = (double *)malloc(NY*sizeof(double));
psi_tmp[i] = (double *)malloc(NY*sizeof(double));
total_energy[i] = (double *)malloc(NY*sizeof(double));
xy_in[i] = (short int *)malloc(NY*sizeof(short int));
}
/* initialise positions and radii of circles */
if (B_DOMAIN == D_CIRCLES) init_circle_config();
if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN == D_CIRCLES_IN_RECT)) init_circle_config();
courant2 = COURANT*COURANT;
courantb2 = COURANTB*COURANTB;
/* initialize wave with a drop at one point, zero elsewhere */
init_circular_wave(0.0, -LAMBDA, phi, psi, xy_in);
// init_circular_wave(0.0, -LAMBDA, phi, psi, xy_in);
/* initialize total energy table */
if ((PLOT == P_MEAN_ENERGY)||(PLOT_B == P_MEAN_ENERGY))
for (i=0; i<NX; i++)
for (j=0; j<NY; j++)
total_energy[i][j] = 0.0;
ratio = (XMAX - XMIN)/8.4; /* for Tokarsky billiard */
isospectral_initial_point(0.25, 0.0, startleft, startright); /* for isospectral billiards */
xy_to_ij(startleft[0], startleft[1], sample_left);
xy_to_ij(startright[0], startright[1], sample_right);
// printf("xleft = (%.3f, %.3f) xright = (%.3f, %.3f)\n", xin_left, yin_left, xin_right, yin_right);
// init_wave_flat(phi, psi, xy_in);
// init_wave_plus(LAMBDA - 0.3*MU, 0.5*MU, phi, psi, xy_in);
// init_wave(LAMBDA - 0.3*MU, 0.5*MU, phi, psi, xy_in);
// init_wave(0.0, 0.0, phi, psi, xy_in);
// init_circular_wave(X_SHOOTER, Y_SHOOTER, phi, psi, xy_in);
init_circular_wave(-LAMBDA, 0.0, phi, psi, xy_in);
// init_circular_wave(0.5, 0.5, phi, psi, xy_in);
// add_circular_wave(-1.0, 0.0, LAMBDA, phi, psi, xy_in);
// add_circular_wave(1.0, -LAMBDA, 0.0, phi, psi, xy_in);
// add_circular_wave(-1.0, 0.0, -LAMBDA, phi, psi, xy_in);
// init_circular_wave_xplusminus(startleft[0], startleft[1], startright[0], startright[1], phi, psi, xy_in);
// init_circular_wave_xplusminus(-0.9, 0.0, 0.81, 0.0, phi, psi, xy_in);
// init_circular_wave(-2.0*ratio, 0.0, phi, psi, xy_in);
// init_planar_wave(XMIN + 0.015, 0.0, phi, psi, xy_in);
// init_planar_wave(XMIN + 0.02, 0.0, phi, psi, xy_in);
// init_planar_wave(XMIN + 0.8, 0.0, phi, psi, xy_in);
@ -364,8 +406,11 @@ void animation()
blank();
glColor3f(0.0, 0.0, 0.0);
draw_wave(phi, psi, xy_in, 1.0, 0, PLOT);
// draw_wave(phi, psi, xy_in, 1.0, 0, PLOT);
draw_wave_e(phi, psi, total_energy, xy_in, 1.0, 0, PLOT);
draw_billiard();
if (DRAW_COLOR_SCHEME) draw_color_scheme(1.7, YMIN + 0.1, 1.9, YMAX - 0.1, PLOT, -12.0, 12.0);
glutSwapBuffers();
@ -385,21 +430,31 @@ void animation()
}
else scale = 1.0;
draw_wave(phi, psi, xy_in, scale, i, PLOT);
// draw_wave(phi, psi, xy_in, scale, i, PLOT);
draw_wave_e(phi, psi, total_energy, xy_in, scale, i, PLOT);
for (j=0; j<NVID; j++)
{
evolve_wave(phi, psi, phi_tmp, psi_tmp, xy_in);
if (SAVE_TIME_SERIES)
{
wave_value = (long int)(phi[sample_left[0]][sample_left[1]]*1.0e16);
fprintf(time_series_left, "%019ld\n", wave_value);
wave_value = (long int)(phi[sample_right[0]][sample_right[1]]*1.0e16);
fprintf(time_series_right, "%019ld\n", wave_value);
if ((j == 0)&&(i%10 == 0)) printf("Frame %i of %i\n", i, NSTEPS);
// fprintf(time_series_right, "%.15f\n", phi[sample_right[0]][sample_right[1]]);
}
// if (i % 10 == 9) oscillate_linear_wave(0.2*scale, 0.15*(double)(i*NVID + j), -1.5, YMIN, -1.5, YMAX, phi, psi);
}
draw_billiard();
if (DRAW_COLOR_SCHEME) draw_color_scheme(1.7, YMIN + 0.1, 1.9, YMAX - 0.1, PLOT, -12.0, 12.0);
/* add oscillating waves */
// if (i%160 == 159)
if (i%150 == 149)
if (i%345 == 344)
{
add_circular_wave(1.0, 0.0, LAMBDA, phi, psi, xy_in);
add_circular_wave(1.0, 0.0, -LAMBDA, phi, psi, xy_in);
add_circular_wave(1.0, -LAMBDA, 0.0, phi, psi, xy_in);
}
glutSwapBuffers();
@ -411,7 +466,9 @@ void animation()
if ((i >= INITIAL_TIME)&&(DOUBLE_MOVIE))
{
draw_wave(phi, psi, xy_in, scale, i, PLOT_B);
// draw_wave(phi, psi, xy_in, scale, i, PLOT_B);
draw_wave_e(phi, psi, total_energy, xy_in, scale, i, PLOT_B);
if (DRAW_COLOR_SCHEME) draw_color_scheme(1.7, YMIN + 0.1, 1.9, YMAX - 0.1, PLOT_B, -12.0, 12.0);
draw_billiard();
glutSwapBuffers();
save_frame_counter(NSTEPS + 21 + counter);
@ -434,14 +491,16 @@ void animation()
{
if (DOUBLE_MOVIE)
{
draw_wave(phi, psi, xy_in, scale, i, PLOT);
// draw_wave(phi, psi, xy_in, scale, i, PLOT);
draw_wave_e(phi, psi, total_energy, xy_in, scale, NSTEPS, PLOT);
draw_billiard();
glutSwapBuffers();
}
for (i=0; i<MID_FRAMES; i++) save_frame();
if (DOUBLE_MOVIE)
{
draw_wave(phi, psi, xy_in, scale, i, PLOT_B);
// draw_wave(phi, psi, xy_in, scale, i, PLOT_B);
draw_wave_e(phi, psi, total_energy, xy_in, scale, NSTEPS, PLOT_B);
draw_billiard();
glutSwapBuffers();
// for (i=0; i<END_FRAMES; i++) save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter + i);
@ -456,8 +515,16 @@ void animation()
free(psi[i]);
free(phi_tmp[i]);
free(psi_tmp[i]);
free(total_energy[i]);
free(xy_in[i]);
}
if (SAVE_TIME_SERIES)
{
fclose(time_series_left);
fclose(time_series_right);
}
}

84
wave_common.c
View File

@ -94,6 +94,28 @@ void init_wave_plus(double x, double y, double *phi[NX], double *psi[NX], short
}
}
void init_circular_wave_xplusminus(double xleft, double yleft, double xright, double yright,
double *phi[NX], double *psi[NX], short int * xy_in[NX])
/* initialise field with two drops, x > 0 and x < 0 do not communicate - phi is wave height, psi is phi at time t-1 */
{
int i, j;
double xy[2], dist2;
for (i=0; i<NX; i++)
for (j=0; j<NY; j++)
{
ij_to_xy(i, j, xy);
if (xy[0] < 0.0) dist2 = (xy[0]-xleft)*(xy[0]-xleft) + (xy[1]-yleft)*(xy[1]-yleft);
else dist2 = (xy[0]-xright)*(xy[0]-xright) + (xy[1]-yright)*(xy[1]-yright);
xy_in[i][j] = xy_in_billiard(xy[0],xy[1]);
if ((xy_in[i][j])||(TWOSPEEDS))
phi[i][j] = INITIAL_AMP*exp(-dist2/INITIAL_VARIANCE)*cos(-sqrt(dist2)/INITIAL_WAVELENGTH);
else phi[i][j] = 0.0;
psi[i][j] = 0.0;
}
}
void init_planar_wave(double x, double y, double *phi[NX], double *psi[NX], short int * xy_in[NX])
/* initialise field with drop at (x,y) - phi is wave height, psi is phi at time t-1 */
/* beta version, works for vertical planar wave only so far */
@ -339,7 +361,10 @@ void draw_wave_e(double *phi[NX], double *psi[NX], double *total_energy[NX], sho
}
case (P_ENERGY):
{
color_scheme(COLOR_SCHEME, compute_energy(phi, psi, xy_in, i, j), scale, time, rgb);
energy = compute_energy(phi, psi, xy_in, i, j);
/* adjust energy to color palette */
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):
@ -352,7 +377,9 @@ void draw_wave_e(double *phi[NX], double *psi[NX], double *total_energy[NX], sho
{
energy = compute_energy(phi, psi, xy_in, i, j);
total_energy[i][j] += energy;
color_scheme(COLOR_SCHEME, total_energy[i][j]/(double)(time+1), scale, time, rgb);
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;
}
}
@ -369,6 +396,59 @@ void draw_wave_e(double *phi[NX], double *psi[NX], double *total_energy[NX], sho
}
void draw_color_scheme(double x1, double y1, double x2, double y2, int plot, double min, double max)
{
int j, ij_botleft[2], ij_topright[2], imin, imax, jmin, jmax;
double y, dy, dy_e, rgb[3], value;
xy_to_ij(x1, y1, ij_botleft);
xy_to_ij(x2, y2, ij_topright);
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));
for (j = jmin; j < jmax; j++)
{
switch (plot) {
case (P_AMPLITUDE):
{
value = min + 1.0*dy*(double)(j - jmin);
color_scheme(COLOR_SCHEME, value, 1.0, 1, rgb);
break;
}
case (P_ENERGY):
{
value = dy_e*(double)(j - jmin);
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_MEAN_ENERGY):
{
value = dy_e*(double)(j - jmin);
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;
}
}
glColor3f(rgb[0], rgb[1], rgb[2]);
glVertex2i(imin, j);
glVertex2i(imax, j);
glVertex2i(imax, j+1);
glVertex2i(imin, j+1);
}
glEnd ();
glColor3f(1.0, 1.0, 1.0);
draw_rectangle(x1, y1, x2, y2);
}

23
wave_comparison.c
View File

@ -81,6 +81,17 @@
#define NGRIDX 20 /* 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
#define X_TARGET 0.4
#define Y_TARGET 0.7 /* shooter and target positions in laser fight */
#define ISO_XSHIFT_LEFT -1.65
#define ISO_XSHIFT_RIGHT 0.4
#define ISO_YSHIFT_LEFT -0.05
#define ISO_YSHIFT_RIGHT -0.05
#define ISO_SCALE 0.85 /* coordinates for isospectral billiards */
/* You can add more billiard tables by adapting the functions */
/* xy_in_billiard and draw_billiard below */
@ -133,9 +144,11 @@
/* Parameters of initial condition */
#define INITIAL_AMP 0.2 /* amplitude of initial condition */
#define INITIAL_VARIANCE 0.002 /* variance of initial condition */
#define INITIAL_WAVELENGTH 0.1 /* wavelength 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.02 /* wavelength of initial condition */
/* Plot type, see list in global_pdes.c */
@ -143,7 +156,7 @@
/* Color schemes */
#define COLOR_PALETTE 0 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE 0 /* Color palette, see list in global_pdes.c */
#define BLACK 1 /* background */
@ -166,8 +179,6 @@
#define VMAX 5.0 /* max value of wave amplitude */
#include "hsluv.c"
#include "global_pdes.c" /* constants and global variables */
#include "sub_wave.c" /* common functions for wave_billiard, heat and schrodinger */
#include "wave_common.c" /* common functions for wave_billiard, wave_comparison, etc */

43
wave_energy.c
View File

@ -50,23 +50,23 @@
#define NX 1280 /* number of grid points on x axis */
#define NY 720 /* number of grid points on y axis */
#define XMIN -1.777777778
#define XMAX 1.777777778 /* x interval */
#define YMIN -1.0
#define YMAX 1.0 /* 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 XMIN -1.777777778
// #define XMAX 1.777777778 /* x interval */
// #define YMIN -1.0
// #define YMAX 1.0 /* 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 JULIA_SCALE 1.0 /* scaling for Julia sets */
/* Choice of the billiard table */
#define B_DOMAIN 15 /* choice of domain shape, see list in global_pdes.c */
#define B_DOMAIN_B 15 /* choice of domain shape, see list in global_pdes.c */
#define B_DOMAIN 20 /* choice of domain shape, see list in global_pdes.c */
#define B_DOMAIN_B 20 /* choice of domain shape, see list in global_pdes.c */
#define CIRCLE_PATTERN 2 /* pattern of circles, see list in global_pdes.c */
#define CIRCLE_PATTERN 1 /* pattern of circles, see list in global_pdes.c */
#define CIRCLE_PATTERN_B 11 /* pattern of circles, see list in global_pdes.c */
#define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */
@ -85,6 +85,17 @@
#define NGRIDX 15 /* 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
#define X_TARGET 0.4
#define Y_TARGET 0.7 /* shooter and target positions in laser fight */
#define ISO_XSHIFT_LEFT -1.65
#define ISO_XSHIFT_RIGHT 0.4
#define ISO_YSHIFT_LEFT -0.05
#define ISO_YSHIFT_RIGHT -0.05
#define ISO_SCALE 0.85 /* coordinates for isospectral billiards */
/* You can add more billiard tables by adapting the functions */
/* xy_in_billiard and draw_billiard below */
@ -98,14 +109,8 @@
#define AMPLITUDE 0.025 /* amplitude of periodic excitation */
#define COURANT 0.02 /* Courant number */
#define COURANTB 0.004 /* Courant number in medium B */
// #define COURANTB 0.005 /* Courant number in medium B */
// #define COURANTB 0.008 /* Courant number in medium B */
#define GAMMA 0.0 /* damping factor in wave equation */
// #define GAMMA 1.0e-8 /* damping factor in wave equation */
#define GAMMAB 1.0e-8 /* damping factor in wave equation */
// #define GAMMAB 1.0e-6 /* damping factor in wave equation */
// #define GAMMAB 2.0e-4 /* damping factor in wave equation */
// #define GAMMAB 2.5e-4 /* damping factor in wave equation */
#define GAMMA_SIDES 1.0e-4 /* damping factor on boundary */
#define GAMMA_TOPBOT 1.0e-6 /* damping factor on boundary */
#define KAPPA 0.0 /* "elasticity" term enforcing oscillations */
@ -122,7 +127,7 @@
/* Parameters for length and speed of simulation */
#define NSTEPS 4500 /* number of frames of movie */
#define NSTEPS 3750 /* 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 200 /* time after which to start saving frames */
@ -170,8 +175,6 @@
#define VMAX 5.0 /* max value of wave amplitude */
#include "hsluv.c"
#include "global_pdes.c" /* constants and global variables */
#include "sub_wave.c" /* common functions for wave_billiard, heat and schrodinger */
#include "wave_common.c" /* common functions for wave_billiard, wave_comparison, etc */