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

Add files via upload

Browse Source
This commit is contained in:
Nils Berglund 2023-07-08 23:28:23 +02:00 committed by GitHub
parent b809ce9e55
commit de54999d98
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
19 changed files with 3393 additions and 678 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
global_3d.c
View File

@ -27,6 +27,7 @@
#define E_SCHRODINGER 5 /* Schrodinger equation */
#define E_EULER_INCOMP 6 /* incompressible Euler equation */
#define E_EULER_COMP 7 /* compressible Euler equation */
#define E_SHALLOW_WATER 8 /* shallow water equation */
/* Choice of potential */
@ -51,7 +52,14 @@
#define GF_ELLIPSE 2 /* repelling ellipse */
#define GF_AIRFOIL 3 /* curved repelling ellipse */
#define GF_WING 4 /* wing shape */
#define GF_COMPUTE_FROM_BC 5 /* compute force field as gradient of bc_field */
#define GF_COMPUTE_FROM_BC 5 /* compute force field as gradient of bc_field2 */
/* Choice of water depth for shallow water equation */
#define SH_CIRCLE 1 /* circular shallower obstacle */
#define SH_CIRCLES 2 /* shallow obstacle specified by CIRCLE_PATTERN */
#define SH_COAST 3 /* depth varying with x-coordinate */
#define SH_COAST_MONOTONE 4 /* depth decreasing with x-coordinate */
/* macros to avoid unnecessary computations in 3D plots */
@ -74,6 +82,7 @@
int global_time = 0;
double max_depth = 1.0;
/* structure used for color and height representations */
/* possible extra fields: zfield, cfield, interpolated coordinates */
@ -111,11 +120,20 @@ typedef struct
double cos_angle; /* cos of angle between normal vector and direction of light */
double log_vorticity; /* logarithm of vorticity (for Euler equation) */
double Lpressure; /* Laplacian of pressure (for Euler equation) */
double height; /* wave height (for shallow wave equation) */
double dxu, dyu, dxv, dyv; /* gradient of velocity field (for compressible Euler equation) */
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) */
double depth; /* water depth */
double cos_depth_angle; /* cos of angle of depth profile */
double gradx, grady; /* gradient of water depth */
} t_rde;
typedef struct
{
double depth; /* water depth */
double gradx, grady; /* gradient of water depth */
} t_swater_depth;

11
global_ljones.c
View File

@ -44,6 +44,7 @@
#define O_GENUS_TWO 2 /* obstacles in corners of L-shape domeain (for genus 2 b.c.) */
#define O_POOL_TABLE 3 /* obstacles around pockets of pool table */
#define O_HLINE_HOLE_SPOKES 181 /* tips of spokes for S_HLINE_HOLE_SPOKES segment pattern */
#define O_CIRCLE 4 /* one circle at the origin */
/* pattern of additional repelling segments */
#define S_RECTANGLE 0 /* segments forming a rectangle */
@ -64,6 +65,7 @@
#define S_DAM_WITH_HOLE 13 /* segments forming a dam in which a hole can open */
#define S_DAM_WITH_HOLE_AND_RAMP 14 /* segments forming a dam in which a hole can open */
#define S_MAZE 15 /* segments forming a maze */
#define S_MAZE_DIAG 151 /* segments forming a maze with diagonally opposed exits */
#define S_EXT_RECTANGLE 16 /* particles outside a rectangle */
#define S_DAM_BRICKS 17 /* dam made of several bricks */
#define S_HLINE_HOLE 18 /* horizontal line with a hole in the bottom */
@ -84,6 +86,7 @@
#define I_VICSEK 8 /* Vicsek-type interaction */
#define I_VICSEK_REPULSIVE 9 /* Vicsek-type interaction with harmonic repulsion */
#define I_VICSEK_SPEED 10 /* Vicsek-type interaction with speed adjustment */
#define I_VICSEK_SHARK 11 /* Vicsek-type interaction with speed adjustment, and one shark */
/* Boundary conditions */
@ -103,6 +106,7 @@
#define BC_GENUS_TWO 14 /* surface of genus 2, obtained by identifying opposite sides of an L shape */
#define BC_ABSORBING 20 /* "no-return" boundary conditions outside BC area */
#define BC_REFLECT_ABS 21 /* reflecting on lower boundary, and "no-return" boundary conditions outside BC area */
#define BC_REFLECT_ABS_BOTTOM 22 /* absorbing on lower boundary, and reflecting elsewhere */
/* Regions for partial thermostat couplings */
@ -122,6 +126,7 @@
#define RCK_DISC 0 /* disc-shaped rocket */
#define RCK_RECT 1 /* rectangular rocket */
#define RCK_RECT_HAT 2 /* rectangular rocket with a hat */
#define RCK_RECT_BAR 3 /* rectangular rocket with a hat and a separating bar */
/* Nozzle shapes */
@ -131,6 +136,7 @@
#define NZ_CONE 3 /* cone-shaped nozzle */
#define NZ_TRUMPET 4 /* trumpet-shaped nozzle */
#define NZ_BROAD 5 /* broad straight nozzle */
#define NZ_DELAVAL 6 /* a type of de Laval nozzle */
#define NZ_NONE 99 /* no nozzle */
/* Types of chemical reactions */
@ -151,6 +157,7 @@
#define CHEM_ABCDABC 13 /* reactions A + B <-> C, A + C <-> D */
#define CHEM_BZ 14 /* simplified Belousov-Zhabotinski reaction with 6 types (Oregonator) */
#define CHEM_BRUSSELATOR 15 /* Brusselator oscillating reaction */
#define CHEM_ABDACBE 16 /* A + B -> D, A + C -> B + E */
/* Initial conditions for chemical reactions */
@ -164,6 +171,7 @@
#define IC_BZ 6 /* initial state for BZ reaction */
#define IC_SIGNX 7 /* type 1 or 2 depending on sign of x */
#define IC_TWOROCKETS 8 /* type 1 or 2 depending on rocket position */
#define IC_TWOROCKETS_TWOFUELS 9 /* type 1 and 2 or 1 and 3 depending on rocket */
/* Plot types */
@ -179,6 +187,7 @@
#define P_DIFF_NEIGHB 9 /* colors represent number of neighbours of different type */
#define P_THERMOSTAT 10 /* colors show which particles are coupled to the thermostat */
#define P_INITIAL_POS 11 /* colors depend on initial position of particle */
#define P_NUMBER 12 /* colors depend on particle number */
/* Color schemes */
@ -206,7 +215,7 @@
#define COL_TURBO_CYCLIC 101 /* TURBO color palette (by Anton Mikhailov) corrected to be cyclic, beta */
#define VICSEK_INT ((INTERACTION == I_VICSEK)||(INTERACTION == I_VICSEK_REPULSIVE)||(INTERACTION == I_VICSEK_SPEED))
#define VICSEK_INT ((INTERACTION == I_VICSEK)||(INTERACTION == I_VICSEK_REPULSIVE)||(INTERACTION == I_VICSEK_SPEED)||(INTERACTION == I_VICSEK_SHARK))
typedef struct
{

2
global_particles.c
View File

@ -58,7 +58,7 @@ double x_shooter = -0.2, y_shooter = -0.6, x_target = 0.4, y_target = 0.7;
#define D_ALT_REU 11 /* alternating between star and Reuleaux */
#define D_ANGLE 12 /* angular sector */
#define D_LSHAPE 13 /* L-shaped billiard for conical singularity */
#define D_GENUSN 14 /* polygon with identifies opposite sides */
#define D_GENUSN 14 /* polygon with identified opposite sides */
#define D_PARABOLAS 15 /* polygon with parabolic sides */
#define D_PENROSE 16 /* Penrose solution to illumination problem */

12
global_pdes.c
View File

@ -69,6 +69,7 @@
#define D_MAZE 53 /* maze */
#define D_MAZE_CLOSED 54 /* closed maze */
#define D_MAZE_CHANNELS 541 /* maze with two channels attached */
#define D_MAZE_CHANNELS_INT 542 /* maze with two channels attached, distance defined from interior of cells */
#define D_CHESSBOARD 55 /* chess board configuration */
#define D_TRIANGLE_TILES 56 /* triangular tiling */
#define D_HEX_TILES 57 /* honeycomb tiling */
@ -128,6 +129,8 @@
#define IOR_RANDOM_WELLS 5 /* random superposition of "wells" */
#define IOR_PERIODIC_WELLS_ROTATING 6 /* periodic superposition rotating in time */
#define IOR_PERIODIC_WELLS_ROTATING_LARGE 7 /* periodic superposition rotating in time, larger area */
#define IOR_POISSON_WELLS 8 /* wells located on a random Poisson disc process */
#define IOR_PPP_WELLS 9 /* wells located on a Poisson point process */
/* Boundary conditions */
@ -229,12 +232,17 @@
#define Z_EULER_LPRESSURE 53 /* Laplacian of pressure */
#define Z_EULER_PRESSURE 54 /* pressure */
/* for Euler compressible Euler equation */
/* for compressible Euler equation */
#define Z_EULER_DENSITY 60 /* density */
#define Z_EULER_SPEED 61 /* norm of velocity */
#define Z_EULERC_VORTICITY 62 /* vorticity of velocity */
#define Z_EULER_DIRECTION 63 /* direction of velocity */
#define Z_EULER_DIRECTION_SPEED 64 /* hut for direction of velocity, luminosity for speed */
#define Z_EULER_DIRECTION_SPEED 64 /* hue for direction of velocity, luminosity for speed */
/* for shallow water equation */
#define Z_SWATER_HEIGHT 70 /* height */
#define Z_SWATER_SPEED 71 /* speed */
#define Z_SWATER_DIRECTION_SPEED 74 /* hue for direction of velocity, luminosity for speed */
/* special boundary conditions for Euler equation */
#define BCE_TOPBOTTOM 1 /* special flow at top and bottom */

246
lennardjones.c
View File

@ -36,15 +36,15 @@
#include <omp.h>
#include <time.h>
#define MOVIE 0 /* set to 1 to generate movie */
#define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */
#define MOVIE 0 /* set to 1 to generate movie */
#define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */
#define SAVE_MEMORY 1 /* set to 1 to save memory while saving frames */
#define NO_EXTRA_BUFFER_SWAP 1 /* some OS require one less buffer swap when recording images */
#define TIME_LAPSE 0 /* set to 1 to add a time-lapse movie at the end */
/* so far incompatible with double movie */
#define TIME_LAPSE_FACTOR 3 /* factor of time-lapse movie */
#define TIME_LAPSE_FIRST 0 /* set to 1 to show time-lapse version first */
#define TIME_LAPSE_FIRST 1 /* set to 1 to show time-lapse version first */
#define SAVE_TIME_SERIES 0 /* set to 1 to save time series of particle positions */
@ -58,10 +58,18 @@
#define YMIN -1.125
#define YMAX 1.125 /* y interval for 9/16 aspect ratio */
#define INITXMIN -1.97
#define INITXMAX 1.97 /* x interval for initial condition */
#define INITXMIN -2.0
#define INITXMAX 2.0 /* x interval for initial condition */
#define INITYMIN -1.1
#define INITYMAX 1.1 /* y interval for initial condition */
// #define INITYMAX 7.0 /* y interval for initial condition */
// #define INITYMAX 9.0 /* y interval for initial condition */
#define ADDXMIN -1.97
#define ADDXMAX -0.8 /* x interval for adding particles */
#define ADDYMIN -1.125
#define ADDYMAX 1.125 /* y interval for adding particles */
#define BCXMIN -2.0
#define BCXMAX 2.0 /* x interval for boundary condition */
@ -71,39 +79,42 @@
#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 CIRCLE_PATTERN 1 /* pattern of circles, see list in global_ljones.c */
#define ADD_FIXED_OBSTACLES 0 /* set to 1 do add fixed circular obstacles */
#define OBSTACLE_PATTERN 181 /* pattern of obstacles, see list in global_ljones.c */
#define ADD_INITIAL_PARTICLES 0 /* set to 1 to add a second type of particles */
#define CIRCLE_PATTERN_B 1 /* pattern of circles for additional particles */
#define ADD_FIXED_OBSTACLES 1 /* set to 1 do add fixed circular obstacles */
#define OBSTACLE_PATTERN 4 /* pattern of obstacles, see list in global_ljones.c */
#define ADD_FIXED_SEGMENTS 0 /* set to 1 to add fixed segments as obstacles */
#define SEGMENT_PATTERN 181 /* pattern of repelling segments, see list in global_ljones.c */
#define ROCKET_SHAPE 2 /* shape of rocket combustion chamber, see list in global_ljones.c */
#define ROCKET_SHAPE_B 2 /* shape of second rocket */
#define NOZZLE_SHAPE 2 /* shape of nozzle, see list in global_ljones.c */
#define NOZZLE_SHAPE_B 4 /* shape of nozzle for second rocket, see list in global_ljones.c */
#define SEGMENT_PATTERN 151 /* pattern of repelling segments, see list in global_ljones.c */
#define ROCKET_SHAPE 3 /* shape of rocket combustion chamber, see list in global_ljones.c */
#define ROCKET_SHAPE_B 3 /* shape of second rocket */
#define NOZZLE_SHAPE 6 /* shape of nozzle, see list in global_ljones.c */
#define NOZZLE_SHAPE_B 6 /* shape of nozzle for second rocket, see list in global_ljones.c */
#define TWO_TYPES 0 /* set to 1 to have two types of particles */
#define TYPE_PROPORTION 0.3 /* proportion of particles of first type */
#define TYPE_PROPORTION 0.5 /* 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 */
#define CENTER_PX 0 /* set to 1 to center horizontal momentum */
#define CENTER_PY 0 /* set to 1 to center vertical momentum */
#define CENTER_PANGLE 0 /* set to 1 to center angular momentum */
#define INTERACTION 10 /* particle interaction, see list in global_ljones.c */
#define INTERACTION_B 1 /* particle interaction for second type of particle, see list in global_ljones.c */
#define SPIN_INTER_FREQUENCY 1.0 /* angular frequency of spin-spin interaction */
#define SPIN_INTER_FREQUENCY_B 2.0 /* angular frequency of spin-spin interaction for second particle type */
#define INTERACTION 2 /* particle interaction, see list in global_ljones.c */
#define INTERACTION_B 2 /* particle interaction for second type of particle, see list in global_ljones.c */
#define SPIN_INTER_FREQUENCY 4.0 /* angular frequency of spin-spin interaction */
#define SPIN_INTER_FREQUENCY_B 4.0 /* angular frequency of spin-spin interaction for second particle type */
#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 1.3 /* minimal distance in Poisson disc process, controls density of particles */
#define PDISC_DISTANCE 2.8 /* minimal distance in Poisson disc process, controls density of particles */
#define PDISC_CANDIDATES 100 /* number of candidates in construction of Poisson disc process */
#define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */
#define LAMBDA 0.5 /* parameter controlling the dimensions of domain */
#define MU 0.035 /* parameter controlling radius of particles */
#define MU_B 0.012 /* parameter controlling radius of particles of second type */
#define LAMBDA 0.8 /* parameter controlling the dimensions of domain */
#define MU 0.015 /* parameter controlling radius of particles */
#define MU_B 0.015 /* parameter controlling radius of particles of second type */
#define NPOLY 25 /* number of sides of polygon */
#define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */
#define MDEPTH 4 /* depth of computation of Menger gasket */
@ -111,8 +122,8 @@
#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 120 /* number of grid point for grid of disks */
#define NGRIDY 51 /* number of grid point for grid of disks */
#define NGRIDX 140 /* number of grid point for grid of disks */
#define NGRIDY 70 /* 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 */
#define TWO_CIRCLES_RADIUS_RATIO 0.8 /* ratio of radii for S_TWO_CIRCLES_EXT segment configuration */
@ -125,11 +136,12 @@
/* Parameters for length and speed of simulation */
#define NSTEPS 3800 /* number of frames of movie */
#define NVID 50 /* number of iterations between images displayed on screen */
#define NSTEPS 2800 /* number of frames of movie */
// #define NSTEPS 200 /* number of frames of movie */
#define NVID 60 /* number of iterations between images displayed on screen */
#define NSEG 250 /* number of segments of boundary */
#define INITIAL_TIME 10 /* time after which to start saving frames */
#define OBSTACLE_INITIAL_TIME 200 /* time after which to start moving obstacle */
#define INITIAL_TIME 0 /* time after which to start saving frames */
#define OBSTACLE_INITIAL_TIME 150 /* time after which to start moving obstacle */
#define BOUNDARY_WIDTH 1 /* width of particle boundary */
#define LINK_WIDTH 2 /* width of links between particles */
#define CONTAINER_WIDTH 4 /* width of container boundary */
@ -148,17 +160,20 @@
/* Plot type, see list in global_ljones.c */
#define PLOT 4
#define PLOT_B 6 /* plot type for second movie */
#define PLOT_B 11 /* plot type for second movie */
#define DRAW_BONDS 1 /* set to 1 to draw bonds between neighbours */
#define COLOR_BONDS 1 /* set to 1 to color bonds according to length */
#define FILL_TRIANGLES 0 /* set to 1 to fill triangles between neighbours */
#define ALTITUDE_LINES 0 /* set to 1 to add horizontal lines to show altitude */
#define COLOR_SEG_GROUPS 0 /* set to 1 to collor segment groups differently */
#define COLOR_SEG_GROUPS 1 /* set to 1 to collor segment groups differently */
#define N_PARTICLE_COLORS 200 /* number of colors for P_NUMBER color scheme */
/* Color schemes */
#define COLOR_PALETTE 0 /* Color palette, see list in global_ljones.c */
#define COLOR_PALETTE 18 /* Color palette, see list in global_ljones.c */
#define COLOR_PALETTE_ANGLE 10 /* Color palette for angle representation */
#define COLOR_PALETTE_INITIAL_POS 10 /* Color palette for initial position representation */
#define BLACK 1 /* background */
@ -175,7 +190,7 @@
#define HUEMEAN 220.0 /* mean value of hue for color scheme C_HUE */
#define HUEAMP -50.0 /* amplitude of variation of hue for color scheme C_HUE */
#define PRINT_PARAMETERS 1 /* set to 1 to print certain parameters */
#define PRINT_PARAMETERS 0 /* set to 1 to print certain parameters */
#define PRINT_TEMPERATURE 0 /* set to 1 to print current temperature */
/* particle properties */
@ -184,54 +199,57 @@
#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.2e3 /* energy of particle with hottest color */
#define HUE_TYPE0 260.0 /* hue of particles of type 0 */
#define HUE_TYPE1 210.0 /* hue of particles of type 1 */
#define HUE_TYPE2 160.0 /* hue of particles of type 2 */
#define HUE_TYPE3 290.0 /* hue of particles of type 3 */
#define PARTICLE_EMAX 6000.0 /* energy of particle with hottest color */
#define HUE_TYPE0 60.0 /* hue of particles of type 0 */
#define HUE_TYPE1 280.0 /* hue of particles of type 1 */
#define HUE_TYPE2 140.0 /* hue of particles of type 2 */
#define HUE_TYPE3 200.0 /* hue of particles of type 3 */
#define RANDOM_RADIUS 0 /* set to 1 for random circle radius */
#define DT_PARTICLE 3.0e-6 /* time step for particle displacement */
#define KREPEL 0.1 /* constant in repelling force between particles */
#define EQUILIBRIUM_DIST 2.0 /* Lennard-Jones equilibrium distance */
#define EQUILIBRIUM_DIST_B 2.0 /* Lennard-Jones equilibrium distance for second type of particle */
#define KREPEL 1.0 /* constant in repelling force between particles */
#define EQUILIBRIUM_DIST 5.0 /* Lennard-Jones equilibrium distance */
#define EQUILIBRIUM_DIST_B 5.0 /* Lennard-Jones equilibrium distance for second type of particle */
#define REPEL_RADIUS 15.0 /* radius in which repelling force acts (in units of particle radius) */
#define DAMPING 0.0 /* damping coefficient of particles */
#define INITIAL_DAMPING 50.0 /* damping coefficient of particles during initial phase */
#define DAMPING 100.0 /* damping coefficient of particles */
#define OMEGA_INITIAL 10.0 /* initial angular velocity range */
#define INITIAL_DAMPING 5.0 /* damping coefficient of particles during initial phase */
#define DAMPING_ROT 100.0 /* dampint coefficient for rotation of particles */
#define PARTICLE_MASS 0.25 /* mass of particle of radius MU */
#define PARTICLE_MASS_B 0.5 /* mass of particle of radius MU */
#define PARTICLE_INERTIA_MOMENT 0.5 /* moment of inertia of particle */
#define PARTICLE_MASS 2.5 /* mass of particle of radius MU */
#define PARTICLE_MASS_B 1.0 /* mass of particle of radius MU */
#define PARTICLE_INERTIA_MOMENT 0.2 /* moment of inertia of particle */
#define PARTICLE_INERTIA_MOMENT_B 0.02 /* moment of inertia of second type of particle */
#define V_INITIAL 20.0 /* initial velocity range */
#define OMEGA_INITIAL 5.0 /* initial angular velocity range */
#define V_INITIAL 0.0 /* initial velocity range */
#define OMEGA_INITIAL 10.0 /* initial angular velocity range */
#define VICSEK_VMIN 1.0 /* minimal speed of particles in Vicsek model */
#define VICSEK_VMAX 40.0 /* minimal speed of particles in Vicsek model */
#define THERMOSTAT 0 /* set to 1 to switch on thermostat */
#define VARY_THERMOSTAT 0 /* set to 1 for time-dependent thermostat schedule */
#define SIGMA 5.0 /* noise intensity in thermostat */
#define BETA 0.002 /* initial inverse temperature */
#define BETA 0.02 /* initial inverse temperature */
#define MU_XI 0.01 /* friction constant in thermostat */
#define KSPRING_BOUNDARY 1.0e7 /* confining harmonic potential outside simulation region */
#define KSPRING_OBSTACLE 1.0e11 /* harmonic potential of obstacles */
#define NBH_DIST_FACTOR 3.0 /* radius in which to count neighbours */
#define NBH_DIST_FACTOR 2.5 /* radius in which to count neighbours */
// #define NBH_DIST_FACTOR 4.0 /* radius in which to count neighbours */
#define GRAVITY 0.0 /* gravity acting on all particles */
#define GRAVITY_X 0.0 /* horizontal gravity acting on all particles */
#define GRAVITY_X 5000.0 /* horizontal gravity acting on all particles */
#define INCREASE_GRAVITY 0 /* set to 1 to increase gravity during the simulation */
#define GRAVITY_SCHEDULE 2 /* type of gravity schedule, see list in global_ljones.c */
#define GRAVITY_FACTOR 100.0 /* factor by which to increase gravity */
#define GRAVITY_INITIAL_TIME 200 /* time at start of simulation with constant gravity */
#define GRAVITY_RESTORE_TIME 700 /* time at end of simulation with gravity restored to initial value */
#define KSPRING_VICSEK 0.2 /* spring constant for I_VICSEK_SPEED interaction */
#define VICSEK_REPULSION 10.0 /* repulsion between particles in Vicsek model */
#define ROTATION 1 /* 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 */
#define COUPLE_ANGLE_TO_THERMOSTAT 1 /* set to 1 to couple angular degrees of freedom to thermostat */
#define DIMENSION_FACTOR 1.0 /* scaling factor taking into account number of degrees of freedom */
#define KTORQUE 2.0e3 /* force constant in angular dynamics */
#define KTORQUE 50.0 /* force constant in angular dynamics */
#define KTORQUE_BOUNDARY 1.0e6 /* constant in torque from the boundary */
#define KTORQUE_B 10.0 /* force constant in angular dynamics */
#define KTORQUE_DIFF 150.0 /* force constant in angular dynamics for different particles */
#define KTORQUE_BOUNDARY 1.0e6 /* constant in torque from the boundary */
#define DRAW_SPIN 0 /* set to 1 to draw spin vectors of particles */
#define DRAW_SPIN_B 0 /* set to 1 to draw spin vectors of particles */
#define DRAW_CROSS 1 /* set to 1 to draw cross on particles of second type */
@ -240,7 +258,7 @@
#define QUADRUPOLE_RATIO 0.6 /* anisotropy in quadrupole potential */
#define INCREASE_BETA 0 /* set to 1 to increase BETA during simulation */
#define BETA_FACTOR 0.5 /* factor by which to change BETA during simulation */
#define BETA_FACTOR 0.02 /* factor by which to change BETA during simulation */
#define N_TOSCILLATIONS 1.5 /* number of temperature oscillations in BETA schedule */
#define NO_OSCILLATION 1 /* set to 1 to have exponential BETA change only */
#define MIDDLE_CONSTANT_PHASE 2000 /* final phase in which temperature is constant */
@ -257,7 +275,7 @@
#define CENTER_VIEW_ON_OBSTACLE 0 /* set to 1 to center display on moving obstacle */
#define RESAMPLE_Y 0 /* set to 1 to resample y coordinate of moved particles (for shock waves) */
#define NTRIALS 2000 /* number of trials when resampling */
#define OBSTACLE_RADIUS 0.12 /* radius of obstacle for circle boundary conditions */
#define OBSTACLE_RADIUS 0.35 /* radius of obstacle for circle boundary conditions */
#define FUNNEL_WIDTH 0.25 /* funnel width for funnel boundary conditions */
#define OBSTACLE_XMIN 0.0 /* initial position of obstacle */
#define OBSTACLE_XMAX 3.0 /* final position of obstacle */
@ -267,7 +285,7 @@
#define N_T_AVERAGE 1 /* 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 some particles to thermostat */
#define PARTIAL_THERMO_REGION 4 /* region for partial thermostat coupling (see list in global_ljones.c) */
#define PARTIAL_THERMO_REGION 1 /* region for partial thermostat coupling (see list in global_ljones.c) */
#define PARTIAL_THERMO_SHIFT 0.2 /* distance from obstacle at the right of which particles are coupled to thermostat */
#define PARTIAL_THERMO_WIDTH 0.5 /* vertical size of partial thermostat coupling */
#define PARTIAL_THERMO_HEIGHT 0.25 /* vertical size of partial thermostat coupling */
@ -279,11 +297,11 @@
#define MASS_PART_BOTTOM 10000.0 /* mass of particles at bottom */
#define NPART_BOTTOM 100 /* number of particles at the bottom */
#define ADD_PARTICLES 0 /* set to 1 to add particles */
#define ADD_PARTICLES 0 /* set to 1 to add particles */
#define ADD_TIME 0 /* time at which to add first particle */
#define ADD_PERIOD 10000 /* time interval between adding further particles */
#define N_ADD_PARTICLES 20 /* number of particles to add */
#define FINAL_NOADD_PERIOD 200 /* final period where no particles are added */
#define ADD_PERIOD 10 /* time interval between adding further particles */
#define N_ADD_PARTICLES 2 /* number of particles to add */
#define FINAL_NOADD_PERIOD 0 /* 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 0 /* set to 1 to have a tracer particle */
@ -306,49 +324,51 @@
#define MOVE_BOUNDARY 0 /* set to 1 to move repelling segments, due to force from particles */
#define SEGMENTS_MASS 40.0 /* mass of collection of segments */
#define DEACTIVATE_SEGMENT 1 /* set to 1 to deactivate last segment after a certain time */
#define SEGMENT_DEACTIVATION_TIME 200 /* time at which to deactivate last segment */
#define RELEASE_ROCKET_AT_DEACTIVATION 1 /* set to 1 to limit segments velocity before segment release */
#define SEGMENT_DEACTIVATION_TIME 20 /* time at which to deactivate last segment */
#define RELEASE_ROCKET_AT_DEACTIVATION 0 /* set to 1 to limit segments velocity before segment release */
#define SEGMENTS_X0 1.5 /* initial position of segments */
#define SEGMENTS_Y0 0.0 /* initial position of segments */
#define SEGMENTS_VX0 0.0 /* initial velocity of segments */
#define SEGMENTS_VY0 0.0 /* initial velocity of segments */
#define DAMP_SEGS_AT_NEGATIVE_Y 0 /* set to 1 to dampen segments when y coordinate is negative */
#define MOVE_SEGMENT_GROUPS 1 /* set to 1 to group segments into moving units */
#define SEGMENT_GROUP_MASS 1000.0 /* mass of segment group */
#define MOVE_SEGMENT_GROUPS 0 /* set to 1 to group segments into moving units */
#define SEGMENT_GROUP_MASS 500.0 /* mass of segment group */
#define SEGMENT_GROUP_I 1000.0 /* moment of inertia of segment group */
#define SEGMENT_GROUP_DAMPING 0.0 /* damping of segment groups */
#define GROUP_REPULSION 1 /* set to 1 for groups of segments to repel each other */
#define KSPRING_GROUPS 1.0e11 /* harmonic potential between segment groups */
#define GROUP_REPULSION 0 /* set to 1 for groups of segments to repel each other */
#define KSPRING_GROUPS 5.0e11 /* harmonic potential between segment groups */
#define GROUP_WIDTH 0.05 /* interaction width of groups */
#define GROUP_G_REPEL 1 /* set to 1 to add repulsion between centers of mass of groups */
#define GROUP_G_REPEL 0 /* set to 1 to add repulsion between centers of mass of groups */
#define GROUP_G_REPEL_RADIUS 1.2 /* radius within which centers of mass of groups repel each other */
#define TRACK_SEGMENT_GROUPS 1 /* set to 1 for view to track group of segments */
#define TRACK_SEGMENT_GROUPS 0 /* set to 1 for view to track group of segments */
#define TRACK_X_PADDING 2.0 /* distance from x boundary where tracking starts */
#define POSITION_DEPENDENT_TYPE 0 /* set to 1 to make particle type depend on initial position */
#define POSITION_Y_DEPENDENCE 0 /* set to 1 for the separation between particles to be horizontal */
#define POSITION_DEP_SIGN -1.0 /* sign in position dependence condition */
#define POSITION_DEP_X -0.625 /* threshold value for position-dependent type */
#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_REACTION 15 /* type of reaction, see list in global_ljones.c */
#define RD_REACTION 16 /* type of reaction, see list in global_ljones.c */
#define RD_TYPES 5 /* number of types in reaction-diffusion equation */
#define RD_INITIAL_COND 2 /* initial condition of particles */
#define REACTION_DIST 3.5 /* maximal distance for reaction to occur */
#define REACTION_PROB 0.5 /* probability controlling reaction term */
#define RD_INITIAL_COND 9 /* initial condition of particles */
#define REACTION_DIST 4.0 /* maximal distance for reaction to occur */
#define REACTION_PROB 1.0 /* probability controlling reaction term */
#define DISSOCIATION_PROB 0.002 /* probability controlling dissociation reaction */
#define CENTER_COLLIDED_PARTICLES 0 /* set to 1 to recenter particles upon reaction (may interfere with thermostat) */
#define EXOTHERMIC 0 /* set to 1 to make reaction exo/endothermic */
#define DELTA_EKIN 500.0 /* change of kinetic energy in reaction */
#define EXOTHERMIC 1 /* set to 1 to make reaction exo/endothermic */
#define DELTA_EKIN 2000.0 /* change of kinetic energy in reaction */
#define COLLISION_TIME 15 /* time during which collisions are shown */
#define PRINT_PARTICLE_NUMBER 0 /* set to 1 to print total number of particles */
#define PLOT_PARTICLE_NUMBER 0 /* set to 1 to make of plot of particle number over time */
#define PARTICLE_NB_PLOT_FACTOR 0.5 /* expected final number of particles over initial number */
#define PRINT_LEFT 1 /* set to 1 to print certain parameters at the top left instead of right */
#define PRINT_LEFT 0 /* set to 1 to print certain parameters at the top left instead of right */
#define PLOT_SPEEDS 0 /* set to 1 to add a plot of obstacle speeds (e.g. for rockets) */
#define PLOT_TRAJECTORIES 0 /* set to 1 to add a plot of obstacle trajectories (e.g. for rockets) */
#define VMAX_PLOT_SPEEDS 0.6 /* vertical scale of plot of obstacle speeds */
#define VMAX_PLOT_SPEEDS 0.25 /* vertical scale of plot of obstacle speeds */
#define EHRENFEST_COPY 0 /* set to 1 to add equal number of larger particles (for Ehrenfest model) */
@ -360,20 +380,20 @@
#define WALL_VMAX 100.0 /* max speed of wall */
#define WALL_TIME 0 /* time during which to keep wall */
#define NXMAZE 10 /* width of maze */
#define NYMAZE 10 /* height of maze */
#define NXMAZE 12 /* width of maze */
#define NYMAZE 12 /* height of maze */
#define MAZE_MAX_NGBH 4 /* max number of neighbours of maze cell */
#define RAND_SHIFT 200 /* seed of random number generator */
#define MAZE_XSHIFT 0.0 /* horizontal shift of maze */
#define RAND_SHIFT 4 /* seed of random number generator */
#define MAZE_XSHIFT 0.5 /* horizontal shift of maze */
#define MAZE_WIDTH 0.01 /* width of maze walls */
#define FLOOR_FORCE 1 /* set to 1 to limit force on particle to FMAX */
#define FMAX 1.0e10 /* maximal force */
#define FMAX 1.0e9 /* maximal force */
#define FLOOR_OMEGA 0 /* set to 1 to limit particle momentum to PMAX */
#define PMAX 1000.0 /* maximal force */
#define HASHX 60 /* size of hashgrid in x direction */
#define HASHY 30 /* size of hashgrid in y direction */
#define HASHX 100 /* size of hashgrid in x direction */
#define HASHY 50 /* 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 */
@ -382,6 +402,8 @@
#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 */
#define LIMIT_ENERGY 0 /* set to 1 to limit energy, when there is no thermostat */
#define NO_WRAP_BC ((BOUNDARY_COND != BC_PERIODIC)&&(BOUNDARY_COND != BC_PERIODIC_CIRCLE)&&(BOUNDARY_COND != BC_PERIODIC_TRIANGLE)&&(BOUNDARY_COND != BC_KLEIN)&&(BOUNDARY_COND != BC_PERIODIC_FUNNEL)&&(BOUNDARY_COND != BC_BOY)&&(BOUNDARY_COND != BC_GENUS_TWO))
#define PERIODIC_BC ((BOUNDARY_COND == BC_PERIODIC)||(BOUNDARY_COND == BC_PERIODIC_CIRCLE)||(BOUNDARY_COND == BC_PERIODIC_FUNNEL)||(BOUNDARY_COND == BC_PERIODIC_TRIANGLE))
#define TWO_OBSTACLES ((SEGMENT_PATTERN == S_TWO_CIRCLES_EXT)||(SEGMENT_PATTERN == S_TWO_ROCKETS))
@ -826,7 +848,8 @@ void evolve_segments(t_segment segment[NMAXSEGMENTS], int time)
if (y < YMIN + padding) fy[group] += KSPRING_BOUNDARY*(YMIN + padding - y);
else if (y > YMAX - padding) fy[group] -= KSPRING_BOUNDARY*(y - YMAX + padding);
}
else if (BOUNDARY_COND == BC_REFLECT_ABS) /* add force from simulation boundary */
else if ((BOUNDARY_COND == BC_REFLECT_ABS)||(BOUNDARY_COND == BC_REFLECT_ABS_BOTTOM))
/* add force from simulation boundary */
{
y = 0.5*(segment[i].y1 + segment[i].y2);
if (y < YMIN) fy[group] += KSPRING_BOUNDARY*(YMIN - y);
@ -899,7 +922,7 @@ void evolve_segment_groups(t_segment segment[NMAXSEGMENTS], int time, t_group_se
/* new version of evolve_segments that takes the group structure into account */
{
double fx[NMAXGROUPS], fy[NMAXGROUPS], torque[NMAXGROUPS], dx[NMAXGROUPS], dy[NMAXGROUPS], dalpha[NMAXGROUPS];
double x, y, dx0, dy0, padding, proj, distance, f, xx[2], yy[2], xmean = 0.0, ymean = 0.0;
double x, y, dx0, dy0, padding, proj, distance, f, xx[2], yy[2], xmean = 0.0, ymean = 0.0, ymax = 0.0;
int i, j, k, group = 0;
static double maxdepth, saturation_depth, xmax;
static int first = 1;
@ -944,7 +967,8 @@ void evolve_segment_groups(t_segment segment[NMAXSEGMENTS], int time, t_group_se
if (y < YMIN + padding) fy[group] += KSPRING_BOUNDARY*(YMIN + padding - y);
else if (y > YMAX - padding) fy[group] -= KSPRING_BOUNDARY*(y - YMAX + padding);
}
else if (BOUNDARY_COND == BC_REFLECT_ABS) /* add force from simulation boundary */
else if ((BOUNDARY_COND == BC_REFLECT_ABS)||(BOUNDARY_COND == BC_REFLECT_ABS_BOTTOM))
/* add force from simulation boundary */
{
y = 0.5*(segment[i].y1 + segment[i].y2);
if (y < YMIN) fy[group] += KSPRING_BOUNDARY*(YMIN - y);
@ -1047,10 +1071,14 @@ void evolve_segment_groups(t_segment segment[NMAXSEGMENTS], int time, t_group_se
{
xmean += segment_group[group].xc;
ymean += segment_group[group].yc;
if (segment_group[group].yc > ymax) ymax = segment_group[group].yc;
}
xmean = xmean/((double)(ngroups-1));
ymean = ymean/((double)(ngroups-1));
/* bias towards ymax */
ymean = 0.75*ymax + 0.25*ymean;
if (ymean > ytrack) ytrack = ymean;
if (xmean > xmax)
xtrack = xmean - xmax;
@ -1063,7 +1091,7 @@ void evolve_segment_groups(t_segment segment[NMAXSEGMENTS], int time, t_group_se
void animation()
{
double time, scale, diss, rgb[3], dissip, gradient[2], x, y, dx, dy, dt, xleft, xright,
a, b, length, fx, fy, force[2], totalenergy = 0.0, krepel = KREPEL, pos[2], prop, vx, beta = BETA, xi = 0.0, xmincontainer = BCXMIN, xmaxcontainer = BCXMAX, torque, torque_ij, fboundary = 0.0, pleft = 0.0, pright = 0.0, entropy[2], mean_energy, gravity = GRAVITY, speed_ratio, ymin, ymax, delta_energy, speed, ratio = 1.0, ratioc;
a, b, length, fx, fy, force[2], totalenergy = 0.0, krepel = KREPEL, pos[2], prop, vx, beta = BETA, xi = 0.0, xmincontainer = BCXMIN, xmaxcontainer = BCXMAX, torque, torque_ij, fboundary = 0.0, pleft = 0.0, pright = 0.0, entropy[2], mean_energy, gravity = GRAVITY, speed_ratio, ymin, ymax, delta_energy, speed, ratio = 1.0, ratioc, cum_etot = 0.0, emean = 0.0;
double *qx, *qy, *px, *py, *qangle, *pangle, *pressure, *obstacle_speeds;
int i, j, k, n, m, s, ij[2], i0, iplus, iminus, j0, jplus, jminus, p, q, p1, q1, p2, q2, total_neighbours = 0,
min_nb, max_nb, close, wrapx = 0, wrapy = 0, nactive = 0, nadd_particle = 0, nmove = 0, nsuccess = 0,
@ -1117,13 +1145,6 @@ void animation()
lj_final_position = fopen("lj_final_position.dat", "w");
}
/* initialise positions and radii of circles */
init_particle_config(particle);
init_hashgrid(hashgrid);
xshift = OBSTACLE_XMIN;
speed_ratio = (double)(25*NVID)*DT_PARTICLE;
if (ADD_FIXED_OBSTACLES) init_obstacle_config(obstacle);
if (ADD_FIXED_SEGMENTS) init_segment_config(segment);
@ -1133,6 +1154,17 @@ void animation()
group_speeds = (t_group_data *)malloc(ngroups*(INITIAL_TIME + NSTEPS)*sizeof(t_group_data));
}
/* initialise positions and radii of circles */
init_particle_config(particle);
/* add some particles, beta */
if (ADD_INITIAL_PARTICLES) add_particle_config(particle, -0.6, 1.6, -1.0, 1.0, MU_B);
init_hashgrid(hashgrid);
xshift = OBSTACLE_XMIN;
speed_ratio = (double)(25*NVID)*DT_PARTICLE;
if (RECORD_PRESSURES) for (i=0; i<N_PRESSURES; i++) pressure[i] = 0.0;
if (PLOT_SPEEDS) obstacle_speeds = (double *)malloc(2*ngroups*(INITIAL_TIME + NSTEPS)*sizeof(double));
if (PLOT_PARTICLE_NUMBER)
@ -1251,15 +1283,19 @@ void animation()
{
speed = module2(particle[j].vx,particle[j].vy);
if ((VICSEK_VMIN > 0.0)&&(speed < VICSEK_VMIN)) speed = VICSEK_VMIN;
if ((INTERACTION == I_VICSEK_SHARK)&&(particle[j].type == 2)) speed *= 1.75;
if (speed > VICSEK_VMAX) speed = 0.5*(speed + VICSEK_VMAX);
particle[j].vx = speed*cos(particle[j].angle);
particle[j].vy = speed*sin(particle[j].angle);
speed = module2(px[j],py[j]);
if ((VICSEK_VMIN > 0.0)&&(speed < VICSEK_VMIN)) speed = VICSEK_VMIN;
if ((INTERACTION == I_VICSEK_SHARK)&&(particle[j].type == 2)) speed *= 1.75;
if (speed > VICSEK_VMAX) speed = 0.5*(speed + VICSEK_VMAX);
px[j] = speed*cos(particle[j].angle);
py[j] = speed*sin(particle[j].angle);
}
/* add gravity */
@ -1284,6 +1320,7 @@ void animation()
/* timestep of thermostat algorithm */
totalenergy = evolve_particles(particle, hashgrid, qx, qy, qangle, px, py, pangle, beta, &nactive, &nsuccess, &nmove, i < INITIAL_TIME);
/* evolution of lid coordinate */
if (BOUNDARY_COND == BC_RECTANGLE_LID) evolve_lid(fboundary);
if (BOUNDARY_COND == BC_RECTANGLE_WALL)
@ -1378,6 +1415,25 @@ void animation()
}
printf("Mean kinetic energy: %.3f\n", totalenergy/(double)ncircles);
if ((!THERMOSTAT)&&(LIMIT_ENERGY))
{
if (cum_etot > 0.0)
{
emean = cum_etot/(double)(i+1);
if (totalenergy > 10.0*emean)
{
reset_energy(particle, px, py, totalenergy, emean);
totalenergy = 0.0;
for (j=0; j<ncircles; j++) if (particle[j].active)
totalenergy += particle[i].energy;
printf("Reset mean kinetic energy: %.3f\n", totalenergy/(double)ncircles);
}
}
printf("Emean: %.3f\n", emean/(double)ncircles);
cum_etot += totalenergy;
}
printf("Boundary force: %.3f\n", fboundary/(double)(ncircles*NVID));
if (RESAMPLE_Y) printf("%i succesful moves out of %i trials\n", nsuccess, nmove);
if (INCREASE_GRAVITY) printf("Gravity: %.3f\n", gravity);
@ -1467,7 +1523,7 @@ void animation()
{
if (i >= INITIAL_TIME)
{
if (TIME_LAPSE_FIRST)
if ((TIME_LAPSE)&&(TIME_LAPSE_FIRST))
{
if ((TIME_LAPSE)&&((i - INITIAL_TIME)%TIME_LAPSE_FACTOR == 0)&&(!DOUBLE_MOVIE))
{

33
particle_billiard.c
View File

@ -35,15 +35,18 @@
#include <omp.h>
#include <time.h>
#define MOVIE 0 /* set to 1 to generate movie */
#define MOVIE 1 /* set to 1 to generate movie */
#define SAVE_MEMORY 1 /* set to 1 to save memory when writing tiff images */
#define INVERT_COUNTER 0 /* set to 1 to save frames in inverse order */
#define WINWIDTH 1280 /* window width */
// #define WINWIDTH 1280 /* window width */
#define WINWIDTH 720 /* window width */
#define WINHEIGHT 720 /* window height */
#define XMIN -1.5
#define XMAX 2.5 /* x interval */
// #define XMIN -1.5
// #define XMAX 2.5 /* x interval */
#define XMIN -1.125
#define XMAX 1.125 /* x interval */
#define YMIN -1.125
#define YMAX 1.125 /* y interval for 9/16 aspect ratio */
@ -51,7 +54,7 @@
/* Choice of the billiard table, see global_particles.c */
#define B_DOMAIN 31 /* choice of domain shape */
#define B_DOMAIN 30 /* choice of domain shape */
#define CIRCLE_PATTERN 1 /* pattern of circles */
#define POLYLINE_PATTERN 10 /* pattern of polyline */
@ -83,7 +86,7 @@
/* Simulation parameters */
// #define NPART 10 /* number of particles */
#define NPART 5000 /* number of particles */
#define NPART 50000 /* number of particles */
#define NPARTMAX 100000 /* maximal number of particles after resampling */
#define LMAX 0.01 /* minimal segment length triggering resampling */
#define DMIN 0.02 /* minimal distance to boundary for triggering resampling */
@ -91,10 +94,10 @@
#define SHOWTRAILS 0 /* set to 1 to keep trails of the particles */
#define HEATMAP 1 /* set to 1 to show heat map of particles */
#define DRAW_FINAL_HEATMAP 1 /* set to 1 to show final heat map of particles */
#define DRAW_HEATMAP_HISTOGRAM 1 /* set to 1 to draw a histogram of particle distribution in heat map */
#define NBIN_FACTOR 8.0 /* constant controlling number of bins in histogram */
#define DRAW_HEATMAP_HISTOGRAM 0 /* set to 1 to draw a histogram of particle distribution in heat map */
#define NBIN_FACTOR 6.0 /* constant controlling number of bins in histogram */
#define DRAW_HEATMAP_PARTICLES 1 /* set to 1 to draw particles in heat map */
#define HEATMAP_MAX_PART_BY_CELL 5 /* set to positive value to draw only limited number of particles in cell */
#define HEATMAP_MAX_PART_BY_CELL 50 /* set to positive value to draw only limited number of particles in cell */
#define PLOT_HEATMAP_AVERAGE 1 /* set to 1 to plot average number of particles in heat map */
#define SHOWZOOM 0 /* set to 1 to show zoom on specific area */
#define PRINT_PARTICLE_NUMBER 0 /* set to 1 to print number of particles */
@ -105,8 +108,8 @@
#define TEST_INITIAL_COND 0 /* set to 1 to allow only initial conditions that pass a test */
#define NSTEPS 22000 /* number of frames of movie */
#define TIME 2000 /* time between movie frames, for fluidity of real-time simulation */
#define NSTEPS 1300 /* number of frames of movie */
#define TIME 3000 /* time between movie frames, for fluidity of real-time simulation */
// #define DPHI 0.000002 /* integration step */
#define DPHI 0.00002 /* integration step */
#define NVID 25 /* number of iterations between images displayed on screen */
@ -120,12 +123,12 @@
/* Colors and other graphical parameters */
#define COLOR_PALETTE 17 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE 11 /* Color palette, see list in global_pdes.c */
#define NCOLORS 500 /* number of colors */
#define COLORSHIFT 0 /* hue of initial color */
#define COLOR_HUEMIN 0 /* minimal color hue */
#define COLOR_HUEMAX 160 /* maximal color hue */
#define COLOR_HUEMAX 300 /* maximal color hue */
#define RAINBOW_COLOR 1 /* set to 1 to use different colors for all particles */
#define FLOWER_COLOR 0 /* set to 1 to adapt initial colors to flower billiard (tracks vs core) */
#define NSEG 100 /* number of segments of boundary */
@ -145,8 +148,8 @@
#define SLEEP1 1 /* initial sleeping time */
#define SLEEP2 1 /* final sleeping time */
#define NXMAZE 36 /* width of maze */
#define NYMAZE 36 /* height of maze */
#define NXMAZE 18 /* width of maze */
#define NYMAZE 18 /* height of maze */
#define MAZE_MAX_NGBH 8 /* max number of neighbours of maze cell */
#define RAND_SHIFT 15 /* seed of random number generator */
#define MAZE_XSHIFT 0.0 /* horizontal shift of maze */

621
rde.c
View File

@ -48,14 +48,21 @@
#define WINWIDTH 1920 /* window width */
#define WINHEIGHT 1150 /* window height */
#define NX 960 /* number of grid points on x axis */
#define NY 575 /* number of grid points on y axis */
// #define NY 575 /* number of grid points on y axis */
// #define NX 960 /* number of grid points on x axis */
// #define NY 575 /* number of grid points on y axis */
// #define WINWIDTH 1280 /* window width */
// #define WINHEIGHT 720 /* window height */
// #define NX 640 /* number of grid points on x axis */
// #define NY 360 /* 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 320 /* number of grid points on x axis */
// #define NY 180 /* number of grid points on y axis */
// #define XMIN -1.2
// #define XMAX 1.2 /* x interval */
// #define YMIN -1.2
// #define YMAX 1.2 /* y interval for 9/16 aspect ratio */
#define XMIN -2.0
#define XMAX 2.0 /* x interval */
#define YMIN -1.197916667
@ -63,7 +70,7 @@
/* Choice of simulated equation */
#define RDE_EQUATION 7 /* choice of reaction term, see list in global_3d.c */
#define RDE_EQUATION 8 /* choice of reaction term, see list in global_3d.c */
#define NFIELDS 3 /* number of fields in reaction-diffusion equation */
#define NLAPLACIANS 0 /* number of fields for which to compute Laplacian */
@ -72,13 +79,16 @@
#define ADD_FORCE_FIELD 0 /* set to 1 to add a foce field (for compressible Euler equation) */
#define POTENTIAL 7 /* type of potential or vector potential, see list in global_3d.c */
#define FORCE_FIELD 5 /* type of force field, see list in global_3d.c */
#define ADD_CORIOLIS_FORCE 1 /* set to 1 to add Coriolis force (quasigeostrophic Euler equations) */
#define ADD_CORIOLIS_FORCE 0 /* set to 1 to add Coriolis force (quasigeostrophic Euler equations) */
#define VARIABLE_DEPTH 1 /* set to 1 for variable depth in shallow water equation */
#define SWATER_DEPTH 4 /* variable depth in shallow water equation */
#define ANTISYMMETRIZE_WAVE_FCT 0 /* set tot 1 to make wave function antisymmetric */
#define ADAPT_STATE_TO_BC 0 /* to smoothly adapt initial state to obstacles */
#define OBSTACLE_GEOMETRY 72 /* geometry of obstacles, as in B_DOMAIN */
#define BC_STIFFNESS 50.0 /* controls region of boundary condition control */
#define OBSTACLE_GEOMETRY 7 /* geometry of obstacles, as in B_DOMAIN */
// #define BC_STIFFNESS 100.0 /* controls region of boundary condition control */
#define BC_STIFFNESS 50.0 /* controls region of boundary condition control */
#define CHECK_INTEGRAL 1 /* set to 1 to check integral of first field */
#define JULIA_SCALE 0.5 /* scaling for Julia sets */
@ -86,14 +96,14 @@
#define B_DOMAIN 999 /* choice of domain shape, see list in global_pdes.c */
#define CIRCLE_PATTERN 99 /* pattern of circles, see list in global_pdes.c */
#define CIRCLE_PATTERN 8 /* pattern of circles, see list in global_pdes.c */
#define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */
#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.6 /* parameter controlling the dimensions of domain */
#define MU 0.08 /* parameter controlling the dimensions of domain */
#define LAMBDA 0.5 /* parameter controlling the dimensions of domain */
#define MU 0.06 /* parameter controlling the dimensions of domain */
#define NPOLY 5 /* number of sides of polygon */
#define APOLY 2.0 /* angle by which to turn polygon, in units of Pi/2 */
#define MDEPTH 7 /* depth of computation of Menger gasket */
@ -101,8 +111,8 @@
#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 NGRIDX 6 /* number of grid point for grid of disks */
#define NGRIDY 8 /* number of grid point for grid of disks */
#define REVERSE_TESLA_VALVE 1 /* set to 1 to orient Tesla valve in blocking configuration */
#define X_SHOOTER -0.2
@ -123,7 +133,12 @@
#define DT 0.00000025
#define VISCOSITY 2.0
// #define VISCOSITY 0.0001
#define VISCOSITY 1.5e-5
// #define VISCOSITY 5.0e-4
// #define VISCOSITY 1.0e-3
#define DISSIPATION 1.0e-8
// #define DISSIPATION 1.0e-7
#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 */
@ -138,9 +153,8 @@
#define BZQ 0.0008 /* parameter in BZ equation */
#define BZF 1.2 /* parameter in BZ equation */
#define B_FIELD 10.0 /* magnetic field */
#define G_FIELD 0.004 /* gravity/constant in repulsive field from obstacles */
// #define G_FIELD 1.0e-4 /* gravity/constant in repulsive field from obstacles */
// #define G_FIELD 2.0e-5 /* gravity/constant in repulsive field from obstacles */
#define G_FIELD 0.75 /* gravity */
#define BC_FIELD 1.0e-5 /* constant in repulsive field from obstacles */
#define AB_RADIUS 0.2 /* radius of region with magnetic field for Aharonov-Bohm effect */
#define K_EULER 50.0 /* constant in stream function integration of Euler equation */
#define K_EULER_INC 0.5 /* constant in incompressible Euler equation */
@ -181,11 +195,22 @@
#define IN_OUT_BC_FACTOR 0.001 /* factor of convex combination between old and new flow */
#define BC_FLOW_TYPE 1 /* type of initial condition */
/* see list in global_pdes.c */
#define IN_OUT_FLOW_MIN_AMP 0.45 /* amplitude of in-flow/out-flow boundary conditions (for Euler equation) - min value */
#define IN_OUT_FLOW_AMP 0.45 /* amplitude of in-flow/out-flow boundary conditions (for Euler equation) - max value */
#define IN_OUT_FLOW_MIN_AMP 0.25 /* amplitude of in-flow/out-flow boundary conditions (for Euler equation) - min value */
#define IN_OUT_FLOW_AMP 0.25 /* amplitude of in-flow/out-flow boundary conditions (for Euler equation) - max value */
#define LAMINAR_FLOW_MODULATION 0.01 /* asymmetry of laminar flow */
#define LAMINAR_FLOW_YPERIOD 1.0 /* period of laminar flow in y direction */
#define PRESSURE_GRADIENT 0.3 /* amplitude of pressure gradient for Euler equation */
#define PRESSURE_GRADIENT 0.2 /* amplitude of pressure gradient for Euler equation */
#define SWATER_MIN_HEIGHT 2.0 /* min height of initial condition for shallow water equation */
// #define DEPTH_FACTOR 0.0125 /* proportion of min height in variable depth */
// #define DEPTH_FACTOR 0.001 /* proportion of min height in variable depth */
// #define DEPTH_FACTOR 0.0012 /* proportion of min height in variable depth */
// #define DEPTH_FACTOR 0.0015 /* proportion of min height in variable depth */
// #define DEPTH_FACTOR 0.002 /* proportion of min height in variable depth */
// #define DEPTH_FACTOR 0.005 /* proportion of min height in variable depth */
// #define DEPTH_FACTOR 0.01 /* proportion of min height in variable depth */
#define DEPTH_FACTOR 0.015 /* proportion of min height in variable depth */
#define TANH_FACTOR 1.0 /* steepness of variable depth */
#define EULER_GRADIENT_YSHIFT 0.0 /* y-shift in computation of gradient in Euler equation */
@ -193,15 +218,22 @@
#define B_COND 1
#define B_COND_LEFT 0
#define B_COND_RIGHT 0
#define B_COND_TOP 1
#define B_COND_BOTTOM 1
/* Parameters for length and speed of simulation */
#define NSTEPS 2250 /* number of frames of movie */
// #define NSTEPS 500 /* number of frames of movie */
#define NVID 100 /* number of iterations between images displayed on screen */
#define NSTEPS 2100 /* number of frames of movie */
// #define NSTEPS 500 /* number of frames of movie */
// #define NVID 100 /* number of iterations between images displayed on screen */
#define NVID 30 /* 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 NSEG 999 /* number of segments of boundary */
#define BOUNDARY_WIDTH 2 /* width of billiard boundary */
#define PAUSE 100 /* number of frames after which to pause */
@ -215,7 +247,7 @@
/* Visualisation */
#define PLOT_3D 0 /* controls whether plot is 2D or 3D */
#define PLOT_3D 1 /* controls whether plot is 2D or 3D */
#define ROTATE_VIEW 0 /* set to 1 to rotate position of observer */
#define ROTATE_ANGLE 360.0 /* total angle of rotation during simulation */
@ -224,22 +256,27 @@
/* Plot type - color scheme */
#define CPLOT 64
#define CPLOT_B 62
#define CPLOT 70
#define CPLOT_B 74
/* Plot type - height of 3D plot */
#define ZPLOT 61 /* z coordinate in 3D plot */
#define ZPLOT_B 64 /* z coordinate in second 3D plot */
#define ZPLOT 70 /* z coordinate in 3D plot */
#define ZPLOT_B 71 /* 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 FADE_WATER_DEPTH 0 /* set to 1 to make wave color depth-dependent */
#define DRAW_OUTSIDE_GRAY 0 /* experimental - draw outside of billiard in gray */
#define ADD_POTENTIAL_TO_Z 1 /* set to 1 to add the external potential to z-coordinate of plot */
#define ADD_POT_CONSTANT 0.35 /* constant in front of added potential */
#define ADD_POTENTIAL_TO_Z 0 /* set to 1 to add the external potential to z-coordinate of plot */
#define ADD_POT_CONSTANT 0.35 /* constant added to wave height */
#define ADD_HEIGHT_CONSTANT -2.0 /* constant added to wave height */
#define DRAW_DEPTH 1 /* set to 1 to draw water depth */
#define DEPTH_SCALE 0.75 /* vertical scaling of depth plot */
#define DEPTH_SHIFT -0.015 /* vertical shift of depth plot */
#define PLOT_SCALE_ENERGY 0.05 /* vertical scaling in energy plot */
@ -269,26 +306,28 @@
/* Color schemes, see list in global_pdes.c */
#define COLOR_PALETTE 17 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 10 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE 11 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 10 /* Color palette, see list in global_pdes.c */
// #define COLOR_PALETTE_B 0 /* 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 0.0 /* phase shift of color scheme, in units of Pi */
#define COLOR_PHASE_SHIFT 0.5 /* 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 15.0 /* additional scaling factor for color scheme P_3D_AMPLITUDE */
#define VSCALE_AMPLITUDE 10.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 SLOPE_SCHROD_LUM 200.0 /* sensitivity of luminosity on module, for color scheme Z_ARGUMENT */
#define MIN_SCHROD_LUM 0.2 /* minimal luminosity in color scheme Z_ARGUMENT*/
#define VSCALE_PRESSURE 2.0 /* additional scaling factor for color scheme Z_EULER_PRESSURE */
#define PRESSURE_SHIFT 10.0 /* shift for color scheme Z_EULER_PRESSURE */
#define PRESSURE_LOG_SHIFT -2.5 /* shift for color scheme Z_EULER_PRESSURE */
#define VSCALE_WATER_HEIGHT 0.4 /* vertical scaling of water height */
#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 */
@ -301,25 +340,28 @@
#define LOG_SCALE 0.5 /* scaling factor for energy log representation */
#define LOG_SHIFT 1.0
#define LOG_MIN 1.0e-3 /* floor value for log vorticity plot */
#define VSCALE_SPEED 10.0 /* additional scaling factor for color scheme Z_EULER_SPEED */
#define VSCALE_SPEED 200.0 /* additional scaling factor for color scheme Z_EULER_SPEED */
#define VMEAN_SPEED 0.0 /* mean value around which to scale for color scheme Z_EULER_SPEED */
#define SHIFT_DENSITY 2.6 /* shift for color scheme Z_EULER_DENSITY */
#define VSCALE_DENSITY 7.5 /* additional scaling factor for color scheme Z_EULER_DENSITY */
#define SHIFT_DENSITY 8.5 /* shift for color scheme Z_EULER_DENSITY */
#define VSCALE_DENSITY 3.0 /* additional scaling factor for color scheme Z_EULER_DENSITY */
#define VSCALE_VORTICITY 20.0 /* additional scaling factor for color scheme Z_EULERC_VORTICITY */
#define VORTICITY_SHIFT 0.0 /* vertical shift of vorticity */
#define ZSCALE_SPEED 1.0 /* additional scaling factor for z-coord Z_EULER_SPEED */
#define ZSCALE_SPEED 0.3 /* additional scaling factor for z-coord Z_EULER_SPEED and Z_SWATER_SPEED */
#define VSCALE_SWATER 300.0 /* additional scaling factor for color scheme Z_EULER_DENSITY */
#define NXMAZE 13 /* width of maze */
#define NYMAZE 13 /* height of maze */
#define NXMAZE 7 /* width of maze */
#define NYMAZE 7 /* height of maze */
#define MAZE_MAX_NGBH 4 /* max number of neighbours of maze cell */
#define RAND_SHIFT 3 /* seed of random number generator */
#define MAZE_XSHIFT 0.0 /* horizontal shift of maze */
#define MAZE_WIDTH 0.03 /* half width of maze walls */
#define MAZE_WIDTH 0.04 /* half width of maze walls */
#define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */
#define COLORBAR_RANGE 2.5 /* scale of color scheme bar */
#define COLORBAR_RANGE_B 2.5 /* scale of color scheme bar for 2nd part */
#define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */
#define CIRC_COLORBAR 0 /* set to 1 to draw circular color scheme */
#define CIRC_COLORBAR_B 1 /* set to 1 to draw circular color scheme */
/* only for compatibility with wave_common.c */
#define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */
@ -355,23 +397,23 @@ 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] = {8.0, 8.0, 10.0}; /* location of observer for REP_PROJ_3D representation */
double observer[3] = {8.0, 8.0, 7.0}; /* location of observer for REP_PROJ_3D representation */
int reset_view = 0; /* switch to reset 3D view parameters (for option ROTATE_VIEW) */
#define Z_SCALING_FACTOR 1.25 /* overall scaling factor of z axis for REP_PROJ_3D representation */
#define Z_SCALING_FACTOR 35.0 /* overall scaling factor of z axis for REP_PROJ_3D representation */
#define XY_SCALING_FACTOR 1.7 /* overall scaling factor for on-screen (x,y) coordinates after projection */
#define ZMAX_FACTOR 1.0 /* max value of z coordinate for REP_PROJ_3D representation */
#define XSHIFT_3D 0.0 /* overall x shift for REP_PROJ_3D representation */
#define YSHIFT_3D -0.1 /* overall y shift for REP_PROJ_3D representation */
#define YSHIFT_3D 0.1 /* overall y shift for REP_PROJ_3D representation */
#define BORDER_PADDING 0 /* distance from boundary at which to plot points, to avoid boundary effects due to gradient */
/* For debugging purposes only */
#define FLOOR 1 /* set to 1 to limit wave amplitude to VMAX */
#define VMAX 1000.0 /* max value of wave amplitude */
#define VMAX 100.0 /* max value of wave amplitude */
#define TEST_GRADIENT 0 /* print norm squared of gradient */
#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)||(cplot == Z_ANGLE_GRADIENTX)||(cplot == Z_ARGUMENT)||(cplot == Z_ANGLE_GRADIENTX)))
#define COMPUTE_WRAP_ANGLE ((WRAP_ANGLE)&&((cplot == Z_ANGLE_GRADIENT)||(cplot == Z_ANGLE_GRADIENTX)||(cplot == Z_ARGUMENT)||(cplot == Z_ANGLE_GRADIENTX)||(cplot == Z_EULER_DIRECTION_SPEED)||(cplot == Z_SWATER_DIRECTION_SPEED)))
#define PRINT_PARAMETERS ((PRINT_TIME)||(PRINT_VISCOSITY)||(PRINT_RPSLZB)||(PRINT_PROBABILITIES)||(PRINT_NOISE)||(PRINT_FLOW_SPEED)||(PRINT_AVERAGE_SPEED))
#define COMPUTE_PRESSURE ((ZPLOT == Z_EULER_PRESSURE)||(CPLOT == Z_EULER_PRESSURE)||(ZPLOT_B == Z_EULER_PRESSURE)||(CPLOT_B == Z_EULER_PRESSURE))
@ -519,16 +561,16 @@ void compute_gfield(int i, int j, double *gx, double *gy)
{
r = module2(x,y) + 1.0e-2;
f = 0.5*(1.0 - tanh(BC_STIFFNESS*(r - LAMBDA)));
*gx = G_FIELD*f*x/r;
*gy = G_FIELD*f*y/r;
*gx = BC_FIELD*f*x/r;
*gy = BC_FIELD*f*y/r;
break;
}
case (GF_ELLIPSE):
{
r = module2(x/LAMBDA,y/MU) + 1.0e-2;
f = 0.5*(1.0 - tanh(BC_STIFFNESS*(r - 1.0)));
*gx = G_FIELD*f*x/(LAMBDA*LAMBDA);
*gy = G_FIELD*f*y/(MU*MU);
*gx = BC_FIELD*f*x/(LAMBDA*LAMBDA);
*gy = BC_FIELD*f*y/(MU*MU);
break;
}
case (GF_AIRFOIL):
@ -536,8 +578,8 @@ void compute_gfield(int i, int j, double *gx, double *gy)
y1 = y + a*x*x;
r = module2(x/LAMBDA,y1/MU) + 1.0e-2;
f = 0.5*(1.0 - tanh(BC_STIFFNESS*(r - 1.0)));
*gx = G_FIELD*f*(x/(LAMBDA*LAMBDA) + a*y1/(MU*MU));
*gy = G_FIELD*f*y1/(MU*MU);
*gx = BC_FIELD*f*(x/(LAMBDA*LAMBDA) + a*y1/(MU*MU));
*gy = BC_FIELD*f*y1/(MU*MU);
break;
}
case (GF_WING):
@ -554,8 +596,8 @@ void compute_gfield(int i, int j, double *gx, double *gy)
y1 = y + a*x*x;
r = module2(x/LAMBDA,y1/(MU*x1)) + 1.0e-2;
f = 0.5*(1.0 - tanh(BC_STIFFNESS*(r - 1.0)));
*gx = G_FIELD*f*(x/(LAMBDA*LAMBDA) + 2.0*a*x*y1/(MU*MU*x1*x1) - y1*y1/(MU*MU*x1*x1*x1));
*gy = G_FIELD*f*y1/(MU*MU*x1*x1);
*gx = BC_FIELD*f*(x/(LAMBDA*LAMBDA) + 2.0*a*x*y1/(MU*MU*x1*x1) - y1*y1/(MU*MU*x1*x1*x1));
*gy = BC_FIELD*f*y1/(MU*MU*x1*x1);
// *gx = 0.1*G_FIELD*f*(x/(LAMBDA*LAMBDA) + 2.0*a*x*y1/(MU*MU*x1*x1) - y1*y1/(MU*MU*x1*x1*x1));
// *gy = 0.1*G_FIELD*f*y1/(MU*MU*x1*x1);
// hx = x/(LAMBDA*LAMBDA) + 2.0*a*x*y1/(MU*MU*x1*x1) - y1*y1/(MU*MU*x1*x1*x1);
@ -613,8 +655,8 @@ void initialize_gfield(double gfield[2*NX*NY], double bc_field[NX*NY], double bc
#pragma omp parallel for private(i,j)
for (i=1; i<NX-1; i++){
for (j=1; j<NY-1; j++){
gfield[i*NY+j] = G_FIELD*(bc_field2[(i+1)*NY+j] - bc_field2[(i-1)*NY+j])/dx;
gfield[NX*NY+i*NY+j] = G_FIELD*(bc_field2[i*NY+j+1] - bc_field2[i*NY+j-1])/dy;
gfield[i*NY+j] = BC_FIELD*(bc_field2[(i+1)*NY+j] - bc_field2[(i-1)*NY+j])/dx;
gfield[NX*NY+i*NY+j] = BC_FIELD*(bc_field2[i*NY+j+1] - bc_field2[i*NY+j-1])/dy;
// printf("gfield at (%i,%i): (%.3lg, %.3lg)\n", i, j, gfield[i*NY+j], gfield[NX*NY+i*NY+j]);
}
}
@ -647,14 +689,232 @@ void initialize_gfield(double gfield[2*NX*NY], double bc_field[NX*NY], double bc
}
}
void evolve_wave_half(double *phi_in[NFIELDS], double *phi_out[NFIELDS], short int xy_in[NX*NY],
double potential_field[NX*NY], double vector_potential_field[2*NX*NY],
double gfield[2*NX*NY], t_rde rde[NX*NY])
void compute_water_depth(int i, int j, t_rde *rde, int ncircles)
/* initialize the vector potential, for Schrodinger equation in a magnetic field */
{
double x, y, xy[2], r0, r, z, h, h0, hh, x1, y1, d;
int n, nx, ny;
ij_to_xy(i, j, xy);
x = xy[0];
y = xy[1];
switch (SWATER_DEPTH) {
case (SH_CIRCLE):
{
r0 = module2(x,y);
r = r0/LAMBDA;
h = tanh(TANH_FACTOR*(r-1.0));
hh = 1.0 - h*h;
z = SWATER_MIN_HEIGHT*DEPTH_FACTOR*TANH_FACTOR*hh/(r0*LAMBDA);
rde->depth = SWATER_MIN_HEIGHT*DEPTH_FACTOR*h;
rde->gradx = x*z;
rde->grady = y*z;
break;
}
case (SH_CIRCLES):
{
h = 0.0;
z = 0.0;
r = 10.0;
/* compute minimal distance to circles */
for (n = 0; n < ncircles; n++)
for (nx = -1; nx < 2; nx++)
for (ny = -1; ny < 2; ny++)
{
x1 = circles[n].xc + (double)nx*(XMAX-XMIN);
y1 = circles[n].yc + (double)ny*(YMAX-YMIN);
r0 = module2(x - x1,y - y1)/circles[n].radius;
if (r0 < r) r = r0;
}
h = tanh(TANH_FACTOR*(r-1.0));
h = 0.5*(h + 1.0);
rde->depth = SWATER_MIN_HEIGHT*DEPTH_FACTOR*h;
break;
}
case (SH_COAST):
{
x1 = PI*(x + 0.1 - 0.1*cos(PI*y/YMAX))/XMAX;
d = 3.0*sin(x1);
h = tanh(-TANH_FACTOR*d);
h = 0.5*(h + 1.0);
rde->depth = SWATER_MIN_HEIGHT*DEPTH_FACTOR*h;
break;
}
case (SH_COAST_MONOTONE):
{
x1 = PI*(x + 0.1 - 0.2*cos(PI*y/YMAX))/XMAX;
h = tanh(-TANH_FACTOR*x1);
h = 0.5*(h + 1.0);
rde->depth = SWATER_MIN_HEIGHT*DEPTH_FACTOR*h;
break;
}
}
}
double initialize_water_depth(t_rde rde[NX*NY])
{
int i, j, ncircles;
double dx, dy, min, max, pscal, norm, vz = 0.01;
if (SWATER_DEPTH == SH_CIRCLES) ncircles = init_circle_config_pattern(circles, CIRCLE_PATTERN);
#pragma omp parallel for private(i,j)
for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
compute_water_depth(i, j, &rde[i*NY + j], ncircles);
if (rde[i*NY + j].depth < 0.0) rde[i*NY + j].depth = 0.0;
}
}
dx = (XMAX - XMIN)/(double)NX;
dy = (YMAX - YMIN)/(double)NY;
/* compute x gradient */
#pragma omp parallel for private(i,j)
for (i=1; i<NX-1; i++){
for (j=0; j<NY; j++){
rde[i*NY + j].gradx = (rde[(i+1)*NY + j].depth - rde[(i-1)*NY + j].depth)/dx;
}
}
/* left boundary */
for (j=0; j<NY; j++){
switch (B_COND_LEFT) {
case (BC_PERIODIC):
{
rde[j].gradx = (rde[NY + j].depth - rde[(NX-1)*NY + j].depth)/dx;
break;
}
case (BC_DIRICHLET):
{
rde[j].gradx = (rde[NY + j].depth - rde[j].depth)*2.0/dx;
break;
}
}
}
/* right boundary */
for (j=0; j<NY; j++){
switch (B_COND_RIGHT) {
case (BC_PERIODIC):
{
rde[(NX-1)*NY + j].gradx = (rde[j].depth - rde[(NX-2)*NY + j].depth)/dx;
break;
}
case (BC_DIRICHLET):
{
rde[(NX-1)*NY + j].gradx = (rde[(NX-1)*NY + j].depth - rde[(NX-2)*NY + j].depth)*2.0/dx;
break;
}
}
}
/* compute y gradient */
#pragma omp parallel for private(i,j)
for (i=0; i<NX; i++){
for (j=1; j<NY-1; j++){
rde[i*NY + j].grady = (rde[i*NY + j+1].depth - rde[i*NY + j-1].depth)/dy;
}
}
/* bottom boundary */
for (i=0; i<NX; i++){
switch (B_COND_BOTTOM) {
case (BC_PERIODIC):
{
rde[i*NY].grady = (rde[i*NY + 1].depth - rde[i*NY + NY-1].depth)/dy;
break;
}
case (BC_DIRICHLET):
{
rde[i*NY].grady = (rde[i*NY + 1].depth - rde[i*NY].depth)*2.0/dy;
break;
}
}
}
/* top boundary */
for (i=0; i<NX; i++){
switch (B_COND_TOP) {
case (BC_PERIODIC):
{
rde[i*NY + NY-1].grady = (rde[i*NY].depth - rde[i*NY + NY-2].depth)/dy;
break;
}
case (BC_DIRICHLET):
{
rde[i*NY + NY-1].grady = (rde[i*NY + NY-1].depth - rde[i*NY + NY-2].depth)*2.0/dy;
break;
}
}
}
/* compute light angle */
#pragma omp parallel for private(i,j)
for (i=0; i<NX; i++){
for (j=1; j<NY-1; j++){
norm = sqrt(vz*vz + rde[i*NY + j].gradx*rde[i*NY + j].gradx + rde[i*NY + j].grady*rde[i*NY + j].grady);
pscal = -rde[i*NY + j].gradx*light[0] - rde[i*NY + j].grady*light[1] + vz;
rde[i*NY+j].cos_depth_angle = pscal/norm;
}
}
/* for monitoring */
min = 0.0;
max = 0.0;
#pragma omp parallel for private(i,j)
for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
if (rde[i*NY + j].gradx > max) max = rde[i*NY + j].gradx;
if (rde[i*NY + j].gradx < min) min = rde[i*NY + j].gradx;
}
}
printf("gradx min = %.3lg, max = %.3lg\n", min, max);
min = 0.0;
max = 0.0;
#pragma omp parallel for private(i,j)
for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
if (rde[i*NY + j].grady > max) max = rde[i*NY + j].grady;
if (rde[i*NY + j].grady < min) min = rde[i*NY + j].grady;
}
}
printf("grady min = %.3lg, max = %.3lg\n", min, max);
min = 0.0;
max = 0.0;
#pragma omp parallel for private(i,j)
for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
if (rde[i*NY + j].depth > max) max = rde[i*NY + j].depth;
if (rde[i*NY + j].depth < min) min = rde[i*NY + j].depth;
}
}
printf("Depth min = %.3lg, max = %.3lg\n", min, max);
// for (j=0; j<NY; j++)
// {
// for (i=0; i<2; i++){
// printf("point (%03i,%03i): depth %.03lg, gradient (%.03lg, %.03lg)\n", i, j, rde[i*NY + j].depth, rde[i*NY + j].gradx, rde[i*NY + j].grady);
// }
// for (i=NX-2; i<NX; i++){
// printf("point (%03i,%03i): depth %.03lg, gradient (%.03lg, %.03lg)\n", i, j, rde[i*NY + j].depth, rde[i*NY + j].gradx, rde[i*NY + j].grady);
// }
// }
return(max);
}
void evolve_wave_half(double *phi_in[NFIELDS], double *phi_out[NFIELDS], short int xy_in[NX*NY], double potential_field[NX*NY], double vector_potential_field[2*NX*NY],
double gfield[2*NX*NY], t_rde rde[NX*NY])
/* time step of field evolution */
{
int i, j, k, iplus, iminus, jplus, jminus, ropening, w;
double x, y, z, deltax, deltay, deltaz, rho, rhox, rhoy, pot, u, v, ux, uy, vx, vy, test = 0.0, dx, dy, xy[2], padding, a;
double *delta_phi[NLAPLACIANS], *nabla_phi, *nabla_psi, *nabla_omega, *delta_vorticity, *delta_pressure, *delta_p, *delta_u, *delta_v, *nabla_rho, *nabla_u, *nabla_v;
double x, y, z, deltax, deltay, deltaz, rho, rhox, rhoy, pot, u, v, ux, uy, vx, vy, test = 0.0, dx, dy, xy[2], padding, a, eta, etax, etay;
double *delta_phi[NLAPLACIANS], *nabla_phi, *nabla_psi, *nabla_omega, *delta_vorticity, *delta_pressure, *delta_p, *delta_u, *delta_v, *nabla_rho, *nabla_u, *nabla_v, *nabla_eta;
// double u_bc[NY], v_bc[NY];
static double invsqr3 = 0.577350269; /* 1/sqrt(3) */
static double stiffness = 2.0; /* stiffness of Poisson equation solver */
@ -668,7 +928,7 @@ void evolve_wave_half(double *phi_in[NFIELDS], double *phi_out[NFIELDS], short i
y = YMIN + 0.02 + dy*((double)ropening);
x = YMAX - padding + MAZE_XSHIFT;
xy_to_pos(x, y, xy);
if ((B_DOMAIN == D_MAZE_CHANNELS)||(OBSTACLE_GEOMETRY == D_MAZE_CHANNELS))
if ((B_DOMAIN == D_MAZE_CHANNELS)||(OBSTACLE_GEOMETRY == D_MAZE_CHANNELS)||(OBSTACLE_GEOMETRY == D_MAZE_CHANNELS_INT))
{
imax = xy[0] + 2;
x = YMIN + padding + MAZE_XSHIFT;
@ -706,64 +966,80 @@ void evolve_wave_half(double *phi_in[NFIELDS], double *phi_out[NFIELDS], short i
compute_gradient_xy(phi_in[1], nabla_psi);
}
/* compute gradients of stream function and vorticity for Euler equation */
if (RDE_EQUATION == E_EULER_INCOMP)
{
nabla_psi = (double *)malloc(2*NX*NY*sizeof(double));
nabla_omega = (double *)malloc(2*NX*NY*sizeof(double));
compute_gradient_euler(phi_in[0], nabla_psi, EULER_GRADIENT_YSHIFT);
compute_gradient_euler(phi_in[1], nabla_omega, 0.0);
if (COMPUTE_PRESSURE) compute_pressure_laplacian(phi_in, delta_p);
dx = (XMAX-XMIN)/((double)NX);
dy = (YMAX-YMIN)/((double)NY);
if (SMOOTHEN_VORTICITY) /* beta: try to reduce formation of ripples */
switch (RDE_EQUATION){
/* compute gradients of stream function and vorticity for Euler equation */
case (E_EULER_INCOMP):
{
if (smooth == 0)
nabla_psi = (double *)malloc(2*NX*NY*sizeof(double));
nabla_omega = (double *)malloc(2*NX*NY*sizeof(double));
compute_gradient_euler(phi_in[0], nabla_psi, EULER_GRADIENT_YSHIFT);
compute_gradient_euler(phi_in[1], nabla_omega, 0.0);
if (COMPUTE_PRESSURE) compute_pressure_laplacian(phi_in, delta_p);
dx = (XMAX-XMIN)/((double)NX);
dy = (YMAX-YMIN)/((double)NY);
if (SMOOTHEN_VORTICITY) /* beta: try to reduce formation of ripples */
{
delta_vorticity = (double *)malloc(NX*NY*sizeof(double));
compute_laplacian_rde(phi_in[1], delta_vorticity, xy_in);
if (smooth == 0)
{
delta_vorticity = (double *)malloc(NX*NY*sizeof(double));
compute_laplacian_rde(phi_in[1], delta_vorticity, xy_in);
// #pragma omp parallel for private(i,delta_vorticity)
for (i=0; i<NX*NY; i++) phi_in[1][i] += intstep*SMOOTH_FACTOR*delta_vorticity[i];
free(delta_vorticity);
for (i=0; i<NX*NY; i++) phi_in[1][i] += intstep*SMOOTH_FACTOR*delta_vorticity[i];
free(delta_vorticity);
}
smooth++;
if (smooth >= SMOOTHEN_PERIOD) smooth = 0;
}
smooth++;
if (smooth >= SMOOTHEN_PERIOD) smooth = 0;
break;
}
}
/* compute gradients of fields for compressible Euler equation */
else if (RDE_EQUATION == E_EULER_COMP)
{
nabla_rho = (double *)malloc(2*NX*NY*sizeof(double));
// nabla_u = (double *)malloc(2*NX*NY*sizeof(double));
// nabla_v = (double *)malloc(2*NX*NY*sizeof(double));
compute_gradient_euler_test(phi_in[0], nabla_rho, xy_in);
compute_velocity_gradients(phi_in, rde);
// compute_gradient_euler_test(phi_in[1], nabla_u, xy_in);
// compute_gradient_euler_test(phi_in[2], nabla_v, xy_in);
if (SMOOTHEN_VELOCITY) /* beta: try to reduce formation of ripples */
/* compute gradients of fields for compressible Euler equation */
case (E_EULER_COMP):
{
if (smooth == 0)
nabla_rho = (double *)malloc(2*NX*NY*sizeof(double));
compute_gradient_euler_test(phi_in[0], nabla_rho, xy_in);
compute_velocity_gradients(phi_in, rde, xy_in);
if (SMOOTHEN_VELOCITY) /* beta: try to reduce formation of ripples */
{
if (smooth == 0)
{
delta_u = (double *)malloc(NX*NY*sizeof(double));
delta_v = (double *)malloc(NX*NY*sizeof(double));
compute_laplacian_rde(phi_in[1], delta_u, xy_in);
compute_laplacian_rde(phi_in[2], delta_v, xy_in);
#pragma omp parallel for private(i)
for (i=0; i<NX*NY; i++) phi_in[1][i] += intstep*SMOOTH_FACTOR*delta_u[i];
#pragma omp parallel for private(i)
for (i=0; i<NX*NY; i++) phi_in[2][i] += intstep*SMOOTH_FACTOR*delta_v[i];
free(delta_u);
free(delta_v);
}
smooth++;
if (smooth >= SMOOTHEN_PERIOD) smooth = 0;
}
break;
}
case (E_SHALLOW_WATER):
{
nabla_eta = (double *)malloc(2*NX*NY*sizeof(double));
compute_gradient_euler_test(phi_in[0], nabla_eta, xy_in);
compute_velocity_gradients(phi_in, rde, xy_in);
if (VISCOSITY > 0.0)
{
delta_u = (double *)malloc(NX*NY*sizeof(double));
delta_v = (double *)malloc(NX*NY*sizeof(double));
compute_laplacian_rde(phi_in[1], delta_u, xy_in);
compute_laplacian_rde(phi_in[2], delta_v, xy_in);
#pragma omp parallel for private(i)
for (i=0; i<NX*NY; i++) phi_in[1][i] += intstep*SMOOTH_FACTOR*delta_u[i];
#pragma omp parallel for private(i)
for (i=0; i<NX*NY; i++) phi_in[2][i] += intstep*SMOOTH_FACTOR*delta_v[i];
free(delta_u);
free(delta_v);
}
smooth++;
if (smooth >= SMOOTHEN_PERIOD) smooth = 0;
break;
}
}
}
if (TEST_GRADIENT) {
test = 0.0;
@ -899,6 +1175,48 @@ void evolve_wave_half(double *phi_in[NFIELDS], double *phi_out[NFIELDS], short i
}
break;
}
case (E_SHALLOW_WATER):
{
eta = phi_in[0][i*NY+j];
if (eta == 0.0) eta = 1.0e-6;
u = phi_in[1][i*NY+j];
v = phi_in[2][i*NY+j];
etax = nabla_eta[i*NY+j];
etay = nabla_eta[NX*NY+i*NY+j];
ux = rde[i*NY+j].dxu;
uy = rde[i*NY+j].dyu;
vx = rde[i*NY+j].dxv;
vy = rde[i*NY+j].dyv;
phi_out[0][i*NY+j] = eta - intstep*(u*etax + v*etay + eta*(ux + vy));
phi_out[1][i*NY+j] = u - intstep*(u*ux + v*uy + G_FIELD*etax);
phi_out[2][i*NY+j] = v - intstep*(u*vx + v*vy + G_FIELD*etay);
if (VISCOSITY > 0.0)
{
phi_out[1][i*NY+j] += intstep*VISCOSITY*delta_u[i*NY+j];
phi_out[2][i*NY+j] += intstep*VISCOSITY*delta_v[i*NY+j];
}
if (DISSIPATION > 0.0)
{
phi_out[1][i*NY+j] -= intstep*DISSIPATION*u;
phi_out[2][i*NY+j] -= intstep*DISSIPATION*v;
}
if (ADD_FORCE_FIELD)
{
phi_out[1][i*NY+j] += intstep*gfield[i*NY+j];
phi_out[2][i*NY+j] += intstep*gfield[NX*NY+i*NY+j];
}
if (VARIABLE_DEPTH)
{
phi_out[0][i*NY+j] -= intstep*rde[i*NY+j].depth*(ux + vy);
phi_out[0][i*NY+j] -= intstep*rde[i*NY+j].gradx*u;
phi_out[0][i*NY+j] -= intstep*rde[i*NY+j].grady*v;
}
break;
}
}
}
}
@ -942,6 +1260,15 @@ void evolve_wave_half(double *phi_in[NFIELDS], double *phi_out[NFIELDS], short i
{
free(nabla_rho);
}
else if (RDE_EQUATION == E_SHALLOW_WATER)
{
free(nabla_eta);
if (VISCOSITY > 0.0)
{
free(delta_u);
free(delta_v);
}
}
if (COMPUTE_PRESSURE)
{
@ -1142,7 +1469,7 @@ void print_parameters(double *phi[NFIELDS], t_rde rde[NX*NY], short int xy_in[NX
}
}
void draw_color_bar_palette(int plot, double range, int palette, int fade, double fade_value)
void draw_color_bar_palette(int plot, double range, int palette, int circular, int fade, double fade_value)
{
double width = 0.14;
// double width = 0.2;
@ -1151,7 +1478,8 @@ void draw_color_bar_palette(int plot, double range, int palette, int fade, doubl
{
if (ROTATE_COLOR_SCHEME)
draw_color_scheme_palette_3d(XMIN + 0.3, YMIN + 0.1, XMAX - 0.3, YMIN + 0.1 + width, plot, -range, range, palette, fade, fade_value);
// draw_color_scheme_palette_3d(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range, palette, fade, fade_value);
else if (circular)
draw_circular_color_scheme_palette_3d(XMAX - 2.0*width, YMAX - 2.0*width, 1.5*width, 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);
}
@ -1159,7 +1487,6 @@ void draw_color_bar_palette(int plot, double range, int palette, int fade, doubl
{
if (ROTATE_COLOR_SCHEME)
draw_color_scheme_palette_fade(XMIN + 0.8, YMIN + 0.1, XMAX - 0.8, YMIN + 0.1 + width, plot, -range, range, palette, fade, fade_value);
// draw_color_scheme_palette_fade(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range, palette, fade, fade_value);
else
draw_color_scheme_palette_fade(XMAX - 1.5*width, YMIN + 0.1, XMAX - 0.5*width, YMAX - 0.1, plot, -range, range, palette, fade, fade_value);
}
@ -1238,12 +1565,13 @@ void viewpoint_schedule(int i)
void animation()
{
double time = 0.0, scale, dx, var, jangle, cosj, sinj, sqrintstep,
intstep0, viscosity_printed, fade_value, noise = NOISE_INTENSITY;
intstep0, viscosity_printed, fade_value, noise = NOISE_INTENSITY, x, y, sign;
double *phi[NFIELDS], *phi_tmp[NFIELDS], *potential_field, *vector_potential_field, *tracers, *gfield, *bc_field, *bc_field2;
short int *xy_in;
int i, j, k, s, nvid, field;
static int counter = 0;
t_rde *rde;
// t_swater_depth *water_depth;
/* Since NX and NY are big, it seemed wiser to use some memory allocation here */
for (i=0; i<NFIELDS; i++)
@ -1274,14 +1602,27 @@ void animation()
initialize_vector_potential(vector_potential_field);
}
if (VARIABLE_DEPTH)
{
// water_depth = (t_swater_depth *)malloc(NX*NY*sizeof(t_swater_depth));
max_depth = initialize_water_depth(rde);
printf("Max depth = %.3lg\n", max_depth);
}
if (ADAPT_STATE_TO_BC)
{
bc_field = (double *)malloc(NX*NY*sizeof(double));
bc_field2 = (double *)malloc(NX*NY*sizeof(double));
initialize_bcfield(bc_field, bc_field2, polyrect);
}
if (ADD_FORCE_FIELD)
{
if (!ADAPT_STATE_TO_BC)
{
printf("Error: if ADD_FORCE_FIELD = 1, then ADAPT_STATE_TO_BC should be 1\n");
exit(1);
}
gfield = (double *)malloc(2*NX*NY*sizeof(double));
initialize_gfield(gfield, bc_field, bc_field2);
}
@ -1311,26 +1652,31 @@ void animation()
// init_fermion_state(-0.5, 0.5, 2.0, 0.0, 0.1, phi, xy_in);
// init_boson_state(-0.5, 0.5, 2.0, 0.0, 0.1, phi, xy_in);
// init_vortex_state(0.1, 0.4, 0.1, 0.3, -0.1, phi, xy_in);
// add_vortex_state(0.15, -0.4, -0.1, 0.4, 0.1, phi, xy_in);
// init_laminar_flow(flow_speed_schedule(0), LAMINAR_FLOW_MODULATION, LAMINAR_FLOW_YPERIOD, 0.0, phi, xy_in);
// init_laminar_flow(IN_OUT_FLOW_AMP, LAMINAR_FLOW_MODULATION, 0.02, 0.1, 1.0, 0.0, 0.1, phi, xy_in);
init_laminar_flow(0.05, LAMINAR_FLOW_MODULATION, 0.015, 0.1, 1.0, 0.0, 0.1, phi, xy_in);
// add_vortex_state(0.2, -0.4, -0.1, 0.3, -2.0, phi, xy_in);
// add_vortex_state(0.2, -1.0, -0.1, 0.3, -0.5, phi, xy_in);
// init_vortex_state(0.1, 1.0, 0.1, 0.3, -0.4, phi, xy_in);
// add_vortex_state(0.15, -1.0, -0.1, 0.4, 0.5, phi, xy_in);
add_vortex_state(0.2, 0.75, 0.1, 0.3, -0.5, phi, xy_in);
add_vortex_state(-0.35, -0.75, -0.1, 0.4, 0.5, phi, xy_in);
add_vortex_state(0.1, -0.3, 0.7, 0.1, -0.5, phi, xy_in);
// init_laminar_flow(0.05, LAMINAR_FLOW_MODULATION, 0.015, 0.1, 1.0, 0.0, 0.1, phi, xy_in);
// for (i=0; i<5; i++)
// for (j=0; j<3; j++)
// {
// x = XMIN + ((double)i - 0.5)*(XMAX-XMIN)/4.0 + 0.15*gaussian();
// y = YMIN + ((double)j - 0.5)*(YMAX-YMIN)/2.0 + 0.15*gaussian();
// sign = (double)((i+j)%2)*2.0 - 1.0;
// add_vortex_state(0.1*sign, x, y, 0.1, -0.35*sign + 0.1*gaussian(), phi, xy_in);
// }
// add_vortex_state(0.2, 0.75, 0.1, 0.3, -0.5, phi, xy_in);
// add_vortex_state(-0.35, -0.75, -0.1, 0.4, 0.5, phi, xy_in);
// add_vortex_state(0.1, -0.3, 0.7, 0.1, -0.5, phi, xy_in);
// init_pressure_gradient_flow(flow_speed_schedule(0), 1.0 + PRESSURE_GRADIENT, 1.0 - PRESSURE_GRADIENT, phi, xy_in, bc_field);
// init_shear_flow(-1.0, 0.1, 0.2, 1, 1, 0.2, phi, xy_in);
// init_shear_flow(1.0, 0.02, 0.15, 1, 1, 0.0, phi, xy_in);
// init_gaussian_wave(-1.0, 0.0, 0.005, 0.25, SWATER_MIN_HEIGHT, phi, xy_in);
init_linear_wave(-1.0, 0.01, 5.0e-8, 0.25, SWATER_MIN_HEIGHT, phi, xy_in);
// init_linear_blob(0.0, -0.75, 0.025, 0.0, 0.001, 0.25, 0.5, SWATER_MIN_HEIGHT, phi, xy_in);
// add_gaussian_wave(-1.6, -0.5, 0.015, 0.25, SWATER_MIN_HEIGHT, phi, xy_in);
if (ADAPT_STATE_TO_BC) adapt_state_to_bc(phi, bc_field, xy_in);
@ -1356,12 +1702,12 @@ void animation()
glutSwapBuffers();
printf("Drawing wave\n");
draw_wave_rde(0, phi, xy_in, rde, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1);
// draw_billiard();
if (PRINT_PARAMETERS) print_parameters(phi, rde, xy_in, time, PRINT_LEFT, VISCOSITY, noise);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 0, 1.0);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, 0, 1.0);
glutSwapBuffers();
sleep(SLEEP1);
@ -1425,7 +1771,7 @@ void animation()
if (ANTISYMMETRIZE_WAVE_FCT) antisymmetrize_wave_function(phi, xy_in);
for (j=0; j<NFIELDS; j++) printf("field[%i] = %.3lg\t", j, phi[j][0]);
for (j=0; j<NFIELDS; j++) printf("field[%i] = %.3lg\t", j, phi[j][NX*NY/2]);
printf("\n");
if (ADD_NOISE == 1)
@ -1460,7 +1806,7 @@ void animation()
// draw_billiard();
if (PRINT_PARAMETERS) print_parameters(phi, rde, xy_in, time, PRINT_LEFT, viscosity_printed, noise);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 0, 1.0);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, 0, 1.0);
// print_level(MDEPTH);
// print_Julia_parameters(i);
@ -1486,7 +1832,7 @@ void animation()
if (ADD_TRACERS) draw_tracers(phi, tracers, i, 0, 1.0);
// draw_billiard();
if (PRINT_PARAMETERS) print_parameters(phi, rde, xy_in, time, PRINT_LEFT, viscosity_printed, noise);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, CIRC_COLORBAR_B, 0, 1.0);
glutSwapBuffers();
// if (NO_EXTRA_BUFFER_SWAP) glutSwapBuffers();
save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter);
@ -1518,7 +1864,7 @@ void animation()
if (ADD_TRACERS) draw_tracers(phi, tracers, NSTEPS, 0, 1.0);
// draw_billiard();
if (PRINT_PARAMETERS) print_parameters(phi, rde, xy_in, time, PRINT_LEFT, viscosity_printed, noise);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 0, 1.0);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, 0, 1.0);
// if (!NO_EXTRA_BUFFER_SWAP) glutSwapBuffers();
glutSwapBuffers();
@ -1530,14 +1876,14 @@ void animation()
if (ADD_TRACERS) draw_tracers(phi, tracers, NSTEPS, 1, fade_value);
// draw_billiard();
if (PRINT_PARAMETERS) print_parameters(phi, rde, xy_in, time, PRINT_LEFT, viscosity_printed, noise);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 1, fade_value);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, 1, fade_value);
if (!NO_EXTRA_BUFFER_SWAP) glutSwapBuffers();
save_frame_counter(NSTEPS + i + 1);
}
draw_wave_rde(1, phi, xy_in, rde, potential_field, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 0, 1.0, REFRESH_B);
if (ADD_TRACERS) draw_tracers(phi, tracers, NSTEPS, 0, 1.0);
if (PRINT_PARAMETERS) print_parameters(phi, rde, xy_in, time, PRINT_LEFT, viscosity_printed, noise);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, CIRC_COLORBAR_B, 0, 1.0);
glutSwapBuffers();
if (!FADE) for (i=0; i<END_FRAMES; i++) save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter + i);
@ -1547,7 +1893,7 @@ void animation()
draw_wave_rde(1, phi, xy_in, rde, potential_field, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 1, fade_value, 0);
if (ADD_TRACERS) draw_tracers(phi, tracers, NSTEPS, 1, fade_value);
if (PRINT_PARAMETERS) print_parameters(phi, rde, xy_in, time, PRINT_LEFT, viscosity_printed, noise);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 1, fade_value);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, CIRC_COLORBAR_B, 1, fade_value);
glutSwapBuffers();
save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter + i);
}
@ -1560,7 +1906,7 @@ void animation()
fade_value = 1.0 - (double)i/(double)END_FRAMES;
draw_wave_rde(0, phi, xy_in, rde, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 1, fade_value, 0);
if (ADD_TRACERS) draw_tracers(phi, tracers, NSTEPS, 1, fade_value);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 1, fade_value);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, 1, fade_value);
glutSwapBuffers();
save_frame_counter(NSTEPS + 1 + counter + i);
}
@ -1581,6 +1927,7 @@ void animation()
free(potential_field);
free(vector_potential_field);
}
// if (VARIABLE_DEPTH) free(water_depth);
if (ADD_TRACERS) free(tracers);
if (ADD_FORCE_FIELD) free(gfield);
if (ADAPT_STATE_TO_BC)

1
sub_hashgrid.c
View File

@ -22,6 +22,7 @@ int bc_grouped(int bc)
case (BC_GENUS_TWO): return(4);
case (BC_ABSORBING): return(0);
case (BC_REFLECT_ABS): return(0);
case (BC_REFLECT_ABS_BOTTOM): return(0);
default:
{
printf("Warning: Hashgrid will not be properly initialised, update bc_grouped()\n\n");

472
sub_lj.c
View File

@ -481,6 +481,9 @@ double type_hue(int type)
if ((RD_REACTION == CHEM_CATALYTIC_A2D)&&(type == 4)) return(HUE_TYPE3);
if ((RD_REACTION == CHEM_ABDACBE)&&(type == 4)) return(HUE_TYPE3);
if ((RD_REACTION == CHEM_ABDACBE)&&(type == 5)) return(280.0);
switch (type) {
case (0): return(HUE_TYPE0);
case (1): return(HUE_TYPE0);
@ -898,6 +901,122 @@ void init_particle_config(t_particle particles[NMAXCIRCLES])
}
}
void add_particle_config(t_particle particles[NMAXCIRCLES], double xmin, double xmax, double ymin, double ymax, double radius)
/* add particles to configuration */
{
int i, j, k, n, ncirc0, n_p_active, ncandidates = PDISC_CANDIDATES, naccepted, newcircles;
double dx, dy, p, phi, r, r0, ra[5], sa[5], height, x, y = 0.0, gamma, dpoisson, xx[4], yy[4];
short int active_poisson[NMAXCIRCLES], far;
dpoisson = PDISC_DISTANCE*radius;
switch (CIRCLE_PATTERN_B) {
case (C_SQUARE):
{
ncircles = NGRIDX*NGRIDY;
dy = (YMAX - YMIN)/((double)NGRIDY);
for (i = 0; i < NGRIDX; i++)
for (j = 0; j < NGRIDY; j++)
{
n = NGRIDY*i + j;
particles[n].xc = ((double)(i-NGRIDX/2) + 0.5)*dy;
particles[n].yc = YMIN + ((double)j + 0.5)*dy;
particles[n].radius = MU;
particles[n].active = 1;
}
break;
}
case (C_HEX):
{
dx = (xmax - xmin)/((double)NGRIDX);
dy = (ymax - ymin)/((double)NGRIDY);
// dx = dy*0.5*sqrt(3.0);
for (i = 0; i < NGRIDX; i++)
for (j = 0; j < NGRIDY+1; j++)
{
n = ncircles + (NGRIDY+1)*i + j;
// particles[n].xc = ((double)(i-NGRIDX/2) + 0.5)*dx; /* is +0.5 needed? */
particles[n].xc = xmin + ((double)i - 0.5)*dx;
particles[n].yc = ymin + ((double)j - 0.5)*dy;
if ((i+NGRIDX)%2 == 1) particles[n].yc += 0.5*dy;
particles[n].radius = radius;
/* activate only circles that intersect the domain */
if ((particles[n].yc < ymax + radius)&&(particles[n].yc > ymin - radius)&&(particles[n].xc < xmax + radius)&&(particles[n].xc > xmin - radius)) particles[n].active = 1;
else particles[n].active = 0;
}
ncircles += NGRIDX*(NGRIDY+1);
break;
}
case (C_POISSON_DISC):
{
ncirc0 = ncircles;
printf("Generating new Poisson disc sample\n");
/* generate first circle */
particles[ncirc0].xc = (xmax - xmin)*(double)rand()/RAND_MAX + xmin;
particles[ncirc0].yc = (ymax - ymin)*(double)rand()/RAND_MAX + ymin;
active_poisson[0] = 1;
// // particles[0].active = 1;
n_p_active = 1;
newcircles = 1;
while ((n_p_active > 0)&&(ncircles < NMAXCIRCLES))
{
/* randomly select an active circle */
i = rand()%(newcircles);
while (!active_poisson[i]) i = rand()%(ncircles);
// printf("Starting from circle %i at (%.3f,%.3f)\n", i, particles[i].xc, particles[i].yc);
/* generate new candidates */
naccepted = 0;
for (j=0; j<ncandidates; j++)
{
r = dpoisson*(2.0*(double)rand()/RAND_MAX + 1.0);
phi = DPI*(double)rand()/RAND_MAX;
x = particles[ncirc0 + i].xc + r*cos(phi);
y = particles[ncirc0 + i].yc + r*sin(phi);
// printf("Testing new circle at (%.3f,%.3f)\t", x, y);
far = 1;
for (k=0; k<ncircles; k++) if ((k!=i))
{
/* new circle is far away from circle k */
far = far*((x - particles[k].xc)*(x - particles[k].xc) + (y - particles[k].yc)*(y - particles[k].yc) >= dpoisson*dpoisson);
/* new circle is in domain */
far = far*(x < xmax)*(x > xmin)*(y < ymax)*(y > ymin);
// far = far*(vabs(x) < LAMBDA)*(y < INITYMAX)*(y > INITYMIN);
}
if (far) /* accept new circle */
{
printf("New circle at (%.3f,%.3f) accepted\n", x, y);
particles[ncircles].xc = x;
particles[ncircles].yc = y;
particles[ncircles].radius = radius;
particles[ncircles].active = 1;
active_poisson[ncircles] = 1;
ncircles++;
ncirc0++;
n_p_active++;
naccepted++;
}
// else printf("Rejected\n");
}
if (naccepted == 0) /* inactivate circle i */
{
// printf("No candidates work, inactivate circle %i\n", i);
active_poisson[i] = 0;
n_p_active--;
}
printf("%i active circles\n", n_p_active);
}
printf("Generated %i circles\n", ncircles);
break;
}
default:
{
printf("Function init_circle_config not defined for this pattern \n");
}
}
}
void init_people_config(t_person people[NMAXCIRCLES])
/* initialise particle configuration */
{
@ -1021,6 +1140,17 @@ void init_obstacle_config(t_obstacle obstacle[NMAXOBSTACLES])
add_obstacle(0.5*(double)i, YMIN + 0.3, radius, obstacle);
break;
}
case (O_CIRCLE):
{
n = 0;
obstacle[n].xc = 0.0;
obstacle[n].yc = 0.0;
obstacle[n].radius = OBSTACLE_RADIUS;
obstacle[n].active = 1;
nobstacles = 1;
break;
}
default:
{
printf("Function init_obstacle_config not defined for this pattern \n");
@ -1162,7 +1292,7 @@ void add_circle_to_segments(double x, double y, double r, int nsegs, double angl
double nozzle_width(double x, double width, int nozzle_shape)
/* width of bell-shaped nozzle */
{
double lam = 0.5*LAMBDA;
double lam = 0.5*LAMBDA, a, b;
if (x >= 0.0) return(width);
else switch (nozzle_shape) {
@ -1176,6 +1306,16 @@ double nozzle_width(double x, double width, int nozzle_shape)
if (-x < 0.1) return(width - (0.5 - width)*x/0.1);
else return(0.5);
}
case (NZ_DELAVAL):
{
a = 1.5;
b = 0.05;
return(sqrt(width*width - 0.5*x) + a*b*x*(1.0 + x)/(b + x*x));
// a = (sqrt(width*width+0.5) - width)/sqrt(0.5);
// c = (a*sqrt(0.5*h) - width)/(h*h);
// if (-x < h) return(width + c*x*x);
// else return(width + a * sqrt(-0.5*x));
}
default: return(0.0);
}
}
@ -1185,7 +1325,7 @@ void add_rocket_to_segments(t_segment segment[NMAXSEGMENTS], double x0, double y
/* add one or several rocket_shaped set of segments */
{
int i, j, cycle = 0, nsegments0;
double angle, dx, x1, y1, x2, y2, nozx, nozy, a, b;
double angle, dx, x1, y1, x2, y2, nozx, nozy, a, b, c, w;
nsegments0 = nsegments;
@ -1234,6 +1374,25 @@ void add_rocket_to_segments(t_segment segment[NMAXSEGMENTS], double x0, double y
add_rotated_angle_to_segments(x0-a, nozy, x0-nozx, nozy, 0.02, segment, 0);
break;
}
case (RCK_RECT_BAR): /* rectangular chamber with a hat and separating bar */
{
a = 0.5*LAMBDA;
b = (0.49*PI-0.25)*LAMBDA;
c = 0.5*a;
w = 0.025;
add_rotated_angle_to_segments(x0+nozx, nozy, x0+nozx, nozy+0.5*c, 0.02, segment, 0);
add_rotated_angle_to_segments(x0+nozx, nozy+0.5*c, x0+nozx+0.5*c, nozy+c, 0.02, segment, 0);
add_rotated_angle_to_segments(x0+nozx+0.5*c, nozy+c, x0+a, nozy+c, 0.02, segment, 0);
add_rotated_angle_to_segments(x0+a, nozy+c, x0+a, nozy+b+c, 0.02, segment, 0);
add_rotated_angle_to_segments(x0+a, nozy+b+c, x0, nozy+b+a+c, 0.02, segment, 0);
add_rotated_angle_to_segments(x0-w, nozy+a+c, x0-w, nozy+b+a+c, 2.0*w, segment, 0);
add_rotated_angle_to_segments(x0, nozy+b+a+c, x0-a, nozy+b+c, 0.02, segment, 0);
add_rotated_angle_to_segments(x0-a, nozy+b+c, x0-a, nozy+c, 0.02, segment, 0);
add_rotated_angle_to_segments(x0-a, nozy+c, x0-nozx-0.5*c, nozy+c, 0.02, segment, 0);
add_rotated_angle_to_segments(x0-nozx-0.5*c, nozy+c, x0-nozx+w, nozy+0.5*c, 0.02, segment, 0);
add_rotated_angle_to_segments(x0-nozx+w, nozy+0.5*c, x0-nozx+w, nozy, 0.02, segment, 0);
break;
}
}
dx = LAMBDA/(double)(nsides);
@ -1275,7 +1434,7 @@ void add_rocket_to_segments(t_segment segment[NMAXSEGMENTS], double x0, double y
for (i=nsegments0; i<nsegments; i++) segment[i].group = group;
}
int init_maze_segments(t_segment segment[NMAXSEGMENTS])
int init_maze_segments(t_segment segment[NMAXSEGMENTS], int diag)
/* init segments forming a maze */
{
t_maze* maze;
@ -1297,7 +1456,11 @@ int init_maze_segments(t_segment segment[NMAXSEGMENTS])
x1 = YMIN + padding + (double)i*dx + MAZE_XSHIFT;
y1 = YMIN + padding + (double)j*dy;
if (((i>0)||(j!=NYMAZE/2))&&(maze[n].west)) add_rectangle_to_segments(x1, y1, x1 - width, y1 + dy, segment, 0);
if (diag)
{
if (((i>0)||(j<NYMAZE-1))&&(maze[n].west)) add_rectangle_to_segments(x1, y1, x1 - width, y1 + dy, segment, 0);
}
else if (((i>0)||(j!=NYMAZE/2))&&(maze[n].west)) add_rectangle_to_segments(x1, y1, x1 - width, y1 + dy, segment, 0);
if (maze[n].south) add_rectangle_to_segments(x1, y1, x1 + dx, y1 - width, segment, 0);
}
@ -1305,15 +1468,29 @@ int init_maze_segments(t_segment segment[NMAXSEGMENTS])
add_rectangle_to_segments(YMIN + padding + MAZE_XSHIFT, YMAX - padding, YMAX - padding + MAZE_XSHIFT, YMAX - padding - width, segment, 0);
/* right side of maze */
y1 = YMIN + padding + dy*((double)NYMAZE/2);
x1 = YMAX - padding + MAZE_XSHIFT;
add_rectangle_to_segments(x1, YMIN - 1.0, x1 - width, y1 - dy, segment, 0);
add_rectangle_to_segments(x1, y1, x1 - width, YMAX + 1.0, segment, 0);
if (diag)
{
y1 = YMIN + padding + dy;
add_rectangle_to_segments(x1, y1, x1 - width, YMAX + 10.0, segment, 0);
}
else
{
y1 = YMIN + padding + dy*((double)NYMAZE/2);
add_rectangle_to_segments(x1, YMIN - 1.0, x1 - width, y1 - dy, segment, 0);
add_rectangle_to_segments(x1, y1, x1 - width, YMAX + 1.0, segment, 0);
}
/* left side of maze */
x1 = YMIN + padding + MAZE_XSHIFT;
add_rectangle_to_segments(x1, YMIN - 1.0, x1 - width, YMIN + padding, segment, 0);
add_rectangle_to_segments(x1, YMAX - padding, x1 - width, YMAX + 1.0, segment, 0);
add_rectangle_to_segments(x1, YMAX - padding, x1 - width, YMAX + 10.0, segment, 0);
if (diag)
{
add_rotated_angle_to_segments(XMIN, YMAX - 0.5*dy, x1, YMAX - dy - 2.0*width, width, segment, 0);
add_rectangle_to_segments(XMIN, YMAX - 0.5*dy, XMIN - width, YMAX + 10.0, segment, 0);
}
free(maze);
}
@ -1822,7 +1999,20 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
}
case (S_MAZE):
{
init_maze_segments(segment);
init_maze_segments(segment, 0);
cycle = 0;
for (i=0; i<nsegments; i++)
{
segment[i].group = 0;
segment[i].inactivate = 0;
}
break;
}
case (S_MAZE_DIAG):
{
init_maze_segments(segment, 1);
cycle = 0;
for (i=0; i<nsegments; i++)
@ -2096,7 +2286,7 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
int in_rocket(double x, double y, int rocket_shape)
/* returns 1 if (x,y) is in rocket chamber, with translated coordinates */
{
double l, y1, a, b;
double l, y1, a, b, c;
switch (rocket_shape) {
case (RCK_DISC) :
@ -2125,6 +2315,18 @@ int in_rocket(double x, double y, int rocket_shape)
return(y1 < a + b - 0.05 - vabs(x));
// return(1);
}
case (RCK_RECT_BAR) :
{
a = 0.5*LAMBDA;
b = (0.49*PI-0.25)*LAMBDA;
c = 0.5*a;
y1 = y - YMIN - LAMBDA;
if (vabs(x) > 0.95*a) return(0);
if (vabs(x) < 0.1*a) return(0);
if (y1 < 0.05 + c) return(0);
if (y1 < b - 0.05 + c) return(1);
return(y1 < a + b - 0.05 + c - vabs(x));
}
}
}
@ -2293,6 +2495,14 @@ int in_segment_region(double x, double y, t_segment segment[NMAXSEGMENTS])
}
return(1);
}
case (S_MAZE_DIAG):
{
for (i=0; i<nsegments; i++)
{
if (distance_to_segment(x, y, segment[i].x1, segment[i].y1, segment[i].x2, segment[i].y2) < 5.0*MAZE_WIDTH) return(0);
}
return(1);
}
default: return(1);
}
}
@ -2911,7 +3121,7 @@ int compute_particle_interaction(int i, int k, double force[2], double *torque,
/* compute repelling force and torque of particle #k on particle #i */
/* returns 1 if distance between particles is smaller than NBH_DIST_FACTOR*MU */
{
double x1, y1, x2, y2, r, f, angle, aniso, fx, fy, ff[2], dist_scaled, spin_f, ck, sk, ck_rel, sk_rel;
double x1, y1, x2, y2, r, f, angle, aniso, fx, fy, ff[2], dist_scaled, spin_f, ck, sk, ck_rel, sk_rel, alpha, amp;
static double dxhalf = 0.5*(BCXMAX - BCXMIN), dyhalf = 0.5*(BCYMAX - BCYMIN);
int wwrapx, wwrapy, twrapx, twrapy;
@ -3010,6 +3220,37 @@ int compute_particle_interaction(int i, int k, double force[2], double *torque,
force[1] = f*KSPRING_VICSEK*(particle[k].vy - particle[i].vy);
break;
}
case (I_VICSEK_SHARK):
{
if (particle[k].type == particle[i].type)
{
f = cos(0.5*(particle[k].angle - particle[i].angle));
force[0] = f*KSPRING_VICSEK*(particle[k].vx - particle[i].vx);
force[1] = f*KSPRING_VICSEK*(particle[k].vy - particle[i].vy);
}
else if (particle[i].type != 2)
{
f = krepel*coulomb_force(distance, particle[k]);
force[0] = f*ca;
force[1] = f*sa;
}
else
{
if (VICSEK_REPULSION > 0.0)
{
// f = VICSEK_REPULSION*harmonic_force(distance, particle[k]);
f = VICSEK_REPULSION*coulomb_force(distance, particle[k]);
force[0] = f*ca;
force[1] = f*sa;
}
else
{
force[0] = 0.0;
force[1] = 0.0;
}
}
break;
}
}
if (ROTATION)
@ -3053,6 +3294,27 @@ int compute_particle_interaction(int i, int k, double force[2], double *torque,
else *torque = sin(particle[k].angle - particle[i].angle);
break;
}
case (I_VICSEK_SHARK):
{
if (dist_scaled > 10.0) *torque = 0.0;
else if (particle[k].type == particle[i].type) /* fish adjusting direction */
{
if (twrapx||twrapy) *torque = sin(-particle[k].angle - particle[i].angle);
else *torque = sin(particle[k].angle - particle[i].angle);
}
else if (particle[k].type == 2) /* fish fleeing a shark */
{
alpha = argument(ca,sa);
particle[i].angle = alpha + PI;
*torque = 0.0;
}
else /* shark tracking fish */
{
*torque = cos(particle[k].angle)*sa - sin(particle[k].angle)*ca;
}
break;
}
default:
{
spin_f = particle[i].spin_freq;
@ -3266,6 +3528,9 @@ int add_particle(double x, double y, double vx, double vy, double mass, short in
// particle[i].spin_freq = SPIN_INTER_FREQUENCY_B;
// }
if ((PLOT == P_NUMBER)||(PLOT_B == P_NUMBER))
particle[i].color_hue = 360.0*(double)(i%N_PARTICLE_COLORS)/(double)N_PARTICLE_COLORS;
ncircles++;
printf("Added particle at (%.3lg, %.3lg)\n", x, y);
@ -3442,6 +3707,13 @@ void compute_particle_colors(t_particle particle, int plot, double rgb[3], doubl
*width = BOUNDARY_WIDTH;
break;
}
case (P_NUMBER):
{
hue = particle.color_hue;
*radius = particle.radius;
*width = BOUNDARY_WIDTH;
break;
}
}
switch (plot) {
@ -3466,6 +3738,13 @@ void compute_particle_colors(t_particle particle, int plot, double rgb[3], doubl
hsl_to_rgb_twilight(hue, 0.9, 0.5, rgby);
break;
}
case (P_ANGLE):
{
hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_ANGLE);
hsl_to_rgb_palette(hue, 0.9, 0.5, rgbx, COLOR_PALETTE_ANGLE);
hsl_to_rgb_palette(hue, 0.9, 0.5, rgby, COLOR_PALETTE_ANGLE);
break;
}
case (P_DIRECT_ENERGY):
{
hsl_to_rgb_twilight(hue, 0.9, 0.5, rgb);
@ -3486,6 +3765,13 @@ void compute_particle_colors(t_particle particle, int plot, double rgb[3], doubl
hsl_to_rgb_twilight(hue, 0.9, 0.5, rgby);
break;
}
case (P_INITIAL_POS):
{
hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_INITIAL_POS);
hsl_to_rgb_palette(hue, 0.9, 0.5, rgbx, COLOR_PALETTE_INITIAL_POS);
hsl_to_rgb_palette(hue, 0.9, 0.5, rgby, COLOR_PALETTE_INITIAL_POS);
break;
}
default:
{
hsl_to_rgb(hue, 0.9, 0.5, rgb);
@ -3988,6 +4274,29 @@ int draw_special_particle(t_particle particle, double xc1, double yc1, double ra
}
break;
}
case (CHEM_ABDACBE):
{
if (particle.type == 4)
{
for (wsign = -1; wsign <= 1; wsign+=2)
{
x1 = xc1 + (double)wsign*0.7*radius*cos(particle.angle);
y1 = yc1 + (double)wsign*0.7*radius*sin(particle.angle);
if (wsign == 1)
{
if (fill) draw_colored_polygon(x1, y1, 1.2*MU_B, nsides, angle + APOLY*PID, rgb);
else draw_polygon(x1, y1, 1.2*MU_B, nsides, angle + APOLY*PID);
}
else
{
if (fill) draw_colored_polygon(x1, y1, 1.2*MU, nsides, angle + APOLY*PID, rgb);
else draw_polygon(x1, y1, 1.2*MU, nsides, angle + APOLY*PID);
}
}
return(0);
}
break;
}
}
return(1);
}
@ -4011,7 +4320,7 @@ void draw_one_particle(t_particle particle, double xc, double yc, double radius,
/* specific shapes for chemical reactions */
if (REACTION_DIFFUSION) cont = draw_special_particle(particle, xc1, yc1, radius, angle, nsides, rgb, 1);
if ((particle.interaction == I_LJ_QUADRUPOLE)||(particle.interaction == I_LJ_DIPOLE)||(particle.interaction == I_VICSEK)||(particle.interaction == I_VICSEK_REPULSIVE)||(particle.interaction == I_VICSEK_SPEED))
if ((particle.interaction == I_LJ_QUADRUPOLE)||(particle.interaction == I_LJ_DIPOLE)||(particle.interaction == I_VICSEK)||(particle.interaction == I_VICSEK_REPULSIVE)||(particle.interaction == I_VICSEK_SPEED)||(particle.interaction == I_VICSEK_SHARK))
draw_colored_rhombus(xc1, yc1, radius, angle + APOLY*PID, rgb);
else if (cont) draw_colored_polygon(xc1, yc1, radius, nsides, angle + APOLY*PID, rgb);
@ -4041,7 +4350,7 @@ void draw_one_particle(t_particle particle, double xc, double yc, double radius,
glColor3f(1.0, 1.0, 1.0);
if (REACTION_DIFFUSION) cont = draw_special_particle(particle, xc1, yc1, radius, angle, nsides, rgb, 0);
if ((particle.interaction == I_LJ_QUADRUPOLE)||(particle.interaction == I_LJ_DIPOLE)||(particle.interaction == I_VICSEK)||(particle.interaction == I_VICSEK_REPULSIVE)||(particle.interaction == I_VICSEK_SPEED))
if ((particle.interaction == I_LJ_QUADRUPOLE)||(particle.interaction == I_LJ_DIPOLE)||(particle.interaction == I_VICSEK)||(particle.interaction == I_VICSEK_REPULSIVE)||(particle.interaction == I_VICSEK_SPEED)||(particle.interaction == I_VICSEK_SHARK))
draw_rhombus(xc1, yc1, radius, angle + APOLY*PID);
else if (cont) draw_polygon(xc1, yc1, radius, nsides, angle + APOLY*PID);
@ -5646,6 +5955,30 @@ double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacl
return(fperp);
}
case (BC_REFLECT_ABS_BOTTOM):
{
/* add harmonic force outside screen */
padding = MU;
x = particle[j].xc;
y = particle[j].yc;
if (y < BCYMIN)
{
particle[j].active = 0;
particle[j].vx = 0.0;
particle[j].vy = 0.0;
particle[j].xc = BCXMAX + 2.0*padding;
particle[j].yc = BCYMIN - 2.0*padding;
}
else
{
if (y > BCYMAX + padding) particle[j].fy -= KSPRING_BOUNDARY*(y - BCYMAX - padding);
if (x > BCXMAX - padding) particle[j].fx -= KSPRING_BOUNDARY*(x - BCXMAX + padding);
else if (x < BCXMIN + padding) particle[j].fx += KSPRING_BOUNDARY*(BCXMIN - x + padding);
}
return(0.0);
}
}
}
@ -5730,9 +6063,8 @@ int reorder_particles(t_particle particle[NMAXCIRCLES], double py[NMAXCIRCLES],
}
int initialize_configuration(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HASHX*HASHY],
t_obstacle obstacle[NMAXOBSTACLES], double px[NMAXCIRCLES], double py[NMAXCIRCLES], double pangle[NMAXCIRCLES], int tracer_n[N_TRACER_PARTICLES],
t_segment segment[NMAXSEGMENTS])
int initialize_configuration(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HASHX*HASHY], t_obstacle obstacle[NMAXOBSTACLES], double px[NMAXCIRCLES], double py[NMAXCIRCLES], double pangle[NMAXCIRCLES], int tracer_n[N_TRACER_PARTICLES],
t_segment segment[NMAXSEGMENTS])
/* initialize all particles, obstacles, and the hashgrid */
{
int i, j, k, n, nactive = 0;
@ -5747,6 +6079,11 @@ int initialize_configuration(t_particle particle[NMAXCIRCLES], t_hashgrid hashgr
particle[i].type = 2;
particle[i].radius = MU_B;
}
if ((INTERACTION == I_VICSEK_SHARK)&&(i==1))
{
particle[i].type = 2;
particle[i].radius = MU_B;
}
particle[i].neighb = 0;
particle[i].diff_neighb = 0;
@ -5781,6 +6118,13 @@ int initialize_configuration(t_particle particle[NMAXCIRCLES], t_hashgrid hashgr
particle[i].vx = V_INITIAL*gaussian();
particle[i].vy = V_INITIAL*gaussian();
if ((INTERACTION == I_VICSEK_SHARK)&&(i==1))
{
particle[i].vx *= 1000.0;
particle[i].vy *= 1000.0;
}
particle[i].energy = (particle[i].vx*particle[i].vx + particle[i].vy*particle[i].vy)*particle[i].mass_inv;
px[i] = particle[i].vx;
@ -5802,6 +6146,8 @@ int initialize_configuration(t_particle particle[NMAXCIRCLES], t_hashgrid hashgr
if ((PLOT == P_INITIAL_POS)||(PLOT_B == P_INITIAL_POS))
particle[i].color_hue = 360.0*(particle[i].yc - INITYMIN)/(INITYMAX - INITYMIN);
else if ((PLOT == P_NUMBER)||(PLOT_B == P_NUMBER))
particle[i].color_hue = 360.0*(double)(i%N_PARTICLE_COLORS)/(double)N_PARTICLE_COLORS;
}
/* initialize dummy values in case particles are added */
for (i=ncircles; i < NMAXCIRCLES; i++)
@ -5892,7 +6238,7 @@ int initialize_configuration(t_particle particle[NMAXCIRCLES], t_hashgrid hashgr
/* position-dependent particle type */
if (POSITION_DEPENDENT_TYPE) for (i=0; i<ncircles; i++)
if (((!POSITION_Y_DEPENDENCE)&&(particle[i].xc < 0))||((POSITION_Y_DEPENDENCE)&&(particle[i].yc < 0)))
if (((!POSITION_Y_DEPENDENCE)&&((particle[i].xc - POSITION_DEP_X)*POSITION_DEP_SIGN < 0.0))||((POSITION_Y_DEPENDENCE)&&(particle[i].yc*POSITION_DEP_SIGN < 0.0)))
{
particle[i].type = 2;
particle[i].mass_inv = 1.0/PARTICLE_MASS_B;
@ -6074,7 +6420,8 @@ int initialize_configuration(t_particle particle[NMAXCIRCLES], t_hashgrid hashgr
particle[i].mass_inv = 1.0/PARTICLE_MASS_B;
}
break;
}case (IC_TWOROCKETS):
}
case (IC_TWOROCKETS):
{
if (vabs(particle[i].xc) < SEGMENTS_X0) particle[i].type = 1;
else
@ -6085,6 +6432,18 @@ int initialize_configuration(t_particle particle[NMAXCIRCLES], t_hashgrid hashgr
}
break;
}
case (IC_TWOROCKETS_TWOFUELS):
{
if (vabs(particle[i].xc) < SEGMENTS_X0) particle[i].type = 1;
else
{
if (particle[i].xc < 0) particle[i].type = 2;
else particle[i].type = 3;
particle[i].radius = MU_B;
particle[i].mass_inv = 1.0/PARTICLE_MASS_B;
}
break;
}
}
}
@ -6136,11 +6495,13 @@ int add_particles(t_particle particle[NMAXCIRCLES], double px[NMAXCIRCLES], doub
// x = INITXMIN + (INITXMAX - INITXMIN)*(double)rand()/(double)RAND_MAX;
// y = INITYMIN + (INITYMAX - INITYMIN)*(double)rand()/(double)RAND_MAX;
// x = BCXMIN + (BCXMAX - BCXMIN)*(double)rand()/(double)RAND_MAX;
// y = YMAX + 0.5*(BCYMAX - YMAX)*(double)rand()/(double)RAND_MAX;
printf("Adding a particle\n\n");
x = BCXMIN + (BCXMAX - BCXMIN)*(double)rand()/(double)RAND_MAX;
y = YMAX + 0.5*(BCYMAX - YMAX)*(double)rand()/(double)RAND_MAX;
x = ADDXMIN + (ADDXMAX - ADDXMIN)*(double)rand()/(double)RAND_MAX;
y = ADDYMIN + 0.5*(ADDYMAX - ADDYMIN)*(double)rand()/(double)RAND_MAX;
add_particle(x, y, 0.0, 0.0, PARTICLE_MASS, 0, particle);
// add_particle(BCXMIN + 0.1, 0.5*(BCYMIN + BCYMAX), 200.0, 0.0, PARTICLE_MASS, 0, particle);
// i++;
@ -6381,6 +6742,13 @@ void compute_inverse_masses(double inv_masses[RD_TYPES+1])
inv_masses[6] = inv_masses[1];
break;
}
case (CHEM_ABDACBE):
{
inv_masses[3] = inv_masses[1];
inv_masses[4] = 1.0/(PARTICLE_MASS + PARTICLE_MASS_B);
inv_masses[5] = inv_masses[1];
break;
}
}
}
@ -6459,6 +6827,13 @@ void compute_radii(double radii[RD_TYPES+1])
radii[6] = MU;
break;
}
case (CHEM_ABDACBE):
{
radii[3] = MU;
radii[4] = 1.2*MU_B;
radii[5] = MU;
break;
}
}
}
@ -6855,6 +7230,12 @@ int chem_transfer(int i, int type2, int newtype1, int newtype2, t_particle parti
particle[k].vy *= newmass_inv2/mass_inv2;
particle[k].mass_inv = newmass_inv2;
if (EXOTHERMIC)
{
adapt_speed_exothermic(&particle[i], DELTA_EKIN);
adapt_speed_exothermic(&particle[k], DELTA_EKIN);
}
collisions[ncollisions].x = 0.5*(particle[i].xc + particle[k].xc);
collisions[ncollisions].y = 0.5*(particle[i].yc + particle[k].yc);
collisions[ncollisions].time = COLLISION_TIME;
@ -6996,7 +7377,6 @@ int update_types(t_particle particle[NMAXCIRCLES], t_collision *collisions, int
printf("%i collisions\n", ncollisions);
delta_n = ncollisions - oldncollisions;
printf("delta_n = %i\n", delta_n);
// if (delta_n > 1) delta_n = 1;
if (EXOTHERMIC) *delta_e = (double)(delta_n)*DELTA_EKIN;
return(ncollisions);
}
@ -7240,11 +7620,30 @@ int update_types(t_particle particle[NMAXCIRCLES], t_collision *collisions, int
}
}
return(ncollisions);
}
}
case (CHEM_ABDACBE):
{
for (i=0; i<ncircles; i++)
{
// oldncollisions = ncollisions;
if ((particle[i].active)&&(particle[i].type == 1))
{
ncollisions = chem_merge(i, 2, 4, particle, collisions, ncollisions, inv_masses, radii);
// if (EXOTHERMIC) *delta_e += (double)(ncollisions - oldncollisions)*DELTA_EKIN;
ncollisions = chem_transfer(i, 3, 2, 5, particle, collisions, ncollisions, inv_masses, radii, REACTION_DIST);
// if (EXOTHERMIC) *delta_e += (double)(ncollisions - oldncollisions)*DELTA_EKIN;
}
}
printf("%i collisions\n", ncollisions);
delta_n = ncollisions - oldncollisions;
printf("delta_n = %i\n", delta_n);
if (EXOTHERMIC) *delta_e = (double)(delta_n)*DELTA_EKIN;
return(ncollisions);
}
}
}
double plot_coord(double x, double xmin, double xmax)
{
return(xmin + x*(xmax - xmin));
@ -7598,3 +7997,30 @@ void init_segment_group(t_segment segment[NMAXSEGMENTS], int group, t_group_segm
printf("Segment group data %i: (%.3lg, %.3lg)\n", group, segment_group[group].xc, segment_group[group].yc);
}
void reset_energy(t_particle particle[NMAXCIRCLES], double px[NMAXCIRCLES], double py[NMAXCIRCLES], double totalenergy, double emean)
/* decrease energy in case of blow-up */
{
int i;
double vratio, emax;
emax = 10.0*emean/(double)ncircles;
// printf("Warning: blow-up, resetting energy from %.5lg to %.5lg\n\n", totalenergy, emean);
printf("Warning: blow-up, resetting energy of some particles\n");
// vratio = sqrt(emean/totalenergy);
for (i=0; i < ncircles; i++) if (particle[i].energy > emax)
{
printf("Particle %i at (%.3lg, %.3lg) has energy %.5lg, resetting to 0\n", i, particle[i].xc, particle[i].yc, particle[i].energy);
particle[i].vx = 0.0;
particle[i].vy = 0.0;
px[i] = 0.0;
py[i] = 0.0;
particle[i].energy = 0.0;
}
printf("\n\n");
}

1350
sub_rde.c
View File

File diff suppressed because it is too large Load Diff

323
sub_wave.c
View File

@ -2000,6 +2000,93 @@ int compute_circular_maze_coordinates(t_rect_rotated polyrectrot[NMAXPOLY], t_ar
free(maze);
}
int compute_interior_maze_coordinates(t_rectangle polyrect[NMAXPOLY], int type)
/* compute positions of complement of maze */
{
t_maze* maze;
int i, j, n, nsides = 0, ropening;
double dx, dy, x1, y1, x0, padding = 0.02, pos[2], width = MAZE_WIDTH;
/* TODO */
maze = (t_maze *)malloc(NXMAZE*NYMAZE*sizeof(t_maze));
ropening = (NYMAZE+1)/2;
init_maze(maze);
/* move the entrance for maze type with two channels */
if (type == 2)
{
n = nmaze(0, ropening-1);
maze[n].west = 0;
n = nmaze(0, ropening);
maze[n].west = 1;
}
/* build walls of maze */
// x0 = LAMBDA - 1.0;
dx = (YMAX - YMIN - 2.0*padding)/(double)(NXMAZE);
dy = (YMAX - YMIN - 2.0*padding)/(double)(NYMAZE);
for (i=0; i<NXMAZE; i++)
for (j=0; j<NYMAZE; j++)
{
n = nmaze(i, j);
x1 = YMIN + padding + (double)i*dx + MAZE_XSHIFT;
y1 = YMIN + padding + (double)j*dy;
if (!maze[n].west)
{
polyrect[nsides].x1 = x1 - 0.5*dx - width;
polyrect[nsides].y1 = y1 + 0.5*dx - width;
polyrect[nsides].x2 = x1 + 0.5*dx + width;
polyrect[nsides].y2 = y1 + 0.5*dx + width;
nsides++;
}
if (!maze[n].south)
{
polyrect[nsides].x1 = x1 + 0.5*dx - width;
polyrect[nsides].y1 = y1 - 0.5*dx - width;
polyrect[nsides].x2 = x1 + 0.5*dx + width;
polyrect[nsides].y2 = y1 + 0.5*dx + width;
nsides++;
}
}
/* right channel of maze */
y1 = YMIN + padding + dy*((double)ropening);
x1 = YMAX - padding + MAZE_XSHIFT;
polyrect[nsides].x1 = x1 - 0.5*dx - width;
polyrect[nsides].y1 = y1 - 0.5*dy - width;
polyrect[nsides].x2 = XMAX;
polyrect[nsides].y2 = y1 - 0.5*dy + width;
nsides++;
/* left channel of maze */
x1 = YMIN + padding + MAZE_XSHIFT;
polyrect[nsides].x1 = XMIN;
polyrect[nsides].y1 = y1 - 0.5*dy - width;
polyrect[nsides].x2 = x1 + 0.5*dx + width;
polyrect[nsides].y2 = y1 - 0.5*dy + width;
nsides++;
for (i=0; i<nsides; i++)
{
xy_to_pos(polyrect[i].x1, polyrect[i].y1, pos);
polyrect[i].posi1 = pos[0];
polyrect[i].posj1 = pos[1];
xy_to_pos(polyrect[i].x2, polyrect[i].y2, pos);
polyrect[i].posi2 = pos[0];
polyrect[i].posj2 = pos[1];
}
free(maze);
return(nsides);
}
int init_polyline(int depth, t_vertex polyline[NMAXPOLY])
/* initialise variable polyline, for certain polygonal domain shapes */
{
@ -2103,6 +2190,10 @@ int init_polyrect_euler(t_rectangle polyrect[NMAXPOLY], int domain)
{
return(compute_maze_coordinates(polyrect, 2));
}
case (D_MAZE_CHANNELS_INT):
{
return(compute_interior_maze_coordinates(polyrect, 2));
}
}
}
@ -5055,6 +5146,158 @@ void draw_color_scheme_palette_fade(double x1, double y1, double x2, double y2,
}
void draw_circular_color_scheme_palette_fade(double x1, double y1, double radius, int plot, double min, double max, int palette, int fade, double fade_value)
{
int j, k;
double x, y, dy, dy_e, dy_phase, rgb[3], value, lum, amp, dphi, pos[2], phi, xy[2];
glBegin(GL_TRIANGLE_FAN);
xy_to_pos(x1, y1, xy);
glVertex2d(xy[0], xy[1]);
dy = (max - min)/360.0;
dy_e = max/360.0;
dy_phase = 1.0/360.0;
dphi = DPI/360.0;
for (j = 0; j <= 360; j++)
{
switch (plot) {
case (P_AMPLITUDE):
{
value = min + 1.0*dy*(double)(j);
color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
break;
}
case (P_ENERGY):
{
value = dy_e*(double)(j)*100.0/E_SCALE;
if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
break;
}
case (P_MEAN_ENERGY):
{
value = dy_e*(double)(j)*100.0/E_SCALE;
if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
break;
}
case (P_LOG_ENERGY):
{
value = LOG_SCALE*log(dy_e*(double)(j)*100.0/E_SCALE);
color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
break;
}
case (P_LOG_MEAN_ENERGY):
{
value = LOG_SHIFT + LOG_SCALE*log(dy_e*(double)(j)*100.0/E_SCALE);
color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
break;
}
case (P_ENERGY_FLUX):
{
value = dy_e*(double)(j);
if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
break;
}
case (P_TOTAL_ENERGY_FLUX):
{
value = dy_phase*(double)(j);
color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb);
break;
}
case (P_PHASE):
{
value = min + 1.0*dy*(double)(j);
amp = 0.5*color_amplitude_linear(value, 1.0, 1);
while (amp > 1.0) amp -= 2.0;
while (amp < -1.0) amp += 2.0;
amp_to_rgb(0.5*(1.0 + amp), rgb);
break;
}
case (Z_EULER_VORTICITY):
{
value = min + 1.0*dy*(double)(j);
color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb);
break;
}
case (Z_EULER_LOG_VORTICITY):
{
value = min + 1.0*dy*(double)(j);
color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb);
break;
}
case (Z_EULER_VORTICITY_ASYM):
{
value = min + 1.0*dy*(double)(j);
color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb);
break;
}
case (Z_EULER_LPRESSURE):
{
value = min + 1.0*dy*(double)(j);
color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb);
break;
}
case (Z_EULER_PRESSURE):
{
value = min + 1.0*dy*(double)(j);
color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb);
break;
}
case (Z_EULER_DENSITY):
{
value = min + 1.0*dy*(double)(j);
color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb);
break;
}
case (Z_EULER_SPEED):
{
value = min + 1.0*dy*(double)(j);
color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb);
break;
}
case (Z_EULERC_VORTICITY):
{
value = min + 1.0*dy*(double)(j);
color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb);
break;
}
default:
{
value = min + 1.0*dy*(double)(j);
color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb);
break;
}
}
if (fade) for (k=0; k<3; k++) rgb[k] *= fade_value;
glColor3f(rgb[0], rgb[1], rgb[2]);
xy_to_pos(x1 + radius*cos(dphi*(double)j), y1 + radius*sin(dphi*(double)j), xy);
glVertex2d(xy[0], xy[1]);
}
xy_to_pos(x1 + radius*cos(dphi), y1 + radius*sin(dphi), xy);
glVertex2d(xy[0], xy[1]);
glEnd ();
if (fade) glColor3f(fade_value, fade_value, fade_value);
else glColor3f(1.0, 1.0, 1.0);
glLineWidth(BOUNDARY_WIDTH*3/2);
glEnable(GL_LINE_SMOOTH);
dphi = DPI/(double)NSEG;
glBegin(GL_LINE_LOOP);
for (j = 0; j < NSEG; j++)
{
xy_to_pos(x1 + radius*cos(dphi*(double)j), y1 + radius*sin(dphi*(double)j), xy);
glVertex2d(xy[0], xy[1]);
glVertex2d(xy[0], xy[1]);
}
glEnd ();
}
void print_speed(double speed, int fade, double fade_value)
{
char message[100];
@ -5550,12 +5793,13 @@ void init_ior_2d(short int *xy_in[NX], double *tcc_table[NX], double ior_angle)
/* compute variable index of refraction */
/* should be at some point merged with 3D version in suv_wave_3d.c */
{
int i, j, k, n, inlens;
int i, j, k, n, inlens, ncircles;
double courant2 = COURANT*COURANT, courantb2 = COURANTB*COURANTB, lambda1, mu1;
double u, v, u1, x, y, xy[2], norm2, speed, r2, c, salpha, h, ll, ca, sa, x1, y1, dx, dy, sum, sigma, x0, y0, rgb[3];
double xc[NGRIDX*NGRIDY], yc[NGRIDX*NGRIDY], height[NGRIDX*NGRIDY];
static double xc_stat[NGRIDX*NGRIDY], yc_stat[NGRIDX*NGRIDY], sigma_stat;
static int first = 1;
t_circle circles[NMAXCIRCLES];
rgb[0] = 1.0;
rgb[1] = 1.0;
@ -5749,7 +5993,9 @@ void init_ior_2d(short int *xy_in[NX], double *tcc_table[NX], double ior_angle)
}
// #pragma omp parallel for private(i,j)
for (i=0; i<NX; i++){
for (i=0; i<NX; i++)
{
if (i%100 == 0) printf("Computing potential for column %i of %i\n", i, NX);
for (j=0; j<NY; j++){
ij_to_xy(i, j, xy);
x = xy[0];
@ -5857,7 +6103,78 @@ void init_ior_2d(short int *xy_in[NX], double *tcc_table[NX], double ior_angle)
}
break;
}
default:
case (IOR_POISSON_WELLS):
{
ncircles = init_circle_config_pattern(circles, C_POISSON_DISC);
for (n = 0; n<ncircles; n++)
{
height[n] = 0.5 + 0.5*gaussian();
if (height[n] > 1.0) height[n] = 1.0;
if (height[n] < 0.0) height[n] = 0.0;
}
for (n = 0; n<ncircles; n++) printf("Circle %i at (%.3lg, %.3lg) height %.3lg\n", n, circles[n].xc, circles[n].yc, height[n]);
sigma = 0.2*(XMAX - XMIN)*(YMAX - YMIN)/(double)ncircles;
// sigma = MU*MU;
for (i=0; i<NX; i++)
{
if (i%100 == 0) printf("Computing potential for column %i of %i\n", i, NX);
for (j=0; j<NY; j++){
ij_to_xy(i, j, xy);
x = xy[0];
y = xy[1];
sum = 0.0;
for (n = 0; n<ncircles; n++)
{
r2 = (x - circles[n].xc)*(x - circles[n].xc) + (y - circles[n].yc)*(y - circles[n].yc);
sum += exp(-r2/(sigma))*height[n];
}
sum = tanh(sum);
// printf("%.3lg\n", sum);
tcc_table[i][j] = COURANT*sum + COURANTB*(1.0-sum);
}
}
break;
}
case (IOR_PPP_WELLS):
{
ncircles = init_circle_config_pattern(circles, C_RAND_POISSON);
for (n = 0; n<ncircles; n++)
{
height[n] = 0.5 + 0.5*gaussian();
if (height[n] > 1.0) height[n] = 1.0;
if (height[n] < 0.0) height[n] = 0.0;
}
for (n = 0; n<ncircles; n++) printf("Circle %i at (%.3lg, %.3lg) height %.3lg\n", n, circles[n].xc, circles[n].yc, height[n]);
sigma = 0.2*(XMAX - XMIN)*(YMAX - YMIN)/(double)ncircles;
// sigma = MU*MU;
for (i=0; i<NX; i++)
{
if (i%100 == 0) printf("Computing potential for column %i of %i\n", i, NX);
for (j=0; j<NY; j++){
ij_to_xy(i, j, xy);
x = xy[0];
y = xy[1];
sum = 0.0;
for (n = 0; n<ncircles; n++)
{
r2 = (x - circles[n].xc)*(x - circles[n].xc) + (y - circles[n].yc)*(y - circles[n].yc);
sum += exp(-r2/(sigma))*height[n];
}
sum = tanh(sum);
// printf("%.3lg\n", sum);
tcc_table[i][j] = COURANT*sum + COURANTB*(1.0-sum);
}
}
break;
}default:
{
for (i=0; i<NX; i++){
for (j=0; j<NY; j++){

78
sub_wave_3d.c
View File

@ -1296,6 +1296,7 @@ void init_speed_dissipation(short int xy_in[NX*NY], double tc[NX*NY], double tcc
double courant2 = COURANT*COURANT, courantb2 = COURANTB*COURANTB, lambda1, mu1;
double u, v, u1, x, y, xy[2], norm2, speed, r2, c, salpha, h, ll, ca, sa, x1, y1, dx, dy, sum, sigma, x0, y0, rgb[3];
double xc[NGRIDX*NGRIDY], yc[NGRIDX*NGRIDY], height[NGRIDX*NGRIDY];
t_circle circles[NMAXCIRCLES];
if (VARIABLE_IOR)
{
@ -1541,6 +1542,82 @@ void init_speed_dissipation(short int xy_in[NX*NY], double tc[NX*NY], double tcc
}
break;
}
case (IOR_POISSON_WELLS):
{
ncircles = init_circle_config_pattern(circles, C_POISSON_DISC);
for (n = 0; n<ncircles; n++)
{
height[n] = 0.5 + 0.5*gaussian();
if (height[n] > 1.0) height[n] = 1.0;
if (height[n] < 0.0) height[n] = 0.0;
}
for (n = 0; n<ncircles; n++) printf("Circle %i at (%.3lg, %.3lg) height %.3lg\n", n, circles[n].xc, circles[n].yc, height[n]);
sigma = 0.2*(XMAX - XMIN)*(YMAX - YMIN)/(double)ncircles;
// sigma = MU*MU;
for (i=0; i<NX; i++)
{
if (i%100 == 0) printf("Computing potential for column %i of %i\n", i, NX);
for (j=0; j<NY; j++){
ij_to_xy(i, j, xy);
x = xy[0];
y = xy[1];
sum = 0.0;
for (n = 0; n<ncircles; n++)
{
r2 = (x - circles[n].xc)*(x - circles[n].xc) + (y - circles[n].yc)*(y - circles[n].yc);
sum += exp(-r2/(sigma))*height[n];
}
sum = tanh(sum);
// printf("%.3lg\n", sum);
tc[i*NY+j] = COURANT*sum + COURANTB*(1.0-sum);
tcc[i*NY+j] = COURANT*sum + COURANTB*(1.0-sum);
tgamma[i*NY+j] = GAMMA;
}
}
break;
}
case (IOR_PPP_WELLS):
{
ncircles = init_circle_config_pattern(circles, C_RAND_POISSON);
for (n = 0; n<ncircles; n++)
{
height[n] = 0.5 + 0.5*gaussian();
if (height[n] > 1.0) height[n] = 1.0;
if (height[n] < 0.0) height[n] = 0.0;
}
for (n = 0; n<ncircles; n++) printf("Circle %i at (%.3lg, %.3lg) height %.3lg\n", n, circles[n].xc, circles[n].yc, height[n]);
sigma = 0.2*(XMAX - XMIN)*(YMAX - YMIN)/(double)ncircles;
// sigma = MU*MU;
for (i=0; i<NX; i++)
{
if (i%100 == 0) printf("Computing potential for column %i of %i\n", i, NX);
for (j=0; j<NY; j++){
ij_to_xy(i, j, xy);
x = xy[0];
y = xy[1];
sum = 0.0;
for (n = 0; n<ncircles; n++)
{
r2 = (x - circles[n].xc)*(x - circles[n].xc) + (y - circles[n].yc)*(y - circles[n].yc);
sum += exp(-r2/(sigma))*height[n];
}
sum = tanh(sum);
// printf("%.3lg\n", sum);
tc[i*NY+j] = COURANT*sum + COURANTB*(1.0-sum);
tcc[i*NY+j] = COURANT*sum + COURANTB*(1.0-sum);
tgamma[i*NY+j] = GAMMA;
}
}
break;
}
default:
{
for (i=0; i<NX; i++){
@ -2121,7 +2198,6 @@ void draw_color_scheme_palette_3d(double x1, double y1, double x2, double y2, in
draw_rectangle_noscale(x1, y1, x2, y2);
}
void print_speed_3d(double speed, int fade, double fade_value)
{
char message[100];

176
sub_wave_3d_rde.c
View File

@ -1529,3 +1529,179 @@ void draw_color_scheme_palette_3d(double x1, double y1, double x2, double y2, in
draw_rectangle_noscale(x1, y1, x2, y2);
}
void draw_circular_color_scheme_palette_3d(double x1, double y1, double radius, int plot,
double min, double max, int palette, int fade, double fade_value)
{
int j, k, ij_center[2], ij_right[2], ic, jc, ir;
double x, y, dy, dy_e, dy_phase, rgb[3], value, lum, amp, dphi, pos[2], phi, xy[2];
// printf("Drawing color bar\n");
xy_to_ij(x1, y1, ij_center);
xy_to_ij(x1 + radius, y1, ij_right);
// rgb[0] = 0.0; rgb[1] = 0.0; rgb[2] = 0.0;
// erase_area_rgb(0.5*(x1 + x2), x2 - x1, 0.5*(y1 + y2), y2 - y1, rgb);
ic = ij_center[0];
jc = ij_center[1];
ir = ij_right[0] - ij_center[0];
// imax = ij_topright[0];
// jmax = ij_topright[1];
glBegin(GL_TRIANGLE_FAN);
draw_vertex_ij(ic, jc);
dy = (max - min)/360.0;
dy_e = max/360.0;
dy_phase = 1.0/360.0;
dphi = DPI/360.0;
for (j = 0; j < 360; j++)
{
switch (plot) {
case (P_3D_AMPLITUDE):
{
value = min + 1.0*dy*(double)(j);
color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb);
break;
}
case (P_3D_ANGLE):
{
value = 1.0*dy*(double)(j);
color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 0, rgb);
break;
}
case (P_3D_AMP_ANGLE):
{
value = min + 1.0*dy*(double)(j);
color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb);
break;
}
case (P_3D_ENERGY):
{
value = dy_e*(double)(j)*100.0/E_SCALE;
if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
break;
}
case (P_3D_LOG_ENERGY):
{
value = LOG_SCALE*dy_e*(double)(j)*100.0/E_SCALE;
color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
break;
}
case (P_3D_TOTAL_ENERGY):
{
value = dy_e*(double)(j)*100.0/E_SCALE;
if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
break;
}
case (P_3D_LOG_TOTAL_ENERGY):
{
value = LOG_SCALE*dy_e*(double)(j)*100.0/E_SCALE;
color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
break;
}
case (P_3D_MEAN_ENERGY):
{
value = dy_e*(double)(j)*100.0/E_SCALE;
if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
break;
}
case (P_3D_LOG_MEAN_ENERGY):
{
value = LOG_SCALE*dy_e*(double)(j)*100.0/E_SCALE;
color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
break;
}
case (P_3D_PHASE):
{
value = dy_phase*(double)(j);
color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb);
break;
}
case (P_3D_FLUX_INTENSITY):
{
value = dy_e*(double)(j)*100.0/E_SCALE;
if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
break;
}
case (P_3D_FLUX_DIRECTION):
{
value = dy_phase*(double)(j);
color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb);
break;
}
case (Z_EULER_VORTICITY):
{
value = min + 1.0*dy*(double)(j);
color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb);
break;
}
case (Z_EULER_LOG_VORTICITY):
{
value = min + 1.0*dy*(double)(j);
color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb);
break;
}
case (Z_EULER_VORTICITY_ASYM):
{
value = min + 1.0*dy*(double)(j);
color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb);
break;
}
case (Z_EULER_LPRESSURE):
{
value = min + 1.0*dy*(double)(j);
color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb);
break;
}
case (Z_EULERC_VORTICITY):
{
value = min + 1.0*dy*(double)(j);
printf("Palette value %.3lg\n", value);
color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb);
break;
}
case (Z_SWATER_DIRECTION_SPEED):
{
value = dy_phase*(double)(j);
color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb);
break;
}
}
if (fade) for (k=0; k<3; k++) rgb[k] *= fade_value;
glColor3f(rgb[0], rgb[1], rgb[2]);
x = cos(dphi*(double)j)*(double)ir;
y = sin(dphi*(double)j)*(double)ir;
ij_to_xy(ic + x, jc + y, xy);
glVertex2d(xy[0], xy[1]);
}
draw_vertex_ij(ic + ir, jc);
glEnd ();
if (fade) glColor3f(fade_value, fade_value, fade_value);
else glColor3f(1.0, 1.0, 1.0);
glLineWidth(BOUNDARY_WIDTH*3/2);
glEnable(GL_LINE_SMOOTH);
dphi = DPI/(double)NSEG;
glBegin(GL_LINE_LOOP);
for (j = 0; j < NSEG; j++)
{
x = cos(dphi*(double)j)*(double)ir;
y = sin(dphi*(double)j)*(double)ir;
xy[0] = XMIN + ((double)ic + x)*(XMAX-XMIN)/((double)NX);
xy[1] = YMIN + ((double)jc + y)*(YMAX-YMIN)/((double)NY);
glVertex2d(xy[0], xy[1]);
}
glEnd ();
}

132
sub_wave_comp.c
View File

@ -483,6 +483,11 @@ int xy_in_billiard_half(double x, double y, int domain, int pattern, int top)
int i, j, k, k1, k2, condition, type;
switch (domain) {
case (D_NOTHING):
{
return(1);
break;
}
case (D_MENGER):
{
x1 = 0.5*(x+1.0);
@ -744,7 +749,12 @@ void draw_billiard_half(int domain, int pattern, int top) /* draws the bill
/* Do nothing */
break;
}
default:
case (D_NOTHING):
{
/* Do nothing */
break;
}
default:
{
printf("Function draw_billiard not defined for this billiard \n");
}
@ -799,6 +809,31 @@ void init_wave_flat_comp( double *phi[NX], double *psi[NX], short int * xy_in[NX
}
void add_circular_wave_comp(double factor, double x, double y, double *phi[NX], double *psi[NX], short int * xy_in[NX], int half)
/* add drop at (x,y) to the field with given prefactor */
{
int i, j;
double xy[2], dist2;
for (i=0; i<NX; i++)
{
if (half == 1) for (j=NY/2; j<NY; j++)
{
ij_to_xy(i, j, xy);
dist2 = (xy[0]-x)*(xy[0]-x) + (xy[1]-y)*(xy[1]-y);
if ((xy_in[i][j])||(TWOSPEEDS))
phi[i][j] += INITIAL_AMP*factor*exp(-dist2/INITIAL_VARIANCE)*cos(-sqrt(dist2)/INITIAL_WAVELENGTH);
}
else for (j=0; j<NY/2; j++)
{
ij_to_xy(i, j, xy);
dist2 = (xy[0]-x)*(xy[0]-x) + (xy[1]-y)*(xy[1]-y);
if ((xy_in[i][j])||(TWOSPEEDS))
phi[i][j] += INITIAL_AMP*factor*exp(-dist2/INITIAL_VARIANCE)*cos(-sqrt(dist2)/INITIAL_WAVELENGTH);
}
}
}
void draw_wave_comp(double *phi[NX], double *psi[NX], short int *xy_in[NX], double scale, int time, int plot)
/* draw the field */
{
@ -893,6 +928,101 @@ void draw_wave_comp(double *phi[NX], double *psi[NX], short int *xy_in[NX], doub
glEnd();
}
void draw_wave_comp_highres_palette(int size, double *phi[NX], double *psi[NX], double *total_energy[NX], short int *xy_in[NX], double scale, int time, int plot, int palette, int fade, double fade_value)
/* draw the field */
{
int i, j, iplus, iminus, jplus, jminus, k;
double rgb[3], xy[2], x1, y1, x2, y2, velocity, energy, gradientx2, gradienty2, pos[2];
static double dtinverse = ((double)NX)/(COURANT*(XMAX-XMIN)), dx = (XMAX-XMIN)/((double)NX);
glBegin(GL_QUADS);
// printf("dtinverse = %.5lg\n", dtinverse);
for (i=0; i<NX; i++)
for (j=0; j<NY; j++)
{
if (((TWOSPEEDS)&&(xy_in[i][j] != 2))||(xy_in[i][j] == 1))
{
switch (plot) {
case (P_AMPLITUDE):
{
color_scheme_palette(COLOR_SCHEME, palette, phi[i][j], scale, time, rgb);
break;
}
case (P_ENERGY):
{
energy = compute_energy(phi, psi, xy_in, i, j);
if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, energy, scale, time, rgb);
else color_scheme_palette(COLOR_SCHEME, palette, energy, scale, time, rgb);
break;
}
case (P_MIXED):
{
if (j > NY/2) color_scheme_palette(COLOR_SCHEME, palette, phi[i][j], scale, time, rgb);
else color_scheme_palette(COLOR_SCHEME, palette, compute_energy(phi, psi, xy_in, i, j), scale, time, rgb);
break;
}
case (P_MEAN_ENERGY):
{
energy = compute_energy(phi, psi, xy_in, i, j);
total_energy[i][j] += energy;
if (COLOR_PALETTE >= COL_TURBO)
color_scheme_asym_palette(COLOR_SCHEME, palette, total_energy[i][j]/(double)(time+1), scale, time, rgb);
else color_scheme_palette(COLOR_SCHEME, palette, total_energy[i][j]/(double)(time+1), scale, time, rgb);
break;
}
case (P_LOG_ENERGY):
{
energy = compute_energy(phi, psi, xy_in, i, j);
color_scheme(COLOR_SCHEME, LOG_SHIFT + LOG_SCALE*log(energy), scale, time, rgb);
break;
}
// case (P_LOG_MEAN_ENERGY):
// {
// energy = compute_energy(phi, psi, xy_in, i, j);
// if (energy == 0.0) energy = 1.0e-20;
// total_energy[i][j] += energy;
// color_scheme(COLOR_SCHEME, LOG_SCALE*log(total_energy[i][j]/(double)(time+1)), scale, time, rgb);
// break;
// }
}
// if (PLOT == P_AMPLITUDE)
// color_scheme(COLOR_SCHEME, phi[i][j], scale, time, rgb);
// else if (PLOT == P_ENERGY)
// {
// energy = compute_energy(phi, psi, xy_in, i, j);
// if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym(COLOR_SCHEME, energy, scale, time, rgb);
// else color_scheme(COLOR_SCHEME, energy, scale, time, rgb);
// }
// else if (PLOT == P_MIXED)
// {
// if (j > NY/2) color_scheme(COLOR_SCHEME, phi[i][j], scale, time, rgb);
// else color_scheme(COLOR_SCHEME, compute_energy(phi, psi, xy_in, i, j), scale, time, rgb);
// }
if (fade) for (k=0; k<3; k++) rgb[k] *= fade_value;
glColor3f(rgb[0], rgb[1], rgb[2]);
glVertex2i(i, j);
glVertex2i(i+size, j);
glVertex2i(i+size, j+size);
glVertex2i(i, j+size);
}
}
glEnd ();
/* draw horizontal mid line */
glColor3f(1.0, 1.0, 1.0);
glBegin(GL_LINE_STRIP);
xy_to_pos(XMIN, 0.5*(YMIN+YMAX), pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(XMAX, 0.5*(YMIN+YMAX), pos);
glVertex2d(pos[0], pos[1]);
glEnd();
}
void compute_energy_tblr(double *phi[NX], double *psi[NX], short int *xy_in[NX], double energies[6])
/* compute total energy in top/bottom left/right boxes */
{

163
wave_3d.c
View File

@ -54,21 +54,19 @@
#define WINHEIGHT 1150 /* window height */
// // // #define NX 2500 /* number of grid points on x axis */
// // // #define NY 1250 /* number of grid points on y axis */
#define NX 1920 /* number of grid points on x axis */
#define NY 1150 /* number of grid points on y axis */
// #define NX 2840 /* number of grid points on x axis */
// #define NY 2300 /* number of grid points on y axis */
// #define NX 2500 /* number of grid points on x axis */
// #define NY 1250 /* number of grid points on y axis */
#define NX 1700 /* number of grid points on x axis */
#define NY 1700 /* number of grid points on y axis */
// #define NX 1700 /* number of grid points on x axis */
// #define NY 1700 /* 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 XMIN -1.5
#define XMAX 1.5 /* x interval */
#define YMIN -1.5
#define YMAX 1.5 /* y interval for 9/16 aspect ratio */
#define XMIN -1.0
#define XMAX 3.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 */
@ -94,7 +92,7 @@
/* Choice of the billiard table */
#define B_DOMAIN 999 /* choice of domain shape, see list in global_pdes.c */
#define B_DOMAIN 10 /* choice of domain shape, see list in global_pdes.c */
#define CIRCLE_PATTERN 2 /* pattern of circles or polygons, see list in global_pdes.c */
@ -103,16 +101,16 @@
#define CIRCLE_PATTERN_B 0 /* second pattern of circles or polygons */
#define VARIABLE_IOR 0 /* set to 1 for a variable index of refraction */
#define IOR 101 /* choice of index of refraction, see list in global_pdes.c */
#define IOR 9 /* choice of index of refraction, see list in global_pdes.c */
#define IOR_TOTAL_TURNS 1.0 /* total angle of rotation for IOR_PERIODIC_WELLS_ROTATING */
#define MANDEL_IOR_SCALE -0.05 /* parameter controlling dependence 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 NPOISSON 1000 /* number of points for Poisson C_RAND_POISSON arrangement */
#define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */
#define LAMBDA 0.5 /* parameter controlling the dimensions of domain */
#define MU 0.5 /* parameter controlling the dimensions of domain */
#define LAMBDA 0.1197916667 /* parameter controlling the dimensions of domain */
#define MU 0.035 /* parameter controlling the dimensions of domain */
#define NPOLY 6 /* number of sides of polygon */
#define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */
#define MDEPTH 7 /* depth of computation of Menger gasket */
@ -120,8 +118,8 @@
#define MANDELLEVEL 2000 /* iteration level for Mandelbrot set */
#define MANDELLIMIT 20.0 /* limit value for approximation of Mandelbrot set */
#define FOCI 1 /* set to 1 to draw focal points of ellipse */
#define NGRIDX 8 /* number of grid point for grid of disks */
#define NGRIDY 8 /* number of grid point for grid of disks */
#define NGRIDX 30 /* number of grid point for grid of disks */
#define NGRIDY 18 /* number of grid point for grid of disks */
#define X_SHOOTER -0.2
#define Y_SHOOTER -0.6
@ -149,8 +147,8 @@
#define AMPLITUDE 0.8 /* amplitude of periodic excitation */
#define ACHIRP 0.2 /* acceleration coefficient in chirp */
#define DAMPING 0.0 /* damping of periodic excitation */
#define COURANT 0.06 /* Courant number */
#define COURANTB 0.06 /* Courant number in medium B */
#define COURANT 0.1 /* Courant number */
#define COURANTB 0.01 /* Courant number in medium B */
#define GAMMA 0.0 /* damping factor in wave equation */
#define GAMMAB 1.0e-6 /* damping factor in wave equation */
#define GAMMA_SIDES 1.0e-4 /* damping factor on boundary */
@ -164,7 +162,7 @@
/* For similar wave forms, COURANT^2*GAMMA should be kept constant */
#define ADD_OSCILLATING_SOURCE 1 /* set to 1 to add an oscillating wave source */
#define OSCILLATING_SOURCE_PERIOD 2 /* period of oscillating source */
#define OSCILLATING_SOURCE_PERIOD 30 /* period of oscillating source */
#define ALTERNATE_OSCILLATING_SOURCE 1 /* set to 1 to alternate sign of oscillating source */
#define ADD_WAVE_PACKET_SOURCES 0 /* set to 1 to add several sources emitting wave packets */
@ -181,9 +179,9 @@
/* Parameters for length and speed of simulation */
#define NSTEPS 2450 /* number of frames of movie */
// #define NSTEPS 500 /* number of frames of movie */
#define NVID 5 /* number of iterations between images displayed on screen */
#define NSTEPS 2300 /* number of frames of movie */
// #define NSTEPS 300 /* number of frames of movie */
#define NVID 4 /* number of iterations between images displayed on screen */
// #define NVID 10 /* 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 */
@ -201,18 +199,19 @@
/* Parameters of initial condition */
// #define INITIAL_AMP 0.02 /* amplitude of initial condition */
#define INITIAL_AMP 0.007 /* amplitude of initial condition */
#define INITIAL_VARIANCE 0.0003 /* variance of initial condition */
#define INITIAL_WAVELENGTH 0.015 /* wavelength of initial condition */
#define INITIAL_AMP 2.0 /* amplitude of initial condition */
// #define INITIAL_VARIANCE 0.000025 /* variance of initial condition */
#define INITIAL_VARIANCE 0.00005 /* variance of initial condition */
#define INITIAL_WAVELENGTH 0.05 /* wavelength of initial condition */
/* Plot type, see list in global_pdes.c */
#define ZPLOT 103 /* wave height */
#define CPLOT 103 /* color scheme */
#define ZPLOT_B 112
#define CPLOT_B 113 /* plot type for second movie */
#define ZPLOT_B 107
#define CPLOT_B 107 /* plot type for second movie */
#define CHANGE_LUMINOSITY 1 /* set to 1 to let luminosity depend on energy flux intensity */
#define FLUX_WINDOW 30 /* size of averaging window of flux intensity */
@ -242,8 +241,8 @@
/* Color schemes */
#define COLOR_PALETTE 15 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 10 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE 11 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 13 /* Color palette, see list in global_pdes.c */
#define BLACK 1 /* background */
@ -252,11 +251,11 @@
#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 2.0 /* additional scaling factor for color scheme P_3D_AMPLITUDE */
#define VSCALE_ENERGY 15.0 /* additional scaling factor for color scheme P_3D_ENERGY */
#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 240.0 /* scaling factor for energy representation */
#define E_SCALE 20.0 /* scaling factor for energy representation */
#define LOG_SCALE 0.75 /* scaling factor for energy log representation */
#define LOG_SHIFT 0.5 /* shift of colors on log scale */
#define LOG_ENERGY_FLOOR -10.0 /* floor value for log of (total) energy */
@ -281,7 +280,7 @@
#define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */
#define COLORBAR_RANGE 2.5 /* scale of color scheme bar */
#define COLORBAR_RANGE_B 5.0 /* scale of color scheme bar for 2nd part */
#define COLORBAR_RANGE_B 0.5 /* 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 */
@ -289,7 +288,7 @@
#define ADD_POTENTIAL 1 /* set to 1 to add potential to z coordinate */
// #define POT_MAZE 7
#define POTENTIAL 10
#define POT_FACT 400.0
#define POT_FACT 20.0
/* end of constants only used by sub_wave and sub_maze */
@ -303,14 +302,14 @@ 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] = {8.0, 8.0, 8.0}; /* location of observer for REP_PROJ_3D representation */
double observer[3] = {8.0, 6.0, 5.0}; /* location of observer for REP_PROJ_3D representation */
int reset_view = 0; /* switch to reset 3D view parameters (for option ROTATE_VIEW) */
#define Z_SCALING_FACTOR 0.5 /* overall scaling factor of z axis for REP_PROJ_3D representation */
#define XY_SCALING_FACTOR 1.8 /* overall scaling factor for on-screen (x,y) coordinates after projection */
#define Z_SCALING_FACTOR 0.1 /* overall scaling factor of z axis for REP_PROJ_3D representation */
#define XY_SCALING_FACTOR 1.7 /* overall scaling factor for on-screen (x,y) coordinates after projection */
#define ZMAX_FACTOR 1.0 /* max value of z coordinate for REP_PROJ_3D representation */
#define XSHIFT_3D -0.1 /* overall x shift for REP_PROJ_3D representation */
#define YSHIFT_3D 0.2 /* overall y shift for REP_PROJ_3D representation */
#define XSHIFT_3D 0.2 /* overall x shift for REP_PROJ_3D representation */
#define YSHIFT_3D 0.5 /* overall y shift for REP_PROJ_3D representation */
#include "global_pdes.c" /* constants and global variables */
@ -1067,8 +1066,8 @@ void animation()
// for (j=1; j<NPOLY; j++)
// add_circular_wave_mod(1.0, lambda1*cos(((double)j+0.5)*angle), lambda1*sin(((double)j+0.5)*angle), phi, psi, xy_in);
init_wave_flat_mod(phi, psi, xy_in);
// init_circular_wave_mod(0.0, 0.0, phi, psi, xy_in);
// init_wave_flat_mod(phi, psi, xy_in);
init_circular_wave_mod(-0.5, 0.0, phi, psi, xy_in);
// add_circular_wave_mod(1.0, 1.0, 0.0, phi, psi, xy_in);
// printf("Wave initialized\n");
@ -1152,47 +1151,47 @@ void animation()
// if ((ADD_OSCILLATING_SOURCE)&&(i%OSCILLATING_SOURCE_PERIOD == OSCILLATING_SOURCE_PERIOD - 1))
if ((ADD_OSCILLATING_SOURCE)&&(i%OSCILLATING_SOURCE_PERIOD == 1))
{
// if (ALTERNATE_OSCILLATING_SOURCE) sign = -sign;
// add_circular_wave_mod(sign, 0.0, 0.0, phi, psi, xy_in);
if (ALTERNATE_OSCILLATING_SOURCE) sign = -sign;
add_circular_wave_mod(sign, -0.5, 0.0, phi, psi, xy_in);
phase_shift = 1/30.0;
// phase_shift = 1/60.0;
if (first_source) for (k=0; k<25; k++)
{
angle = 0.0;
omega = DPI/25.0;
wave_source[k].xc = 0.05*cos((double)k*omega);;
wave_source[k].yc = 0.05*sin((double)k*omega);
wave_source[k].phase = 0.99 - 1.4*sin(0.7*(1.0 + wave_source[k].xc/0.05));
wave_source[k].amp = 1.0;
// if (wave_source[k].phase > 0.0) wave_source[k].sign = 1;
// else wave_source[k].sign = -1;
wave_source[k].sign = 1;
}
first_source = 0;
for (k=0; k<25; k++)
wave_source[k].phase += 1.4*sin(0.7*(1.0 + wave_source[k].xc*cos(angle)/0.05 + wave_source[k].yc*sin(angle)/0.05));
angle = DPI*(double)i/(double)NSTEPS;
// phase_shift = 0.02 + 0.08*(double)i/(double)NSTEPS;
// printf("Phase shift = %.3lg\n", phase_shift);
for (k=0; k<25; k++)
{
// wave_source[k].phase += 0.07;
wave_source[k].phase += phase_shift;
wave_source[k].phase -= 1.4*sin(0.7*(1.0 + wave_source[k].xc*cos(angle)/0.05 + wave_source[k].yc*sin(angle)/0.05));
if (wave_source[k].phase > 1.0)
{
add_circular_wave_mod((double)wave_source[k].sign*wave_source[k].amp, wave_source[k].xc, wave_source[k].yc, phi, psi, xy_in);
printf("Adding wave at (%.2lg, %.2lg)\n", wave_source[k].xc, wave_source[k].yc);
wave_source[k].phase -= 1.0;
wave_source[k].sign *= -1;
}
}
// phase_shift = 1/30.0;
// // phase_shift = 1/60.0;
//
// if (first_source) for (k=0; k<25; k++)
// {
// angle = 0.0;
// omega = DPI/25.0;
// wave_source[k].xc = 0.05*cos((double)k*omega);;
// wave_source[k].yc = 0.05*sin((double)k*omega);
// wave_source[k].phase = 0.99 - 1.4*sin(0.7*(1.0 + wave_source[k].xc/0.05));
// wave_source[k].amp = 1.0;
// // if (wave_source[k].phase > 0.0) wave_source[k].sign = 1;
// // else wave_source[k].sign = -1;
// wave_source[k].sign = 1;
// }
// first_source = 0;
//
// for (k=0; k<25; k++)
// wave_source[k].phase += 1.4*sin(0.7*(1.0 + wave_source[k].xc*cos(angle)/0.05 + wave_source[k].yc*sin(angle)/0.05));
//
// angle = DPI*(double)i/(double)NSTEPS;
// // phase_shift = 0.02 + 0.08*(double)i/(double)NSTEPS;
// // printf("Phase shift = %.3lg\n", phase_shift);
//
// for (k=0; k<25; k++)
// {
// // wave_source[k].phase += 0.07;
// wave_source[k].phase += phase_shift;
// wave_source[k].phase -= 1.4*sin(0.7*(1.0 + wave_source[k].xc*cos(angle)/0.05 + wave_source[k].yc*sin(angle)/0.05));
//
// if (wave_source[k].phase > 1.0)
// {
// add_circular_wave_mod((double)wave_source[k].sign*wave_source[k].amp, wave_source[k].xc, wave_source[k].yc, phi, psi, xy_in);
// printf("Adding wave at (%.2lg, %.2lg)\n", wave_source[k].xc, wave_source[k].yc);
// wave_source[k].phase -= 1.0;
// wave_source[k].sign *= -1;
// }
// }
// for (j=0; j<NPOLY; j++)
// add_circular_wave_mod(sign, lambda1*cos(((double)j+0.5)*angle), lambda1*sin(((double)j+0.5)*angle), phi, psi, xy_in);

150
wave_billiard.c
View File

@ -49,7 +49,7 @@
#define NO_EXTRA_BUFFER_SWAP 1 /* some OS require one less buffer swap when recording images */
#define VARIABLE_IOR 0 /* set to 1 for a variable index of refraction */
#define IOR 7 /* choice of index of refraction, see list in global_pdes.c */
#define IOR 9 /* choice of index of refraction, see list in global_pdes.c */
#define IOR_TOTAL_TURNS 1.5 /* total angle of rotation for IOR_PERIODIC_WELLS_ROTATING */
#define MANDEL_IOR_SCALE -0.05 /* parameter controlling dependence of IoR on Mandelbrot escape speed */
@ -61,8 +61,8 @@
#define NX 3840 /* number of grid points on x axis */
#define NY 2300 /* number of grid points on y axis */
#define XMIN -2.0
#define XMAX 2.0 /* x interval */
#define XMIN -1.0
#define XMAX 3.0 /* x interval */
#define YMIN -1.197916667
#define YMAX 1.197916667 /* y interval for 9/16 aspect ratio */
@ -87,7 +87,7 @@
/* Choice of the billiard table */
#define B_DOMAIN 999 /* choice of domain shape, see list in global_pdes.c */
#define B_DOMAIN 10 /* choice of domain shape, see list in global_pdes.c */
#define CIRCLE_PATTERN 1 /* pattern of circles or polygons, see list in global_pdes.c */
@ -96,11 +96,11 @@
#define CIRCLE_PATTERN_B 0 /* second pattern of circles or polygons */
#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 NPOISSON 1000 /* number of points for Poisson C_RAND_POISSON arrangement */
#define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */
#define LAMBDA 0.5 /* parameter controlling the dimensions of domain */
#define MU 0.5 /* parameter controlling the dimensions of domain */
#define LAMBDA 0.1197916667 /* parameter controlling the dimensions of domain */
#define MU 0.035 /* parameter controlling the dimensions of domain */
#define NPOLY 6 /* number of sides of polygon */
#define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */
#define MDEPTH 6 /* depth of computation of Menger gasket */
@ -108,8 +108,9 @@
#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 12 /* number of grid point for grid of disks */
#define NGRIDY 12 /* number of grid point for grid of disks */
#define NGRIDX 30 /* number of grid point for grid of disks */
#define NGRIDY 30 /* number of grid point for grid of disks */
// #define NGRIDY 18 /* number of grid point for grid of disks */
#define X_SHOOTER -0.2
#define Y_SHOOTER -0.6
@ -130,14 +131,14 @@
#define TWOSPEEDS 0 /* 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 OSCILLATION_SCHEDULE 1 /* oscillation schedule, see list in global_pdes.c */
#define OSCILLATION_SCHEDULE 0 /* oscillation schedule, see list in global_pdes.c */
#define OMEGA 0.0005 /* frequency of periodic excitation */
#define AMPLITUDE 0.8 /* amplitude of periodic excitation */
#define OMEGA 0.004 /* frequency of periodic excitation */
#define AMPLITUDE 1.0 /* amplitude of periodic excitation */
#define ACHIRP 0.25 /* acceleration coefficient in chirp */
#define DAMPING 0.0 /* damping of periodic excitation */
#define COURANT 0.06 /* Courant number */
#define COURANTB 0.0 /* Courant number in medium B */
#define COURANT 0.1 /* Courant number */
#define COURANTB 0.01 /* Courant number in medium B */
#define GAMMA 0.0 /* damping factor in wave equation */
#define GAMMAB 0.0 /* damping factor in wave equation */
#define GAMMA_SIDES 1.0e-4 /* damping factor on boundary */
@ -151,7 +152,7 @@
/* For similar wave forms, COURANT^2*GAMMA should be kept constant */
#define ADD_OSCILLATING_SOURCE 1 /* set to 1 to add an oscillating wave source */
#define OSCILLATING_SOURCE_PERIOD 2 /* period of oscillating source */
#define OSCILLATING_SOURCE_PERIOD 30 /* period of oscillating source */
#define ALTERNATE_OSCILLATING_SOURCE 1 /* set to 1 to alternate sign of oscillating source */
#define ADD_WAVE_PACKET_SOURCES 0 /* set to 1 to add several sources emitting wave packets */
@ -165,9 +166,9 @@
/* Parameters for length and speed of simulation */
#define NSTEPS 2750 /* number of frames of movie */
// #define NSTEPS 2100 /* number of frames of movie */
#define NVID 10 /* number of iterations between images displayed on screen */
#define NSTEPS 2400 /* number of frames of movie */
// #define NSTEPS 1300 /* number of frames of movie */
#define NVID 6 /* number of iterations between images displayed on screen */
#define NSEG 1000 /* number of segments of boundary */
#define INITIAL_TIME 0 /* time after which to start saving frames */
#define BOUNDARY_WIDTH 2 /* width of billiard boundary */
@ -184,22 +185,23 @@
/* Parameters of initial condition */
#define INITIAL_AMP 0.007 /* amplitude of initial condition */
#define INITIAL_VARIANCE 0.0003 /* variance of initial condition */
#define INITIAL_WAVELENGTH 0.015 /* wavelength of initial condition */
#define INITIAL_AMP 2.0 /* amplitude of initial condition */
// #define INITIAL_VARIANCE 0.000015 /* variance of initial condition */
#define INITIAL_VARIANCE 0.000025 /* variance of initial condition */
#define INITIAL_WAVELENGTH 0.05 /* wavelength of initial condition */
/* Plot type, see list in global_pdes.c */
#define PLOT 0
// #define PLOT 7
#define PLOT_B 7 /* plot type for second movie */
#define PLOT_B 5 /* plot type for second movie */
/* Color schemes */
// #define COLOR_PALETTE 15 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE 11 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 0 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE 17 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 13 /* Color palette, see list in global_pdes.c */
#define BLACK 1 /* background */
@ -212,7 +214,7 @@
#define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */
#define E_SCALE 100.0 /* scaling factor for energy representation */
#define LOG_SCALE 1.0 /* scaling factor for energy log representation */
#define LOG_SHIFT 3.5 /* shift of colors on log scale */
#define LOG_SHIFT 1.0 /* shift of colors on log scale */
#define FLUX_SCALE 5.0e3 /* scaling factor for enegy flux represtnation */
#define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */
@ -224,9 +226,11 @@
#define HUEAMP -180.0 /* amplitude of variation of hue for color scheme C_HUE */
#define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */
#define COLORBAR_RANGE 2.5 /* scale of color scheme bar */
#define COLORBAR_RANGE_B 0.1 /* scale of color scheme bar for 2nd part */
#define COLORBAR_RANGE 2.0 /* scale of color scheme bar */
#define COLORBAR_RANGE_B 2.0 /* scale of color scheme bar for 2nd part */
#define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */
#define CIRC_COLORBAR 0 /* set to 1 to draw circular color scheme */
#define CIRC_COLORBAR_B 0 /* set to 1 to draw circular color scheme */
#define SAVE_TIME_SERIES 0 /* set to 1 to save wave time series at a point */
@ -562,19 +566,15 @@ void draw_color_bar(int plot, double range)
// else draw_color_scheme(1.7, YMIN + 0.25, 1.9, YMAX - 0.25, plot, -range, range);
}
// void draw_color_bar_palette(int plot, double range, int palette)
// {
// if (ROTATE_COLOR_SCHEME) draw_color_scheme_palette(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range, palette);
// else draw_color_scheme_palette(XMAX - 0.3, YMIN + 0.1, XMAX - 0.1, YMAX - 0.1, plot, -range, range, palette);
// }
void draw_color_bar_palette(int plot, double range, int palette, int fade, double fade_value)
void draw_color_bar_palette(int plot, double range, int palette, int circular, int fade, double fade_value)
{
double width = 0.14;
// double width = 0.2;
if (ROTATE_COLOR_SCHEME)
draw_color_scheme_palette_fade(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range, palette, fade, fade_value);
else if (circular)
draw_circular_color_scheme_palette_fade(XMAX - 2.0*width, YMIN + 2.0*width, 1.5*width, plot, -range, range, palette, fade, fade_value);
else
draw_color_scheme_palette_fade(XMAX - 1.5*width, YMIN + 0.1, XMAX - 0.5*width, YMAX - 0.1, plot, -range, range, palette, fade, fade_value);
}
@ -676,9 +676,9 @@ void animation()
// 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_flat(phi, psi, xy_in);
// init_circular_wave(0.0, 0.0, phi, psi, xy_in);
init_circular_wave(-0.5, 0.0, phi, psi, xy_in);
// x = XMIN + (XMAX - XMIN)*rand()/RAND_MAX;
// y = YMIN + (YMAX - YMIN)*rand()/RAND_MAX;
// init_circular_wave(0.0, -0.8, phi, psi, xy_in);
@ -733,7 +733,7 @@ void animation()
draw_billiard(0, 1.0);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, fade, fade_value);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, fade, fade_value);
if (PRINT_SPEED)
{
a = 0.0075;
@ -784,48 +784,48 @@ void animation()
draw_billiard(0, 1.0);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, fade, fade_value);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, fade, fade_value);
/* add oscillating waves */
if ((ADD_OSCILLATING_SOURCE)&&(i%OSCILLATING_SOURCE_PERIOD == 1))
{
// if (ALTERNATE_OSCILLATING_SOURCE) sign = -sign;
// add_circular_wave(sign, 0.0, 0.0, phi, psi, xy_in);
if (ALTERNATE_OSCILLATING_SOURCE) sign = -sign;
add_circular_wave(sign, -0.5, 0.0, phi, psi, xy_in);
// p = phased_array_schedule(i);
// phase_shift = 0.02 + 0.06*(double)i/(double)NSTEPS;
if (first_source) for (k=0; k<25; k++)
{
omega = DPI/25.0;
wave_source[k].xc = 0.05*cos((double)k*omega);;
wave_source[k].yc = 0.05*sin((double)k*omega);
wave_source[k].phase = 0.99 - 1.4*sin(0.7*(1.0 + wave_source[k].xc/0.05));
wave_source[k].amp = 1.0;
if (wave_source[k].phase) wave_source[k].sign = 1;
else wave_source[k].sign = -1;
first_source = 0;
}
for (k=0; k<25; k++)
wave_source[k].phase += 1.4*sin(0.7*(1.0 + wave_source[k].xc*cos(angle)/0.05 + wave_source[k].yc*sin(angle)/0.05));
angle = DPI*(double)i/(double)NSTEPS;
for (k=0; k<25; k++)
{
wave_source[k].phase += 0.07;
wave_source[k].phase -= 1.4*sin(0.7*(1.0 + wave_source[k].xc*cos(angle)/0.05 + wave_source[k].yc*sin(angle)/0.05));
if (wave_source[k].phase > 1.0)
{
add_circular_wave((double)wave_source[k].sign*wave_source[k].amp, wave_source[k].xc, wave_source[k].yc, phi, psi, xy_in);
printf("Adding wave at (%.2lg, %.2lg)\n", wave_source[k].xc, wave_source[k].yc);
wave_source[k].phase -= 1.0;
wave_source[k].sign *= -1;
}
}
// if (first_source) for (k=0; k<25; k++)
// {
// omega = DPI/25.0;
// wave_source[k].xc = 0.05*cos((double)k*omega);;
// wave_source[k].yc = 0.05*sin((double)k*omega);
// wave_source[k].phase = 0.99 - 1.4*sin(0.7*(1.0 + wave_source[k].xc/0.05));
// wave_source[k].amp = 1.0;
// if (wave_source[k].phase) wave_source[k].sign = 1;
// else wave_source[k].sign = -1;
// first_source = 0;
// }
//
// for (k=0; k<25; k++)
// wave_source[k].phase += 1.4*sin(0.7*(1.0 + wave_source[k].xc*cos(angle)/0.05 + wave_source[k].yc*sin(angle)/0.05));
//
// angle = DPI*(double)i/(double)NSTEPS;
//
// for (k=0; k<25; k++)
// {
// wave_source[k].phase += 0.07;
// wave_source[k].phase -= 1.4*sin(0.7*(1.0 + wave_source[k].xc*cos(angle)/0.05 + wave_source[k].yc*sin(angle)/0.05));
//
// if (wave_source[k].phase > 1.0)
// {
// add_circular_wave((double)wave_source[k].sign*wave_source[k].amp, wave_source[k].xc, wave_source[k].yc, phi, psi, xy_in);
// printf("Adding wave at (%.2lg, %.2lg)\n", wave_source[k].xc, wave_source[k].yc);
// wave_source[k].phase -= 1.0;
// wave_source[k].sign *= -1;
// }
// }
// p = 3;
// y = -1.0;
@ -893,7 +893,7 @@ void animation()
draw_wave_highres_palette(2, phi, psi, total_energy, total_flux, xy_in, scale, i, PLOT_B, COLOR_PALETTE_B, 0, 1.0);
else draw_wave_epalette(phi, psi, total_energy, total_flux, color_scale, xy_in, scale, i, PLOT_B, COLOR_PALETTE_B, 0, 1.0);
draw_billiard(0, 1.0);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, CIRC_COLORBAR_B, 0, 1.0);
if (PRINT_SPEED) print_speed(speed, 0, 1.0);
if (PRINT_FREQUENCY) print_frequency(phase_shift, 0, 1.0);
glutSwapBuffers();
@ -922,7 +922,7 @@ void animation()
if (HIGHRES) draw_wave_highres_palette(2, phi, psi, total_energy, total_flux, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 0, 1.0);
else draw_wave_epalette(phi, psi, total_energy, total_flux, color_scale, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 0, 1.0);
draw_billiard(0, 1.0);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, 0, 1.0);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, 0, 1.0);
if (PRINT_SPEED) print_speed(speed, 0, 1.0);
if (PRINT_FREQUENCY) print_frequency(phase_shift, 0, 1.0);
glutSwapBuffers();
@ -935,7 +935,7 @@ void animation()
draw_wave_highres_palette(2, phi, psi, total_energy, total_flux, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 1, fade_value);
else draw_wave_epalette(phi, psi, total_energy, total_flux, color_scale, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 1, fade_value);
draw_billiard(1, fade_value);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, 1, fade_value);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, 1, fade_value);
if (PRINT_SPEED) print_speed(speed, 1, fade_value);
if (PRINT_FREQUENCY) print_frequency(phase_shift, 1, fade_value);
if (!NO_EXTRA_BUFFER_SWAP) glutSwapBuffers();
@ -948,7 +948,7 @@ void animation()
draw_wave_highres_palette(2, phi, psi, total_energy, total_flux, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B, 0, 1.0);
else draw_wave_epalette(phi, psi, total_energy, total_flux, color_scale, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B, 0, 1.0);
draw_billiard(0, 1.0);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, CIRC_COLORBAR_B, 0, 1.0);
if (PRINT_SPEED) print_speed(speed, 0, 1.0);
if (PRINT_FREQUENCY) print_frequency(phase_shift, 0, 1.0);
glutSwapBuffers();
@ -961,7 +961,7 @@ void animation()
draw_wave_highres_palette(2, phi, psi, total_energy, total_flux, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B, 1, fade_value);
else draw_wave_epalette(phi, psi, total_energy, total_flux, color_scale, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B, 1, fade_value);
draw_billiard(1, fade_value);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 1, fade_value);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, CIRC_COLORBAR_B, 1, fade_value);
if (PRINT_SPEED) print_speed(speed, 1, fade_value);
if (PRINT_FREQUENCY) print_frequency(phase_shift, 0, 1.0);
glutSwapBuffers();

3
wave_common.c
View File

@ -386,8 +386,7 @@ void draw_wave(double *phi[NX], double *psi[NX], short int *xy_in[NX], double sc
glEnd ();
}
void draw_wave_e(double *phi[NX], double *psi[NX], double *total_energy[NX], double *color_scale[NX], short int *xy_in[NX],
double scale, int time, int plot)
void draw_wave_e(double *phi[NX], double *psi[NX], double *total_energy[NX], double *color_scale[NX], short int *xy_in[NX], double scale, int time, int plot)
/* draw the field, new version with total energy option */
{
int i, j, iplus, iminus, jplus, jminus;

272
wave_comparison.c
View File

@ -41,40 +41,47 @@
#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 SAVE_MEMORY 0 /* set to 1 to save memory when writing tiff images */
#define MOVIE 1 /* set to 1 to generate movie */
#define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */
#define SAVE_MEMORY 1 /* set to 1 to save memory when writing tiff images */
#define NO_EXTRA_BUFFER_SWAP 1 /* some OS require one less buffer swap when recording images */
#define TIME_LAPSE 0 /* set to 1 to add a time-lapse movie at the end */
#define TIME_LAPSE_FACTOR 4 /* factor of time-lapse movie */
// #define WINWIDTH 1920 /* window width */
// #define WINHEIGHT 1000 /* window height */
#define WINWIDTH 1920 /* window width */
#define WINHEIGHT 1150 /* window height */
// #define NX 1920 /* number of grid points on x axis */
// #define NY 1000 /* number of grid points on y axis */
// #define YMID 500 /* mid point of display */
// #define XMIN -0.5
// #define XMAX 3.5 /* x interval */
// #define YMIN -1.041666667
// #define YMAX 1.041666667 /* y interval for 9/16 aspect ratio */
#define NX 3840 /* number of grid points on x axis */
#define NY 2300 /* number of grid points on y axis */
#define YMID 1150 /* mid point of display */
#define XMIN -2.0
#define XMAX 2.0 /* x interval */
#define YMIN -1.197916667
#define YMAX 1.197916667 /* y interval for 9/16 aspect ratio */
#define WINWIDTH 1280 /* window width */
#define WINHEIGHT 720 /* window height */
#define NX 1280 /* number of grid points on x axis */
#define NY 720 /* number of grid points on y axis */
#define YMID 360 /* mid point of display */
#define XMIN -0.5
#define XMAX 3.5 /* x interval */
#define YMIN -1.125
#define YMAX 1.125 /* y interval for 9/16 aspect ratio */
// #define WINWIDTH 1280 /* window width */
// #define WINHEIGHT 720 /* window height */
// #define NX 1280 /* number of grid points on x axis */
// #define NY 720 /* number of grid points on y axis */
// #define YMID 360 /* mid point of display */
// #define XMIN -2.0
// #define XMAX 2.0 /* x interval */
// #define YMIN -1.125
// #define YMAX 1.125 /* y interval for 9/16 aspect ratio */
#define HIGHRES 1 /* set to 1 if resolution of grid is double that of displayed image */
#define JULIA_SCALE 1.0 /* scaling for Julia sets */
/* Choice of the billiard table */
#define B_DOMAIN 43 /* choice of domain shape, see list in global_pdes.c */
#define B_DOMAIN_B 47 /* choice of domain shape, see list in global_pdes.c */
#define B_DOMAIN 999 /* choice of domain shape, see list in global_pdes.c */
#define B_DOMAIN_B 999 /* choice of domain shape, see list in global_pdes.c */
#define CIRCLE_PATTERN 13 /* pattern of circles, see list in global_pdes.c */
#define CIRCLE_PATTERN_B 13 /* pattern of circles, see list in global_pdes.c */
@ -123,13 +130,13 @@
/* Physical parameters of wave equation */
#define TWOSPEEDS 1 /* set to 1 to replace hardcore boundary by medium with different speed */
#define OSCILLATE_LEFT 1 /* set to 1 to add oscilating boundary condition on the left */
#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.004 /* frequency of periodic excitation */
#define AMPLITUDE 1.0 /* amplitude of periodic excitation */
#define COURANT 0.02 /* Courant number */
#define COURANTB 0.01154 /* Courant number in medium B */
#define COURANT 0.045 /* Courant number */
#define COURANTB 0.045 /* Courant number in medium B */
// #define COURANTB 0.005 /* Courant number in medium B */
// #define COURANTB 0.008 /* Courant number in medium B */
#define GAMMA 0.0 /* damping factor in wave equation */
@ -148,17 +155,22 @@
/* Increasing COURANT speeds up the simulation, but decreases accuracy */
/* For similar wave forms, COURANT^2*GAMMA should be kept constant */
#define ADD_OSCILLATING_SOURCE 1 /* set to 1 to add an oscillating wave source */
#define OSCILLATING_SOURCE_PERIOD 2 /* period of oscillating source */
#define ALTERNATE_OSCILLATING_SOURCE 1 /* set to 1 to alternate sign of oscillating source */
#define NSOURCES 48 /* number of sources */
/* Boundary conditions, see list in global_pdes.c */
#define B_COND 0
#define B_COND 4
/* Parameters for length and speed of simulation */
#define NSTEPS 2500 /* number of frames of movie */
// #define NSTEPS 3300 /* number of frames of movie */
// #define NSTEPS 500 /* number of frames of movie */
#define NSTEPS 2600 /* number of frames of movie */
#define NVID 25 /* number of iterations between images displayed on screen */
#define NSEG 100 /* number of segments of boundary */
#define INITIAL_TIME 20 /* time after which to start saving frames */
#define INITIAL_TIME 0 /* time after which to start saving frames */
#define COMPUTE_ENERGIES 0 /* set to 1 to compute and print energies */
#define BOUNDARY_WIDTH 2 /* width of billiard boundary */
@ -168,24 +180,25 @@
#define SLEEP2 1 /* final sleeping time */
#define MID_FRAMES 20 /* number of still frames between movies */
#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.75 /* amplitude of initial condition */
// #define INITIAL_VARIANCE 0.0003 /* variance of initial condition */
// #define INITIAL_WAVELENGTH 0.015 /* wavelength of initial condition */
#define INITIAL_VARIANCE 0.0003 /* variance of initial condition */
#define INITIAL_WAVELENGTH 0.02 /* wavelength of initial condition */
#define INITIAL_AMP 0.002 /* amplitude of initial condition */
// #define INITIAL_AMP 0.002 /* amplitude of initial condition */
#define INITIAL_VARIANCE 0.0003 /* variance of initial condition */
#define INITIAL_WAVELENGTH 0.015 /* wavelength of initial condition */
/* Plot type, see list in global_pdes.c */
#define PLOT 0
#define PLOT_B 1
#define PLOT_B 3
/* Color schemes */
#define COLOR_PALETTE 18 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 13 /* Color palette, see list in global_pdes.c */
#define BLACK 1 /* background */
#define BLACK_TEXT 1 /* set to 1 to write text in black instead of white */
@ -198,7 +211,7 @@
#define PHASE_FACTOR 1.0 /* factor in computation of phase in color scheme P_3D_PHASE */
#define PHASE_SHIFT 0.0 /* shift of phase in color scheme P_3D_PHASE */
#define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */
#define E_SCALE 200.0 /* scaling factor for energy representation */
#define E_SCALE 500.0 /* scaling factor for energy representation */
#define LOG_SCALE 1.5 /* scaling factor for energy log representation */
#define LOG_SHIFT 1.0 /* shift of colors on log scale */
#define FLUX_SCALE 1.0e4 /* scaling factor for enegy flux represtnation */
@ -213,7 +226,7 @@
#define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */
#define COLORBAR_RANGE 1.5 /* scale of color scheme bar */
#define COLORBAR_RANGE_B 2.5 /* scale of color scheme bar for 2nd part */
#define COLORBAR_RANGE_B 12.5 /* scale of color scheme bar for 2nd part */
#define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */
@ -238,16 +251,16 @@
#define POT_MAZE 7
#define POTENTIAL 0
#define MAZE_WIDTH 0.02 /* half width of maze walls */
#define VARIABLE_IOR 0 /* set to 1 for a variable index of refraction */
#define VARIABLE_IOR 1 /* set to 1 for a variable index of refraction */
#define IOR 7 /* choice of index of refraction, see list in global_pdes.c */
#define IOR_TOTAL_TURNS 1.5 /* total angle of rotation for IOR_PERIODIC_WELLS_ROTATING */
#define MANDEL_IOR_SCALE -0.05 /* parameter controlling dependence of IoR on Mandelbrot escape speed */
#define COURANT 0.04 /* Courant number */
#define COURANTB 0.0 /* Courant number in medium B */
#define INITIAL_AMP 0.5 /* amplitude of initial condition */
#define INITIAL_VARIANCE 0.0003 /* variance of initial condition */
#define INITIAL_WAVELENGTH 0.015 /* wavelength of initial condition */
#define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */
// #define COURANT 0.04 /* Courant number */
// #define COURANTB 0.0 /* Courant number in medium B */
// #define INITIAL_AMP 0.5 /* amplitude of initial condition */
// #define INITIAL_VARIANCE 0.0003 /* variance of initial condition */
// #define INITIAL_WAVELENGTH 0.015 /* wavelength of initial condition */
// #define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */
#define WAVE_PACKET_SOURCE_TYPE 1 /* type of wave packet sources */
#define N_WAVE_PACKETS 15 /* number of wave packets */
/* end of constants only used by sub_wave and sub_maze */
@ -261,6 +274,7 @@
double courant2, courantb2; /* Courant parameters squared */
FILE *monitor_sources;
/*********************/
/* animation part */
@ -677,20 +691,32 @@ void evolve_wave(double *phi[NX], double *psi[NX], double *tmp[NX], short int *x
}
void draw_color_bar(int plot, double range)
// void draw_color_bar(int plot, double range)
// {
// if (ROTATE_COLOR_SCHEME) draw_color_scheme(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range);
// else draw_color_scheme(XMAX - 0.3, YMIN + 0.1, XMAX - 0.1, YMAX - 0.1, plot, -range, range);
// }
void draw_color_bar_palette(int plot, double range, int palette, int fade, double fade_value)
{
if (ROTATE_COLOR_SCHEME) draw_color_scheme(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range);
else draw_color_scheme(XMAX - 0.3, YMIN + 0.1, XMAX - 0.1, YMAX - 0.1, plot, -range, range);
double width = 0.14;
// double width = 0.2;
if (ROTATE_COLOR_SCHEME)
draw_color_scheme_palette_fade(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range, palette, fade, fade_value);
else
draw_color_scheme_palette_fade(XMAX - 1.5*width, YMIN + 0.1, XMAX - 0.5*width, YMAX - 0.1, plot, -range, range, palette, fade, fade_value);
}
void animation()
{
double time, scale, energies[6], top_energy, bottom_energy;
double *phi[NX], *psi[NX], *tmp[NX];
double time, scale, energies[6], top_energy, bottom_energy, omega, angle, fade_value;
double *phi[NX], *psi[NX], *tmp[NX], *total_energy[NX];
short int *xy_in[NX];
int i, j, s, counter = 0;
int i, j, s, counter = 0, k, first_source = 1, fade, resol = HIGHRES + 1;
t_wave_source wave_source[NSOURCES];
monitor_sources = fopen("monitor_sources.dat", "w");
/* Since NX and NY are big, it seemed wiser to use some memory allocation here */
for (i=0; i<NX; i++)
@ -698,6 +724,7 @@ void animation()
phi[i] = (double *)malloc(NY*sizeof(double));
psi[i] = (double *)malloc(NY*sizeof(double));
tmp[i] = (double *)malloc(NY*sizeof(double));
total_energy[i] = (double *)malloc(NY*sizeof(double));
xy_in[i] = (short int *)malloc(NY*sizeof(short int));
}
@ -711,6 +738,11 @@ void animation()
npolyline = init_polyline(MDEPTH, polyline);
for (i=0; i<npolyline; i++) printf("vertex %i: (%.3f, %.3f)\n", i, polyline[i].x, polyline[i].y);
/* initialize total energy table */
if ((PLOT == P_MEAN_ENERGY)||(PLOT_B == P_MEAN_ENERGY)||(PLOT == P_LOG_MEAN_ENERGY)||(PLOT_B == P_LOG_MEAN_ENERGY))
for (i=0; i<NX; i++)
for (j=0; j<NY; j++)
total_energy[i][j] = 0.0;
courant2 = COURANT*COURANT;
courantb2 = COURANTB*COURANTB;
@ -749,7 +781,7 @@ void animation()
printf("drawing billiard\n");
draw_billiard_comp();
if (DRAW_COLOR_SCHEME) draw_color_bar(PLOT, COLORBAR_RANGE);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, fade, fade_value);
glutSwapBuffers();
@ -769,11 +801,68 @@ void animation()
}
else scale = 1.0;
draw_wave_comp(phi, psi, xy_in, scale, i, PLOT);
// draw_wave_comp(phi, psi, xy_in, scale, i, PLOT);
draw_wave_comp_highres_palette(resol, phi, psi, total_energy, xy_in, scale, i, PLOT, COLOR_PALETTE, 0, 1.0);
for (j=0; j<NVID; j++)
{
evolve_wave(phi, psi, tmp, xy_in);
// if (i % 10 == 9) oscillate_linear_wave(0.2*scale, 0.15*(double)(i*NVID + j), -1.5, YMIN, -1.5, YMAX, phi, psi);
}
/* add oscillating waves */
if ((ADD_OSCILLATING_SOURCE)&&(i%OSCILLATING_SOURCE_PERIOD == 1))
{
// if (ALTERNATE_OSCILLATING_SOURCE) sign = -sign;
// add_circular_wave(sign, 0.0, 0.0, phi, psi, xy_in);
// p = phased_array_schedule(i);
// phase_shift = 0.02 + 0.06*(double)i/(double)NSTEPS;
if (first_source) for (k=0; k<NSOURCES; k++)
{
omega = DPI/(double)NSOURCES;
wave_source[k].xc = 0.05*cos(((double)k + 0.5)*omega);;
wave_source[k].yc = 0.05*sin(((double)k + 0.5)*omega);
if (wave_source[k].yc < 0.0)
{
wave_source[k].xc *= 0.5;
wave_source[k].yc *= 0.5;
}
// wave_source[k].phase = 0.99 - 1.4*sin(0.35*(1.0 + wave_source[k].xc/0.1));
// wave_source[k].phase = 0.99 - 1.4*sin(0.7*(1.0 + wave_source[k].xc/0.1));
// wave_source[k].phase = 0.99 - 1.4*sin(0.7*(1.0 + wave_source[k].xc/0.05));
wave_source[k].phase = 0.99 - 1.4*sin(0.35*(1.0 + wave_source[k].xc/0.05));
wave_source[k].amp = 1.0;
// if (wave_source[k].phase) wave_source[k].sign = 1;
// else wave_source[k].sign = -1;
wave_source[k].sign = 1;
first_source = 0;
}
fprintf(monitor_sources, "Frame %i\n\n", i);
for (k=0; k<NSOURCES; k++) /*if (wave_source[k].xc > 0.0) */
{
wave_source[k].phase += 0.06;
// if (k==1) printf("x = %.3lg, phase = %.3lg\n", wave_source[k].xc, wave_source[k].phase);
// fprintf(monitor_sources, "x = %.3lg, y = %.3lg, phase = %.3lg\n", wave_source[k].xc, wave_source[k].yc, wave_source[k].phase);
if (wave_source[k].phase > 1.0)
{
add_circular_wave_comp((double)wave_source[k].sign*wave_source[k].amp, wave_source[k].xc, wave_source[k].yc, phi, psi, xy_in, (wave_source[k].yc > 0));
fprintf(monitor_sources, "x = %.3lg, y = %.3lg, phase = %.3lg\n", wave_source[k].xc, wave_source[k].yc, wave_source[k].phase);
printf("Adding wave at (%.2lg, %.2lg)\n", wave_source[k].xc, wave_source[k].yc);
fprintf(monitor_sources, "Adding wave at (%.2lg, %.2lg)\n", wave_source[k].xc, wave_source[k].yc);
wave_source[k].phase -= 1.0;
wave_source[k].sign *= -1;
}
}
}
draw_billiard_comp();
if (COMPUTE_ENERGIES)
{
compute_energy_tblr(phi, psi, xy_in, energies);
@ -785,37 +874,39 @@ void animation()
print_energies(energies, top_energy, bottom_energy);
}
for (j=0; j<NVID; j++)
{
evolve_wave(phi, psi, tmp, xy_in);
// if (i % 10 == 9) oscillate_linear_wave(0.2*scale, 0.15*(double)(i*NVID + j), -1.5, YMIN, -1.5, YMAX, phi, psi);
}
if (DRAW_COLOR_SCHEME) draw_color_bar(PLOT, COLORBAR_RANGE);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, fade, fade_value);
glutSwapBuffers();
if (!((NO_EXTRA_BUFFER_SWAP)&&(MOVIE))) glutSwapBuffers();
if (MOVIE)
{
if (i >= INITIAL_TIME)
if (i >= INITIAL_TIME) save_frame();
if ((i >= INITIAL_TIME)&&(DOUBLE_MOVIE))
{
save_frame();
// save_frame();
if ((TIME_LAPSE)&&((i - INITIAL_TIME)%TIME_LAPSE_FACTOR == 0))
{
save_frame_counter(NSTEPS + END_FRAMES + (i - INITIAL_TIME)/TIME_LAPSE_FACTOR);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0);
counter++;
}
else if (DOUBLE_MOVIE)
else /*if (DOUBLE_MOVIE)*/
{
draw_wave_comp(phi, psi, xy_in, scale, i, PLOT_B);
// draw_wave_comp(phi, psi, xy_in, scale, i, PLOT_B);
draw_wave_comp_highres_palette(resol, phi, psi, total_energy, xy_in, scale, i, PLOT_B, COLOR_PALETTE_B, 0, 1.0);
draw_billiard_comp();
if (DRAW_COLOR_SCHEME) draw_color_bar(PLOT_B, COLORBAR_RANGE_B);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0);
glutSwapBuffers();
save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter);
counter++;
}
}
else printf("Initial phase time %i of %i\n", i, INITIAL_TIME);
else if (NO_EXTRA_BUFFER_SWAP) glutSwapBuffers();
// else printf("Initial phase time %i of %i\n", i, INITIAL_TIME);
@ -835,20 +926,42 @@ void animation()
{
if (DOUBLE_MOVIE)
{
draw_wave_comp(phi, psi, xy_in, scale, i, PLOT);
// draw_wave_comp(phi, psi, xy_in, scale, i, PLOT);
draw_wave_comp_highres_palette(resol, phi, psi, total_energy, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 0, 1.0);
draw_billiard_comp();
if (DRAW_COLOR_SCHEME) draw_color_bar(PLOT, COLORBAR_RANGE);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, 0, 1.0);
glutSwapBuffers();
}
for (i=0; i<MID_FRAMES; i++) save_frame();
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_comp_highres_palette(resol, phi, psi, total_energy, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 1, fade_value);
draw_billiard_comp();
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, 1, fade_value);
if (!NO_EXTRA_BUFFER_SWAP) glutSwapBuffers();
save_frame_counter(NSTEPS + i + 1);
}
if (DOUBLE_MOVIE)
{
draw_wave_comp(phi, psi, xy_in, scale, i, PLOT_B);
// draw_wave_comp(phi, psi, xy_in, scale, i, PLOT_B);
draw_wave_comp_highres_palette(resol, phi, psi, total_energy, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B, 0, 1.0);
draw_billiard_comp();
if (DRAW_COLOR_SCHEME) draw_color_bar(PLOT_B, COLORBAR_RANGE_B);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_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_comp_highres_palette(resol, phi, psi, total_energy, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B, 1, fade_value);
draw_billiard_comp();
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 1, fade_value);
glutSwapBuffers();
save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter + i);
}
}
for (i=0; i<END_FRAMES; i++) save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter + i);
if (TIME_LAPSE) for (i=0; i<END_FRAMES; i++) save_frame_counter(NSTEPS + END_FRAMES + NSTEPS/TIME_LAPSE_FACTOR + i);
s = system("mv wave*.tif tif_wave/");
@ -858,6 +971,7 @@ void animation()
free(phi[i]);
free(psi[i]);
free(tmp[i]);
free(total_energy[i]);
free(xy_in[i]);
}
@ -866,6 +980,12 @@ void animation()
void display(void)
{
time_t rawtime;
struct tm * timeinfo;
time(&rawtime);
timeinfo = localtime(&rawtime);
glPushMatrix();
blank();
@ -879,7 +999,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));
}

6
wave_energy.c
View File

@ -219,12 +219,6 @@
#define IOR 7 /* choice of index of refraction, see list in global_pdes.c */
#define IOR_TOTAL_TURNS 1.5 /* total angle of rotation for IOR_PERIODIC_WELLS_ROTATING */
#define MANDEL_IOR_SCALE -0.05 /* parameter controlling dependence of IoR on Mandelbrot escape speed */
#define COURANT 0.04 /* Courant number */
#define COURANTB 0.0 /* Courant number in medium B */
#define INITIAL_AMP 0.5 /* amplitude of initial condition */
#define INITIAL_VARIANCE 0.0003 /* variance of initial condition */
#define INITIAL_WAVELENGTH 0.015 /* wavelength of initial condition */
#define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */
#define WAVE_PACKET_SOURCE_TYPE 1 /* type of wave packet sources */
#define N_WAVE_PACKETS 15 /* number of wave packets */
/* end of constants only used by sub_wave and sub_maze */