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

Add files via upload

This commit is contained in:
nilsberglund-orleans 2022-05-14 20:58:08 +02:00 committed by GitHub
parent 1112516d3d
commit 701f09b2ca
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
13 changed files with 5074 additions and 606 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

103
global_3d.c
View File

@ -1,39 +1,106 @@
/* global variables and definition used by sub_wave_3d.c */
/* global variables and definitions used by sub_wave_3d.c */
/* plot types used by wave_3d */
#define P_3D_AMPLITUDE 101 /* color depends on amplitude */
#define P_3D_ANGLE 102 /* color depends on angle with fixed direction */
#define P_3D_AMP_ANGLE 103 /* color depends on amplitude, luminosity depends on angle */
#define P_3D_ENERGY 104 /* color depends on energy, luminosity depends on angle */
#define P_3D_LOG_ENERGY 105 /* color depends on logarithm of energy, luminosity depends on angle */
#define P_3D_AMPLITUDE 101 /* height/color depends on amplitude - DEPRECATED, instead use set SHADE_3D to 0 */
#define P_3D_ANGLE 102 /* height/color depends on angle with fixed direction - TO DO */
#define P_3D_AMP_ANGLE 103 /* height/color depends on amplitude, luminosity depends on angle */
#define P_3D_ENERGY 104 /* height/color depends on energy, luminosity depends on angle */
#define P_3D_LOG_ENERGY 105 /* height/color depends on logarithm of energy, luminosity depends on angle */
#define P_3D_TOTAL_ENERGY 106 /* height/color depends on total energy over time, luminosity depends on angle */
#define P_3D_LOG_TOTAL_ENERGY 107 /* height/color depends on log on total energy over time, luminosity depends on angle */
#define P_3D_MEAN_ENERGY 108 /* height/color depends on energy averaged over time, luminosity depends on angle */
#define P_3D_LOG_MEAN_ENERGY 109 /* height/color depends on log on energy averaged over time, luminosity depends on angle */
#define P_3D_PHASE 111 /* phase of wave */
/* Choice of simulated reaction-diffusion equation in rde.c */
#define E_HEAT 0 /* heat equation */
#define E_ALLEN_CAHN 1 /* Allen-Cahn equation */
#define E_CAHN_HILLIARD 2 /* Cahn-Hilliard equation */
#define E_FHN 3 /* FitzHugh-Nagumo equation */
#define E_RPS 4 /* rock-paper-scissors equation */
#define E_RPSLZ 41 /* rock-paper-scissors-lizard-Spock equation */
#define E_SCHRODINGER 5 /* Schrodinger equation */
/* Choice of potential */
#define POT_HARMONIC 1 /* harmonic oscillator */
/* plot types used by rde */
#define Z_AMPLITUDE 0 /* amplitude of first field */
#define Z_NORM_GRADIENT 22 /* gradient of polar angle */
#define Z_NORM_GRADIENTX 23 /* direction of gradient of u */
#define Z_RGB 20 /* RGB plot */
#define Z_POLAR 21 /* polar angle associated with RBG plot */
#define Z_NORM_GRADIENT 22 /* gradient of polar angle */
#define Z_ANGLE_GRADIENT 221 /* direction of polar angle */
#define Z_NORM_GRADIENTX 23 /* norm of gradient of u */
#define Z_ANGLE_GRADIENTX 231 /* direction of gradient of u */
#define Z_NORM_GRADIENT_INTENSITY 24 /* gradient and intensity of polar angle */
#define Z_VORTICITY 25 /* curl of polar angle */
/* for Schrodinger equation */
#define Z_MODULE 30 /* module squared of first two fields */
#define Z_ARGUMENT 31 /* argument of first two fields, with luminosity depending on module */
#define COMPUTE_THETA ((plot == P_POLAR)||(plot == P_GRADIENT)||(plot == P_GRADIENTX)||(plot == P_GRADIENT_INTENSITY)||(plot == P_VORTICITY))
#define COMPUTE_THETAZ ((zplot == P_POLAR)||(zplot == P_GRADIENT)||(zplot == P_GRADIENTX)||(zplot == P_GRADIENT_INTENSITY)||(zplot == P_VORTICITY))
/* for RPSLZ equation */
#define Z_MAXTYPE_RPSLZ 40 /* color of type with maximal density */
#define Z_THETA_RPSLZ 41 /* polar angle */
#define Z_NORM_GRADIENT_RPSLZ 42 /* gradient of polar angle */
/* macros to avoid unnecessary computations in 3D plots */
#define COMPUTE_THETA ((cplot == Z_POLAR)||(cplot == Z_NORM_GRADIENT)||(cplot == Z_ANGLE_GRADIENT)||(cplot == Z_NORM_GRADIENT_INTENSITY)||(cplot == Z_VORTICITY))
#define COMPUTE_THETAZ ((zplot == Z_POLAR)||(zplot == Z_NORM_GRADIENT)||(zplot == Z_ANGLE_GRADIENT)||(zplot == Z_NORM_GRADIENT_INTENSITY)||(zplot == Z_VORTICITY))
#define COMPUTE_ENERGY ((zplot == P_3D_ENERGY)||(cplot == P_3D_ENERGY)||(zplot == P_3D_LOG_ENERGY)||(cplot == P_3D_LOG_ENERGY)||(zplot == P_3D_TOTAL_ENERGY)||(cplot == P_3D_TOTAL_ENERGY)||(zplot == P_3D_LOG_TOTAL_ENERGY)||(cplot == P_3D_LOG_TOTAL_ENERGY)||(zplot == P_3D_MEAN_ENERGY)||(cplot == P_3D_MEAN_ENERGY)||(zplot == P_3D_LOG_MEAN_ENERGY)||(cplot == P_3D_LOG_MEAN_ENERGY))
#define COMPUTE_LOG_TOTAL_ENERGY ((ZPLOT == P_3D_LOG_TOTAL_ENERGY)||(CPLOT == P_3D_LOG_TOTAL_ENERGY)||(ZPLOT_B == P_3D_LOG_TOTAL_ENERGY)||(CPLOT_B == P_3D_LOG_TOTAL_ENERGY))
#define COMPUTE_LOG_MEAN_ENERGY ((ZPLOT == P_3D_LOG_MEAN_ENERGY)||(CPLOT == P_3D_LOG_MEAN_ENERGY)||(ZPLOT_B == P_3D_LOG_MEAN_ENERGY)||(CPLOT_B == P_3D_LOG_MEAN_ENERGY))
#define COMPUTE_LOG_ENERGY ((ZPLOT == P_3D_LOG_TOTAL_ENERGY)||(CPLOT == P_3D_LOG_TOTAL_ENERGY)||(ZPLOT_B == P_3D_LOG_TOTAL_ENERGY)||(CPLOT_B == P_3D_LOG_TOTAL_ENERGY)||(ZPLOT == P_3D_LOG_MEAN_ENERGY)||(CPLOT == P_3D_LOG_MEAN_ENERGY)||(ZPLOT_B == P_3D_LOG_MEAN_ENERGY)||(CPLOT_B == P_3D_LOG_MEAN_ENERGY))
#define COMPUTE_MEAN_ENERGY ((ZPLOT == P_3D_MEAN_ENERGY)||(CPLOT == P_3D_MEAN_ENERGY)||(ZPLOT_B == P_3D_MEAN_ENERGY)||(CPLOT_B == P_3D_MEAN_ENERGY)||(ZPLOT == P_3D_LOG_MEAN_ENERGY)||(CPLOT == P_3D_LOG_MEAN_ENERGY)||(ZPLOT_B == P_3D_LOG_MEAN_ENERGY)||(CPLOT_B == P_3D_LOG_MEAN_ENERGY))
#define COMPUTE_TOTAL_ENERGY ((ZPLOT == P_3D_TOTAL_ENERGY)||(CPLOT == P_3D_TOTAL_ENERGY)||(ZPLOT == P_3D_LOG_TOTAL_ENERGY)||(CPLOT == P_3D_LOG_TOTAL_ENERGY)||(ZPLOT == P_3D_MEAN_ENERGY)||(CPLOT == P_3D_MEAN_ENERGY)||(ZPLOT == P_3D_LOG_MEAN_ENERGY)||(CPLOT == P_3D_LOG_MEAN_ENERGY)||(ZPLOT_B == P_3D_TOTAL_ENERGY)||(CPLOT_B == P_3D_TOTAL_ENERGY)||(ZPLOT_B == P_3D_LOG_TOTAL_ENERGY)||(CPLOT_B == P_3D_LOG_TOTAL_ENERGY)||(ZPLOT_B == P_3D_MEAN_ENERGY)||(CPLOT_B == P_3D_MEAN_ENERGY)||(ZPLOT_B == P_3D_LOG_MEAN_ENERGY)||(CPLOT_B == P_3D_LOG_MEAN_ENERGY))
int global_time = 0;
/* structure used for color and height representations */
/* possible extra fields: zfield, cfield, interpolated coordinates */
typedef struct
{
double energy; /* wave energy */
double phase; /* wave phase */
double log_energy; /* log of wave energy */
double total_energy; /* total energy since beginning of simulation */
double cos_angle; /* cos of angle between normal vector and direction of light */
double rgb[3]; /* RGB color code */
double *p_zfield; /* pointer to z field */
double *p_cfield; /* pointer to color field */
double energy; /* wave energy */
double phase; /* wave phase */
double log_energy; /* log of wave energy */
double total_energy; /* total energy since beginning of simulation */
double log_total_energy; /* log of total energy since beginning of simulation */
double mean_energy; /* energy averaged since beginning of simulation */
double log_mean_energy; /* log of energy averaged since beginning of simulation */
double cos_angle; /* cos of angle between normal vector and direction of light */
double rgb[3]; /* RGB color code */
double *p_zfield[2]; /* pointers to z field (second pointer for option DOUBLE_MOVIE) */
double *p_cfield[2]; /* pointers to color field (second pointer for option DOUBLE_MOVIE) */
} t_wave;
typedef struct
{
double theta; /* angle for Rock-Paper-Scissors equation */
double nablax; /* gradient of first field */
double nablay; /* gradient of first field */
double field_norm; /* norm of field or gradient */
double field_arg; /* argument of field or gradient */
double curl; /* curl of field */
double cos_angle; /* cos of angle between normal vector and direction of light */
double rgb[3]; /* RGB color code */
double *p_zfield[2]; /* pointers to z field (second pointer for option DOUBLE_MOVIE) */
double *p_cfield[2]; /* pointers to color field (second pointer for option DOUBLE_MOVIE) */
} t_rde;

36
global_ljones.c
View File

@ -14,6 +14,7 @@
#define NMAXCIRCLES 20000 /* total number of circles/polygons (must be at least NCX*NCY for square grid) */
#define MAXNEIGH 20 /* max number of neighbours kept in memory */
#define NMAXOBSTACLES 100 /* max number of obstacles */
#define NMAXSEGMENTS 100 /* max number of repelling segments */
#define C_SQUARE 0 /* square grid of circles */
#define C_HEX 1 /* hexagonal/triangular grid of circles */
@ -38,6 +39,14 @@
#define O_GALTON_BOARD 1 /* Galton board pattern */
#define O_GENUS_TWO 2 /* obstacles in corners of L-shape domeain (for genus 2 b.c.) */
/* pattern of additional repelling segments */
#define S_RECTANGLE 0 /* segments forming a rectangle */
#define S_CUP 1 /* segments forming a cup (for increasing gravity) */
#define S_HOURGLASS 2 /* segments forming an hour glass */
#define S_PENTA 3 /* segments forming a pentagon with 3 angles of 120° and 2 right angles */
#define S_CENTRIFUGE 4 /* segments forming "centrifuge" (polygon with radial segments) */
#define S_POLY_ELLIPSE 5 /* segments forming a polygonal approximation of an ellipse */
/* particle interaction */
#define I_COULOMB 0 /* Coulomb force */
@ -66,6 +75,11 @@
#define BC_BOY 13 /* Boy surface/projective plane (periodic with twisted horizontal and vertical parts) */
#define BC_GENUS_TWO 14 /* surface of genus 2, obtained by identifying opposite sides of an L shape */
/* Regions for partial thermostat couplings */
#define TH_VERTICAL 0 /* only particles at the right of x = PARTIAL_THERMO_SHIFT are coupled */
#define TH_INSEGMENT 1 /* only particles in region defined by segments are coupled */
/* Plot types */
#define P_KINETIC 0 /* colors represent kinetic energy of particles */
@ -76,13 +90,15 @@
#define P_TYPE 5 /* colors represent type of particle */
#define P_DIRECTION 6 /* colors represent direction of velocity */
#define P_ANGULAR_SPEED 7 /* colors represent angular speed */
#define P_DIRECT_ENERGY 8 /* hues represent direction, luminosity represents energy */
#define P_DIFF_NEIGHB 9 /* colors represent number of neighbours of different type */
/* Color schemes */
#define C_LUM 0 /* color scheme modifies luminosity (with slow drift of hue) */
#define C_HUE 1 /* color scheme modifies hue */
#define C_PHASE 2 /* color scheme shows phase */
#define C_ONEDIM 3 /* use preset 1d color scheme (for Turbo, Viridis, Magma, Inferno, Plasma) */
#define C_ONEDIM 3 /* use preset 1d color scheme (for Turbo, Viridis, Magma, Inferno, Plasma, Twilight) */
#define C_ONEDIM_LINEAR 4 /* use preset 1d color scheme with linear scale */
/* Color palettes */
@ -121,6 +137,7 @@ typedef struct
short int thermostat; /* whether particle is coupled to thermostat */
int hashcell; /* hash cell in which particle is located */
int neighb; /* number of neighbours within given distance */
int diff_neighb; /* number of neighbours of different type */
int hash_nneighb; /* number of neighbours in hashgrid */
int hashneighbour[9*HASHMAX]; /* particle numbers of neighbours in hashgrid */
double deltax[9*HASHMAX]; /* relative position of neighbours */
@ -164,11 +181,26 @@ typedef struct
short int active; /* circle is active */
} t_obstacle;
typedef struct
{
double x1, x2, y1, y2; /* extremities of segment */
double nx, ny; /* normal vector */
double c; /* constant term in cartesian eq nx*x + ny*y = c */
double length; /* length of segment */
short int concave; /* corner is concave, to add extra repelling force */
double angle1, angle2; /* angles in which concave corners repel */
short int active; /* segment is active */
double x01, x02, y01, y02; /* initial values of extremities, in case of rotation/translation */
double nx0, ny0; /* initial normal vector */
double angle01, angle02; /* initial values of angles in which concave corners repel */
double fx, fy; /* x and y-components of force on segment */
} t_segment;
typedef struct
{
double xc, yc; /* center of circle */
} t_tracer;
int ncircles, nobstacles, counter = 0;
int ncircles, nobstacles, nsegments, counter = 0;

6
global_pdes.c
View File

@ -11,6 +11,7 @@
/* shape of domain */
#define D_NOTHING 999 /* no boundaries */
#define D_RECTANGLE 0 /* rectangular domain */
#define D_ELLIPSE 1 /* elliptical domain */
#define D_STADIUM 2 /* stadium-shaped domain */
@ -94,6 +95,9 @@
#define D_MANDELBROT_CIRCLE 26 /* Mandelbrot set with circular conductor */
#define D_VONKOCH_HEATED 27 /* von Koch snowflake in larger circle */
/* Variable index of refraction */
#define IOR_MANDELBROT 1 /* index of refraction depends on escape time in Mandelbrot set */
/* Boundary conditions */
#define BC_DIRICHLET 0 /* Dirichlet boundary conditions */
@ -114,7 +118,7 @@
#define P_ENERGY 1 /* plot energy of wave */
#define P_MIXED 2 /* plot amplitude in upper half, energy in lower half */
#define P_MEAN_ENERGY 3 /* energy averaged over time */
#define P_LOG_ENERGY 4 /* log of energy averaged over time */
#define P_LOG_ENERGY 4 /* log of energy averaged over time */
#define P_LOG_MEAN_ENERGY 5 /* log of energy averaged over time */
/* For Schrodinger equation */

286
lennardjones.c
View File

@ -34,11 +34,12 @@
#include <sys/types.h>
#include <tiffio.h> /* Sam Leffler's libtiff library. */
#include <omp.h>
#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 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 */
@ -53,24 +54,27 @@
#define YMIN -1.125
#define YMAX 1.125 /* y interval for 9/16 aspect ratio */
#define INITXMIN -1.9
#define INITXMAX 1.9 /* x interval for initial condition */
#define INITYMIN -1.0
#define INITYMAX 1.0 /* y interval for initial condition */
#define INITXMIN -0.7
#define INITXMAX 0.7 /* x interval for initial condition */
#define INITYMIN 0.1
#define INITYMAX 0.6 /* y interval for initial 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 BCXMIN -0.7
#define BCXMAX 0.7 /* x interval for boundary condition */
#define BCYMIN -0.85
#define BCYMAX 0.85 /* y interval for boundary condition */
#define OBSXMIN -2.0
#define OBSXMAX 2.0 /* x interval for motion of obstacle */
#define CIRCLE_PATTERN 8 /* pattern of circles, see list in global_ljones.c */
#define ADD_FIXED_OBSTACLES 1 /* set to 1 do add fixed circular obstacles */
#define ADD_FIXED_OBSTACLES 0 /* set to 1 do add fixed circular obstacles */
#define OBSTACLE_PATTERN 2 /* 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 2 /* pattern of repelling segments, see list in global_ljones.c */
#define TWO_TYPES 0 /* set to 1 to have two types of particles */
#define TPYE_PROPORTION 0.8 /* proportion of particles of first type */
#define SYMMETRIZE_FORCE 1 /* set to 1 to symmetrize two-particle interaction, only needed if particles are not all the same */
@ -85,23 +89,22 @@
#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 5.75 /* minimal distance in Poisson disc process, controls density of particles */
// #define PDISC_DISTANCE 2.25 /* minimal distance in Poisson disc process, controls density of particles */
#define PDISC_DISTANCE 2.75 /* minimal distance in Poisson disc process, controls density of particles */
#define PDISC_CANDIDATES 100 /* number of candidates in construction of Poisson disc process */
#define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */
#define LAMBDA 2.0 /* parameter controlling the dimensions of domain */
#define MU 0.015 /* parameter controlling radius of particles */
#define MU 0.013 /* parameter controlling radius of particles */
#define MU_B 0.0254 /* parameter controlling radius of particles of second type */
#define NPOLY 3 /* number of sides of polygon */
#define APOLY 0.125 /* angle by which to turn polygon, in units of Pi/2 */
#define APOLY 0.666666666 /* angle by which to turn polygon, in units of Pi/2 */
#define MDEPTH 4 /* depth of computation of Menger gasket */
#define MRATIO 3 /* ratio defining Menger gasket */
#define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */
#define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */
#define FOCI 1 /* set to 1 to draw focal points of ellipse */
#define NGRIDX 30 /* number of grid point for grid of disks */
#define NGRIDY 20 /* number of grid point for grid of disks */
#define NGRIDX 46 /* number of grid point for grid of disks */
#define NGRIDY 24 /* number of grid point for grid of disks */
#define EHRENFEST_RADIUS 0.9 /* radius of container for Ehrenfest urn configuration */
#define EHRENFEST_WIDTH 0.035 /* width of tube for Ehrenfest urn configuration */
@ -112,11 +115,10 @@
/* Parameters for length and speed of simulation */
#define NSTEPS 5100 /* number of frames of movie */
// #define NSTEPS 2750 /* number of frames of movie */
#define NVID 200 /* number of iterations between images displayed on screen */
#define NSTEPS 5500 /* number of frames of movie */
#define NVID 650 /* number of iterations between images displayed on screen */
#define NSEG 250 /* number of segments of boundary */
#define INITIAL_TIME 20 /* time after which to start saving frames */
#define INITIAL_TIME 10 /* time after which to start saving frames */
#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 */
@ -130,12 +132,12 @@
/* Boundary conditions, see list in global_ljones.c */
#define BOUNDARY_COND 14
#define BOUNDARY_COND 0
/* Plot type, see list in global_ljones.c */
#define PLOT 5
#define PLOT_B 4 /* plot type for second movie */
#define PLOT 0
#define PLOT_B 8 /* plot type for second movie */
#define COLOR_BONDS 1 /* set to 1 to color bonds according to length */
@ -164,7 +166,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 1.0e3 /* energy of particle with hottest color */
#define PARTICLE_EMAX 2.0e2 /* energy of particle with hottest color */
#define HUE_TYPE0 280.0 /* hue of particles of type 0 */
#define HUE_TYPE1 135.0 /* hue of particles of type 1 */
#define HUE_TYPE2 70.0 /* hue of particles of type 1 */
@ -173,7 +175,7 @@
#define RANDOM_RADIUS 0 /* set to 1 for random circle radius */
#define DT_PARTICLE 5.0e-7 /* time step for particle displacement */
#define KREPEL 12.0 /* constant in repelling force between particles */
#define EQUILIBRIUM_DIST 5.0 /* Lennard-Jones equilibrium distance */
#define EQUILIBRIUM_DIST 4.0 /* Lennard-Jones equilibrium distance */
#define EQUILIBRIUM_DIST_B 5.0 /* 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 0.0 /* damping coefficient of particles */
@ -185,16 +187,18 @@
#define OMEGA_INITIAL 10.0 /* initial angular velocity range */
#define THERMOSTAT 1 /* 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.0001 /* initial inverse temperature */
#define BETA 0.01 /* initial inverse temperature */
#define MU_XI 0.01 /* friction constant in thermostat */
#define KSPRING_BOUNDARY 5.0e9 /* confining harmonic potential outside simulation region */
#define KSPRING_OBSTACLE 5.0e8 /* harmonic potential of obstacles */
#define NBH_DIST_FACTOR 4.0 /* radius in which to count neighbours */
#define GRAVITY 0.0 /* gravity acting on all particles */
#define KSPRING_OBSTACLE 5.0e9 /* harmonic potential of obstacles */
#define NBH_DIST_FACTOR 4.5 /* radius in which to count neighbours */
#define GRAVITY 10000.0 /* gravity acting on all particles */
#define INCREASE_GRAVITY 0 /* set to 1 to increase gravity during the simulation */
#define GRAVITY_FACTOR 100.0 /* factor by which to increase gravity */
#define GRAVITY_RESTORE_TIME 750 /* time at end of simulation with gravity restored to initial value */
#define GRAVITY_INITIAL_TIME 500 /* time at start of simulation with constant gravity */
#define GRAVITY_RESTORE_TIME 1000 /* time at end of simulation with gravity restored to initial value */
#define ROTATION 0 /* set to 1 to include rotation of particles */
#define COUPLE_ANGLE_TO_THERMOSTAT 0 /* set to 1 to couple angular degrees of freedom to thermostat */
@ -212,8 +216,8 @@
#define INCREASE_BETA 0 /* set to 1 to increase BETA during simulation */
#define BETA_FACTOR 20.0 /* factor by which to change BETA during simulation */
#define N_TOSCILLATIONS 1.5 /* number of temperature oscillations in BETA schedule */
#define NO_OSCILLATION 0 /* set to 1 to have exponential BETA change only */
#define FINAL_CONSTANT_PHASE 0 /* final phase in which temperature is constant */
#define NO_OSCILLATION 1 /* set to 1 to have exponential BETA change only */
#define FINAL_CONSTANT_PHASE 2000 /* final phase in which temperature is constant */
#define DECREASE_CONTAINER_SIZE 0 /* set to 1 to decrease size of container */
#define SYMMETRIC_DECREASE 0 /* set tp 1 to decrease container symmetrically */
@ -235,6 +239,7 @@
#define N_T_AVERAGE 50 /* size of temperature averaging window */
#define MAX_PRESSURE 3.0e10 /* pressure shown in "hottest" color */
#define PARTIAL_THERMO_COUPLING 0 /* set to 1 to couple only particles to the right of obstacle to thermostat */
#define PARTIAL_THERMO_REGION 1 /* region for partial thermostat coupling (see list in global_ljones.c) */
#define PARTIAL_THERMO_SHIFT 0.5 /* distance from obstacle at the right of which particles are coupled to thermostat */
#define INCREASE_KREPEL 0 /* set to 1 to increase KREPEL during simulation */
@ -245,22 +250,40 @@
#define NPART_BOTTOM 100 /* number of particles at the bottom */
#define ADD_PARTICLES 0 /* set to 1 to add particles */
#define ADD_TIME 25 /* time at which to add first particle */
#define ADD_PERIOD 20 /* time interval between adding further particles */
#define FINAL_NOADD_PERIOD 250 /* final period where no particles are added */
#define ADD_TIME 500 /* time at which to add first particle */
#define ADD_PERIOD 250 /* time interval between adding further particles */
#define FINAL_NOADD_PERIOD 200 /* final period where no particles are added */
#define SAFETY_FACTOR 2.0 /* no particles are added at distance less than MU*SAFETY_FACTOR of other particles */
#define TRACER_PARTICLE 1 /* set to 1 to have a tracer particle */
#define TRACER_PARTICLE 0 /* set to 1 to have a tracer particle */
#define N_TRACER_PARTICLES 3 /* number of tracer particles */
#define TRAJECTORY_LENGTH 8000 /* length of recorded trajectory */
#define TRACER_PARTICLE_MASS 4.0 /* relative mass of tracer particle */
#define TRAJECTORY_WIDTH 3 /* width of tracer particle trajectory */
#define ROTATE_BOUNDARY 1 /* set to 1 to rotate the repelling segments */
#define SMOOTH_ROTATION 1 /* set to 1 to update segments at each time step (rather than at each movie frame) */
#define PERIOD_ROTATE_BOUNDARY 2500 /* period of rotating boundary */
#define ROTATE_INITIAL_TIME 0 /* initial time without rotation */
#define ROTATE_FINAL_TIME 500 /* final time without rotation */
#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 MOVE_BOUNDARY 0 /* set to 1 to move repelling segments, due to force from particles */
#define SEGMENTS_MASS 100.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 1000 /* time at which to deactivate last segment */
#define POSITION_DEPENDENT_TYPE 0 /* set to 1 to make particle type depend on initial position */
#define POSITION_Y_DEPENDENCE 0 /* set to 1 for the separation between particles to be vertical */
#define PRINT_ENTROPY 0 /* set to 1 to compute entropy */
#define REACTION_DIFFUSION 0 /* set to 1 to simulate a chemical reaction (particles may change type) */
#define RD_TYPES 3 /* number of types in reaction-diffusion equation */
#define REACTION_PROB 0.03 /* probability controlling reaction term */
#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 EHRENFEST_COPY 0 /* set to 1 to add equal number of larger particles (for Ehrenfest model) */
@ -277,8 +300,8 @@
#define FLOOR_OMEGA 1 /* set to 1 to limit particle momentum to PMAX */
#define PMAX 1000.0 /* maximal force */
#define HASHX 34 /* size of hashgrid in x direction */
#define HASHY 18 /* size of hashgrid in y direction */
#define HASHX 40 /* size of hashgrid in x direction */
#define HASHY 20 /* 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 */
@ -296,6 +319,14 @@ double ylid = 0.9; /* y coordinate of lid (for BC_RECTANGLE_LID b.c.) */
double vylid = 0.0; /* y speed coordinate of lid (for BC_RECTANGLE_LID b.c.) */
double xwall = 0.0; /* x coordinate of wall (for BC_RECTANGLE_WALL b.c.) */
double vxwall = 0.0; /* x speed of wall (for BC_RECTANGLE_WALL b.c.) */
double angular_speed = 0.0; /* angular speed of rotating segments */
double xsegments = 0.0; /* x coordinate of segments (for option MOVE_BOUNDARY) */
double ysegments = 0.0; /* y coordinate of segments (for option MOVE_BOUNDARY) */
double vxsegments = 0.0; /* vx coordinate of segments (for option MOVE_BOUNDARY) */
double vysegments = 0.0; /* vy coordinate of segments (for option MOVE_BOUNDARY) */
int thermostat_on = 1; /* thermostat switch used when VARY_THERMOSTAT is on */
#define THERMOSTAT_ON ((THERMOSTAT)&&((!VARY_THERMOSTAT)||(thermostat_on)))
#include "global_ljones.c"
#include "sub_lj.c"
@ -332,8 +363,8 @@ double temperature_schedule(int i)
if (first)
{
bexponent = log(BETA_FACTOR)/(double)(INITIAL_TIME + NSTEPS - FINAL_CONSTANT_PHASE);
omega = N_TOSCILLATIONS*DPI/(double)(INITIAL_TIME + NSTEPS - FINAL_CONSTANT_PHASE);
bexponent = log(BETA_FACTOR)/(double)(NSTEPS - FINAL_CONSTANT_PHASE);
omega = N_TOSCILLATIONS*DPI/(double)(NSTEPS - FINAL_CONSTANT_PHASE);
first = 0;
}
if (i < INITIAL_TIME) beta = BETA;
@ -411,16 +442,73 @@ double gravity_schedule(int i, int j)
{
double time, gravity;
if ((i < INITIAL_TIME)||(i > NSTEPS + INITIAL_TIME - GRAVITY_RESTORE_TIME)) return(GRAVITY);
if ((i < INITIAL_TIME + GRAVITY_INITIAL_TIME)||(i > NSTEPS + INITIAL_TIME - GRAVITY_RESTORE_TIME)) return(GRAVITY);
else
{
time = ((double)(i - INITIAL_TIME) + (double)j/(double)NVID)/(double)(NSTEPS - GRAVITY_RESTORE_TIME);
time = ((double)(i - INITIAL_TIME - GRAVITY_INITIAL_TIME)
+ (double)j/(double)NVID)/(double)(NSTEPS - GRAVITY_RESTORE_TIME - GRAVITY_INITIAL_TIME);
gravity = GRAVITY*(1.0 + time*(GRAVITY_FACTOR - 1.0));
// printf("i = %i, time = %.3lg, Gravity = %.3lg\n", i, time, gravity);
return(gravity);
}
}
double rotation_angle(double phase)
{
double omegamax = 15.0;
/* case of rotating hourglass */
while (phase > DPI) phase -= DPI;
return(phase - 0.5*sin(2.0*phase));
/* case of centrifuge */
// while (phase > DPI) phase -= DPI;
// angular_speed = 0.5*omegamax*(1.0 - cos(phase));
// return(0.5*omegamax*(phase - sin(phase)));
}
double rotation_schedule(int i)
{
double phase;
static int imin = INITIAL_TIME + ROTATE_INITIAL_TIME, imax = INITIAL_TIME + NSTEPS - ROTATE_FINAL_TIME;
if (i < imin)
{
angular_speed = 0.0;
return(0.0);
}
else
{
if (i > imax) i = imax;
phase = (DPI/(double)PERIOD_ROTATE_BOUNDARY)*(double)(i - imin);
return(rotation_angle(phase));
}
}
double rotation_schedule_smooth(int i, int j)
{
double phase;
static int imin = INITIAL_TIME + ROTATE_INITIAL_TIME, imax = INITIAL_TIME + NSTEPS - ROTATE_FINAL_TIME;
if (i < imin)
{
angular_speed = 0.0;
return(0.0);
}
else
{
if (i > imax) i = imax;
phase = (DPI/(double)PERIOD_ROTATE_BOUNDARY)*((double)(i - imin) + (double)j/(double)NVID);
return(rotation_angle(phase));
}
}
int thermostat_schedule(int i)
{
if (i < INITIAL_TIME) return(1);
else return(0);
}
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],
@ -449,7 +537,7 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA
qy[j] = particle[j].yc + 0.5*DT_PARTICLE*py[j]*particle[j].mass_inv;
qangle[j] = particle[j].angle + 0.5*DT_PARTICLE*pangle[j]*particle[j].inertia_moment_inv;
if ((THERMOSTAT)&&(particle[j].thermostat))
if ((THERMOSTAT_ON)&&(particle[j].thermostat))
{
px[j] *= exp(- 0.5*DT_PARTICLE*xi);
py[j] *= exp(- 0.5*DT_PARTICLE*xi);
@ -467,7 +555,7 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA
// *nactive++;
}
totalenergy *= DIMENSION_FACTOR; /* normalize energy to take number of degrees of freedom into account */
if (THERMOSTAT)
if (THERMOSTAT_ON)
{
a = DT_PARTICLE*(totalenergy - (double)*nactive/beta)/MU_XI;
a += SIGMA*sqrt(DT_PARTICLE)*gaussian();
@ -477,7 +565,7 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA
move = 0;
for (j=0; j<ncircles; j++) if (particle[j].active)
{
if ((THERMOSTAT)&&(particle[j].thermostat))
if ((THERMOSTAT_ON)&&(particle[j].thermostat))
{
px[j] *= exp(- 0.5*DT_PARTICLE*xi);
py[j] *= exp(- 0.5*DT_PARTICLE*xi);
@ -570,6 +658,45 @@ void evolve_wall(double fboundary)
}
void evolve_segments(t_segment segment[NMAXSEGMENTS])
{
int i, nactive = 0;
double fx = 0.0, fy = 0.0;
for (i=0; i<nsegments; i++) if (segment[i].active)
{
fx += segment[i].fx;
fy += segment[i].fy;
nactive++;
}
if (nactive > 0)
{
fx = fx/(double)nactive;
fy = fy/(double)nactive;
}
if (FLOOR_FORCE)
{
if (fx > FMAX) fx = FMAX;
else if (fx < -FMAX) fx = -FMAX;
if (fy > FMAX) fy = FMAX;
else if (fy < -FMAX) fy = -FMAX;
}
vxsegments += fx*DT_PARTICLE/(SEGMENTS_MASS);
vysegments += fy*DT_PARTICLE/(SEGMENTS_MASS);
xsegments += vxsegments*DT_PARTICLE;
ysegments += vysegments*DT_PARTICLE;
/* to avoid numerical instabilities */
if (xsegments + 1.0 > BCXMAX)
{
xsegments = BCXMAX - 1.0;
vxsegments = 0.0;
}
translate_segments(segment, xsegments, ysegments);
}
void animation()
{
double time, scale, diss, rgb[3], dissip, gradient[2], x, y, dx, dy, dt, xleft, xright, a, b,
@ -584,13 +711,15 @@ void animation()
static short int first = 1;
t_particle *particle;
t_obstacle *obstacle;
t_segment *segment;
t_tracer *trajectory;
t_hashgrid *hashgrid;
char message[100];
particle = (t_particle *)malloc(NMAXCIRCLES*sizeof(t_particle)); /* particles */
if (ADD_FIXED_OBSTACLES) obstacle = (t_obstacle *)malloc(NMAXOBSTACLES*sizeof(t_obstacle)); /* obstacles */
if (ADD_FIXED_OBSTACLES) obstacle = (t_obstacle *)malloc(NMAXOBSTACLES*sizeof(t_obstacle)); /* circular obstacles */
if (ADD_FIXED_SEGMENTS) segment = (t_segment *)malloc(NMAXSEGMENTS*sizeof(t_segment)); /* linear obstacles */
if (TRACER_PARTICLE) trajectory = (t_tracer *)malloc(TRAJECTORY_LENGTH*N_TRACER_PARTICLES*sizeof(t_tracer));
@ -611,10 +740,11 @@ void animation()
xshift = OBSTACLE_XMIN;
if (ADD_FIXED_OBSTACLES) init_obstacle_config(obstacle);
if (ADD_FIXED_SEGMENTS) init_segment_config(segment);
if (RECORD_PRESSURES) for (i=0; i<N_PRESSURES; i++) pressure[i] = 0.0;
printf("1\n");
// printf("1\n");
nactive = initialize_configuration(particle, hashgrid, obstacle, px, py, pangle, tracer_n);
@ -647,6 +777,14 @@ void animation()
xmincontainer = container_size_schedule(i);
if (SYMMETRIC_DECREASE) xmaxcontainer = -container_size_schedule(i);
}
if ((ROTATE_BOUNDARY)&&(!SMOOTH_ROTATION)) rotate_segments(segment, rotation_schedule(i));
if (VARY_THERMOSTAT)
{
thermostat_on = thermostat_schedule(i);
printf("Termostat: %i\n", thermostat_on);
}
if ((DEACTIVATE_SEGMENT)&&(i > INITIAL_TIME + SEGMENT_DEACTIVATION_TIME))
segment[nsegments-1].active = 0;
blank();
@ -662,10 +800,18 @@ void animation()
xmincontainer = obstacle_schedule_smooth(i, n);
xshift = xmincontainer;
}
if ((ROTATE_BOUNDARY)&&(SMOOTH_ROTATION)) rotate_segments(segment, rotation_schedule_smooth(i,n));
if (INCREASE_GRAVITY) gravity = gravity_schedule(i,n);
if ((BOUNDARY_COND == BC_RECTANGLE_WALL)&&(i < INITIAL_TIME + WALL_TIME)) wall = 1;
else wall = 0;
if (MOVE_BOUNDARY) for (j=0; j<nsegments; j++)
{
segment[j].fx = 0.0;
segment[j].fy = 0.0;
}
compute_relative_positions(particle, hashgrid);
update_hashgrid(particle, hashgrid, 0);
@ -680,7 +826,7 @@ void animation()
compute_particle_force(j, krepel, particle, hashgrid);
/* take care of boundary conditions */
fboundary += compute_boundary_force(j, particle, obstacle, xmincontainer, xmaxcontainer, &pleft, &pright, pressure, wall);
fboundary += compute_boundary_force(j, particle, obstacle, segment, xmincontainer, xmaxcontainer, &pleft, &pright, pressure, wall);
/* add gravity */
if (INCREASE_GRAVITY) particle[j].fy -= gravity;
@ -707,8 +853,11 @@ void animation()
if (i < INITIAL_TIME + WALL_TIME) evolve_wall(fboundary);
else xwall = 0.0;
}
if (MOVE_BOUNDARY) evolve_segments(segment);
} /* end of for (n=0; n<NVID; n++) */
if (MOVE_BOUNDARY) printf("segments position (%.3lg, %.3lg), speed (%.3lg, %.3lg)\n", xsegments, ysegments, vxsegments, vysegments);
// if ((PARTIAL_THERMO_COUPLING))
if ((PARTIAL_THERMO_COUPLING)&&(i>N_T_AVERAGE))
{
@ -770,16 +919,19 @@ void animation()
if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length);
draw_particles(particle, PLOT);
draw_container(xmincontainer, xmaxcontainer, obstacle, wall);
draw_container(xmincontainer, xmaxcontainer, obstacle, segment, wall);
/* add a particle */
if ((ADD_PARTICLES)&&((i - INITIAL_TIME - ADD_TIME + 1)%ADD_PERIOD == 0)&&(i < NSTEPS - FINAL_NOADD_PERIOD))
if ((ADD_PARTICLES)&&(i > ADD_TIME)&&((i - INITIAL_TIME - ADD_TIME + 1)%ADD_PERIOD == 0)&&(i < NSTEPS - FINAL_NOADD_PERIOD))
nadd_particle = add_particles(particle, px, py, nadd_particle);
/* case of reaction-diffusion equation */
if (REACTION_DIFFUSION) update_types(particle);
update_hashgrid(particle, hashgrid, 1);
print_parameters(beta, mean_energy, krepel, xmaxcontainer - xmincontainer,
fboundary/(double)(ncircles*NVID), 0, pressure, gravity);
fboundary/(double)(ncircles*NVID), PRINT_LEFT, pressure, gravity);
if ((BOUNDARY_COND == BC_EHRENFEST)||(BOUNDARY_COND == BC_RECTANGLE_WALL))
print_ehrenfest_parameters(particle, pleft, pright);
else if (PRINT_PARTICLE_NUMBER) print_particle_number(ncircles);
@ -790,6 +942,9 @@ void animation()
printf("Entropy 1 = %.5lg, Entropy 2 = %.5lg\n", entropy[0], entropy[1]);
print_entropy(entropy);
}
if (PRINT_OMEGA) print_omega(angular_speed);
else if (PRINT_PARTICLE_SPEEDS) print_particles_speeds(particle);
glutSwapBuffers();
@ -811,11 +966,13 @@ void animation()
{
if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length);
draw_particles(particle, PLOT_B);
draw_container(xmincontainer, xmaxcontainer, obstacle, wall);
draw_container(xmincontainer, xmaxcontainer, obstacle, segment, wall);
print_parameters(beta, mean_energy, krepel, xmaxcontainer - xmincontainer,
fboundary/(double)(ncircles*NVID), 0, pressure, gravity);
fboundary/(double)(ncircles*NVID), PRINT_LEFT, pressure, gravity);
if (BOUNDARY_COND == BC_EHRENFEST) print_ehrenfest_parameters(particle, pleft, pright);
else if (PRINT_PARTICLE_NUMBER) print_particle_number(ncircles);
if (PRINT_OMEGA) print_omega(angular_speed);
else if (PRINT_PARTICLE_SPEEDS) print_particles_speeds(particle);
glutSwapBuffers();
save_frame_lj_counter(NSTEPS + MID_FRAMES + 1 + counter);
counter++;
@ -839,11 +996,13 @@ void animation()
{
if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length);
draw_particles(particle, PLOT);
draw_container(xmincontainer, xmaxcontainer, obstacle, wall);
draw_container(xmincontainer, xmaxcontainer, obstacle, segment, wall);
print_parameters(beta, mean_energy, krepel, xmaxcontainer - xmincontainer,
fboundary/(double)(ncircles*NVID), 0, pressure, gravity);
fboundary/(double)(ncircles*NVID), PRINT_LEFT, pressure, gravity);
if (BOUNDARY_COND == BC_EHRENFEST) print_ehrenfest_parameters(particle, pleft, pright);
else if (PRINT_PARTICLE_NUMBER) print_particle_number(ncircles);
if (PRINT_OMEGA) print_omega(angular_speed);
else if (PRINT_PARTICLE_SPEEDS) print_particles_speeds(particle);
glutSwapBuffers();
}
for (i=0; i<MID_FRAMES; i++) save_frame_lj();
@ -851,11 +1010,13 @@ void animation()
{
if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length);
draw_particles(particle, PLOT_B);
draw_container(xmincontainer, xmaxcontainer, obstacle, wall);
draw_container(xmincontainer, xmaxcontainer, obstacle, segment, wall);
print_parameters(beta, mean_energy, krepel, xmaxcontainer - xmincontainer,
fboundary/(double)(ncircles*NVID), 0, pressure, gravity);
fboundary/(double)(ncircles*NVID), PRINT_LEFT, pressure, gravity);
if (BOUNDARY_COND == BC_EHRENFEST) print_ehrenfest_parameters(particle, pleft, pright);
else if (PRINT_PARTICLE_NUMBER) print_particle_number(ncircles);
if (PRINT_OMEGA) print_omega(angular_speed);
else if (PRINT_PARTICLE_SPEEDS) print_particles_speeds(particle);
glutSwapBuffers();
}
for (i=0; i<END_FRAMES; i++) save_frame_lj_counter(NSTEPS + MID_FRAMES + 1 + counter + i);
@ -869,6 +1030,7 @@ void animation()
free(particle);
if (ADD_FIXED_OBSTACLES) free(obstacle);
if (ADD_FIXED_SEGMENTS) free(segment);
if (TRACER_PARTICLE) free(trajectory);
free(hashgrid);
free(qx);
@ -883,6 +1045,12 @@ void animation()
void display(void)
{
time_t rawtime;
struct tm * timeinfo;
time(&rawtime);
timeinfo = localtime(&rawtime);
glPushMatrix();
blank();
@ -897,6 +1065,10 @@ void display(void)
glutDestroyWindow(glutGetWindow());
printf("Start local time and date: %s", asctime(timeinfo));
time(&rawtime);
timeinfo = localtime(&rawtime);
printf("Current local time and date: %s", asctime(timeinfo));
}

827
rde.c Normal file
View File

@ -0,0 +1,827 @@
/*********************************************************************************/
/* */
/* Animation of reaction-diffusion equation in a planar domain */
/* */
/* N. Berglund, January 2022 */
/* */
/* Feel free to reuse, but if doing so it would be nice to drop a */
/* line to nils.berglund@univ-orleans.fr - Thanks! */
/* */
/* compile with */
/* gcc -o rde rde.c */
/* -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lXmu -lglut -O3 -fopenmp */
/* */
/* OMP acceleration may be more effective after executing */
/* export OMP_NUM_THREADS=2 in the shell before running the program */
/* */
/* To make a video, set MOVIE to 1 and create subfolder tif_bz */
/* It may be possible to increase parameter PAUSE */
/* */
/* create movie using */
/* ffmpeg -i wave.%05d.tif -vcodec libx264 wave.mp4 */
/* */
/*********************************************************************************/
/*********************************************************************************/
/* */
/* NB: The algorithm used to simulate the wave equation is highly paralellizable */
/* One could make it much faster by using a GPU */
/* */
/*********************************************************************************/
#include <math.h>
#include <string.h>
#include <GL/glut.h>
#include <GL/glu.h>
#include <unistd.h>
#include <sys/types.h>
#include <tiffio.h> /* Sam Leffler's libtiff library. */
#include <omp.h>
#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 */
/* General geometrical parameters */
// #define WINWIDTH 1920 /* window width */
// #define WINHEIGHT 1000 /* window height */
// #define NX 480 /* number of grid points on x axis */
// #define NY 250 /* number of grid points on y axis */
// // #define NX 1920 /* number of grid points on x axis */
// // #define NY 1000 /* number of grid points on y axis */
//
// #define XMIN -2.0
// #define XMAX 2.0 /* x interval */
// #define YMIN -1.041666667
// #define YMAX 1.041666667 /* y interval for 9/16 aspect ratio */
#define WINWIDTH 1280 /* window width */
#define WINHEIGHT 720 /* window height */
// #define NX 160 /* number of grid points on x axis */
// #define NY 90 /* number of grid points on y axis */
#define NX 320 /* number of grid points on x axis */
#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 1280 /* number of grid points on x axis */
// #define NX 720 /* number of grid points on x axis */
// #define NY 720 /* number of grid points on y axis */
#define XMIN -2.0
#define XMAX 2.0 /* x interval */
#define YMIN -1.125
#define YMAX 1.125 /* y interval for 9/16 aspect ratio */
/* Choice of simulated equation */
#define RDE_EQUATION 5 /* choice of reaction term, see list in global_3d.c */
#define NFIELDS 2 /* number of fields in reaction-diffusion equation */
#define NLAPLACIANS 2 /* number of fields for which to compute Laplacian */
#define ADD_POTENTIAL 0 /* set to 1 to add a potential (for Schrodiner equation) */
#define POTENTIAL 1 /* type of potential, see list in global_3d.c */
#define JULIA_SCALE 0.5 /* scaling for Julia sets */
/* Choice of the billiard table */
#define B_DOMAIN 9 /* choice of domain shape, see list in global_pdes.c */
#define CIRCLE_PATTERN 99 /* pattern of circles, see list in global_pdes.c */
#define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */
#define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */
#define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */
#define LAMBDA 0.4 /* parameter controlling the dimensions of domain */
#define MU 0.06 /* parameter controlling the dimensions of domain */
#define NPOLY 6 /* number of sides of polygon */
#define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */
#define MDEPTH 5 /* depth of computation of Menger gasket */
#define MRATIO 5 /* ratio defining Menger gasket */
#define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */
#define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */
#define FOCI 1 /* set to 1 to draw focal points of ellipse */
#define NGRIDX 15 /* number of grid point for grid of disks */
#define NGRIDY 20 /* number of grid point for grid of disks */
#define X_SHOOTER -0.2
#define Y_SHOOTER -0.6
#define X_TARGET 0.4
#define Y_TARGET 0.7 /* shooter and target positions in laser fight */
#define ISO_XSHIFT_LEFT -1.65
#define ISO_XSHIFT_RIGHT 0.4
#define ISO_YSHIFT_LEFT -0.05
#define ISO_YSHIFT_RIGHT -0.05
#define ISO_SCALE 0.85 /* coordinates for isospectral billiards */
/* You can add more billiard tables by adapting the functions */
/* xy_in_billiard and draw_billiard in sub_wave.c */
/* Physical patameters of wave equation */
#define DT 0.00000001
#define VISCOSITY 0.1
#define RPSA 0.75 /* parameter in Rock-Paper-Scissors-type interaction */
#define RPSLZB 0.75 /* second parameter in Rock-Paper-Scissors-Lizard-Spock type interaction */
#define EPSILON 0.8 /* time scale separation */
#define DELTA 0.1 /* time scale separation */
#define FHNA 1.0 /* parameter in FHN equation */
#define FHNC -0.01 /* parameter in FHN equation */
#define K_HARMONIC 1.5 /* spring constant of harmonic potential */
#define BZQ 0.0008 /* parameter in BZ equation */
#define BZF 1.2 /* parameter in BZ equation */
#define T_OUT 2.0 /* outside temperature */
#define T_IN 0.0 /* inside temperature */
#define SPEED 0.0 /* speed of drift to the right */
#define ADD_NOISE 0 /* set to 1 to add noise, set to 2 to add noise in right half */
#define NOISE_INTENSITY 0.1 /* noise intensity */
#define CHANGE_VISCOSITY 0 /* set to 1 to change the viscosity in the course of the simulation */
#define ADJUST_INTSTEP 1 /* set to 1 to decrease integration step when viscosity increases */
#define VISCOSITY_INITIAL_TIME 10 /* initial time during which viscosity remains constant */
#define VISCOSITY_FACTOR 100.0 /* factor by which to change viscosity */
#define VISCOSITY_MAX 2.0 /* max value of viscosity beyond which NVID is increased and integration step is decrase,
for numerical stability */
#define CHANGE_RPSLZB 0 /* set to 1 to change second parameter in Rock-Paper-Scissors-Lizard-Spock equation */
#define RPSLZB_CHANGE 0.75 /* factor by which to rpslzb parameter */
#define RPSLZB_INITIAL_TIME 0 /* initial time during which rpslzb remains constant */
#define RPSLZB_FINAL_TIME 500 /* final time during which rpslzb remains constant */
/* Boundary conditions, see list in global_pdes.c */
#define B_COND 1
/* Parameters for length and speed of simulation */
#define NSTEPS 1500 /* number of frames of movie */
#define NVID 900 /* 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 */
#define MAX_DT 0.024 /* maximal value of integration step */
#define NSEG 100 /* number of segments of boundary */
#define BOUNDARY_WIDTH 1 /* width of billiard boundary */
#define PAUSE 100 /* number of frames after which to pause */
#define PSLEEP 2 /* sleep time during pause */
#define SLEEP1 2 /* initial sleeping time */
#define SLEEP2 1 /* final sleeping time */
#define INITIAL_TIME 0 /* initial still time */
#define MID_FRAMES 50 /* number of still frames between parts of two-part movie */
#define END_FRAMES 100 /* number of still frames at end of movie */
#define FADE 1 /* set to 1 to fade at end of movie */
#define PLOT_3D 1 /* controls whether plot is 2D or 3D */
/* Plot type - color scheme */
#define CPLOT 30
#define CPLOT_B 31
/* Plot type - height of 3D plot */
#define ZPLOT 30 /* z coordinate in 3D plot */
#define ZPLOT_B 30 /* z coordinate in second 3D plot */
#define AMPLITUDE_HIGH_RES 1 /* set to 1 to increase resolution of P_3D_AMPLITUDE plot */
#define SHADE_3D 1 /* set to 1 to change luminosity according to normal vector */
#define NON_DIRICHLET_BC 0 /* set to 1 to draw only facets in domain, if field is not zero on boundary */
#define WRAP_ANGLE 1 /* experimental: wrap angle to [0, 2Pi) for interpolation in angle schemes */
#define FADE_IN_OBSTACLE 0 /* set to 1 to fade color inside obstacles */
#define DRAW_OUTSIDE_GRAY 1 /* experimental - draw outside of billiard in gray */
#define PLOT_SCALE_ENERGY 0.05 /* vertical scaling in energy plot */
#define PRINT_TIME 0 /* set to 1 to print running time */
#define PRINT_VISCOSITY 0 /* set to 1 to print viscosity */
#define PRINT_RPSLZB 0 /* set to 1 to print rpslzb parameter */
#define DRAW_FIELD_LINES 0 /* set to 1 to draw field lines */
#define FIELD_LINE_WIDTH 1 /* width of field lines */
#define N_FIELD_LINES 120 /* number of field lines */
#define FIELD_LINE_FACTOR 120 /* factor controlling precision when computing origin of field lines */
#define DRAW_BILLIARD 1 /* set to 1 to draw boundary */
#define DRAW_BILLIARD_FRONT 1 /* set to 1 to draw boundary */
#define FILL_BILLIARD_COMPLEMENT 1 /* set to 1 to fill complement of billiard (for certain shapes only) */
/* 3D representation */
#define REPRESENTATION_3D 1 /* choice of 3D representation */
#define REP_AXO_3D 0 /* linear projection (axonometry) */
#define REP_PROJ_3D 1 /* projection on plane orthogonal to observer line of sight */
/* 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 BLACK 1 /* black background */
#define COLOR_SCHEME 3 /* choice of color scheme */
#define COLOR_PHASE_SHIFT 1.0 /* phase shift of color scheme, in units of Pi */
#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 20.0 /* 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) */
#define SLOPE_SCHROD_LUM 50.0 /* sensitivity of luminosity on module, for color scheme Z_ARGUMENT */
#define COLORHUE 260 /* initial hue of water color for scheme C_LUM */
#define COLORDRIFT 0.0 /* how much the color hue drifts during the whole simulation */
#define LUMMEAN 0.5 /* amplitude of luminosity variation for scheme C_LUM */
#define LUMAMP 0.3 /* amplitude of luminosity variation for scheme C_LUM */
#define HUEMEAN 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 1.0 /* scaling factor for energy log representation */
#define LOG_SHIFT 0.0
#define DRAW_COLOR_SCHEME 0 /* set to 1 to plot the color scheme */
#define COLORBAR_RANGE 0.55 /* scale of color scheme bar */
#define COLORBAR_RANGE_B 12.0 /* scale of color scheme bar for 2nd part */
#define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */
/* only for compatibility with wave_common.c */
#define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */
#define OMEGA 0.005 /* frequency of periodic excitation */
#define COURANT 0.08 /* Courant number */
#define COURANTB 0.03 /* Courant number in medium B */
#define INITIAL_AMP 0.5 /* amplitude of initial condition */
#define INITIAL_VARIANCE 0.0002 /* variance of initial condition */
#define INITIAL_WAVELENGTH 0.1 /* wavelength of initial condition */
#define VSCALE_ENERGY 200.0 /* additional scaling factor for color scheme P_3D_ENERGY */
#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 */
/* end of constants added only for compatibility with wave_common.c */
double u_3d[2] = {0.75, -0.45}; /* projections of basis vectors for REP_AXO_3D representation */
double v_3d[2] = {-0.75, -0.45};
double w_3d[2] = {0.0, 0.015};
double light[3] = {0.816496581, -0.40824829, 0.40824829}; /* vector of "light" direction for P_3D_ANGLE color scheme */
double observer[3] = {8.0, 8.0, 8.0}; /* location of observer for REP_PROJ_3D representation */
#define Z_SCALING_FACTOR 0.75 /* overall scaling factor of z axis for REP_PROJ_3D representation */
#define XY_SCALING_FACTOR 2.0 /* 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.0 /* overall x shift for REP_PROJ_3D representation */
#define YSHIFT_3D 0.05 /* overall y shift for REP_PROJ_3D representation */
/* For debugging purposes only */
#define FLOOR 1 /* set to 1 to limit wave amplitude to VMAX */
#define VMAX 10.0 /* max value of wave amplitude */
#define REFRESH_B (ZPLOT_B != ZPLOT)||(CPLOT_B != CPLOT) /* to save computing time, to be improved */
#define COMPUTE_WRAP_ANGLE ((WRAP_ANGLE)&&((cplot == Z_ANGLE_GRADIENT)||(Z_ANGLE_GRADIENTX)||(cplot == Z_ARGUMENT)||(Z_ANGLE_GRADIENTX)))
#include "global_pdes.c"
#include "sub_wave.c"
#include "wave_common.c" /* common functions for wave_billiard, wave_comparison, etc */
#include "global_3d.c" /* constants and global variables */
#include "sub_wave_3d_rde.c" /* should be later replaced by sub_wave_rde.c */
#include "sub_rde.c"
double potential(int i, int j)
/* compute potential (e.g. for Schrödinger equation) */
{
double x, y, xy[2];
ij_to_xy(i, j, xy);
x = xy[0];
y = xy[1];
switch (POTENTIAL) {
case (POT_HARMONIC):
{
return (K_HARMONIC*(x*x + y*y));
}
default:
{
return(0.0);
}
}
}
void initialize_potential(double potential_field[NX*NY])
/* initialize the potential field, e.f. for the Schrödinger equation */
{
int i, j;
#pragma omp parallel for private(i,j)
for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
potential_field[i*NY+j] = potential(i,j);
}
}
}
void evolve_wave_half(double *phi_in[NFIELDS], double *phi_out[NFIELDS], short int xy_in[NX*NY], double potential_field[NX*NY])
/* time step of field evolution */
{
int i, j, k, iplus, iminus, jplus, jminus;
double x, y, z, deltax, deltay, deltaz, rho, pot;
double *delta_phi[NLAPLACIANS];
static double invsqr3 = 0.577350269; /* 1/sqrt(3) */
for (i=0; i<NLAPLACIANS; i++) delta_phi[i] = (double *)malloc(NX*NY*sizeof(double));
/* compute the Laplacian of phi */
for (i=0; i<NLAPLACIANS; i++) compute_laplacian(phi_in[i], delta_phi[i], xy_in);
#pragma omp parallel for private(i,j,k,x,y,z,deltax,deltay,deltaz,rho)
for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
if (xy_in[i*NY+j]) switch (RDE_EQUATION){
case (E_HEAT):
{
deltax = viscosity*delta_phi[0][i*NY+j];
phi_out[0][i*NY+j] = phi_in[0][i*NY+j] + intstep*deltax;
break;
}
case (E_ALLEN_CAHN):
{
x = phi_in[0][i*NY+j];
deltax = viscosity*delta_phi[0][i*NY+j];
phi_out[0][i*NY+j] = phi_in[0][i*NY+j] + intstep*(deltax + x*(1.0-x*x));
break;
}
case (E_CAHN_HILLIARD):
{
/* TO DO */
break;
}
case (E_FHN):
{
x = phi_in[0][i*NY+j];
y = phi_in[1][i*NY+j];
deltax = viscosity*delta_phi[0][i*NY+j];
phi_out[0][i*NY+j] = phi_in[0][i*NY+j] + intstep*(deltax + x*(1.0-x*x) + y);
phi_out[1][i*NY+j] = phi_in[0][i*NY+j] + intstep*EPSILON*(- invsqr3 - FHNC - FHNA*x);
break;
}
case (E_RPS):
{
x = phi_in[0][i*NY+j];
y = phi_in[1][i*NY+j];
z = phi_in[2][i*NY+j];
rho = x + y + z;
deltax = viscosity*delta_phi[0][i*NY+j];
deltay = viscosity*delta_phi[1][i*NY+j];
deltaz = viscosity*delta_phi[2][i*NY+j];
phi_out[0][i*NY+j] = x + intstep*(deltax + x*(1.0 - rho - RPSA*y));
phi_out[1][i*NY+j] = y + intstep*(deltay + y*(1.0 - rho - RPSA*z));
phi_out[2][i*NY+j] = z + intstep*(deltaz + z*(1.0 - rho - RPSA*x));
break;
}
case (E_RPSLZ):
{
rho = 0.0;
for (k=0; k<5; k++) rho += phi_in[k][i*NY+j];
for (k=0; k<5; k++)
{
x = phi_in[k][i*NY+j];
y = phi_in[(k+1)%5][i*NY+j];
z = phi_in[(k+3)%5][i*NY+j];
phi_out[k][i*NY+j] = x + intstep*(delta_phi[k][i*NY+j] + x*(1.0 - rho - RPSA*y - rpslzb*z));
}
break;
}
case (E_SCHRODINGER):
{
phi_out[0][i*NY+j] = phi_in[0][i*NY+j] - intstep*delta_phi[1][i*NY+j];
phi_out[1][i*NY+j] = phi_in[1][i*NY+j] + intstep*delta_phi[0][i*NY+j];
if (ADD_POTENTIAL)
{
pot = potential_field[i*NY+j];
phi_out[0][i*NY+j] += intstep*pot*phi_in[1][i*NY+j];
phi_out[1][i*NY+j] -= intstep*pot*phi_in[0][i*NY+j];
}
}
}
}
}
if (FLOOR) for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
if (xy_in[i*NY+j] != 0) for (k=0; k<NFIELDS; k++)
{
if (phi_out[k][i*NY+j] > VMAX) phi_out[k][i*NY+j] = VMAX;
if (phi_out[k][i*NY+j] < -VMAX) phi_out[k][i*NY+j] = -VMAX;
}
}
}
for (i=0; i<NLAPLACIANS; i++) free(delta_phi[i]);
}
void evolve_wave(double *phi[NFIELDS], double *phi_tmp[NFIELDS], short int xy_in[NX*NY], double potential_field[NX*NY])
/* time step of field evolution */
{
evolve_wave_half(phi, phi_tmp, xy_in, potential_field);
evolve_wave_half(phi_tmp, phi, xy_in, potential_field);
}
void print_level(int level)
{
double pos[2];
char message[50];
glColor3f(1.0, 1.0, 1.0);
sprintf(message, "Level %i", level);
xy_to_pos(XMIN + 0.1, YMAX - 0.2, pos);
write_text(pos[0], pos[1], message);
}
void print_parameters(double time, short int left, double viscosity)
{
char message[100];
double density, hue, rgb[3], logratio, y, pos[2];
static double xbox, xtext, boxwidth, boxheight;
static int first = 1;
if (first)
{
if (WINWIDTH > 1280)
{
boxheight = 0.035;
boxwidth = 0.21;
if (left)
{
xbox = XMIN + 0.4;
xtext = XMIN + 0.2;
}
else
{
xbox = XMAX - 0.39;
xtext = XMAX - 0.55;
}
}
else
{
boxwidth = 0.3;
boxheight = 0.05;
if (left)
{
xbox = XMIN + 0.4;
xtext = XMIN + 0.1;
}
else
{
xbox = XMAX - 0.39;
xtext = XMAX - 0.61;
}
}
first = 0;
}
y = YMAX - 0.1;
erase_area_hsl(xbox, y + 0.02, boxwidth, boxheight, 0.0, 0.9, 0.0);
glColor3f(1.0, 1.0, 1.0);
if (PRINT_TIME) sprintf(message, "Time %.3f", time);
else if (PRINT_VISCOSITY) sprintf(message, "Viscosity %.3f", viscosity);
else if (PRINT_RPSLZB) sprintf(message, "b = %.3f", rpslzb);
if (PLOT_3D) write_text(xtext, y, message);
else
{
xy_to_pos(xtext, y, pos);
write_text(pos[0], pos[1], message);
}
}
void draw_color_bar_palette(int plot, double range, int palette, int fade, double fade_value)
{
double width = 0.14;
// double width = 0.2;
if (ROTATE_COLOR_SCHEME)
draw_color_scheme_palette_3d(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range, palette, fade, fade_value);
else
draw_color_scheme_palette_3d(XMAX - 1.5*width, YMIN + 0.1, XMAX - 0.5*width, YMAX - 0.1, plot, -range, range, palette, fade, fade_value);
}
double viscosity_schedule(int i)
{
double ratio;
if (i < VISCOSITY_INITIAL_TIME) return (VISCOSITY);
else
{
ratio = (double)(i - VISCOSITY_INITIAL_TIME)/(double)(NSTEPS - VISCOSITY_INITIAL_TIME);
return (VISCOSITY*(1.0 + ratio*(VISCOSITY_FACTOR - 1.0)));
}
}
double rpslzb_schedule(int i)
{
double ratio;
if (i < RPSLZB_INITIAL_TIME) return (RPSLZB);
else if (i > NSTEPS - RPSLZB_FINAL_TIME) return(RPSLZB - RPSLZB_CHANGE);
else
{
ratio = (double)(i - RPSLZB_INITIAL_TIME)/(double)(NSTEPS - RPSLZB_INITIAL_TIME - RPSLZB_FINAL_TIME);
return (RPSLZB - ratio*RPSLZB_CHANGE);
}
}
void animation()
{
double time = 0.0, scale, dx, var, jangle, cosj, sinj, sqrintstep, intstep0, viscosity_printed, fade_value;
double *phi[NFIELDS], *phi_tmp[NFIELDS], *potential_field;
short int *xy_in;
int i, j, k, s, nvid, field;
static int counter = 0;
t_rde *rde;
/* Since NX and NY are big, it seemed wiser to use some memory allocation here */
for (i=0; i<NFIELDS; i++)
{
phi[i] = (double *)malloc(NX*NY*sizeof(double));
phi_tmp[i] = (double *)malloc(NX*NY*sizeof(double));
}
xy_in = (short int *)malloc(NX*NY*sizeof(short int));
rde = (t_rde *)malloc(NX*NY*sizeof(t_rde));
if (ADD_POTENTIAL) potential_field = (double *)malloc(NX*NY*sizeof(double));
npolyline = init_polyline(MDEPTH, polyline);
for (i=0; i<npolyline; i++) printf("vertex %i: (%.3f, %.3f)\n", i, polyline[i].x, polyline[i].y);
dx = (XMAX-XMIN)/((double)NX);
intstep = DT/(dx*dx);
intstep0 = intstep;
intstep1 = DT/dx;
viscosity = VISCOSITY;
sqrintstep = sqrt(intstep*(double)NVID);
printf("Integration step %.3lg\n", intstep);
/* initialize field */
// init_random(0.5, 0.4, phi, xy_in);
// init_random(0.0, 0.4, phi, xy_in);
// init_gaussian(x, y, mean, amplitude, scalex, phi, xy_in)
init_coherent_state(-0.75, 0.0, 5.0, -7.0, 0.2, phi, xy_in);
init_cfield_rde(phi, xy_in, CPLOT, rde, 0);
if (PLOT_3D) init_zfield_rde(phi, xy_in, ZPLOT, rde, 0);
if (DOUBLE_MOVIE)
{
init_cfield_rde(phi, xy_in, CPLOT_B, rde, 1);
if (PLOT_3D) init_zfield_rde(phi, xy_in, ZPLOT_B, rde, 1);
}
blank();
glColor3f(0.0, 0.0, 0.0);
glutSwapBuffers();
printf("Drawing wave\n");
draw_wave_rde(0, phi, xy_in, rde, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1);
// draw_billiard();
if ((PRINT_TIME)||(PRINT_VISCOSITY)||(PRINT_RPSLZB)) print_parameters(time, 0, VISCOSITY);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 0, 1.0);
glutSwapBuffers();
sleep(SLEEP1);
// printf("Saving frame %i\n", i);
if (MOVIE) for (i=0; i<INITIAL_TIME; i++) save_frame();
for (i=0; i<=NSTEPS; i++)
{
nvid = NVID;
if (CHANGE_VISCOSITY)
{
viscosity = viscosity_schedule(i);
viscosity_printed = viscosity;
printf("Viscosity = %.3lg\n", viscosity);
if ((ADJUST_INTSTEP)&&(viscosity > VISCOSITY_MAX))
{
nvid = (int)((double)NVID*viscosity/VISCOSITY_MAX);
// viscosity = VISCOSITY_MAX;
intstep = intstep0*VISCOSITY_MAX/viscosity;
printf("Nvid = %i, intstep = %.3lg\n", nvid, intstep);
}
}
if (CHANGE_RPSLZB) rpslzb = rpslzb_schedule(i);
printf("Drawing wave %i\n", i);
draw_wave_rde(0, phi, xy_in, rde, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1);
// nvid = (int)((double)NVID*(1.0 + (ACCELERATION_FACTOR - 1.0)*(double)i/(double)NSTEPS));
/* increase integration step */
// intstep = intstep0*exp(log(DT_ACCELERATION_FACTOR)*(double)i/(double)NSTEPS);
// if (intstep > MAX_DT)
// {
// nvid *= intstep/MAX_DT;
// intstep = MAX_DT;
// }
// printf("Steps per frame: %i\n", nvid);
// printf("Integration step %.5lg\n", intstep);
printf("Evolving wave\n");
for (j=0; j<nvid; j++) evolve_wave(phi, phi_tmp, xy_in, potential_field);
for (j=0; j<NFIELDS; j++) printf("field[%i] = %.3lg\t", j, phi[j][0]);
printf("\n");
if (ADD_NOISE == 1)
{
for (field=0; field<=NFIELDS; field++)
for (j=0; j<NX; j++)
for (k=0; k<NY; k++)
phi[field][j*NY+k] += sqrintstep*NOISE_INTENSITY*gaussian();
}
else if (ADD_NOISE == 2)
{
for (field=0; field<=NFIELDS; field++)
for (j=NX/2; j<NX; j++)
for (k=0; k<NY; k++)
phi[field][j*NY+k] += sqrintstep*NOISE_INTENSITY*gaussian();
}
time += nvid*intstep;
// draw_billiard();
if ((PRINT_TIME)||(PRINT_VISCOSITY)||(PRINT_RPSLZB)) print_parameters(time, 0, viscosity_printed);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 0, 1.0);
// print_level(MDEPTH);
// print_Julia_parameters(i);
glutSwapBuffers();
/* modify Julia set */
// set_Julia_parameters(i, phi, xy_in);
if (MOVIE)
{
printf("Saving frame %i\n", i);
save_frame();
if ((i >= INITIAL_TIME)&&(DOUBLE_MOVIE))
{
draw_wave_rde(1, phi, xy_in, rde, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 0, 1.0, REFRESH_B);
// draw_billiard();
if ((PRINT_TIME)||(PRINT_VISCOSITY)||(PRINT_RPSLZB)) print_parameters(time, 0, viscosity_printed);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0);
glutSwapBuffers();
save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter);
counter++;
}
/* it seems that saving too many files too fast can cause trouble with the file system */
/* so this is to make a pause from time to time - parameter PAUSE may need adjusting */
if (i % PAUSE == PAUSE - 1)
{
printf("Making a short pause\n");
sleep(PSLEEP);
s = system("mv wave*.tif tif_bz/");
}
}
else printf("Computing frame %i\n", i);
}
if (MOVIE)
{
if (DOUBLE_MOVIE)
{
draw_wave_rde(0, phi, xy_in, rde, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1);
// draw_billiard();
if ((PRINT_TIME)||(PRINT_VISCOSITY)||(PRINT_RPSLZB)) print_parameters(time, 0, viscosity_printed);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 0, 1.0);
glutSwapBuffers();
if (!FADE) for (i=0; i<MID_FRAMES; i++) save_frame();
else for (i=0; i<MID_FRAMES; i++)
{
fade_value = 1.0 - (double)i/(double)MID_FRAMES;
draw_wave_rde(0, phi, xy_in, rde, ZPLOT, CPLOT, COLOR_PALETTE, 1, fade_value, 0);
// draw_billiard();
if ((PRINT_TIME)||(PRINT_VISCOSITY)||(PRINT_RPSLZB)) print_parameters(time, 0, viscosity_printed);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 1, fade_value);
glutSwapBuffers();
save_frame_counter(NSTEPS + i + 1);
}
draw_wave_rde(1, phi, xy_in, rde, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 0, 1.0, REFRESH_B);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0);
glutSwapBuffers();
if (!FADE) for (i=0; i<END_FRAMES; i++) save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter + i);
else for (i=0; i<END_FRAMES; i++)
{
fade_value = 1.0 - (double)i/(double)END_FRAMES;
draw_wave_rde(1, phi, xy_in, rde, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 1, fade_value, 0);
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);
}
}
else
{
if (!FADE) for (i=0; i<END_FRAMES; i++) save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter + i);
else for (i=0; i<END_FRAMES; i++)
{
fade_value = 1.0 - (double)i/(double)END_FRAMES;
draw_wave_rde(0, phi, xy_in, rde, ZPLOT, CPLOT, COLOR_PALETTE, 1, fade_value, 0);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 1, fade_value);
glutSwapBuffers();
save_frame_counter(NSTEPS + 1 + counter + i);
}
}
s = system("mv wave*.tif tif_bz/");
}
for (i=0; i<NFIELDS; i++)
{
free(phi[i]);
free(phi_tmp[i]);
}
free(xy_in);
if (ADD_POTENTIAL) free(potential_field);
printf("Time %.5lg\n", time);
}
void display(void)
{
time_t rawtime;
struct tm * timeinfo;
time(&rawtime);
timeinfo = localtime(&rawtime);
glPushMatrix();
blank();
glutSwapBuffers();
blank();
glutSwapBuffers();
animation();
sleep(SLEEP2);
glPopMatrix();
glutDestroyWindow(glutGetWindow());
printf("Start local time and date: %s", asctime(timeinfo));
time(&rawtime);
timeinfo = localtime(&rawtime);
printf("Current local time and date: %s", asctime(timeinfo));
}
int main(int argc, char** argv)
{
glutInit(&argc, argv);
glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH);
glutInitWindowSize(WINWIDTH,WINHEIGHT);
glutCreateWindow("FitzHugh-Nagumo equation in a planar domain");
if (PLOT_3D) init_3d();
else init();
glutDisplayFunc(display);
glutMainLoop();
return 0;
}

32
sub_hashgrid.c
View File

@ -349,22 +349,22 @@ void init_hashgrid(t_hashgrid hashgrid[HASHX*HASHY])
default: /* do nothing */;
}
for (i=0; i<HASHX; i++)
{
for (j=0; j<HASHY; j++)
{
for (k=0; k<hashgrid[mhash(i,j)].nneighb; k++)
{
m = hashgrid[mhash(i,j)].neighbour[k];
p = m/HASHY;
q = m%HASHY;
if (vabs((double)(p-i)) + vabs((double)(q-j)) > 2.0)
printf("Grid cell (%i, %i) - neighbour %i = %i = (%i, %i)\n", i, j, k, m, p, q);
}
// sleep(1);
}
// sleep(1);
}
// for (i=0; i<HASHX; i++)
// {
// for (j=0; j<HASHY; j++)
// {
// for (k=0; k<hashgrid[mhash(i,j)].nneighb; k++)
// {
// m = hashgrid[mhash(i,j)].neighbour[k];
// p = m/HASHY;
// q = m%HASHY;
// if (vabs((double)(p-i)) + vabs((double)(q-j)) > 2.0)
// printf("Grid cell (%i, %i) - neighbour %i = %i = (%i, %i)\n", i, j, k, m, p, q);
// }
// // sleep(1);
// }
// // sleep(1);
// }
sleep(1);
}

734
sub_lj.c
View File

@ -837,7 +837,7 @@ void init_people_config(t_person people[NMAXCIRCLES])
}
void init_obstacle_config(t_obstacle obstacle[NMAXOBSTACLES])
/* initialise particle configuration */
/* initialise circular obstacle configuration */
{
int i, j, n;
double x, y, dx, dy;
@ -898,6 +898,354 @@ void init_obstacle_config(t_obstacle obstacle[NMAXOBSTACLES])
}
void init_segment_config(t_segment segment[NMAXSEGMENTS])
/* initialise linear obstacle configuration */
{
int i, cycle = 0, iminus, iplus;
double angle, angle2, dx, width, height, a, b, length;
switch (SEGMENT_PATTERN) {
case (S_RECTANGLE):
{
segment[0].x1 = BCXMIN;
segment[0].y1 = BCYMAX;
segment[1].x1 = BCXMIN;
segment[1].y1 = BCYMIN;
segment[2].x1 = BCXMAX;
segment[2].y1 = BCYMIN;
segment[3].x1 = BCXMAX;
segment[3].y1 = BCYMAX;
cycle = 1;
nsegments = 4;
for (i=0; i<nsegments; i++) segment[i].concave = 0;
break;
}
case (S_CUP):
{
angle = APOLY*PID;
dx = (BCYMAX - BCYMIN)/tan(angle);
segment[0].x1 = BCXMIN;
segment[0].y1 = BCYMAX;
segment[1].x1 = BCXMIN + dx;
segment[1].y1 = BCYMIN;
segment[2].x1 = BCXMAX - dx;
segment[2].y1 = BCYMIN;
segment[3].x1 = BCXMAX;
segment[3].y1 = BCYMAX;
cycle = 1;
nsegments = 4;
for (i=0; i<nsegments; i++) segment[i].concave = 0;
break;
}
case (S_HOURGLASS):
{
angle = APOLY*PID;
width = 2.5*MU;
height = 2.5*MU;
segment[0].x1 = BCXMIN;
segment[0].y1 = BCYMAX;
segment[0].concave = 0;
segment[1].x1 = -width;
segment[1].y1 = height;
segment[1].concave = 1;
segment[2].x1 = -width;
segment[2].y1 = -height;
segment[2].concave = 1;
segment[3].x1 = BCXMIN;
segment[3].y1 = BCYMIN;
segment[3].concave = 0;
segment[4].x1 = BCXMAX;
segment[4].y1 = BCYMIN;
segment[4].concave = 0;
segment[5].x1 = width;
segment[5].y1 = -height;
segment[5].concave = 1;
segment[6].x1 = width;
segment[6].y1 = height;
segment[6].concave = 1;
segment[7].x1 = BCXMAX;
segment[7].y1 = BCYMAX;
segment[7].concave = 0;
cycle = 1;
nsegments = 8;
break;
}
case (S_PENTA):
{
height = 0.5*(BCYMAX - BCYMIN);
width = height/sqrt(3.0);
segment[0].x1 = BCXMIN;
segment[0].y1 = 0.5*(BCYMIN + BCYMAX);
segment[1].x1 = BCXMIN + width;
segment[1].y1 = BCYMIN;
segment[2].x1 = BCXMAX;
segment[2].y1 = BCYMIN;
segment[3].x1 = BCXMAX;
segment[3].y1 = BCYMAX;
segment[4].x1 = BCXMIN + width;
segment[4].y1 = BCYMAX;
cycle = 1;
nsegments = 5;
for (i=0; i<nsegments; i++) segment[i].concave = 0;
break;
}
case (S_CENTRIFUGE):
{
angle = DPI/(double)NPOLY;
for (i=0; i<NPOLY; i++)
{
segment[i*4].x1 = LAMBDA*cos(((double)i + 0.02)*angle);
segment[i*4].y1 = LAMBDA*sin(((double)i + 0.02)*angle);
segment[i*4].concave = 1;
segment[i*4 + 1].x1 = cos(((double)i + 0.05)*angle);
segment[i*4 + 1].y1 = sin(((double)i + 0.05)*angle);
segment[i*4 + 1].concave = 0;
segment[i*4 + 2].x1 = cos(((double)i + 0.95)*angle);
segment[i*4 + 2].y1 = sin(((double)i + 0.95)*angle);
segment[i*4 + 2].concave = 0;
segment[i*4 + 3].x1 = LAMBDA*cos(((double)i + 0.98)*angle);
segment[i*4 + 3].y1 = LAMBDA*sin(((double)i + 0.98)*angle);
segment[i*4 + 3].concave = 1;
}
cycle = 1;
nsegments = 4*NPOLY;
break;
}
case (S_POLY_ELLIPSE):
{
angle = DPI/(double)NPOLY;
for (i=0; i<NPOLY; i++)
{
segment[i].x1 = -cos(((double)i + 0.5)*angle);
segment[i].y1 = -LAMBDA*sin(((double)i + 0.5)*angle);
segment[i].concave = 0;
}
segment[0].concave = 1;
segment[NPOLY-1].concave = 1;
for (i=0; i<NPOLY; i++)
{
segment[NPOLY+i].x1 = 1.05*segment[NPOLY-1-i].x1;
segment[NPOLY+i].y1 = 1.05*segment[NPOLY-1-i].y1;
segment[NPOLY+i].concave = 1;
}
cycle = 1;
nsegments = 2*NPOLY;
break;
}
default:
{
printf("Function init_segment_config not defined for this pattern \n");
}
}
if (cycle) for (i=0; i<nsegments; i++)
{
segment[i].x2 = segment[(i+1)%(nsegments)].x1;
segment[i].y2 = segment[(i+1)%(nsegments)].y1;
}
/* add one segment for S_POLY_ELLIPSE configuration */
if (SEGMENT_PATTERN == S_POLY_ELLIPSE)
{
segment[nsegments].x1 = -cos(((double)i + 0.5)*angle);
segment[nsegments].y1 = LAMBDA*sin(((double)i + 0.5)*angle);
segment[nsegments].x2 = -cos(((double)i + 0.5)*angle);
segment[nsegments].y2 = -LAMBDA*sin(((double)i + 0.5)*angle);
nsegments++;
}
/* activate all segments */
for (i=0; i<nsegments; i++) segment[i].active = 1;
/* compute parameters for slope and normal of segments */
for (i=0; i<nsegments; i++)
{
a = segment[i].y1 - segment[i].y2;
b = segment[i].x2 - segment[i].x1;
length = module2(a, b);
segment[i].nx = a/length;
segment[i].ny = b/length;
segment[i].c = segment[i].nx*segment[i].x1 + segment[i].ny*segment[i].y1;
segment[i].length = length;
segment[i].fx = 0.0;
segment[i].fy = 0.0;
}
/* deal with concave corners */
for (i=0; i<nsegments; i++) if (segment[i].concave)
{
iminus = i-1; if (iminus == 0) iminus = nsegments - 1;
iplus = i+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;
}
/* make copy of initial values in case of rotation/translation */
if ((ROTATE_BOUNDARY)||(MOVE_BOUNDARY)) for (i=0; i<nsegments; i++)
{
segment[i].x01 = segment[i].x1;
segment[i].x02 = segment[i].x2;
segment[i].y01 = segment[i].y1;
segment[i].y02 = segment[i].y2;
segment[i].nx0 = segment[i].nx;
segment[i].ny0 = segment[i].ny;
segment[i].angle01 = segment[i].angle1;
segment[i].angle02 = segment[i].angle2;
}
for (i=0; i<nsegments; i++)
{
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);
}
}
int in_segment_region(double x, double y)
/* returns 1 if (x,y) is inside region delimited by obstacle segments */
{
double angle, dx, height, width, theta;
if (x >= BCXMAX) return(0);
if (x <= BCXMIN) return(0);
if (y >= BCYMAX) return(0);
if (y <= BCYMIN) return(0);
switch (SEGMENT_PATTERN) {
case (S_CUP):
{
angle = APOLY*PID;
dx = (BCYMAX - BCYMIN)/tan(angle);
if (y < BCYMAX - (BCYMAX - BCYMIN)*(x - BCXMIN)/dx) return(0);
if (y < BCYMAX - (BCYMAX - BCYMIN)*(BCXMAX - x)/dx) return(0);
}
case (S_HOURGLASS):
{
angle = APOLY*PID;
width = 2.5*MU;
height = 2.5*MU;
x = vabs(x);
y = vabs(y);
if ((x >= width)&&(x - width >= (y - height)*(BCXMAX - width)/(BCYMAX - height))) return(0);
return(1);
}
case (S_PENTA):
{
height = 0.5*(BCYMAX - BCYMIN);
width = height/sqrt(3.0);
if (y < BCYMIN + height*(1.0 - (x - BCXMIN)/width)) return(0);
if (y > BCYMAX - height*(1.0 - (x - BCXMIN)/width)) return(0);
return(1);
}
case (S_CENTRIFUGE):
{
angle = argument(x,y);
theta = DPI/(double)NPOLY;
while (angle > theta) angle -= theta;
while (angle < 0.0) angle += theta;
if (angle < 0.1) return(0);
if (angle > 0.9) return(0);
return(1);
}
case (S_POLY_ELLIPSE):
{
if (x*x + y*y/(LAMBDA*LAMBDA) < 0.95) return(1);
else return(0);
}
default: return(1);
}
}
void rotate_segments(t_segment segment[NMAXSEGMENTS], double angle)
/* rotates the repelling segments by given angle */
{
int i;
double ca, sa;
ca = cos(angle);
sa = sin(angle);
for (i=0; i<nsegments; i++)
{
segment[i].x1 = ca*segment[i].x01 - sa*segment[i].y01;
segment[i].y1 = sa*segment[i].x01 + ca*segment[i].y01;
segment[i].x2 = ca*segment[i].x02 - sa*segment[i].y02;
segment[i].y2 = sa*segment[i].x02 + ca*segment[i].y02;
segment[i].nx = ca*segment[i].nx0 - sa*segment[i].ny0;
segment[i].ny = sa*segment[i].nx0 + ca*segment[i].ny0;
if (segment[i].concave)
{
segment[i].angle1 = segment[i].angle01 + angle;
segment[i].angle2 = segment[i].angle02 + angle;
while (segment[i].angle1 > DPI)
{
segment[i].angle1 -= DPI;
segment[i].angle2 -= DPI;
}
}
}
}
void translate_segments(t_segment segment[NMAXSEGMENTS], double deltax, double deltay)
/* rotates the repelling segments by given angle */
{
int i;
for (i=0; i<nsegments; i++)
{
segment[i].x1 = segment[i].x01 + deltax;
segment[i].y1 = segment[i].y01 + deltay;
segment[i].x2 = segment[i].x02 + deltax;
segment[i].y2 = segment[i].y02 + deltay;
segment[i].c = segment[i].nx*segment[i].x1 + segment[i].ny*segment[i].y1;
}
}
/* Computation of interaction force */
double lennard_jones_force(double r, t_particle particle)
@ -1637,6 +1985,7 @@ int add_particle(double x, double y, double vx, double vy, double mass, short in
particle[i].radius = MU*sqrt(mass);
particle[i].active = 1;
particle[i].neighb = 0;
particle[i].diff_neighb = 0;
particle[i].thermostat = 1;
particle[i].energy = 0.0;
@ -1754,97 +2103,6 @@ void draw_one_particle_links(t_particle particle)
if (length < 1.0*NBH_DIST_FACTOR) draw_line(x1, y1, x2, y2);
}
// {
// draw_line(x1, y1, x2, y2);
/* in case of periodic b.c., draw translates of particles */
// if (PERIODIC_BC)
// {
// for (i=-1; i<2; i++)
// for (j=-1; j<2; j++)
// draw_line(x1 + (double)i*(BCXMAX - BCXMIN), y1 + (double)j*(BCYMAX - BCYMIN), x2 + (double)i*(BCXMAX - BCXMIN), y2 + (double)j*(BCYMAX - BCYMIN));
// }
// else if (BOUNDARY_COND == BC_KLEIN)
// {
// x1 = particle[j].xc;
// y1 = particle[j].yc;
//
// for (i=-2; i<3; i++)
// {
// if (vabs(i) == 1) sign = -1.0;
// else sign = 1.0;
// angle1 = angle*sign;
// for (k=-1; k<2; k++)
// draw_one_particle(particle[j], x1 + (double)i*(BCXMAX - BCXMIN), sign*(y1 + (double)k*(BCYMAX - BCYMIN)),
// radius, angle1, nsides, width, rgb);
// }
// }
// else if (BOUNDARY_COND == BC_BOY)
// {
// x1 = particle[j].xc;
// y1 = particle[j].yc;
//
// for (i=-1; i<2; i++) for (k=-1; k<2; k++)
// {
// if (vabs(i) == 1) sign = -1.0;
// else sign = 1.0;
// if (vabs(k) == 1) signy = -1.0;
// else signy = 1.0;
// if (signy == 1.0) angle1 = angle*sign;
// else angle1 = PI - angle;
// if (sign == -1.0) draw_one_particle(particle[j], signy*(x1 + (double)i*(BCXMAX - BCXMIN)),
// sign*(y1 + (double)k*(BCYMAX - BCYMIN)), radius, angle1, nsides, width, rgbx);
// else if (signy == -1.0) draw_one_particle(particle[j], signy*(x1 + (double)i*(BCXMAX - BCXMIN)),
// sign*(y1 + (double)k*(BCYMAX - BCYMIN)), radius, angle1, nsides, width, rgby);
// else draw_one_particle(particle[j], signy*(x1 + (double)i*(BCXMAX - BCXMIN)),
// sign*(y1 + (double)k*(BCYMAX - BCYMIN)), radius, angle1, nsides, width, rgb);
// }
// }
// else if (BOUNDARY_COND == BC_GENUS_TWO)
// {
// if (x1 < 0.0) periody = BCYMAX - BCYMIN;
// else periody = 0.5*(BCYMAX - BCYMIN);
//
// if (y1 < 0.0) periodx = BCXMAX - BCXMIN;
// else periodx = 0.5*(BCXMAX - BCXMIN);
// if ((x1 < 0.0)&&(y1 < 0.0))
// for (i=-1; i<2; i++)
// for (j=-1; j<2; j++)
// {
// xt1 = x1 + (double)i*periodx;
// yt1 = y1 + (double)j*periody;
// xt2 = x2 + (double)i*periodx;
// yt2 = y2 + (double)j*periody;
// draw_line(xt1, yt1, xt2, yt2);
// }
// else if ((x1 < 0.0)&&(y1 >= 0.0))
// for (i=-1; i<2; i++)
// for (j=-1; j<2; j++)
// {
// xt1 = x1 + (double)i*periodx;
// yt1 = y1 + (double)j*periody;
// xt2 = x2 + (double)i*periodx;
// yt2 = y2 + (double)j*periody;
// draw_line(xt1, yt1, xt2, yt2);
// }
// else if ((x1 >= 0.0)&&(y1 < 0.0))
// for (i=-1; i<2; i++)
// for (j=-1; j<2; j++)
// {
// xt1 = x1 + (double)i*periodx;
// yt1 = y1 + (double)j*periody;
// xt2 = x2 + (double)i*periodx;
// yt2 = y2 + (double)j*periody;
// draw_line(xt1, yt1, xt2, yt2);
// }
// }
// }
// }
// sprintf(message, "%i - %i", particle[j].hash_nneighb, particle[j].hashcell);
// write_text(particle[j].xc, particle[j].yc, message);
// }
}
void draw_one_particle(t_particle particle, double xc, double yc, double radius, double angle, int nsides, double width, double rgb[3])
@ -1969,7 +2227,7 @@ void draw_trajectory(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTICLES],
void draw_particles(t_particle particle[NMAXCIRCLES], int plot)
{
int i, j, k, m, width, nnbg, nsides;
double ej, hue, huex, huey, rgb[3], rgbx[3], rgby[3], radius, x1, y1, x2, y2, angle, ca, sa, length, linkcolor, sign = 1.0, angle1, signy = 1.0, periodx, periody, x, y;
double ej, hue, huex, huey, rgb[3], rgbx[3], rgby[3], radius, x1, y1, x2, y2, angle, ca, sa, length, linkcolor, sign = 1.0, angle1, signy = 1.0, periodx, periody, x, y, lum;
char message[100];
if (!TRACER_PARTICLE) blank();
@ -1980,32 +2238,6 @@ void draw_particles(t_particle particle[NMAXCIRCLES], int plot)
{
glLineWidth(LINK_WIDTH);
for (j=0; j<ncircles; j++) if (particle[j].active) draw_one_particle_links(particle[j]);
// if (0) {
// // radius = particle[j].radius;
// for (k = 0; k < particle[j].hash_nneighb; k++)
// {
// x1 = particle[j].xc;
// if (CENTER_VIEW_ON_OBSTACLE) x1 -= xshift;
// y1 = particle[j].yc;
// x2 = x1 + particle[j].deltax[k];
// y2 = y1 + particle[j].deltay[k];
//
// length = module2(particle[j].deltax[k], particle[j].deltay[k])/particle[j].radius;
//
// // if (COLOR_BONDS)
// // {
// // if (length < 1.5) linkcolor = 1.0;
// // else linkcolor = 1.0 - 0.75*(length - 1.5)/(NBH_DIST_FACTOR - 1.5);
// // glColor3f(linkcolor, linkcolor, linkcolor);
// // }
//
// if (length < 2.0*NBH_DIST_FACTOR)
// draw_line(x1, y1, x2, y2);
// }
//
// // sprintf(message, "%i - %i", particle[j].hash_nneighb, particle[j].hashcell);
// // write_text(particle[j].xc, particle[j].yc, message);
// }
}
/* determine particle color and size */
@ -2082,6 +2314,18 @@ void draw_particles(t_particle particle[NMAXCIRCLES], int plot)
width = BOUNDARY_WIDTH;
break;
}
case (P_DIRECT_ENERGY):
{
hue = argument(particle[j].vx, particle[j].vy);
if (hue > DPI) hue -= DPI;
if (hue < 0.0) hue += DPI;
hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*hue/DPI;
if (particle[j].energy < 0.1*PARTICLE_EMAX) lum = 10.0*particle[j].energy/PARTICLE_EMAX;
else lum = 1.0;
radius = particle[j].radius;
width = BOUNDARY_WIDTH;
break;
}
case (P_ANGULAR_SPEED):
{
hue = 160.0*(1.0 + tanh(SLOPE*particle[j].omega));
@ -2090,6 +2334,15 @@ void draw_particles(t_particle particle[NMAXCIRCLES], int plot)
width = BOUNDARY_WIDTH;
break;
}
case (P_DIFF_NEIGHB):
{
hue = (double)(particle[j].diff_neighb+1)/(double)(particle[j].neighb+1);
// hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*hue;
hue = 180.0*(1.0 + hue);
radius = particle[j].radius;
width = BOUNDARY_WIDTH;
break;
}
}
switch (particle[j].interaction) {
@ -2139,6 +2392,26 @@ void draw_particles(t_particle particle[NMAXCIRCLES], int plot)
hsl_to_rgb_twilight(hue, 0.9, 0.5, rgby);
break;
}
case (P_DIRECT_ENERGY):
{
hsl_to_rgb_twilight(hue, 0.9, 0.5, rgb);
hsl_to_rgb_twilight(hue, 0.9, 0.5, rgbx);
hsl_to_rgb_twilight(hue, 0.9, 0.5, rgby);
for (i=0; i<3; i++)
{
rgb[i] *= lum;
rgbx[i] *= lum;
rgby[i] *= lum;
}
break;
}
case (P_DIFF_NEIGHB):
{
hsl_to_rgb_twilight(hue, 0.9, 0.5, rgb);
hsl_to_rgb_twilight(hue, 0.9, 0.5, rgbx);
hsl_to_rgb_twilight(hue, 0.9, 0.5, rgby);
break;
}
default:
{
hsl_to_rgb(hue, 0.9, 0.5, rgb);
@ -2257,7 +2530,7 @@ void draw_particles(t_particle particle[NMAXCIRCLES], int plot)
}
void draw_container(double xmin, double xmax, t_obstacle obstacle[NMAXOBSTACLES], int wall)
void draw_container(double xmin, double xmax, t_obstacle obstacle[NMAXOBSTACLES], t_segment segment[NMAXSEGMENTS], int wall)
/* draw the container, for certain boundary conditions */
{
int i, j;
@ -2275,6 +2548,13 @@ void draw_container(double xmin, double xmax, t_obstacle obstacle[NMAXOBSTACLES]
for (i = 0; i < nobstacles; i++)
draw_circle(obstacle[i].xc, obstacle[i].yc, obstacle[i].radius, NSEG);
}
if (ADD_FIXED_SEGMENTS)
{
glLineWidth(CONTAINER_WIDTH);
glColor3f(1.0, 1.0, 1.0);
for (i = 0; i < nsegments; i++) if (segment[i].active)
draw_line(segment[i].x1, segment[i].y1, segment[i].x2, segment[i].y2);
}
switch (BOUNDARY_COND) {
case (BC_SCREEN):
@ -2462,20 +2742,8 @@ void draw_container(double xmin, double xmax, t_obstacle obstacle[NMAXOBSTACLES]
hsl_to_rgb(300.0, 0.1, 0.5, rgb);
draw_colored_rectangle(0.0, 0.0, BCXMAX, BCYMAX, rgb);
}
}
// /* draw fixed obstacles */
// if (ADD_FIXED_OBSTACLES)
// {
// glLineWidth(CONTAINER_WIDTH);
// hsl_to_rgb(300.0, 0.1, 0.5, rgb);
// for (i = 0; i < nobstacles; i++)
// draw_colored_circle(obstacle[i].xc, obstacle[i].yc, obstacle[i].radius, NSEG, rgb);
// glColor3f(1.0, 1.0, 1.0);
// for (i = 0; i < nobstacles; i++)
// draw_circle(obstacle[i].xc, obstacle[i].yc, obstacle[i].radius, NSEG);
// }
}
void print_parameters(double beta, double temperature, double krepel, double lengthcontainer, double boundary_force,
@ -2628,7 +2896,7 @@ void print_parameters(double beta, double temperature, double krepel, double len
}
}
if (PARTIAL_THERMO_COUPLING)
if ((PARTIAL_THERMO_COUPLING)&&(!INCREASE_BETA))
{
printf("Temperature %i in average: %.3lg\n", i_temp, temperature);
temp[i_temp] = temperature;
@ -2828,15 +3096,77 @@ void print_entropy(double entropy[2])
}
void print_omega(double angular_speed)
{
char message[100];
double rgb[3];
static double xleftbox, xlefttext, xrightbox, xrighttext, y = YMAX - 0.1, ymin = YMIN + 0.05;
static int first = 1;
if (first)
{
// xleftbox = XMIN + 0.4;
// xlefttext = xleftbox - 0.55;
xrightbox = XMAX - 0.39;
xrighttext = xrightbox - 0.55;
first = 0;
}
erase_area_hsl(xrightbox, y + 0.025, 0.35, 0.05, 0.0, 0.9, 0.0);
glColor3f(1.0, 1.0, 1.0);
sprintf(message, "Angular speed = %.4f", DPI*angular_speed*25.0/(double)(PERIOD_ROTATE_BOUNDARY));
write_text(xrighttext + 0.1, y, message);
}
void print_particles_speeds(t_particle particle[NMAXCIRCLES])
{
char message[100];
double y = YMAX - 0.1, vx = 0.0, vy = 0.0;
int i, nactive = 0;
static double xleftbox, xlefttext, xrightbox, xrighttext, ymin = YMIN + 0.05;
static int first = 1;
if (first)
{
xrightbox = XMAX - 0.39;
xrighttext = xrightbox - 0.55;
first = 0;
}
for (i=0; i<NMAXCIRCLES; i++) if (particle[i].active)
{
nactive++;
vx += particle[i].vx;
vy += particle[i].vy;
}
vx = vx/(double)nactive;
vy = vy/(double)nactive;
erase_area_hsl(xrightbox, y + 0.025, 0.35, 0.05, 0.0, 0.9, 0.0);
glColor3f(1.0, 1.0, 1.0);
sprintf(message, "Average vx = %.4f", vx);
write_text(xrighttext + 0.1, y, message);
y -= 0.1;
erase_area_hsl(xrightbox, y + 0.025, 0.35, 0.05, 0.0, 0.9, 0.0);
glColor3f(1.0, 1.0, 1.0);
sprintf(message, "Average vy = %.4f", vy);
write_text(xrighttext + 0.1, y, message);
}
double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NMAXOBSTACLES],
t_segment segment[NMAXSEGMENTS],
double xleft, double xright, double *pleft, double *pright, double pressure[N_PRESSURES], int wall)
{
int i, k;
double xmin, xmax, ymin, ymax, padding, r, rp, r2, cphi, sphi, f, fperp = 0.0, x, y, xtube, distance, dx, dy,
width, ybin, angle, x1, x2, h, ytop, norm, dleft, dplus, dminus, tmp_pleft = 0.0, tmp_pright = 0.0;
width, ybin, angle, x1, x2, h, ytop, norm, dleft, dplus, dminus, tmp_pleft = 0.0, tmp_pright = 0.0, proj;
/* compute force from fixed obstacles */
/* compute force from fixed circular obstacles */
if (ADD_FIXED_OBSTACLES) for (i=0; i<nobstacles; i++)
{
x = particle[j].xc - obstacle[i].xc;
@ -2854,15 +3184,67 @@ double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacl
particle[j].fy += f*sphi;
}
}
/* compute force from fixed linear obstacles */
if (ADD_FIXED_SEGMENTS) for (i=0; i<nsegments; i++) if (segment[i].active)
{
x = particle[j].xc;
y = particle[j].yc;
proj = (segment[i].ny*(x - segment[i].x1) - segment[i].nx*(y - segment[i].y1))/segment[i].length;
if ((proj > 0.0)&&(proj < 1.0))
{
distance = segment[i].nx*x + segment[i].ny*y - segment[i].c;
r = 1.5*particle[j].radius;
if (vabs(distance) < r)
{
f = KSPRING_OBSTACLE*(r - distance);
particle[j].fx += f*segment[i].nx;
particle[j].fy += f*segment[i].ny;
if (MOVE_BOUNDARY)
{
segment[i].fx -= f*segment[i].nx;
segment[i].fy -= f*segment[i].ny;
// segment[i].fx += f*segment[i].nx;
// segment[i].fy += f*segment[i].ny;
}
}
}
/* compute force from concave corners */
if (segment[i].concave)
{
distance = module2(x - segment[i].x1, y - segment[i].y1);
angle = argument(x - segment[i].x1, y - segment[i].y1);
if (angle < segment[i].angle1) angle += DPI;
r = 1.5*particle[j].radius;
if ((distance < r)&&(angle > segment[i].angle1)&&(angle < segment[i].angle2))
{
f = KSPRING_OBSTACLE*(r - distance);
particle[j].fx += f*cos(angle);
particle[j].fy += f*sin(angle);
if (MOVE_BOUNDARY)
{
segment[i].fx -= f*cos(angle);
segment[i].fy -= f*sin(angle);
// segment[i].fx += f*cos(angle);
// segment[i].fy += f*sin(angle);
}
}
}
}
switch(BOUNDARY_COND){
case (BC_SCREEN):
{
/* add harmonic force outside screen */
if (particle[j].xc > BCXMAX) particle[j].fx -= KSPRING_BOUNDARY*(particle[j].xc - BCXMAX);
else if (particle[j].xc < BCXMIN) particle[j].fx += KSPRING_BOUNDARY*(BCXMIN - particle[j].xc);
if (particle[j].yc > BCYMAX) particle[j].fy -= KSPRING_BOUNDARY*(particle[j].yc - BCYMAX);
else if (particle[j].yc < BCYMIN) particle[j].fy += KSPRING_BOUNDARY*(BCYMIN - particle[j].yc);
if (particle[j].xc > XMAX) particle[j].fx -= KSPRING_BOUNDARY*(particle[j].xc - XMAX);
else if (particle[j].xc < XMIN) particle[j].fx += KSPRING_BOUNDARY*(XMIN - particle[j].xc);
if (particle[j].yc > YMAX) particle[j].fy -= KSPRING_BOUNDARY*(particle[j].yc - YMAX);
else if (particle[j].yc < YMIN) particle[j].fy += KSPRING_BOUNDARY*(YMIN - particle[j].yc);
// if (particle[j].xc > BCXMAX) particle[j].fx -= KSPRING_BOUNDARY*(particle[j].xc - BCXMAX);
// else if (particle[j].xc < BCXMIN) particle[j].fx += KSPRING_BOUNDARY*(BCXMIN - particle[j].xc);
// if (particle[j].yc > BCYMAX) particle[j].fy -= KSPRING_BOUNDARY*(particle[j].yc - BCYMAX);
// else if (particle[j].yc < BCYMIN) particle[j].fy += KSPRING_BOUNDARY*(BCYMIN - particle[j].yc);
return(fperp);
}
case (BC_RECTANGLE):
@ -3235,6 +3617,7 @@ void compute_particle_force(int j, double krepel, t_particle particle[NMAXCIRCLE
double fx = 0.0, fy = 0.0, force[2], torque = 0.0, torque_ij, x, y;
particle[j].neighb = 0;
if (REACTION_DIFFUSION) particle[j].diff_neighb = 0;
for (k=0; k<particle[j].hash_nneighb; k++)
{
@ -3242,7 +3625,12 @@ void compute_particle_force(int j, double krepel, t_particle particle[NMAXCIRCLE
fx += force[0];
fy += force[1];
torque += torque_ij;
if (close) particle[j].neighb++;
if (close)
{
particle[j].neighb++;
if (REACTION_DIFFUSION&&(particle[j].type != particle[particle[j].hashneighbour[k]].type))
particle[j].diff_neighb++;
}
}
particle[j].fx += fx;
@ -3264,11 +3652,12 @@ int initialize_configuration(t_particle particle[NMAXCIRCLES], t_hashgrid hashgr
particle[i].type = 0;
if ((TWO_TYPES)&&((double)rand()/RAND_MAX > TPYE_PROPORTION))
{
particle[i].type = 1;
particle[i].type = 2;
particle[i].radius = MU_B;
}
particle[i].neighb = 0;
particle[i].diff_neighb = 0;
particle[i].thermostat = 1;
// particle[i].energy = 0.0;
@ -3462,17 +3851,29 @@ int initialize_configuration(t_particle particle[NMAXCIRCLES], t_hashgrid hashgr
}
if (ADD_FIXED_OBSTACLES)
{
for (i=0; i< ncircles; i++) for (j=0; j < nobstacles; j++)
if (module2(particle[i].xc - obstacle[j].xc, particle[i].yc - obstacle[j].yc) < OBSTACLE_RADIUS + particle[i].radius)
particle[i].active = 0;
}
/* case of segment obstacles */
if (ADD_FIXED_SEGMENTS) for (i=0; i< ncircles; i++)
if (!in_segment_region(particle[i].xc, particle[i].yc))
particle[i].active = 0;
/* case of reaction-diffusion equation */
if (REACTION_DIFFUSION) for (i=0; i< ncircles; i++)
{
particle[i].type = 1 + (int)(RD_TYPES*(double)rand()/(double)RAND_MAX);
}
/* count number of active particles */
for (i=0; i< ncircles; i++) nactive += particle[i].active;
printf("%i active particles\n", nactive);
printf("%i active particles\n", nactive);
for (i=0; i<ncircles; i++) printf("particle %i of type %i\n", i, particle[i].type);
return(nactive);
}
@ -3481,6 +3882,7 @@ int initialize_configuration(t_particle particle[NMAXCIRCLES], t_hashgrid hashgr
int add_particles(t_particle particle[NMAXCIRCLES], double px[NMAXCIRCLES], double py[NMAXCIRCLES], int nadd_particle)
/* add several particles to the system */
{
static int i = 0;
// add_particle(XMIN + 0.1, 0.0, 50.0, 0.0, 3.0, 0, particle);
// px[ncircles - 1] = particle[ncircles - 1].vx;
// py[ncircles - 1] = particle[ncircles - 1].vy;
@ -3504,11 +3906,14 @@ int add_particles(t_particle particle[NMAXCIRCLES], double px[NMAXCIRCLES], doub
// particle[ncircles - 1].thermostat = 0;
// px[ncircles - 1] = particle[ncircles - 1].vx;
// py[ncircles - 1] = particle[ncircles - 1].vy;
// add_particle(MU*(2.0*rand()/RAND_MAX - 1.0), YMAX + 2.0*MU, 0.0, 0.0, PARTICLE_MASS, 0, particle);
printf("Adding a particle\n\n");
add_particle(MU*(2.0*rand()/RAND_MAX - 1.0), YMAX + 2.0*MU, 0.0, 0.0, PARTICLE_MASS, 0, particle);
particle[ncircles - 1].radius = MU;
add_particle(XMIN - 0.5*MU, 0.0, 50.0 + 5.0*(double)i, 0.0, 2.0*PARTICLE_MASS, 0, particle);
i++;
particle[ncircles - 1].radius = 0.5*MU;
particle[ncircles - 1].eq_dist = EQUILIBRIUM_DIST;
particle[ncircles - 1].thermostat = 0;
px[ncircles - 1] = particle[ncircles - 1].vx;
@ -3553,11 +3958,24 @@ int floor_momentum(double p[NMAXCIRCLES])
int partial_thermostat_coupling(t_particle particle[NMAXCIRCLES], double xmin)
/* only couple particles with x > xmin to thermostat */
{
int i, nthermo = 0;
int condition, i, nthermo = 0;
for (i=0; i<ncircles; i++)
{
if (particle[i].xc > xmin)
switch (PARTIAL_THERMO_REGION) {
case (TH_VERTICAL):
{
condition = (particle[i].xc > xmin);
break;
}
case (TH_INSEGMENT):
{
condition = (in_segment_region(particle[i].xc - xsegments, particle[i].yc - ysegments));
break;
}
default: condition = 1;
}
if (condition)
{
particle[i].thermostat = 1;
nthermo++;
@ -3578,4 +3996,24 @@ double compute_mean_energy(t_particle particle[NMAXCIRCLES])
nactive++;
}
return(total_energy/(double)nactive);
}
}
void update_types(t_particle particle[NMAXCIRCLES])
/* update the types in case of reaction-diffusion equation */
{
int i, j, k;
double distance;
for (i=0; i<ncircles; i++)
for (j=0; j<particle[i].hash_nneighb; j++)
{
k = particle[i].hashneighbour[j];
if ((particle[k].type == particle[i].type + 1)||((particle[i].type == RD_TYPES)&&(particle[k].type == 1)))
{
distance = module2(particle[i].deltax[j], particle[i].deltay[j]);
if ((distance < EQUILIBRIUM_DIST)&&((double)rand()/RAND_MAX < REACTION_PROB))
particle[k].type = particle[i].type;
}
}
}

1123
sub_rde.c Normal file
View File

File diff suppressed because it is too large Load Diff

13
sub_wave.c
View File

@ -1434,10 +1434,14 @@ int compute_noisepanel_coordinates(t_vertex polyline[NMAXPOLY])
if (even) y = YMIN + 0.1;
else y = YMIN + 0.1 + MU;
polyline[i].x = x;
x1 = x;
if (x1 > XMAX) x1 = XMAX;
else if (x1 < XMIN) x1 = XMIN;
polyline[i].x = x1;
polyline[i].y = y;
xy_to_pos(x, y, pos);
xy_to_pos(x1, y, pos);
polyline[i].posi = pos[0];
polyline[i].posj = pos[1];
@ -1612,6 +1616,11 @@ int xy_in_billiard(double x, double y)
static int first = 1, nsides;
switch (B_DOMAIN) {
case (D_NOTHING):
{
return(1);
break;
}
case (D_RECTANGLE):
{
if ((vabs(x) <LAMBDA)&&(vabs(y) < 1.0)) return(1);

811
sub_wave_3d.c
View File

File diff suppressed because it is too large Load Diff

1430
sub_wave_3d_rde.c Normal file
View File

File diff suppressed because it is too large Load Diff

257
wave_3d.c
View File

@ -41,34 +41,38 @@
#include <sys/types.h>
#include <tiffio.h> /* Sam Leffler's libtiff library. */
#include <omp.h>
#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 */
/* General geometrical parameters */
/* uncomment for higher resolution */
// #define WINWIDTH 1920 /* window width */
// #define WINHEIGHT 1000 /* window height */
// #define NX 1920 /* number of grid points on x axis */
// #define NY 1000 /* number of grid points on y axis */
// // #define NX 3840 /* number of grid points on x axis */
// // #define NY 2000 /* number of grid points on y axis */
// // #define NX 900 /* number of grid points on x axis */
// // #define NY 450 /* number of grid points on y axis */
// #define NX 1800 /* number of grid points on x axis */
// #define NY 900 /* number of grid points on y axis */
// // #define NX 1920 /* number of grid points on x axis */
// // #define NY 1000 /* number of grid points on y axis */
//
// #define XMIN -2.0
// #define XMAX 2.0 /* x interval */
// // #define 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 0 /* set to 1 if resolution of grid is double that of displayed image */
/* comment out for higher resolution */
#define WINWIDTH 1280 /* window width */
#define WINHEIGHT 720 /* window height */
#define NX 1280 /* number of grid points on x axis */
#define NY 720 /* number of grid points on y axis */
// #define NX 640 /* number of grid points on x axis */
// #define NY 360 /* number of grid points on y axis */
#define XMIN -2.0
#define XMAX 2.0 /* x interval */
#define YMIN -1.125
@ -78,25 +82,29 @@
/* Choice of the billiard table */
#define B_DOMAIN 16 /* choice of domain shape, see list in global_pdes.c */
#define B_DOMAIN 44 /* choice of domain shape, see list in global_pdes.c */
#define CIRCLE_PATTERN 201 /* pattern of circles or polygons, see list in global_pdes.c */
#define CIRCLE_PATTERN 202 /* pattern of circles or polygons, see list in global_pdes.c */
#define VARIABLE_IOR 0 /* set to 1 for a variable index of refraction */
#define IOR 1 /* choice of index of refraction, see list in global_pdes.c */
#define MANDEL_IOR_SCALE -0.05 /* parameter controlling dependece of IoR on Mandelbrot escape speed */
#define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */
#define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */
#define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */
#define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */
#define LAMBDA 0.6 /* parameter controlling the dimensions of domain */
#define MU 0.6 /* 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 LAMBDA 0.2 /* parameter controlling the dimensions of domain */
#define MU 0.5 /* 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 */
#define MDEPTH 3 /* depth of computation of Menger gasket */
#define MRATIO 3 /* ratio defining Menger gasket */
#define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */
#define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */
#define FOCI 1 /* set to 1 to draw focal points of ellipse */
#define NGRIDX 36 /* number of grid point for grid of disks */
#define NGRIDY 6 /* number of grid point for grid of disks */
#define NGRIDY 6 /* number of grid point for grid of disks */
#define X_SHOOTER -0.2
#define Y_SHOOTER -0.6
@ -115,18 +123,16 @@
/* Physical parameters of wave equation */
// #define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */
#define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */
#define TWOSPEEDS 1 /* set to 1 to replace hardcore boundary by medium with different speed */
#define OSCILLATE_LEFT 0 /* set to 1 to add oscilating boundary condition on the left */
#define OSCILLATE_TOPBOT 0 /* set to 1 to enforce a planar wave on top and bottom boundary */
#define OMEGA 0.005 /* frequency of periodic excitation */
#define AMPLITUDE 0.8 /* amplitude of periodic excitation */
#define COURANT 0.06 /* Courant number */
#define COURANTB 0.03 /* Courant number in medium B */
// #define COURANTB 0.016363636 /* Courant number in medium B */
#define OMEGA 0.017 /* frequency of periodic excitation */
#define AMPLITUDE 0.9 /* amplitude of periodic excitation */
#define COURANT 0.1 /* Courant number */
#define COURANTB 0.05 /* Courant number in medium B */
#define GAMMA 0.0 /* damping factor in wave equation */
#define GAMMAB 1.0e-7 /* damping factor in wave equation */
#define GAMMAB 2.0e-5 /* 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 */
@ -139,53 +145,54 @@
#define ADD_OSCILLATING_SOURCE 0 /* set to 1 to add an oscillating wave source */
#define OSCILLATING_SOURCE_PERIOD 30 /* period of oscillating source */
// #define OSCILLATING_SOURCE_PERIOD 14 /* period of oscillating source */
/* Boundary conditions, see list in global_pdes.c */
#define B_COND 2
// #define B_COND 2
/* Parameters for length and speed of simulation */
#define NSTEPS 2500 /* number of frames of movie */
#define NVID 10 /* number of iterations between images displayed on screen */
#define NSTEPS 1500 /* number of frames of movie */
#define NVID 8 /* 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 3 /* width of billiard boundary */
#define BOUNDARY_WIDTH 1 /* width of billiard boundary */
#define PAUSE 200 /* number of frames after which to pause */
#define PAUSE 100 /* number of frames after which to pause */
#define PSLEEP 2 /* sleep time during pause */
#define SLEEP1 1 /* initial sleeping time */
#define SLEEP2 1 /* final sleeping time */
#define MID_FRAMES 200 /* number of still frames between parts of two-part movie */
#define MID_FRAMES 100 /* number of still frames between parts of two-part movie */
#define END_FRAMES 100 /* number of still frames at end of movie */
#define FADE 1 /* set to 1 to fade at end of movie */
/* Parameters of initial condition */
#define INITIAL_AMP 0.5 /* amplitude of initial condition */
#define INITIAL_VARIANCE 0.0005 /* variance of initial condition */
#define INITIAL_WAVELENGTH 0.1 /* wavelength of initial condition */
#define INITIAL_AMP 0.75 /* amplitude of initial condition */
#define INITIAL_VARIANCE 0.0003 /* variance of initial condition */
#define INITIAL_WAVELENGTH 0.02 /* wavelength of initial condition */
/* Plot type, see list in global_pdes.c */
#define ZPLOT 103 /* wave height */
#define CPLOT 103 /* color scheme */
#define ZPLOT 108 /* wave height */
#define CPLOT 108 /* color scheme */
#define ZPLOT_B 104
#define CPLOT_B 104 /* plot type for second movie */
#define ZPLOT_B 109
#define CPLOT_B 109 /* 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 */
#define NON_DIRICHLET_BC 0 /* set to 1 to draw only facets in domain, if field is not zero on boundary */
#define FLOOR_ZCOORD 1 /* set to 1 to draw only facets with z not too negative */
#define DRAW_BILLIARD 1 /* set to 1 to draw boundary */
#define DRAW_BILLIARD_FRONT 1 /* set to 1 to draw front of boundary after drawing wave */
#define DRAW_BILLIARD_FRONT 0 /* set to 1 to draw front of boundary after drawing wave */
#define DRAW_CONSTRUCTION_LINES 0 /* set to 1 to draw construction lines of certain domains */
#define FADE_IN_OBSTACLE 1 /* set to 1 to fade color inside obstacles */
#define DRAW_OUTSIDE_GRAY 0 /* experimental, draw outside of billiard in gray */
#define PLOT_SCALE_ENERGY 0.05 /* vertical scaling in energy plot */
#define PLOT_SCALE_LOG_ENERGY 0.6 /* vertical scaling in log energy plot */
#define PLOT_SCALE_ENERGY 0.02 /* vertical scaling in energy plot */
#define PLOT_SCALE_LOG_ENERGY 1.0 /* vertical scaling in log energy plot */
/* 3D representation */
@ -194,10 +201,9 @@
#define REP_AXO_3D 0 /* linear projection (axonometry) */
#define REP_PROJ_3D 1 /* projection on plane orthogonal to observer line of sight */
/* Color schemes */
#define COLOR_PALETTE 14 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE 14 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 11 /* Color palette, see list in global_pdes.c */
#define BLACK 1 /* background */
@ -205,15 +211,17 @@
#define COLOR_SCHEME 3 /* choice of color scheme, see list in global_pdes.c */
#define SCALE 0 /* set to 1 to adjust color scheme to variance of field */
#define SLOPE 1.0 /* sensitivity of color on wave amplitude */
#define VSCALE_AMPLITUDE 0.2 /* additional scaling factor for color scheme P_3D_AMPLITUDE */
#define VSCALE_ENERGY 0.35 /* additional scaling factor for color scheme P_3D_ENERGY */
#define SLOPE 1.0 /* sensitivity of color on wave amplitude */
#define VSCALE_AMPLITUDE 1.0 /* additional scaling factor for color scheme P_3D_AMPLITUDE */
#define VSCALE_ENERGY 10.0 /* additional scaling factor for color scheme P_3D_ENERGY */
#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 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 1000.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 LOG_ENERGY_FLOOR -1000.0 /* floor value for log of (total) energy */
#define LOG_MEAN_ENERGY_SHIFT 5.0 /* additional shift for log of mean energy */
#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 */
@ -223,9 +231,9 @@
#define HUEMEAN 240.0 /* mean value of hue for color scheme C_HUE */
#define HUEAMP -200.0 /* amplitude of variation of hue for color scheme C_HUE */
#define DRAW_COLOR_SCHEME 0 /* set to 1 to plot the color scheme */
#define COLORBAR_RANGE 3.0 /* scale of color scheme bar */
#define COLORBAR_RANGE_B 5.0 /* scale of color scheme bar for 2nd part */
#define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */
#define COLORBAR_RANGE 20.0 /* scale of color scheme bar */
#define COLORBAR_RANGE_B 100.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 */
@ -240,13 +248,13 @@ double u_3d[2] = {0.75, -0.45}; /* projections of basis vectors for REP_AXO_
double v_3d[2] = {-0.75, -0.45};
double w_3d[2] = {0.0, 0.015};
double light[3] = {0.816496581, -0.40824829, 0.40824829}; /* vector of "light" direction for P_3D_ANGLE color scheme */
double observer[3] = {10.0, 6.0, 8.5}; /* location of observer for REP_PROJ_3D representation */
double observer[3] = {6.0, 8.0, 6.0}; /* location of observer for REP_PROJ_3D representation */
#define Z_SCALING_FACTOR 0.018 /* overall scaling factor of z axis for REP_PROJ_3D representation */
#define XY_SCALING_FACTOR 3.75 /* overall scaling factor for on-screen (x,y) coordinates after projection */
#define Z_SCALING_FACTOR 0.05 /* overall scaling factor of z axis for REP_PROJ_3D representation */
#define XY_SCALING_FACTOR 2.0 /* 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.0 /* overall x shift for REP_PROJ_3D representation */
#define YSHIFT_3D 0.0 /* overall y shift for REP_PROJ_3D representation */
#define XSHIFT_3D -0.1 /* overall x shift for REP_PROJ_3D representation */
#define YSHIFT_3D 0.15 /* overall y shift for REP_PROJ_3D representation */
#include "global_pdes.c" /* constants and global variables */
@ -500,10 +508,15 @@ void evolve_wave(double phi[NX*NY], double psi[NX*NY], double phi_tmp[NX*NY], do
}
void draw_color_bar_palette(int plot, double range, int palette)
void draw_color_bar_palette(int plot, double range, int palette, int fade, double fade_value)
{
if (ROTATE_COLOR_SCHEME) draw_color_scheme_palette_3d(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range, palette);
else draw_color_scheme_palette_3d(XMAX - 0.3, YMIN + 0.1, XMAX - 0.1, YMAX - 0.1, plot, -range, range, palette);
double width = 0.14;
// double width = 0.2;
if (ROTATE_COLOR_SCHEME)
draw_color_scheme_palette_3d(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range, palette, fade, fade_value);
else
draw_color_scheme_palette_3d(XMAX - 1.5*width, YMIN + 0.1, XMAX - 0.5*width, YMAX - 0.1, plot, -range, range, palette, fade, fade_value);
}
@ -550,9 +563,6 @@ void animation()
courant2 = COURANT*COURANT;
courantb2 = COURANTB*COURANTB;
/* initialize color scale, for option RESCALE_COLOR_IN_CENTER */
if (RESCALE_COLOR_IN_CENTER)
{
@ -576,84 +586,33 @@ void animation()
ratio = (XMAX - XMIN)/8.4; /* for Tokarsky billiard */
// isospectral_initial_point(0.2, 0.0, startleft, startright); /* for isospectral billiards */
// homophonic_initial_point(0.5, -0.25, 1.5, -0.25, startleft, startright);
// homophonic_initial_point(0.5, -0.25, 1.5, -0.25, startleft, startright);
// printf("xleft = (%.3f, %.3f) xright = (%.3f, %.3f)\n", startleft[0], startleft[1], startright[0], startright[1]);
// xy_to_ij(startleft[0], startleft[1], sample_left);
// xy_to_ij(startright[0], startright[1], sample_right);
// printf("xleft = (%.3f, %.3f) xright = (%.3f, %.3f)\n", xin_left, yin_left, xin_right, yin_right);
// init_wave_flat(phi, psi, xy_in);
// init_wave_plus(LAMBDA - 0.3*MU, 0.5*MU, phi, psi, xy_in);
// init_wave(LAMBDA - 0.3*MU, 0.5*MU, phi, psi, xy_in);
// init_circular_wave(X_SHOOTER, Y_SHOOTER, phi, psi, xy_in);
// printf("Initializing wave\n");
// init_circular_wave_mod(polyline[85].x, polyline[85].y, phi, psi, xy_in);
// init_circular_wave_mod(0.0, 0.0, phi, psi, xy_in);
init_circular_wave_mod(0.2, 0.4, phi, psi, xy_in);
add_circular_wave_mod(-1.0, -0.2, -0.4, phi, psi, xy_in);
init_circular_wave_mod(0.0, 0.0, phi, psi, xy_in);
// add_circular_wave(-1.0, -0.2, -0.4, phi, psi, xy_in);
// printf("Wave initialized\n");
// init_circular_wave(0.6*cos((double)(period)*DPI/3.0), 0.6*sin((double)(period)*DPI/3.0), phi, psi, xy_in);
// period++;
// for (i=0; i<3; i++)
// {
// add_circular_wave(-1.0, 0.6*cos(PID + (double)(i)*DPI/3.0), 0.6*sin(PID + (double)(i)*DPI/3.0), phi, psi, xy_in);
// }
// add_circular_wave(1.0, -LAMBDA, 0.0, phi, psi, xy_in);
// add_circular_wave(-1.0, 0.0, -LAMBDA, phi, psi, xy_in);
// init_circular_wave_xplusminus(startleft[0], startleft[1], startright[0], startright[1], phi, psi, xy_in);
// init_circular_wave_xplusminus(-0.9, 0.0, 0.81, 0.0, phi, psi, xy_in);
// init_circular_wave(-2.0*ratio, 0.0, phi, psi, xy_in);
// init_planar_wave(XMIN + 0.015, 0.0, phi, psi, xy_in);
// init_planar_wave(XMIN + 0.02, 0.0, phi, psi, xy_in);
// init_planar_wave(XMIN + 0.5, 0.0, phi, psi, xy_in);
// init_wave(-1.5, 0.0, phi, psi, xy_in);
// init_wave(0.0, 0.0, phi, psi, xy_in);
/* add a drop at another point */
// add_drop_to_wave(1.0, 0.7, 0.0, phi, psi);
// add_drop_to_wave(1.0, -0.7, 0.0, phi, psi);
// add_drop_to_wave(1.0, 0.0, -0.7, phi, psi);
/* initialize table of wave speeds/dissipation */
for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
if (xy_in[i*NY+j] != 0)
{
tc[i*NY+j] = COURANT;
tcc[i*NY+j] = courant2;
if (xy_in[i*NY+j] == 1) tgamma[i*NY+j] = GAMMA;
else tgamma[i*NY+j] = GAMMAB;
}
else if (TWOSPEEDS)
{
tc[i*NY+j] = COURANTB;
tcc[i*NY+j] = courantb2;
tgamma[i*NY+j] = GAMMAB;
}
}
init_speed_dissipation(xy_in, tc, tcc, tgamma);
init_zfield(phi, psi, xy_in, ZPLOT, wave, 0);
init_cfield(phi, psi, xy_in, CPLOT, wave, 0);
if (DOUBLE_MOVIE)
{
init_zfield(phi, psi, xy_in, ZPLOT_B, wave, 1);
init_cfield(phi, psi, xy_in, CPLOT_B, wave, 1);
}
blank();
glColor3f(0.0, 0.0, 0.0);
draw_wave_3d(phi, psi, xy_in, wave, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0);
draw_wave_3d(0, phi, psi, xy_in, wave, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1);
// draw_billiard();
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 0, 1.0);
glutSwapBuffers();
@ -663,6 +622,8 @@ void animation()
for (i=0; i<=INITIAL_TIME + NSTEPS; i++)
{
global_time++;
//printf("%d\n",i);
/* compute the variance of the field to adjust color scheme */
/* the color depends on the field divided by sqrt(1 + variance) */
@ -673,7 +634,7 @@ void animation()
}
else scale = 1.0;
draw_wave_3d(phi, psi, xy_in, wave, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0);
draw_wave_3d(0, phi, psi, xy_in, wave, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1);
for (j=0; j<NVID; j++)
{
evolve_wave(phi, psi, phi_tmp, psi_tmp, xy_in, tc, tcc, tgamma);
@ -691,7 +652,7 @@ void animation()
// draw_billiard();
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, fade, fade_value);
/* add oscillating waves */
if ((ADD_OSCILLATING_SOURCE)&&(i%OSCILLATING_SOURCE_PERIOD == OSCILLATING_SOURCE_PERIOD - 1))
@ -711,11 +672,10 @@ void animation()
if ((i >= INITIAL_TIME)&&(DOUBLE_MOVIE))
{
draw_wave_3d(phi, psi, xy_in, wave, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 0, 1.0);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B);
draw_wave_3d(1, phi, psi, xy_in, wave, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 0, 1.0, 1);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0);
glutSwapBuffers();
save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter);
// save_frame_counter(NSTEPS + 21 + counter);
counter++;
}
@ -735,27 +695,29 @@ void animation()
{
if (DOUBLE_MOVIE)
{
draw_wave_3d(phi, psi, xy_in, wave, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE);
draw_wave_3d(0, phi, psi, xy_in, wave, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 0, 1.0);
glutSwapBuffers();
if (!FADE) for (i=0; i<MID_FRAMES; i++) save_frame();
else for (i=0; i<MID_FRAMES; i++)
{
draw_wave_3d(phi, psi, xy_in, wave, ZPLOT, CPLOT, COLOR_PALETTE, 1, 1.0 - (double)i/(double)MID_FRAMES);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE);
fade_value = 1.0 - (double)i/(double)MID_FRAMES;
draw_wave_3d(0, phi, psi, xy_in, wave, ZPLOT, CPLOT, COLOR_PALETTE, 1, fade_value, 0);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 1, fade_value);
glutSwapBuffers();
save_frame_counter(NSTEPS + i + 1);
}
draw_wave_3d(phi, psi, xy_in, wave, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 0, 1.0);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B);
draw_wave_3d(1, phi, psi, xy_in, wave, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 0, 1.0, 1);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0);
glutSwapBuffers();
if (!FADE) for (i=0; i<END_FRAMES; i++) save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter + i);
else for (i=0; i<END_FRAMES; i++)
{
draw_wave_3d(phi, psi, xy_in, wave, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 1, 1.0 - (double)i/(double)END_FRAMES);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B);
fade_value = 1.0 - (double)i/(double)END_FRAMES;
draw_wave_3d(1, phi, psi, xy_in, wave, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 1, fade_value, 0);
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);
}
@ -765,8 +727,9 @@ void animation()
if (!FADE) for (i=0; i<END_FRAMES; i++) save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter + i);
else for (i=0; i<END_FRAMES; i++)
{
draw_wave_3d(phi, psi, xy_in, wave, ZPLOT, CPLOT, COLOR_PALETTE, 1, 1.0 - (double)i/(double)END_FRAMES);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE);
fade_value = 1.0 - (double)i/(double)END_FRAMES;
draw_wave_3d(0, phi, psi, xy_in, wave, ZPLOT, CPLOT, COLOR_PALETTE, 1, fade_value, 0);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 1, fade_value);
glutSwapBuffers();
save_frame_counter(NSTEPS + 1 + counter + i);
}
@ -801,6 +764,12 @@ void animation()
void display(void)
{
time_t rawtime;
struct tm * timeinfo;
time(&rawtime);
timeinfo = localtime(&rawtime);
glPushMatrix();
blank();
@ -814,7 +783,11 @@ void display(void)
glPopMatrix();
glutDestroyWindow(glutGetWindow());
printf("Start local time and date: %s", asctime(timeinfo));
time(&rawtime);
timeinfo = localtime(&rawtime);
printf("Current local time and date: %s", asctime(timeinfo));
}

22
wave_common.c
View File

@ -856,7 +856,7 @@ void init_circular_wave_mod(double x, double y, double phi[NX*NY], double psi[NX
double xy[2], dist2;
printf("Initializing wave\n");
#pragma omp parallel for private(i,j,xy,dist2)
// #pragma omp parallel for private(i,j,xy,dist2)
for (i=0; i<NX; i++)
{
if (i%100 == 0) printf("Initializing column %i of %i\n", i, NX);
@ -892,6 +892,26 @@ void add_circular_wave_mod(double factor, double x, double y, double phi[NX*NY],
}
}
void init_wave_flat_mod(double phi[NX*NY], double psi[NX*NY], short int xy_in[NX*NY])
/* initialise flat field - phi is wave height, psi is phi at time t-1 */
{
int i, j;
double xy[2];
#pragma omp parallel for private(i,j,xy)
for (i=0; i<NX; i++) {
if (i%100 == 0) printf("Wave and table xy_in - Initialising column %i of %i\n", i, NX);
for (j=0; j<NY; j++)
{
ij_to_xy(i, j, xy);
xy_in[i*NY+j] = xy_in_billiard(xy[0],xy[1]);
phi[i*NY+j] = 0.0;
psi[i*NY+j] = 0.0;
}
}
}
double compute_variance_mod(double phi[NX*NY], double psi[NX*NY], short int xy_in[NX*NY])
/* compute the variance of the field, to adjust color scheme */
{