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:
Nils Berglund 2022-11-20 23:17:39 +01:00 committed by GitHub
parent 30fc3e5a8b
commit cb2bcdb13d
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
21 changed files with 1625 additions and 288 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
global_ljones.c
View File

@ -64,6 +64,8 @@
#define S_MAZE 15 /* segments forming a maze */
#define S_EXT_RECTANGLE 16 /* particles outside a rectangle */
#define S_DAM_BRICKS 17 /* dam made of several bricks */
#define S_HLINE_HOLE 18 /* horizontal line with a hole in the bottom */
#define S_EXT_CIRCLE_RECT 19 /* particles outside a circle and a rectangle */
/* particle interaction */

1
global_particles.c
View File

@ -79,6 +79,7 @@ double x_shooter = -0.2, y_shooter = -0.6, x_target = 0.4, y_target = 0.7;
#define P_TREE 7 /* pine tree */
#define P_TOKA_NONSELF 8 /* Tokarsky non-self-unilluminable room */
#define P_MAZE 10 /* maze */
#define P_MAZE_DIAG 11 /* maze with 45 degrees angles */
/* Color palettes */

8
global_pdes.c
View File

@ -14,6 +14,9 @@
#define D_NOTHING 999 /* no boundaries */
#define D_RECTANGLE 0 /* rectangular domain */
#define D_ELLIPSE 1 /* elliptical domain */
#define D_EXT_ELLIPSE 199 /* exterior of elliptical domain */
#define D_EXT_ELLIPSE_CURVED 198 /* exterior of curved elliptical domain */
#define D_EXT_ELLIPSE_CURVED_BDRY 197 /* exterior of curved elliptical domain, with horizontal boundaries */
#define D_STADIUM 2 /* stadium-shaped domain */
#define D_SINAI 3 /* Sinai billiard */
#define D_DIAMOND 4 /* diamond-shaped billiard */
@ -65,6 +68,9 @@
#define D_WAVEGUIDE_W 52 /* W-shaped wave guide */
#define D_MAZE 53 /* maze */
#define D_MAZE_CLOSED 54 /* closed maze */
#define D_CHESSBOARD 55 /* chess board configuration */
#define D_TRIANGLE_TILES 56 /* triangular tiling */
#define D_HEX_TILES 57 /* honeycomb tiling */
#define NMAXCIRCLES 10000 /* total number of circles/polygons (must be at least NCX*NCY for square grid) */
#define NMAXPOLY 50000 /* maximal number of vertices of polygonal lines (for von Koch et al) */
@ -140,6 +146,8 @@
#define P_MEAN_ENERGY 3 /* energy averaged over time */
#define P_LOG_ENERGY 4 /* log of energy averaged over time */
#define P_LOG_MEAN_ENERGY 5 /* log of energy averaged over time */
#define P_ENERGY_FLUX 6 /* energy flux */
#define P_TOTAL_ENERGY_FLUX 7 /* energy flux averaged over time */
/* For Schrodinger equation */
#define P_MODULE 10 /* plot module of wave function squared */

1
heat.c
View File

@ -36,6 +36,7 @@
#include <omp.h>
#define MOVIE 0 /* set to 1 to generate movie */
#define SAVE_MEMORY 0 /* set to 1 to save memory when writing tiff images */
/* General geometrical parameters */

115
lennardjones.c
View File

@ -37,42 +37,32 @@
#include <time.h>
#define MOVIE 0 /* set to 1 to generate movie */
#define DOUBLE_MOVIE 0 /* set to 1 to produce movies for wave height and energy simultaneously */
#define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */
#define TIME_LAPSE 1 /* set to 1 to add a time-lapse movie at the end */
#define TIME_LAPSE 0 /* set to 1 to add a time-lapse movie at the end */
/* so far incompatible with double movie */
#define TIME_LAPSE_FACTOR 3 /* factor of time-lapse movie */
#define TIME_LAPSE_FIRST 1 /* set to 1 to show time-lapse version first */
#define TIME_LAPSE_FIRST 0 /* set to 1 to show time-lapse version first */
/* General geometrical parameters */
#define WINWIDTH 1280 /* window width */
#define WINHEIGHT 720 /* window height */
// #define XMIN -2.3
// #define XMAX 3.7 /* x interval */
// #define YMIN -1.6875
// #define YMAX 1.6875 /* 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 -3.3
#define XMAX 4.7 /* x interval */
#define YMIN -2.25
#define YMAX 2.25 /* y interval for 9/16 aspect ratio */
#define INITXMIN -1.95
#define INITXMAX 1.95 /* x interval for initial condition */
#define INITYMIN -1.1
#define INITYMAX 0.15 /* y interval for initial condition */
#define INITXMIN -2.5
#define INITXMAX 2.5 /* x interval for initial condition */
#define INITYMIN -1.7
#define INITYMAX 0.7 /* y interval for initial condition */
// #define BCXMIN -3.1
// #define BCXMAX 3.1 /* x interval for boundary condition */
// #define BCYMIN -4.5
// #define BCYMAX 4.5 /* y interval for boundary condition */
#define BCXMIN -5.1
#define BCXMAX 6.1 /* x interval for boundary condition */
#define BCYMIN -6.5
#define BCYMAX 44.0 /* y interval for boundary condition */
#define BCXMIN -2.0
#define BCXMAX 2.0 /* x interval for boundary condition */
#define BCYMIN -1.125
#define BCYMAX 1.125 /* y interval for boundary condition */
#define OBSXMIN -2.0
#define OBSXMAX 2.0 /* x interval for motion of obstacle */
@ -83,7 +73,7 @@
#define OBSTACLE_PATTERN 3 /* pattern of obstacles, see list in global_ljones.c */
#define ADD_FIXED_SEGMENTS 1 /* set to 1 to add fixed segments as obstacles */
#define SEGMENT_PATTERN 102 /* pattern of repelling segments, see list in global_ljones.c */
#define SEGMENT_PATTERN 19 /* pattern of repelling segments, see list in global_ljones.c */
#define ROCKET_SHAPE 2 /* shape of rocket combustion chamber, see list in global_ljones.c */
#define ROCKET_SHAPE_B 2 /* shape of second rocket */
#define NOZZLE_SHAPE 1 /* shape of nozzle, see list in global_ljones.c */
@ -103,18 +93,16 @@
#define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */
#define NPOISSON 100 /* number of points for Poisson C_RAND_POISSON arrangement */
#define PDISC_DISTANCE 2.7 /* minimal distance in Poisson disc process, controls density of particles */
#define PDISC_DISTANCE 3.0 /* minimal distance in Poisson disc process, controls density of particles */
#define PDISC_CANDIDATES 100 /* number of candidates in construction of Poisson disc process */
#define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */
#define LAMBDA 0.8 /* parameter controlling the dimensions of domain */
// #define MU 0.02 /* parameter controlling radius of particles */
// #define MU 0.015 /* parameter controlling radius of particles */
#define LAMBDA 0.3 /* parameter controlling the dimensions of domain */
#define MU 0.009 /* parameter controlling radius of particles */
// #define MU 0.012 /* parameter controlling radius of particles */
#define MU_B 0.018 /* parameter controlling radius of particles of second type */
#define NPOLY 25 /* number of sides of polygon */
#define APOLY 0.666666666 /* angle by which to turn polygon, in units of Pi/2 */
// #define NPOLY 4 /* number of sides of polygon */
#define NPOLY 3 /* number of sides of polygon */
#define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */
#define MDEPTH 4 /* depth of computation of Menger gasket */
#define MRATIO 3 /* ratio defining Menger gasket */
#define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */
@ -134,13 +122,11 @@
/* Parameters for length and speed of simulation */
#define NSTEPS 3300 /* number of frames of movie */
// #define NSTEPS 2000 /* number of frames of movie */
#define NSTEPS 1500 /* number of frames of movie */
#define NVID 100 /* number of iterations between images displayed on screen */
#define NSEG 250 /* number of segments of boundary */
#define INITIAL_TIME 10 /* time after which to start saving frames */
// #define OBSTACLE_INITIAL_TIME 10 /* time after which to start moving obstacle */
#define OBSTACLE_INITIAL_TIME 200 /* time after which to start moving obstacle */
#define INITIAL_TIME 250 /* time after which to start saving frames */
#define OBSTACLE_INITIAL_TIME 350 /* time after which to start moving obstacle */
#define BOUNDARY_WIDTH 1 /* width of particle boundary */
#define LINK_WIDTH 2 /* width of links between particles */
#define CONTAINER_WIDTH 4 /* width of container boundary */
@ -154,7 +140,7 @@
/* Boundary conditions, see list in global_ljones.c */
#define BOUNDARY_COND 20
#define BOUNDARY_COND 0
/* Plot type, see list in global_ljones.c */
@ -164,6 +150,7 @@
#define DRAW_BONDS 0 /* set to 1 to draw bonds between neighbours */
#define COLOR_BONDS 1 /* set to 1 to color bonds according to length */
#define FILL_TRIANGLES 1 /* set to 1 to fill triangles between neighbours */
#define ALTITUDE_LINES 0 /* set to 1 to add horizontal lines to show altitude */
#define COLOR_SEG_GROUPS 1 /* set to 1 to collor segment groups differently */
/* Color schemes */
@ -191,7 +178,7 @@
#define ENERGY_HUE_MAX 50.0 /* color of saturated particle */
#define PARTICLE_HUE_MIN 359.0 /* color of original particle */
#define PARTICLE_HUE_MAX 0.0 /* color of saturated particle */
#define PARTICLE_EMAX 2.0e2 /* energy of particle with hottest color */
#define PARTICLE_EMAX 3.0e2 /* energy of particle with hottest color */
#define HUE_TYPE0 70.0 /* hue of particles of type 0 */
#define HUE_TYPE1 280.0 /* hue of particles of type 1 */
#define HUE_TYPE2 .0 /* hue of particles of type 2 */
@ -200,10 +187,11 @@
#define RANDOM_RADIUS 0 /* set to 1 for random circle radius */
#define DT_PARTICLE 3.0e-6 /* time step for particle displacement */
#define KREPEL 12.0 /* constant in repelling force between particles */
#define EQUILIBRIUM_DIST 4.5 /* Lennard-Jones equilibrium distance */
#define EQUILIBRIUM_DIST 4.0 /* Lennard-Jones equilibrium distance */
#define EQUILIBRIUM_DIST_B 3.5 /* Lennard-Jones equilibrium distance for second type of particle */
#define REPEL_RADIUS 20.0 /* radius in which repelling force acts (in units of particle radius) */
#define DAMPING 5.0 /* damping coefficient of particles */
#define DAMPING 10.0 /* damping coefficient of particles */
#define INITIAL_DAMPING 35.0 /* damping coefficient of particles during initial phase */
#define PARTICLE_MASS 1.0 /* mass of particle of radius MU */
#define PARTICLE_MASS_B 1.0 /* mass of particle of radius MU */
#define PARTICLE_INERTIA_MOMENT 0.2 /* moment of inertia of particle */
@ -211,7 +199,7 @@
#define V_INITIAL 0.0 /* initial velocity range */
#define OMEGA_INITIAL 10.0 /* initial angular velocity range */
#define THERMOSTAT 1 /* set to 1 to switch on thermostat */
#define THERMOSTAT 0 /* set to 1 to switch on thermostat */
#define VARY_THERMOSTAT 0 /* set to 1 for time-dependent thermostat schedule */
#define SIGMA 5.0 /* noise intensity in thermostat */
#define BETA 0.02 /* initial inverse temperature */
@ -219,7 +207,7 @@
#define KSPRING_BOUNDARY 1.0e7 /* confining harmonic potential outside simulation region */
#define KSPRING_OBSTACLE 1.0e11 /* harmonic potential of obstacles */
#define NBH_DIST_FACTOR 7.5 /* radius in which to count neighbours */
#define GRAVITY 15.0 /* gravity acting on all particles */
#define GRAVITY 1500.0 /* gravity acting on all particles */
#define GRAVITY_X 0.0 /* horizontal gravity acting on all particles */
#define INCREASE_GRAVITY 0 /* set to 1 to increase gravity during the simulation */
#define GRAVITY_SCHEDULE 2 /* type of gravity schedule, see list in global_ljones.c */
@ -240,7 +228,7 @@
#define SPIN_RANGE_B 5.0 /* range of spin-spin interaction for second type of particle */
#define QUADRUPOLE_RATIO 0.6 /* anisotropy in quadrupole potential */
#define INCREASE_BETA 1 /* set to 1 to increase BETA during simulation */
#define INCREASE_BETA 0 /* set to 1 to increase BETA during simulation */
#define BETA_FACTOR 0.025 /* factor by which to change BETA during simulation */
#define N_TOSCILLATIONS 1.5 /* number of temperature oscillations in BETA schedule */
#define NO_OSCILLATION 1 /* set to 1 to have exponential BETA change only */
@ -301,29 +289,29 @@
#define OMEGAMAX 100.0 /* maximal rotation speed */
#define PRINT_OMEGA 0 /* set to 1 to print angular speed */
#define PRINT_PARTICLE_SPEEDS 0 /* set to 1 to print average speeds/momenta of particles */
#define PRINT_SEGMENTS_SPEEDS 1 /* set to 1 to print velocity of moving segments */
#define PRINT_SEGMENTS_SPEEDS 0 /* set to 1 to print velocity of moving segments */
#define MOVE_BOUNDARY 0 /* set to 1 to move repelling segments, due to force from particles */
#define SEGMENTS_MASS 40.0 /* mass of collection of segments */
#define DEACTIVATE_SEGMENT 1 /* set to 1 to deactivate last segment after a certain time */
#define SEGMENT_DEACTIVATION_TIME 200 /* time at which to deactivate last segment */
#define RELEASE_ROCKET_AT_DEACTIVATION 1 /* set to 1 to limit segments velocity before segment release */
#define SEGMENTS_X0 1.5 /* initial position of segments */
#define SEGMENTS_Y0 0.0 /* initial position of segments */
#define SEGMENTS_X0 1.0 /* initial position of segments */
#define SEGMENTS_Y0 0.8 /* initial position of segments */
#define SEGMENTS_VX0 0.0 /* initial velocity of segments */
#define SEGMENTS_VY0 0.0 /* initial velocity of segments */
#define DAMP_SEGS_AT_NEGATIVE_Y 0 /* set to 1 to dampen segments when y coordinate is negative */
#define MOVE_SEGMENT_GROUPS 1 /* set to 1 to group segments into moving units */
#define SEGMENT_GROUP_MASS 1000.0 /* mass of segment group */
#define SEGMENT_GROUP_I 1000.0 /* moment of inertia of segment group */
#define SEGMENT_GROUP_MASS 750.0 /* mass of segment group */
#define SEGMENT_GROUP_I 500.0 /* moment of inertia of segment group */
#define SEGMENT_GROUP_DAMPING 0.0 /* damping of segment groups */
#define GROUP_REPULSION 1 /* set to 1 for groups of segments to repel each other */
#define KSPRING_GROUPS 1.0e11 /* harmonic potential between segment groups */
#define GROUP_WIDTH 0.05 /* interaction width of groups */
#define GROUP_G_REPEL 1 /* set to 1 to add repulsion between centers of mass of groups */
#define GROUP_G_REPEL 0 /* set to 1 to add repulsion between centers of mass of groups */
#define GROUP_G_REPEL_RADIUS 1.2 /* radius within which centers of mass of groups repel each other */
#define TRACK_SEGMENT_GROUPS 1 /* set to 1 for view to track group of segments */
#define TRACK_SEGMENT_GROUPS 0 /* set to 1 for view to track group of segments */
#define TRACK_X_PADDING 2.0 /* distance from x boundary where tracking starts */
#define POSITION_DEPENDENT_TYPE 0 /* set to 1 to make particle type depend on initial position */
@ -336,8 +324,8 @@
#define PRINT_PARTICLE_NUMBER 0 /* set to 1 to print total number of particles */
#define PRINT_LEFT 0 /* set to 1 to print certain parameters at the top left instead of right */
#define PLOT_SPEEDS 1 /* set to 1 to add a plot of obstacle speeds (e.g. for rockets) */
#define PLOT_TRAJECTORIES 1 /* set to 1 to add a plot of obstacle trajectories (e.g. for rockets) */
#define PLOT_SPEEDS 0 /* set to 1 to add a plot of obstacle speeds (e.g. for rockets) */
#define PLOT_TRAJECTORIES 0 /* set to 1 to add a plot of obstacle trajectories (e.g. for rockets) */
#define VMAX_PLOT_SPEEDS 0.6 /* vertical scale of plot of obstacle speeds */
#define EHRENFEST_COPY 0 /* set to 1 to add equal number of larger particles (for Ehrenfest model) */
@ -362,7 +350,7 @@
#define PMAX 1000.0 /* maximal force */
#define HASHX 120 /* size of hashgrid in x direction */
#define HASHY 450 /* size of hashgrid in y direction */
#define HASHY 60 /* size of hashgrid in y direction */
#define HASHMAX 100 /* maximal number of particles per hashgrid cell */
#define HASHGRID_PADDING 0.1 /* padding of hashgrid outside simulation window */
@ -636,12 +624,15 @@ int thermostat_schedule(int i)
double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HASHX*HASHY],
double qx[NMAXCIRCLES], double qy[NMAXCIRCLES], double qangle[NMAXCIRCLES],
double px[NMAXCIRCLES], double py[NMAXCIRCLES], double pangle[NMAXCIRCLES],
double beta, int *nactive, int *nsuccess, int *nmove)
double beta, int *nactive, int *nsuccess, int *nmove, int initial_phase)
{
double a, totalenergy = 0.0;
double a, totalenergy = 0.0, damping;
static double b = 0.25*SIGMA*SIGMA*DT_PARTICLE/MU_XI, xi = 0.0;
int j, move;
if (initial_phase) damping = INITIAL_DAMPING;
else damping = DAMPING;
#pragma omp parallel for private(j,xi,totalenergy,a,move)
for (j=0; j<ncircles; j++) if (particle[j].active)
{
@ -696,10 +687,11 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA
}
else
{
px[j] *= exp(- DT_PARTICLE*DAMPING);
py[j] *= exp(- DT_PARTICLE*DAMPING);
px[j] *= exp(- DT_PARTICLE*damping);
py[j] *= exp(- DT_PARTICLE*damping);
pangle[j] *= exp(- DT_PARTICLE*damping);
}
if ((COUPLE_ANGLE_TO_THERMOSTAT)&&(particle[j].thermostat))
if ((THERMOSTAT_ON)&&(COUPLE_ANGLE_TO_THERMOSTAT)&&(particle[j].thermostat))
pangle[j] *= exp(- 0.5*DT_PARTICLE*xi);
particle[j].xc = qx[j] + 0.5*DT_PARTICLE*px[j]*particle[j].mass_inv;
@ -888,6 +880,7 @@ void evolve_segment_groups(t_segment segment[NMAXSEGMENTS], int time, t_group_se
maxdepth = 0.5*GROUP_WIDTH;
saturation_depth = 0.1*GROUP_WIDTH;
padding = 0.1;
for (group=0; group<ngroups; group++)
{
@ -1206,7 +1199,7 @@ void animation()
}
/* timestep of thermostat algorithm */
totalenergy = evolve_particles(particle, hashgrid, qx, qy, qangle, px, py, pangle, beta, &nactive, &nsuccess, &nmove);
totalenergy = evolve_particles(particle, hashgrid, qx, qy, qangle, px, py, pangle, beta, &nactive, &nsuccess, &nmove, i < INITIAL_TIME);
/* evolution of lid coordinate */
if (BOUNDARY_COND == BC_RECTANGLE_LID) evolve_lid(fboundary);

2
mangrove.c
View File

@ -43,6 +43,7 @@
#include <omp.h>
#define MOVIE 0 /* set to 1 to generate movie */
#define SAVE_MEMORY 0 /* set to 1 to save memory when writing tiff images */
/* General geometrical parameters */
@ -179,6 +180,7 @@
#define E_SCALE 2500.0 /* scaling factor for energy representation */
#define LOG_SCALE 1.0 /* scaling factor for energy log representation */
#define LOG_SHIFT 0.0 /* shift of colors on log scale */
#define FLUX_SCALE 1.0e4 /* scaling factor for enegy flux represtnation */
#define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */
#define COLORHUE 260 /* initial hue of water color for scheme C_LUM */

59
particle_billiard.c
View File

@ -31,7 +31,7 @@
#include <tiffio.h> /* Sam Leffler's libtiff library. */
#include <time.h>
#define MOVIE 0 /* set to 1 to generate movie */
#define MOVIE 1 /* set to 1 to generate movie */
#define WINWIDTH 1280 /* window width */
#define WINHEIGHT 720 /* window height */
@ -57,9 +57,9 @@
#define B_DOMAIN 30 /* choice of domain shape */
#define CIRCLE_PATTERN 1 /* pattern of circles */
#define POLYLINE_PATTERN 10 /* pattern of polyline */
#define POLYLINE_PATTERN 11 /* pattern of polyline */
#define ABSORBING_CIRCLES 1 /* set to 1 for circular scatterers to be absorbing */
#define ABSORBING_CIRCLES 0 /* set to 1 for circular scatterers to be absorbing */
#define NMAXCIRCLES 100000 /* total number of circles (must be at least NCX*NCY for square grid) */
#define NMAXPOLY 100000 /* total number of sides of polygonal line */
@ -87,7 +87,7 @@
/* Simulation parameters */
#define NPART 1000 /* number of particles */
#define NPART 10000 /* number of particles */
// #define NPART 2000 /* number of particles */
#define NPARTMAX 100000 /* maximal number of particles after resampling */
#define LMAX 0.01 /* minimal segment length triggering resampling */
@ -100,9 +100,9 @@
#define PRINT_COLLISION_NUMBER 0 /* set to 1 to print number of collisions */
#define TEST_ACTIVE 1 /* set to 1 to test whether particle is in billiard */
#define TEST_INITIAL_COND 1 /* set to 1 to allow only initial conditions that pass a test */
#define TEST_INITIAL_COND 0 /* set to 1 to allow only initial conditions that pass a test */
#define NSTEPS 6000 /* number of frames of movie */
#define NSTEPS 4500 /* number of frames of movie */
#define TIME 1500 /* time between movie frames, for fluidity of real-time simulation */
#define DPHI 0.00002 /* integration step */
#define NVID 25 /* number of iterations between images displayed on screen */
@ -117,17 +117,17 @@
/* Colors and other graphical parameters */
#define COLOR_PALETTE 11 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE 10 /* Color palette, see list in global_pdes.c */
#define NCOLORS 16 /* number of colors */
#define NCOLORS 128 /* number of colors */
#define COLORSHIFT 0 /* hue of initial color */
#define RAINBOW_COLOR 1 /* set to 1 to use different colors for all particles */
#define RAINBOW_COLOR 0 /* 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.04 /* length of velocity vectors */
#define BILLIARD_WIDTH 3 /* width of billiard */
#define PARTICLE_WIDTH 3 /* width of particles */
#define BILLIARD_WIDTH 2 /* width of billiard */
#define PARTICLE_WIDTH 2 /* width of particles */
#define FRONT_WIDTH 3 /* width of wave front */
#define BLACK 1 /* set to 1 for black background */
@ -141,10 +141,12 @@
#define SLEEP1 1 /* initial sleeping time */
#define SLEEP2 1 /* final sleeping time */
#define NXMAZE 8 /* width of maze */
#define NYMAZE 8 /* height of maze */
#define NXMAZE 16 /* width of maze */
#define NYMAZE 16 /* height of maze */
// #define NXMAZE 10 /* width of maze */
// #define NYMAZE 10 /* height of maze */
#define MAZE_MAX_NGBH 4 /* max number of neighbours of maze cell */
#define RAND_SHIFT 58 /* seed of random number generator */
#define RAND_SHIFT 1 /* seed of random number generator */
#define MAZE_XSHIFT 0.0 /* horizontal shift of maze */
@ -660,16 +662,25 @@ void print_left_right_part_number(double *configs[NPARTMAX], int active[NPARTMAX
{
char message[50];
int i, nleft = 0, nright = 0;
double rgb[3], x1, cosphi;
double rgb[3], x1, y1, cosphi, sinphi;
static int maxnleft = 0, maxnright = 0;
/* count active particles, using the fact that absorbed particles have been given dummy coordinates */
for (i=0; i<nparticles; i++) /*if (active[i])*/
for (i=0; i<nparticles; i++) if (configs[i][0] >= DUMMY_ABSORBING)
{
cosphi = (configs[i][6] - configs[i][4])/configs[i][3];
sinphi = (configs[i][7] - configs[i][5])/configs[i][3];
x1 = configs[i][4] + configs[i][2]*cosphi;
y1 = configs[i][5] + configs[i][2]*cosphi;
if (x1 < xmin) nleft++;
else if (x1 > xmax) nright++;
else if (x1 > xmax)
{
nright++;
printf("Detected leaving particle %i at (%.2f, %2f)\n", i, x1, y1);
}
}
if (nleft > maxnleft) maxnleft = nleft;
if (nright > maxnright) maxnright = nright;
hsl_to_rgb(0.0, 0.0, 0.0, rgb);
@ -677,11 +688,11 @@ void print_left_right_part_number(double *configs[NPARTMAX], int active[NPARTMAX
erase_area(xr, yr - 0.03, 0.4, 0.12, rgb);
glColor3f(1.0, 1.0, 1.0);
if (nleft > 1) sprintf(message, "%i particles", nleft);
else sprintf(message, "%i particle", nleft);
if (nleft > 1) sprintf(message, "%i particles", maxnleft);
else sprintf(message, "%i particle", maxnleft);
write_text(xl, yl, message);
if (nright > 1) sprintf(message, "%i particles", nright);
else sprintf(message, "%i particle", nright);
if (nright > 1) sprintf(message, "%i particles", maxnright);
else sprintf(message, "%i particle", maxnright);
write_text(xr, yr, message);
}
@ -770,7 +781,7 @@ void animation()
// init_drop_config(-0.5, 0.0, 0.2, 0.4, configs);
init_drop_config(-1.3, 0.2, -0.25*PI, 0.25*PI, configs);
init_drop_config(-1.15, 0.01, 0.00*PI, 0.3*PI, configs);
// init_drop_config(-1.3, -0.1, 0.0, DPI, configs);
// init_drop_config(1.4, 0.1, 0.0, DPI, configs);
@ -792,7 +803,7 @@ void animation()
if (DRAW_BILLIARD) draw_billiard();
if (PRINT_PARTICLE_NUMBER) print_part_number(configs, active, XMIN + 0.1, YMIN + 0.1);
else if (PRINT_LEFT_RIGHT_PARTICLE_NUMBER)
print_left_right_part_number(configs, active, XMIN + 0.1, YMIN + 0.1, XMAX - 0.45, YMIN + 0.1, YMAX + MAZE_XSHIFT, YMAX + MAZE_XSHIFT);
print_left_right_part_number(configs, active, XMIN + 0.1, YMIN + 0.05, XMAX - 0.45, YMIN + 0.05, XMIN + MAZE_XSHIFT, XMAX + MAZE_XSHIFT);
else if (PRINT_COLLISION_NUMBER) print_collision_number(ncollisions, XMIN + 0.1, YMIN + 0.1);
glutSwapBuffers();
@ -850,7 +861,7 @@ void animation()
if (DRAW_BILLIARD) draw_billiard();
if (PRINT_PARTICLE_NUMBER) print_part_number(configs, active, XMIN + 0.1, YMIN + 0.1);
else if (PRINT_LEFT_RIGHT_PARTICLE_NUMBER)
print_left_right_part_number(configs, active, XMIN + 0.1, YMIN + 0.1, XMAX - 0.45, YMIN + 0.1, YMAX + MAZE_XSHIFT, YMAX + MAZE_XSHIFT);
print_left_right_part_number(configs, active, XMIN + 0.1, YMIN + 0.05, XMAX - 0.45, YMIN + 0.05, YMIN + MAZE_XSHIFT, YMAX + MAZE_XSHIFT);
// print_left_right_part_number(configs, XMIN + 0.1, YMIN + 0.1, XMAX - 0.45, YMIN + 0.1, YMIN + MAZE_XSHIFT, YMAX + MAZE_XSHIFT);
else if (PRINT_COLLISION_NUMBER) print_collision_number(ncollisions, XMIN + 0.1, YMIN + 0.1);

158
rde.c
View File

@ -41,6 +41,7 @@
#define MOVIE 0 /* set to 1 to generate movie */
#define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */
#define SAVE_MEMORY 0 /* set to 1 to save memory when writing tiff images */
/* General geometrical parameters */
@ -63,6 +64,8 @@
// // #define NY 180 /* 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 960 /* number of grid points on x axis */
// // #define NY 540 /* 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 */
@ -88,7 +91,7 @@
/* Choice of the billiard table */
#define B_DOMAIN 999 /* choice of domain shape, see list in global_pdes.c */
#define B_DOMAIN 197 /* choice of domain shape, see list in global_pdes.c */
#define CIRCLE_PATTERN 99 /* pattern of circles, see list in global_pdes.c */
@ -96,8 +99,8 @@
#define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */
#define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */
#define LAMBDA 1.0 /* parameter controlling the dimensions of domain */
#define MU 1.0 /* parameter controlling the dimensions of domain */
#define LAMBDA 0.7 /* parameter controlling the dimensions of domain */
#define MU 0.15 /* parameter controlling the dimensions of domain */
#define NPOLY 5 /* number of sides of polygon */
#define APOLY 2.0 /* angle by which to turn polygon, in units of Pi/2 */
#define MDEPTH 7 /* depth of computation of Menger gasket */
@ -150,7 +153,13 @@
#define SMOOTHEN_VORTICITY 1 /* set to 1 to smoothen vorticity field in Euler equation */
#define SMOOTHEN_PERIOD 10 /* period between smoothenings */
#define SMOOTH_FACTOR 0.015 /* factor by which to smoothen */
// #define SMOOTH_FACTOR 0.05 /* factor by which to smoothen */
#define SMOOTH_FACTOR 0.03 /* factor by which to smoothen */
// #define SMOOTH_FACTOR 0.015 /* factor by which to smoothen */
// #define SMOOTH_FACTOR 0.01 /* factor by which to smoothen */
#define ADD_TRACERS 1 /* set to 1 to add tracer particles (for Euler equations) */
#define N_TRACERS 1000 /* number of tracer particles */
#define T_OUT 2.0 /* outside temperature */
#define T_IN 0.0 /* inside temperature */
@ -179,12 +188,14 @@
#define B_COND 1
#define EULER_GRADIENT_YSHIFT 0.0 /* y-shift in computation of gradient in Euler equation */
/* Parameters for length and speed of simulation */
#define NSTEPS 4000 /* number of frames of movie */
// #define NSTEPS 2500 /* number of frames of movie */
#define NVID 50 /* number of iterations between images displayed on screen */
// #define NVID 100 /* number of iterations between images displayed on screen */
#define NSTEPS 2250 /* number of frames of movie */
// #define NSTEPS 500 /* number of frames of movie */
// #define NVID 90 /* number of iterations between images displayed on screen */
#define NVID 120 /* number of iterations between images displayed on screen */
// #define NVID 1100 /* number of iterations between images displayed on screen */
#define ACCELERATION_FACTOR 1.0 /* factor by which to increase NVID in course of simulation */
#define DT_ACCELERATION_FACTOR 1.0 /* factor by which to increase time step in course of simulation */
@ -208,6 +219,8 @@
#define ROTATE_VIEW 0 /* set to 1 to rotate position of observer */
#define ROTATE_ANGLE 360.0 /* total angle of rotation during simulation */
#define DRAW_PERIODICISED 1 /* set to 1 to repeat wave periodically in x and y directions */
/* Plot type - color scheme */
#define CPLOT 52
@ -253,8 +266,8 @@
/* Color schemes, see list in global_pdes.c */
#define COLOR_PALETTE 11 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 17 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE 14 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 13 /* Color palette, see list in global_pdes.c */
#define BLACK 1 /* black background */
@ -264,7 +277,7 @@
#define SCALE 0 /* set to 1 to adjust color scheme to variance of field */
#define SLOPE 1.0 /* sensitivity of color on wave amplitude */
#define VSCALE_AMPLITUDE 0.5 /* additional scaling factor for color scheme P_3D_AMPLITUDE */
#define VSCALE_AMPLITUDE 1.5 /* additional scaling factor for color scheme P_3D_AMPLITUDE */
#define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */
#define CURL_SCALE 0.000015 /* scaling factor for curl representation */
#define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */
@ -278,18 +291,20 @@
#define HUEMEAN 359.0 /* mean value of hue for color scheme C_HUE */
#define HUEAMP -359.0 /* amplitude of variation of hue for color scheme C_HUE */
#define E_SCALE 100.0 /* scaling factor for energy representation */
#define LOG_SCALE 0.75 /* scaling factor for energy log representation */
#define FLUX_SCALE 100.0 /* scaling factor for energy representation */
#define LOG_SCALE 0.5 /* scaling factor for energy log representation */
#define LOG_SHIFT 1.0
#define LOG_MIN 1.0e-3 /* floor value for log vorticity plot */
#define NXMAZE 7 /* width of maze */
#define NYMAZE 7 /* height of maze */
#define MAZE_MAX_NGBH 4 /* max number of neighbours of maze cell */
#define RAND_SHIFT 24 /* seed of random number generator */
#define RAND_SHIFT 0 /* seed of random number generator */
#define MAZE_XSHIFT 0.0 /* horizontal shift of maze */
#define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */
#define COLORBAR_RANGE 3.0 /* scale of color scheme bar */
#define COLORBAR_RANGE_B 2.5 /* scale of color scheme bar for 2nd part */
#define COLORBAR_RANGE_B 3.0 /* scale of color scheme bar for 2nd part */
#define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */
/* only for compatibility with wave_common.c */
@ -320,12 +335,12 @@ double light[3] = {0.816496581, -0.40824829, 0.40824829}; /* vector of "lig
double observer[3] = {8.0, 8.0, 8.0}; /* location of observer for REP_PROJ_3D representation */
int reset_view = 0; /* switch to reset 3D view parameters (for option ROTATE_VIEW) */
#define Z_SCALING_FACTOR 0.25 /* overall scaling factor of z axis for REP_PROJ_3D representation */
#define XY_SCALING_FACTOR 1.8 /* overall scaling factor for on-screen (x,y) coordinates after projection */
#define Z_SCALING_FACTOR 0.08 /* overall scaling factor of z axis for REP_PROJ_3D representation */
#define XY_SCALING_FACTOR 1.7 /* overall scaling factor for on-screen (x,y) coordinates after projection */
#define ZMAX_FACTOR 1.0 /* max value of z coordinate for REP_PROJ_3D representation */
#define XSHIFT_3D -0.1 /* overall x shift for REP_PROJ_3D representation */
#define YSHIFT_3D 0.0 /* overall y shift for REP_PROJ_3D representation */
#define BORDER_PADDING 2 /* distance from boundary at which to plot points, to avoid boundary effects due to gradient */
#define YSHIFT_3D 0.1 /* overall y shift for REP_PROJ_3D representation */
#define BORDER_PADDING 0 /* distance from boundary at which to plot points, to avoid boundary effects due to gradient */
/* For debugging purposes only */
#define FLOOR 1 /* set to 1 to limit wave amplitude to VMAX */
@ -514,8 +529,8 @@ void evolve_wave_half(double *phi_in[NFIELDS], double *phi_out[NFIELDS], short i
{
nabla_psi = (double *)malloc(2*NX*NY*sizeof(double));
nabla_omega = (double *)malloc(2*NX*NY*sizeof(double));
compute_gradient_euler(phi_in[0], nabla_psi);
compute_gradient_euler(phi_in[1], nabla_omega);
compute_gradient_euler(phi_in[0], nabla_psi, EULER_GRADIENT_YSHIFT);
compute_gradient_euler(phi_in[1], nabla_omega, 0.0);
dx = (XMAX-XMIN)/((double)NX);
if (SMOOTHEN_VORTICITY) /* beta: try to reduce formation of ripples */
@ -622,9 +637,15 @@ void evolve_wave_half(double *phi_in[NFIELDS], double *phi_out[NFIELDS], short i
}
case (E_EULER_INCOMP):
{
phi_out[0][i*NY+j] = phi_in[0][i*NY+j] + intstep*stiffness*(delta_phi[0][i*NY+j] + phi_in[1][i*NY+j]*dx*dx);
phi_out[1][i*NY+j] = phi_in[1][i*NY+j] - intstep*K_EULER*(nabla_omega[i*NY+j]*nabla_psi[NX*NY+i*NY+j]);
phi_out[1][i*NY+j] += intstep*K_EULER*(nabla_omega[NX*NY+i*NY+j]*nabla_psi[i*NY+j]);
// if ((j > 1)&&(j < NY - 1))
{
phi_out[0][i*NY+j] = phi_in[0][i*NY+j] + intstep*stiffness*(delta_phi[0][i*NY+j] + phi_in[1][i*NY+j]*dx*dx);
// phi_out[0][i*NY+j] += intstep*EULER_GRADIENT_YSHIFT;
phi_out[1][i*NY+j] = phi_in[1][i*NY+j] - intstep*K_EULER*(nabla_omega[i*NY+j]*nabla_psi[NX*NY+i*NY+j]);
phi_out[1][i*NY+j] += intstep*K_EULER*(nabla_omega[NX*NY+i*NY+j]*nabla_psi[i*NY+j]);
// if ((i == 0)&&(j%10 == 0)) printf("j = %i, psi = %.5lg\n", j, phi_out[0][i*NY+j]);
}
break;
}
}
@ -664,7 +685,7 @@ void evolve_wave_half(double *phi_in[NFIELDS], double *phi_out[NFIELDS], short i
}
}
void evolve_wave(double *phi[NFIELDS], double *phi_tmp[NFIELDS], short int xy_in[NX*NY], double potential_field[NX*NY], double vector_potential_field[2*NX*NY])
void evolve_wave(double *phi[NFIELDS], double *phi_tmp[NFIELDS], short int xy_in[NX*NY], double potential_field[NX*NY], double vector_potential_field[2*NX*NY])
/* time step of field evolution */
{
evolve_wave_half(phi, phi_tmp, xy_in, potential_field, vector_potential_field);
@ -672,6 +693,63 @@ void evolve_wave(double *phi[NFIELDS], double *phi_tmp[NFIELDS], short int xy_in
}
void evolve_tracers(double *phi[NFIELDS], double tracers[2*N_TRACERS*NSTEPS], int time, int nsteps, double step)
/* time steps of tracer particle evolution (for Euler equation) */
{
int tracer, i, j, t, ij[2], iplus, jplus;
double x, y, xy[2], vx, vy;
step = 0.2;
for (tracer = 0; tracer < N_TRACERS; tracer++)
{
x = tracers[time*2*N_TRACERS + 2*tracer];
y = tracers[time*2*N_TRACERS + 2*tracer + 1];
// printf("Tracer %i position (%.2f, %.2f)\n", tracer, x, y);
for (t=0; t<nsteps; t++)
{
xy_to_ij_safe(x, y, ij);
i = ij[0];
j = ij[1];
iplus = i + 1; if (iplus == NX) iplus = 0;
jplus = j + 1; if (jplus == NY) jplus = 0;
vx = phi[0][i*NY+jplus] - phi[0][i*NY+j];
vy = -(phi[0][iplus*NY+j] - phi[0][i*NY+j]);
if (j == 0) vx += EULER_GRADIENT_YSHIFT;
else if (j == NY-1) vx -= EULER_GRADIENT_YSHIFT;
// v = module2(vx, vy);
// if ((v > 0.0)&&(v < 0.1))
// {
// vx = vx*0.1/v;
// vy = vy*0.1/v;
// }
// printf("(i, j) = (%i, %i), Tracer %i velocity (%.6f, %.6f)\n", i, j, tracer, vx, vy);
x += vx*step;
y += vy*step;
}
// printf("Tracer %i velocity (%.2f, %.2f)\n", tracer, vx, vy);
if (x > XMAX) x += (XMIN - XMAX);
if (x < XMIN) x += (XMAX - XMIN);
if (y > YMAX) y += (YMIN - YMAX);
if (y < YMIN) y += (YMAX - YMIN);
if (time+1 < NSTEPS)
{
tracers[(time+1)*2*N_TRACERS + 2*tracer] = x;
tracers[(time+1)*2*N_TRACERS + 2*tracer + 1] = y;
}
}
}
void print_level(int level)
{
double pos[2];
@ -852,7 +930,7 @@ void animation()
{
double time = 0.0, scale, dx, var, jangle, cosj, sinj, sqrintstep,
intstep0, viscosity_printed, fade_value, noise = NOISE_INTENSITY;
double *phi[NFIELDS], *phi_tmp[NFIELDS], *potential_field, *vector_potential_field;
double *phi[NFIELDS], *phi_tmp[NFIELDS], *potential_field, *vector_potential_field, *tracers;
short int *xy_in;
int i, j, k, s, nvid, field;
static int counter = 0;
@ -886,6 +964,8 @@ void animation()
initialize_potential(potential_field);
initialize_vector_potential(vector_potential_field);
}
if (ADD_TRACERS) tracers = (double *)malloc(2*NSTEPS*N_TRACERS*sizeof(double));
dx = (XMAX-XMIN)/((double)NX);
intstep = DT/(dx*dx);
@ -908,10 +988,10 @@ void animation()
// init_fermion_state(-0.5, 0.5, 2.0, 0.0, 0.1, phi, xy_in);
// init_boson_state(-0.5, 0.5, 2.0, 0.0, 0.1, phi, xy_in);
// init_vortex_state(0.4, 0.0, 0.1, phi, xy_in);
// add_vortex_state(-0.4, 0.0, 0.1, phi, xy_in);
// init_shear_flow(1.0, 0.02, 0.15, 1, 1, phi, xy_in);
// init_laminar_flow(1.0, 0.1, 0.5, 0.0, phi, xy_in);
init_shear_flow(1.0, 0.02, 0.03, 1, 1, phi, xy_in);
init_shear_flow(-1.0, 0.0, 0.1, 1, 1, 0.0, phi, xy_in);
init_cfield_rde(phi, xy_in, CPLOT, rde, 0);
if (PLOT_3D) init_zfield_rde(phi, xy_in, ZPLOT, rde, 0);
@ -922,6 +1002,12 @@ void animation()
if (PLOT_3D) init_zfield_rde(phi, xy_in, ZPLOT_B, rde, 1);
}
if (ADD_TRACERS) for (i=0; i<N_TRACERS; i++)
{
tracers[2*i] = XMIN + 0.05 + (XMAX - XMIN - 0.1)*rand()/RAND_MAX;
tracers[2*i+1] = YMIN + 0.05 + (YMAX - YMIN - 0.1)*rand()/RAND_MAX;
}
blank();
glColor3f(0.0, 0.0, 0.0);
@ -981,6 +1067,15 @@ void animation()
printf("Evolving wave\n");
for (j=0; j<nvid; j++) evolve_wave(phi, phi_tmp, xy_in, potential_field, vector_potential_field);
if (ADD_TRACERS)
{
printf("Evolving tracer particles\n");
evolve_tracers(phi, tracers, i, 10, 0.1);
// for (j=0; j<N_TRACERS; j++)
// printf("Tracer %i position (%.2f, %.2f)\n", j, tracers[2*N_TRACERS*i + 2*j], tracers[2*N_TRACERS*i + 2*j + 1]);
draw_tracers(phi, tracers, i, 0, 1.0);
}
if (ANTISYMMETRIZE_WAVE_FCT) antisymmetrize_wave_function(phi, xy_in);
for (j=0; j<NFIELDS; j++) printf("field[%i] = %.3lg\t", j, phi[j][0]);
@ -1036,6 +1131,7 @@ void animation()
if ((i >= INITIAL_TIME)&&(DOUBLE_MOVIE))
{
draw_wave_rde(1, phi, xy_in, rde, potential_field, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 0, 1.0, REFRESH_B);
if (ADD_TRACERS) draw_tracers(phi, tracers, i, 0, 1.0);
// draw_billiard();
if (PRINT_PARAMETERS) print_parameters(rde, xy_in, time, 0, viscosity_printed, noise);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0);
@ -1062,6 +1158,7 @@ void animation()
if (DOUBLE_MOVIE)
{
draw_wave_rde(0, phi, xy_in, rde, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1);
if (ADD_TRACERS) draw_tracers(phi, tracers, NSTEPS, 0, 1.0);
// draw_billiard();
if (PRINT_PARAMETERS) print_parameters(rde, xy_in, time, 0, viscosity_printed, noise);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 0, 1.0);
@ -1072,6 +1169,7 @@ void animation()
{
fade_value = 1.0 - (double)i/(double)MID_FRAMES;
draw_wave_rde(0, phi, xy_in, rde, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 1, fade_value, 0);
if (ADD_TRACERS) draw_tracers(phi, tracers, NSTEPS, 1, fade_value);
// draw_billiard();
if (PRINT_PARAMETERS) print_parameters(rde, xy_in, time, 0, viscosity_printed, noise);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 1, fade_value);
@ -1079,6 +1177,7 @@ void animation()
save_frame_counter(NSTEPS + i + 1);
}
draw_wave_rde(1, phi, xy_in, rde, potential_field, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 0, 1.0, REFRESH_B);
if (ADD_TRACERS) draw_tracers(phi, tracers, NSTEPS, 0, 1.0);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0);
glutSwapBuffers();
@ -1087,6 +1186,7 @@ void animation()
{
fade_value = 1.0 - (double)i/(double)END_FRAMES;
draw_wave_rde(1, phi, xy_in, rde, potential_field, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 1, fade_value, 0);
if (ADD_TRACERS) draw_tracers(phi, tracers, NSTEPS, 1, fade_value);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 1, fade_value);
glutSwapBuffers();
save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter + i);
@ -1099,6 +1199,7 @@ void animation()
{
fade_value = 1.0 - (double)i/(double)END_FRAMES;
draw_wave_rde(0, phi, xy_in, rde, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 1, fade_value, 0);
if (ADD_TRACERS) draw_tracers(phi, tracers, NSTEPS, 1, fade_value);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 1, fade_value);
glutSwapBuffers();
save_frame_counter(NSTEPS + 1 + counter + i);
@ -1120,6 +1221,7 @@ void animation()
free(potential_field);
free(vector_potential_field);
}
if (ADD_TRACERS) free(tracers);
printf("Time %.5lg\n", time);

1
schrodinger.c
View File

@ -36,6 +36,7 @@
#include <omp.h>
#define MOVIE 0 /* set to 1 to generate movie */
#define SAVE_MEMORY 0 /* set to 1 to save memory when writing tiff images */
/* General geometrical parameters */

220
sub_lj.c
View File

@ -955,7 +955,7 @@ void init_obstacle_config(t_obstacle obstacle[NMAXOBSTACLES])
}
}
void add_rotated_angle_to_segments(double x1, double y1, double x2, double y2, double width, t_segment segment[NMAXSEGMENTS])
void add_rotated_angle_to_segments(double x1, double y1, double x2, double y2, double width, t_segment segment[NMAXSEGMENTS], int group)
/* add four segments forming a rectangle, specified by two adjacent corners and width */
{
double tx, ty, ux, uy, norm, x3, y3, x4, y4;
@ -993,7 +993,11 @@ void add_rotated_angle_to_segments(double x1, double y1, double x2, double y2, d
segment[n+3].x2 = x1;
segment[n+3].y2 = y1;
for (i=0; i<4; i++) segment[n+i].concave = 1;
for (i=0; i<4; i++)
{
segment[n+i].concave = 1;
segment[n+i].group = group;
}
nsegments += 4;
}
else printf("Warning: NMAXSEGMENTS too small\n");
@ -1050,6 +1054,38 @@ void add_rectangle_to_segments(double x1, double y1, double x2, double y2, t_seg
}
void add_circle_to_segments(double x, double y, double r, int nsegs, double angle0, t_segment segment[NMAXSEGMENTS], int group)
/* add segments forming a circle to linear obstacle configuration */
{
int i, n = nsegments, nplus, nminus;
double angle;
angle = DPI/(double)nsegs;
if (nsegments + nsegs < NMAXSEGMENTS)
{
for (i=0; i<nsegs; i++)
{
segment[n+i].x1 = -SEGMENTS_X0 + r*cos(((double)i)*angle + angle0*PID);
segment[n+i].y1 = SEGMENTS_Y0 - r*sin(((double)i)*angle + angle0*PID);
segment[n+i].x2 = -SEGMENTS_X0 + r*cos(((double)(i+1))*angle + angle0*PID);
segment[n+i].y2 = SEGMENTS_Y0 - r*sin(((double)(i+1))*angle + angle0*PID);
segment[n+i].angle1 = -((double)i + 0.5)*angle - angle0*PID;
segment[n+i].angle2 = -((double)i - 0.5)*angle - angle0*PID;
while (segment[n+i].angle1 < 0.0) segment[n+i].angle1 += DPI;
while (segment[n+i].angle2 < segment[n+i].angle1) segment[n+i].angle2 += DPI;
segment[n+i].concave = 1;
segment[n+i].group = group;
}
nsegments += nsegs;
}
else printf("Warning: NMAXSEGMENTS too small\n");
}
double nozzle_width(double x, double width, int nozzle_shape)
/* width of bell-shaped nozzle */
{
@ -1097,7 +1133,7 @@ void add_rocket_to_segments(t_segment segment[NMAXSEGMENTS], double x0, double y
angle = -PID + (double)(i+1)*DPI/(double)NPOLY;
x2 = x0 + 0.7*LAMBDA*cos(angle);
y2 = y0 + YMIN + LAMBDA*(1.7 + 0.7*sin(angle));
add_rotated_angle_to_segments(x1, y1, x2, y2, 0.02, segment);
add_rotated_angle_to_segments(x1, y1, x2, y2, 0.02, segment, 0);
}
break;
}
@ -1106,23 +1142,23 @@ void add_rocket_to_segments(t_segment segment[NMAXSEGMENTS], double x0, double y
/* dimensions chosen to have same area as circular chamber */
a = 0.5*LAMBDA;
b = 0.49*PI*LAMBDA;
add_rotated_angle_to_segments(x0+nozx, nozy, x0+a, nozy, 0.02, segment);
add_rotated_angle_to_segments(x0+a, nozy, x0+a, nozy+b, 0.02, segment);
add_rotated_angle_to_segments(x0+a, nozy+b, x0-a, nozy+b, 0.02, segment);
add_rotated_angle_to_segments(x0-a, nozy+b, x0-a, nozy, 0.02, segment);
add_rotated_angle_to_segments(x0-a, nozy, x0-nozx, nozy, 0.02, segment);
add_rotated_angle_to_segments(x0+nozx, nozy, x0+a, nozy, 0.02, segment, 0);
add_rotated_angle_to_segments(x0+a, nozy, x0+a, nozy+b, 0.02, segment, 0);
add_rotated_angle_to_segments(x0+a, nozy+b, x0-a, nozy+b, 0.02, segment, 0);
add_rotated_angle_to_segments(x0-a, nozy+b, x0-a, nozy, 0.02, segment, 0);
add_rotated_angle_to_segments(x0-a, nozy, x0-nozx, nozy, 0.02, segment, 0);
break;
}
case (RCK_RECT_HAT): /* rectangular chamber with a hat */
{
a = 0.5*LAMBDA;
b = 0.49*PI*LAMBDA;
add_rotated_angle_to_segments(x0+nozx, nozy, x0+a, nozy, 0.02, segment);
add_rotated_angle_to_segments(x0+a, nozy, x0+a, nozy+b, 0.02, segment);
add_rotated_angle_to_segments(x0+a, nozy+b, x0, nozy+b+a, 0.02, segment);
add_rotated_angle_to_segments(x0, nozy+b+a, x0-a, nozy+b, 0.02, segment);
add_rotated_angle_to_segments(x0-a, nozy+b, x0-a, nozy, 0.02, segment);
add_rotated_angle_to_segments(x0-a, nozy, x0-nozx, nozy, 0.02, segment);
add_rotated_angle_to_segments(x0+nozx, nozy, x0+a, nozy, 0.02, segment, 0);
add_rotated_angle_to_segments(x0+a, nozy, x0+a, nozy+b, 0.02, segment, 0);
add_rotated_angle_to_segments(x0+a, nozy+b, x0, nozy+b+a, 0.02, segment, 0);
add_rotated_angle_to_segments(x0, nozy+b+a, x0-a, nozy+b, 0.02, segment, 0);
add_rotated_angle_to_segments(x0-a, nozy+b, x0-a, nozy, 0.02, segment, 0);
add_rotated_angle_to_segments(x0-a, nozy, x0-nozx, nozy, 0.02, segment, 0);
break;
}
}
@ -1138,18 +1174,18 @@ void add_rocket_to_segments(t_segment segment[NMAXSEGMENTS], double x0, double y
x1 = x0 + nozzle_width(y1 - y0, nozx, nozzle_shape);
y2 = y1 + dx;
x2 = x0 + nozzle_width(y2 - y0, nozx, nozzle_shape);
add_rotated_angle_to_segments(x1, y1 + YMIN + LAMBDA, x2, y2 + YMIN + LAMBDA, 0.02, segment);
add_rotated_angle_to_segments(x1, y1 + YMIN + LAMBDA, x2, y2 + YMIN + LAMBDA, 0.02, segment, 0);
}
add_rotated_angle_to_segments(x2, y2 + YMIN + LAMBDA, x0 + nozx, nozy, 0.02, segment);
add_rotated_angle_to_segments(x2, y2 + YMIN + LAMBDA, x0 + nozx, nozy, 0.02, segment, 0);
for (i=0; i<nsides; i++)
{
y1 = y0 - LAMBDA + dx*(double)(i-1);
x1 = x0 - nozzle_width(y1 - y0, nozx, nozzle_shape);
y2 = y1 + dx;
x2 = x0 - nozzle_width(y2 - y0, nozx, nozzle_shape);
add_rotated_angle_to_segments(x1, y1 + YMIN + LAMBDA, x2, y2 + YMIN + LAMBDA, 0.02, segment);
add_rotated_angle_to_segments(x1, y1 + YMIN + LAMBDA, x2, y2 + YMIN + LAMBDA, 0.02, segment, 0);
}
add_rotated_angle_to_segments(x2, y2 + YMIN + LAMBDA, x0 - nozx, nozy, 0.02, segment);
add_rotated_angle_to_segments(x2, y2 + YMIN + LAMBDA, x0 - nozx, nozy, 0.02, segment, 0);
}
for (i=nsegments0; i<nsegments; i++) segment[i].inactivate = 0;
@ -1491,7 +1527,7 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
y1 = LAMBDA*sin(angle2);
x2 = 0.7*cos(angle2);
y2 = 0.7*sin(angle2);
add_rotated_angle_to_segments(x1, y1, x2, y2, MU, segment);
add_rotated_angle_to_segments(x1, y1, x2, y2, MU, segment, 0);
for (j=0; j<nsides; j++) if (j!=nsides/2)
{
@ -1509,7 +1545,7 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
x2 = r*cos(angle2);
y2 = r*sin(angle2);
add_rotated_angle_to_segments(x1, y1, x2, y2, 0.5*MU, segment);
add_rotated_angle_to_segments(x1, y1, x2, y2, 0.5*MU, segment, 0);
}
}
@ -1537,10 +1573,11 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
cycle = 1;
nsegments = NPOLY;
ngroups = 2;
for (i=0; i<nsegments; i++)
{
segment[i].group = 0;
segment[i].group = 1;
segment[i].inactivate = 0;
}
@ -1557,7 +1594,7 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
angle = -PI + (double)(i+1)*DPI/(double)NPOLY;
x2 = 0.7*LAMBDA*(1.0 + cos(angle));
y2 = 0.5*LAMBDA*sin(angle);
add_rotated_angle_to_segments(x1, y1, x2, y2, 0.02, segment);
add_rotated_angle_to_segments(x1, y1, x2, y2, 0.02, segment, 0);
}
/* compute intersection point of nozzle and ellipse */
@ -1575,18 +1612,18 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
y1 = nozzle_width(x1, nozy, NOZZLE_SHAPE);
x2 = x1 + dx;
y2 = nozzle_width(x2, nozy, NOZZLE_SHAPE);
add_rotated_angle_to_segments(x1, y1, x2, y2, 0.02, segment);
add_rotated_angle_to_segments(x1, y1, x2, y2, 0.02, segment, 0);
}
add_rotated_angle_to_segments(x2, y2, nozx, nozy, 0.02, segment);
add_rotated_angle_to_segments(x2, y2, nozx, nozy, 0.02, segment, 0);
for (i=0; i<nsides; i++)
{
x1 = -LAMBDA + dx*(double)(i-1);
y1 = -nozzle_width(x1, nozy, NOZZLE_SHAPE);
x2 = x1 + dx;
y2 = -nozzle_width(x2, nozy, NOZZLE_SHAPE);
add_rotated_angle_to_segments(x1, y1, x2, y2, 0.02, segment);
add_rotated_angle_to_segments(x1, y1, x2, y2, 0.02, segment, 0);
}
add_rotated_angle_to_segments(x2, y2, nozx, -nozy, 0.02, segment);
add_rotated_angle_to_segments(x2, y2, nozx, -nozy, 0.02, segment, 0);
/* closing segment */
segment[nsegments].x1 = nozx;
@ -1696,7 +1733,7 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
y1 = YMIN + r - (r + MU)*sin(angle);
x2 = XMAX - r + (r + MU)*cos(angle + dangle);
y2 = YMIN + r - (r + MU)*sin(angle + dangle);
add_rotated_angle_to_segments(x1, y1, x2, y2, MU, segment);
add_rotated_angle_to_segments(x1, y1, x2, y2, MU, segment, 0);
}
cycle = 0;
@ -1727,17 +1764,17 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
{
width = 0.1*LAMBDA;
segment[0].x1 = -LAMBDA;
segment[0].y1 = -width;
segment[0].x1 = SEGMENTS_X0 - LAMBDA;
segment[0].y1 = SEGMENTS_Y0 - width;
segment[1].x1 = -LAMBDA;
segment[1].y1 = width;
segment[1].x1 = SEGMENTS_X0 - LAMBDA;
segment[1].y1 = SEGMENTS_Y0 + width;
segment[2].x1 = LAMBDA;
segment[2].y1 = width;
segment[2].x1 = SEGMENTS_X0 + LAMBDA;
segment[2].y1 = SEGMENTS_Y0 + width;
segment[3].x1 = LAMBDA;
segment[3].y1 = -width;
segment[3].x1 = SEGMENTS_X0 + LAMBDA;
segment[3].y1 = SEGMENTS_Y0 - width;
cycle = 1;
nsegments = 4;
@ -1753,7 +1790,7 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
}
case (S_DAM_BRICKS):
{
add_rectangle_to_segments(DAM_WIDTH, BCYMIN - 0.5, -DAM_WIDTH, BCYMIN, segment, i);
add_rectangle_to_segments(DAM_WIDTH, BCYMIN - 0.5, -DAM_WIDTH, BCYMIN, segment, 0);
dy = 0.1*(LAMBDA - BCYMIN);
for (i=1; i<11; i++)
@ -1773,6 +1810,55 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
break;
}
case (S_HLINE_HOLE):
{
x = 0.15;
x1 = XMAX + 1.0;
y1 = 0.0;
width = 0.05;
add_rectangle_to_segments(x1, y1 - width, x, y1, segment, 0);
add_rectangle_to_segments(-x, y1 - width, -x1, y1, segment, 0);
/* closing segment */
segment[nsegments].x1 = -x;
segment[nsegments].y1 = y1;
segment[nsegments].x2 = x;
segment[nsegments].y2 = y1;
nsegments++;
cycle = 0;
concave = 0; /* add_rectangle_to_segments already deals with concave corners */
for (i=0; i<nsegments; i++)
{
segment[i].inactivate = 0;
}
segment[nsegments-1].inactivate = 1;
break;
}
case (S_EXT_CIRCLE_RECT):
{
width = 0.1*LAMBDA;
add_rectangle_to_segments(SEGMENTS_X0 + LAMBDA, SEGMENTS_Y0 - width, SEGMENTS_X0 - LAMBDA, SEGMENTS_Y0 + width, segment, 1);
add_circle_to_segments(-SEGMENTS_X0, SEGMENTS_Y0, 0.5*LAMBDA, NPOLY, APOLY, segment, 2);
cycle = 0;
concave = 0;
nsegments = NPOLY + 4;
ngroups = 3;
for (i=0; i<nsegments; i++)
{
segment[i].concave = 1;
segment[i].inactivate = 0;
}
break;
}
default:
{
printf("Function init_segment_config not defined for this pattern \n");
@ -1840,14 +1926,16 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
}
else
{
if (iminus == 0) iminus = nsegments - 1;
if (iplus == nsegments - 1) iplus = 0;
if (iminus < 0) iminus = nsegments - 1;
if (iplus > nsegments - 1) iplus = 0;
}
angle = argument(segment[iplus].x1 - segment[i].x1, segment[iplus].y1 - segment[i].y1) + PID;
angle2 = argument(segment[i].x1 - segment[iminus].x1, segment[i].y1 - segment[iminus].y1) + PID;
if (angle2 < angle) angle2 += DPI;
segment[i].angle1 = angle;
segment[i].angle2 = angle2;
printf("i = %i, iplus = %i, iminus = %i, angle1 = %.0f, angle2 = %.0f\n", i, iplus, iminus, angle*360.0/DPI, angle2*360.0/DPI);
}
/* make copy of initial values in case of rotation/translation */
@ -1868,6 +1956,7 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
// printf("Segment %i: (x1, y1) = (%.3lg,%.3lg), (x2, y2) = (%.3lg,%.3lg)\n (nx, ny) = (%.3lg,%.3lg), c = %.3lg, length = %.3lg\n", i, segment[i].x1, segment[i].y1, segment[i].x2, segment[i].y2, segment[i].nx, segment[i].ny, segment[i].c, segment[i].length);
// if (segment[i].concave) printf("Concave with angles %.3lg Pi, %.3lg Pi\n", segment[i].angle1/PI, segment[i].angle2/PI);
// }
// sleep(4);
}
int in_rocket(double x, double y, int rocket_shape)
@ -2044,6 +2133,13 @@ int in_segment_region(double x, double y)
else if (vabs(y) > 0.1*LAMBDA + padding) return(1);
else return(0);
}
case (S_EXT_CIRCLE_RECT):
{
padding = 0.1;
if ((vabs(x - SEGMENTS_X0) < LAMBDA + padding)&&(vabs(y - SEGMENTS_Y0) < 0.1*LAMBDA + padding)) return(0);
else if (module2(x + SEGMENTS_X0, y - SEGMENTS_Y0) < 0.5*LAMBDA + padding) return(0);
else return(1);
}
default: return(1);
}
}
@ -3192,6 +3288,33 @@ void set_segment_group_color(int group, double lum, double rgb[3])
}
void draw_altitude_lines()
{
int i, i1, i2;
double x, y;
glColor3f(0.5, 0.5, 0.5);
glLineWidth(1);
i1 = (int)(YMIN + ytrack) - 1.0;
i2 = (int)(YMAX + ytrack) + 1.0;
for (i = i1; i < i2; i++)
{
y = (double)i - ytrack;
draw_line(BCXMIN, y, XMAX - 1.8, y);
}
i1 = (int)(XMIN + xtrack) - 1.0;
i2 = (int)(XMAX - 1.8 + xtrack) + 1.0;
for (i = i1; i < i2; i++)
{
x = (double)i - xtrack;
if (x < XMAX - 1.8) draw_line(x, YMIN, x, YMAX);
}
}
void draw_one_triangle(t_particle particle, int same_table[9*HASHMAX], int p, int q, int nsame)
{
double x, y, dx, dy;
@ -3450,6 +3573,9 @@ void draw_particles(t_particle particle[NMAXCIRCLES], int plot, double beta)
}
/* draw "altitude lines" */
if (ALTITUDE_LINES) draw_altitude_lines();
/* draw the bonds first */
if ((DRAW_BONDS)||(plot == P_BONDS))
{
@ -4996,7 +5122,7 @@ int initialize_configuration(t_particle particle[NMAXCIRCLES], t_hashgrid hashgr
pangle[i] = particle[i].omega;
if ((PLOT == P_INITIAL_POS)||(PLOT_B == P_INITIAL_POS))
particle[i].color_hue = 360.0*(particle[i].yc - YMIN)/(YMAX - YMIN);
particle[i].color_hue = 360.0*(particle[i].yc - INITYMIN)/(INITYMAX - INITYMIN);
}
/* initialize dummy values in case particles are added */
for (i=ncircles; i < NMAXCIRCLES; i++)
@ -5469,6 +5595,17 @@ void draw_trajectory_plot(t_group_data *group_speeds, int i)
first = 0;
}
if (ALTITUDE_LINES)
{
glLineWidth(1);
glColor3f(0.5, 0.5, 0.5);
for (j=1; j<(int)scaley + 1; j++)
{
y1 = plot_coord((double)j/scaley, plotymin, plotymax);
if (y1 < plotymax) draw_line(plotxmin, y1, plotxmax + 0.05, y1);
}
}
glLineWidth(2);
/* plot trajectories */
@ -5514,7 +5651,7 @@ void draw_trajectory_plot(t_group_data *group_speeds, int i)
// sprintf(message, "%.1f", VMAX_PLOT_SPEEDS);
sprintf(message, "y");
write_text_fixedwidth(plotxmin - 0.1, plotymax - 0.01, message);
write_text_fixedwidth(plotxmin - 0.1, plotymax - 0.005, message);
sprintf(message, "x");
write_text_fixedwidth(plotxmax - 0.1, plotymin - 0.1, message);
@ -5532,6 +5669,9 @@ void init_segment_group(t_segment segment[NMAXSEGMENTS], int group, t_group_segm
xc += 0.5*(segment[i].x1 + segment[i].x2);
yc += 0.5*(segment[i].y1 + segment[i].y2);
// printf("Segment group data %i: z1 = (%.3lg, %.3lg)\n z2 = (%.3lg, %.3lg)\n zc = (%.3lg, %.3lg)\n\n", group,
// segment[i].x1, segment[i].y1, segment[i].x2, segment[i].y2, xc, yc);
nseg_group++;
}

121
sub_maze.c
View File

@ -12,6 +12,8 @@ typedef struct
short int north, east, south, west; /* closed walls */
short int active; /* takes value 1 if currently active in RW path */
short int tested; /* takes value 1 if tested */
short int connected; /* takes value 1 if connected to exit */
short int closed; /* takes value 1 if no untested neighbours */
} t_maze;
@ -122,6 +124,8 @@ void init_maze_graph(t_maze maze[NXMAZE*NYMAZE])
n = nmaze(i, j);
maze[n].active = 0;
maze[n].tested = 0;
maze[n].connected = 0;
maze[n].closed = 0;
maze[n].north = 1;
maze[n].east = 1;
maze[n].south = 1;
@ -129,14 +133,17 @@ void init_maze_graph(t_maze maze[NXMAZE*NYMAZE])
}
}
int find_maze_path(t_maze maze[NXMAZE*NYMAZE], int n0)
int find_maze_path(t_maze maze[NXMAZE*NYMAZE], int n0, int *path, int *pathlength)
/* find a random walk path in the maze */
/* returns 0 or 1 depending on whether path reaches a tested cell or a deadend */
{
int active_counter = 0, i, n = n0, npaths, inext, nextcell, trial, nnext;
int active_counter = 0, i, n = n0, npaths, inext, nextcell, trial, nnext, deadend = 1, length = 0;
int next_table[4];
/* contruct random walk */
npaths = maze[n].nneighb;
path[0] = n0;
// while ((npaths > 0)&&(!maze[n].tested))
while ((npaths > 0))
{
@ -148,7 +155,7 @@ int find_maze_path(t_maze maze[NXMAZE*NYMAZE], int n0)
for (i=0; i<npaths; i++)
{
nextcell = maze[n].neighb[i];
if (!maze[nextcell].active)
if ((!maze[nextcell].active)&&((maze[nextcell].connected)||(!maze[nextcell].tested)))
{
next_table[nnext] = i;
nnext++;
@ -157,12 +164,15 @@ int find_maze_path(t_maze maze[NXMAZE*NYMAZE], int n0)
if (nnext == 0)
{
deadend = 1;
printf("Ended path\n");
// sleep(5);
npaths = 0;
maze[n].closed = 1;
}
else
{
deadend = 0;
inext = next_table[rand()%nnext];
nextcell = maze[n].neighb[inext];
switch(maze[n].directions[inext]){
@ -200,8 +210,18 @@ int find_maze_path(t_maze maze[NXMAZE*NYMAZE], int n0)
if (maze[n].tested) npaths = 0;
else npaths = maze[n].nneighb;
active_counter++;
if (length < NXMAZE*NYMAZE)
{
length++;
path[length] = n;
}
deadend = 0;
}
}
printf("Reached tested cell (%i, %i)\n", n%NXMAZE, n/NXMAZE);
if (!maze[n].connected) deadend = 1;
/* update cell status */
for (n=0; n<NXMAZE*NYMAZE; n++) if (maze[n].active)
@ -211,21 +231,108 @@ int find_maze_path(t_maze maze[NXMAZE*NYMAZE], int n0)
}
printf("Ended path\n");
if (deadend) printf("Deadend\n");
*pathlength = length;
printf("Path length %i \n", length);
return(active_counter);
return(deadend);
// return(active_counter);
}
void init_maze(t_maze maze[NXMAZE*NYMAZE])
void init_maze_old(t_maze maze[NXMAZE*NYMAZE])
/* init a maze */
{
int i;
int i, pathlength, *path;
init_maze_graph(maze);
for (i=0; i<RAND_SHIFT; i++) rand();
for (i=0; i<NXMAZE*NYMAZE; i++) if (!maze[i].tested) find_maze_path(maze, i);
for (i=0; i<NXMAZE*NYMAZE; i++) if (!maze[i].tested) find_maze_path(maze, i, path, &pathlength);
}
void init_maze(t_maze maze[NXMAZE*NYMAZE])
/* init a maze with exit at (nx, ny) */
{
int i, j, n, deadend, pathlength, newpathlength;
int *path, *newpath;
path = (int *)malloc(2*NXMAZE*NYMAZE*sizeof(short int));
newpath = (int *)malloc(2*NXMAZE*NYMAZE*sizeof(short int));
init_maze_graph(maze);
for (i=0; i<RAND_SHIFT; i++) rand();
find_maze_path(maze, 0, path, &pathlength);
for (n=0; n<pathlength; n++) maze[path[n]].connected = 1;
for (i=0; i<NXMAZE*NYMAZE; i++) if ((!maze[i].tested)&&(!maze[i].connected))
{
deadend = find_maze_path(maze, i, path, &pathlength);
if (!deadend) for (n=0; n<pathlength; n++) maze[path[n]].connected = 1;
j = 0;
printf("deadend = %i, pathlength = %i\n", deadend, pathlength);
// while ((deadend)&&(j < pathlength))
while (deadend)
{
j++;
if (j > pathlength) j = 0;
printf("j = %i\n", j);
// while (deadend)
if (!maze[path[j]].connected) deadend = find_maze_path(maze, path[j], newpath, &newpathlength);
if (!deadend) for (n=0; n<newpathlength; n++) maze[newpath[n]].connected = 1;
}
// for (n=0; n<newpathlength; n++) maze[newpath[n]].connected = 1;
for (j=0; j<NXMAZE*NYMAZE; j++) if (maze[j].tested) maze[j].connected = 1;
}
free(path);
free(newpath);
}
void init_maze_exit(int nx, int ny, t_maze maze[NXMAZE*NYMAZE])
/* init a maze with exit at (nx, ny) */
{
int i, j, n, deadend, pathlength, newpathlength;
int *path, *newpath;
path = (int *)malloc(2*NXMAZE*NYMAZE*sizeof(short int));
newpath = (int *)malloc(2*NXMAZE*NYMAZE*sizeof(short int));
init_maze_graph(maze);
for (i=0; i<RAND_SHIFT; i++) rand();
find_maze_path(maze, nmaze(nx, ny), path, &pathlength);
for (n=0; n<pathlength; n++) maze[path[n]].connected = 1;
for (i=0; i<NXMAZE*NYMAZE; i++) if ((!maze[i].tested)&&(!maze[i].connected))
{
deadend = find_maze_path(maze, i, path, &pathlength);
if (!deadend) for (n=0; n<pathlength; n++) maze[path[n]].connected = 1;
j = 0;
printf("deadend = %i, pathlength = %i\n", deadend, pathlength);
// while ((deadend)&&(j < pathlength))
while (deadend)
{
j++;
if (j > pathlength) j = 0;
printf("j = %i\n", j);
// while (deadend)
if (!maze[path[j]].connected) deadend = find_maze_path(maze, path[j], newpath, &newpathlength);
if (!deadend) for (n=0; n<newpathlength; n++) maze[newpath[n]].connected = 1;
}
// for (n=0; n<newpathlength; n++) maze[newpath[n]].connected = 1;
for (j=0; j<NXMAZE*NYMAZE; j++) if (maze[j].tested) maze[j].connected = 1;
}
free(path);
free(newpath);
}

164
sub_part_billiard.c
View File

@ -1,8 +1,8 @@
#include "colormaps.c"
#define DUMMY_ABSORBING -1000.0 /* dummy value of config[0] for absorbing circles */
#define DUMMY_ABSORBING 100000.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 */
#define DUMMY_SIDE_ABS -100000 /* 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;
@ -4617,7 +4617,7 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */
/* position varies between 0 and nsides */
/* returns number of hit setment */
{
double len, rlarge = 1000.0;
double len, rlarge = 1.0e8;
int nside;
nside = (int)conf[0];
@ -4647,7 +4647,7 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */
int vpolyline_xy(double config[8], double alpha, double pos[2])
/* determine initial configuration for start at point pos = (x,y) */
{
double c0, s0, a, b, c, t, dx, delta, s, xi, yi, margin = 1.0e-12, tmin, rlarge = 1000.0;
double c0, s0, a, b, c, t, dx, delta, s, xi, yi, margin = 1.0e-12, tmin, rlarge = 1.0e8;
double tval[nsides + ncircles], xint[nsides + ncircles], yint[nsides + ncircles], sint[nsides + ncircles];
int i, nt = 0, nsegment[nsides + ncircles], ntmin;
@ -5378,6 +5378,7 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */
{
/* not easy to implement for non-convex polygons */
if (POLYLINE_PATTERN == P_MAZE) return ((vabs(x) < 1.1*XMAX)&&(vabs(y) < 1.1*YMAX));
else if (POLYLINE_PATTERN == P_MAZE_DIAG) return ((vabs(x) < 1.1*XMAX)&&(vabs(y) < 1.1*YMAX));
else return(1);
break;
}
@ -6179,7 +6180,7 @@ void init_polyline(t_segment polyline[NMAXPOLY], t_circle circles[NMAXCIRCLES])
{
maze = (t_maze *)malloc(NXMAZE*NYMAZE*sizeof(t_maze));
init_maze(maze);
init_maze_exit(0, NYMAZE/2, maze);
/* build walls of maze */
dx = (YMAX - YMIN - 2.0*padding)/(double)(NXMAZE);
@ -6233,7 +6234,7 @@ void init_polyline(t_segment polyline[NMAXPOLY], t_circle circles[NMAXCIRCLES])
y1 = YMIN + padding + dy*((double)NYMAZE/2);
x1 = YMAX - padding + MAZE_XSHIFT;
polyline[nsides].x1 = x1;
polyline[nsides].y1 = YMIN - 1.0;
polyline[nsides].y1 = YMIN - 1000.0;
polyline[nsides].x2 = x1;
polyline[nsides].y2 = y1 - dy;
polyline[nsides].angle = PID;
@ -6242,14 +6243,14 @@ void init_polyline(t_segment polyline[NMAXPOLY], t_circle circles[NMAXCIRCLES])
polyline[nsides].x1 = x1;
polyline[nsides].y1 = y1;
polyline[nsides].x2 = x1;
polyline[nsides].y2 = YMAX + 1.0;
polyline[nsides].y2 = YMAX + 1000.0;
polyline[nsides].angle = PID;
nsides++;
/* left side of maze */
x1 = YMIN + padding + MAZE_XSHIFT;
polyline[nsides].x1 = x1;
polyline[nsides].y1 = YMIN - 1.0;
polyline[nsides].y1 = YMIN - 1000.0;
polyline[nsides].x2 = x1;
polyline[nsides].y2 = YMIN + padding;
polyline[nsides].angle = PID;
@ -6258,7 +6259,143 @@ void init_polyline(t_segment polyline[NMAXPOLY], t_circle circles[NMAXCIRCLES])
polyline[nsides].x1 = x1;
polyline[nsides].y1 = YMAX - padding;
polyline[nsides].x2 = x1;
polyline[nsides].y2 = YMAX + 1.0;
polyline[nsides].y2 = YMAX + 1000.0;
polyline[nsides].angle = PID;
nsides++;
free(maze);
break;
}
case (P_MAZE_DIAG):
{
maze = (t_maze *)malloc(NXMAZE*NYMAZE*sizeof(t_maze));
init_maze_exit(0, NYMAZE/2, maze);
/* build walls of maze */
dx = (YMAX - YMIN - 2.0*padding)/(double)(NXMAZE);
dy = (YMAX - YMIN - 2.0*padding)/(double)(NYMAZE);
nsides = 0;
ncircles = 0;
for (i=0; i<NXMAZE; i++)
for (j=0; j<NYMAZE; j++)
{
n = nmaze(i, j);
x1 = YMIN + padding + (double)i*dx + MAZE_XSHIFT;
y1 = YMIN + padding + (double)j*dy;
if (((i>0)||(j!=NYMAZE/2))&&(maze[n].west))
{
polyline[nsides].x1 = x1;
polyline[nsides].y1 = y1;
polyline[nsides].x2 = x1;
polyline[nsides].y2 = y1 + dy;
polyline[nsides].angle = PID;
nsides++;
}
if (maze[n].south)
{
polyline[nsides].x1 = x1;
polyline[nsides].y1 = y1;
polyline[nsides].x2 = x1 + dx;
polyline[nsides].y2 = y1;
polyline[nsides].angle = 0.0;
nsides++;
}
}
/* 45 degrees parts */
for (i=0; i<NXMAZE; i++)
for (j=0; j<NYMAZE; j++)
{
n = nmaze(i, j);
x1 = YMIN + padding + (double)i*dx + MAZE_XSHIFT;
y1 = YMIN + padding + (double)j*dy;
if ((maze[n].south)&&(maze[n].west))
{
polyline[nsides].x1 = x1 + 0.25*dx;
polyline[nsides].y1 = y1;
polyline[nsides].x2 = x1;
polyline[nsides].y2 = y1 + 0.25*dy;
polyline[nsides].angle = 1.5*PID;
nsides++;
}
if ((maze[n].south)&&(maze[n].east))
{
polyline[nsides].x1 = x1 + 0.75*dx;
polyline[nsides].y1 = y1;
polyline[nsides].x2 = x1 + dx;
polyline[nsides].y2 = y1 + 0.25*dy;
polyline[nsides].angle = 0.5*PID;
nsides++;
}
if ((maze[n].north)&&(maze[n].east))
{
polyline[nsides].x1 = x1 + dx;
polyline[nsides].y1 = y1 + 0.75*dy;
polyline[nsides].x2 = x1 + 0.75*dx;
polyline[nsides].y2 = y1 + dy;
polyline[nsides].angle = 1.5*PID;
nsides++;
}
if ((maze[n].north)&&(maze[n].west))
{
polyline[nsides].x1 = x1;
polyline[nsides].y1 = y1 + 0.75*dy;
polyline[nsides].x2 = x1 + 0.25*dx;
polyline[nsides].y2 = y1 + dy;
polyline[nsides].angle = 0.5*PID;
nsides++;
}
}
/* top side of maze */
polyline[nsides].x1 = YMIN + padding + MAZE_XSHIFT;
polyline[nsides].y1 = YMAX - padding;
polyline[nsides].x2 = YMAX - padding + MAZE_XSHIFT;
polyline[nsides].y2 = YMAX - padding;
polyline[nsides].angle = 0.0;
nsides++;
/* right side of maze */
y1 = YMIN + padding + dy*((double)NYMAZE/2);
x1 = YMAX - padding + MAZE_XSHIFT;
polyline[nsides].x1 = x1;
polyline[nsides].y1 = YMIN - 1000.0;
polyline[nsides].x2 = x1;
polyline[nsides].y2 = y1 - dy;
polyline[nsides].angle = PID;
nsides++;
polyline[nsides].x1 = x1;
polyline[nsides].y1 = y1;
polyline[nsides].x2 = x1;
polyline[nsides].y2 = YMAX + 1000.0;
polyline[nsides].angle = PID;
nsides++;
/* left side of maze */
x1 = YMIN + padding + MAZE_XSHIFT;
polyline[nsides].x1 = x1;
polyline[nsides].y1 = YMIN - 1000.0;
polyline[nsides].x2 = x1;
polyline[nsides].y2 = YMIN + padding;
polyline[nsides].angle = PID;
nsides++;
polyline[nsides].x1 = x1;
polyline[nsides].y1 = YMAX - padding;
polyline[nsides].x2 = x1;
polyline[nsides].y2 = YMAX + 1000.0;
polyline[nsides].angle = PID;
nsides++;
@ -6378,7 +6515,8 @@ int test_initial_condition(double *configs[NPARTMAX], int active[NPARTMAX], int
for (j=0; j<8; j++) newconf[j] = configs[i][j];
time = 0;
while ((time < 100000)&&(newconf[4] < 1000.0))
// while ((time < 100000)&&(newconf[4] < 1000.0))
while ((time < 100000)&&(newconf[0] < DUMMY_ABSORBING))
{
for (j=0; j<8; j++) conf[j] = newconf[j];
vbilliard(newconf);
@ -6390,9 +6528,11 @@ int test_initial_condition(double *configs[NPARTMAX], int active[NPARTMAX], int
printf("tmax = %i\n", time);
cosphi = (conf[6] - conf[4])/conf[3];
x2 = conf[4] + 10.0*cosphi;
x2 = conf[4] + conf[2]*cosphi;
// x2 = conf[4] + 10.0*cosphi;
if (x2 > 0.0)
if ((conf[0] >= DUMMY_ABSORBING)&&(x2 > XMAX))
// if (conf[0] >= DUMMY_ABSORBING)
{
active[i] = 1;
if (time > tmaxmax) tmaxmax = time;

185
sub_perco_3d.c Normal file
View File

@ -0,0 +1,185 @@
/* routines for 3D representation of percolation clusters, taken from sub_wave_3d_rde.c */
void xyz_to_xy(double x, double y, double z, double xy_out[2])
{
int i;
double s, t, xinter[3];
static double n2, m2, d, sm2, sn2, v[3], h[2], plane_ratio = 0.5;
static int first = 1;
if (((first)&&(REPRESENTATION_3D == REP_PROJ_3D))||(reset_view))
{
m2 = observer[0]*observer[0] + observer[1]*observer[1];
n2 = m2 + observer[2]*observer[2];
d = plane_ratio*n2;
sm2 = sqrt(m2);
sn2 = sqrt(n2);
h[0] = observer[1]/sm2;
h[1] = -observer[0]/sm2;
v[0] = -observer[0]*observer[2]/(sn2*sm2);
v[1] = -observer[1]*observer[2]/(sn2*sm2);
v[2] = m2/(sn2*sm2);
first = 0;
reset_view = 0;
// printf("h = (%.3lg, %.3lg)\n", h[0], h[1]);
// printf("v = (%.3lg, %.3lg, %.3lg)\n", v[0], v[1], v[2]);
}
switch (REPRESENTATION_3D) {
case (REP_AXO_3D):
{
for (i=0; i<2; i++)
xy_out[i] = x*u_3d[i] + y*v_3d[i] + z*w_3d[i];
break;
}
case (REP_PROJ_3D):
{
if (z > ZMAX_FACTOR*n2) z = ZMAX_FACTOR*n2;
z *= Z_SCALING_FACTOR;
s = observer[0]*x + observer[1]*y + observer[2]*z;
t = (d - s)/(n2 - s);
xinter[0] = t*observer[0] + (1.0-t)*x;
xinter[1] = t*observer[1] + (1.0-t)*y;
xinter[2] = t*observer[2] + (1.0-t)*z;
xy_out[0] = XSHIFT_3D + XY_SCALING_FACTOR*(xinter[0]*h[0] + xinter[1]*h[1]);
xy_out[1] = YSHIFT_3D + XY_SCALING_FACTOR*(xinter[0]*v[0] + xinter[1]*v[1] + xinter[2]*v[2]);
break;
}
}
}
// void draw_vertex_ij(int i, int j)
// {
// double xy[2];
//
// ij_to_xy(i, j, xy);
// // if (xy[1] > 0.0) printf("(i,j) = (%i,%i), (x,y) = (%.2lg,%.2lg)\n", i, j, xy[0], xy[1]);
// glVertex2d(xy[0], xy[1]);
// }
void draw_vertex_xyz(double xy[2], double z)
{
double xy_screen[2];
xyz_to_xy(xy[0], xy[1], z, xy_screen);
glVertex2d(xy_screen[0], xy_screen[1]);
}
void draw_vertex_x_y_z(double x, double y, double z)
{
double xy_screen[2];
xyz_to_xy(x, y, z, xy_screen);
glVertex2d(xy_screen[0], xy_screen[1]);
}
void draw_line_3d(double x1, double y1, double z1, double x2, double y2, double z2)
{
glBegin(GL_LINE_LOOP);
draw_vertex_x_y_z(x1, y1, z1);
draw_vertex_x_y_z(x2, y2, z2);
glEnd();
}
double angle_lum(double cosangle)
{
double mid_lum = 0.6;
if (cosangle > 0.0) return(mid_lum + (1.0 - mid_lum)*cosangle);
else return (mid_lum + 0.8*(1.0 - mid_lum)*cosangle);
}
void draw_cube(double x, double y, double z, double size, double rgb[3])
/* draw a cube */
{
double lum_factor;
/* front/back face */
if (observer[0] > 0.0)
{
lum_factor = angle_lum(light[0]);
glColor3f(rgb[0]*lum_factor, rgb[1]*lum_factor, rgb[2]*lum_factor);
glBegin(GL_TRIANGLE_FAN);
draw_vertex_x_y_z(x + size, y, z);
draw_vertex_x_y_z(x + size, y + size, z);
draw_vertex_x_y_z(x + size, y + size, z + size);
draw_vertex_x_y_z(x + size, y, z + size);
glEnd();
}
else
{
// cosangle = -light[0];
lum_factor = angle_lum(-light[0]);
glColor3f(rgb[0]*lum_factor, rgb[1]*lum_factor, rgb[2]*lum_factor);
glBegin(GL_TRIANGLE_FAN);
draw_vertex_x_y_z(x, y, z);
draw_vertex_x_y_z(x, y + size, z);
draw_vertex_x_y_z(x, y + size, z + size);
draw_vertex_x_y_z(x, y, z + size);
glEnd();
}
/* right/left face */
if (observer[1] > 0.0)
{
// cosangle = light[1];
lum_factor = angle_lum(light[1]);
glColor3f(rgb[0]*lum_factor, rgb[1]*lum_factor, rgb[2]*lum_factor);
glBegin(GL_TRIANGLE_FAN);
draw_vertex_x_y_z(x + size, y + size, z);
draw_vertex_x_y_z(x, y + size, z);
draw_vertex_x_y_z(x, y + size, z + size);
draw_vertex_x_y_z(x + size, y + size, z + size);
glEnd();
}
else
{
// cosangle = -light[1];
lum_factor = angle_lum(-light[1]);
glColor3f(rgb[0]*lum_factor, rgb[1]*lum_factor, rgb[2]*lum_factor);
glBegin(GL_TRIANGLE_FAN);
draw_vertex_x_y_z(x + size, y, z);
draw_vertex_x_y_z(x, y, z);
draw_vertex_x_y_z(x, y, z + size);
draw_vertex_x_y_z(x + size, y, z + size);
glEnd();
}
/* top face */
// cosangle = light[2];
lum_factor = angle_lum(light[2]);
glColor3f(rgb[0]*lum_factor, rgb[1]*lum_factor, rgb[2]*lum_factor);
glBegin(GL_TRIANGLE_FAN);
draw_vertex_x_y_z(x + size, y, z + size);
draw_vertex_x_y_z(x + size, y + size, z + size);
draw_vertex_x_y_z(x, y + size, z + size);
draw_vertex_x_y_z(x, y, z + size);
glEnd();
}
void viewpoint_schedule(int i)
/* change position of observer */
{
int j;
double angle, ca, sa;
static double observer_initial[3];
static int first = 1;
if (first)
{
for (j=0; j<3; j++) observer_initial[j] = observer[j];
first = 0;
}
angle = (ROTATE_ANGLE*DPI/360.0)*(double)i/(double)NSTEPS;
ca = cos(angle);
sa = sin(angle);
observer[0] = ca*observer_initial[0] - sa*observer_initial[1];
observer[1] = sa*observer_initial[0] + ca*observer_initial[1];
printf("Angle %.3lg, Observer position (%.3lg, %.3lg, %.3lg)\n", angle, observer[0], observer[1], observer[2]);
}

206
sub_rde.c
View File

@ -318,7 +318,7 @@ void init_vortex_state(double x, double y, double scale, double *phi[NFIELDS], s
/* phi[0] is stream function, phi[1] is vorticity */
{
int i, j;
double xy[2], dist2, module, phase, scale2;
double xy[2], dist2, module, phase, scale2, sign;
// scale2 = scale*scale;
for (i=0; i<NX; i++)
@ -330,10 +330,13 @@ void init_vortex_state(double x, double y, double scale, double *phi[NFIELDS], s
if (xy_in[i*NY+j])
{
dist2 = (xy[0]-x)*(xy[0]-x) + (xy[1]-y)*(xy[1]-y);
module = exp(-dist2/scale);
module = exp(-dist2/vabs(scale));
phi[1][i*NY+j] = module;
phi[0][i*NY+j] = -module; /* approximate, stream function should solve Poisson equation */
if (scale > 0.0) sign = 1.0;
else sign = -1.0;
phi[1][i*NY+j] += sign*module;
phi[0][i*NY+j] -= sign*module; /* approximate, stream function should solve Poisson equation */
}
else
{
@ -348,7 +351,7 @@ void add_vortex_state(double x, double y, double scale, double *phi[NFIELDS], sh
/* phi[0] is stream function, phi[1] is vorticity */
{
int i, j;
double xy[2], dist2, module, phase, scale2;
double xy[2], dist2, module, phase, scale2, sign;
// scale2 = scale*scale;
for (i=0; i<NX; i++)
@ -360,10 +363,13 @@ void add_vortex_state(double x, double y, double scale, double *phi[NFIELDS], sh
if (xy_in[i*NY+j])
{
dist2 = (xy[0]-x)*(xy[0]-x) + (xy[1]-y)*(xy[1]-y);
module = exp(-dist2/scale);
module = exp(-dist2/vabs(scale));
phi[1][i*NY+j] += module;
phi[0][i*NY+j] -= module; /* approximate, stream function should solve Poisson equation */
if (scale > 0.0) sign = 1.0;
else sign = -1.0;
phi[1][i*NY+j] += sign*module;
phi[0][i*NY+j] -= sign*module; /* approximate, stream function should solve Poisson equation */
}
else
{
@ -374,7 +380,7 @@ void add_vortex_state(double x, double y, double scale, double *phi[NFIELDS], sh
}
void init_shear_flow(double amp, double delta, double rho, int nx, int ny, double *phi[NFIELDS], short int xy_in[NX*NY])
void init_shear_flow(double amp, double delta, double rho, int nx, int ny, double yshift, double *phi[NFIELDS], short int xy_in[NX*NY])
/* initialise field with a shear flow */
/* phi[0] is stream function, phi[1] is vorticity */
/* amp is global amplitude */
@ -395,7 +401,7 @@ void init_shear_flow(double amp, double delta, double rho, int nx, int ny, doubl
ij_to_xy(i, j, xy);
xy_in[i*NY+j] = xy_in_billiard(xy[0],xy[1]);
y1 = xy[1]*(double)ny/YMAX;
y1 = xy[1]*(double)ny/YMAX - yshift;
while (y1 > 1.0) y1 -= 2.0;
while (y1 < -1.0) y1 += 2.0;
@ -424,6 +430,39 @@ void init_shear_flow(double amp, double delta, double rho, int nx, int ny, doubl
}
}
void init_laminar_flow(double amp, double modulation, double period, double yshift, double *phi[NFIELDS], short int xy_in[NX*NY])
/* initialise field with a laminar flow in x direction */
/* phi[0] is stream function, phi[1] is vorticity */
/* amp is global amplitude */
{
int i, j;
double xy[2], y1, a, b, f, cplus, cminus;
a = period*PI/YMAX;
// a = PID/YMAX;
for (i=0; i<NX; i++)
for (j=0; j<NY; j++)
{
ij_to_xy(i, j, xy);
xy_in[i*NY+j] = xy_in_billiard(xy[0],xy[1]);
y1 = xy[1] + yshift;
if (xy_in[i*NY+j])
{
phi[0][i*NY+j] = amp*(y1 + modulation*sin(a*y1)/a);
phi[1][i*NY+j] = amp*modulation*a*sin(a*y1);
// phi[0][i*NY+j] = amp*sin(a*xy[1]);
// phi[1][i*NY+j] = a*a*amp*sin(a*xy[1]);
}
else
{
phi[0][i*NY+j] = 0.0;
phi[1][i*NY+j] = 0.0;
}
}
}
/*********************/
/* animation part */
@ -604,7 +643,7 @@ void compute_gradient_xy(double phi[NX*NY], double gradient[2*NX*NY])
}
}
void compute_gradient_euler(double phi[NX*NY], double gradient[2*NX*NY])
void compute_gradient_euler(double phi[NX*NY], double gradient[2*NX*NY], double yshift)
/* compute the gradient of the field */
{
int i, j, iplus, iminus, jplus, jminus, padding = 0;
@ -640,14 +679,16 @@ void compute_gradient_euler(double phi[NX*NY], double gradient[2*NX*NY])
jminus = NY-1;
gradient[i*NY+j] = 0.5*(phi[iplus*NY+j] - phi[iminus*NY+j]);
gradient[NX*NY+i*NY+j] = 0.5*(phi[i*NY+jplus] - phi[i*NY+jminus]);
gradient[NX*NY+i*NY+j] = 0.5*(phi[i*NY+jplus] - phi[i*NY+jminus] - yshift);
// if (i == 1) printf("psi+ = %.5lg, psi- = %.5lg, gradient = %.5lg\n", phi[i*NY+jplus], phi[i*NY+jminus], gradient[NX*NY+i*NY+j]);
j = NY-1;
jplus = 0;
jminus = NY-2;
gradient[i*NY+j] = 0.5*(phi[iplus*NY+j] - phi[iminus*NY+j]);
gradient[NX*NY+i*NY+j] = 0.5*(phi[i*NY+jplus] - phi[i*NY+jminus]);
gradient[NX*NY+i*NY+j] = 0.5*(phi[i*NY+jplus] - phi[i*NY+jminus] + yshift);
}
for (j=1; j<NY-1; j++)
@ -661,6 +702,9 @@ void compute_gradient_euler(double phi[NX*NY], double gradient[2*NX*NY])
gradient[i*NY+j] = 0.5*(phi[iplus*NY+j] - phi[iminus*NY+j]);
gradient[NX*NY+i*NY+j] = 0.5*(phi[i*NY+jplus] - phi[i*NY+jminus]);
// printf("j = %i, psi+ = %.5lg, psi- = %.5lg, gradient = %.5lg\n", j, phi[i*NY+jplus], phi[i*NY+jminus], gradient[NX*NY+i*NY+j]);
i = NX-1;
iplus = 0;
@ -675,21 +719,21 @@ void compute_gradient_euler(double phi[NX*NY], double gradient[2*NX*NY])
j = 0; jplus = 1; jminus = NY-1;
gradient[i*NY+j] = 0.5*(phi[iplus*NY+j] - phi[iminus*NY+j]);
gradient[NX*NY+i*NY+j] = 0.5*(phi[i*NY+jplus] - phi[i*NY+jminus]);
gradient[NX*NY+i*NY+j] = 0.5*(phi[i*NY+jplus] - phi[i*NY+jminus] + yshift);
j = NY-1; jplus = 0; jminus = NY-2;
gradient[i*NY+j] = 0.5*(phi[iplus*NY+j] - phi[iminus*NY+j]);
gradient[NX*NY+i*NY+j] = 0.5*(phi[i*NY+jplus] - phi[i*NY+jminus]);
gradient[NX*NY+i*NY+j] = 0.5*(phi[i*NY+jplus] - phi[i*NY+jminus] - yshift);
i = NX-1; iplus = 0; iminus = NX-2;
j = 0; jplus = 1; jminus = NY-1;
gradient[i*NY+j] = 0.5*(phi[iplus*NY+j] - phi[iminus*NY+j]);
gradient[NX*NY+i*NY+j] = 0.5*(phi[i*NY+jplus] - phi[i*NY+jminus]);
gradient[NX*NY+i*NY+j] = 0.5*(phi[i*NY+jplus] - phi[i*NY+jminus] + yshift);
j = NY-1; jplus = 0; jminus = NY-2;
gradient[i*NY+j] = 0.5*(phi[iplus*NY+j] - phi[iminus*NY+j]);
gradient[NX*NY+i*NY+j] = 0.5*(phi[i*NY+jplus] - phi[i*NY+jminus]);
gradient[NX*NY+i*NY+j] = 0.5*(phi[i*NY+jplus] - phi[i*NY+jminus] - yshift);
}
void compute_gradient_rde(double phi[NX*NY], t_rde rde[NX*NY])
@ -844,7 +888,7 @@ void compute_field_log(double *phi[NFIELDS], t_rde rde[NX*NY])
for (j=0; j<NY; j++)
{
value = vabs(phi[1][i*NY+j]);
if (value < 1.0e-5) value = 1.0e-5;
if (value < LOG_MIN) value = LOG_MIN;
rde[i*NY+j].log_vorticity = LOG_SHIFT + LOG_SCALE*log(value);
}
}
@ -1164,7 +1208,7 @@ void compute_field_color_rde(double value, int cplot, int palette, double rgb[3]
// if (value < 0.0) value = -value;
// if (value < 1.0e-10) value = 1.0e-10;
// color_scheme_palette(COLOR_SCHEME, palette, LOG_SCALE*value + LOG_SHIFT, 1.0, 0, rgb);
color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 0, rgb);
color_scheme_palette(COLOR_SCHEME, palette, LOG_SCALE*value, 1.0, 0, rgb);
break;
}
case (Z_EULER_VORTICITY_ASYM):
@ -1688,6 +1732,70 @@ void draw_wave_3d_ij_rde(int i, int j, int movie, double *phi[NFIELDS], short in
}
}
void draw_wave_3d_ij_rde_periodic(int shiftx, int shifty, int i, int j, int movie, double *phi[NFIELDS],
short int xy_in[NX*NY], t_rde rde[NX*NY], double potential[NX*NY], int zplot, int cplot, int palette, int fade, double fade_value)
{
int k, l, draw = 1;
double xy[2], xy_screen[2], rgb[3], pos[2], ca, rgb_e[3], rgb_w[3], rgb_n[3], rgb_s[3];
double z, z_sw, z_se, z_nw, z_ne, z_mid, zw, ze, zn, zs, min = 1000.0, max = 0.0;
double xy_sw[2], xy_se[2], xy_nw[2], xy_ne[2], xy_mid[2];
double energy;
if (NON_DIRICHLET_BC)
draw = (xy_in[i*NY+j])&&(xy_in[(i+1)*NY+j])&&(xy_in[i*NY+j+1])&&(xy_in[(i+1)*NY+j+1]);
else draw = (TWOSPEEDS)||(xy_in[i*NY+j]);
if (draw)
{
if (AMPLITUDE_HIGH_RES > 0)
{
z_mid = compute_interpolated_colors_rde(i, j, rde, potential, palette, cplot,
rgb_e, rgb_w, rgb_n, rgb_s, &z_sw, &z_se, &z_nw, &z_ne,
fade, fade_value, movie);
ij_to_xy(i, j, xy_sw);
ij_to_xy(i+1, j, xy_se);
ij_to_xy(i, j+1, xy_nw);
ij_to_xy(i+1, j+1, xy_ne);
for (k=0; k<2; k++) xy_mid[k] = 0.25*(xy_sw[k] + xy_se[k] + xy_nw[k] + xy_ne[k]);
if (AMPLITUDE_HIGH_RES == 1)
{
glBegin(GL_TRIANGLE_FAN);
glColor3f(rgb_w[0], rgb_w[1], rgb_w[2]);
draw_vertex_xyz_shift(xy_mid, z_mid, shiftx, shifty);
draw_vertex_xyz_shift(xy_nw, z_nw, shiftx, shifty);
draw_vertex_xyz_shift(xy_sw, z_sw, shiftx, shifty);
glColor3f(rgb_s[0], rgb_s[1], rgb_s[2]);
draw_vertex_xyz_shift(xy_se, z_se, shiftx, shifty);
glColor3f(rgb_e[0], rgb_e[1], rgb_e[2]);
draw_vertex_xyz_shift(xy_ne, z_ne, shiftx, shifty);
glColor3f(rgb_n[0], rgb_n[1], rgb_n[2]);
draw_vertex_xyz_shift(xy_nw, z_nw, shiftx, shifty);
glEnd ();
}
}
else
{
glColor3f(rde[i*NY+j].rgb[0], rde[i*NY+j].rgb[1], rde[i*NY+j].rgb[2]);
glBegin(GL_TRIANGLE_FAN);
ij_to_xy(i, j, xy);
draw_vertex_xyz_shift(xy, adjust_field(*rde[i*NY+j].p_zfield[movie], potential[i*NY+j]), shiftx, shifty);
ij_to_xy(i+1, j, xy);
draw_vertex_xyz_shift(xy, adjust_field(*rde[(i+1)*NY+j].p_zfield[movie], potential[(i+1)*NY+j]), shiftx, shifty);
ij_to_xy(i+1, j+1, xy);
draw_vertex_xyz_shift(xy, adjust_field(*rde[(i+1)*NY+j+1].p_zfield[movie], potential[(i+1)*NY+j+1]), shiftx, shifty);
ij_to_xy(i, j+1, xy);
draw_vertex_xyz_shift(xy, adjust_field(*rde[i*NY+j+1].p_zfield[movie], potential[i*NY+j+1]), shiftx, shifty);
glEnd ();
}
}
}
void draw_wave_3d_rde(int movie, double *phi[NFIELDS], short int xy_in[NX*NY], t_rde rde[NX*NY], double potential[NX*NY],
@ -1742,6 +1850,27 @@ void draw_wave_3d_rde(int movie, double *phi[NFIELDS], short int xy_in[NX*NY], t
}
void draw_periodicised_wave_3d(int movie, double *phi[NFIELDS], short int xy_in[NX*NY], t_rde rde[NX*NY], double potential[NX*NY],
int zplot, int cplot, int palette, int fade, double fade_value)
/* a variant where the wave is repeated periodically in x and y directions (beta) */
{
int i, j, shiftx, shifty;
double observer_angle;
blank();
if (DRAW_BILLIARD) draw_billiard_3d(fade, fade_value);
if (!ROTATE_VIEW)
{
for (shiftx = -1; shiftx < 2; shiftx++)
for (shifty = -2; shifty < 2; shifty++)
for (i=BORDER_PADDING; i<NX-1-BORDER_PADDING; i++)
for (j=BORDER_PADDING; j<NY-1-BORDER_PADDING; j++)
draw_wave_3d_ij_rde_periodic(shiftx, shifty, i, j, movie, phi, xy_in, rde, potential, zplot, cplot, palette, fade, fade_value);
}
}
void draw_wave_rde(int movie, double *phi[NFIELDS], short int xy_in[NX*NY], t_rde rde[NX*NY], double potential[NX*NY],
int zplot, int cplot, int palette, int fade, double fade_value, int refresh)
{
@ -1760,8 +1889,43 @@ void draw_wave_rde(int movie, double *phi[NFIELDS], short int xy_in[NX*NY], t_rd
}
compute_cfield_rde(xy_in, cplot, palette, rde, fade, fade_value, movie);
if (PLOT_3D) draw_wave_3d_rde(movie, phi, xy_in, rde, potential, zplot, cplot, palette, fade, fade_value);
if (PLOT_3D)
{
if (DRAW_PERIODICISED)
draw_periodicised_wave_3d(movie, phi, xy_in, rde, potential, zplot, cplot, palette, fade, fade_value);
else draw_wave_3d_rde(movie, phi, xy_in, rde, potential, zplot, cplot, palette, fade, fade_value);
}
else draw_wave_2d_rde(xy_in, rde);
}
void draw_tracers(double *phi[NFIELDS], double tracers[2*N_TRACERS*NSTEPS], int time, int fade, double fade_value)
/* draw trajectories of tracers */
{
int tracer, t, t1, length = 50;
double x1, y1, x2, y2, lum;
glColor3f(1.0, 1.0, 1.0);
glLineWidth(1);
t1 = time - length;
if (t1 < 1) t1 = 1;
for (t = t1 + 1; t < time; t++)
for (tracer = 0; tracer < N_TRACERS; tracer++)
{
x1 = tracers[2*(t-1)*N_TRACERS + 2*tracer];
y1 = tracers[2*(t-1)*N_TRACERS + 2*tracer + 1];
x2 = tracers[2*t*N_TRACERS + 2*tracer];
y2 = tracers[2*t*N_TRACERS + 2*tracer + 1];
lum = 1.0 - 0.75*(double)(time - t)/(double)length;
if (fade) lum *= fade_value;
glColor3f(lum, lum, lum);
if (module2(x2 - x1, y2 - y1) < 0.2) draw_line(x1, y1, x2, y2);
// printf("time = %i, tracer = %i, coord = %i, x1 = %.2lg, y1 = %.2lg, x2 = %.2lg, y2 = %.2lg\n", t, tracer,2*t*N_TRACERS + 2*tracer, x1, y1, x2, y2);
}
}

351
sub_wave.c
View File

@ -3,6 +3,7 @@
/*********************/
#include "colors_waves.c"
#define TIFF_FREE_PERIOD 10
int writetiff_new(char *filename, char *description, int x, int y, int width, int height, int compression)
{
@ -63,6 +64,7 @@ int writetiff(char *filename, char *description, int x, int y, int width, int he
TIFF *file;
GLubyte *image, *p;
int i;
static int counter = 0;
file = TIFFOpen(filename, "w");
if (file == NULL)
@ -106,7 +108,16 @@ int writetiff(char *filename, char *description, int x, int y, int width, int he
p += width * sizeof(GLubyte) * 3;
}
/* added 9/9/22 and removed again, since it produces an unwanted "band" on the right */
// free(image); /* prevents RAM consumption*/
/* readded 5/11/22 */
if (SAVE_MEMORY)
{
counter++;
if (counter%TIFF_FREE_PERIOD == 0)
{
free(image); /* prevents RAM consumption*/
counter = 0;
}
}
TIFFClose(file);
return 0;
}
@ -346,6 +357,24 @@ void xy_to_ij(double x, double y, int ij[2])
}
void xy_to_ij_safe(double x, double y, int ij[2])
/* convert (x,y) position to (i,j) in table representing wave, making sure (i,j) are between 0 and NX or NY */
{
double x1, y1;
x1 = (x - XMIN)/(XMAX - XMIN);
y1 = (y - YMIN)/(YMAX - YMIN);
ij[0] = (int)(x1 * (double)NX);
ij[1] = (int)(y1 * (double)NY);
if (ij[0] < 0) ij[0] = 0;
if (ij[0] > NX-1) ij[0] = NX-1;
if (ij[1] < 0) ij[1] = 0;
if (ij[1] > NY-1) ij[1] = NY-1;
}
void xy_to_pos(double x, double y, double pos[2])
/* convert (x,y) position to double-valued position in table representing wave */
{
@ -1606,13 +1635,14 @@ int compute_maze_coordinates(t_rectangle polyrect[NMAXPOLY], int closed)
{
t_maze* maze;
int i, j, n, nsides = 0, ropening;
double dx, dy, x1, y1, padding = 0.02, pos[2], width = 0.02;
double dx, dy, x1, y1, x0, padding = 0.02, pos[2], width = 0.02;
maze = (t_maze *)malloc(NXMAZE*NYMAZE*sizeof(t_maze));
init_maze(maze);
/* build walls of maze */
// x0 = LAMBDA - 1.0;
dx = (YMAX - YMIN - 2.0*padding)/(double)(NXMAZE);
dy = (YMAX - YMIN - 2.0*padding)/(double)(NYMAZE);
@ -1918,6 +1948,7 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_
double l2, r2, r2mu, omega, b, c, angle, z, x1, y1, x2, y2, u, v, u1, v1, dx, dy, width, alpha, s, a, r, height;
int i, j, k, k1, k2, condition = 0, m;
static int first = 1, nsides;
static double h, hh, ra, rb;
switch (b_domain) {
case (D_NOTHING):
@ -1937,6 +1968,28 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_
else return(0);
break;
}
case (D_EXT_ELLIPSE):
{
if (x*x/(LAMBDA*LAMBDA) + y*y/(MU*MU) > 1.0) return(1);
else return(0);
break;
}
case (D_EXT_ELLIPSE_CURVED):
{
y1 = y + 0.4*x*x;
if (x*x/(LAMBDA*LAMBDA) + y1*y1/(MU*MU) > 1.0) return(1);
else return(0);
break;
}
case (D_EXT_ELLIPSE_CURVED_BDRY):
{
if (y > YMAX - 0.05) return(0);
if (y < YMIN + 0.05) return(0);
y1 = y + 0.4*x*x;
if (x*x/(LAMBDA*LAMBDA) + y1*y1/(MU*MU) > 1.0) return(1);
else return(0);
break;
}
case (D_STADIUM):
{
if ((x > -0.5*LAMBDA)&&(x < 0.5*LAMBDA)&&(y > -1.0)&&(y < 1.0)) return(1);
@ -2485,6 +2538,49 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_
if ((x > polyrect[i].x1)&&(x < polyrect[i].x2)&&(y > polyrect[i].y1)&&(y < polyrect[i].y2)) return(0);
return(1);
}
case (D_CHESSBOARD):
{
i = (int)(vabs(x)/LAMBDA + 0.5);
j = (int)(vabs(y)/LAMBDA + 0.5);
if ((i+j)%2 == 0) return(1);
else return(0);
}
case (D_TRIANGLE_TILES):
{
if (first)
{
h = LAMBDA/(2.0*sqrt(3.0));
hh = h*3.0;
first = 0;
}
i = (int)((y + h)/hh + 10.0);
y1 = sin(DPI/3.0)*x - 0.5*y;
j = (int)((y1 + h)/hh + 10.0);
y1 = sin(-DPI/3.0)*x -0.5*y;
k = (int)((y1 + h)/hh + 10.0);
if ((i+j+k)%2 == 0) return(1);
else return(0);
}
case (D_HEX_TILES):
{
if (first)
{
ra = -1.0/sqrt(3.0);
rb = -2.0*ra;
first = 0;
}
x1 = (x + ra*y)/LAMBDA + 30.0;
y1 = rb*y/LAMBDA + 30.0;
x1 = x1 - (double)(3*(int)(x1/3.0));
y1 = y1 - (double)(3*(int)(y1/3.0));
if ((x1 > 2.0)&&(y1 < 1.0)) return(1);
if ((x1 < 1.0)&&(y1 > 2.0)) return(1);
if (x1 + y1 < 1.0) return(1);
if (x1 + y1 > 5.0) return(1);
return(0);
}
case (D_MENGER):
{
x1 = 0.5*(x+1.0);
@ -2675,11 +2771,30 @@ void tvertex_lineto(t_vertex z)
}
void hex_transfo(double u, double v, double *x, double *y)
/* linear transformation of plane used for hex tiles */
{
static double ra, rb;
static int first = 1;
if (first)
{
ra = 0.5;
rb = 0.5*sqrt(3.0);
first = 0;
}
*x = u + ra*v;
*y = rb*v;
}
void draw_billiard(int fade, double fade_value) /* draws the billiard boundary */
{
double x0, x, y, x1, y1, x2, y2, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z, l, width, a, b, c, ymax, height;
int i, j, k, k1, k2, mr2;
double x0, y0, x, y, x1, y1, x2, y2, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z, l, width, a, b, c, ymax, height;
int i, j, k, k1, k2, mr2, ntiles;
static int first = 1, nsides;
static double h, hh, sqr3;
if (fade)
{
@ -2738,6 +2853,67 @@ void draw_billiard(int fade, double fade_value) /* draws the billiard bound
}
break;
}
case (D_EXT_ELLIPSE):
{
glBegin(GL_LINE_LOOP);
for (i=0; i<=NSEG; i++)
{
phi = (double)i*DPI/(double)NSEG;
x = LAMBDA*cos(phi);
y = MU*sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
/* draw foci */
if (FOCI)
{
if (fade) glColor3f(0.3*fade_value, 0.3*fade_value, 0.3*fade_value);
else glColor3f(0.3, 0.3, 0.3);
x0 = sqrt(LAMBDA*LAMBDA-MU*MU);
glLineWidth(2);
glEnable(GL_LINE_SMOOTH);
draw_circle(x0, 0.0, r, NSEG);
draw_circle(-x0, 0.0, r, NSEG);
}
break;
}
case (D_EXT_ELLIPSE_CURVED):
{
glBegin(GL_LINE_LOOP);
for (i=0; i<=NSEG; i++)
{
phi = (double)i*DPI/(double)NSEG;
x = LAMBDA*cos(phi);
y = MU*sin(phi) - 0.4*x*x;
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
break;
}
case (D_EXT_ELLIPSE_CURVED_BDRY):
{
glBegin(GL_LINE_LOOP);
for (i=0; i<=NSEG; i++)
{
phi = (double)i*DPI/(double)NSEG;
x = LAMBDA*cos(phi);
y = MU*sin(phi) - 0.4*x*x;
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
draw_line(XMIN, YMAX - 0.05, XMAX, YMAX - 0.05);
draw_line(XMIN, YMIN + 0.05, XMAX, YMIN + 0.05);
break;
}
case (D_STADIUM):
{
glBegin(GL_LINE_LOOP);
@ -3631,6 +3807,109 @@ void draw_billiard(int fade, double fade_value) /* draws the billiard bound
draw_filled_rectangle(polyrect[i].x1, polyrect[i].y1, polyrect[i].x2, polyrect[i].y2);
break;
}
case (D_CHESSBOARD):
{
glLineWidth(BOUNDARY_WIDTH);
x = 0.5*LAMBDA;
while (x < XMAX)
{
draw_line(x, YMIN, x, YMAX);
x += LAMBDA;
}
x = -0.5*LAMBDA;
while (x > XMIN)
{
draw_line(x, YMIN, x, YMAX);
x -= LAMBDA;
}
y = 0.5*LAMBDA;
while (y < YMAX)
{
draw_line(XMIN, y, XMAX, y);
y += LAMBDA;
}
y = -0.5*LAMBDA;
while (y > YMIN)
{
draw_line(XMIN, y, XMAX, y);
y -= LAMBDA;
}
break;
}
case (D_TRIANGLE_TILES):
{
if (first)
{
h = LAMBDA/(2.0*sqrt(3.0));
hh = h*3.0;
sqr3 = sqrt(3.0);
first = 0;
}
glLineWidth(BOUNDARY_WIDTH);
y = -h;
while (y < YMAX)
{
draw_line(XMIN, y, XMAX, y);
y += hh;
}
y = -h;
while (y > YMIN)
{
draw_line(XMIN, y, XMAX, y);
y -= hh;
}
x = -0.5*LAMBDA;
y = -h;
while (x < 1.5*XMAX)
{
draw_line(x - 10.0, y - 10.0*sqr3, x + 10.0, y + 10.0*sqr3);
draw_line(x - 10.0, y + 10.0*sqr3, x + 10.0, y - 10.0*sqr3);
x += LAMBDA;
}
x = -0.5*LAMBDA;
while (x > 1.5*XMIN)
{
draw_line(x - 10.0, y - 10.0*sqr3, x + 10.0, y + 10.0*sqr3);
draw_line(x - 10.0, y + 10.0*sqr3, x + 10.0, y - 10.0*sqr3);
x -= LAMBDA;
}
break;
}
case (D_HEX_TILES):
{
ntiles = (int)(XMAX/LAMBDA) + 1;
for (i=-ntiles; i<ntiles; i++)
for (j=-ntiles; j<ntiles; j++)
{
x0 = 3.0*LAMBDA*(double)i;
y0 = 3.0*LAMBDA*(double)j;
hex_transfo(x0, y0 + LAMBDA, &x, &y);
hex_transfo(x0, y0 + 2.0*LAMBDA, &x1, &y1);
draw_line(x, y, x1, y1);
hex_transfo(x0 + LAMBDA, y0 + 2.0*LAMBDA, &x2, &y2);
draw_line(x1, y1, x2, y2);
hex_transfo(x0 + LAMBDA, y0 + 3.0*LAMBDA, &x1, &y1);
draw_line(x1, y1, x2, y2);
hex_transfo(x0 + 2.0*LAMBDA, y0 + LAMBDA, &x1, &y1);
draw_line(x1, y1, x2, y2);
hex_transfo(x0 + 3.0*LAMBDA, y0 + LAMBDA, &x2, &y2);
draw_line(x1, y1, x2, y2);
hex_transfo(x0 + 2.0*LAMBDA, y0, &x2, &y2);
draw_line(x1, y1, x2, y2);
hex_transfo(x0 + LAMBDA, y0, &x1, &y1);
draw_line(x1, y1, x2, y2);
draw_line(x1, y1, x, y);
hex_transfo(x0 + 2.0*LAMBDA, y0 + 3.0*LAMBDA, &x, &y);
hex_transfo(x0 + 3.0*LAMBDA, y0 + 2.0*LAMBDA, &x1, &y1);
draw_line(x, y, x1, y1);
}
break;
}
case (D_MENGER):
{
glLineWidth(3);
@ -3936,6 +4215,27 @@ void draw_color_scheme(double x1, double y1, double x2, double y2, int plot, dou
color_scheme(COLOR_SCHEME, value, 1.0, 1, rgb);
break;
}
case (P_ENERGY_FLUX):
{
value = dy_e*(double)(j - jmin)*100.0/E_SCALE;
if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, COLOR_PALETTE, value, 1.0, 1, rgb);
else color_scheme_palette(COLOR_SCHEME, COLOR_PALETTE, value, 1.0, 1, rgb);
// value = min + 1.0*dy*(double)(j - jmin);
// amp = 0.7*color_amplitude_linear(value, 1.0, 1);
// while (amp > 1.0) amp -= 2.0;
// while (amp < -1.0) amp += 2.0;
// amp_to_rgb(0.5*(1.0 + amp), rgb);
break;
}
case (P_TOTAL_ENERGY_FLUX):
{
value = min + 1.0*dy*(double)(j - jmin);
amp = 0.7*color_amplitude_linear(value, 1.0, 1);
while (amp > 1.0) amp -= 2.0;
while (amp < -1.0) amp += 2.0;
amp_to_rgb(0.5*(1.0 + amp), rgb);
break;
}
case (P_PHASE):
{
value = min + 1.0*dy*(double)(j - jmin);
@ -4042,6 +4342,27 @@ void draw_color_scheme_palette(double x1, double y1, double x2, double y2, int p
color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
break;
}
case (P_ENERGY_FLUX):
{
value = dy_e*(double)(j - jmin)*100.0/E_SCALE;
if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
// value = min + 1.0*dy*(double)(j - jmin);
// amp = 0.7*color_amplitude_linear(value, 1.0, 1);
// while (amp > 1.0) amp -= 2.0;
// while (amp < -1.0) amp += 2.0;
// amp_to_rgb(0.5*(1.0 + amp), rgb);
break;
}
case (P_TOTAL_ENERGY_FLUX):
{
value = min + 1.0*dy*(double)(j - jmin);
amp = 0.7*color_amplitude_linear(value, 1.0, 1);
while (amp > 1.0) amp -= 2.0;
while (amp < -1.0) amp += 2.0;
amp_to_rgb(0.5*(1.0 + amp), rgb);
break;
}
case (P_PHASE):
{
value = min + 1.0*dy*(double)(j - jmin);
@ -4137,6 +4458,7 @@ void draw_color_scheme_palette_fade(double x1, double y1, double x2, double y2,
case (P_LOG_ENERGY):
{
value = LOG_SHIFT + LOG_SCALE*log(dy_e*(double)(j - jmin)*100.0/E_SCALE);
// printf("value = %.2lg\n", value);
// if (value <= 0.0) value = 0.0;
color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
break;
@ -4148,6 +4470,27 @@ void draw_color_scheme_palette_fade(double x1, double y1, double x2, double y2,
color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
break;
}
case (P_ENERGY_FLUX):
{
value = dy_e*(double)(j - jmin);
if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
// value = min + 1.0*dy*(double)(j - jmin);
// amp = 0.7*color_amplitude_linear(value, 1.0, 1);
// while (amp > 1.0) amp -= 2.0;
// while (amp < -1.0) amp += 2.0;
// amp_to_rgb(0.5*(1.0 + amp), rgb);
break;
}
case (P_TOTAL_ENERGY_FLUX):
{
value = min + 1.0*dy*(double)(j - jmin);
amp = 0.7*color_amplitude_linear(value, 1.0, 1);
while (amp > 1.0) amp -= 2.0;
while (amp < -1.0) amp += 2.0;
amp_to_rgb(0.5*(1.0 + amp), rgb);
break;
}
case (P_PHASE):
{
value = min + 1.0*dy*(double)(j - jmin);

8
sub_wave_3d_rde.c
View File

@ -82,6 +82,14 @@ void draw_vertex_xyz(double xy[2], double z)
glVertex2d(xy_screen[0], xy_screen[1]);
}
void draw_vertex_xyz_shift(double xy[2], double z, int shiftx, int shifty)
{
double xy_screen[2];
xyz_to_xy(xy[0] + (double)shiftx*(XMAX - XMIN), xy[1] + (double)shifty*(YMAX - YMIN), z, xy_screen);
glVertex2d(xy_screen[0], xy_screen[1]);
}
void draw_vertex_x_y_z(double x, double y, double z)
{
double xy_screen[2];

101
wave_3d.c
View File

@ -45,44 +45,47 @@
#define MOVIE 0 /* set to 1 to generate movie */
#define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */
#define SAVE_MEMORY 0 /* set to 1 to save memory when writing tiff images */
/* General geometrical parameters */
#define WINWIDTH 1920 /* window width */
#define WINHEIGHT 1000 /* window height */
#define NX 2000 /* number of grid points on x axis */
#define NY 2000 /* number of grid points on y axis */
#define NX 2500 /* number of grid points on x axis */
#define NY 1250 /* number of grid points on y axis */
// #define NX 2700 /* number of grid points on x axis */
// #define NY 1350 /* number of grid points on y axis */
// // #define YMIN -1.041666667
// // #define YMAX 1.041666667 /* y interval for 9/16 aspect ratio */
// #define XMIN -2.0
// #define XMAX 2.0 /* x interval */
// #define YMIN -1.041666667
// #define YMAX 1.041666667 /* y interval for 9/16 aspect ratio */
#define XMIN -2.0
#define XMAX 2.0 /* x interval */
#define YMIN -1.041666667
#define YMAX 1.041666667 /* y interval for 9/16 aspect ratio */
#define HIGHRES 1 /* set to 1 if resolution of grid is double that of displayed image */
// #define WINWIDTH 1280 /* window width */
// #define WINHEIGHT 720 /* window height */
// // //
// #define NX 1500 /* number of grid points on x axis */
// #define NY 1500 /* 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 NX 2560 /* number of grid points on x axis */
// #define NY 1440 /* 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 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 JULIA_SCALE 0.8 /* scaling for Julia sets */
/* Choice of the billiard table */
#define B_DOMAIN 7 /* choice of domain shape, see list in global_pdes.c */
#define B_DOMAIN 57 /* choice of domain shape, see list in global_pdes.c */
#define CIRCLE_PATTERN 2 /* pattern of circles or polygons, see list in global_pdes.c */
@ -98,7 +101,7 @@
#define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */
#define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */
#define LAMBDA 0.7 /* parameter controlling the dimensions of domain */
#define LAMBDA 0.25 /* parameter controlling the dimensions of domain */
#define MU 0.0 /* parameter controlling the dimensions of domain */
#define NPOLY 3 /* number of sides of polygon */
#define APOLY 2.0 /* angle by which to turn polygon, in units of Pi/2 */
@ -137,10 +140,10 @@
#define AMPLITUDE 0.8 /* amplitude of periodic excitation */
#define ACHIRP 0.2 /* acceleration coefficient in chirp */
#define DAMPING 0.0 /* damping of periodic excitation */
#define COURANT 0.06 /* Courant number */
#define COURANTB 0.15 /* Courant number in medium B */
#define COURANT 0.05 /* Courant number */
#define COURANTB 0.1 /* Courant number in medium B */
#define GAMMA 0.0 /* damping factor in wave equation */
#define GAMMAB 0.015 /* damping factor in wave equation */
#define GAMMAB 0.0 /* 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 */
@ -151,8 +154,9 @@
/* Increasing COURANT speeds up the simulation, but decreases accuracy */
/* For similar wave forms, COURANT^2*GAMMA should be kept constant */
#define ADD_OSCILLATING_SOURCE 0 /* set to 1 to add an oscillating wave source */
#define OSCILLATING_SOURCE_PERIOD 30 /* period of oscillating source */
#define ADD_OSCILLATING_SOURCE 1 /* set to 1 to add an oscillating wave source */
#define OSCILLATING_SOURCE_PERIOD 100 /* period of oscillating source */
#define ALTERNATE_OSCILLATING_SOURCE 1 /* set to 1 to alternate sign of oscillating source */
/* Boundary conditions, see list in global_pdes.c */
@ -163,9 +167,9 @@
/* Parameters for length and speed of simulation */
// #define NSTEPS 200 /* number of frames of movie */
#define NSTEPS 1500 /* number of frames of movie */
#define NVID 10 /* number of iterations between images displayed on screen */
#define NSTEPS 1400 /* number of frames of movie */
// #define NSTEPS 250 /* number of frames of movie */
#define NVID 5 /* number of iterations between images displayed on screen */
#define NSEG 1000 /* number of segments of boundary */
#define INITIAL_TIME 0 /* time after which to start saving frames */
#define BOUNDARY_WIDTH 2 /* width of billiard boundary */
@ -182,12 +186,9 @@
/* Parameters of initial condition */
// #define INITIAL_AMP 0.75 /* amplitude of initial condition */
// #define INITIAL_VARIANCE 0.001 /* variance of initial condition */
// #define INITIAL_WAVELENGTH 0.05 /* wavelength of initial condition */
#define INITIAL_AMP 0.75 /* amplitude of initial condition */
#define INITIAL_VARIANCE 0.0002 /* variance of initial condition */
#define INITIAL_WAVELENGTH 0.01 /* wavelength of initial condition */
#define INITIAL_AMP 0.25 /* 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 */
@ -195,8 +196,8 @@
#define ZPLOT 103 /* wave height */
#define CPLOT 103 /* color scheme */
#define ZPLOT_B 105
#define CPLOT_B 105 /* plot type for second movie */
#define ZPLOT_B 107
#define CPLOT_B 107 /* plot type for second movie */
#define AMPLITUDE_HIGH_RES 1 /* set to 1 to increase resolution of plot */
#define SHADE_3D 1 /* set to 1 to change luminosity according to normal vector */
@ -218,13 +219,13 @@
#define REP_AXO_3D 0 /* linear projection (axonometry) */
#define REP_PROJ_3D 1 /* projection on plane orthogonal to observer line of sight */
#define ROTATE_VIEW 0 /* set to 1 to rotate position of observer */
#define ROTATE_VIEW 1 /* set to 1 to rotate position of observer */
#define ROTATE_ANGLE 360.0 /* total angle of rotation during simulation */
/* Color schemes */
#define COLOR_PALETTE 10 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 16 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 14 /* Color palette, see list in global_pdes.c */
#define BLACK 1 /* background */
@ -237,11 +238,12 @@
#define PHASE_FACTOR 20.0 /* factor in computation of phase in color scheme P_3D_PHASE */
#define PHASE_SHIFT 0.0 /* shift of phase in color scheme P_3D_PHASE */
#define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */
#define E_SCALE 100.0 /* scaling factor for energy representation */
#define LOG_SCALE 0.5 /* scaling factor for energy log representation */
#define E_SCALE 100.0 /* scaling factor for energy representation */
#define LOG_SCALE 0.75 /* scaling factor for energy log representation */
#define LOG_SHIFT 0.5 /* shift of colors on log scale */
#define LOG_ENERGY_FLOOR -10.0 /* floor value for log of (total) energy */
#define LOG_MEAN_ENERGY_SHIFT 1.0 /* additional shift for log of mean energy */
#define FLUX_SCALE 1.0e4 /* scaling factor for enegy flux represtnation */
#define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */
#define COLORHUE 260 /* initial hue of water color for scheme C_LUM */
@ -255,13 +257,11 @@
#define NYMAZE 7 /* height of maze */
#define MAZE_MAX_NGBH 4 /* max number of neighbours of maze cell */
#define RAND_SHIFT 24 /* seed of random number generator */
// 0, 9, 18, 20, 22
// #define RAND_SHIFT 58 /* seed of random number generator */
#define MAZE_XSHIFT 0.0 /* horizontal shift of maze */
#define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */
#define COLORBAR_RANGE 2.0 /* scale of color scheme bar */
#define COLORBAR_RANGE_B 5.0 /* scale of color scheme bar for 2nd part */
#define COLORBAR_RANGE 3.5 /* scale of color scheme bar */
#define COLORBAR_RANGE_B 3.0 /* scale of color scheme bar for 2nd part */
#define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */
#define SAVE_TIME_SERIES 0 /* set to 1 to save wave time series at a point */
@ -286,7 +286,7 @@ double light[3] = {0.816496581, -0.40824829, 0.40824829}; /* vector of "lig
double observer[3] = {8.0, 8.0, 8.0}; /* location of observer for REP_PROJ_3D representation */
int reset_view = 0; /* switch to reset 3D view parameters (for option ROTATE_VIEW) */
#define Z_SCALING_FACTOR 0.35 /* overall scaling factor of z axis for REP_PROJ_3D representation */
#define Z_SCALING_FACTOR 0.3 /* overall scaling factor of z axis for REP_PROJ_3D representation */
#define XY_SCALING_FACTOR 2.4 /* overall scaling factor for on-screen (x,y) coordinates after projection */
#define ZMAX_FACTOR 1.0 /* max value of z coordinate for REP_PROJ_3D representation */
#define XSHIFT_3D -0.1 /* overall x shift for REP_PROJ_3D representation */
@ -1028,7 +1028,7 @@ void animation()
// init_circular_wave_mod(polyline[85].x, polyline[85].y, phi, psi, xy_in);
init_circular_wave_mod(-0.85, 0.0, phi, psi, xy_in);
init_circular_wave_mod(0.0, 0.0, phi, psi, xy_in);
// init_wave_flat_mod(phi, psi, xy_in);
// add_circular_wave_mod(1.0, 1.0, 0.0, phi, psi, xy_in);
@ -1108,15 +1108,16 @@ void animation()
// add_circular_wave_mod(1.0, -1.0, 0.0, phi, psi, xy_in);
// add_circular_wave(1.0, -1.5*LAMBDA, 0.0, phi, psi, xy_in);
// add_circular_wave(-1.0, 0.6*cos((double)(period)*DPI/3.0), 0.6*sin((double)(period)*DPI/3.0), phi, psi, xy_in);
yshift = (double)period*a + (double)(period*period)*b;
add_circular_wave_mod(sign, -1.5 + yshift, 0.0, phi, psi, xy_in);
// yshift = (double)period*a + (double)(period*period)*b;
if (ALTERNATE_OSCILLATING_SOURCE) sign = -sign;
add_circular_wave_mod(sign, 0.0, 0.0, phi, psi, xy_in);
// speed = (a + 2.0*(double)(period)*b)/((double)(NVID));
// speed = 0.55*(a + 2.0*(double)(period)*b)/((double)(NVID*OSCILLATING_SOURCE_PERIOD));
speed = (a + 2.0*(double)(period)*b)/((double)(3*NVID*OSCILLATING_SOURCE_PERIOD));
printf("v = %.3lg, c = %.3lg\n", speed, c);
speed = speed/c;
sign = -sign;
period++;
// speed = (a + 2.0*(double)(period)*b)/((double)(3*NVID*OSCILLATING_SOURCE_PERIOD));
// printf("v = %.3lg, c = %.3lg\n", speed, c);
// speed = speed/c;
// sign = -sign;
// period++;
}
if (PRINT_SPEED) print_speed_3d(speed, 0, 1.0);

84
wave_billiard.c
View File

@ -45,6 +45,7 @@
#define MOVIE 0 /* set to 1 to generate movie */
#define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */
#define SAVE_MEMORY 0 /* set to 1 to save memory when writing tiff images */
/* General geometrical parameters */
@ -81,7 +82,7 @@
/* Choice of the billiard table */
#define B_DOMAIN 53 /* choice of domain shape, see list in global_pdes.c */
#define B_DOMAIN 57 /* choice of domain shape, see list in global_pdes.c */
#define CIRCLE_PATTERN 1 /* pattern of circles or polygons, see list in global_pdes.c */
@ -93,8 +94,8 @@
#define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */
#define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */
#define LAMBDA 1.0 /* parameter controlling the dimensions of domain */
#define MU 0.3 /* parameter controlling the dimensions of domain */
#define LAMBDA 0.25 /* parameter controlling the dimensions of domain */
#define MU 0.0 /* parameter controlling the dimensions of domain */
#define NPOLY 6 /* number of sides of polygon */
#define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */
#define MDEPTH 6 /* depth of computation of Menger gasket */
@ -130,10 +131,10 @@
#define AMPLITUDE 0.8 /* amplitude of periodic excitation */
#define ACHIRP 0.25 /* acceleration coefficient in chirp */
#define DAMPING 0.0 /* damping of periodic excitation */
#define COURANT 0.1 /* Courant number */
#define COURANTB 0.05 /* Courant number in medium B */
#define COURANT 0.05 /* Courant number */
#define COURANTB 0.1 /* Courant number in medium B */
#define GAMMA 0.0 /* damping factor in wave equation */
#define GAMMAB 5.0e-3 /* damping factor in wave equation */
#define GAMMAB 0.0 /* 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 */
@ -145,7 +146,8 @@
/* For similar wave forms, COURANT^2*GAMMA should be kept constant */
#define ADD_OSCILLATING_SOURCE 1 /* set to 1 to add an oscillating wave source */
#define OSCILLATING_SOURCE_PERIOD 50 /* period of oscillating source */
#define OSCILLATING_SOURCE_PERIOD 50 /* period of oscillating source */
#define ALTERNATE_OSCILLATING_SOURCE 1 /* set to 1 to alternate sign of oscillating source */
/* Boundary conditions, see list in global_pdes.c */
@ -153,8 +155,8 @@
/* Parameters for length and speed of simulation */
// #define NSTEPS 1000 /* number of frames of movie */
#define NSTEPS 2400 /* number of frames of movie */
#define NSTEPS 1700 /* number of frames of movie */
// #define NSTEPS 3500 /* number of frames of movie */
#define NVID 12 /* number of iterations between images displayed on screen */
#define NSEG 1000 /* number of segments of boundary */
#define INITIAL_TIME 0 /* time after which to start saving frames */
@ -172,8 +174,7 @@
/* Parameters of initial condition */
// #define INITIAL_AMP 0.75 /* amplitude of initial condition */
#define INITIAL_AMP 1.0 /* amplitude of initial condition */
#define INITIAL_AMP 0.35 /* amplitude of initial condition */
#define INITIAL_VARIANCE 0.0002 /* variance of initial condition */
#define INITIAL_WAVELENGTH 0.01 /* wavelength of initial condition */
@ -186,8 +187,8 @@
/* Color schemes */
#define COLOR_PALETTE 17 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 13 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE 13 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 18 /* Color palette, see list in global_pdes.c */
#define BLACK 1 /* background */
@ -198,9 +199,10 @@
#define PHASE_FACTOR 1.0 /* factor in computation of phase in color scheme P_3D_PHASE */
#define PHASE_SHIFT 0.0 /* shift of phase in color scheme P_3D_PHASE */
#define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */
#define E_SCALE 200.0 /* scaling factor for energy representation */
#define LOG_SCALE 1.0 /* scaling factor for energy log representation */
#define LOG_SHIFT 1.0 /* shift of colors on log scale */
#define E_SCALE 100.0 /* scaling factor for energy representation */
#define LOG_SCALE 0.25 /* scaling factor for energy log representation */
#define LOG_SHIFT 0.0 /* shift of colors on log scale */
#define FLUX_SCALE 1.0e4 /* scaling factor for enegy flux represtnation */
#define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */
#define COLORHUE 260 /* initial hue of water color for scheme C_LUM */
@ -211,8 +213,8 @@
#define HUEAMP -180.0 /* amplitude of variation of hue for color scheme C_HUE */
#define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */
#define COLORBAR_RANGE 1.5 /* scale of color scheme bar */
#define COLORBAR_RANGE_B 1.5 /* scale of color scheme bar for 2nd part */
#define COLORBAR_RANGE 2.0 /* scale of color scheme bar */
#define COLORBAR_RANGE_B 7.0 /* scale of color scheme bar for 2nd part */
#define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */
#define SAVE_TIME_SERIES 0 /* set to 1 to save wave time series at a point */
@ -234,6 +236,8 @@
#define FLOOR 0 /* set to 1 to limit wave amplitude to VMAX */
#define VMAX 10.0 /* max value of wave amplitude */
#define MEAN_FLUX (PLOT == P_TOTAL_ENERGY_FLUX)||(PLOT_B == P_TOTAL_ENERGY_FLUX)
#include "global_pdes.c" /* constants and global variables */
#include "sub_maze.c" /* support for generating mazes */
#include "sub_wave.c" /* common functions for wave_billiard, heat and schrodinger */
@ -568,7 +572,7 @@ void draw_color_bar_palette(int plot, double range, int palette, int fade, doubl
void animation()
{
double time, scale, ratio, startleft[2], startright[2], sign = 1.0, r2, xy[2], fade_value, yshift, speed = 0.0, a, b, c;
double *phi[NX], *psi[NX], *tmp[NX], *total_energy[NX], *color_scale[NX];
double *phi[NX], *psi[NX], *tmp[NX], *total_energy[NX], *color_scale[NX], *total_flux;
short int *xy_in[NX];
int i, j, s, sample_left[2], sample_right[2], period = 0, fade;
static int counter = 0;
@ -591,6 +595,8 @@ void animation()
color_scale[i] = (double *)malloc(NY*sizeof(double));
}
if (MEAN_FLUX) total_flux = (double *)malloc(4*NX*NY*sizeof(double));
/* initialise positions and radii of circles */
if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN == D_CIRCLES_IN_RECT)) ncircles = init_circle_config(circles);
else if (B_DOMAIN == D_POLYGONS) ncircles = init_polygon_config(polygons);
@ -627,6 +633,10 @@ void animation()
for (i=0; i<NX; i++)
for (j=0; j<NY; j++)
total_energy[i][j] = 0.0;
if (MEAN_FLUX)
for (i=0; i<4*NX*NY; i++)
total_flux[i] = 0.0;
ratio = (XMAX - XMIN)/8.4; /* for Tokarsky billiard */
@ -642,7 +652,7 @@ void animation()
// init_wave_flat(phi, psi, xy_in);
// init_circular_wave(sqrt(LAMBDA*LAMBDA - 1.0), 0.0, phi, psi, xy_in);
init_circular_wave(0.0, 0.3, phi, psi, xy_in);
init_circular_wave(0.0, 0.0, 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);
@ -680,8 +690,8 @@ void animation()
blank();
glColor3f(0.0, 0.0, 0.0);
// draw_wave(phi, psi, xy_in, 1.0, 0, PLOT);
if (HIGHRES) draw_wave_highres_palette(2, phi, psi, total_energy, xy_in, 1.0, 0, PLOT, COLOR_PALETTE, 0, 1.0);
else draw_wave_epalette(phi, psi, total_energy, color_scale, xy_in, 1.0, 0, PLOT, COLOR_PALETTE, 0, 1.0);
if (HIGHRES) draw_wave_highres_palette(2, phi, psi, total_energy, total_flux, xy_in, 1.0, 0, PLOT, COLOR_PALETTE, 0, 1.0);
else draw_wave_epalette(phi, psi, total_energy, total_flux, color_scale, xy_in, 1.0, 0, PLOT, COLOR_PALETTE, 0, 1.0);
draw_billiard(0, 1.0);
@ -717,8 +727,8 @@ void animation()
else scale = 1.0;
// draw_wave(phi, psi, xy_in, scale, i, PLOT);
if (HIGHRES) draw_wave_highres_palette(2, phi, psi, total_energy, xy_in, scale, i, PLOT, COLOR_PALETTE, 0, 1.0);
else draw_wave_epalette(phi, psi, total_energy, color_scale, xy_in, scale, i, PLOT, COLOR_PALETTE, 0, 1.0);
if (HIGHRES) draw_wave_highres_palette(2, phi, psi, total_energy, total_flux, xy_in, scale, i, PLOT, COLOR_PALETTE, 0, 1.0);
else draw_wave_epalette(phi, psi, total_energy, total_flux, color_scale, xy_in, scale, i, PLOT, COLOR_PALETTE, 0, 1.0);
for (j=0; j<NVID; j++)
{
// evolve_wave(phi, psi, phi_tmp, psi_tmp, xy_in);
@ -742,8 +752,8 @@ void animation()
/* add oscillating waves */
if ((ADD_OSCILLATING_SOURCE)&&(i%OSCILLATING_SOURCE_PERIOD == OSCILLATING_SOURCE_PERIOD - 1))
{
sign = -sign;
add_circular_wave(sign, 0.0, 0.3, phi, psi, xy_in);
if (ALTERNATE_OSCILLATING_SOURCE) sign = -sign;
add_circular_wave(sign, 0.0, 0.0, phi, psi, xy_in);
// add_circular_wave(1.0, -1.5*LAMBDA, 0.0, phi, psi, xy_in);
// add_circular_wave(-1.0, 0.6*cos((double)(period)*DPI/3.0), 0.6*sin((double)(period)*DPI/3.0), phi, psi, xy_in);
@ -771,8 +781,8 @@ void animation()
{
// draw_wave(phi, psi, xy_in, scale, i, PLOT_B);
if (HIGHRES)
draw_wave_highres_palette(2, phi, psi, total_energy, xy_in, scale, i, PLOT_B, COLOR_PALETTE_B, 0, 1.0);
else draw_wave_epalette(phi, psi, total_energy, color_scale, xy_in, scale, i, PLOT_B, COLOR_PALETTE_B, 0, 1.0);
draw_wave_highres_palette(2, phi, psi, total_energy, total_flux, xy_in, scale, i, PLOT_B, COLOR_PALETTE_B, 0, 1.0);
else draw_wave_epalette(phi, psi, total_energy, total_flux, color_scale, xy_in, scale, i, PLOT_B, COLOR_PALETTE_B, 0, 1.0);
draw_billiard(0, 1.0);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0);
if (PRINT_SPEED) print_speed(speed, 0, 1.0);
@ -799,8 +809,8 @@ void animation()
if (DOUBLE_MOVIE)
{
// draw_wave(phi, psi, xy_in, scale, i, PLOT);
if (HIGHRES) draw_wave_highres_palette(2, phi, psi, total_energy, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 0, 1.0);
else draw_wave_epalette(phi, psi, total_energy, color_scale, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 0, 1.0);
if (HIGHRES) draw_wave_highres_palette(2, phi, psi, total_energy, total_flux, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 0, 1.0);
else draw_wave_epalette(phi, psi, total_energy, total_flux, color_scale, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 0, 1.0);
draw_billiard(0, 1.0);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, 0, 1.0);
if (PRINT_SPEED) print_speed(speed, 0, 1.0);
@ -811,8 +821,8 @@ void animation()
{
fade_value = 1.0 - (double)i/(double)MID_FRAMES;
if (HIGHRES)
draw_wave_highres_palette(2, phi, psi, total_energy, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 1, fade_value);
else draw_wave_epalette(phi, psi, total_energy, color_scale, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 1, fade_value);
draw_wave_highres_palette(2, phi, psi, total_energy, total_flux, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 1, fade_value);
else draw_wave_epalette(phi, psi, total_energy, total_flux, color_scale, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 1, fade_value);
draw_billiard(1, fade_value);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, 1, fade_value);
if (PRINT_SPEED) print_speed(speed, 1, fade_value);
@ -823,8 +833,8 @@ void animation()
{
// draw_wave(phi, psi, xy_in, scale, i, PLOT_B);
if (HIGHRES)
draw_wave_highres_palette(2, phi, psi, total_energy, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B, 0, 1.0);
else draw_wave_epalette(phi, psi, total_energy, color_scale, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B, 0, 1.0);
draw_wave_highres_palette(2, phi, psi, total_energy, total_flux, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B, 0, 1.0);
else draw_wave_epalette(phi, psi, total_energy, total_flux, color_scale, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B, 0, 1.0);
draw_billiard(0, 1.0);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0);
if (PRINT_SPEED) print_speed(speed, 0, 1.0);
@ -835,8 +845,8 @@ void animation()
{
fade_value = 1.0 - (double)i/(double)END_FRAMES;
if (HIGHRES)
draw_wave_highres_palette(2, phi, psi, total_energy, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B, 1, fade_value);
else draw_wave_epalette(phi, psi, total_energy, color_scale, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B, 1, fade_value);
draw_wave_highres_palette(2, phi, psi, total_energy, total_flux, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B, 1, fade_value);
else draw_wave_epalette(phi, psi, total_energy, total_flux, color_scale, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B, 1, fade_value);
draw_billiard(1, fade_value);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 1, fade_value);
if (PRINT_SPEED) print_speed(speed, 1, fade_value);
@ -859,6 +869,8 @@ void animation()
free(color_scale[i]);
}
if (MEAN_FLUX) free(total_flux);
if (SAVE_TIME_SERIES)
{
fclose(time_series_left);

122
wave_common.c
View File

@ -273,6 +273,45 @@ double compute_energy(double *phi[NX], double *psi[NX], short int *xy_in[NX], in
}
void compute_energy_flux(double *phi[NX], double *psi[NX], short int *xy_in[NX], int i, int j, double *gx, double *gy, double *arg, double *module)
/* computes energy flux given by c^2 norm(nabla u) du/dt*/
{
double velocity, energy, gradientx, gradienty;
int iplus, iminus, jplus, jminus;
velocity = vabs(phi[i][j] - psi[i][j]);
iplus = (i+1); if (iplus == NX) iplus = NX-1;
iminus = (i-1); if (iminus == -1) iminus = 0;
jplus = (j+1); if (jplus == NY) jplus = NY-1;
jminus = (j-1); if (jminus == -1) jminus = 0;
gradientx = (phi[iplus][j] - phi[iminus][j]);
gradienty = (phi[i][jplus] - phi[i][jminus]);
*arg = argument(gradientx,gradienty);
if (*arg < 0.0) *arg += DPI;
if (*arg > DPI) *arg -= DPI;
if ((xy_in[i][j])||(TWOSPEEDS))
{
*module = velocity*module2(gradientx, gradienty);
*gx = velocity*gradientx;
*gy = velocity*gradienty;
}
else
{
*module = 0.0;
*gx = 0.0;
*gy = 0.0;
}
// if (xy_in[i][j]) return(E_SCALE*E_SCALE*(velocity*COURANT*module2(gradientx,gradienty)));
// else if (TWOSPEEDS) return(E_SCALE*E_SCALE*(velocity*COURANTB*module2(gradientx,gradienty)));
}
double compute_variance(double *phi[NX], double *psi[NX], short int *xy_in[NX])
/* compute the variance of the field, to adjust color scheme */
{
@ -682,12 +721,12 @@ void draw_wave_highres_diss(int size, double *phi[NX], double *psi[NX], double *
}
void draw_wave_epalette(double *phi[NX], double *psi[NX], double *total_energy[NX], double *color_scale[NX], short int *xy_in[NX],
double scale, int time, int plot, int palette, int fade, double fade_value)
void draw_wave_epalette(double *phi[NX], double *psi[NX], double *total_energy[NX], double *total_flux, double *color_scale[NX],
short int *xy_in[NX], double scale, int time, int plot, int palette, int fade, double fade_value)
/* same as draw_wave_e, but with color scheme specification */
{
int i, j, k, iplus, iminus, jplus, jminus;
double rgb[3], xy[2], x1, y1, x2, y2, velocity, field_value, energy, gradientx2, gradienty2, r2;
double rgb[3], xy[2], x1, y1, x2, y2, velocity, field_value, energy, gradientx2, gradienty2, r2, arg, mod, flux_factor, gx, gy, mgx, mgy, ffactor;
static double dtinverse = ((double)NX)/(COURANT*(XMAX-XMIN)), dx = (XMAX-XMIN)/((double)NX);
glBegin(GL_QUADS);
@ -755,6 +794,50 @@ void draw_wave_epalette(double *phi[NX], double *psi[NX], double *total_energy[N
color_scheme_palette(COLOR_SCHEME, palette, energy, scale, time, rgb);
break;
}
case (P_ENERGY_FLUX):
{
compute_energy_flux(phi, psi, xy_in, i, j, &gx, &gy, &arg, &mod);
// color_scheme_palette(C_ONEDIM_LINEAR, palette, arg/DPI, 1.0, 1, rgb);
// flux_factor = tanh(mod*E_SCALE);
// for (k=0; k<3; k++) rgb[k] *= flux_factor;
color_scheme_asym_palette(COLOR_SCHEME, palette, mod*FLUX_SCALE, scale, time, rgb);
break;
}
case (P_TOTAL_ENERGY_FLUX):
{
// ffactor = 1.0;
compute_energy_flux(phi, psi, xy_in, i, j, &gx, &gy, &arg, &mod);
total_flux[2*NX*NY + 2*j*NX + 2*i] += gx;
total_flux[2*NX*NY + 2*j*NX + 2*i + 1] += gy;
total_flux[2*j*NX + 2*i] += total_flux[2*NX*NY + 2*j*NX + 2*i];
total_flux[2*j*NX + 2*i + 1] += total_flux[2*NX*NY + 2*j*NX + 2*i + 1];
// total_flux[2*j*NX + 2*i] *= 1.0 + 1.0/(double)(time+1);
// total_flux[2*j*NX + 2*i + 1] *= 1.0 + 1.0/(double)(time+1);
// total_flux[2*j*NX + 2*i] *= ffactor;
// total_flux[2*j*NX + 2*i + 1] *= ffactor;
// total_flux[2*j*NX + 2*i] += gx;
// total_flux[2*j*NX + 2*i + 1] += gy;
// total_flux[2*j*NX + 2*i] *= 1.0/ffactor;
// total_flux[2*j*NX + 2*i + 1] *= 1.0/ffactor;
// total_flux[2*j*NX + 2*i] += 0.1*gx;
// total_flux[2*j*NX + 2*i + 1] += 0.1*gy;
mgx = total_flux[2*j*NX + 2*i];
mgy = total_flux[2*j*NX + 2*i + 1];
// mgx = total_flux[2*j*NX + 2*i]/sqrt((double)(time+1));
// mgy = total_flux[2*j*NX + 2*i + 1]/sqrt((double)(time+1));
// mgx = total_flux[2*j*NX + 2*i]/(1.0 + 0.1*log((double)(time+2)));
// mgy = total_flux[2*j*NX + 2*i + 1]/(1.0 + 0.1*log((double)(time+2)));
mod = module2(mgx, mgy);
arg = argument(mgx, mgy);
if (arg < 0.0) arg += DPI;
color_scheme_palette(C_ONEDIM_LINEAR, palette, arg/DPI, 1.0, 1, rgb);
flux_factor = tanh(mod*FLUX_SCALE);
for (k=0; k<3; k++) rgb[k] *= flux_factor;
// color_scheme_asym_palette(COLOR_SCHEME, palette, mod, scale, time, rgb);
break;
}
}
if (fade) for (k=0; k<3; k++) rgb[k] *= fade_value;
glColor3f(rgb[0], rgb[1], rgb[2]);
@ -770,11 +853,11 @@ void draw_wave_epalette(double *phi[NX], double *psi[NX], double *total_energy[N
}
void draw_wave_highres_palette(int size, double *phi[NX], double *psi[NX], double *total_energy[NX], short int *xy_in[NX], double scale, int time, int plot, int palette, int fade, double fade_value)
void draw_wave_highres_palette(int size, double *phi[NX], double *psi[NX], double *total_energy[NX], double *total_flux, short int *xy_in[NX], double scale, int time, int plot, int palette, int fade, double fade_value)
/* same as draw_wave_highres, but with color scheme option */
{
int i, j, k, iplus, iminus, jplus, jminus;
double rgb[3], xy[2], x1, y1, x2, y2, velocity, energy, gradientx2, gradienty2;
double rgb[3], xy[2], x1, y1, x2, y2, velocity, energy, gradientx2, gradienty2, arg, mod, flux_factor, gx, gy, mgx, mgy;
static double dtinverse = ((double)NX)/(COURANT*(XMAX-XMIN)), dx = (XMAX-XMIN)/((double)NX);
glBegin(GL_QUADS);
@ -834,6 +917,35 @@ void draw_wave_highres_palette(int size, double *phi[NX], double *psi[NX], doubl
color_scheme_palette(COLOR_SCHEME, palette, LOG_SHIFT + LOG_SCALE*log(total_energy[i][j]/(double)(time+1)), scale, time, rgb);
break;
}
case (P_ENERGY_FLUX):
{
compute_energy_flux(phi, psi, xy_in, i, j, &gx, &gy, &arg, &mod);
color_scheme_palette(C_ONEDIM_LINEAR, palette, arg/DPI, 1.0, 1, rgb);
flux_factor = tanh(mod*E_SCALE);
for (k=0; k<3; k++) rgb[k] *= flux_factor;
break;
}
case (P_TOTAL_ENERGY_FLUX):
{
compute_energy_flux(phi, psi, xy_in, i, j, &gx, &gy, &arg, &mod);
total_flux[2*j*NX + 2*i] *= 0.99;
total_flux[2*j*NX + 2*i + 1] *= 0.99;
total_flux[2*j*NX + 2*i] += gx;
total_flux[2*j*NX + 2*i + 1] += gy;
// mgx = total_flux[2*j*NX + 2*i]/(double)(time+1);
// mgy = total_flux[2*j*NX + 2*i + 1]/(double)(time+1);
mgx = total_flux[2*j*NX + 2*i];
mgy = total_flux[2*j*NX + 2*i + 1];
// mgx = total_flux[2*j*NX + 2*i]/log((double)(time+2));
// mgy = total_flux[2*j*NX + 2*i + 1]/log((double)(time+2));
mod = module2(mgx, mgy);
arg = argument(mgx, mgy);
if (arg < 0.0) arg += DPI;
color_scheme_palette(C_ONEDIM_LINEAR, palette, arg/DPI, 1.0, 1, rgb);
flux_factor = tanh(mod*E_SCALE);
for (k=0; k<3; k++) rgb[k] *= flux_factor;
break;
}
}
if (fade) for (k=0; k<3; k++) rgb[k] *= fade_value;

2
wave_comparison.c
View File

@ -44,6 +44,7 @@
#define MOVIE 0 /* set to 1 to generate movie */
#define DOUBLE_MOVIE 0 /* set to 1 to produce movies for wave height and energy simultaneously */
#define SAVE_MEMORY 0 /* set to 1 to save memory when writing tiff images */
#define TIME_LAPSE 0 /* set to 1 to add a time-lapse movie at the end */
#define TIME_LAPSE_FACTOR 4 /* factor of time-lapse movie */
@ -200,6 +201,7 @@
#define E_SCALE 200.0 /* scaling factor for energy representation */
#define LOG_SCALE 1.5 /* scaling factor for energy log representation */
#define LOG_SHIFT 1.0 /* shift of colors on log scale */
#define FLUX_SCALE 1.0e4 /* scaling factor for enegy flux represtnation */
#define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */
#define COLORHUE 260 /* initial hue of water color for scheme C_LUM */

2
wave_energy.c
View File

@ -43,6 +43,7 @@
#include <omp.h>
#define MOVIE 0 /* set to 1 to generate movie */
#define SAVE_MEMORY 0 /* set to 1 to save memory when writing tiff images */
#define WINWIDTH 1280 /* window width */
#define WINHEIGHT 720 /* window height */
@ -178,6 +179,7 @@
#define E_SCALE 500.0 /* scaling factor for energy representation */
#define LOG_SCALE 1.5 /* scaling factor for energy log representation */
#define LOG_SHIFT 1.0 /* shift of colors on log scale */
#define FLUX_SCALE 1.0e4 /* scaling factor for enegy flux represtnation */
#define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */
#define COLORHUE 260 /* initial hue of water color for scheme C_LUM */