From 6d0d707fcc958539036a295b0d31c48ec0949baf Mon Sep 17 00:00:00 2001 From: Nils Berglund <83530463+nilsberglund-orleans@users.noreply.github.com> Date: Sat, 25 Mar 2023 19:56:19 +0100 Subject: [PATCH] Add files via upload --- drop_billiard.c | 1 + global_3d.c | 6 + global_ljones.c | 14 +- global_particles.c | 1 + global_pdes.c | 26 +- heat.c | 14 + lennardjones.c | 110 ++--- mangrove.c | 15 +- particle_billiard.c | 267 +++++++++--- particle_pinball.c | 1 + rde.c | 312 +++++++++----- schrodinger.c | 13 + sub_lj.c | 997 +++++++++++++++++++++++++++++++++++++++----- sub_maze.c | 295 ++++++++++++- sub_part_billiard.c | 404 +++++++++++++++++- sub_rde.c | 433 ++++++++++++++++--- sub_wave.c | 644 +++++++++++++++++++++++++++- sub_wave_3d.c | 121 +++++- wave_3d.c | 181 ++++---- wave_billiard.c | 195 +++++---- wave_comparison.c | 15 +- wave_energy.c | 15 +- 22 files changed, 3491 insertions(+), 589 deletions(-) diff --git a/drop_billiard.c b/drop_billiard.c index 53ba3ad..e82ec1f 100644 --- a/drop_billiard.c +++ b/drop_billiard.c @@ -132,6 +132,7 @@ #define MAZE_XSHIFT 0.0 /* horizontal shift of maze */ #define MAZE_RANDOM_FACTOR 0.1 /* randomization factor for S_MAZE_RANDOM */ #define MAZE_CORNER_RADIUS 0.5 /* radius of tounded corners in maze */ +#define CLOSE_MAZE 0 /* set to 1 to close maze exits */ diff --git a/global_3d.c b/global_3d.c index 859b3bf..203ef5c 100644 --- a/global_3d.c +++ b/global_3d.c @@ -37,6 +37,7 @@ #define POT_FERMIONS 5 /* two interacting 1D fermions */ #define POT_FERMIONS_PERIODIC 6 /* two interacting 1D fermions on the circle */ #define POT_MAZE 7 /* higher potential on walls of a maze */ +#define POT_IOR 10 /* index of refraction, for z coordinate of wave equation */ /* Choice of vector potential */ @@ -47,6 +48,10 @@ #define GF_VERTICAL 0 /* gravity acting vertically */ #define GF_CIRCLE 1 /* repelling circle */ +#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 */ /* macros to avoid unnecessary computations in 3D plots */ @@ -88,6 +93,7 @@ typedef struct double flux_int_table[FLUX_WINDOW]; /* table of energy flux intensities (for averaging) */ short int flux_counter; /* counter for averaging of energy flux */ double rgb[3]; /* RGB color code */ + double *potential; /* pointer to "potential" to add to z-coordinate */ double *p_zfield[2]; /* pointers to z field (second pointer for option DOUBLE_MOVIE) */ double *p_cfield[4]; /* pointers to color field (second pointer for option DOUBLE_MOVIE) */ /* third and fourth pointer for color luminosity (for energy flux) */ diff --git a/global_ljones.c b/global_ljones.c index 8366ffc..ae2b099 100644 --- a/global_ljones.c +++ b/global_ljones.c @@ -11,11 +11,12 @@ #define D_CIRCLES 20 /* several circles */ #define D_CIRCLES_IN_RECT 201 /* several circles in a rectangle */ -#define NMAXCIRCLES 20000 /* total number of circles/polygons (must be at least NCX*NCY for square grid) */ +#define NMAXCIRCLES 100000 /* total number of circles/polygons (must be at least NCX*NCY for square grid) */ #define MAXNEIGH 20 /* max number of neighbours kept in memory */ #define NMAXOBSTACLES 100 /* max number of obstacles */ #define NMAXSEGMENTS 1000 /* max number of repelling segments */ #define NMAXGROUPS 50 /* max number of groups of segments */ +#define NMAXCOLLISIONS 200000 /* max number of collisions */ #define C_SQUARE 0 /* square grid of circles */ #define C_HEX 1 /* hexagonal/triangular grid of circles */ @@ -105,6 +106,8 @@ #define TH_VERTICAL 0 /* only particles at the right of x = PARTIAL_THERMO_SHIFT are coupled */ #define TH_INSEGMENT 1 /* only particles in region defined by segments are coupled */ #define TH_INBOX 2 /* only particles in a given box are coupled */ +#define TH_LAYER 3 /* only particles above -LAMBDA are coupled */ +#define TH_LAYER_TYPE2 4 /* only particles above highest type 2 particle are coupled */ /* Gravity schedules */ @@ -138,6 +141,13 @@ #define CHEM_AABAA 6 /* reaction A + A <-> B (reversible) */ #define CHEM_POLYMER 7 /* reaction A + B -> C, A + C -> D, etc */ #define CHEM_POLYMER_DISS 8 /* polimerisation with dissociation */ +#define CHEM_POLYMER_STEP 9 /* step growth polimerisation with dissociation */ +#define CHEM_AUTOCATALYSIS 10 /* autocatalytic reaction 2A + B -> 2B */ +#define CHEM_CATALYTIC_A2D 11 /* catalytic reaction A + B -> C, A + C -> B + D */ +#define CHEM_ABCAB 12 /* reaction A + B <-> C (reversible) */ +#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 */ /* Initial conditions for chemical reactions */ @@ -147,6 +157,8 @@ #define IC_RANDOM_TWO 2 /* particle type chosen randomly between 1 and 2, with TYPE_PROPORTION */ #define IC_CIRCLE 3 /* type 1 in a disc */ #define IC_CATALYSIS 4 /* mix of 1 and 2 in left half, only 1 in right half */ +#define IC_LAYERS 5 /* layer of 2 below 1 */ +#define IC_BZ 6 /* initial state for BZ reaction */ /* Plot types */ diff --git a/global_particles.c b/global_particles.c index 010340c..41b1499 100644 --- a/global_particles.c +++ b/global_particles.c @@ -98,6 +98,7 @@ double x_shooter = -0.2, y_shooter = -0.6, x_target = 0.4, y_target = 0.7; #define P_MAZE_CIRCULAR 13 /* circular maze */ #define P_MAZE_CIRC_SCATTERER 14 /* circular maze with scatterers */ #define P_MAZE_HEX 15 /* hexagonal maze */ +#define P_MAZE_OCT 16 /* maze with octagonal and square cells */ /* Color palettes */ diff --git a/global_pdes.c b/global_pdes.c index fde056a..fdcdd1a 100644 --- a/global_pdes.c +++ b/global_pdes.c @@ -77,6 +77,9 @@ #define D_LENSES_RING 59 /* several lenses forming a ring */ #define D_MAZE_CIRCULAR 60 /* circular maze */ +#define D_WING 70 /* complement of wing-shaped domain */ +#define D_TESLA 71 /* Tesla valve */ + #define NMAXCIRCLES 10000 /* total number of circles/polygons (must be at least NCX*NCY for square grid) */ #define NMAXPOLY 50000 /* maximal number of vertices of polygonal lines (for von Koch et al) */ @@ -119,6 +122,10 @@ #define IOR_MANDELBROT_LIN 100 /* index of refraction depends on escape time in Mandelbrot set (linear) */ #define IOR_EARTH 2 /* index of refraction models speed of seismic waves */ #define IOR_EXPLO_LENSING 3 /* explosive lensing */ +#define IOR_PERIODIC_WELLS 4 /* periodic superposition of "wells" */ +#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 */ /* Boundary conditions */ @@ -221,15 +228,18 @@ #define Z_EULER_PRESSURE 54 /* pressure */ /* for Euler compressible Euler equation */ -#define Z_EULER_DENSITY 60 /* density */ -#define Z_EULER_SPEED 61 /* norm of velocity */ +#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 */ /* special boundary conditions for Euler equation */ #define BCE_TOPBOTTOM 1 /* laminar flow at top and bottom */ #define BCE_TOPBOTTOMLEFT 2 /* laminar flow at top, bottom and left side */ #define BCE_CHANNELS 3 /* laminar flow in channels at left and right */ #define BCE_MIDDLE_STRIP 4 /* laminar flow in horizontal strip in the middle */ +#define BCE_LEFT 5 /* laminar flow at left side */ typedef struct { @@ -280,6 +290,18 @@ typedef struct } t_laplacian; +typedef struct +{ + double xc, yc; /* (x,y) coordinates of center */ + int ix, iy; /* lattice coordinates of center */ + double period, amp; /* period and amplitude */ + double phase; /* phase shift */ + double var_envelope; /* variance of Gaussian envelope */ + int time_shift; /* time shift */ +} t_wave_packet; + + + int ncircles = NMAXCIRCLES; /* actual number of circles, can be decreased e.g. for random patterns */ int npolyline = NMAXPOLY; /* actual length of polyline */ int npolyrect = NMAXPOLY; /* actual number of polyrect */ diff --git a/heat.c b/heat.c index dd81957..57e2093 100644 --- a/heat.c +++ b/heat.c @@ -195,9 +195,23 @@ #define MAZE_MAX_NGBH 4 /* max number of neighbours of maze cell */ #define RAND_SHIFT 24 /* seed of random number generator */ #define MAZE_XSHIFT 0.0 /* horizontal shift of maze */ +#define MAZE_WIDTH 0.02 /* half width of maze walls */ #define ADD_POTENTIAL 0 #define POT_MAZE 7 #define POTENTIAL 0 +#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 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 */ #include "global_pdes.c" diff --git a/lennardjones.c b/lennardjones.c index b4aaa4d..cd20d96 100644 --- a/lennardjones.c +++ b/lennardjones.c @@ -58,10 +58,10 @@ #define YMIN -1.125 #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ -#define INITXMIN -1.92 -#define INITXMAX 1.92 /* x interval for initial condition */ -#define INITYMIN -1.05 -#define INITYMAX 1.05 /* y interval for initial condition */ +#define INITXMIN -1.97 +#define INITXMAX 1.97 /* x interval for initial condition */ +#define INITYMIN -1.1 +#define INITYMAX 1.1 /* y interval for initial condition */ #define BCXMIN -2.0 #define BCXMAX 2.0 /* x interval for boundary condition */ @@ -84,7 +84,7 @@ #define NOZZLE_SHAPE_B 4 /* 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.66 /* proportion of particles of first type */ +#define TYPE_PROPORTION 0.3 /* 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 */ @@ -97,11 +97,11 @@ #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 4.7 /* minimal distance in Poisson disc process, controls density of particles */ +#define PDISC_DISTANCE 4.2 /* minimal distance in Poisson disc process, controls density of particles */ #define PDISC_CANDIDATES 100 /* number of candidates in construction of Poisson disc process */ #define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */ -#define LAMBDA 0.8 /* parameter controlling the dimensions of domain */ +#define LAMBDA 0.5 /* parameter controlling the dimensions of domain */ #define MU 0.008 /* parameter controlling radius of particles */ #define MU_B 0.012 /* parameter controlling radius of particles of second type */ #define NPOLY 25 /* number of sides of polygon */ @@ -111,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 90 /* number of grid point for grid of disks */ -#define NGRIDY 100 /* number of grid point for grid of disks */ +#define NGRIDX 120 /* number of grid point for grid of disks */ +#define NGRIDY 51 /* 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,10 +125,11 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 1000 /* number of frames of movie */ -#define NVID 150 /* number of iterations between images displayed on screen */ +#define NSTEPS 5000 /* number of frames of movie */ +// #define NSTEPS 3000 /* number of frames of movie */ +#define NVID 175 /* number of iterations between images displayed on screen */ #define NSEG 250 /* number of segments of boundary */ -#define INITIAL_TIME 25 /* time after which to start saving frames */ +#define INITIAL_TIME 20 /* time after which to start saving frames */ #define OBSTACLE_INITIAL_TIME 200 /* time after which to start moving obstacle */ #define BOUNDARY_WIDTH 1 /* width of particle boundary */ #define LINK_WIDTH 2 /* width of links between particles */ @@ -143,7 +144,7 @@ /* Boundary conditions, see list in global_ljones.c */ -#define BOUNDARY_COND 3 +#define BOUNDARY_COND 0 /* Plot type, see list in global_ljones.c */ @@ -175,6 +176,9 @@ #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_TEMPERATURE 0 /* set to 1 to print current temperature */ + /* particle properties */ #define ENERGY_HUE_MIN 330.0 /* color of original particle */ @@ -182,24 +186,21 @@ #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 70.0 /* hue of particles of type 0 */ -#define HUE_TYPE1 150.0 /* hue of particles of type 1 */ -#define HUE_TYPE2 190.0 /* hue of particles of type 2 */ -#define HUE_TYPE3 220.0 /* hue of particles of type 3 */ -// #define HUE_TYPE1 310.0 /* hue of particles of type 1 */ -// #define HUE_TYPE2 150.0 /* hue of particles of type 2 */ -// #define HUE_TYPE3 180.0 /* hue of particles of type 3 */ +#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 RANDOM_RADIUS 0 /* set to 1 for random circle radius */ #define DT_PARTICLE 3.0e-6 /* time step for particle displacement */ #define KREPEL 12.0 /* constant in repelling force between particles */ -#define EQUILIBRIUM_DIST 3.5 /* Lennard-Jones equilibrium distance */ -#define EQUILIBRIUM_DIST_B 3.5 /* Lennard-Jones equilibrium distance for second type of particle */ +#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 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 0.0 /* damping coefficient of particles during initial phase */ +#define DAMPING 200.0 /* damping coefficient of particles */ +#define INITIAL_DAMPING 1000.0 /* damping coefficient of particles during initial phase */ #define PARTICLE_MASS 1.0 /* mass of particle of radius MU */ -#define PARTICLE_MASS_B 2.0 /* mass of particle of radius MU */ +#define PARTICLE_MASS_B 0.5 /* mass of particle of radius MU */ #define PARTICLE_INERTIA_MOMENT 0.02 /* moment of inertia of particle */ #define PARTICLE_INERTIA_MOMENT_B 0.02 /* moment of inertia of second type of particle */ #define V_INITIAL 0.0 /* initial velocity range */ @@ -208,12 +209,11 @@ #define THERMOSTAT 1 /* set to 1 to switch on thermostat */ #define VARY_THERMOSTAT 0 /* set to 1 for time-dependent thermostat schedule */ #define SIGMA 5.0 /* noise intensity in thermostat */ -#define BETA 0.01 /* initial inverse temperature */ +#define BETA 0.002 /* 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 10.0 /* radius in which to count neighbours */ -// #define NBH_DIST_FACTOR 7.5 /* 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 INCREASE_GRAVITY 0 /* set to 1 to increase gravity during the simulation */ @@ -236,7 +236,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.025 /* factor by which to change BETA during simulation */ +#define BETA_FACTOR 0.5 /* 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 */ @@ -260,13 +260,13 @@ #define RECORD_PRESSURES 0 /* set to 1 to record pressures on obstacle */ #define N_PRESSURES 100 /* number of intervals to record pressure */ #define N_P_AVERAGE 100 /* size of pressure averaging window */ -#define N_T_AVERAGE 50 /* size of temperature averaging window */ +#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 particles to the right of obstacle to thermostat */ -#define PARTIAL_THERMO_REGION 1 /* region for partial thermostat coupling (see list in global_ljones.c) */ +#define PARTIAL_THERMO_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_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.2 /* vertical size of partial thermostat coupling */ +#define PARTIAL_THERMO_HEIGHT 0.25 /* vertical size of partial thermostat coupling */ #define INCREASE_KREPEL 0 /* set to 1 to increase KREPEL during simulation */ #define KREPEL_FACTOR 1000.0 /* factor by which to change KREPEL during simulation */ @@ -327,19 +327,21 @@ #define PRINT_ENTROPY 0 /* set to 1 to compute entropy */ #define REACTION_DIFFUSION 1 /* set to 1 to simulate a chemical reaction (particles may change type) */ -#define RD_REACTION 8 /* type of reaction, see list in global_ljones.c */ -#define RD_TYPES 9 /* number of types in reaction-diffusion equation */ +#define RD_REACTION 15 /* 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 REACION_DIST 3.5 /* maximal distance for reaction to occur */ +#define REACTION_DIST 3.5 /* maximal distance for reaction to occur */ #define REACTION_PROB 0.5 /* probability controlling reaction term */ -#define DISSOCIATION_PROB 0.005 /* probability controlling dissociation reaction */ +#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 COLLISION_TIME 25 /* time during which collisions are shown */ +#define EXOTHERMIC 0 /* set to 1 to make reaction exo/endothermic */ +#define DELTA_EKIN 500.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 1 /* set to 1 to make of plot of particle number over time */ -#define PARTICLE_NB_PLOT_FACTOR 1.0 /* expected final number of particles over initial number */ -#define PRINT_LEFT 0 /* set to 1 to print certain parameters at the top left instead of right */ +#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 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 */ @@ -365,8 +367,8 @@ #define FLOOR_OMEGA 0 /* set to 1 to limit particle momentum to PMAX */ #define PMAX 1000.0 /* maximal force */ -#define HASHX 90 /* size of hashgrid in x direction */ -#define HASHY 45 /* 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 */ @@ -1049,8 +1051,7 @@ 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; + fboundary = 0.0, pleft = 0.0, pright = 0.0, entropy[2], mean_energy, gravity = GRAVITY, speed_ratio, ymin, ymax, delta_energy; 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, @@ -1093,8 +1094,8 @@ void animation() if (REACTION_DIFFUSION) { - collisions = (t_collision *)malloc(2*NMAXCIRCLES*sizeof(t_collision)); - for (i=0; i<2*NMAXCIRCLES; i++) collisions[i].time = 0; + collisions = (t_collision *)malloc(2*NMAXCOLLISIONS*sizeof(t_collision)); + for (i=0; i<2*NMAXCOLLISIONS; i++) collisions[i].time = 0; } if (SAVE_TIME_SERIES) @@ -1364,7 +1365,8 @@ void animation() /* case of reaction-diffusion equation */ if ((i > INITIAL_TIME)&&(REACTION_DIFFUSION)) { - ncollisions = update_types(particle, collisions, ncollisions); + ncollisions = update_types(particle, collisions, ncollisions, particle_numbers, i - INITIAL_TIME - 1, &delta_energy); + if (EXOTHERMIC) beta *= 1.0/(1.0 + delta_energy/totalenergy); nactive = 0; for (j=0; j #define MOVIE 0 /* set to 1 to generate movie */ -#define SAVE_MEMORY 1 /* set to 1 to save memory when writing tiff images */ +#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 WINHEIGHT 720 /* window height */ -#define XMIN -2.0 -#define XMAX 2.0 /* x interval */ +#define XMIN -1.5 +#define XMAX 2.5 /* x interval */ #define YMIN -1.125 #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ -// #define WINWIDTH 1920 /* window width */ -// #define WINHEIGHT 1000 /* window height */ -// -// #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 SCALING_FACTOR 1.0 /* scaling factor of drawing, needed for flower billiards, otherwise set to 1.0 */ /* Choice of the billiard table, see global_particles.c */ @@ -59,14 +54,12 @@ #define B_DOMAIN 31 /* choice of domain shape */ #define CIRCLE_PATTERN 1 /* pattern of circles */ -#define POLYLINE_PATTERN 15 /* pattern of polyline */ +#define POLYLINE_PATTERN 10 /* pattern of polyline */ #define ABSORBING_CIRCLES 0 /* set to 1 for circular scatterers to be absorbing */ #define NMAXCIRCLES 100000 /* total number of circles (must be at least NCX*NCY for square grid) */ #define NMAXPOLY 100000 /* total number of sides of polygonal line */ -// #define NCX 10 /* number of circles in x direction */ -// #define NCY 10 /* number of circles in y direction */ #define NCX 30 /* number of circles in x direction */ #define NCY 20 /* number of circles in y direction */ #define NPOISSON 500 /* number of points for Poisson C_RAND_POISSON arrangement */ @@ -89,37 +82,35 @@ /* Simulation parameters */ -// #define NPART 1 /* number of particles */ -#define NPART 20000 /* number of particles */ -// #define NPART 10000 /* number of particles */ +// #define NPART 10 /* number of particles */ +#define NPART 5000 /* number of particles */ #define NPARTMAX 100000 /* maximal number of particles after resampling */ #define LMAX 0.01 /* minimal segment length triggering resampling */ #define DMIN 0.02 /* minimal distance to boundary for triggering resampling */ #define CYCLE 1 /* set to 1 for closed curve (start in all directions) */ #define SHOWTRAILS 0 /* set to 1 to keep trails of the particles */ #define 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_PARTICLES 1 /* set to 1 to draw particles in heat map */ -#define HEATMAP_MAX_PART_BY_CELL 0 /* to draw only limited number of particles in cell */ -#define PLOT_HEATMAP_AVERAGE 0 /* set to 1 to plot average number of 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 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 */ -#define PRINT_LEFT_RIGHT_PARTICLE_NUMBER 1 /* set to 1 to print number of particles on left and right side */ +#define PRINT_LEFT_RIGHT_PARTICLE_NUMBER 0 /* set to 1 to print number of particles on left and right side */ #define PRINT_CIRCLE_PARTICLE_NUMBER 0 /* set to 1 to print number of particles outside circular maze */ #define PRINT_COLLISION_NUMBER 0 /* set to 1 to print number of collisions */ #define TEST_ACTIVE 1 /* set to 1 to test whether particle is in billiard */ #define TEST_INITIAL_COND 0 /* set to 1 to allow only initial conditions that pass a test */ -#define NSTEPS 12300 /* number of frames of movie */ -// #define NSTEPS 6500 /* number of frames of movie */ -#define TIME 1500 /* time between movie frames, for fluidity of real-time simulation */ -// #define TIME 750 /* time between movie frames, for fluidity of real-time simulation */ +#define NSTEPS 22000 /* number of frames of movie */ +#define TIME 2000 /* time between movie frames, for fluidity of real-time simulation */ +// #define DPHI 0.000002 /* integration step */ #define DPHI 0.00002 /* integration step */ -// #define DPHI 0.00005 /* integration step */ -// #define DPHI 0.00001 /* integration step */ #define NVID 25 /* number of iterations between images displayed on screen */ -// #define NVID 100 /* number of iterations between images displayed on screen */ -#define END_FRAMES 100 /* number of still frames at the end of the movie */ +#define END_FRAMES 50 /* number of still frames at the end of the movie */ /* Decreasing TIME accelerates the animation and the movie */ /* For constant speed of movie, TIME*DPHI should be kept constant */ @@ -131,15 +122,14 @@ #define COLOR_PALETTE 17 /* Color palette, see list in global_pdes.c */ -#define NCOLORS 1000 /* number of colors */ +#define NCOLORS 500 /* number of colors */ #define COLORSHIFT 0 /* hue of initial color */ #define COLOR_HUEMIN 0 /* minimal color hue */ -#define COLOR_HUEMAX 150 /* maximal color hue */ -#define RAINBOW_COLOR 0 /* set to 1 to use different colors for all particles */ +#define COLOR_HUEMAX 160 /* 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 */ #define LENGTH 0.025 /* length of velocity vectors */ -// #define LENGTH 0.04 /* length of velocity vectors */ #define BILLIARD_WIDTH 2 /* width of billiard */ #define PARTICLE_WIDTH 2 /* width of particles */ #define FRONT_WIDTH 3 /* width of wave front */ @@ -155,14 +145,14 @@ #define SLEEP1 1 /* initial sleeping time */ #define SLEEP2 1 /* final sleeping time */ -#define NXMAZE 24 /* width of maze */ -#define NYMAZE 20 /* height of maze */ -#define MAZE_MAX_NGBH 6 /* max number of neighbours of maze cell */ -#define RAND_SHIFT 11 /* seed of random number generator */ -// #define RAND_SHIFT 0 /* seed of random number generator */ +#define NXMAZE 36 /* width of maze */ +#define NYMAZE 36 /* 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 */ #define MAZE_RANDOM_FACTOR 0.1 /* randomization factor for S_MAZE_RANDOM */ #define MAZE_CORNER_RADIUS 0.5 /* radius of tounded corners in maze */ +#define CLOSE_MAZE 1 /* set to 1 to close maze exits */ #include "global_particles.c" #include "sub_maze.c" @@ -563,7 +553,8 @@ void draw_config_heatmap(double *configs[NPARTMAX], int active[NPARTMAX], int he { drawtable[i] = ((n == -2)||((n >= -1)&&(heatmap_number[n] <= HEATMAP_MAX_PART_BY_CELL))); } - else drawtable[i] = (n >= -1); + else drawtable[i] = 1; +// else drawtable[i] = (n >= -1); } // printf("Particle %i is in maze cell %i\n", i, n); @@ -596,6 +587,156 @@ void draw_config_heatmap(double *configs[NPARTMAX], int active[NPARTMAX], int he free(drawtable); } +double plot_coord(double x, double xmin, double xmax) +{ + return(xmin + x*(xmax - xmin)); +} + +void draw_chosen_heatmap_histogram(int heatmap_number[NXMAZE*NYMAZE+1], int normalisation) +{ + int i, j, k, n, bin, nbins, binwidth, maxpart, part_number, maxhist = 0, *histo; + double xmin, xmax, ymin, ymax, xmid, ymid, dx, dy, plotxmin, plotxmax, plotymin, plotymax, c; + double x1, x2, y, y1, y2, hue, rgb[3]; + static int first = 1, prevmaxhist, prevmaxpart; + static double minprop = 0.01; + char message[100]; + + if (first) + { + xmin = XMAX - 1.25; + xmax = XMAX - 0.05; + ymin = YMAX - 1.25; + ymax = YMAX - 0.05; + + xmid = 0.5*(xmin + xmax); + ymid = 0.5*(ymin + ymax); + + dx = 0.5*(xmax - xmin); + dy = 0.5*(ymax - ymin); + + plotxmin = xmin + 0.15; + plotxmax = xmax - 0.1; + plotymin = ymin + 0.07; + plotymax = ymax - 0.1; + + prevmaxhist = 1010; + prevmaxpart = 10; + + first = 0; + } + +// nbins = NXMAZE; +// nbins = (int)(2.0*sqrt((double)NPART/(double)NXMAZE)); + nbins = (int)(NBIN_FACTOR*sqrt((double)NPART/(double)(NXMAZE*NYMAZE))); + + histo = (int *)malloc(nbins*sizeof(int)); + + for (i=0; i maxpart) maxpart = heatmap_number[j]/normalisation; + + if (maxpart < 1010) maxpart = 1010; + + c = log((double)maxpart)/(double)nbins; + + for (j=0; j 0) + { + i = (int)(log((double)n)/c) - 2; + if (i < 0) i = 0; + histo[i]++; + } + } + +// for (i=0; i maxhist) maxhist = histo[i]; + if (maxhist < 10) maxhist = 10; + + /* some smoothing in time */ + maxpart = (maxpart + 3*prevmaxpart)/4; + prevmaxpart = maxpart; + maxhist = (maxhist + 3*prevmaxhist)/4; + prevmaxhist = maxhist; + + glColor3f(0.0, 0.0, 0.0); + glLineWidth(1); + + x1 = plotxmin; + y1 = plotymin; + + for (i=0; i 1) sprintf(message, "%i particles", nleft); - else sprintf(message, "%i particle", nleft); - write_text(xl, yl, message); +// if (nleft > 1) sprintf(message, "%i particles", nleft); +// else sprintf(message, "%i particle", nleft); + if (nleft > 1) sprintf(message, "%i part.", nleft); + else sprintf(message, "%i part.", nleft); + write_text_fixedwidth(xl, yl, message); if (nright > 1) sprintf(message, "%i particles", nright); else sprintf(message, "%i particle", nright); - write_text(xr, yr, message); + write_text_fixedwidth(xr, yr, message); } void print_circle_part_number(double *configs[NPARTMAX], int active[NPARTMAX], double xr, double yr) @@ -864,6 +1008,7 @@ void animation() int i, j, resamp = 1, s, i1, i2, c, lengthmax; int *color, *newcolor, *active, *heatmap_number, *heatmap_total; short int *heatmap_visited; + char message[100]; t_exit *exits; // t_circle *circles; /* experimental */ @@ -948,13 +1093,9 @@ void animation() // alphamax = 2.50949; // init_drop_config(x_shooter, y_shooter, alphamax, alphamax + DPI, configs); - init_drop_config(-0.05, 0.05, 0.0, 0.3*PID, configs); -// init_drop_config(-0.5, 0.0, 0.2, 0.4, configs); -// init_drop_config(0.0, 0.05, 0.0, DPI, configs); + init_drop_config(0.05, 0.05, 0.0, DPI, configs); +// init_drop_config(-0.95, 0.95, 0.0, DPI, configs); -// init_drop_config(-1.3, -0.1, 0.0, DPI, configs); -// init_drop_config(1.4, 0.1, 0.0, DPI, configs); -// init_drop_config(0.5, 0.5, -1.0, 1.0, configs); // init_sym_drop_config(-1.0, 0.5, -PID, PID, configs); // init_drop_config(-0.999, 0.0, -alpha, alpha, configs); @@ -972,7 +1113,7 @@ void animation() if (DRAW_BILLIARD) draw_billiard(); if (PRINT_PARTICLE_NUMBER) print_part_number(configs, active, XMIN + 0.1, YMIN + 0.1); else if (PRINT_LEFT_RIGHT_PARTICLE_NUMBER) - print_left_right_part_number(configs, active, XMIN + 0.1, YMIN + 0.05, XMAX - 0.45, YMIN + 0.05, exits); + print_left_right_part_number(configs, active, XMIN + 0.05, YMIN + 0.05, XMAX - 0.35, YMIN + 0.05, exits); else if (PRINT_CIRCLE_PARTICLE_NUMBER) print_circle_part_number(configs, active, XMAX - 0.45, YMIN + 0.05); else if (PRINT_COLLISION_NUMBER) print_collision_number(ncollisions, XMIN + 0.1, YMIN + 0.1); @@ -1029,6 +1170,7 @@ void animation() else if (HEATMAP) { draw_config_heatmap(configs, active, heatmap_number, heatmap_total, heatmap_visited, DRAW_HEATMAP_PARTICLES); + if (DRAW_HEATMAP_HISTOGRAM) draw_heatmap_histogram(heatmap_number, heatmap_total, i); // draw_config(newcolor, configs, active); } else draw_config(newcolor, configs, active); @@ -1036,7 +1178,7 @@ void animation() if (DRAW_BILLIARD) draw_billiard(); if (PRINT_PARTICLE_NUMBER) print_part_number(configs, active, XMIN + 0.1, YMIN + 0.1); else if (PRINT_LEFT_RIGHT_PARTICLE_NUMBER) - print_left_right_part_number(configs, active, XMIN + 0.1, YMIN + 0.05, XMAX - 0.45, YMIN + 0.05, exits); + print_left_right_part_number(configs, active, XMIN + 0.05, YMIN + 0.05, XMAX - 0.35, YMIN + 0.05, exits); else if (PRINT_CIRCLE_PARTICLE_NUMBER) print_circle_part_number(configs, active, XMAX - 0.45, YMIN + 0.05); // print_left_right_part_number(configs, XMIN + 0.1, YMIN + 0.1, XMAX - 0.45, YMIN + 0.1, YMIN + MAZE_XSHIFT, YMAX + MAZE_XSHIFT); else if (PRINT_COLLISION_NUMBER) print_collision_number(ncollisions, XMIN + 0.1, YMIN + 0.1); @@ -1050,7 +1192,8 @@ void animation() if (MOVIE) { - save_frame(); + if (INVERT_COUNTER) save_frame_counter(NSTEPS+1-i); + else save_frame(); /* it seems that saving too many files too fast can cause trouble with the file system */ /* so this is to make a pause from time to time - parameter PAUSE may need adjusting */ @@ -1066,9 +1209,23 @@ void animation() if (MOVIE) { - if (HEATMAP) + if (DRAW_FINAL_HEATMAP) draw_config_heatmap(configs, active, heatmap_number, heatmap_total, heatmap_visited, 0); - for (i=0; i= LAMBDA) + { + *gx = 0.0; + *gy = 0.0; + } + else + { + x1 = 1.0 - x/LAMBDA; + if (x1 < 0.1) x1 = 0.1; + 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 = 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); +// hy = y1/(MU*MU*x1*x1); +// h = module2(hx, hy) + 1.0e-2; +// *gx = G_FIELD*f*hx/h; +// *gy = G_FIELD*f*hy/h; + } + break; + } default: { *gx = 0.0; @@ -565,15 +609,50 @@ void initialize_vector_potential(double vpotential_field[2*NX*NY]) } } -void initialize_gfield(double gfield[2*NX*NY]) +void initialize_gfield(double gfield[2*NX*NY], double bc_field[NX*NY]) /* initialize the exterior field, e.g. for the compressible Euler equation */ { int i, j; + double dx, dy; - #pragma omp parallel for private(i,j) - for (i=0; i 0)) + if (((RDE_EQUATION == E_EULER_INCOMP)||(RDE_EQUATION == E_EULER_COMP))&&(IN_OUT_FLOW_BC > 0)) { switch (IN_OUT_FLOW_BC) { + case (BCE_LEFT): + { + set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, 0.02, 0.1, 1.0, 0.0, 0.1, phi_out, xy_in, 0, 5, 0, NY, IN_OUT_BC_FACTOR); + break; + } case (BCE_TOPBOTTOM): { - set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, LAMINAR_FLOW_YPERIOD, -0.1, phi_out, xy_in, 0, NX, 0, 10); - set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, LAMINAR_FLOW_YPERIOD, -0.1, phi_out, xy_in, 0, NX, NY-10, NY); + set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, 0.02, LAMINAR_FLOW_YPERIOD, 1.0, -0.1, 0.1, phi_out, xy_in, 0, NX, 0, 10, IN_OUT_BC_FACTOR); + set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, 0.02, LAMINAR_FLOW_YPERIOD, 1.0, -0.1, 0.1, phi_out, xy_in, 0, NX, NY-10, NY, IN_OUT_BC_FACTOR); break; } case (BCE_TOPBOTTOMLEFT): { - set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, LAMINAR_FLOW_YPERIOD, -0.1, phi_out, xy_in, 0, NX, 0, 10); - set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, LAMINAR_FLOW_YPERIOD, -0.1, phi_out, xy_in, 0, NX, NY-10, NY); - set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, LAMINAR_FLOW_YPERIOD, -0.1, phi_out, xy_in, 0, 2, 0, NY); + set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, 0.02, LAMINAR_FLOW_YPERIOD, 1.0, -0.1, 0.1, phi_out, xy_in, 0, NX, 0, 10, IN_OUT_BC_FACTOR); + set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, 0.02, LAMINAR_FLOW_YPERIOD, 1.0, -0.1, 0.1, phi_out, xy_in, 0, NX, NY-10, NY, IN_OUT_BC_FACTOR); + set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, 0.02, LAMINAR_FLOW_YPERIOD, 1.0, -0.1, 0.1, phi_out, xy_in, 0, 2, 0, NY, IN_OUT_BC_FACTOR); break; } case (BCE_CHANNELS): { - set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, LAMINAR_FLOW_YPERIOD, 0.0, phi_out, xy_in, imin-5, imin+10, NY - y_maze_entry, y_maze_entry); - set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, LAMINAR_FLOW_YPERIOD, 0.0, phi_out, xy_in, imax-10, imax+5, NY- y_maze_entry, y_maze_entry); + set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, 0.02, LAMINAR_FLOW_YPERIOD, 1.0, 0.0, 0.1, phi_out, xy_in, 0, imin+5, NY - y_channels, y_channels, IN_OUT_BC_FACTOR); + set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, 0.02, LAMINAR_FLOW_YPERIOD, 1.0, 0.0, 0.1, phi_out, xy_in, imax-5, NX - 1, NY- y_channels, y_channels, IN_OUT_BC_FACTOR); +// set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, 0.02, LAMINAR_FLOW_YPERIOD, 1.0, 0.0, 0.1, phi_out, xy_in, imin-5, imin+10, NY - y_channels, y_channels); +// set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, 0.02, LAMINAR_FLOW_YPERIOD, 1.0, 0.0, 0.1, phi_out, xy_in, imax-10, imax+5, NY- y_channels, y_channels); break; } case (BCE_MIDDLE_STRIP): { - set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, LAMINAR_FLOW_YPERIOD, 0.0, phi_out, xy_in, 0, NX, NY/2 - 10, NY/2 + 10); - set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, LAMINAR_FLOW_YPERIOD, 0.0, phi_out, xy_in, 0, 2, 0, NY); - set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, LAMINAR_FLOW_YPERIOD, 0.0, phi_out, xy_in, NX-2, NX, 0, NY); + set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, 0.02, LAMINAR_FLOW_YPERIOD, 1.0, 0.0, 0.1, phi_out, xy_in, 0, NX, NY/2 - 10, NY/2 + 10, IN_OUT_BC_FACTOR); + set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, 0.02, LAMINAR_FLOW_YPERIOD, 1.0, 0.0, 0.1, phi_out, xy_in, 0, 2, 0, NY, IN_OUT_BC_FACTOR); + set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, 0.02, LAMINAR_FLOW_YPERIOD, 1.0, 0.0, 0.1, phi_out, xy_in, NX-2, NX, 0, NY, IN_OUT_BC_FACTOR); break; } } @@ -927,7 +1019,7 @@ void evolve_tracers(double *phi[NFIELDS], double tracers[2*N_TRACERS*NSTEPS], in int tracer, i, j, t, ij[2], iplus, jplus; double x, y, xy[2], vx, vy; - step = 0.01; + step = TRACERS_STEP; for (tracer = 0; tracer < N_TRACERS; tracer++) { @@ -1221,15 +1313,15 @@ void animation() initialize_vector_potential(vector_potential_field); } - if (ADD_FORCE_FIELD) - { - gfield = (double *)malloc(2*NX*NY*sizeof(double)); - initialize_gfield(gfield); - } if (ADAPT_STATE_TO_BC) { bc_field = (double *)malloc(NX*NY*sizeof(double)); - initialize_bcfield(bc_field); + initialize_bcfield(bc_field, polyrect); + } + if (ADD_FORCE_FIELD) + { + gfield = (double *)malloc(2*NX*NY*sizeof(double)); + initialize_gfield(gfield, bc_field); } @@ -1262,7 +1354,8 @@ void animation() // init_shear_flow(1.0, 0.02, 0.15, 1, 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(IN_OUT_FLOW_AMP, LAMINAR_FLOW_MODULATION, 0.02, 0.1, 1.0, 0.0, 0.1, phi, xy_in); + init_laminar_flow(flow_speed_schedule(0), LAMINAR_FLOW_MODULATION, 0.02, 0.1, 1.0, 0.0, 0.1, phi, xy_in); // 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); @@ -1271,7 +1364,7 @@ void animation() init_cfield_rde(phi, xy_in, CPLOT, rde, 0); if (PLOT_3D) init_zfield_rde(phi, xy_in, ZPLOT, rde, 0); - + if (DOUBLE_MOVIE) { init_cfield_rde(phi, xy_in, CPLOT_B, rde, 1); @@ -1296,6 +1389,7 @@ void animation() if (PRINT_PARAMETERS) print_parameters(rde, xy_in, time, 0, VISCOSITY, noise); if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 0, 1.0); + glutSwapBuffers(); sleep(SLEEP1); diff --git a/schrodinger.c b/schrodinger.c index 9a1cafa..e2d02a5 100644 --- a/schrodinger.c +++ b/schrodinger.c @@ -175,9 +175,22 @@ #define MAZE_MAX_NGBH 4 /* max number of neighbours of maze cell */ #define RAND_SHIFT 24 /* seed of random number generator */ #define MAZE_XSHIFT 0.0 /* horizontal shift of maze */ +#define MAZE_WIDTH 0.02 /* half width of maze walls */ #define ADD_POTENTIAL 0 #define POT_MAZE 7 #define POTENTIAL 0 +#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 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 */ diff --git a/sub_lj.c b/sub_lj.c index 19ac849..11350d4 100644 --- a/sub_lj.c +++ b/sub_lj.c @@ -479,6 +479,8 @@ double type_hue(int type) first = 0; } + if ((RD_REACTION == CHEM_CATALYTIC_A2D)&&(type == 4)) return(HUE_TYPE3); + switch (type) { case (0): return(HUE_TYPE0); case (1): return(HUE_TYPE0); @@ -486,13 +488,14 @@ double type_hue(int type) case (3): return(HUE_TYPE2); default: { + if (RD_REACTION == CHEM_BZ) + { + if (type == 7) return(HUE_TYPE2); + if (type == 8) type = 5; + } + else if ((RD_REACTION == CHEM_BRUSSELATOR)&&(type >= 5)) return(70.0); t2 = (double)(type*type); hue = (hmax*t2 - b)/t2; -// b = 4.0*hmax - 4.0*HUE_TYPE3; -// hue = (hmax*t - b)/t; -// hue = (hmax*t - b)/(t+1.0); -// hue = HUE_TYPE3 + 30.0*(double)(type - 4); -// if (hue > 360.0) hue = 360.0; return(hue); } } @@ -2835,8 +2838,7 @@ double water_force(double r, double ca, double sa, double ca_rel, double sa_rel, return(torque); } -int compute_particle_interaction(int i, int k, double force[2], double *torque, t_particle* particle, - double distance, double krepel, double ca, double sa, double ca_rel, double sa_rel) +int compute_particle_interaction(int i, int k, double force[2], double *torque, t_particle* particle, double distance, double krepel, double ca, double sa, double ca_rel, double sa_rel) /* compute repelling force and torque of particle #k on particle #i */ /* returns 1 if distance between particles is smaller than NBH_DIST_FACTOR*MU */ { @@ -3142,13 +3144,13 @@ int add_particle(double x, double y, double vx, double vy, double mass, short in particle[i].angle = DPI*(double)rand()/RAND_MAX; particle[i].omega = 0.0; - if (particle[i].type == 1) - { - particle[i].interaction = INTERACTION_B; - particle[i].eq_dist = EQUILIBRIUM_DIST_B; - particle[i].spin_range = SPIN_RANGE_B; - particle[i].spin_freq = SPIN_INTER_FREQUENCY_B; - } +// if (particle[i].type == 1) +// { +// particle[i].interaction = INTERACTION_B; +// particle[i].eq_dist = EQUILIBRIUM_DIST_B; +// particle[i].spin_range = SPIN_RANGE_B; +// particle[i].spin_freq = SPIN_INTER_FREQUENCY_B; +// } ncircles++; @@ -3465,7 +3467,7 @@ void draw_triangles(t_particle particle[NMAXCIRCLES], int plot) int i, j, k, p, q, t0, tj, tmax, nsame = 0, same_table[9*HASHMAX], width; double rgb[3], hue, dx, dy, x, y, radius, rgbx[3], rgby[3]; - printf("Number of nbs: "); +// printf("Number of nbs: "); for (i=0; i= 2) + { + omega = DPI/(double)(particle.type - 1); + draw_colored_polygon(xc1, yc1, MU, nsides, angle + APOLY*PID, rgb); + for (i=0; i= 2)&&(particle.type <= 4)) + { + omega = DPI/(double)(particle.type - 1); + if (fill) draw_colored_polygon(xc1, yc1, MU, nsides, angle + APOLY*PID, rgb); + else draw_polygon(xc1, yc1, MU, nsides, angle + APOLY*PID); + for (i=0; i 1.0) logratio = 1.0; @@ -4375,7 +4584,7 @@ void print_parameters(double beta, double temperature, double krepel, double len } } - if ((PARTIAL_THERMO_COUPLING)&&(!INCREASE_BETA)) + if ((PARTIAL_THERMO_COUPLING)&&(!INCREASE_BETA)&&(!EXOTHERMIC)) { printf("Temperature %i in average: %.3lg\n", i_temp, temperature); temp[i_temp] = temperature; @@ -5375,7 +5584,7 @@ int initialize_configuration(t_particle particle[NMAXCIRCLES], t_hashgrid hashgr /* initialize all particles, obstacles, and the hashgrid */ { int i, j, k, n, nactive = 0; - double x, y, h, xx, yy; + double x, y, h, xx, yy, rnd; for (i=0; i < ncircles; i++) { @@ -5656,6 +5865,51 @@ int initialize_configuration(t_particle particle[NMAXCIRCLES], t_hashgrid hashgr } break; } + case (IC_LAYERS): + { + if (particle[i].yc > 0.0) + { + particle[i].type = 1; + particle[i].radius = MU; + particle[i].eq_dist = EQUILIBRIUM_DIST; + particle[i].mass_inv = 1.0/PARTICLE_MASS; + } + else if (particle[i].yc < -LAMBDA) + { + particle[i].type = 2; + particle[i].radius = MU_B; + particle[i].eq_dist = EQUILIBRIUM_DIST_B; + particle[i].mass_inv = 1.0/PARTICLE_MASS_B; + } + else particle[i].active = 0; + break; + } + case (IC_BZ): + { + rnd = (double)rand()/(double)RAND_MAX; + if (rnd < 60.0/180.0) + { + particle[i].type = 4; + particle[i].radius = MU; + particle[i].eq_dist = EQUILIBRIUM_DIST; + particle[i].mass_inv = 1.0/(PARTICLE_MASS + 3.0*PARTICLE_MASS_B); + } + else if (rnd < 100.0/180.0) + { + particle[i].type = 5; + particle[i].radius = MU_B; + particle[i].eq_dist = EQUILIBRIUM_DIST; + particle[i].mass_inv = 1.0/(3.5*PARTICLE_MASS); + } + else + { + particle[i].type = 6; + particle[i].radius = 1.5*MU_B; + particle[i].eq_dist = EQUILIBRIUM_DIST; + particle[i].mass_inv = 0.2/(PARTICLE_MASS); + } + break; + } } } @@ -5759,10 +6013,31 @@ int floor_momentum(double p[NMAXCIRCLES]) } int partial_thermostat_coupling(t_particle particle[NMAXCIRCLES], double xmin) -/* only couple particles with x > xmin to thermostat */ +/* only couple particles satisfying condition PARTIAL_THERMO_REGION to thermostat */ { int condition, i, nthermo = 0; - double x, y; + double x, y, height; + static double maxheight; + static int first = 1; + + if (first) + { + maxheight = YMIN + PARTIAL_THERMO_HEIGHT*(YMAX - YMIN); + first = 0; + } + + if (PARTIAL_THERMO_REGION == TH_LAYER_TYPE2) + { + height = YMIN; + for (i=0; i height)&&(y <= maxheight)) + height = y; + } + if (height > maxheight) height = maxheight; + printf("Thermostat region y > %.3lg, max height = %.3lg\n", height, maxheight); + } for (i=0; i PARTIAL_THERMO_HEIGHT); + break; + } + case (TH_LAYER_TYPE2): + { + y = particle[i].yc; + condition = (y > height); + break; + } default: condition = 1; } if (condition) @@ -5873,10 +6160,146 @@ void compute_inverse_masses(double inv_masses[RD_TYPES+1]) } break; } + case (CHEM_POLYMER_STEP): + { + for (type = 2; type < RD_TYPES+1; type++) + inv_masses[type] = 1.0/((double)type*PARTICLE_MASS); + break; + } + case (CHEM_CATALYTIC_A2D): + { + inv_masses[3] = 1.0/(1.0/PARTICLE_MASS + 1.0/PARTICLE_MASS_B); + inv_masses[4] = 0.5/PARTICLE_MASS; + break; + } + case (CHEM_ABCAB): + { + inv_masses[3] = 1.0/(PARTICLE_MASS + PARTICLE_MASS_B); + break; + } + case (CHEM_ABCDABC): + { + inv_masses[3] = 1.0/(PARTICLE_MASS + PARTICLE_MASS_B); + inv_masses[4] = 1.0/(2.0*PARTICLE_MASS + PARTICLE_MASS_B); + break; + } + case (CHEM_BZ): + { + for (type = 2; type <= 4; type++) + { + mass = PARTICLE_MASS + (double)(type-2)*PARTICLE_MASS_B; + inv_masses[type] = 1.0/(mass); + } +// inv_masses[5] = 1.0/(3.5*PARTICLE_MASS); +// inv_masses[6] = 0.2/PARTICLE_MASS; + inv_masses[5] = 0.5/(PARTICLE_MASS); + inv_masses[6] = 0.5/PARTICLE_MASS; + inv_masses[7] = 0.5*inv_masses[3]; + inv_masses[8] = 0.5*inv_masses[5]; + break; + } + case (CHEM_BRUSSELATOR): + { + inv_masses[3] = inv_masses[1]; + inv_masses[4] = 1.0/(PARTICLE_MASS + PARTICLE_MASS_B); + inv_masses[5] = inv_masses[1]; + inv_masses[6] = inv_masses[1]; + break; + } } } -int chem_merge_AAB(int i, int newtype, t_particle particle[NMAXCIRCLES], t_collision *collisions, int ncollisions, double inv_masses[RD_TYPES+1]) +void compute_radii(double radii[RD_TYPES+1]) +/* compute radii of molecules in chemical reactions */ +{ + int type; + + /* default values that apply in most cases */ + radii[1] = MU; + radii[2] = MU_B; + + switch (RD_REACTION) { + case (CHEM_ABC): + { + radii[3] = 1.5*MU; + break; + } + case (CHEM_A2BC): + { + radii[3] = 1.5*MU; + break; + } + case (CHEM_CATALYSIS): + { + radii[3] = MU_B; + break; + } + case (CHEM_POLYMER): + { + for (type = 3; type < RD_TYPES+1; type++) radii[type] = 1.5*MU; + break; + } + case (CHEM_POLYMER_DISS): + { + for (type = 3; type < RD_TYPES+1; type++) radii[type] = 1.5*MU; + break; + } + case (CHEM_POLYMER_STEP): + { + for (type = 2; type < RD_TYPES+1; type++) radii[type] = 1.2*MU; + break; + } + case (CHEM_CATALYTIC_A2D): + { + radii[3] = MU_B; + radii[4] = MU*1.2; + break; + } + case (CHEM_ABCAB): + { + radii[3] = 1.5*MU; + break; + } + case (CHEM_ABCDABC): + { + radii[3] = 1.5*MU; + radii[4] = 1.5*MU; + break; + } + case (CHEM_BZ): + { +// for (type = 2; type <= 8; type++) radii[type] = MU; + for (type = 2; type <= 4; type++) radii[type] = MU; + radii[5] = MU_B; + radii[6] = 1.5*MU_B; + radii[7] = 1.2*radii[3]; + radii[8] = 1.2*radii[5]; + break; + } + case (CHEM_BRUSSELATOR): + { + radii[3] = MU; + radii[4] = 1.2*MU_B; + radii[5] = MU; + radii[6] = MU; + break; + } + } +} + +void adapt_speed_exothermic(t_particle *particle, double delta_ekin) +/* change the particle speed in case of exothermic reaction */ +{ + double old_energy, e_ratio; + + old_energy = (particle->vx*(particle->vx) + particle->vy*(particle->vy))*(particle->mass_inv); + if (old_energy + delta_ekin > 0.0) e_ratio = sqrt((old_energy + delta_ekin)/old_energy); + else e_ratio = 0.0; + particle->vx *= e_ratio; + particle->vy *= e_ratio; +} + +int chem_merge_AAB(int i, int newtype, t_particle particle[NMAXCIRCLES], t_collision *collisions, int ncollisions, double inv_masses[RD_TYPES+1], double radii[RD_TYPES+1]) /* merging particle i with a particle of same type into a particle of type newtype */ /* particular case of chem_merge */ { @@ -5895,10 +6318,10 @@ int chem_merge_AAB(int i, int newtype, t_particle particle[NMAXCIRCLES], t_colli if ((search)&&(particle[k].active)&&(particle[k].type == type1)) { distance = module2(particle[i].deltax[j], particle[i].deltay[j]); - if ((distance < REACION_DIST*MU)&&((double)rand()/RAND_MAX < REACTION_PROB)) + if ((distance < REACTION_DIST*MU)&&((double)rand()/RAND_MAX < REACTION_PROB)) { particle[i].type = newtype; - particle[i].radius *= sqrt(2.0); + particle[i].radius = radii[newtype]; rx = particle[i].xc - particle[k].xc; ry = particle[i].yc - particle[k].yc; particle[i].angle = argument(rx, ry); @@ -5920,7 +6343,7 @@ int chem_merge_AAB(int i, int newtype, t_particle particle[NMAXCIRCLES], t_colli collisions[ncollisions].time = COLLISION_TIME; collisions[ncollisions].color = 0.0; - if (ncollisions < NMAXCIRCLES - 1) ncollisions++; + if (ncollisions < NMAXCOLLISIONS - 1) ncollisions++; else printf("Too many collisions\n"); search = 0; @@ -5930,12 +6353,12 @@ int chem_merge_AAB(int i, int newtype, t_particle particle[NMAXCIRCLES], t_colli return(ncollisions); } -int chem_merge(int i, int type2, int newtype, t_particle particle[NMAXCIRCLES], t_collision *collisions, int ncollisions, double inv_masses[RD_TYPES+1]) +int chem_merge(int i, int type2, int newtype, t_particle particle[NMAXCIRCLES], t_collision *collisions, int ncollisions, double inv_masses[RD_TYPES+1], double radii[RD_TYPES+1]) /* merging of particle i and a particle of type type2 into a particle of type newtype */ { int j, k, type1; short int search = 1; - double distance, rx, ry, m2, mr1, mr2, newmass_inv; + double distance, rx, ry, m2, mr1, mr2, newmass_inv, old_energy, e_ratio; type1 = particle[i].type; newmass_inv = inv_masses[newtype]; @@ -5948,10 +6371,10 @@ int chem_merge(int i, int type2, int newtype, t_particle particle[NMAXCIRCLES], if ((particle[k].active)&&(particle[k].type == type2)) { distance = module2(particle[i].deltax[j], particle[i].deltay[j]); - if ((distance < REACION_DIST*MU)&&((double)rand()/RAND_MAX < REACTION_PROB)) + if ((distance < REACTION_DIST*MU)&&((double)rand()/RAND_MAX < REACTION_PROB)) { particle[i].type = newtype; - particle[i].radius *= 1.5; + particle[i].radius = radii[newtype]; rx = particle[i].xc - particle[k].xc; ry = particle[i].yc - particle[k].yc; particle[i].angle = argument(rx, ry); @@ -5966,6 +6389,8 @@ int chem_merge(int i, int type2, int newtype, t_particle particle[NMAXCIRCLES], particle[i].vy = mr1*particle[i].vy + mr2*particle[k].vy; particle[i].mass_inv = newmass_inv; + if (EXOTHERMIC) adapt_speed_exothermic(&particle[i], DELTA_EKIN); + particle[k].active = 0; collisions[ncollisions].x = particle[i].xc; @@ -5973,7 +6398,7 @@ int chem_merge(int i, int type2, int newtype, t_particle particle[NMAXCIRCLES], collisions[ncollisions].time = COLLISION_TIME; collisions[ncollisions].color = 0.0; - if (ncollisions < 2*NMAXCIRCLES - 1) ncollisions++; + if (ncollisions < 2*NMAXCOLLISIONS - 1) ncollisions++; else printf("Too many collisions\n"); } } @@ -5982,7 +6407,7 @@ int chem_merge(int i, int type2, int newtype, t_particle particle[NMAXCIRCLES], } -int chem_multi_merge(int i, int ni, int type2, int newtype, t_particle particle[NMAXCIRCLES], t_collision *collisions, int ncollisions, double inv_masses[RD_TYPES+1]) +int chem_multi_merge(int i, int ni, int type2, int newtype, t_particle particle[NMAXCIRCLES], t_collision *collisions, int ncollisions, double inv_masses[RD_TYPES+1], double radii[RD_TYPES+1]) /* merging of ni particles of the same type as particle i (including particle i) */ /* and one particle of type type2 into a particle of type newtype */ { @@ -6009,7 +6434,7 @@ int chem_multi_merge(int i, int ni, int type2, int newtype, t_particle particle[ if ((particle[k].active)&&(particle[k].type == type1)) { distance = module2(particle[i].deltax[j], particle[i].deltay[j]); - if ((distance < REACION_DIST*MU)&&((double)rand()/RAND_MAX < REACTION_PROB)) + if ((distance < REACTION_DIST*MU)&&((double)rand()/RAND_MAX < REACTION_PROB)) { k1[n1] = k; n1++; @@ -6018,7 +6443,7 @@ int chem_multi_merge(int i, int ni, int type2, int newtype, t_particle particle[ else if ((particle[k].active)&&(particle[k].type == type2)) { distance = module2(particle[i].deltax[j], particle[i].deltay[j]); - if ((distance < REACION_DIST*MU)&&((double)rand()/RAND_MAX < REACTION_PROB)) + if ((distance < REACTION_DIST*MU)&&((double)rand()/RAND_MAX < REACTION_PROB)) { k2 = k; n2 = 1; @@ -6028,7 +6453,7 @@ int chem_multi_merge(int i, int ni, int type2, int newtype, t_particle particle[ if ((n1 == ni-1)&&(n2 == 1)) { - particle[i].type = newtype; + particle[i].type = radii[newtype]; particle[i].radius *= 1.5; xg = particle[i].xc; yg = particle[i].yc; @@ -6077,7 +6502,7 @@ int chem_multi_merge(int i, int ni, int type2, int newtype, t_particle particle[ collisions[ncollisions].time = COLLISION_TIME; collisions[ncollisions].color = 0.0; - if (ncollisions < 2*NMAXCIRCLES - 1) ncollisions++; + if (ncollisions < 2*NMAXCOLLISIONS - 1) ncollisions++; else printf("Too many collisions\n"); } @@ -6085,7 +6510,7 @@ int chem_multi_merge(int i, int ni, int type2, int newtype, t_particle particle[ } -int chem_catalytic_merge(int i, int type_catalyst, int newtype, t_particle particle[NMAXCIRCLES], t_collision *collisions, int ncollisions, double inv_masses[RD_TYPES+1]) +int chem_catalytic_merge(int i, int type_catalyst, int newtype, t_particle particle[NMAXCIRCLES], t_collision *collisions, int ncollisions, double inv_masses[RD_TYPES+1], double radii[RD_TYPES+1]) /* merging of 2 particles of the same type as particle i (including particle i) */ /* into a particle of type newtype, provided a particle of type type_catalyst is present */ { @@ -6105,7 +6530,7 @@ int chem_catalytic_merge(int i, int type_catalyst, int newtype, t_particle parti if ((particle[k].active)&&(particle[k].type == type1)) { distance = module2(particle[i].deltax[j], particle[i].deltay[j]); - if ((distance < REACION_DIST*MU)&&((double)rand()/RAND_MAX < REACTION_PROB)) + if ((distance < REACTION_DIST*MU)&&((double)rand()/RAND_MAX < REACTION_PROB)) { k1 = k; n1 = 1; @@ -6114,7 +6539,7 @@ int chem_catalytic_merge(int i, int type_catalyst, int newtype, t_particle parti else if ((particle[k].active)&&(particle[k].type == type_catalyst)) { distance = module2(particle[i].deltax[j], particle[i].deltay[j]); - if ((distance < REACION_DIST*MU)&&((double)rand()/RAND_MAX < REACTION_PROB)) + if ((distance < REACTION_DIST*MU)&&((double)rand()/RAND_MAX < REACTION_PROB)) { k2 = k; n2 = 1; @@ -6125,7 +6550,7 @@ int chem_catalytic_merge(int i, int type_catalyst, int newtype, t_particle parti if ((n1 == 1)&&(n2 == 1)) { particle[i].type = newtype; - particle[i].radius *= 1.2; + particle[i].radius = radii[newtype]; xg = 0.5*(particle[i].xc + particle[k1].xc + particle[k2].xc); yg = 0.5*(particle[i].yc + particle[k1].yc + particle[k2].yc); rx = particle[i].xc - xg; @@ -6151,7 +6576,7 @@ int chem_catalytic_merge(int i, int type_catalyst, int newtype, t_particle parti collisions[ncollisions].time = COLLISION_TIME; collisions[ncollisions].color = 0.0; - if (ncollisions < 2*NMAXCIRCLES - 1) ncollisions++; + if (ncollisions < 2*NMAXCOLLISIONS - 1) ncollisions++; else printf("Too many collisions\n"); } @@ -6159,7 +6584,7 @@ int chem_catalytic_merge(int i, int type_catalyst, int newtype, t_particle parti } -int chem_split(int i, int newtype1, int newtype2, t_particle particle[NMAXCIRCLES], t_collision *collisions, int ncollisions, double inv_masses[RD_TYPES+1]) +int chem_split(int i, int newtype1, int newtype2, t_particle particle[NMAXCIRCLES], t_collision *collisions, int ncollisions, double inv_masses[RD_TYPES+1], double radii[RD_TYPES+1]) /* split particle i into particles of type newtype1 and newtype2 */ { int j, k, oldtype; @@ -6182,16 +6607,24 @@ int chem_split(int i, int newtype1, int newtype2, t_particle particle[NMAXCIRCLE { particle[i].type = newtype2; particle[i].mass_inv = inv_masses[newtype2]; - particle[i].radius = MU; + particle[i].radius = radii[newtype2]; - particle[ncircles-1].radius = MU; + particle[ncircles-1].type = newtype1; + particle[ncircles-1].mass_inv = inv_masses[newtype1]; + particle[ncircles-1].radius = radii[newtype1]; + if (EXOTHERMIC) + { + adapt_speed_exothermic(&particle[i], -0.5*DELTA_EKIN); + adapt_speed_exothermic(&particle[ncircles-1], -0.5*DELTA_EKIN); + } + collisions[ncollisions].x = particle[i].xc; collisions[ncollisions].y = particle[i].yc; collisions[ncollisions].time = COLLISION_TIME; collisions[ncollisions].color = type_hue(oldtype); - if (ncollisions < NMAXCIRCLES - 1) ncollisions++; + if (ncollisions < NMAXCOLLISIONS - 1) ncollisions++; else printf("Too many collisions\n"); } else @@ -6203,20 +6636,160 @@ int chem_split(int i, int newtype1, int newtype2, t_particle particle[NMAXCIRCLE return(ncollisions); } -int update_types(t_particle particle[NMAXCIRCLES], t_collision *collisions, int ncollisions) +int chem_transfer(int i, int type2, int newtype1, int newtype2, t_particle particle[NMAXCIRCLES], t_collision *collisions, int ncollisions, double inv_masses[RD_TYPES+1], double radii[RD_TYPES+1], double reac_dist) +/* reaction of particle i and a particle of type type2 into a particles of types newtype1 and newtype2 */ +{ + int j, k, type1; + short int search = 1; + double distance, rx, ry, mass_inv1, mass_inv2, newmass_inv1, newmass_inv2, old_energy, e_ratio, omega; + + type1 = particle[i].type; + mass_inv1 = inv_masses[type1]; + mass_inv2 = inv_masses[type2]; + newmass_inv1 = inv_masses[newtype1]; + newmass_inv2 = inv_masses[newtype2]; + + for (j=0; j 2)&&((double)rand()/RAND_MAX < DISSOCIATION_PROB)) - ncollisions = chem_split(i, 1, particle[i].type - 1, particle, collisions, ncollisions, inv_masses); + ncollisions = chem_split(i, 1, particle[i].type - 1, particle, collisions, ncollisions, inv_masses, radii); } printf("%i collisions\n", ncollisions); return(ncollisions); } + case (CHEM_POLYMER_STEP): + { + for (i=0; i 1)&&((double)rand()/RAND_MAX < DISSOCIATION_PROB)) + { + k = rand()%(type-1) + 1; +// printf("Splitting type %i into type %i and type %i\n", type, k, type - k); + ncollisions = chem_split(i, k, type - k, particle, collisions, ncollisions, inv_masses, radii); + } + } + + printf("%i collisions\n", ncollisions); + return(ncollisions); + } + case (CHEM_CATALYTIC_A2D): + { + for (i=0; i X */ + { + ncollisions = chem_convert(i, 3, particle, collisions, ncollisions, inv_masses, radii, DISSOCIATION_PROB); + break; + } + case (2): /* B + X -> Y + D */ + { + ncollisions = chem_transfer(i, 3, 4, 5, particle, collisions, ncollisions, inv_masses, radii, REACTION_DIST); + break; + } + case (3): + { + /* X -> B, modified from Brusselator */ + ncollisions = chem_convert(i, 2, particle, collisions, ncollisions, inv_masses, radii, 0.5*DISSOCIATION_PROB); + /* 2X + Y -> 3X */ + ncollisions = chem_catalytic_convert(i, 4, 3, particle, collisions, ncollisions, inv_masses, radii, 3.0*REACTION_DIST, 1.0); + break; + } + case (5): /* D -> A or B if concentration of X or Y is small */ + { + n = time*(RD_TYPES+1); + n3 = particle_numbers[n+3]; + n4 = particle_numbers[n+4]; + p1 = 1.0/((double)(n3*n3/10+1)); + p2 = 5.0*DISSOCIATION_PROB + 1.0/((double)(n4*n4/200+1)); + rnd = (double)rand()/RAND_MAX; + if (rnd < p1/(p1+p2)) + ncollisions = chem_convert(i, 1, particle, collisions, ncollisions, inv_masses, radii, p1); + else /*if (n4 < 1000)*/ + ncollisions = chem_convert(i, 2, particle, collisions, ncollisions, inv_masses, radii, p2); + } + } + } + return(ncollisions); + } } } @@ -6514,6 +7272,13 @@ void draw_trajectory_plot(t_group_data *group_speeds, int i) write_text_fixedwidth(plotxmax - 0.1, plotymin - 0.1, message); } + +void write_size_text(double x, double y, char *message, int largewin) +{ + if (largewin) write_text(x, y, message); + else write_text_fixedwidth(x, y, message); +} + void draw_particle_nb_plot(int *particle_numbers, int i) /* draw plot of number of particles as a function of time */ { @@ -6521,7 +7286,7 @@ void draw_particle_nb_plot(int *particle_numbers, int i) char message[100]; static double xmin, xmax, ymin, ymax, xmid, ymid, dx, dy, plotxmin, plotxmax, plotymin, plotymax; double pos[2], x1, y1, x2, y2, rgb[3]; - static int first = 1, gshift = INITIAL_TIME + NSTEPS, nmax, ygrad; + static int first = 1, gshift = INITIAL_TIME + NSTEPS, nmax, ygrad, largewin; if (first) { @@ -6550,6 +7315,8 @@ void draw_particle_nb_plot(int *particle_numbers, int i) power = (int)(log((double)nmax)/log(10.0)) - 1.0; ygrad = (int)ipow(10.0, power); + largewin = (WINWIDTH > 1280); + first = 0; } @@ -6563,19 +7330,23 @@ void draw_particle_nb_plot(int *particle_numbers, int i) set_type_color(type, 1.0, rgb); x1 = plotxmin; y1 = plotymin; + glBegin(GL_LINE_STRIP); for (j=1; j= 1000) write_text_fixedwidth(plotxmin - 0.17, y1 - 0.015, message); - else write_text_fixedwidth(plotxmin - 0.14, y1 - 0.015, message); +// if (j*ygrad >= 1000) write_size_text(plotxmin - 0.15, y1 - 0.015, message, largewin); +// else write_size_text(plotxmin - 0.12, y1 - 0.015, message, largewin); + if (j*ygrad >= 1000) write_text_fixedwidth(plotxmin - 0.15, y1 - 0.015, message); + else write_text_fixedwidth(plotxmin - 0.12, y1 - 0.015, message); } j++; } sprintf(message, "time"); +// write_size_text(plotxmax - 0.1, plotymin - 0.05, message, largewin); write_text_fixedwidth(plotxmax - 0.1, plotymin - 0.05, message); } diff --git a/sub_maze.c b/sub_maze.c index 8c0e755..3f27a62 100644 --- a/sub_maze.c +++ b/sub_maze.c @@ -4,6 +4,11 @@ /* Change constant RAND_SHIFT to change the maze */ +#define MAZE_TYPE_SQUARE 0 /* maze with square cells */ +#define MAZE_TYPE_CIRCLE 1 /* circular maze */ +#define MAZE_TYPE_HEX 2 /* honeycomb maze */ +#define MAZE_TYPE_OCT 3 /* maze with octagonal and square cells */ + typedef struct { short int nneighb; /* number of neighbours */ @@ -488,6 +493,208 @@ void init_hex_maze_graph(t_maze maze[NXMAZE*NYMAZE]) } } +void init_oct_maze_graph(t_maze maze[NXMAZE*NYMAZE]) +/* initialise graph of maze made of octagons and squares */ +{ + int i, j, k, n, p, q; + + printf("Initializing maze\n"); + if (MAZE_MAX_NGBH < 8) + { + printf("Error: MAZE_MAX_NGBH should be at least 8 for circular maze\n"); + exit(0); + } + + /* initialize neighbours */ + /* in the bulk */ + for (i=1; i 0) printf("Cell (%i, %i)\n", i, j); +// for (k = 0; k = COLOR_HUEMAX) hue -= delta_hue; + + hsl_to_rgb(hue, r, lum, rgb); + /* saturation = r */ +} + void blank() { double rgb[3]; @@ -403,6 +421,54 @@ void draw_colored_circle(double x, double y, double r, int nseg, double rgb[3]) glEnd(); } +void draw_colored_octahedron(double x, double y, double r, double rgb[3]) +{ + int i, ij[2]; + double pos[2], alpha, dalpha; + + dalpha = DPI*0.125; + + glColor3f(rgb[0], rgb[1], rgb[2]); + glBegin(GL_TRIANGLE_FAN); + glVertex2d(x,y); + for (i=0; i<=8; i++) + { + alpha = ((double)i+0.5)*dalpha; + glVertex2d(x + r*cos(alpha), y + r*sin(alpha)); + } + + glEnd(); +} + +void draw_colored_rectangle_rgb(double x1, double y1, double x2, double y2, double rgb[3]) +{ + glColor3f(rgb[0], rgb[1], rgb[2]); + + glBegin(GL_QUADS); + glVertex2d(x1, y1); + glVertex2d(x2, y1); + glVertex2d(x2, y2); + glVertex2d(x1, y2); + glEnd(); + + glColor3f(1.0, 1.0, 1.0); + glBegin(GL_LINE_LOOP); + glVertex2d(x1, y1); + glVertex2d(x2, y1); + glVertex2d(x2, y2); + glVertex2d(x1, y2); + glEnd(); +} + +void draw_colored_rectangle(double x1, double y1, double x2, double y2, double hue) +{ + double rgb[3]; + + hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE); + draw_colored_rectangle_rgb(x1, y1, x2, y2, rgb); +} + + void draw_filled_sector(double x, double y, double rmin, double rmax, double phi1, double dphi, int nseg, double rgb[3]) { @@ -5765,7 +5831,7 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ if (POLYLINE_PATTERN == P_MAZE) return ((vabs(x) < 1.1*XMAX)&&(vabs(y) < 1.1*YMAX)); else if ((POLYLINE_PATTERN == P_MAZE_DIAG)||(POLYLINE_PATTERN == P_MAZE_RANDOM)) return ((vabs(x) < 1.1*XMAX)&&(vabs(y) < 1.1*YMAX)); - else if ((POLYLINE_PATTERN == P_MAZE_CIRCULAR)||(POLYLINE_PATTERN == P_MAZE_HEX)) + else if ((POLYLINE_PATTERN == P_MAZE_CIRCULAR)||(POLYLINE_PATTERN == P_MAZE_HEX)||(POLYLINE_PATTERN == P_MAZE_OCT)) return ((vabs(x) < 1.1*XMAX)&&(vabs(y) < 1.1*YMAX)); else return(1); break; @@ -6309,8 +6375,16 @@ void init_polyline_circular_maze(t_segment* polyline, t_circle* circles, t_arc* arcs[narcs].xc = MAZE_XSHIFT; arcs[narcs].yc = 0.0; arcs[narcs].radius = rmax; - arcs[narcs].angle1 = dphi; - arcs[narcs].dangle = DPI - dphi; + if (CLOSE_MAZE) + { + arcs[narcs].angle1 = 0.0; + arcs[narcs].dangle = DPI; + } + else + { + arcs[narcs].angle1 = dphi; + arcs[narcs].dangle = DPI - dphi; + } narcs++; } @@ -6500,7 +6574,191 @@ void init_polyline_hex_maze(t_segment* polyline, t_circle* circles, t_maze* maze } } -// (()||())&& +void init_polyline_oct_maze(t_segment* polyline, t_circle* circles, t_maze* maze) +{ + int i, j, n; + double a, b, a2, c, dx, x1, y1, padding = 0.02; + + dx = (YMAX - YMIN - 2.0*padding)/((double)NYMAZE+0.5); + a = dx*(2.0 - sqrt(2.0)); + b = a/sqrt(2.0); + a2 = 0.5*a; + c = a2 + b; + + for (i = 0; i < NXMAZE; i++) + for (j = 0; j < NYMAZE; j++) + { + n = nmaze(i, j); + x1 = YMIN + padding + ((double)i + 0.5)*dx + MAZE_XSHIFT; + y1 = YMIN + padding + ((double)j + 0.75)*dx; + + if ((i+j)%2 == 0) + { + if (((i>0)||(j!=NYMAZE/2))&&(maze[n].west)) + { + polyline[nsides].x1 = x1 - c; + polyline[nsides].y1 = y1 - a2; + polyline[nsides].x2 = x1 - c; + polyline[nsides].y2 = y1 + a2; + polyline[nsides].angle = PID; + nsides++; + } + + if (maze[n].south) + { + polyline[nsides].x1 = x1 - a2; + polyline[nsides].y1 = y1 - c; + polyline[nsides].x2 = x1 + a2; + polyline[nsides].y2 = y1 - c; + polyline[nsides].angle = 0.0; + nsides++; + } + + if (maze[n].southwest) + { + polyline[nsides].x1 = x1 - a2; + polyline[nsides].y1 = y1 - c; + polyline[nsides].x2 = x1 - c; + polyline[nsides].y2 = y1 - a2; + polyline[nsides].angle = 1.5*PID; + nsides++; + } + + if (maze[n].northwest) + { + polyline[nsides].x1 = x1 - c; + polyline[nsides].y1 = y1 + a2; + polyline[nsides].x2 = x1 - a2; + polyline[nsides].y2 = y1 + c; + polyline[nsides].angle = 0.5*PID; + nsides++; + } + } + else + { + if (((i>0)||(j!=NYMAZE/2))&&(maze[n].west)) + { + polyline[nsides].x1 = x1 - a2; + polyline[nsides].y1 = y1 - a2; + polyline[nsides].x2 = x1 - a2; + polyline[nsides].y2 = y1 + a2; + polyline[nsides].angle = PID; + nsides++; + } + + if (maze[n].south) + { + polyline[nsides].x1 = x1 - a2; + polyline[nsides].y1 = y1 - a2; + polyline[nsides].x2 = x1 + a2; + polyline[nsides].y2 = y1 - a2; + polyline[nsides].angle = 0.0; + nsides++; + } + + + } + } + + /* bottom of maze */ + for (i=0; i0)||(j!=NYMAZE/2))&&(maze[n].west)) + if (((i>0)||(j!=NYMAZE/2)||(CLOSE_MAZE))&&(maze[n].west)) { polyline[nsides].x1 = x1; polyline[nsides].y1 = y1; @@ -7019,7 +7277,8 @@ void init_polyline(t_segment polyline[NMAXPOLY], t_circle circles[NMAXCIRCLES], polyline[nsides].x1 = x1; polyline[nsides].y1 = YMIN - 1000.0; polyline[nsides].x2 = x1; - polyline[nsides].y2 = y1 - dy; + if (CLOSE_MAZE) polyline[nsides].y2 = y1; + else polyline[nsides].y2 = y1 - dy; polyline[nsides].angle = PID; nsides++; @@ -7488,7 +7747,21 @@ void init_polyline(t_segment polyline[NMAXPOLY], t_circle circles[NMAXCIRCLES], break; } - + case (P_MAZE_OCT): + { + maze = (t_maze *)malloc(NXMAZE*NYMAZE*sizeof(t_maze)); + + init_oct_maze(maze); + + nsides = 0; + ncircles = 0; + + init_polyline_oct_maze(polyline, circles, maze); + + free(maze); + + break; + } } } @@ -7662,6 +7935,7 @@ int maze_type(int polyline_pattern) case (P_MAZE_CIRCULAR): return(1); case (P_MAZE_CIRC_SCATTERER): return(1); case (P_MAZE_HEX): return(2); + case (P_MAZE_OCT): return(3); default: return(0); } @@ -7671,8 +7945,8 @@ int find_maze_cell(double x, double y) /* return maze cell number for coordinates (x, y) */ { int i, j, n, block, q, w, k, l; - double padding = 0.02, x1, y1, u, v, r, phi, phi1, angle1, dphi, r1, tolerance = 0.025; - static double dx, dy, rmin, rmax, angle, rr, rrtol, h, dr, dphi_table[5], x0, y0; + double padding = 0.02, x1, y1, u, v, r, phi, phi1, angle1, dphi, r1, tolerance = 0.025, x2, y2; + static double dx, dy, rmin, rmax, angle, rr, rrtol, h, dr, dphi_table[5], x0, y0, a, b, d, a2, a2tol; static int first = 1, nblocks; if (first) @@ -7708,6 +7982,22 @@ int find_maze_cell(double x, double y) x0 = YMIN + padding + MAZE_XSHIFT; y0 = YMIN + padding + h; + break; + } + case (3): /* octahedron-square maze */ + { + dx = (YMAX - YMIN - 2.0*padding)/((double)NYMAZE+0.5); + a = dx*(2.0 - sqrt(2.0)); + a2 = 0.5*a; + b = a/sqrt(2.0); + d = 2.0*dx; + + x0 = YMIN + padding + MAZE_XSHIFT - 0.5*b; + y0 = YMIN + padding + 0.25*dx - 0.5*b; + rr = sqrt((b+a2)*(b+a2) + a2*a2); + rrtol = rr*(1.0 - tolerance); + a2tol = a2*(1.0 - tolerance); + break; } } @@ -7826,6 +8116,69 @@ int find_maze_cell(double x, double y) } return(-10); } + case (3): /* octahedron-square maze */ + { + i = 2*(int)((x-x0)/d); + j = 2*(int)((y-y0)/d); + + x1 = x - x0 - (double)(i)*dx; + y1 = y - y0 - (double)(j)*dx; + +// printf("(x,y) = (%.3lg,%.3lg), (i,j) = (%i,%i), (x1,y1) = (%.3lg,%.3lg)\n", x, y, i, j, x1/d, y1/d); + + if (i < 0) return(-1); + if (i > NXMAZE-1) return(-1); + if (j < 0) return(-1); + if (j > NYMAZE-1) return(-1); + + if (x1 < b) + { + if (x1 + y1 < b) i--; + else if (y1 - x1 > dx) i--; + } + else if ((x1 > dx)&&(x1 < a + 2.0*b)) + { + if (y1 - x1 < - dx) i++; + else if (x1 + y1 > 3.0*a + 2.0*b) i++; + } + else if (x1 >= a + 2.0*b) i++; + + if (y1 < b) + { + if (x1 + y1 < b) j--; + else if (x1 - y1 > a + b) j--; + } + else if ((y1 > a + b)&&(y1 < a + 2.0*b)) + { + if (x1 - y1 < - dx) j++; + else if (x1 + y1 > 3.0*a + 2.0*b) j++; + } + else if (y1 >= a + 2.0*b) j++; + + if (i < 0) return(-1); + if (i > NXMAZE-1) return(-1); + if (j < 0) return(-1); + if (j > NYMAZE-1) return(-1); + + + /* avoid finding wrong cell for particles close to wall */ + x2 = x - x0 - (double)i*dx - b - a2; + y2 = y - y0 - (double)j*dx - b - a2; + if ((i+j)%2 == 0) + { + if (!in_polygon(x2, y2, rrtol, 8, 0.25)) return(-10); + } + else + { + if (vabs(x2) > a2tol) return(-10); + if (vabs(y2) > a2tol) return(-10); + } + + n = nmaze(i, j); + + if (n < NXMAZE*NYMAZE) return(n); + else return(-1); + } } } @@ -7834,7 +8187,7 @@ void draw_maze_cell(int n, int part_number, double minprop) { int i, j, block, q; double x, y, padding = 0.02, rgb[3], w; - static double dx, dy, rmin, rmax, angle, r, r1, dr, phi1, dphi, h, x0, y0; + static double dx, dy, rmin, rmax, angle, r, r1, dr, phi1, dphi, h, x0, y0, a, b, a2, square_prop; static int first = 1, nblocks; if (first) switch (maze_type(POLYLINE_PATTERN)) { @@ -7865,6 +8218,21 @@ void draw_maze_cell(int n, int part_number, double minprop) y0 = YMIN + padding + h; break; } + case (3): /* octahedron-square maze */ + { + dx = (YMAX - YMIN - 2.0*padding)/((double)NYMAZE+0.5); + a = dx*(2.0 - sqrt(2.0)); + b = a/sqrt(2.0); + a2 = 0.5*a; + + r = sqrt((b+a2)*(b+a2) + a2*a2); + square_prop = 2.0*(sqrt(2.0) + 1.0); + + x0 = YMIN + padding + MAZE_XSHIFT; + y0 = YMIN + padding + 0.25*dx; + + break; + } first = 0; } @@ -7925,6 +8293,22 @@ void draw_maze_cell(int n, int part_number, double minprop) draw_colored_circle(x, y, r, 6, rgb); break; } + case (3): /* octahedron-square maze */ + { + i = n%NXMAZE; + j = n/NXMAZE; + + x = x0 + ((double)i + 0.5)*dx; + y = y0 + ((double)j + 0.5)*dx; + + /* adapt particle number in square to have constant density */ + if ((i+j)%2 == 1) part_number = (int)(square_prop*(double)part_number); + + rgb_color_scheme_density(part_number, rgb, minprop); + if ((i+j)%2 == 0) draw_colored_octahedron(x, y, r,rgb); + else erase_rectangle(x-a2, y-a2, x+a2, y+a2, rgb); + break; + } } } diff --git a/sub_rde.c b/sub_rde.c index 16dc790..dc1e7d0 100644 --- a/sub_rde.c +++ b/sub_rde.c @@ -6,32 +6,32 @@ double viscosity; /* viscosity (parameter in front of Laplacian) */ double rpslzb; /* second parameter in Rock-Paper-Scissors-Lizard-Spock equation */ double flow_speed; /* flow speed for laminar boundary conditions in Euler equation */ -double gaussian() -/* returns standard normal random variable, using Box-Mueller algorithm */ -{ - static double V1, V2, S; - static int phase = 0; - double X; - - if (phase == 0) - { - do - { - double U1 = (double)rand() / RAND_MAX; - double U2 = (double)rand() / RAND_MAX; - V1 = 2 * U1 - 1; - V2 = 2 * U2 - 1; - S = V1 * V1 + V2 * V2; - } - while(S >= 1 || S == 0); - X = V1 * sqrt(-2 * log(S) / S); - } - else X = V2 * sqrt(-2 * log(S) / S); - - phase = 1 - phase; - - return X; -} +// double gaussian() +// /* returns standard normal random variable, using Box-Mueller algorithm */ +// { +// static double V1, V2, S; +// static int phase = 0; +// double X; +// +// if (phase == 0) +// { +// do +// { +// double U1 = (double)rand() / RAND_MAX; +// double U2 = (double)rand() / RAND_MAX; +// V1 = 2 * U1 - 1; +// V2 = 2 * U2 - 1; +// S = V1 * V1 + V2 * V2; +// } +// while(S >= 1 || S == 0); +// X = V1 * sqrt(-2 * log(S) / S); +// } +// else X = V2 * sqrt(-2 * log(S) / S); +// +// phase = 1 - phase; +// +// return X; +// } void init_random(double mean, double amplitude, double *phi[NFIELDS], short int xy_in[NX*NY]) /* initialise field with gaussian at position (x,y) */ @@ -510,28 +510,63 @@ void init_shear_flow(double amp, double delta, double rho, int nx, int ny, doubl } } -void set_boundary_laminar_flow(double amp, double modulation, double period, double yshift, double *phi[NFIELDS], short int xy_in[NX*NY], int imin, int imax, int jmin, int jmax) +void set_boundary_laminar_flow(double amp, double xmodulation, double ymodulation, double xperiod, double yperiod, double yshift, double density_mod, double *phi[NFIELDS], short int xy_in[NX*NY], int imin, int imax, int jmin, int jmax, double factor) /* enfoce laminar flow in x direction on top and bottom boundaries */ /* phi[0] is stream function, phi[1] is vorticity */ /* amp is global amplitude */ { int i, j; - double xy[2], y1, a, b, f, cplus, cminus; + double xy[2], y1, a, b, f, cplus, cminus, comp_factor; - a = period*PI/YMAX; + a = xperiod*PI/YMAX; + b = yperiod*PI/XMAX; + comp_factor = 1.0 - factor; - for (i=imin; i= 0)&&(xp <= length)) return(vabs(yp)); + else if (xp < 0) return(module2(xp, yp)); + else return(module2(xp-length, yp)); +} + + +double tesla_distance(double x, double y, double a, double l, double theta) +/* distance to center of Tesla valve */ +{ + double dmin, dist, ct, st, tt, b, c, d, l1, l2, angle; + double xa, ya, xb, yb, xc, yc, xd, yd, xe, ye; + + ct = cos(theta); + st = sin(theta); + tt = st/ct; + + b = a*ct; +// c = (l*st - a)*ct; + d = 0.5*a*tt; + + l1 = l*cos(2.0*theta); +// l2 = l - a/st; + l2 = 0.3*l; /* TODO */ + + /* upper segment */ + xa = l1*ct + 0.5*b*st; + ya = a + l1*st - 0.5*b*ct; + dmin = distance_to_segment(x, y, -d, 0.0, xa, ya); + + /* lower segment */ + xb = l*ct; + yb = -l*st - 0.5*a; + dist = distance_to_segment(x, y, -d, 0.0, xb, yb); + if (dist < dmin) dmin = dist; + + /* small segment */ +// xc = xb; +// yc = -l*st + 0.5*a; +// dist = distance_to_segment(x, y, xc - 0.5*a/tt, -l*st, xc, yc); +// if (dist < dmin) dmin = dist; + + /* middle segment */ + xd = l*ct - 1.0*a*st; + yd = -l*st + a + 1.0*a*ct; + dist = distance_to_segment(x, y, xd - l2*ct, yd - l2*st, xd, yd); + if (dist < dmin) dmin = dist; + + /* circular part */ + xe = 0.5*(xa + xd); + ye = 0.5*(ya + yd); + c = module2(xd - xe, yd - ye) - 0.5*b; + dist = vabs(module2(x - xe, y - ye) - c - 0.5*b*ct); + angle = argument(x - xe, y - ye); + angle += PID - theta; + if (angle > DPI) angle -= DPI; + if ((angle > 0.0)&&(angle < PI)&&(dist < dmin)) dmin = dist; + + return(dmin); +} + +void initialize_bcfield(double bc_field[NX*NY], t_rectangle polyrect[NMAXPOLY]) /* apply smooth modulation to adapt initial state to obstacles */ { - int i, j; - double xy[2], r, f; + int i, j, nsides, s; + double xy[2], x, y, r, f, a, l, theta, x1, x2, y1, y2, distance, d, d0, length, height, mid, fmin, ct, st; switch (OBSTACLE_GEOMETRY) { case (D_SINAI): @@ -615,6 +724,171 @@ void initialize_bcfield(double bc_field[NX*NY]) } break; } + case (D_EXT_ELLIPSE): + { + for (i=0; i l*ct) + { + d = vabs(xy[1]); + if (d < d0) d0 = d; + } + + r = a - d0; + + f = 0.5*(1.0 + tanh(BC_STIFFNESS*r)); + bc_field[i*NY+j] = f; + } + break; + } + } + + if ((OBSTACLE_GEOMETRY == D_MAZE)||(OBSTACLE_GEOMETRY == D_MAZE_CHANNELS)) + { + d0 = 0.25*MAZE_WIDTH; + fmin = 0.5*(1.0 - tanh(d0*BC_STIFFNESS)); + for (i=0; i 4.0*MAZE_WIDTH)&&(height > 4.0*MAZE_WIDTH)) + { + if (x < x1) + { + if (y < y1) d = module2(x - x1, y - y1); + else if (y > y2) d = module2(x - x1, y - y2); + else d = x1 - x; + } + else if (x > x2) + { + if (y < y1) d = module2(x - x2, y - y1); + else if (y > y2) d = module2(x - x2, y - y2); + else d = x - x2; + } + else + { + if (y < y1) d = y1 - y; + else if (y > y2) d = y - y2; + else d = 0.0; + } + } + else if (length > height) + { + mid = 0.5*(polyrect[s].y1 + polyrect[s].y2); + if ((x > x1)&&(x < x2)) d = vabs(y - mid); + else if (x <= x1) d = module2(x - x1, y - mid); + else d = module2(x - x2, y - mid); + } + else + { + mid = 0.5*(polyrect[s].x1 + polyrect[s].x2); + if ((y > y1)&&(y < y2)) d = vabs(x - mid); + else if (y <= y1) d = module2(x - mid, y - y1); + else d = module2(x - mid, y - y2); + } + if (d < distance) distance = d; + } + if (distance < d0) f = fmin*distance/d0; + else f = 0.5*(1.0 + tanh(BC_STIFFNESS*(distance - 1.25*MAZE_WIDTH))); + bc_field[i*NY+j] = f; +// printf("distance = %.5lg, bcfield = %.5lg\n", distance, f); + } + } } } @@ -628,10 +902,10 @@ void adapt_state_to_bc(double *phi[NFIELDS], double bc_field[NX*NY], short int x #pragma omp parallel for private(i,j) for (i=0; i DPI) value -= DPI; + rde[i*NY+j].field_arg = value; + } +} + void compute_vorticity(t_rde rde[NX*NY]) /* compute the log of a field */ { @@ -1255,7 +1546,7 @@ void compute_vorticity(t_rde rde[NX*NY]) for (i=0; i= 1 || S == 0); + X = V1 * sqrt(-2 * log(S) / S); + } + else X = V2 * sqrt(-2 * log(S) / S); + + phase = 1 - phase; + + return X; +} + + + + // int in_polygon(double x, double y, double r, int npoly, double apoly) // /* test whether (x,y) is in regular polygon of npoly sides inscribed in circle of radious r, turned by apoly Pi/2 */ // { @@ -1635,7 +1665,7 @@ int compute_maze_coordinates(t_rectangle polyrect[NMAXPOLY], int type) { t_maze* maze; int i, j, n, nsides = 0, ropening; - double dx, dy, x1, y1, x0, padding = 0.02, pos[2], width = 0.02; + double dx, dy, x1, y1, x0, padding = 0.02, pos[2], width = MAZE_WIDTH; maze = (t_maze *)malloc(NXMAZE*NYMAZE*sizeof(t_maze)); @@ -1756,7 +1786,7 @@ int compute_maze_coordinates(t_rectangle polyrect[NMAXPOLY], int type) polyrect[nsides].x2 = XMAX + padding; polyrect[nsides].y2 = YMAX + padding; nsides++; - + /* left channel */ x1 = YMIN + padding + MAZE_XSHIFT; polyrect[nsides].x1 = XMIN - padding; @@ -2057,6 +2087,26 @@ int init_polyrect(t_rectangle polyrect[NMAXPOLY]) } } +int init_polyrect_euler(t_rectangle polyrect[NMAXPOLY], int domain) +/* initialise variable polyrect, for certain polygonal domain shapes */ +{ + switch (domain) { + case (D_MAZE): + { + return(compute_maze_coordinates(polyrect, 0)); + } + case (D_MAZE_CLOSED): + { + return(compute_maze_coordinates(polyrect, 1)); + } + case (D_MAZE_CHANNELS): + { + return(compute_maze_coordinates(polyrect, 2)); + } + } +} + + void init_polyrect_arc(t_rect_rotated polyrectrot[NMAXPOLY], t_arc polyarc[NMAXPOLY], int *npolyrect, int *npolyarc) /* initialise variables polyrectrot and polyarc, for certain domain shapes */ { @@ -4744,7 +4794,7 @@ void draw_color_scheme_palette(double x1, double y1, double x2, double y2, int p } case (P_LOG_ENERGY): { - value = LOG_SHIFT + LOG_SCALE*log(dy_e*(double)(j - jmin)*100.0/E_SCALE); + value = LOG_SCALE*log(dy_e*(double)(j - jmin)*100.0/E_SCALE); // if (value <= 0.0) value = 0.0; color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); break; @@ -4874,7 +4924,7 @@ void draw_color_scheme_palette_fade(double x1, double y1, double x2, double y2, } case (P_LOG_ENERGY): { - value = LOG_SHIFT + LOG_SCALE*log(dy_e*(double)(j - jmin)*100.0/E_SCALE); + value = LOG_SCALE*log(dy_e*(double)(j - jmin)*100.0/E_SCALE); // printf("value = %.2lg\n", value); // if (value <= 0.0) value = 0.0; color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); @@ -5475,3 +5525,589 @@ double oscillating_bc(int time) } } +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; + 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; + + rgb[0] = 1.0; + rgb[1] = 1.0; + rgb[2] = 1.0; + + if (VARIABLE_IOR) + { + switch (IOR) { + case (IOR_MANDELBROT): + { + #pragma omp parallel for private(i,j) + for (i=0; i= MANDELLEVEL) + { +// tc[i*NY+j] = COURANT; + tcc_table[i][j] = courant2; +// tgamma[i*NY+j] = GAMMA; + } + else + { + speed = (double)k/(double)MANDELLEVEL; + if (speed < 1.0e-10) speed = 1.0e-10; + else if (speed > 10.0) speed = 10.0; + tcc_table[i][j] = courant2*speed; +// tc[i*NY+j] = COURANT*sqrt(speed); +// tgamma[i*NY+j] = GAMMA; + } + } + } + break; + } + case (IOR_EARTH): + { + for (i=0; i 1.0) c = 0.0; + else if (r2 < 0.25*0.25) c = 0.8*COURANT; + else if (r2 < 0.58*0.58) c = COURANT*(0.68 - 0.55*r2); + else c = COURANT*(1.3 - 0.9*r2); +// tc[i*NY+j] = c; + tcc_table[i][j] = c; +// tgamma[i*NY+j] = GAMMA; + } + } + break; + } + case (IOR_EXPLO_LENSING): + { + salpha = DPI/(double)NPOLY; +// lambda1 = LAMBDA; +// mu1 = LAMBDA; + lambda1 = 0.5*LAMBDA; + mu1 = 0.5*LAMBDA; + h = lambda1*tan(PI/(double)NPOLY); + if (h < mu1) ll = sqrt(mu1*mu1 - h*h); + else ll = 0.0; + +// #pragma omp parallel for private(i,j) + for (i=0; i 1.0) height[n] = 1.0; + if (height[n] < 0.0) height[n] = 0.0; + } + +// #pragma omp parallel for private(i,j) + for (i=0; i rmin2)) + { + r2 = (double)(irad2); + wave_height1 = wave_height*exp(-r2/variance); + phi[packet[i].ix + j][packet[i].iy + k] = wave_height1; + phi[packet[i].ix - j][packet[i].iy + k] = wave_height1; + phi[packet[i].ix + j][packet[i].iy - k] = wave_height1; + phi[packet[i].ix - j][packet[i].iy - k] = wave_height1; + } + else if (irad2 <= rmin2) + { + phi[packet[i].ix + j][packet[i].iy + k] = 0.0; + phi[packet[i].ix - j][packet[i].iy + k] = 0.0; + phi[packet[i].ix + j][packet[i].iy - k] = 0.0; + phi[packet[i].ix - j][packet[i].iy - k] = 0.0; + } + + } +// printf("Adding wave packet %i with value %.3lg\n", i, wave_height); + +// if (add==0) + { +// t -= 1.0/(double)NVID; + t -= 0.1/(double)NVID; + wave_height = wave_packet_height(t, packet[i], type_envelope); + for (j=0; j rmin2)) + { + r2 = (double)(irad2); + wave_height1 = wave_height*exp(-r2/variance); + psi[packet[i].ix + j][packet[i].iy + k] = wave_height1; + psi[packet[i].ix - j][packet[i].iy + k] = wave_height1; + psi[packet[i].ix + j][packet[i].iy - k] = wave_height1; + psi[packet[i].ix - j][packet[i].iy - k] = wave_height1; + } + else if (irad2 <= rmin2) + { + psi[packet[i].ix + j][packet[i].iy + k] = 0.0; + psi[packet[i].ix - j][packet[i].iy + k] = 0.0; + psi[packet[i].ix + j][packet[i].iy - k] = 0.0; + psi[packet[i].ix - j][packet[i].iy - k] = 0.0; + } + } + } + } +} + +void add_circular_wave_loc(double factor, double x, double y, double *phi[NX], double *psi[NX], short int * xy_in[NX], int xmin, int xmax, int ymin, int ymax) +/* add drop at (x,y) to the field with given prefactor */ +{ + int i, j; + double xy[2], dist2; + +// for (i=0; i (double)packet[i].var_envelope) amp *= exp(-0.1*(t - (double)packet[i].var_envelope)); + if (t < (double)packet[i].var_envelope + 100.0) + add_circular_wave_loc(amp, packet[i].xc, packet[i].yc, phi, psi, xy_in, xmin, xmax, ymin, ymax); + } +} + + +void add_wave_packets(double *phi[NX], double *psi[NX], short int * xy_in[NX], t_wave_packet *packet, int time, int radius, int local, int add_period, int type_envelope) +/* add some wave packet sources */ +{ + if (local) add_wave_packets_locally(phi, psi, packet, time, radius, add_period, type_envelope); + else add_wave_packets_globally(phi, psi, xy_in, packet, time); +} diff --git a/sub_wave_3d.c b/sub_wave_3d.c index 47edb7a..27c8256 100644 --- a/sub_wave_3d.c +++ b/sub_wave_3d.c @@ -1106,6 +1106,13 @@ void compute_light_angle(short int xy_in[NX*NY], t_wave wave[NX*NY], int movie) { gradx = (*wave[(i+1)*NY+j].p_zfield[movie] - *wave[(i-1)*NY+j].p_zfield[movie])/dx; grady = (*wave[i*NY+j+1].p_zfield[movie] - *wave[i*NY+j-1].p_zfield[movie])/dy; + + if (ADD_POTENTIAL) + { + gradx -= (*wave[(i+1)*NY+j].potential - *wave[(i-1)*NY+j].potential)*POT_FACT/dx; + grady -= (*wave[i*NY+j+1].potential - *wave[i*NY+j-1].potential)*POT_FACT/dy; + } + norm = sqrt(1.0 + gradx*gradx + grady*grady); pscal = -gradx*light[0] - grady*light[1] + 1.0; @@ -1185,14 +1192,14 @@ double compute_interpolated_colors_wave(int i, int j, short int xy_in[NX*NY], t_ { int k; double cw, ce, cn, cs, c_sw, c_se, c_nw, c_ne, c_mid, ca, z_mid; - double cw2, ce2, cn2, cs2; + double cw2, ce2, cn2, cs2, factor; double *z_sw, *z_se, *z_nw, *z_ne; z_sw = wave[i*NY+j].p_zfield[movie]; z_se = wave[(i+1)*NY+j].p_zfield[movie]; z_nw = wave[i*NY+j+1].p_zfield[movie]; z_ne = wave[(i+1)*NY+j+1].p_zfield[movie]; - + z_mid = 0.25*(*z_sw + *z_se + *z_nw + *z_ne); c_sw = *wave[i*NY+j].p_cfield[movie]; @@ -1251,6 +1258,16 @@ double compute_interpolated_colors_wave(int i, int j, short int xy_in[NX*NY], t_ rgb_s[k] *= fade_value; } +// if (ADD_POTENTIAL) +// { +// factor = 0.25*POT_FACT; +// z_mid += *wave[i*NY+j].potential*factor; +// z_mid += *wave[(i+1)*NY+j].potential*factor; +// z_mid += *wave[i*NY+j+1].potential*factor; +// z_mid += *wave[(i+1)*NY+j+1].potential*factor; +// } + + return(z_mid); } @@ -1275,9 +1292,10 @@ void compute_wave_fields(double phi[NX*NY], double psi[NX*NY], short int xy_in[N void init_speed_dissipation(short int xy_in[NX*NY], double tc[NX*NY], double tcc[NX*NY], double tgamma[NX*NY]) /* initialise fields for wave speed and dissipation */ { - int i, j, k, inlens; + int i, j, k, n, inlens; 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; + 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]; if (VARIABLE_IOR) { @@ -1414,6 +1432,44 @@ void init_speed_dissipation(short int xy_in[NX*NY], double tc[NX*NY], double tcc } break; } + case (IOR_RANDOM_WELLS): + { + dx = (XMAX - XMIN)/(double)NGRIDX; + dy = (YMAX - YMIN)/(double)NGRIDY; + sigma = 0.2*dx*dx; + for (i=0; i 1.0) height[n] = 1.0; + if (height[n] < 0.0) height[n] = 0.0; + } + +// #pragma omp parallel for private(i,j) + for (i=0; i #include -#define MOVIE 1 /* set to 1 to generate movie */ +#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 when writing tiff images */ #define NO_EXTRA_BUFFER_SWAP 1 /* some OS require one less buffer swap when recording images */ @@ -54,12 +54,8 @@ // #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 1800 /* number of grid points on x axis */ -// #define NY 1800 /* number of grid points on y axis */ -// // #define NX 2700 /* number of grid points on x axis */ -// // #define NY 1350 /* number of grid points on y axis */ -// // #define YMIN -1.041666667 -// // #define YMAX 1.041666667 /* y interval for 9/16 aspect ratio */ +// #define NX 2840 /* 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 */ @@ -72,30 +68,25 @@ #define WINHEIGHT 720 /* window height */ // #define NX 1280 /* number of grid points on x axis */ -#define NX 720 /* number of grid points on x axis */ -#define NY 720 /* number of grid points on y axis */ -// #define NX 2560 /* number of grid points on x axis */ +// #define NX 720 /* number of grid points on x axis */ +// #define NY 720 /* number of grid points on y axis */ +#define NX 2560 /* number of grid points on x axis */ // #define NX 1440 /* number of grid points on x axis */ -// #define NY 1440 /* number of grid points on y axis */ +#define NY 1440 /* number of grid points on y axis */ // #define NX 360 /* number of grid points on x axis */ // #define NY 360 /* number of grid points on y axis */ -// #define XMIN -2.0 -// #define XMAX 2.0 /* x interval */ -// #define YMIN -1.125 -// #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ -#define XMIN -1.1 -#define XMAX 1.1 /* x interval */ -#define YMIN -1.1 -#define YMAX 1.1 /* y interval */ +#define XMIN -2.0 +#define XMAX 2.0 /* x interval */ +#define YMIN -1.125 +#define YMAX 1.125 /* y interval for 9/16 aspect ratio */ #define JULIA_SCALE 0.8 /* scaling for Julia sets */ /* Choice of the billiard table */ -#define B_DOMAIN 32 /* 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 999 /* 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 */ @@ -104,15 +95,16 @@ #define CIRCLE_PATTERN_B 0 /* second pattern of circles or polygons */ #define VARIABLE_IOR 1 /* set to 1 for a variable index of refraction */ -#define IOR 3 /* choice of index of refraction, see list in global_pdes.c */ +#define IOR 5 /* 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 RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */ -#define LAMBDA 0.75 /* parameter controlling the dimensions of domain */ -#define MU 0.25 /* parameter controlling the dimensions of domain */ +#define LAMBDA 0.5 /* parameter controlling the dimensions of domain */ +#define MU 0.5 /* 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 +112,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 12 /* number of grid point for grid of disks */ -#define NGRIDY 16 /* number of grid point for grid of disks */ +#define NGRIDX 14 /* number of grid point for grid of disks */ +#define NGRIDY 8 /* number of grid point for grid of disks */ #define X_SHOOTER -0.2 #define Y_SHOOTER -0.6 @@ -150,11 +142,10 @@ #define AMPLITUDE 0.8 /* amplitude of periodic excitation */ #define ACHIRP 0.2 /* acceleration coefficient in chirp */ #define DAMPING 0.0 /* damping of periodic excitation */ -#define COURANT 0.1 /* Courant number */ -// #define COURANTB 0.0 /* Courant number in medium B */ -#define COURANTB 0.04658753 /* Courant number in medium B */ +#define COURANT 0.05 /* Courant number */ +#define COURANTB 0.0 /* Courant number in medium B */ #define GAMMA 0.0 /* damping factor in wave equation */ -#define GAMMAB 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 */ #define GAMMA_TOPBOT 1.0e-7 /* damping factor on boundary */ #define KAPPA 0.0 /* "elasticity" term enforcing oscillations */ @@ -166,8 +157,13 @@ /* 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 100 /* period of oscillating source */ -#define ALTERNATE_OSCILLATING_SOURCE 0 /* set to 1 to alternate sign of oscillating source */ +#define OSCILLATING_SOURCE_PERIOD 75 /* 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 */ +#define WAVE_PACKET_SOURCE_TYPE 1 /* type of wave packet sources */ +#define N_WAVE_PACKETS 15 /* number of wave packets */ +#define WAVE_PACKET_RADIUS 20 /* radius of wave packets */ /* Boundary conditions, see list in global_pdes.c */ @@ -178,9 +174,9 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 2800 /* number of frames of movie */ -// #define NSTEPS 200 /* number of frames of movie */ -#define NVID 6 /* number of iterations between images displayed on screen */ +#define NSTEPS 2200 /* number of frames of movie */ +// #define NSTEPS 500 /* number of frames of movie */ +#define NVID 4 /* 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 */ @@ -197,27 +193,23 @@ /* Parameters of initial condition */ -// #define INITIAL_AMP 0.25 /* amplitude of initial condition */ -// #define INITIAL_VARIANCE 0.001 /* variance of initial condition */ -// #define INITIAL_WAVELENGTH 0.1 /* wavelength of initial condition */ -#define INITIAL_AMP 0.3 /* amplitude of initial condition */ -#define INITIAL_VARIANCE 0.0004 /* variance of initial condition */ -#define INITIAL_WAVELENGTH 0.02 /* wavelength of initial condition */ -// #define INITIAL_AMP 0.35 /* amplitude of initial condition */ -// #define INITIAL_VARIANCE 0.0002 /* variance of initial condition */ -// #define INITIAL_WAVELENGTH 0.01 /* wavelength of initial condition */ +#define INITIAL_AMP 0.25 /* amplitude of initial condition */ +#define INITIAL_VARIANCE 0.00015 /* variance of initial condition */ +#define INITIAL_WAVELENGTH 0.0075 /* wavelength of initial condition */ /* Plot type, see list in global_pdes.c */ -// #define ZPLOT 103 /* wave height */ -// #define CPLOT 103 /* color scheme */ -#define ZPLOT 104 /* wave height */ -#define CPLOT 104 /* color scheme */ +#define ZPLOT 103 /* wave height */ +#define CPLOT 103 /* color scheme */ +// #define ZPLOT 104 /* wave height */ +// #define CPLOT 104 /* color scheme */ #define ZPLOT_B 108 #define CPLOT_B 108 /* 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 */ #define AMPLITUDE_HIGH_RES 1 /* set to 1 to increase resolution of plot */ @@ -230,7 +222,7 @@ #define FADE_IN_OBSTACLE 1 /* set to 1 to fade color inside obstacles */ #define DRAW_OUTSIDE_GRAY 0 /* experimental, draw outside of billiard in gray */ -#define PLOT_SCALE_ENERGY 0.01 /* vertical scaling in energy plot */ +#define PLOT_SCALE_ENERGY 0.1 /* vertical scaling in energy plot */ #define PLOT_SCALE_LOG_ENERGY 0.2 /* vertical scaling in log energy plot */ /* 3D representation */ @@ -241,26 +233,26 @@ #define REP_PROJ_3D 1 /* projection on plane orthogonal to observer line of sight */ #define ROTATE_VIEW 1 /* set to 1 to rotate position of observer */ -#define ROTATE_ANGLE 540.0 /* total angle of rotation during simulation */ +#define ROTATE_ANGLE 360.0 /* total angle of rotation during simulation */ // #define ROTATE_ANGLE 45.0 /* total angle of rotation during simulation */ /* Color schemes */ #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 COLOR_PALETTE_B 14 /* Color palette, see list in global_pdes.c */ #define BLACK 1 /* background */ #define COLOR_SCHEME 3 /* choice of color scheme, see list in global_pdes.c */ #define SCALE 0 /* set to 1 to adjust color scheme to variance of field */ -#define SLOPE 0.5 /* sensitivity of color on wave amplitude */ -#define VSCALE_AMPLITUDE 1.0 /* additional scaling factor for color scheme P_3D_AMPLITUDE */ -#define VSCALE_ENERGY 30.0 /* additional scaling factor for color scheme P_3D_ENERGY */ +#define SLOPE 1.0 /* sensitivity of color on wave amplitude */ +#define VSCALE_AMPLITUDE 3.0 /* additional scaling factor for color scheme P_3D_AMPLITUDE */ +#define VSCALE_ENERGY 25.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 100.0 /* scaling factor for energy representation */ +#define E_SCALE 240.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,18 +273,19 @@ #define MAZE_MAX_NGBH 5 /* max number of neighbours of maze cell */ #define RAND_SHIFT 5 /* seed of random number generator */ #define MAZE_XSHIFT 0.0 /* horizontal shift of maze */ +#define MAZE_WIDTH 0.02 /* half width of maze walls */ #define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */ -#define COLORBAR_RANGE 6.0 /* scale of color scheme bar */ -#define COLORBAR_RANGE_B 6.0 /* scale of color scheme bar for 2nd part */ +#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 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 */ -/* for compatibility with sub_wave and sub_maze */ -#define ADD_POTENTIAL 0 +#define ADD_POTENTIAL 0 /* set to 1 to add potential to z coordinate */ // #define POT_MAZE 7 -#define POTENTIAL 0 +#define POTENTIAL 10 +#define POT_FACT 30.0 /* end of constants only used by sub_wave and sub_maze */ @@ -306,11 +299,11 @@ 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, 12.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 0.25 /* overall scaling factor of z axis for REP_PROJ_3D representation */ -#define XY_SCALING_FACTOR 2.0 /* overall scaling factor for on-screen (x,y) coordinates after projection */ +#define Z_SCALING_FACTOR 0.15 /* 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 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.1 /* overall y shift for REP_PROJ_3D representation */ @@ -952,11 +945,11 @@ void viewpoint_schedule(int i) void animation() { - double time, scale, ratio, startleft[2], startright[2], sign = 1.0, r2, xy[2], fade_value, yshift, speed = 0.0, a, b, c, angle, lambda1; + double time, scale, ratio, startleft[2], startright[2], sign = 1.0, r2, xy[2], fade_value, yshift, speed = 0.0, a, b, c, angle, lambda1, y, x1, sign1; double *phi, *psi, *tmp, *color_scale, *tc, *tcc, *tgamma; // double *total_energy; short int *xy_in; - int i, j, s, sample_left[2], sample_right[2], period = 0, fade; + int i, j, s, sample_left[2], sample_right[2], period = 0, fade, source_counter = 0, k, p, q; static int counter = 0; long int wave_value; t_wave *wave; @@ -1065,11 +1058,11 @@ void animation() // init_circular_wave_mod(LAMBDA*cos(APOLY*PID), LAMBDA*sin(APOLY*PID), phi, psi, xy_in); lambda1 = LAMBDA; angle = DPI/(double)NPOLY; - init_circular_wave_mod(lambda1*cos(0.5*angle), lambda1*sin(0.5*angle), phi, psi, xy_in); - for (j=1; j 0.083333333*XMIN)&&(x1 < 0.083333333*XMAX)) +// { +// add_circular_wave_mod(sign1, x1, y, phi, psi, xy_in); +// printf("Adding wave at (%.2lg, %.2lg)\n", x1, y); +// } +// sign1 = -sign1; +// } +// source_counter++; +// if (p > 0) q = p; +// else q = -p; +// if (source_counter >= q) +// { +// source_counter = 0; +// sign = -sign; +// } + } if (PRINT_SPEED) print_speed_3d(speed, 0, 1.0); diff --git a/wave_billiard.c b/wave_billiard.c index 428c4bd..9b18cd6 100644 --- a/wave_billiard.c +++ b/wave_billiard.c @@ -43,43 +43,45 @@ #include #include -#define MOVIE 1 /* set to 1 to generate movie */ +#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 when writing tiff images */ #define NO_EXTRA_BUFFER_SWAP 1 /* some OS require one less buffer swap when recording images */ +#define VARIABLE_IOR 1 /* set to 1 for a variable index of refraction */ +#define IOR 5 /* 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 */ + + /* General geometrical parameters */ -#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 x axis */ -#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 YMIN -1.197916667 -#define YMAX 1.197916667 /* y interval for 9/16 aspect ratio */ -// #define YMIN -1.041666667 -// #define YMAX 1.041666667 /* y interval for 9/16 aspect ratio */ -// -#define HIGHRES 1 /* set to 1 if resolution of grid is double that of displayed image */ - -// #define WINWIDTH 1280 /* window width */ -// #define WINHEIGHT 720 /* window height */ -// -// // #define NX 640 /* number of grid points on x axis */ -// // #define NY 360 /* number of grid points on y axis */ -// // #define NX 1280 /* number of grid points on x axis */ -// // #define NY 720 /* number of grid points on y axis */ -// #define NX 2560 /* number of grid points on x axis */ -// #define NY 1440 /* number of grid points on y axis */ +// #define WINWIDTH 1920 /* window width */ +// #define WINHEIGHT 1150 /* window height */ +// #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 YMIN -1.125 -// #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ +// #define XMAX 2.0 /* x interval */ +// #define YMIN -1.197916667 +// #define YMAX 1.197916667 /* y interval for 9/16 aspect ratio */ + +#define HIGHRES 1 /* set to 1 if resolution of grid is double that of displayed image */ + +#define WINWIDTH 1280 /* window width */ +#define WINHEIGHT 720 /* window height */ + +// #define NX 640 /* number of grid points on x axis */ +// #define NY 360 /* number of grid points on y axis */ +// #define NX 1280 /* number of grid points on x axis */ +// #define NY 720 /* number of grid points on y axis */ +#define NX 2560 /* number of grid points on x axis */ +#define NY 1440 /* number of grid points on y axis */ + +#define XMIN -2.0 +#define XMAX 2.0 /* x interval */ +#define YMIN -1.125 +#define YMAX 1.125 /* y interval for 9/16 aspect ratio */ #define JULIA_SCALE 1.0 /* scaling for Julia sets */ @@ -106,8 +108,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 12 /* number of grid point for grid of disks */ -#define NGRIDY 15 /* number of grid point for grid of disks */ +#define NGRIDX 12 /* number of grid point for grid of disks */ +#define NGRIDY 12 /* number of grid point for grid of disks */ #define X_SHOOTER -0.2 #define Y_SHOOTER -0.6 @@ -134,8 +136,8 @@ #define AMPLITUDE 0.8 /* amplitude of periodic excitation */ #define ACHIRP 0.25 /* acceleration coefficient in chirp */ #define DAMPING 0.0 /* damping of periodic excitation */ -#define COURANT 0.08 /* Courant number */ -#define COURANTB 0.04658753 /* Courant number in medium B */ +#define COURANT 0.04 /* Courant number */ +#define COURANTB 0.0 /* 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 */ @@ -149,8 +151,13 @@ /* 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 15 /* period of oscillating source */ -#define ALTERNATE_OSCILLATING_SOURCE 0 /* set to 1 to alternate sign of oscillating source */ +#define OSCILLATING_SOURCE_PERIOD 50 /* period of oscillating source */ +#define ALTERNATE_OSCILLATING_SOURCE 1 /* set to 1 to alternate sign of oscillating source */ + +#define ADD_WAVE_PACKET_SOURCES 0 /* set to 1 to add several sources emitting wave packets */ +#define WAVE_PACKET_SOURCE_TYPE 1 /* type of wave packet sources */ +#define N_WAVE_PACKETS 15 /* number of wave packets */ +#define WAVE_PACKET_RADIUS 20 /* radius of wave packets */ /* Boundary conditions, see list in global_pdes.c */ @@ -158,13 +165,11 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 2500 /* number of frames of movie */ +#define NSTEPS 1800 /* number of frames of movie */ // #define NSTEPS 500 /* number of frames of movie */ -#define NVID 12 /* number of iterations between images displayed on screen */ -// #define NVID 9 /* 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 */ -// #define INITIAL_TIME 400 /* time after which to start saving frames */ #define BOUNDARY_WIDTH 2 /* width of billiard boundary */ #define PRINT_SPEED 0 /* print speed of moving source */ @@ -178,24 +183,22 @@ /* Parameters of initial condition */ -#define INITIAL_AMP 0.1 /* amplitude of initial condition */ +#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 INITIAL_VARIANCE 0.00015 /* variance of initial condition */ -// #define INITIAL_WAVELENGTH 0.0075 /* wavelength of initial condition */ /* Plot type, see list in global_pdes.c */ -#define PLOT 7 +#define PLOT 0 // #define PLOT 7 -#define PLOT_B 6 /* 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 17 /* Color palette, see list in global_pdes.c */ -#define COLOR_PALETTE_B 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 */ @@ -206,10 +209,10 @@ #define PHASE_FACTOR 1.0 /* factor in computation of phase in color scheme P_3D_PHASE */ #define PHASE_SHIFT 0.0 /* shift of phase in color scheme P_3D_PHASE */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ -#define E_SCALE 140.0 /* scaling factor for energy representation */ +#define E_SCALE 60.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 FLUX_SCALE 2.0e3 /* scaling factor for enegy flux represtnation */ +#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) */ #define COLORHUE 260 /* initial hue of water color for scheme C_LUM */ @@ -219,9 +222,9 @@ #define HUEMEAN 180.0 /* mean value of hue for color scheme C_HUE */ #define HUEAMP -180.0 /* amplitude of variation of hue for color scheme C_HUE */ -#define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */ -#define COLORBAR_RANGE 1.5 /* scale of color scheme bar */ -#define COLORBAR_RANGE_B 2.5 /* scale of color scheme bar for 2nd part */ +#define DRAW_COLOR_SCHEME 0 /* set to 1 to plot the color scheme */ +#define COLORBAR_RANGE 2.0 /* scale of color scheme bar */ +#define COLORBAR_RANGE_B 0.1 /* 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 */ @@ -231,6 +234,7 @@ #define MAZE_MAX_NGBH 5 /* max number of neighbours of maze cell */ #define RAND_SHIFT 0 /* seed of random number generator */ #define MAZE_XSHIFT 0.0 /* horizontal shift of maze */ +#define MAZE_WIDTH 0.02 /* half width of maze walls */ /* for compatibility with sub_wave and sub_maze */ #define ADD_POTENTIAL 0 @@ -243,6 +247,7 @@ #define VMAX 10.0 /* max value of wave amplitude */ #define MEAN_FLUX (PLOT == P_TOTAL_ENERGY_FLUX)||(PLOT_B == P_TOTAL_ENERGY_FLUX) +#define REFRESH_IOR ((IOR == IOR_PERIODIC_WELLS_ROTATING)||(IOR == IOR_PERIODIC_WELLS_ROTATING_LARGE)) #include "global_pdes.c" /* constants and global variables */ #include "sub_maze.c" /* support for generating mazes */ @@ -257,11 +262,10 @@ double courant2, courantb2; /* Courant parameters squared */ /* animation part */ /*********************/ - // void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX], double *psi_out[NX], // short int *xy_in[NX]) void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX], - short int *xy_in[NX]) + short int *xy_in[NX], double *tcc[NX]) /* time step of field evolution */ /* phi is value of field at time t, psi at time t-1 */ /* this version of the function has been rewritten in order to minimize the number of if-branches */ @@ -269,7 +273,7 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX int i, j, iplus, iminus, jplus, jminus; double delta, x, y, c, cc, gamma, tb_shift; static long time = 0; - static double tc[NX][NY], tcc[NX][NY], tgamma[NX][NY]; + static double tc[NX][NY], tgamma[NX][NY]; static short int first = 1; time++; @@ -285,7 +289,7 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX if (xy_in[i][j] != 0) { tc[i][j] = COURANT; - tcc[i][j] = courant2; + if (!VARIABLE_IOR) tcc[i][j] = courant2; if (xy_in[i][j] == 1) tgamma[i][j] = GAMMA; else tgamma[i][j] = GAMMAB; } @@ -529,7 +533,7 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX } -void evolve_wave(double *phi[NX], double *psi[NX], double *tmp[NX], short int *xy_in[NX]) +void evolve_wave(double *phi[NX], double *psi[NX], double *tmp[NX], short int *xy_in[NX], double *tcc_table[NX]) /* time step of field evolution */ /* phi is value of field at time t, psi at time t-1 */ { @@ -539,11 +543,11 @@ void evolve_wave(double *phi[NX], double *psi[NX], double *tmp[NX], short int *x // At the beginning w[t] is saved in phi, w[t-1] in psi and tmp is space // for the next wave state w[t+1]. Take w[t] and w[t-1] to calculate the // next wave state. Write this new state in temp - evolve_wave_half(phi, psi, tmp, xy_in); + evolve_wave_half(phi, psi, tmp, xy_in, tcc_table); // now w[t] is saved in tmp, w[t-1] in phi and the result is written to psi - evolve_wave_half(tmp, phi, psi, xy_in); + evolve_wave_half(tmp, phi, psi, xy_in, tcc_table); // now w[t] is saved in psi, w[t-1] in tmp and the result is written to phi - evolve_wave_half(psi, tmp, phi, xy_in); + evolve_wave_half(psi, tmp, phi, xy_in, tcc_table); // now w[t] is saved in phi, w[t-1] in psi and tmp is free again to take // the new wave state w[t+1] in the next call to this function, thus // matching the given parameter names again @@ -577,12 +581,13 @@ void draw_color_bar_palette(int plot, double range, int palette, int fade, doubl void animation() { - double time, scale, ratio, startleft[2], startright[2], sign = 1.0, r2, xy[2], fade_value, yshift, speed = 0.0, a, b, c, x, y, angle, x1, sign1; - double *phi[NX], *psi[NX], *tmp[NX], *total_energy[NX], *color_scale[NX], *total_flux; + double time, scale, ratio, startleft[2], startright[2], sign = 1.0, r2, xy[2], fade_value, yshift, speed = 0.0, a, b, c, x, y, angle, x1, sign1, ior_angle = 0.0; + double *phi[NX], *psi[NX], *tmp[NX], *total_energy[NX], *color_scale[NX], *total_flux, *tcc_table[NX]; short int *xy_in[NX]; - int i, j, k, s, sample_left[2], sample_right[2], period = 0, fade, source_counter = 0; + int i, j, k, s, sample_left[2], sample_right[2], period = 0, fade, source_counter = 0, p, q; static int counter = 0; long int wave_value; + t_wave_packet *packet; if (SAVE_TIME_SERIES) { @@ -599,10 +604,17 @@ void animation() total_energy[i] = (double *)malloc(NY*sizeof(double)); xy_in[i] = (short int *)malloc(NY*sizeof(short int)); color_scale[i] = (double *)malloc(NY*sizeof(double)); + tcc_table[i] = (double *)malloc(NX*sizeof(double)); } if (MEAN_FLUX) total_flux = (double *)malloc(4*NX*NY*sizeof(double)); + if (ADD_WAVE_PACKET_SOURCES) + { + packet = (t_wave_packet *)malloc(N_WAVE_PACKETS*sizeof(t_wave_packet)); + init_wave_packets(packet, WAVE_PACKET_RADIUS); + } + /* initialise positions and radii of circles */ if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN == D_CIRCLES_IN_RECT)) ncircles = init_circle_config(circles); else if (B_DOMAIN == D_POLYGONS) ncircles = init_polygon_config(polygons); @@ -649,6 +661,8 @@ void animation() for (i=0; i<4*NX*NY; i++) total_flux[i] = 0.0; + if (VARIABLE_IOR) init_ior_2d(xy_in, tcc_table, ior_angle); + ratio = (XMAX - XMIN)/8.4; /* for Tokarsky billiard */ // isospectral_initial_point(0.2, 0.0, startleft, startright); /* for isospectral billiards */ @@ -662,7 +676,7 @@ void animation() init_wave_flat(phi, psi, xy_in); -// init_circular_wave(sqrt(LAMBDA*LAMBDA - 1.0), 0.0, phi, psi, xy_in); +// init_circular_wave(0.0, 0.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); @@ -753,7 +767,7 @@ void animation() for (j=0; j XMIN)&&(x1 < XMAX)) - { - add_circular_wave(sign1, x1, y, phi, psi, xy_in); - printf("Adding wave at (%.2lg, %.2lg)\n", x1, y); - } - sign1 = -sign1; - } - source_counter++; - if (source_counter == 4) - { - source_counter = 0; - sign = -sign; - } +// p = phased_array_schedule(i); +// p = 3; +// y = -1.0; +// sign1 = sign; +// printf("p = %i\n", p); +// for (k=-12; k<13; k++) +// { +// x1 = 0.05*((double)source_counter/(double)p + (double)k); +// if ((x1 > 0.1*XMIN)&&(x1 < 0.1*XMAX)) +// { +// add_circular_wave(sign1, x1, y, phi, psi, xy_in); +// printf("Adding wave at (%.2lg, %.2lg)\n", x1, y); +// } +// sign1 = -sign1; +// } +// source_counter++; +// if (p > 0) q = p; +// else q = -p; +// if (source_counter >= q) +// { +// source_counter = 0; +// sign = -sign; +// } // for (j=0; j