Compare commits
3 Commits
b809ce9e55
...
2bf029d952
| Author | SHA1 | Date |
|---|---|---|
Nils Berglund
|
2bf029d952 | |
|
Nils Berglund
|
a925cf2ca9 | |
|
Nils Berglund
|
de54999d98 |
8846
Parameters.md
8846
Parameters.md
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
2424
Parameters_May23.md
2424
Parameters_May23.md
File diff suppressed because it is too large
Load Diff
20
global_3d.c
20
global_3d.c
|
|
@ -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;
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
{
|
||||
|
|
|
|||
|
|
@ -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 */
|
||||
|
||||
|
|
|
|||
|
|
@ -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
246
lennardjones.c
|
|
@ -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))
|
||||
{
|
||||
|
|
|
|||
|
|
@ -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
621
rde.c
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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
472
sub_lj.c
|
|
@ -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");
|
||||
|
||||
}
|
||||
|
|
|
|||
323
sub_wave.c
323
sub_wave.c
|
|
@ -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++){
|
||||
|
|
|
|||
|
|
@ -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];
|
||||
|
|
|
|||
|
|
@ -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
132
sub_wave_comp.c
|
|
@ -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
163
wave_3d.c
|
|
@ -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
150
wave_billiard.c
|
|
@ -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();
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
|
|
|||
|
|
@ -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));
|
||||
}
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -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 */
|
||||
|
|
|
|||
Loading…
Reference in New Issue