From f570f6885ee7d1cb0ed3d105daaa3a3d8a69ff1e Mon Sep 17 00:00:00 2001 From: nilsberglund-orleans <83530463+nilsberglund-orleans@users.noreply.github.com> Date: Sun, 13 Mar 2022 15:29:50 +0100 Subject: [PATCH] Add files via upload --- colormaps.c | 207 ++++ colors_waves.c | 89 ++ global_ljones.c | 39 +- global_pdes.c | 2 + lennardjones.c | 1925 ++++++++--------------------------- sub_hashgrid.c | 594 +++++++++++ sub_lj.c | 2461 +++++++++++++++++++++++++++++++++++++++++++++ sub_wave.c | 247 ++++- sub_wave_comp.c | 74 +- wave_common.c | 362 ++++++- wave_comparison.c | 157 ++- 11 files changed, 4564 insertions(+), 1593 deletions(-) create mode 100644 sub_hashgrid.c diff --git a/colormaps.c b/colormaps.c index a420008..e15f6ff 100644 --- a/colormaps.c +++ b/colormaps.c @@ -2339,3 +2339,210 @@ void hsl_to_rgb_twilight(double h, double s, double l, double rgb[3]) else if (color_hue < 0) color_hue = 0; for (i=0; i<3; i++) rgb[i] = 2.0*l*(double)twilight_srgb_floats[color_hue][i]; } + +void hsl_to_rgb_palette(double h, double s, double l, double rgb[3], int palette) +/* color conversion from HSL to RGB */ +{ + int color_hue, i; + double r, g, b, ratio = 0.711111111, interpolate; /* ratio equals 256/360 */ + + switch (palette) { + case (COL_JET): + { + hsl_to_rgb_jet(h, s, l, rgb); + break; + } + case (COL_HSLUV): + { + hsluv2rgb(h, 100.0*s, l, &r, &g, &b); + rgb[0] = 2.0*l*r*100.0; + rgb[1] = 2.0*l*g*100.0; + rgb[2] = 2.0*l*b*100.0; + break; + } + case (COL_GRAY): + { + color_hue = h/360.0; + rgb[0] = color_hue; + rgb[1] = color_hue; + rgb[2] = color_hue; + break; + } + case (COL_TURBO): + { + color_hue = 255 - (int)(ratio*h); + for (i=0; i<3; i++) rgb[i] = 2.0*l*(double)turbo_srgb_floats[color_hue][i]; + break; + } + case (COL_VIRIDIS): + { + color_hue = 255 - (int)(ratio*h); + for (i=0; i<3; i++) rgb[i] = 2.0*l*(double)viridis_srgb_floats[color_hue][i]; + break; + } + case (COL_MAGMA): + { + color_hue = 255 - (int)(ratio*h); + for (i=0; i<3; i++) rgb[i] = 2.0*l*(double)magma_srgb_floats[color_hue][i]; + break; + } + case (COL_INFERNO): + { + color_hue = 255 - (int)(ratio*h); + for (i=0; i<3; i++) rgb[i] = 2.0*l*(double)inferno_srgb_floats[color_hue][i]; + break; + } + case (COL_PLASMA): + { + color_hue = 255 - (int)(ratio*h); + for (i=0; i<3; i++) rgb[i] = 2.0*l*(double)plasma_srgb_floats[color_hue][i]; + break; + } + case (COL_CIVIDIS): + { + color_hue = 255 - (int)(ratio*h); + for (i=0; i<3; i++) rgb[i] = 2.0*l*(double)cividis_srgb_floats[color_hue][i]; + break; + } + case (COL_PARULA): + { + color_hue = 255 - (int)(ratio*h); + for (i=0; i<3; i++) rgb[i] = 2.0*l*(double)parula_srgb_floats[color_hue][i]; + break; + } + case (COL_TWILIGHT): + { + color_hue = 510 - (int)(2.0*ratio*h); + if (color_hue > 510) color_hue = 510; + else if (color_hue < 0) color_hue = 0; + for (i=0; i<3; i++) rgb[i] = 2.0*l*(double)twilight_srgb_floats[color_hue][i]; + break; + } + case (COL_TWILIGHT_SHIFTED): + { +// color_hue = 255 - (int)(2.0*ratio*h); + color_hue = (int)(2.0*ratio*h) - 255; + if (color_hue > 510) color_hue -= 510; + else if (color_hue < 0) color_hue += 510; + for (i=0; i<3; i++) rgb[i] = 2.0*l*(double)twilight_srgb_floats[color_hue][i]; + break; + } + case (COL_TURBO_CYCLIC): + { + if (h < 330.0) + { + color_hue = 255 - (int)(ratio*h*12.0/11.0); + for (i=0; i<3; i++) rgb[i] = 2.0*l*(double)turbo_srgb_floats[color_hue][i]; + } + else /* use a convex combination to interpolate between extremal Turbo values */ + { + if (h > 360.0) h = 360.0; + interpolate = (h - 330.0)/30.0; + for (i=0; i<3; i++) + rgb[i] = (1.0 - interpolate)*(double)turbo_srgb_floats[255][i] + interpolate*(double)turbo_srgb_floats[0][i]; + } + break; + } + } +} + +void amp_to_rgb_palette(double h, double rgb[3], int palette) /* color conversion from amplitude in [0,1] to RGB */ +{ + int color_hue, i; + double r, g, b, interpolate; + + color_hue = (int)(256.0*h); + if (color_hue > 255) color_hue = 255; + + switch (palette) { + case (COL_JET): + { + hsl_to_rgb_jet(360.0*(1.0 - h), 0.9, 0.5, rgb); + break; + } + case (COL_HSLUV): + { + hsluv2rgb(360.0*(1.0 - h), 90.0, 0.5, &r, &g, &b); + rgb[0] = r*100.0; + rgb[1] = g*100.0; + rgb[2] = b*100.0; + break; + } + case (COL_GRAY): + { + rgb[0] = h; + rgb[1] = h; + rgb[2] = h; + break; + } + case (COL_TURBO): + { + for (i=0; i<3; i++) rgb[i] = (double)turbo_srgb_floats[color_hue][i]; + break; + } + case (COL_VIRIDIS): + { + for (i=0; i<3; i++) rgb[i] = (double)viridis_srgb_floats[color_hue][i]; + break; + } + case (COL_MAGMA): + { + for (i=0; i<3; i++) rgb[i] = (double)magma_srgb_floats[color_hue][i]; + break; + } + case (COL_INFERNO): + { + for (i=0; i<3; i++) rgb[i] = (double)inferno_srgb_floats[color_hue][i]; + break; + } + case (COL_PLASMA): + { + for (i=0; i<3; i++) rgb[i] = (double)plasma_srgb_floats[color_hue][i]; + break; + } + case (COL_CIVIDIS): + { + for (i=0; i<3; i++) rgb[i] = (double)cividis_srgb_floats[color_hue][i]; + break; + } + case (COL_PARULA): + { + for (i=0; i<3; i++) rgb[i] = (double)parula_srgb_floats[color_hue][i]; + break; + } + case (COL_TWILIGHT): + { + color_hue = (int)(511.0*h); + if (color_hue > 510) color_hue = 510; + else if (color_hue < 0) color_hue = 0; + for (i=0; i<3; i++) rgb[i] = (double)twilight_srgb_floats[color_hue][i]; + break; + } + case (COL_TWILIGHT_SHIFTED): + { +// color_hue = (int)(511.0*h) + 255; + color_hue = 255 - (int)(511.0*h); + if (color_hue > 510) color_hue -= 510; + else if (color_hue < 0) color_hue += 510; + for (i=0; i<3; i++) rgb[i] = (double)twilight_srgb_floats[color_hue][i]; + break; + } + case (COL_TURBO_CYCLIC): + { + if (h < 0.9) + { + color_hue = (int)(256.0*h/0.9); + for (i=0; i<3; i++) rgb[i] = (double)turbo_srgb_floats[color_hue][i]; + } + else /* use a convex combination to interpolate between extremal Turbo values */ + { + if (h > 1.0) h = 1.0; + interpolate = (h - 0.9)/0.1; + for (i=0; i<3; i++) + rgb[i] = (1.0 - interpolate)*(double)turbo_srgb_floats[255][i] + interpolate*(double)turbo_srgb_floats[0][i]; + } + break; + } + } +} + diff --git a/colors_waves.c b/colors_waves.c index 6c88d67..dc0ad4c 100644 --- a/colors_waves.c +++ b/colors_waves.c @@ -141,3 +141,92 @@ void color_scheme_asym(int scheme, double value, double scale, int time, double } } +void color_scheme_palette(int scheme, int palette, double value, double scale, int time, double rgb[3]) /* color scheme */ +{ + double hue, y, r, amplitude; + int intpart; + + /* saturation = r, luminosity = y */ + switch (scheme) { + case (C_LUM): + { + hue = COLORHUE + (double)time*COLORDRIFT/(double)NSTEPS; + if (hue < 0.0) hue += 360.0; + if (hue >= 360.0) hue -= 360.0; + r = 0.9; + amplitude = color_amplitude(value, scale, time); + y = LUMMEAN + amplitude*LUMAMP; + intpart = (int)y; + y -= (double)intpart; + hsl_to_rgb_palette(hue, r, y, rgb, palette); + break; + } + case (C_HUE): + { + r = 0.9; + amplitude = color_amplitude(value, scale, time); + y = 0.5; + hue = HUEMEAN + amplitude*HUEAMP; + if (hue < 0.0) hue += 360.0; + if (hue >= 360.0) hue -= 360.0; + hsl_to_rgb_palette(hue, r, y, rgb, palette); + break; + } + case (C_ONEDIM): + { + amplitude = color_amplitude(value, scale, time); + amp_to_rgb_palette(0.5*(1.0 + amplitude), rgb, palette); + break; + } + case (C_ONEDIM_LINEAR): + { + amplitude = color_amplitude_linear(value, scale, time); + if (amplitude > 1.0) amplitude -= 1.0; + else if (amplitude < 0.0) amplitude += 1.0; + amp_to_rgb_palette(amplitude, rgb, palette); + break; + } + } +} + +void color_scheme_asym_palette(int scheme, int palette, double value, double scale, int time, double rgb[3]) /* color scheme */ +{ + double hue, y, r, amplitude; + int intpart; + + /* saturation = r, luminosity = y */ + switch (scheme) { + case (C_LUM): + { + hue = COLORHUE + (double)time*COLORDRIFT/(double)NSTEPS; + if (hue < 0.0) hue += 360.0; + if (hue >= 360.0) hue -= 360.0; + r = 0.9; + amplitude = color_amplitude(value, scale, time); + y = LUMMEAN + amplitude*LUMAMP; + intpart = (int)y; + y -= (double)intpart; + hsl_to_rgb_palette(hue, r, y, rgb, palette); + break; + } + case (C_HUE): + { + r = 0.9; + amplitude = color_amplitude_asym(value, scale, time); + y = 0.5; + hue = HUEMEAN + 0.8*amplitude*HUEAMP; + if (hue < 0.0) hue += 360.0; + if (hue >= 360.0) hue -= 360.0; + hsl_to_rgb_palette(hue, r, y, rgb, palette); + break; + } + case (C_ONEDIM): + { + amplitude = color_amplitude(value, scale, time); + amp_to_rgb_palette(amplitude, rgb, palette); + break; + } + } +} + + diff --git a/global_ljones.c b/global_ljones.c index 228b1e7..ef2fef1 100644 --- a/global_ljones.c +++ b/global_ljones.c @@ -13,6 +13,7 @@ #define NMAXCIRCLES 20000 /* total number of circles/polygons (must be at least NCX*NCY for square grid) */ #define MAXNEIGH 20 /* max number of neighbours kept in memory */ +#define NMAXOBSTACLES 100 /* max number of obstacles */ #define C_SQUARE 0 /* square grid of circles */ #define C_HEX 1 /* hexagonal/triangular grid of circles */ @@ -32,6 +33,10 @@ #define C_TWO 98 /* two concentric circles of different type */ #define C_NOTHING 99 /* no circle at all, for comparisons */ +/* pattern of additional obstacles */ +#define O_CORNERS 0 /* obstacles in the corners (for Boy b.c.) */ +#define O_GALTON_BOARD 1 /* Galton board pattern */ + /* particle interaction */ #define I_COULOMB 0 /* Coulomb force */ @@ -41,6 +46,7 @@ #define I_GOLDENRATIO 4 /* Lennard-Jones type with equilibria related by golden ratio */ #define I_LJ_DIPOLE 5 /* Lennard-Jones with a dipolar angle dependence */ #define I_LJ_QUADRUPOLE 6 /* Lennard-Jones with a quadropolar angle dependence */ +#define I_LJ_WATER 7 /* model for water molecule */ /* Boundary conditions */ @@ -50,7 +56,12 @@ #define BC_PERIODIC 3 /* periodic boundary conditions */ #define BC_PERIODIC_CIRCLE 4 /* periodic boundary conditions and harmonic b.c. outside moving circle */ #define BC_EHRENFEST 5 /* Ehrenfest urn-type configuration */ +#define BC_PERIODIC_FUNNEL 6 /* funnel with periodic boundary conditions */ +#define BC_RECTANGLE_LID 7 /* rectangular container with moving lid */ +#define BC_PERIODIC_TRIANGLE 8 /* periodic boundary conditions and harmonic b.c. outside moving triangle */ #define BC_KLEIN 11 /* Klein bottle (periodic with twisted vertical parts) */ +#define BC_SCREEN_BINS 12 /* harmonic boundary conditions outside screen area plus "bins" (for Galton board) */ +#define BC_BOY 13 /* Boy surface/projective plane (periodic with twisted horizontal and vertical parts) */ /* Plot types */ @@ -61,6 +72,7 @@ #define P_ANGLE 4 /* colors represent angle/spin of particle */ #define P_TYPE 5 /* colors represent type of particle */ #define P_DIRECTION 6 /* colors represent direction of velocity */ +#define P_ANGULAR_SPEED 7 /* colors represent angular speed */ /* Color schemes */ @@ -104,11 +116,12 @@ typedef struct double fy; /* y component of force on particle */ double torque; /* torque on particle */ short int thermostat; /* whether particle is coupled to thermostat */ - int hashx; /* hash grid positions of particles */ - int hashy; /* hash grid positions of particles */ - int neighb; /* number of neighbours */ - int neighbours[MAXNEIGH]; /* coordinates of neighbours */ - double nghangle[MAXNEIGH]; /* angles of neighbours */ + int hashcell; /* hash cell in which particle is located */ + int neighb; /* number of neighbours within given distance */ + int hash_nneighb; /* number of neighbours in hashgrid */ + int hashneighbour[9*HASHMAX]; /* particle numbers of neighbours in hashgrid */ + double deltax[9*HASHMAX]; /* relative position of neighbours */ + double deltay[9*HASHMAX]; /* relative position of neighbours */ short int type; /* type of particle, for mixture simulations */ short int interaction; /* type of interaction */ double eq_dist; /* equilibrium distance */ @@ -120,6 +133,8 @@ typedef struct { int number; /* total number of particles in cell */ int particles[HASHMAX]; /* numbers of particles in cell */ + int nneighb; /* number of neighbouring cells */ + int neighbour[9]; /* numbers of neighbouring cells */ } t_hashgrid; typedef struct @@ -135,12 +150,22 @@ typedef struct int hashx; /* hash grid positions of particles */ int hashy; /* hash grid positions of particles */ int neighb; /* number of neighbours */ - int health; /* 0 = sane, 1 = infected, 2 = recovered */ + int health; /* 0 = healthy, 1 = infected, 2 = recovered */ double infected_time; /* time since infected */ int protected; /* 0 = not protected, 1 = protected */ } t_person; +typedef struct +{ + double xc, yc, radius; /* center and radius of circle */ + short int active; /* circle is active */ +} t_obstacle; + +typedef struct +{ + double xc, yc; /* center of circle */ +} t_tracer; -int ncircles, counter = 0; +int ncircles, nobstacles, counter = 0; diff --git a/global_pdes.c b/global_pdes.c index 4ce3b2b..9e4d78b 100644 --- a/global_pdes.c +++ b/global_pdes.c @@ -53,6 +53,8 @@ #define D_FRESNEL 43 /* Fresnel lens */ #define D_NOISEPANEL 44 /* zigzag noise insulating panel */ #define D_DOUBLE_FRESNEL 45 /* two facing Fresnel lenses */ +#define D_QRD 46 /* quadratic resonance diffuser */ +#define D_CIRCLE_SEGMENT 47 /* lens-shaped circular segment */ #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) */ diff --git a/lennardjones.c b/lennardjones.c index b4f94cf..f81c773 100644 --- a/lennardjones.c +++ b/lennardjones.c @@ -38,6 +38,11 @@ #define MOVIE 0 /* set to 1 to generate movie */ #define DOUBLE_MOVIE 0 /* set to 1 to produce movies for wave height and energy simultaneously */ +#define TIME_LAPSE 0 /* set to 1 to add a time-lapse movie at the end */ + /* so far incompatible with double movie */ +#define TIME_LAPSE_FACTOR 3 /* factor of time-lapse movie */ + + /* General geometrical parameters */ #define WINWIDTH 1280 /* window width */ @@ -48,40 +53,56 @@ #define YMIN -1.125 #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ -#define INITXMIN -1.95 -#define INITXMAX 1.95 /* x interval for initial condition */ -#define INITYMIN -1.05 -#define INITYMAX 0.95 /* y interval for initial condition */ +#define INITXMIN -1.85 +#define INITXMAX 1.85 /* x interval for initial condition */ +#define INITYMIN -0.9 +#define INITYMAX 0.9 /* y interval for initial condition */ + +#define BCXMIN -2.0 +#define BCXMAX 2.0 /* x interval for boundary condition */ +#define BCYMIN -1.125 +#define BCYMAX 1.125 /* y interval for boundary condition */ + +#define OBSXMIN -2.0 +#define OBSXMAX 2.0 /* x interval for motion of obstacle */ #define CIRCLE_PATTERN 8 /* pattern of circles, see list in global_ljones.c */ +#define ADD_FIXED_OBSTACLES 1 /* set to 1 do add fixed circular obstacles */ +#define OBSTACLE_PATTERN 0 /* pattern of obstacles, see list in global_ljones.c */ + #define TWO_TYPES 0 /* set to 1 to have two types of particles */ #define TPYE_PROPORTION 0.8 /* proportion of particles of first type */ -#define SYMMETRIZE_FORCE 1 /* set to 1 to symmetrize two-particle interaction, only needed if particles are not all the same */ +#define SYMMETRIZE_FORCE 0 /* set to 1 to symmetrize two-particle interaction, only needed if particles are not all the same */ +#define CENTER_PX 0 /* set to 1 to center horizontal momentum */ +#define CENTER_PY 0 /* set to 1 to center vertical momentum */ +#define CENTER_PANGLE 0 /* set to 1 to center angular momentum */ -#define INTERACTION 1 /* particle interaction, see list in global_ljones.c */ +#define INTERACTION 3 /* particle interaction, see list in global_ljones.c */ #define INTERACTION_B 1 /* particle interaction for second type of particle, see list in global_ljones.c */ -#define SPIN_INTER_FREQUENCY 1.0 /* angular frequency of spin-spin interaction */ +#define SPIN_INTER_FREQUENCY 5.0 /* angular frequency of spin-spin interaction */ #define SPIN_INTER_FREQUENCY_B 2.0 /* angular frequency of spin-spin interaction for second particle type */ #define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */ #define NPOISSON 100 /* number of points for Poisson C_RAND_POISSON arrangement */ -#define PDISC_DISTANCE 5.0 /* minimal distance in Poisson disc process, controls density of particles */ +#define PDISC_DISTANCE 1.8 /* minimal distance in Poisson disc process, controls density of particles */ #define PDISC_CANDIDATES 100 /* number of candidates in construction of Poisson disc process */ #define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */ #define LAMBDA 2.0 /* parameter controlling the dimensions of domain */ -#define MU 0.015 /* parameter controlling radius of particles */ +#define MU 0.045 /* parameter controlling radius of particles */ #define MU_B 0.02427051 /* parameter controlling radius of particles of second type */ #define NPOLY 3 /* number of sides of polygon */ -#define APOLY 0.5 /* angle by which to turn polygon, in units of Pi/2 */ +#define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */ #define MDEPTH 4 /* depth of computation of Menger gasket */ #define MRATIO 3 /* ratio defining Menger gasket */ #define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ #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 80 /* number of grid point for grid of disks */ -#define NGRIDY 25 /* number of grid point for grid of disks */ +#define NGRIDX 10 /* number of grid point for grid of disks */ +#define NGRIDY 3 /* 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 X_SHOOTER -0.2 #define Y_SHOOTER -0.6 @@ -90,13 +111,13 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 3100 /* number of frames of movie */ -#define NVID 150 /* number of iterations between images displayed on screen */ +#define NSTEPS 3750 /* number of frames of movie */ +#define NVID 100 /* number of iterations between images displayed on screen */ #define NSEG 150 /* number of segments of boundary */ -#define INITIAL_TIME 200 /* time after which to start saving frames */ +#define INITIAL_TIME 0 /* time after which to start saving frames */ #define BOUNDARY_WIDTH 1 /* width of particle boundary */ #define LINK_WIDTH 2 /* width of links between particles */ -#define CONTAINER_WIDTH 4 /* width of container boundary */ +#define CONTAINER_WIDTH 4 /* width of container boundary */ #define PAUSE 1000 /* number of frames after which to pause */ #define PSLEEP 1 /* sleep time during pause */ @@ -107,18 +128,18 @@ /* Boundary conditions, see list in global_ljones.c */ -#define BOUNDARY_COND 1 +#define BOUNDARY_COND 13 /* Plot type, see list in global_ljones.c */ -#define PLOT 0 -#define PLOT_B 6 /* plot type for second movie */ +#define PLOT 4 +#define PLOT_B 3 /* plot type for second movie */ -#define COLOR_BONDS 0 /* set to 1 to color bonds according to length */ +#define COLOR_BONDS 1 /* set to 1 to color bonds according to length */ /* Color schemes */ -#define COLOR_PALETTE 10 /* Color palette, see list in global_ljones.c */ +#define COLOR_PALETTE 0 /* Color palette, see list in global_ljones.c */ #define BLACK 1 /* background */ @@ -137,52 +158,55 @@ /* particle properties */ -#define PARTICLE_HUE_MIN 330.0 /* color of original particle */ -#define PARTICLE_HUE_MAX 50.0 /* color of saturated particle */ -#define PARTICLE_EMAX 1.0e3 /* max energy for particle to survive */ +#define ENERGY_HUE_MIN 330.0 /* color of original particle */ +#define ENERGY_HUE_MAX 50.0 /* color of saturated particle */ +#define PARTICLE_HUE_MIN 359.0 /* color of original particle */ +#define PARTICLE_HUE_MAX 0.0 /* color of saturated particle */ +#define PARTICLE_EMAX 5.0e2 /* energy of particle with hottest color */ #define RANDOM_RADIUS 0 /* set to 1 for random circle radius */ -#define DT_PARTICLE 2.0e-6 /* time step for particle displacement */ +#define DT_PARTICLE 1.0e-6 /* time step for particle displacement */ #define KREPEL 10.0 /* constant in repelling force between particles */ -#define EQUILIBRIUM_DIST 4.0 /* Lennard-Jones equilibrium distance for second type of particle */ +#define EQUILIBRIUM_DIST 6.0 /* Lennard-Jones equilibrium distance for second type of particle */ #define EQUILIBRIUM_DIST_B 5.0 /* Lennard-Jones equilibrium distance for second type of particle */ #define REPEL_RADIUS 20.0 /* radius in which repelling force acts (in units of particle radius) */ #define DAMPING 0.0 /* damping coefficient of particles */ #define PARTICLE_MASS 1.0 /* mass of particle of radius MU */ #define PARTICLE_MASS_B 0.1 /* mass of particle of radius MU */ -#define PARTICLE_INERTIA_MOMENT 0.05 /* moment of inertia of particle */ +#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.5 /* initial velocity range */ -#define OMEGA_INITIAL 2.0 /* initial angular velocity range */ +#define V_INITIAL 10.0 /* initial velocity range */ +#define OMEGA_INITIAL 10.0 /* initial angular velocity range */ -#define THERMOSTAT 0 /* set to 1 to switch on thermostat */ +#define THERMOSTAT 1 /* set to 1 to switch on thermostat */ #define SIGMA 5.0 /* noise intensity in thermostat */ -#define BETA 0.3 /* initial inverse temperature */ +#define BETA 0.002 /* initial inverse temperature */ #define MU_XI 0.05 /* friction constant in thermostat */ -#define KSPRING_BOUNDARY 2.0e7 /* confining harmonic potential outside simulation region */ -#define NBH_DIST_FACTOR 5.0 /* radius in which to count neighbours */ +#define KSPRING_BOUNDARY 5.0e8 /* confining harmonic potential outside simulation region */ +#define KSPRING_OBSTACLE 5.0e8 /* harmonic potential of obstacles */ +#define NBH_DIST_FACTOR 3.5 /* radius in which to count neighbours */ #define GRAVITY 0.0 /* gravity acting on all particles */ -#define ROTATION 0 /* set to 1 to include rotation of particles */ +#define ROTATION 1 /* set to 1 to include rotation of particles */ #define COUPLE_ANGLE_TO_THERMOSTAT 1 /* set to 1 to couple angular degrees of freedom to thermostat */ #define DIMENSION_FACTOR 1.0 /* scaling factor taking into account number of degrees of freedom */ -#define KTORQUE 7.0 /* force constant in angular dynamics */ +#define KTORQUE 600.0 /* force constant in angular dynamics */ #define KTORQUE_B 10.0 /* force constant in angular dynamics */ #define KTORQUE_DIFF 150.0 /* force constant in angular dynamics for different particles */ #define DRAW_SPIN 0 /* set to 1 to draw spin vectors of particles */ #define DRAW_SPIN_B 0 /* set to 1 to draw spin vectors of particles */ #define DRAW_CROSS 1 /* set to 1 to draw cross on particles of second type */ -#define SPIN_RANGE 5.0 /* range of spin-spin interaction */ +#define SPIN_RANGE 7.0 /* range of spin-spin interaction */ #define SPIN_RANGE_B 5.0 /* range of spin-spin interaction for second type of particle */ #define QUADRUPOLE_RATIO 0.6 /* anisotropy in quadrupole potential */ -#define INCREASE_BETA 0 /* set to 1 to increase BETA during simulation */ -#define BETA_FACTOR 0.001 /* factor by which to change BETA during simulation */ -#define N_TOSCILLATIONS 0.5 /* number of temperature oscillations in BETA schedule */ -#define NO_OSCILLATION 1 /* set to 1 to have exponential BETA change only */ -#define FINAL_CONSTANT_PHASE 200 /* final phase in which temperature is constant */ +#define INCREASE_BETA 1 /* set to 1 to increase BETA during simulation */ +#define BETA_FACTOR 500.0 /* factor by which to change BETA during simulation */ +#define N_TOSCILLATIONS 1.5 /* number of temperature oscillations in BETA schedule */ +#define NO_OSCILLATION 0 /* set to 1 to have exponential BETA change only */ +#define FINAL_CONSTANT_PHASE 0 /* final phase in which temperature is constant */ -#define DECREASE_CONTAINER_SIZE 1 /* set to 1 to decrease size of container */ +#define DECREASE_CONTAINER_SIZE 0 /* set to 1 to decrease size of container */ #define SYMMETRIC_DECREASE 0 /* set tp 1 to decrease container symmetrically */ #define COMPRESSION_RATIO 0.3 /* final size of container */ #define RESTORE_CONTAINER_SIZE 1 /* set to 1 to restore container to initial size at end of simulation */ @@ -190,9 +214,18 @@ #define MOVE_OBSTACLE 0 /* set to 1 to have a moving obstacle */ #define CENTER_VIEW_ON_OBSTACLE 0 /* set to 1 to center display on moving obstacle */ -#define OBSTACLE_RADIUS 0.27 /* radius of obstacle for circle boundary conditions */ -#define OBSTACLE_XMIN 0.0 /* initial position of obstacle */ -#define OBSTACLE_XMAX 3.0 /* final position of obstacle */ +#define RESAMPLE_Y 0 /* set to 1 to resample y coordinate of moved particles (for shock waves) */ +#define NTRIALS 2000 /* number of trials when resampling */ +#define OBSTACLE_RADIUS 0.15 /* radius of obstacle for circle boundary conditions */ +#define FUNNEL_WIDTH 0.25 /* funnel width for funnel boundary conditions */ +#define OBSTACLE_XMIN 0.0 /* initial position of obstacle */ +#define OBSTACLE_XMAX 3.0 /* final position of obstacle */ +#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 MAX_PRESSURE 3.0e10 /* pressure shown in "hottest" color */ +#define N_T_AVERAGE 50 /* size of temperature averaging window */ +#define PARTIAL_THERMO_COUPLING 1 /* set to 1 to couple only particles to the right of obstacle to thermostat */ #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 */ @@ -202,15 +235,33 @@ #define NPART_BOTTOM 100 /* number of particles at the bottom */ #define ADD_PARTICLES 0 /* set to 1 to add particles */ -#define ADD_TIME 100 /* time at which to add first particle */ -#define ADD_PERIOD 500 /* time interval between adding further particles */ -#define SAFETY_FACTOR 3.0 /* no particles are added at distance less than MU*SAFETY_FACTOR of other particles */ +#define ADD_TIME 25 /* time at which to add first particle */ +#define ADD_PERIOD 20 /* time interval between adding further particles */ +#define FINAL_NOADD_PERIOD 250 /* final period where no particles are added */ +#define SAFETY_FACTOR 2.0 /* no particles are added at distance less than MU*SAFETY_FACTOR of other particles */ -#define FLOOR_FORCE 1 /* set to 1 to limit force on particle to FMAX */ -#define FMAX 1.0e10 /* maximal force */ +#define TRACER_PARTICLE 0 /* set to 1 to have a tracer particle */ +#define TRAJECTORY_LENGTH 6000 /* length of recorded trajectory */ +#define TRACER_PARTICLE_MASS 0.1 /* relative mass of tracer particle */ +#define TRAJECTORY_WIDTH 3 /* width of tracer particle trajectory */ -#define HASHX 23 /* size of hashgrid in x direction */ -#define HASHY 14 /* size of hashgrid in y direction */ +#define POSITION_DEPENDENT_TYPE 0 /* set to 1 to make particle type depend on initial position */ +#define PRINT_ENTROPY 0 /* set to 1 to compute entropy */ + +#define PRINT_PARTICLE_NUMBER 0 /* set to 1 to print total number of particles */ + +#define EHRENFEST_COPY 0 /* set to 1 to add equal number of larger particles (for Ehrenfest model) */ + +#define LID_MASS 1000.0 /* mass of lid for BC_RECTANGLE_LID b.c. */ +#define LID_WIDTH 0.1 /* width of lid for BC_RECTANGLE_LID b.c. */ + +#define FLOOR_FORCE 0 /* set to 1 to limit force on particle to FMAX */ +#define FMAX 1.0e9 /* maximal force */ +#define FLOOR_OMEGA 0 /* set to 1 to limit particle momentum to PMAX */ +#define PMAX 10.0 /* maximal force */ + +#define HASHX 30 /* size of hashgrid in x direction */ +#define HASHY 20 /* size of hashgrid in y direction */ #define HASHMAX 100 /* maximal number of particles per hashgrid cell */ #define HASHGRID_PADDING 0.1 /* padding of hashgrid outside simulation window */ @@ -219,38 +270,17 @@ #define COLORBAR_RANGE_B 12.0 /* scale of color scheme bar for 2nd part */ #define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ -#include "global_ljones.c" -#include "sub_lj.c" +#define NO_WRAP_BC ((BOUNDARY_COND != BC_PERIODIC)&&(BOUNDARY_COND != BC_PERIODIC_CIRCLE)&&(BOUNDARY_COND != BC_PERIODIC_TRIANGLE)&&(BOUNDARY_COND != BC_KLEIN)&&(BOUNDARY_COND != BC_PERIODIC_FUNNEL)&&(BOUNDARY_COND != BC_BOY)) +#define PERIODIC_BC ((BOUNDARY_COND == BC_PERIODIC)||(BOUNDARY_COND == BC_PERIODIC_CIRCLE)||(BOUNDARY_COND == BC_PERIODIC_FUNNEL)||(BOUNDARY_COND == BC_PERIODIC_TRIANGLE)) double xshift = 0.0; /* x shift of shown window */ double xspeed = 0.0; /* x speed of obstacle */ +double ylid = 0.9; /* y coordinate of lid (for BC_RECTANGLE_LID b.c.) */ +double vylid = 0.0; /* y speed coordinate of lid (for BC_RECTANGLE_LID b.c.) */ -double 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; -} +#include "global_ljones.c" +#include "sub_lj.c" +#include "sub_hashgrid.c" /*********************/ @@ -258,998 +288,6 @@ double gaussian() /*********************/ - -void hash_xy_to_ij(double x, double y, int ij[2]) -{ - static int first = 1; - static double lx, ly, padding; - int i, j; - - if (first) - { - if ((BOUNDARY_COND == BC_PERIODIC)||(BOUNDARY_COND == BC_PERIODIC_CIRCLE)) padding = 0.0; - else padding = HASHGRID_PADDING; - lx = XMAX - XMIN + 2.0*padding; - ly = YMAX - YMIN + 2.0*padding; - first = 0; - } - - i = (int)((double)HASHX*(x - XMIN + padding)/lx); - j = (int)((double)HASHY*(y - YMIN + padding)/ly); - - if (i<0) i = 0; - else if (i>=HASHX) i = HASHX-1; - if (j<0) j = 0; - else if (j>=HASHY) j = HASHY-1; - - ij[0] = i; - ij[1] = j; - -// printf("Mapped (%.3f,%.3f) to (%i, %i)\n", x, y, ij[0], ij[1]); -} - - -double lennard_jones_force(double r, t_particle particle) -{ - int i; - double rmin = 0.01, rplus, ratio = 1.0; - - if (r > REPEL_RADIUS*particle.radius) return(0.0); - else - { - if (r > rmin) rplus = r; - else rplus = rmin; - -// ratio = ipow(particle.eq_dist*particle.radius/rplus, 6); - ratio = ipow(particle.eq_dist*particle.radius/rplus, 3); - ratio = ratio*ratio; - - return((ratio - 2.0*ratio*ratio)/rplus); - } -} - -void aniso_lj_force(double r, double ca, double sa, double ca_rel, double sa_rel, double force[2], t_particle particle) -{ - int i; - double rmin = 0.01, rplus, ratio = 1.0, c2, s2, c4, s4, a, aprime, f1, f2; - - if (r > REPEL_RADIUS*particle.radius) - { - force[0] = 0.0; - force[1] = 0.0; - } - else - { - if (r > rmin) rplus = r; - else rplus = rmin; - -// ratio = ipow(particle.eq_dist*particle.radius/rplus, 6); - ratio = ipow(particle.eq_dist*particle.radius/rplus, 3); - ratio = ratio*ratio; - - /* cos(2phi) and sin(2phi) */ - c2 = ca_rel*ca_rel - sa_rel*sa_rel; - s2 = 2.0*ca_rel*sa_rel; - - /* cos(4phi) and sin(4phi) */ - c4 = c2*c2 - s2*s2; - s4 = 2.0*c2*s2; - - a = 0.5*(9.0 - 7.0*c4); - aprime = 14.0*s4; - - f1 = ratio*(a - ratio)/rplus; - f2 = ratio*aprime/rplus; - - force[0] = f1*ca - f2*sa; - force[1] = f1*sa + f2*ca; - } -} - -void penta_lj_force(double r, double ca, double sa, double ca_rel, double sa_rel, double force[2], t_particle particle) -{ - int i; - double rmin = 0.01, rplus, ratio = 1.0, c2, s2, c4, s4, c5, s5, a, aprime, f1, f2; - static double a0, b0; - static int first = 1; - - if (first) - { - a0 = cos(0.1*PI) + 0.5; - b0 = a0 - 1.0; - first = 0; - } - - if (r > REPEL_RADIUS*particle.radius) - { - force[0] = 0.0; - force[1] = 0.0; - } - else - { - if (r > rmin) rplus = r; - else rplus = rmin; - -// ratio = ipow(particle.eq_dist*particle.radius/rplus, 6); - ratio = ipow(particle.eq_dist*particle.radius/rplus, 3); - ratio = ratio*ratio; - - - /* cos(2phi) and sin(2phi) */ - c2 = ca_rel*ca_rel - sa_rel*sa_rel; - s2 = 2.0*ca_rel*sa_rel; - - /* cos(4phi) and sin(4phi) */ - c4 = c2*c2 - s2*s2; - s4 = 2.0*c2*s2; - - /* cos(5phi) and sin(5phi) */ - c5 = ca_rel*c4 - sa_rel*s4; - s5 = sa_rel*c4 + ca_rel*s4; - - a = a0 - b0*c5; - aprime = 5.0*b0*s5; - - f1 = ratio*(a - ratio)/rplus; - f2 = ratio*aprime/rplus; - - force[0] = f1*ca - f2*sa; - force[1] = f1*sa + f2*ca; - } -} - -double old_golden_ratio_force(double r, t_particle particle) -/* potential with two minima at distances whose ratio is the golden ratio Phi */ -/* old version that does not work very well */ -{ - int i; - double x, y, z, rplus, ratio = 1.0, phi, a, phi3; - static int first = 1; - static double rmin, b, c, d; - - if (first) - { - rmin = 0.5*particle.radius; - phi = 0.5*(1.0 + sqrt(5.0)); - phi3 = 1.0/(phi*phi*phi); - a = 0.66; - b = 1.0 + phi3 + a; - d = phi3*a; - c = phi3 + a + d; -// b = 7.04; -// c = 13.66; -// d = 6.7; - first = 0; - printf("a = %.4lg, b = %.4lg, c = %.4lg, d = %.4lg\n", a, b, c, d); - } - - if (r > REPEL_RADIUS*particle.radius) return(0.0); - else - { - if (r > rmin) rplus = r; - else rplus = rmin; - - x = particle.eq_dist*particle.radius/rplus; - y = x*x*x; - z = d - c*y + b*y*y - y*y*y; - return(x*z/rplus); - } -} - -double golden_ratio_force(double r, t_particle particle) -/* potential with two minima at distances whose ratio is the golden ratio Phi */ -/* piecewise polynomial/LJ version */ -{ - int i; - double x, rplus, xm6, y1; - static int first = 1; - static double rmin, phi, a, h1, h2, phi6; - - if (first) - { - rmin = 0.5*particle.radius; - phi = 0.5*(1.0 + sqrt(5.0)); - a = 1.2; - - h1 = 1.0; /* inner potential well depth */ - h2 = 10.0; /* outer potential well depth */ - phi6 = ipow(phi, 6); - - first = 0; - } - - if (r > REPEL_RADIUS*particle.radius) return(0.0); - else - { - if (r > rmin) rplus = r; - else rplus = rmin; - - x = rplus/(particle.eq_dist*particle.radius); -// xm6 = 1.0/ipow(x, 6); - xm6 = 1.0/ipow(x, 3); - xm6 = xm6*xm6; - - if (x <= 1.0) return(12.0*h1*xm6*(1.0 - xm6)/x); - else if (x <= a) - { - y1 = ipow(a - 1.0, 3); - return(6.0*h1*(x - 1.0)*(a - x)/y1); - } - else if (x <= phi) - { - y1 = ipow(phi - a, 3); - return(6.0*h2*(x - a)*(x - phi)/y1); - } - else return(12.0*h2*phi6*(1.0 - phi6*xm6)*xm6/x); - } -} - -void dipole_lj_force(double r, double ca, double sa, double ca_rel, double sa_rel, double force[2], t_particle particle) -{ - int i; - double rmin = 0.01, rplus, ratio = 1.0, a, aprime, f1, f2; - - if (r > REPEL_RADIUS*MU) - { - force[0] = 0.0; - force[1] = 0.0; - } - else - { - if (r > rmin) rplus = r; - else rplus = rmin; - -// ratio = ipow(particle.eq_dist*particle.radius/rplus, 6); - ratio = ipow(particle.eq_dist*particle.radius/rplus, 3); - ratio = ratio*ratio; - - a = 1.0 + 0.25*ca_rel; - aprime = -0.25*sa_rel; - - f1 = ratio*(a - ratio)/rplus; - f2 = ratio*aprime/rplus; - - force[0] = f1*ca - f2*sa; - force[1] = f1*sa + f2*ca; - } -} - -void quadrupole_lj_force(double r, double ca, double sa, double ca_rel, double sa_rel, double force[2], t_particle particle) -{ - int i; - double rmin = 0.01, rplus, ratio = 1.0, a, aprime, f1, f2, ca2, sa2, x, y, dplus, dminus; - static int first = 1; - static double a0, b0, aplus, aminus; - - if (first) - { - dplus = cos(0.2*PI)*cos(0.1*PI); -// dminus = 0.8*dplus; - dminus = QUADRUPOLE_RATIO*dplus; - aplus = ipow(1.0/dplus, 6); - aminus = ipow(1.0/dminus, 6); -// aminus = ipow(cos(0.2*PI)*(0.25 + 0.5*sin(0.1*PI)), 6); - a0 = 0.5*(aplus + aminus); - b0 = 0.5*(aplus - aminus); - first = 0; - } - - if (r > REPEL_RADIUS*particle.radius) - { - force[0] = 0.0; - force[1] = 0.0; - } - else - { - if (r > rmin) rplus = r; - else rplus = rmin; - -// ratio = ipow(particle.eq_dist*particle.radius/rplus, 6); - ratio = ipow(particle.eq_dist*particle.radius/rplus, 3); - ratio = ratio*ratio; - - /* cos(2*phi) and sin(2*phi) */ - ca2 = ca_rel*ca_rel - sa_rel*sa_rel; - sa2 = 2.0*ca_rel*sa_rel; - - a = a0 + b0*ca2; -// if (a == 0.0) a = 1.0e-10; - aprime = -2.0*b0*sa2; - - f1 = ratio*(a - ratio)/rplus; - f2 = ratio*aprime/rplus; - - force[0] = f1*ca - f2*sa; - force[1] = f1*sa + f2*ca; - } -} - - -void quadrupole_lj_force2(double r, double ca, double sa, double ca_rel, double sa_rel, double force[2], t_particle particle) -{ - int i; - double rmin = 0.01, rplus, ratio = 1.0, a, aprime, f1, f2, ca2, sa2, x, y, eqdist; - static int first = 1; - static double aplus, aminus, a0, b0; - - if (first) - { - aplus = ipow(cos(0.2*PI)*cos(0.1*PI), 6); - aminus = 0.1*aplus; -// aminus = 0.0; -// aminus = -2.0*ipow(cos(0.2*PI)*(0.5*sin(0.1*PI)), 6); -// aminus = ipow(cos(0.2*PI)*(0.25 + 0.5*sin(0.1*PI)), 6); - a0 = 0.5*(aplus + aminus); - b0 = 0.5*(aplus - aminus); - first = 0; - } - - if (r > REPEL_RADIUS*particle.radius) - { - force[0] = 0.0; - force[1] = 0.0; - } - else - { - if (r > rmin) rplus = r; - else rplus = rmin; - - /* correct distance */ -// ratio = ipow(particle.eq_dist*particle.radius/rplus, 6); - ratio = ipow(particle.eq_dist*particle.radius/rplus, 3); - ratio = ratio*ratio; - - /* cos(2*phi) and sin(2*phi) */ - - ca2 = ca_rel*ca_rel - sa_rel*sa_rel; - sa2 = 2.0*ca_rel*sa_rel; - - a = a0 + b0*ca2; - if (a == 0.0) a = 1.0e-10; - aprime = -2.0*b0*sa2; - -// f1 = ratio*(a - ratio)/rplus; -// f2 = ratio*aprime/rplus; - - f1 = ratio*(aplus - ratio)/(rplus); - f2 = ratio*(aminus - ratio)/(rplus); - -// force[0] = f1*ca_rel - f2*sa_rel; -// force[1] = f1*sa_rel + f2*ca_rel; - - force[0] = f1*ca - f2*sa; - force[1] = f1*sa + f2*ca; - } -} - - - -int compute_repelling_force(int i, int j, double force[2], double *torque, t_particle* particle, - double krepel, int wrapx, int wrapy) -/* compute repelling force of particle j on particle i */ -/* returns 1 if distance between particles is smaller than NBH_DIST_FACTOR*MU */ -{ - double x1, y1, x2, y2, x3, y3, distance, r, f, angle, ca, sa, aniso, fx, fy, ff[2], cj, sj, ca_rel, sa_rel, dist_scaled, spin_f; - double xtemp, ytemp; - static double deltax = XMAX - XMIN, deltay = YMAX - YMIN, dxhalf = 0.5*(XMAX - XMIN), dyhalf = 0.5*(YMAX - YMIN); - int periodx, periody, wwrapx; -// static double factor = 0.5*(XMAX - XMIN); -// int interaction; - - x1 = particle[i].xc; - y1 = particle[i].yc; - x2 = particle[j].xc; - y2 = particle[j].yc; - - /* for monitoring purposes only */ - xtemp = x2; - ytemp = y2; - wwrapx = (BOUNDARY_COND == BC_KLEIN)&&(vabs(x2 - x1) > dxhalf); - - /* case of periodic boundary conditions */ - if ((BOUNDARY_COND == BC_PERIODIC)||(BOUNDARY_COND == BC_PERIODIC_CIRCLE)) - { - if (wrapy != 0) - { - if (y2 > 0.0) y2 -= deltay; - else y2 += deltay; - } - if (wrapx != 0) - { - if (x2 > 0.0) x2 -= deltax; - else x2 += deltax; - } - } - distance = module2(x2 - x1, y2 - y1); - - if ((distance == 0.0)||(i == j)) - { - force[0] = 0.0; - force[1] = 0.0; - *torque = 0.0; - return(1); - } - else if (distance > REPEL_RADIUS*particle[i].radius) - { - force[0] = 0.0; - force[1] = 0.0; - *torque = 0.0; - return(0); - } - else - { - ca = (x2 - x1)/distance; - sa = (y2 - y1)/distance; - - /* compute relative angle in case particles can rotate */ - if (ROTATION) - { - cj = cos(particle[j].angle); - sj = sin(particle[j].angle); - ca_rel = ca*cj + sa*sj; - sa_rel = sa*cj - ca*sj; - } - else - { - ca_rel = ca; - sa_rel = sa; - } - - switch (particle[j].interaction) { - case (I_COULOMB): - { - f = -krepel/(1.0e-8 + distance*distance); - force[0] = f*ca; - force[1] = f*sa; - break; - } - case (I_LENNARD_JONES): - { - f = krepel*lennard_jones_force(distance, particle[j]); - force[0] = f*ca; - force[1] = f*sa; - break; - } - case (I_LJ_DIRECTIONAL): - { - aniso_lj_force(distance, ca, sa, ca_rel, sa_rel, ff, particle[j]); - force[0] = krepel*ff[0]; - force[1] = krepel*ff[1]; - break; - } - case (I_LJ_PENTA): - { - penta_lj_force(distance, ca, sa, ca_rel, sa_rel, ff, particle[j]); - force[0] = krepel*ff[0]; - force[1] = krepel*ff[1]; - break; - } - case (I_GOLDENRATIO): - { - f = krepel*golden_ratio_force(distance, particle[j]); - force[0] = f*ca; - force[1] = f*sa; - break; - } - case (I_LJ_DIPOLE): - { - dipole_lj_force(distance, ca, sa, ca_rel, sa_rel, ff, particle[j]); - force[0] = krepel*ff[0]; - force[1] = krepel*ff[1]; - break; - } - case (I_LJ_QUADRUPOLE): - { - quadrupole_lj_force(distance, ca, sa, ca_rel, sa_rel, ff, particle[j]); - force[0] = krepel*ff[0]; - force[1] = krepel*ff[1]; - break; - } - } - } - - if (ROTATION) - { - dist_scaled = distance/(particle[i].spin_range*particle[i].radius); - spin_f = particle[i].spin_freq; - *torque = sin(spin_f*(particle[j].angle - particle[i].angle))/(1.0e-8 + dist_scaled*dist_scaled); - - if (particle[i].type == particle[j].type) - { - if (particle[i].type == 0) *torque *= KTORQUE; - else *torque *= KTORQUE_B; - } - else *torque *= KTORQUE_DIFF; - } -// if (ROTATION) *torque = -sin(particle[j].angle - particle[i].angle)/(1.0e-8 + distance*distance); -// if (ROTATION) *torque = ff[0]*sin(particle[i].angle) - ff[1]*cos(particle[i].angle); -// if (ROTATION) *torque = ff[0]*sin(particle[j].angle) - ff[1]*cos(particle[j].angle); - else *torque = 0.0; - - if ((distance < NBH_DIST_FACTOR*particle[i].radius)&&(j != i)) return(1); - else return(0); -} - - - -void update_hashgrid(t_particle* particle, int* hashgrid_number, int* hashgrid_particles) -{ - int i, j, k, n, m, ij[2], max = 0; - -// printf("Updating hashgrid_number\n"); - for (i=0; i max) max = n; -// printf("Placed particle %i at (%i,%i) in hashgrid\n", k, ij[0], ij[1]); -// printf("%i particles at (%i,%i)\n", hashgrid_number[ij[0]][ij[1]], ij[0], ij[1]); - } - - printf("Maximal number of particles per hash cell: %i\n", max); -} - -int add_particle(double x, double y, double vx, double vy, double mass, short int type, t_particle particle[NMAXCIRCLES]) -{ - int i, closeby = 0; - double dist; - - /* test distance to other particles */ - for (i=0; i= NMAXCIRCLES)) - { - printf("Cannot add particle at (%.3lg, %.3lg)\n", x, y); - return(0); - } - else - { - i = ncircles; - - particle[i].type = type; - - particle[i].xc = x; - particle[i].yc = y; - particle[i].radius = MU*sqrt(mass); - particle[i].active = 1; - particle[i].neighb = 0; - particle[i].thermostat = 1; - - particle[i].energy = 0.0; - - if (RANDOM_RADIUS) particle[i].radius = particle[i].radius*(0.75 + 0.5*((double)rand()/RAND_MAX)); - - particle[i].mass_inv = 1.0/mass; - if (particle[i].type == 0) particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT; - else particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT; - - particle[i].vx = vx; - particle[i].vy = vy; - particle[i].energy = (particle[i].vx*particle[i].vx + particle[i].vy*particle[i].vy)*particle[i].mass_inv; - - 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; - } - - ncircles++; - - printf("Added particle at (%.3lg, %.3lg)\n", x, y); - printf("Number of particles: %i\n", ncircles); - - return(1); - } -} - -double neighbour_color(int nnbg) -{ - if (nnbg > 7) nnbg = 7; - switch(nnbg){ - case (7): return(340.0); - case (6): return(310.0); - case (5): return(260.0); - case (4): return(200.0); - case (3): return(140.0); - case (2): return(100.0); - case (1): return(70.0); - default: return(30.0); - } -} - -void draw_one_particle(t_particle particle, double xc, double yc, double radius, double angle, int nsides, double width, double rgb[3]) -/* draw one of the particles */ -{ - double ca, sa, x1, x2, y1, y2, xc1; - - if (CENTER_VIEW_ON_OBSTACLE) - { - xc1 = xc - xshift; - if (xc1 < XMIN) xc1 += XMAX - XMIN; - else if (xc1 > XMAX) xc1 -= XMAX - XMIN; - } - else xc1 = xc; - glColor3f(rgb[0], rgb[1], rgb[2]); - if ((particle.interaction == I_LJ_QUADRUPOLE)||(particle.interaction == I_LJ_DIPOLE)) - draw_colored_rhombus(xc1, yc, radius, angle + APOLY*PID, rgb); - else draw_colored_polygon(xc1, yc, radius, nsides, angle + APOLY*PID, rgb); - - /* draw crosses on particles of second type */ - if ((TWO_TYPES)&&(DRAW_CROSS)) - if (particle.type == 1) - { - if (ROTATION) angle = angle + APOLY*PID; - else angle = APOLY*PID; - ca = cos(angle); - sa = sin(angle); - glLineWidth(3); - glColor3f(0.0, 0.0, 0.0); - x1 = xc1 - MU_B*ca; - y1 = yc - MU_B*sa; - x2 = xc1 + MU_B*ca; - y2 = yc + MU_B*sa; - draw_line(x1, y1, x2, y2); - x1 = xc1 - MU_B*sa; - y1 = yc + MU_B*ca; - x2 = xc1 + MU_B*sa; - y2 = yc - MU_B*ca; - draw_line(x1, y1, x2, y2); - } - - glLineWidth(width); - glColor3f(1.0, 1.0, 1.0); - if ((particle.interaction == I_LJ_QUADRUPOLE)||(particle.interaction == I_LJ_DIPOLE)) - draw_rhombus(xc1, yc, radius, angle + APOLY*PID); - else draw_polygon(xc1, yc, radius, nsides, angle + APOLY*PID); -} - -void draw_particles(t_particle particle[NMAXCIRCLES], int plot) -{ - int i, j, k, m, width, nnbg, nsides; - double ej, hue, rgb[3], radius, x1, y1, x2, y2, angle, ca, sa, length, linkcolor; - - blank(); - glColor3f(1.0, 1.0, 1.0); - - /* draw the bonds first */ - if (plot == P_BONDS) - { - glLineWidth(LINK_WIDTH); - for (j=0; j 0.5*(XMAX - XMIN)) x1 += XMAX - XMIN; - else if (x2 - x1 < -0.5*(XMAX - XMIN)) x1 -= XMAX - XMIN; - if (y2 - y1 > 0.5*(YMAX - YMIN)) y1 += YMAX - YMIN; - else if (y2 - y1 < -0.5*(YMAX - YMIN)) y1 -= YMAX - YMIN; - } - - if (COLOR_BONDS) - { - length = module2(x1 - x2, y1 - y2)/particle[j].radius; - if (length < 1.5) linkcolor = 1.0; - else linkcolor = 1.0 - 0.75*(length - 1.5)/(NBH_DIST_FACTOR - 1.5); -// printf("length = %.3lg\t color = %.3lg\n", length, linkcolor); - glColor3f(linkcolor, linkcolor, linkcolor); - } - - draw_line(x1, y1, x2, y2); - } - } - } - - /* determine particle color and size */ - for (j=0; j 0.0) - { - hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*ej/PARTICLE_EMAX; - if (hue > PARTICLE_HUE_MIN) hue = PARTICLE_HUE_MIN; - if (hue < PARTICLE_HUE_MAX) hue = PARTICLE_HUE_MAX; - } - radius = particle[j].radius; - width = BOUNDARY_WIDTH; - break; - } - case (P_NEIGHBOURS): - { - hue = neighbour_color(particle[j].neighb); - radius = particle[j].radius; - width = BOUNDARY_WIDTH; - break; - } - case (P_BONDS): - { -// if (particle[j].type == 1) hue = 70.0; /* to make second particle type more visible */ -// if (particle[j].type == 1) hue = neighbour_color(7 - particle[j].neighb); -// else - hue = neighbour_color(particle[j].neighb); - radius = particle[j].radius; - width = 1; - break; - } - case (P_ANGLE): - { - hue = particle[j].angle*particle[j].spin_freq/DPI; - hue -= (double)((int)hue); - hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*hue; - radius = particle[j].radius; - width = BOUNDARY_WIDTH; - break; - } - case (P_TYPE): - { - if (particle[j].type == 0) hue = 310.0; - else hue = 70.0; - radius = particle[j].radius; - width = BOUNDARY_WIDTH; - break; - } - case (P_DIRECTION): - { - hue = argument(particle[j].vx, particle[j].vy); - if (hue > DPI) hue -= DPI; - if (hue < 0.0) hue += DPI; - hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*hue/DPI; - radius = particle[j].radius; - width = BOUNDARY_WIDTH; - break; - } - } - - switch (particle[j].interaction) { - case (I_LJ_DIRECTIONAL): - { - nsides = 4; - break; - } - case (I_LJ_PENTA): - { - nsides = 5; - break; - } - case (I_LJ_QUADRUPOLE): - { - nsides = 4; - break; - } - default: nsides = NSEG; - } - - if (plot == P_BONDS) hsl_to_rgb_turbo(hue, 0.9, 0.5, rgb); - else if (plot == P_DIRECTION) hsl_to_rgb_twilight(hue, 0.9, 0.5, rgb); - else hsl_to_rgb(hue, 0.9, 0.5, rgb); - angle = particle[j].angle + APOLY*DPI; - - draw_one_particle(particle[j], particle[j].xc, particle[j].yc, radius, angle, nsides, width, rgb); - - /* in case of periodic b.c., draw translates of particles */ - if ((BOUNDARY_COND == BC_PERIODIC)||(BOUNDARY_COND == BC_PERIODIC_CIRCLE)) - { - x1 = particle[j].xc; - y1 = particle[j].yc; - - for (i=-2; i<3; i++) - for (k=-1; k<2; k++) - draw_one_particle(particle[j], x1 + (double)i*(XMAX - XMIN), y1 + (double)k*(YMAX - YMIN), radius, angle, nsides, width, rgb); - } - } - - /* draw spin vectors */ - if ((DRAW_SPIN)||(DRAW_SPIN_B)) - { - glLineWidth(width); - for (j=0; j 1.0) logratio = 1.0; - if (BETA_FACTOR > 1.0) hue = PARTICLE_HUE_MAX - (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*logratio; - else hue = PARTICLE_HUE_MIN - (PARTICLE_HUE_MIN - PARTICLE_HUE_MAX)*logratio; - erase_area_hsl_turbo(xbox, y + 0.025, 0.37, 0.05, hue, 0.9, 0.5); - if ((hue < 90)||(hue > 270)) glColor3f(1.0, 1.0, 1.0); - else glColor3f(0.0, 0.0, 0.0); - sprintf(message, "Temperature %.3f", 1.0/beta); - write_text(xtext, y, message); - y -= 0.1; - } - if (DECREASE_CONTAINER_SIZE) /* print density */ - { - density = (double)ncircles/((lengthcontainer)*(INITYMAX - INITYMIN)); - erase_area_hsl(xbox, y + 0.025, 0.37, 0.05, 0.0, 0.9, 0.0); - glColor3f(1.0, 1.0, 1.0); - sprintf(message, "Density %.3f", density); - write_text(xtext, y, message); - - erase_area_hsl(xmid, y + 0.025, 0.37, 0.05, 0.0, 0.9, 0.0); - glColor3f(1.0, 1.0, 1.0); - sprintf(message, "Temperature %.3f", temperature); - write_text(xmidtext, y, message); - - erase_area_hsl(xxbox, y + 0.025, 0.37, 0.05, 0.0, 0.9, 0.0); - glColor3f(1.0, 1.0, 1.0); - sprintf(message, "Pressure %.3f", meanpressure); - write_text(xxtext, y, message); - - } - else if (INCREASE_KREPEL) /* print force constant */ - { - erase_area_hsl(xbox, y + 0.025, 0.22, 0.05, 0.0, 0.9, 0.0); - glColor3f(1.0, 1.0, 1.0); - sprintf(message, "Force %.0f", krepel); - write_text(xtext + 0.28, y, message); - } -} - - double repel_schedule(int i) { static double kexponent; @@ -1316,7 +354,7 @@ double obstacle_schedule_old(int i) double obstacle_schedule(int i) { - double time; + double time, acceleration = 40.0; double x; // static double t1 = 0.5, t2 = 0.75, t3 = 0.875; @@ -1325,272 +363,218 @@ double obstacle_schedule(int i) { time = (double)(i-INITIAL_TIME)/(double)(NSTEPS); - x = OBSTACLE_XMIN + 50.0*time*time; - xspeed = 100.0*time; - while (x > XMAX) x += XMIN - XMAX; + x = OBSTACLE_XMIN + 0.5*acceleration*time*time; + xspeed = acceleration*time; +// while (x > OBSXMAX) x += OBSXMIN - OBSXMAX; return(x); } } -double compute_boundary_force(t_particle particle, double *fx, double *fy, double xleft, double xright) +double obstacle_schedule_smooth(int i, int j) { - int i; - double xmin, xmax, ymin, ymax, padding, r, cphi, sphi, f, fperp = 0.0, x; + double time, acceleration = 50.0; + double x; +// static double t1 = 0.5, t2 = 0.75, t3 = 0.875; - switch(BOUNDARY_COND){ - case (BC_SCREEN): - { - /* add harmonic force outside screen */ - if (particle.xc > XMAX) *fx -= KSPRING_BOUNDARY*(particle.xc - XMAX); - else if (particle.xc < XMIN) *fx += KSPRING_BOUNDARY*(XMIN - particle.xc); - if (particle.yc > YMAX) *fy -= KSPRING_BOUNDARY*(particle.yc - YMAX); - else if (particle.yc < YMIN) *fy += KSPRING_BOUNDARY*(YMIN - particle.yc); - return(fperp); - } - case (BC_RECTANGLE): - { - /* add harmonic force outside rectangular box */ - padding = MU + 0.01; - xmin = xleft + padding; - xmax = xright - padding; - ymin = INITYMIN + padding; - ymax = INITYMAX - padding; - - if (particle.xc > xmax) - { - fperp = KSPRING_BOUNDARY*(particle.xc - xmax); - *fx -= fperp; - } - else if (particle.xc < xmin) - { - fperp = KSPRING_BOUNDARY*(xmin - particle.xc); - *fx += fperp; - } - if (particle.yc > ymax) - { - fperp = KSPRING_BOUNDARY*(particle.yc - ymax); - *fy -= fperp; - } - else if (particle.yc < ymin) - { - fperp = KSPRING_BOUNDARY*(ymin - particle.yc); - *fy += fperp; - } -// if (particle.xc > xmax) *fx -= KSPRING_BOUNDARY*(particle.xc - xmax); -// else if (particle.xc < xmin) *fx += KSPRING_BOUNDARY*(xmin - particle.xc); -// if (particle.yc > ymax) *fy -= KSPRING_BOUNDARY*(particle.yc - ymax); -// else if (particle.yc < ymin) *fy += KSPRING_BOUNDARY*(ymin - particle.yc); - - return(fperp); - } - case (BC_CIRCLE): - { - /* add harmonic force outside screen */ - if (particle.xc > XMAX) *fx -= KSPRING_BOUNDARY*(particle.xc - XMAX); - else if (particle.xc < XMIN) *fx += KSPRING_BOUNDARY*(XMIN - particle.xc); - if (particle.yc > YMAX) *fy -= KSPRING_BOUNDARY*(particle.yc - YMAX); - else if (particle.yc < YMIN) *fy += KSPRING_BOUNDARY*(YMIN - particle.yc); - - /* add harmonic force from obstacle */ - for (i=-1; i<2; i++) - { - x = xleft + (double)i*(XMAX - XMIN); - if (vabs(particle.xc - x) < 1.1*OBSTACLE_RADIUS) - { - r = module2(particle.xc - x, particle.yc); - if (r < 1.0e-5) r = 1.0e-05; - cphi = (particle.xc - x)/r; - sphi = particle.yc/r; - padding = MU + 0.03; - if (r < OBSTACLE_RADIUS + padding) - { - f = KSPRING_BOUNDARY*(OBSTACLE_RADIUS + padding - r); - *fx += f*cphi; - *fy += f*sphi; - } - } - } - return(fperp); - } - case (BC_PERIODIC_CIRCLE): - { - /* add harmonic force from obstacle */ - for (i=-1; i<2; i++) - { - x = xleft + (double)i*(XMAX - XMIN); - if (vabs(particle.xc - x) < 1.1*OBSTACLE_RADIUS) - { - r = module2(particle.xc - x, particle.yc); - if (r < 1.0e-5) r = 1.0e-05; - cphi = (particle.xc - x)/r; - sphi = particle.yc/r; - padding = MU + 0.03; - if (r < OBSTACLE_RADIUS + padding) - { - f = KSPRING_BOUNDARY*(OBSTACLE_RADIUS + padding - r); - *fx += f*cphi; - *fy += f*sphi; - } - } - } - return(fperp); - } + if (i < INITIAL_TIME) return(OBSTACLE_XMIN); + else + { + time = ((double)(i-INITIAL_TIME) + (double)j/(double)NVID)/(double)(NSTEPS); + + x = OBSTACLE_XMIN + 0.5*acceleration*time*time; + xspeed = acceleration*time; +// while (x > OBSXMAX) x += OBSXMIN - OBSXMAX; + return(x); } } +double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HASHX*HASHY], + double qx[NMAXCIRCLES], double qy[NMAXCIRCLES], double qangle[NMAXCIRCLES], + double px[NMAXCIRCLES], double py[NMAXCIRCLES], double pangle[NMAXCIRCLES], + double beta, int *nactive, int *nsuccess, int *nmove) +{ + double a, totalenergy = 0.0; + static double b = 0.25*SIGMA*SIGMA*DT_PARTICLE/MU_XI, xi = 0.0; + int j, move; + + #pragma omp parallel for private(j,xi,totalenergy,a,move) + for (j=0; j BCXMAX) + else if (particle[j].xc > xshift + BCXMAX) + { + particle[j].xc += BCXMIN - BCXMAX; + } + if (particle[j].yc > BCYMAX) particle[j].yc += BCYMIN - BCYMAX; + else if (particle[j].yc < BCYMIN) particle[j].yc += BCYMAX - BCYMIN; + } + else if (!NO_WRAP_BC) + { + move += wrap_particle(&particle[j], &px[j], &py[j]); + } +// if (move > 0) +// { +// compute_relative_positions(particle, hashgrid); +// update_hashgrid(particle, hashgrid, 0); /* REDUNDANT ? */ +// } + } + + return(totalenergy); +} + + +void evolve_lid(double fboundary) +{ + double force; + + force = fboundary - GRAVITY*LID_MASS; + if (ylid > BCYMAX + LID_WIDTH) force -= KSPRING_BOUNDARY*(ylid - BCYMAX - LID_WIDTH); +// printf("Force on lid = %.3lg\n", force); + vylid += force*DT_PARTICLE/LID_MASS; + ylid += vylid*DT_PARTICLE; +// if (ylid > BCYMAX + LID_WIDTH) ylid = BCYMAX; +} + 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 = INITXMIN, xmaxcontainer = INITXMAX, torque, torque_ij, - fboundary = 0.0; - double *qx, *qy, *px, *py, *qangle, *pangle; - int i, j, k, n, m, s, ij[2], i0, iplus, iminus, j0, jplus, jminus, p, q, p1, q1, total_neighbours = 0, - min_nb, max_nb, close, wrapx = 0, wrapy = 0, nactive = 0, nadd_particle = 0; + fboundary = 0.0, pleft = 0.0, pright = 0.0, entropy[2], mean_energy; + double *qx, *qy, *px, *py, *qangle, *pangle, *pressure; + 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, + tracer_n, traj_position = 0, traj_length = 0, move = 0, old, m0, floor, nthermo; static int imin, imax; static short int first = 1; t_particle *particle; - int *hashgrid_number, *hashgrid_particles; + t_obstacle *obstacle; + t_tracer *trajectory; t_hashgrid *hashgrid; char message[100]; particle = (t_particle *)malloc(NMAXCIRCLES*sizeof(t_particle)); /* particles */ - + if (ADD_FIXED_OBSTACLES) obstacle = (t_obstacle *)malloc(NMAXOBSTACLES*sizeof(t_obstacle)); /* obstacles */ + + if (TRACER_PARTICLE) trajectory = (t_tracer *)malloc(TRAJECTORY_LENGTH*sizeof(t_tracer)); + hashgrid = (t_hashgrid *)malloc(HASHX*HASHY*sizeof(t_hashgrid)); /* hashgrid */ - - hashgrid_number = (int *)malloc(HASHX*HASHY*sizeof(int)); /* total number of particles in each hash grid cell */ - hashgrid_particles = (int *)malloc(HASHX*HASHY*HASHMAX*sizeof(int)); /* numbers of particles in each hash grid cell */ - + qx = (double *)malloc(NMAXCIRCLES*sizeof(double)); qy = (double *)malloc(NMAXCIRCLES*sizeof(double)); px = (double *)malloc(NMAXCIRCLES*sizeof(double)); py = (double *)malloc(NMAXCIRCLES*sizeof(double)); qangle = (double *)malloc(NMAXCIRCLES*sizeof(double)); pangle = (double *)malloc(NMAXCIRCLES*sizeof(double)); + pressure = (double *)malloc(N_PRESSURES*sizeof(double)); /* initialise positions and radii of circles */ init_particle_config(particle); - - - /* initialise particles */ - for (i=0; i < NMAXCIRCLES; i++) - { - /* set particle type */ - particle[i].type = 0; - if ((TWO_TYPES)&&((double)rand()/RAND_MAX > TPYE_PROPORTION)) - { - particle[i].type = 1; - particle[i].radius = MU_B; - } - - particle[i].neighb = 0; - particle[i].thermostat = 1; - -// particle[i].energy = 0.0; - y = particle[i].yc; - if (y >= YMAX) y -= particle[i].radius; - if (y <= YMIN) y += particle[i].radius; - - if (RANDOM_RADIUS) particle[i].radius = particle[i].radius*(0.75 + 0.5*((double)rand()/RAND_MAX)); - - if (particle[i].type == 0) - { - particle[i].interaction = INTERACTION; - particle[i].eq_dist = EQUILIBRIUM_DIST; - particle[i].spin_range = SPIN_RANGE; - particle[i].spin_freq = SPIN_INTER_FREQUENCY; - particle[i].mass_inv = 1.0/PARTICLE_MASS; - particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT; - } - else - { - 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; - particle[i].mass_inv = 1.0/PARTICLE_MASS_B; - particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT_B; - } - - particle[i].vx = V_INITIAL*gaussian(); - particle[i].vy = V_INITIAL*gaussian(); - particle[i].energy = (particle[i].vx*particle[i].vx + particle[i].vy*particle[i].vy)*particle[i].mass_inv; - - px[i] = particle[i].vx; - py[i] = particle[i].vy; - - if (ROTATION) - { - particle[i].angle = DPI*(double)rand()/RAND_MAX; - particle[i].omega = OMEGA_INITIAL*gaussian(); - if (COUPLE_ANGLE_TO_THERMOSTAT) - particle[i].energy += particle[i].omega*particle[i].omega*particle[i].inertia_moment_inv; - } - else - { - particle[i].angle = 0.0; - particle[i].omega = 0.0; - } - pangle[i] = particle[i].omega; - } - /* initialize dummy values in case particles are added */ - for (i=ncircles; i < NMAXCIRCLES; i++) - { - particle[i].type = 0; - particle[i].active = 0; - particle[i].neighb = 0; - particle[i].thermostat = 0; - particle[i].energy = 0.0; - particle[i].mass_inv = 1.0/PARTICLE_MASS; - particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT; - particle[i].vx = 0.0; - particle[i].vy = 0.0; - px[i] = 0.0; - py[i] = 0.0; - particle[i].angle = DPI*(double)rand()/RAND_MAX; - particle[i].omega = 0.0; - pangle[i] = 0.0; - particle[i].interaction = INTERACTION; - particle[i].eq_dist = EQUILIBRIUM_DIST; - particle[i].spin_range = SPIN_RANGE; - particle[i].spin_freq = SPIN_INTER_FREQUENCY; - } + init_hashgrid(hashgrid); - /* add particles at the bottom as seed */ - if (PART_AT_BOTTOM) for (i=0; i<=NPART_BOTTOM; i++) - { - x = XMIN + (double)i*(XMAX - XMIN)/(double)NPART_BOTTOM; - y = YMIN + 2.0*MU; - add_particle(x, y, 0.0, 0.0, MASS_PART_BOTTOM, 0, particle); - } - if (PART_AT_BOTTOM) for (i=0; i<=NPART_BOTTOM; i++) - { - x = XMIN + (double)i*(XMAX - XMIN)/(double)NPART_BOTTOM; - y = YMIN + 4.0*MU; - add_particle(x, y, 0.0, 0.0, MASS_PART_BOTTOM, 0, particle); - } + xshift = OBSTACLE_XMIN; + + if (ADD_FIXED_OBSTACLES) init_obstacle_config(obstacle); - /* inactivate particles in obstacle */ - if ((BOUNDARY_COND == BC_CIRCLE)||(BOUNDARY_COND == BC_PERIODIC_CIRCLE)) - for (i=0; i< ncircles; i++) - if ((module2(particle[i].xc - OBSTACLE_XMIN, particle[i].yc) < 1.2*OBSTACLE_RADIUS)) - particle[i].active = 0; - - /* count number of active particles */ - for (i=0; i< ncircles; i++) nactive += particle[i].active; + if (RECORD_PRESSURES) for (i=0; i XMAX) xshift -= XMAX - XMIN; - } blank(); fboundary = 0.0; + pleft = 0.0; + pright = 0.0; + if (RECORD_PRESSURES) for (j=0; j= HASHX) iplus = HASHX-1; - - if (jminus < 0) jminus = 0; - else if (jplus >= HASHY) jplus = HASHY-1; - } - - - fx = 0.0; - fy = 0.0; - torque = 0.0; - for (p=iminus; p<= iplus; p++) - for (q=jminus; q<= jplus; q++) - { - p1 = p; - q1 = q; - wrapx = 0; - wrapy = 0; - /* case of periodic boundary conditions */ - if ((BOUNDARY_COND == BC_PERIODIC)||(BOUNDARY_COND == BC_PERIODIC_CIRCLE)) - { - if (p1 < 0) - { - p1 = HASHX-1; - wrapx = -1; - } - else if (p1 >= HASHX) - { - p1 = 0; - wrapx = 1; - } - if (q1 < 0) - { - q1 = HASHY-1; - wrapy = -1; - } - else if (q1 >= HASHY) - { - q1 = 0; - wrapy = 1; - } - } - for (k=0; k FMAX) fx = FMAX; - if (fx < -FMAX) fx = -FMAX; - if (fy > FMAX) fy = FMAX; - if (fy < -FMAX) fy = -FMAX; + if (particle[j].fx > FMAX) particle[j].fx = FMAX; + if (particle[j].fx < -FMAX) particle[j].fx = -FMAX; + if (particle[j].fy > FMAX) particle[j].fy = FMAX; + if (particle[j].fy < -FMAX) particle[j].fy = -FMAX; + if (particle[j].torque > FMAX) particle[j].torque = FMAX; + if (particle[j].torque < -FMAX) particle[j].torque = -FMAX; } - - particle[j].fx = fx; - particle[j].fy = fy; - particle[j].torque = torque; - -// if (j%500 == 1) printf("Particle %i angle = %.5lg omega = %.5lg torque = %.5lg\n", j, particle[j].angle, particle[j].omega, particle[j].torque); -// if (j%500 == 1) printf("Particle %i x = %.5lg vx = %.5lg fx = %.5lg\n", j, particle[j].xc, particle[j].vx, particle[j].fx); } /* timestep of thermostat algorithm */ - for (j=0; jN_T_AVERAGE)) + { + nthermo = partial_thermostat_coupling(particle, xshift + OBSTACLE_RADIUS); + printf("%i particles coupled to thermostat out of %i active\n", nthermo, nactive); + mean_energy = compute_mean_energy(particle); } + else mean_energy = totalenergy/(double)ncircles; + + if (CENTER_PX) center_momentum(px); + if (CENTER_PY) center_momentum(py); + if (CENTER_PANGLE) center_momentum(pangle); + if (FLOOR_OMEGA) floor = floor_momentum(pangle); + +// printf("pressure left %.5lg, pressure right %.5lg\n", pleft, pright); +// for (j=0; j XMAX) particle[j].xc += XMIN - XMAX; - else if (particle[j].xc < XMIN) particle[j].xc += XMAX - XMIN; - if (particle[j].yc > YMAX) particle[j].yc += YMIN - YMAX; - else if (particle[j].yc < YMIN) particle[j].yc += YMAX - YMIN; - } + + /* update tracer particle trajectory */ + if ((TRACER_PARTICLE)&&(i > INITIAL_TIME)) + { + trajectory[traj_position].xc = particle[tracer_n].xc; + trajectory[traj_position].yc = particle[tracer_n].yc; + traj_position++; + if (traj_position >= TRAJECTORY_LENGTH) traj_position = 0; + traj_length++; + if (traj_length >= TRAJECTORY_LENGTH) traj_length = TRAJECTORY_LENGTH - 1; +// for (j=0; j max_nb) max_nb = particle[j].neighb; if (particle[j].neighb < min_nb) min_nb = particle[j].neighb; } - printf("Mean number of neighbours: %.3f\n", (double)total_neighbours/(double)ncircles); - printf("Min number of neighbours: %i\n", min_nb); - printf("Max number of neighbours: %i\n", max_nb); +// printf("Mean number of neighbours: %.3f\n", (double)total_neighbours/(double)ncircles); +// printf("Min number of neighbours: %i\n", min_nb); +// printf("Max number of neighbours: %i\n", max_nb); + if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length); draw_particles(particle, PLOT); - draw_container(xmincontainer, xmaxcontainer); + draw_container(xmincontainer, xmaxcontainer, obstacle); /* add a particle */ - if ((ADD_PARTICLES)&&((i - INITIAL_TIME - ADD_TIME + 1)%ADD_PERIOD == 0)) + if ((ADD_PARTICLES)&&((i - INITIAL_TIME - ADD_TIME + 1)%ADD_PERIOD == 0)&&(i < NSTEPS - FINAL_NOADD_PERIOD)) + nadd_particle = add_particles(particle, px, py, nadd_particle); + + update_hashgrid(particle, hashgrid, 1); + + print_parameters(beta, mean_energy, krepel, xmaxcontainer - xmincontainer, + fboundary/(double)(ncircles*NVID), 1, pressure); + if (BOUNDARY_COND == BC_EHRENFEST) print_ehrenfest_parameters(particle, pleft, pright); + else if (PRINT_PARTICLE_NUMBER) print_particle_number(ncircles); + + if (PRINT_ENTROPY) { -// add_particle(XMIN + 0.1, 0.0, 50.0, 0.0, 3.0, 0, particle); -// px[ncircles - 1] = particle[ncircles - 1].vx; -// py[ncircles - 1] = particle[ncircles - 1].vy; -// particle[ncircles - 1].radius = 1.5*MU; -// j = 0; -// while (module2(particle[j].xc,particle[j].yc) > 0.7) j = rand()%ncircles; -// x = particle[j].xc + 2.5*MU; -// y = particle[j].yc; - -// x = XMIN + (XMAX - XMIN)*rand()/RAND_MAX; -// y = YMAX + 0.01*rand()/RAND_MAX; -// add_particle(x, y, 0.0, 0.0, 1.0, 0, particle); - - x = XMIN + 0.25*(XMAX - XMIN); - y = YMAX + 0.01; - prop = 1.0 - (double)nadd_particle/5.0; - vx = 100.0*prop; - add_particle(x, y, vx, -10.0, 5.0*prop, 0, particle); - particle[ncircles - 1].radius = 10.0*MU*prop; - particle[ncircles - 1].eq_dist = 2.0; - particle[ncircles - 1].thermostat = 0; - px[ncircles - 1] = particle[ncircles - 1].vx; - py[ncircles - 1] = particle[ncircles - 1].vy; - nadd_particle++; + compute_entropy(particle, entropy); + printf("Entropy 1 = %.5lg, Entropy 2 = %.5lg\n", entropy[0], entropy[1]); + print_entropy(entropy); } - - update_hashgrid(particle, hashgrid_number, hashgrid_particles); - - print_parameters(beta, totalenergy/(double)ncircles, krepel, xmaxcontainer - xmincontainer, - fboundary/(double)(ncircles*NVID), 0); glutSwapBuffers(); if (MOVIE) { - if (i >= INITIAL_TIME) save_frame_lj(); + if (i >= INITIAL_TIME) + { + save_frame_lj(); + if ((TIME_LAPSE)&&((i - INITIAL_TIME)%TIME_LAPSE_FACTOR == 0)&&(!DOUBLE_MOVIE)) + { + save_frame_lj_counter(NSTEPS + END_FRAMES + (i - INITIAL_TIME)/TIME_LAPSE_FACTOR); + } + } else printf("Initial phase time %i of %i\n", i, INITIAL_TIME); + if ((i >= INITIAL_TIME)&&(DOUBLE_MOVIE)) { + if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length); draw_particles(particle, PLOT_B); - draw_container(xmincontainer, xmaxcontainer); - print_parameters(beta, totalenergy/(double)ncircles, krepel, xmaxcontainer - xmincontainer, - fboundary/(double)(ncircles*NVID), 0); + draw_container(xmincontainer, xmaxcontainer, obstacle); + print_parameters(beta, mean_energy, krepel, xmaxcontainer - xmincontainer, + fboundary/(double)(ncircles*NVID), 0, pressure); + if (BOUNDARY_COND == BC_EHRENFEST) print_ehrenfest_parameters(particle, pleft, pright); + else if (PRINT_PARTICLE_NUMBER) print_particle_number(ncircles); glutSwapBuffers(); - save_frame_lj_counter(NSTEPS + 21 + counter); + save_frame_lj_counter(NSTEPS + MID_FRAMES + 1 + counter); counter++; } @@ -1903,35 +770,47 @@ void animation() { if (DOUBLE_MOVIE) { + if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length); draw_particles(particle, PLOT); - draw_container(xmincontainer, xmaxcontainer); - print_parameters(beta, totalenergy/(double)ncircles, krepel, xmaxcontainer - xmincontainer, - fboundary/(double)(ncircles*NVID), 0); + draw_container(xmincontainer, xmaxcontainer, obstacle); + print_parameters(beta, mean_energy, krepel, xmaxcontainer - xmincontainer, + fboundary/(double)(ncircles*NVID), 0, pressure); + if (BOUNDARY_COND == BC_EHRENFEST) print_ehrenfest_parameters(particle, pleft, pright); + else if (PRINT_PARTICLE_NUMBER) print_particle_number(ncircles); glutSwapBuffers(); } for (i=0; i=HASHX) i = HASHX-1; + if (j<0) j = 0; + else if (j>=HASHY) j = HASHY-1; + + return(mhash(i,j)); +// printf("Mapped (%.3f,%.3f) to (%i, %i)\n", x, y, ij[0], ij[1]); +} + + +void init_hashgrid(t_hashgrid hashgrid[HASHX*HASHY]) +/* initialise table of neighbouring cells for each hashgrid cell, depending on boundary condition */ +{ + int i, j, k, p, q, m, i1, j1; + + printf("Initializing hash grid\n"); + + /* bulk of the table */ + for (i=0; i= 0))||((i == HASHX-1)&&(p <= 0))) j1 = (j+q+HASHY)%HASHY; + else j1 = (2*HASHY-1-j-q)%HASHY; + hashgrid[mhash(i,j)].neighbour[3*(p+1)+q+1] = mhash(i1, j1); + } + } + + /* top/bottom boundaries */ + for (i=1; i= 0))||((i == HASHX-1)&&(p <= 0))) j1 = (j+q+HASHY)%HASHY; + else j1 = (2*HASHY-1-j-q)%HASHY; + hashgrid[mhash(i,j)].neighbour[3*(p+1)+q+1] = mhash(i1, j1); + } + } + + /* top/bottom boundaries */ + for (i=1; i= 0))||((j == HASHY-1)&&(q <= 0))) i1 = (i+p+HASHX)%HASHX; + else i1 = (2*HASHX-1-i-p)%HASHX; + j1 = (j+q+HASHY)%HASHY; + hashgrid[mhash(i,j)].neighbour[3*(p+1)+q+1] = mhash(i1, j1); + } + } + + /* corners */ + for (i=0; i max) max = n; + } + + if(verbose) printf("Maximal number of particles per hash cell: %i\n", max); +} + + +int wrap_particle(t_particle* particle, double *px, double *py) +/* relocate particles in case of periodic and similar boundary conditions */ +{ + double x, y, x1, y1, x2, y2; + int move = 0; + + x = particle->xc; + y = particle->yc; + x1 = x; + y1 = y; + + switch (bc_grouped(BOUNDARY_COND)) { + case (0): /* rectangular b.c. */ + { + /* do nothing */ + return(0); + break; + } + case (1): /* periodic b.c. */ + { + if (x < BCXMIN) + { + x1 += BCXMAX - BCXMIN; + move++; + } + else if (x > BCXMAX) + { + x1 += BCXMIN - BCXMAX; + move++; + } + if (y > BCYMAX) + { + y1 += BCYMIN - BCYMAX; + move++; + } + else if (y < BCYMIN) + { + y1 += BCYMAX - BCYMIN; + move++; + } + particle->xc = x1; + particle->yc = y1; + return(move); + break; + } + case (2): /* Klein bottle b.c. */ + { + if (y > BCYMAX) + { + y1 += BCYMIN - BCYMAX; + move++; + } + else if (y < BCYMIN) + { + y1 += BCYMAX - BCYMIN; + move++; + } + if (x < BCXMIN) + { + x1 += BCXMAX - BCXMIN; + y1 = -y1; + *py *= -1.0; + move++; + } + else if (x > BCXMAX) + { + x1 += BCXMIN - BCXMAX; + y1 = -y1; + *py *= -1.0; + move++; + } + particle->xc = x1; + particle->yc = y1; + return(move); + break; + } + case (3): /* Boy surface b.c. */ + { + if ((y < BCYMAX)&&(y > BCYMIN)) + { + if (x > BCXMAX) + { + x1 += BCXMIN - BCXMAX; + y1 = -y1; + particle->angle *= -1.0; + particle->vy *= -1.0; + *py *= -1.0; + move++; + } + else if (x < BCXMIN) + { + x1 += BCXMAX - BCXMIN; + y1 = -y1; + particle->angle *= -1.0; + particle->vy *= -1.0; + *py *= -1.0; + move++; + } + } +// x = x1; +// y = y1; + if ((x < BCXMAX)&&(x > BCXMIN)) + { + if (y > BCYMAX) + { + y1 += BCYMIN - BCYMAX; + x1 = -x1; + particle->angle = PI - particle->angle; + particle->vx *= -1.0; + *px *= -1.0; + move++; + } + else if (y < BCYMIN) + { + y1 += BCYMAX - BCYMIN; + x1 = -x1; + particle->angle = PI - particle->angle; + particle->vx *= -1.0; + *px *= -1.0; + move++; + } + } + if (((x >= BCXMAX)||(x <= BCXMIN))&&((y >= BCYMAX)||(y <= BCYMIN))) + { + /* This case can lead to numerical instabilities, and can be avoided by putting obstacles in the corners */ + printf("Double wrap!\n"); + if (x >= BCXMAX) x1 = BCXMAX - BCXMIN - x1; + else x1 = -BCXMAX + BCXMIN - x1; + if (y >= BCXMAX) y1 = BCYMAX - BCYMIN - y1; + else y1 = -BCYMAX + BCYMIN - y1; + + particle->vx *= -1.0; + *px *= -1.0; + particle->vy *= -1.0; + *py *= -1.0; + + +// if (x1 >= BCXMAX) x1 = BCXMAX - 1.0e-5; +// else if (x1 <= BCXMIN) x1 = BCXMIN + 1.0e-5; +// if (y1 >= BCXMAX) y1 = BCYMAX - 1.0e-5; +// else if (y1 <= BCYMIN) y1 = BCYMIN + 1.0e-5; + move++; + } + + particle->xc = x1; + particle->yc = y1; + return(move); + break; + } + default: + { + /* do nothing */ + break; + } + } + +} + + +void wrap_relative_positions(double x1, double y1, double *x2, double *y2) +/* computes relative positions of particles, taking boundary conditions into account */ +{ +double dxhalf = 0.5*(BCXMAX - BCXMIN), dyhalf = 0.5*(BCYMAX - BCYMIN); +double dx, dy, x3, y3, x4, y4, dh, dv, wrap_x, wrap_y; +int verbose = 0; + + dx = *x2 - x1; + dy = *y2 - y1; + + switch (bc_grouped(BOUNDARY_COND)) { + case (0): /* rectangular b.c. */ + { + /* do nothing */ + break; + } + case (1): /* periodic b.c. */ + { + if (dx > dxhalf) *x2 -= BCXMAX - BCXMIN; + else if (-dx > dxhalf) *x2 += BCXMAX - BCXMIN; + if (dy > dyhalf) *y2 -= BCYMAX - BCYMIN; + else if (-dy > dyhalf) *y2 += BCYMAX - BCYMIN; +// printf("(x2,y2) = (%.3lg, %.3lg)\n", *x2, *y2); + break; + } + case (2): /* Klein bottle b.c. */ + { + if (dx > dxhalf) + { + *x2 -= BCXMAX - BCXMIN; + *y2 *= -1.0; + } + else if (-dx > dxhalf) + { + *x2 += BCXMAX - BCXMIN; + *y2 *= -1.0; + } + dy = *y2 - y1; + if (dy > dyhalf) *y2 -= BCYMAX - BCYMIN; + else if (-dy > dyhalf) *y2 += BCYMAX - BCYMIN; + break; + } + case (3): /* Boy surface - TO FIX */ + { +// wrap_x = 2.0*(BCXMAX - BCXMIN)/(double)HASHX; +// wrap_y = 2.0*(BCYMAX - BCYMIN)/(double)HASHY; + wrap_x = dxhalf; + wrap_y = dyhalf; + + /* find which wrapped point is closest to (x1, y1) */ + dh = 100.0; + if (dx > wrap_x) + { + x3 = *x2 - BCXMAX + BCXMIN; + y3 = -*y2; + dh = (x3-x1)*(x3-x1) + (y3-y1)*(y3-y1); + } + else if (-dx > wrap_x) + { + x3 = *x2 + BCXMAX - BCXMIN; + y3 = -*y2; + dh = (x3-x1)*(x3-x1) + (y3-y1)*(y3-y1); + } + if ((verbose)&&(dh < 100.0)&&(dh > 0.0)) + { + printf("Case 1: (x1,y1) = (%.3lg, %.3lg)\t", x1, y1); + printf("(x2,y2) = (%.3lg, %.3lg) -> (%.3lg, %.3lg)\n", *x2, *y2, x3, y3); + } + + dv = 100.0; + if (dy > wrap_y) + { + x4 = -*x2; + y4 = *y2 - BCYMAX + BCYMIN; + dv = (x4-x1)*(x4-x1) + (y4-y1)*(y4-y1); + } + else if (-dy > wrap_y) + { + x4 = -*x2; + y4 = *y2 + BCYMAX - BCYMIN; + dv = (x4-x1)*(x4-x1) + (y4-y1)*(y4-y1); + } + if ((verbose)&&(dv < 100.0)&&(dv > 0.0)) + { + printf("Case 2: (x1,y1) = (%.3lg, %.3lg)\t", x1, y1); + printf("(x2,y2) = (%.3lg, %.3lg) -> (%.3lg, %.3lg)\n", *x2, *y2, x4, y4); + } + + if (dh < 100.0) + { + *x2 = x3; + *y2 = y3; + } + if (dv < dh) + { + *x2 = x4; + *y2 = y4; + } + break; + } + } +} + + +void compute_relative_positions(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HASHX*HASHY]) +{ + int j, i0, j0, m0, k, m, p, q, n = 0; + double x1, x2, y1, y2, xtemp, ytemp; + + for (j=0; j < ncircles; j++) if (particle[j].active) + { +// i0 = particle[j].hashx; +// j0 = particle[j].hashy; +// m0 = mhash(i0, j0); + m0 = particle[j].hashcell; + x1 = particle[j].xc; + y1 = particle[j].yc; + n = 0; + + for (q=0; q1.0)||(vabs(y2-y1)>1.0))) +// if (((vabs(x2-x1)>0.7)||(vabs(y2-y1)>0.7))) +// { +// printf("(x1, y1) = (%.3lg, %.3lg), (x2, y2) = (%.3lg, %.3lg) -> (%.3lg, %.3lg)\n", x1, y1, xtemp, ytemp, x2, y2); +// printf("particle[%i].delta[%i] = (%.4lg, %.4lg)\n", j, n, particle[j].deltax[n], particle[j].deltay[n]); +// } + n++; + } + else printf("Not enough memory in particle.deltax, particle.deltay\n"); + } + } + particle[j].hash_nneighb = n; + } +} + diff --git a/sub_lj.c b/sub_lj.c index 3eba398..1e74ba3 100644 --- a/sub_lj.c +++ b/sub_lj.c @@ -4,6 +4,8 @@ #include "colors_waves.c" +#define HUE_TYPE0 300.0 /* hue of particles of type 0 */ +#define HUE_TYPE1 90.0 /* hue of particles of type 1 */ int writetiff(char *filename, char *description, int x, int y, int width, int height, int compression) { @@ -191,6 +193,33 @@ void write_text( double x, double y, char *st) } +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; +} + /*********************/ /* drawing routines */ @@ -262,6 +291,26 @@ void draw_rectangle(double x1, double y1, double x2, double y2) glEnd(); } +void draw_colored_rectangle(double x1, double y1, double x2, double y2, double rgb[3]) +{ + glColor3f(rgb[0], rgb[1], rgb[2]); + glBegin(GL_TRIANGLE_FAN); + glVertex2d(x1, y1); + glVertex2d(x2, y1); + glVertex2d(x2, y2); + glVertex2d(x1, y2); + glEnd(); +} + +void draw_triangle(double x1, double y1, double x2, double y2, double x3, double y3) +{ + glBegin(GL_LINE_LOOP); + glVertex2d(x1, y1); + glVertex2d(x2, y2); + glVertex2d(x3, y3); + glEnd(); +} + void draw_colored_triangle(double x1, double y1, double x2, double y2, double x3, double y3, double rgb[3]) { glColor3f(rgb[0], rgb[1], rgb[2]); @@ -386,6 +435,34 @@ void draw_colored_rhombus(double x, double y, double r, double angle, double rgb glEnd(); } +void draw_colored_sector(double xc, double yc, double r1, double r2, double angle1, double angle2, double rgb[3], int nsides) +{ + int i; + double angle, dangle; + + dangle = (angle2 - angle1)/(double)nsides; + + glColor3f(rgb[0], rgb[1], rgb[2]); + glBegin(GL_TRIANGLE_FAN); + glVertex2d(xc + r1*cos(angle1), yc + r1*sin(angle1)); + for (i = 0; i < nsides+1; i++) + { + angle = angle1 + dangle*(double)i; + glVertex2d(xc + r2*cos(angle), yc + r2*sin(angle)); + } + glEnd(); + glBegin(GL_TRIANGLE_FAN); + glVertex2d(xc + r2*cos(angle2), yc + r2*sin(angle2)); + for (i = 0; i < nsides+1; i++) + { + angle = angle1 + dangle*(double)i; + glVertex2d(xc + r1*cos(angle), yc + r1*sin(angle)); + } + glEnd(); +} + + + void init_particle_config(t_particle particles[NMAXCIRCLES]) /* initialise particle configuration */ { @@ -758,4 +835,2388 @@ void init_people_config(t_person people[NMAXCIRCLES]) } +void init_obstacle_config(t_obstacle obstacle[NMAXOBSTACLES]) +/* initialise particle configuration */ +{ + int i, j, n; + double x, y, dx, dy; + + switch (OBSTACLE_PATTERN) { + case (O_CORNERS): + { + n = 0; + for (i = 0; i < 2; i++) + for (j = 0; j < 2; j++) + { + obstacle[n].xc = BCXMIN + ((double)i)*(BCXMAX - BCXMIN); + obstacle[n].yc = BCYMIN + ((double)j)*(BCYMAX - BCYMIN); + obstacle[n].radius = OBSTACLE_RADIUS; + obstacle[n].active = 1; + n++; + } + nobstacles = n; + break; + } + case (O_GALTON_BOARD): + { + dy = (YMAX - YMIN)/((double)NGRIDX + 3); + dx = dy/cos(PI/6.0); + n = 0; + for (i = 0; i < NGRIDX + 1; i++) + for (j = 0; j < i; j++) + { + obstacle[n].yc = YMAX - ((double)i)*dy; + obstacle[n].xc = ((double)j - 0.5*(double)i + 0.5)*dx; + obstacle[n].radius = OBSTACLE_RADIUS; + obstacle[n].active = 1; + n++; + } + nobstacles = n; + break; + } + default: + { + printf("Function init_obstacle_config not defined for this pattern \n"); + } + } +} + + +/* Computation of interaction force */ +double lennard_jones_force(double r, t_particle particle) +{ + int i; + double rmin = 0.01, rplus, ratio = 1.0; + + if (r > REPEL_RADIUS*particle.radius) return(0.0); + else + { + if (r > rmin) rplus = r; + else rplus = rmin; + +// ratio = ipow(particle.eq_dist*particle.radius/rplus, 6); + ratio = ipow(particle.eq_dist*particle.radius/rplus, 3); + ratio = ratio*ratio; + + return((ratio - 2.0*ratio*ratio)/rplus); + } +} + +void aniso_lj_force(double r, double ca, double sa, double ca_rel, double sa_rel, double force[2], t_particle particle) +{ + int i; + double rmin = 0.01, rplus, ratio = 1.0, c2, s2, c4, s4, a, aprime, f1, f2; + + if (r > REPEL_RADIUS*particle.radius) + { + force[0] = 0.0; + force[1] = 0.0; + } + else + { + if (r > rmin) rplus = r; + else rplus = rmin; + +// ratio = ipow(particle.eq_dist*particle.radius/rplus, 6); + ratio = ipow(particle.eq_dist*particle.radius/rplus, 3); + ratio = ratio*ratio; + + /* cos(2phi) and sin(2phi) */ + c2 = ca_rel*ca_rel - sa_rel*sa_rel; + s2 = 2.0*ca_rel*sa_rel; + + /* cos(4phi) and sin(4phi) */ + c4 = c2*c2 - s2*s2; + s4 = 2.0*c2*s2; + + a = 0.5*(9.0 - 7.0*c4); + aprime = 14.0*s4; + + f1 = ratio*(a - ratio)/rplus; + f2 = ratio*aprime/rplus; + + force[0] = f1*ca - f2*sa; + force[1] = f1*sa + f2*ca; + } +} + +void penta_lj_force(double r, double ca, double sa, double ca_rel, double sa_rel, double force[2], t_particle particle) +{ + int i; + double rmin = 0.01, rplus, ratio = 1.0, c2, s2, c4, s4, c5, s5, a, aprime, f1, f2; + static double a0, b0; + static int first = 1; + + if (first) + { + a0 = cos(0.1*PI) + 0.5; + b0 = a0 - 1.0; + first = 0; + } + + if (r > REPEL_RADIUS*particle.radius) + { + force[0] = 0.0; + force[1] = 0.0; + } + else + { + if (r > rmin) rplus = r; + else rplus = rmin; + +// ratio = ipow(particle.eq_dist*particle.radius/rplus, 6); + ratio = ipow(particle.eq_dist*particle.radius/rplus, 3); + ratio = ratio*ratio; + + + /* cos(2phi) and sin(2phi) */ + c2 = ca_rel*ca_rel - sa_rel*sa_rel; + s2 = 2.0*ca_rel*sa_rel; + + /* cos(4phi) and sin(4phi) */ + c4 = c2*c2 - s2*s2; + s4 = 2.0*c2*s2; + + /* cos(5phi) and sin(5phi) */ + c5 = ca_rel*c4 - sa_rel*s4; + s5 = sa_rel*c4 + ca_rel*s4; + + a = a0 - b0*c5; + aprime = 5.0*b0*s5; + + f1 = ratio*(a - ratio)/rplus; + f2 = ratio*aprime/rplus; + + force[0] = f1*ca - f2*sa; + force[1] = f1*sa + f2*ca; + } +} + +double old_golden_ratio_force(double r, t_particle particle) +/* potential with two minima at distances whose ratio is the golden ratio Phi */ +/* old version that does not work very well */ +{ + int i; + double x, y, z, rplus, ratio = 1.0, phi, a, phi3; + static int first = 1; + static double rmin, b, c, d; + + if (first) + { + rmin = 0.5*particle.radius; + phi = 0.5*(1.0 + sqrt(5.0)); + phi3 = 1.0/(phi*phi*phi); + a = 0.66; + b = 1.0 + phi3 + a; + d = phi3*a; + c = phi3 + a + d; +// b = 7.04; +// c = 13.66; +// d = 6.7; + first = 0; + printf("a = %.4lg, b = %.4lg, c = %.4lg, d = %.4lg\n", a, b, c, d); + } + + if (r > REPEL_RADIUS*particle.radius) return(0.0); + else + { + if (r > rmin) rplus = r; + else rplus = rmin; + + x = particle.eq_dist*particle.radius/rplus; + y = x*x*x; + z = d - c*y + b*y*y - y*y*y; + return(x*z/rplus); + } +} + +double golden_ratio_force(double r, t_particle particle) +/* potential with two minima at distances whose ratio is the golden ratio Phi */ +/* piecewise polynomial/LJ version */ +{ + int i; + double x, rplus, xm6, y1; + static int first = 1; + static double rmin, phi, a, h1, h2, phi6; + + if (first) + { + rmin = 0.5*particle.radius; + phi = 0.5*(1.0 + sqrt(5.0)); + a = 1.2; + + h1 = 1.0; /* inner potential well depth */ + h2 = 10.0; /* outer potential well depth */ + phi6 = ipow(phi, 6); + + first = 0; + } + + if (r > REPEL_RADIUS*particle.radius) return(0.0); + else + { + if (r > rmin) rplus = r; + else rplus = rmin; + + x = rplus/(particle.eq_dist*particle.radius); +// xm6 = 1.0/ipow(x, 6); + xm6 = 1.0/ipow(x, 3); + xm6 = xm6*xm6; + + if (x <= 1.0) return(12.0*h1*xm6*(1.0 - xm6)/x); + else if (x <= a) + { + y1 = ipow(a - 1.0, 3); + return(6.0*h1*(x - 1.0)*(a - x)/y1); + } + else if (x <= phi) + { + y1 = ipow(phi - a, 3); + return(6.0*h2*(x - a)*(x - phi)/y1); + } + else return(12.0*h2*phi6*(1.0 - phi6*xm6)*xm6/x); + } +} + +void dipole_lj_force(double r, double ca, double sa, double ca_rel, double sa_rel, double force[2], t_particle particle) +{ + int i; + double rmin = 0.01, rplus, ratio = 1.0, a, aprime, f1, f2; + + if (r > REPEL_RADIUS*MU) + { + force[0] = 0.0; + force[1] = 0.0; + } + else + { + if (r > rmin) rplus = r; + else rplus = rmin; + +// ratio = ipow(particle.eq_dist*particle.radius/rplus, 6); + ratio = ipow(particle.eq_dist*particle.radius/rplus, 3); + ratio = ratio*ratio; + + a = 1.0 + 0.25*ca_rel; + aprime = -0.25*sa_rel; + + f1 = ratio*(a - ratio)/rplus; + f2 = ratio*aprime/rplus; + + force[0] = f1*ca - f2*sa; + force[1] = f1*sa + f2*ca; + } +} + +void quadrupole_lj_force(double r, double ca, double sa, double ca_rel, double sa_rel, double force[2], t_particle particle) +{ + int i; + double rmin = 0.01, rplus, ratio = 1.0, a, aprime, f1, f2, ca2, sa2, x, y, dplus, dminus; + static int first = 1; + static double a0, b0, aplus, aminus; + + if (first) + { + dplus = cos(0.2*PI)*cos(0.1*PI); +// dminus = 0.8*dplus; + dminus = QUADRUPOLE_RATIO*dplus; + aplus = ipow(1.0/dplus, 6); + aminus = ipow(1.0/dminus, 6); +// aminus = ipow(cos(0.2*PI)*(0.25 + 0.5*sin(0.1*PI)), 6); + a0 = 0.5*(aplus + aminus); + b0 = 0.5*(aplus - aminus); + first = 0; + } + + if (r > REPEL_RADIUS*particle.radius) + { + force[0] = 0.0; + force[1] = 0.0; + } + else + { + if (r > rmin) rplus = r; + else rplus = rmin; + +// ratio = ipow(particle.eq_dist*particle.radius/rplus, 6); + ratio = ipow(particle.eq_dist*particle.radius/rplus, 3); + ratio = ratio*ratio; + + /* cos(2*phi) and sin(2*phi) */ + ca2 = ca_rel*ca_rel - sa_rel*sa_rel; + sa2 = 2.0*ca_rel*sa_rel; + + a = a0 + b0*ca2; +// if (a == 0.0) a = 1.0e-10; + aprime = -2.0*b0*sa2; + + f1 = ratio*(a - ratio)/rplus; + f2 = ratio*aprime/rplus; + + force[0] = f1*ca - f2*sa; + force[1] = f1*sa + f2*ca; + } +} + + +void quadrupole_lj_force2(double r, double ca, double sa, double ca_rel, double sa_rel, double force[2], t_particle particle) +{ + int i; + double rmin = 0.01, rplus, ratio = 1.0, a, aprime, f1, f2, ca2, sa2, x, y, eqdist; + static int first = 1; + static double aplus, aminus, a0, b0; + + if (first) + { + aplus = ipow(cos(0.2*PI)*cos(0.1*PI), 6); + aminus = 0.1*aplus; +// aminus = 0.0; +// aminus = -2.0*ipow(cos(0.2*PI)*(0.5*sin(0.1*PI)), 6); +// aminus = ipow(cos(0.2*PI)*(0.25 + 0.5*sin(0.1*PI)), 6); + a0 = 0.5*(aplus + aminus); + b0 = 0.5*(aplus - aminus); + first = 0; + } + + if (r > REPEL_RADIUS*particle.radius) + { + force[0] = 0.0; + force[1] = 0.0; + } + else + { + if (r > rmin) rplus = r; + else rplus = rmin; + + /* correct distance */ +// ratio = ipow(particle.eq_dist*particle.radius/rplus, 6); + ratio = ipow(particle.eq_dist*particle.radius/rplus, 3); + ratio = ratio*ratio; + + /* cos(2*phi) and sin(2*phi) */ + + ca2 = ca_rel*ca_rel - sa_rel*sa_rel; + sa2 = 2.0*ca_rel*sa_rel; + + a = a0 + b0*ca2; + if (a == 0.0) a = 1.0e-10; + aprime = -2.0*b0*sa2; + +// f1 = ratio*(a - ratio)/rplus; +// f2 = ratio*aprime/rplus; + + f1 = ratio*(aplus - ratio)/(rplus); + f2 = ratio*(aminus - ratio)/(rplus); + +// force[0] = f1*ca_rel - f2*sa_rel; +// force[1] = f1*sa_rel + f2*ca_rel; + + force[0] = f1*ca - f2*sa; + force[1] = f1*sa + f2*ca; + } +} + +double water_torque(double r, double ca, double sa, double ca_rel, double sa_rel, double ck_rel, double sk_rel) +/* compute torque of water molecule #k on water molecule #j (for interaction I_LJ_WATER) - OLD VERSION */ +{ + double c1p, c1m, c2p, c2m, s2p, s2m, s21, s21p, s21m, c21, c21p, c21m, torque; + double r2, rd, rd2, rr[3][3]; + static double cw = -0.5, sw = 0.866025404, delta = 1.5*MU, d2 = 2.25*MU*MU; + int i, j; + + c1p = ck_rel*cw - sk_rel*sw; + c1m = ck_rel*cw + sk_rel*sw; + c2p = ca_rel*cw - sa_rel*sw; + c2m = ca_rel*cw + sa_rel*sw; + s2p = sa_rel*cw + ca_rel*sw; + s2m = sa_rel*cw - ca_rel*sw; + + s21 = sa_rel*ck_rel - ca_rel*sk_rel; + c21 = ca_rel*ck_rel + sa_rel*sk_rel; + + s21p = s21*cw - c21*sw; + s21m = s21*cw + c21*sw; + c21p = c21*cw + s21*sw; + c21m = c21*cw - s21*sw; + + r2 = r*r; + rd = 2.0*r*delta; + rd2 = r2 + d2; + + rr[0][0] = r; + rr[0][1] = sqrt(rd2 + rd*c2p); + rr[0][2] = sqrt(rd2 + rd*c2m); + rr[1][0] = sqrt(rd2 - rd*c1p); + rr[2][0] = sqrt(rd2 - rd*c1m); + + rr[1][1] = sqrt(r2 + rd*(c2p - c1p) + 2.0*d2*(1.0 - c21)); + rr[1][2] = sqrt(r2 + rd*(c2m - c1p) + 2.0*d2*(1.0 - c21m)); + rr[2][1] = sqrt(r2 + rd*(c2p - c1m) + 2.0*d2*(1.0 - c21p)); + rr[2][2] = sqrt(r2 + rd*(c2m - c1m) + 2.0*d2*(1.0 - c21)); + + for (i=0; i<3; i++) for (j=0; j<3; j++) + { + if (rr[i][j] < 1.0e-4) rr[i][j] = 1.0e-4; + rr[i][j] = rr[i][j]*rr[i][j]*rr[i][j]; + } + + torque = rd*(s2p/rr[0][1] + s2m/rr[0][2]); + torque += -0.5*rd*(s2p/rr[1][1] + s2p/rr[2][1] + s2m/rr[1][2] + s2m/rr[2][2]); + torque += d2*(s21/rr[1][1] + s21/rr[2][2] + s21m/rr[1][2] + s21p/rr[2][1]); + + return(torque); + +} + +double water_force(double r, double ca, double sa, double ca_rel, double sa_rel, double ck_rel, double sk_rel, double f[2]) +/* compute force and torque of water molecule #k on water molecule #j (for interaction I_LJ_WATER) */ +{ + double x1[3], y1[3], x2[3], y2[3], rr[3][3], dx[3][3], dy[3][3], fx[3][3], fy[3][3], m[3][3], torque = 0.0; + static double cw[3], sw[3], q[3], d[3], delta = 1.25*MU, dmin = 0.5*MU, fscale = 1.0; + int i, j; + static int first = 1; + + if (first) + { + cw[0] = 1.0; cw[1] = -0.5; cw[2] = -0.5; + sw[0] = 0.0; sw[1] = 866025404; sw[2] = -866025404; /* sines and cosines of angles */ + q[0] = -2.0; q[1] = 1.0; q[2] = 1.0; /* charges */ + d[0] = 0.5*delta; d[1] = delta; d[2] = delta; /* distances to center */ + first = 0; + } + + /* positions of O and H atoms */ + for (i=0; i<3; i++) + { + x1[i] = d[i]*(ca_rel*cw[i] - sa_rel*sw[i]); + y1[i] = d[i]*(ca_rel*sw[i] + sa_rel*cw[i]); + x2[i] = r + d[i]*(ck_rel*cw[i] - sk_rel*sw[i]); + y2[i] = d[i]*(ck_rel*sw[i] + sk_rel*cw[i]); + } + + /* relative positions */ + for (i=0; i<3; i++) for (j=0; j<3; j++) + { + dx[i][j] = x2[j] - x1[i]; + dy[i][j] = y2[j] - y1[i]; + rr[i][j] = module2(dx[i][j], dy[i][j]); + if (rr[i][j] < dmin) rr[i][j] = dmin; + rr[i][j] = ipow(rr[i][j],3); +// rr[i][j] = rr[i][j]*rr[i][j]*rr[i][j]; + } + + /* forces between particles */ + for (i=0; i<3; i++) for (j=0; j<3; j++) + { + fx[i][j] = -q[i]*q[j]*dx[i][j]/rr[i][j]; + fy[i][j] = -q[i]*q[j]*dy[i][j]/rr[i][j]; + } + + /* torques between particles */ + for (i=0; i<3; i++) for (j=0; j<3; j++) + { + m[i][j] = x1[i]*fy[i][j] - y1[i]*fx[i][j]; + } + + /* total force */ + f[0] = 0.0; + f[1] = 0.0; + for (i=0; i<3; i++) for (j=0; j<3; j++) + { + f[0] += fscale*fx[i][j]; + f[1] += fscale*fy[i][j]; + torque += fscale*m[i][j]; + } + + 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) +/* compute repelling force and torque of particle #k on particle #i */ +/* returns 1 if distance between particles is smaller than NBH_DIST_FACTOR*MU */ +{ + double x1, y1, x2, y2, r, f, angle, aniso, fx, fy, ff[2], dist_scaled, spin_f, ck, sk, ck_rel, sk_rel; + static double dxhalf = 0.5*(BCXMAX - BCXMIN), dyhalf = 0.5*(BCYMAX - BCYMIN); + int wwrapx, wwrapy; + + x1 = particle[i].xc; + y1 = particle[i].yc; + x2 = particle[k].xc; + y2 = particle[k].yc; + + wwrapx = ((BOUNDARY_COND == BC_KLEIN)||(BOUNDARY_COND == BC_BOY))&&(vabs(x2 - x1) > dxhalf); + wwrapy = (BOUNDARY_COND == BC_BOY)&&(vabs(y2 - y1) > dyhalf); + + switch (particle[k].interaction) { + case (I_COULOMB): + { + f = -krepel/(1.0e-8 + distance*distance); + force[0] = f*ca; + force[1] = f*sa; + break; + } + case (I_LENNARD_JONES): + { + f = krepel*lennard_jones_force(distance, particle[k]); + force[0] = f*ca; + force[1] = f*sa; + break; + } + case (I_LJ_DIRECTIONAL): + { + aniso_lj_force(distance, ca, sa, ca_rel, sa_rel, ff, particle[k]); + force[0] = krepel*ff[0]; + force[1] = krepel*ff[1]; + break; + } + case (I_LJ_PENTA): + { + penta_lj_force(distance, ca, sa, ca_rel, sa_rel, ff, particle[k]); + force[0] = krepel*ff[0]; + force[1] = krepel*ff[1]; + break; + } + case (I_GOLDENRATIO): + { + f = krepel*golden_ratio_force(distance, particle[k]); + force[0] = f*ca; + force[1] = f*sa; + break; + } + case (I_LJ_DIPOLE): + { + dipole_lj_force(distance, ca, sa, ca_rel, sa_rel, ff, particle[k]); + force[0] = krepel*ff[0]; + force[1] = krepel*ff[1]; + break; + } + case (I_LJ_QUADRUPOLE): + { + quadrupole_lj_force(distance, ca, sa, ca_rel, sa_rel, ff, particle[k]); + force[0] = krepel*ff[0]; + force[1] = krepel*ff[1]; + break; + } + case (I_LJ_WATER): + { + f = krepel*lennard_jones_force(distance, particle[k]); + force[0] = f*ca; + force[1] = f*sa; + break; + } + } + + if (ROTATION) + { + dist_scaled = distance/(particle[i].spin_range*particle[i].radius); + switch (particle[k].interaction) { + case (I_LJ_WATER): + { + ck = cos(particle[k].angle); + sk = sin(particle[k].angle); + ck_rel = ca*ck + sa*sk; + sk_rel = sa*ck - ca*sk; +// *torque = (-3.0*ca_rel*sk_rel + 2.0*sa_rel*ck_rel)/(1.0e-12 + dist_scaled*dist_scaled*dist_scaled); +// *torque = water_torque(distance, ca, sa, ca_rel, sa_rel, ck_rel, sk_rel); +// *torque = (0.5*sin(angle) + 0.5*sin(2.0*angle) - 0.45*sin(3.0*angle))/(1.0e-12 + dist_scaled*dist_scaled*dist_scaled); + + *torque = water_force(distance, ca, sa, ca_rel, sa_rel, ck_rel, sk_rel, ff); + force[0] += ff[0]; + force[1] += ff[1]; +// printf("force = (%.3lg, %.3lg)\n", ff[0], ff[1]); + break; + } + default: + { + spin_f = particle[i].spin_freq; + if (wwrapx||wwrapy) *torque = sin(spin_f*(-particle[k].angle - particle[i].angle))/(1.0e-8 + dist_scaled*dist_scaled); + else + *torque = sin(spin_f*(particle[k].angle - particle[i].angle))/(1.0e-8 + dist_scaled*dist_scaled); + } + } + + if (particle[i].type == particle[k].type) + { + if (particle[i].type == 0) *torque *= KTORQUE; + else *torque *= KTORQUE_B; + } + else *torque *= KTORQUE_DIFF; + } + else *torque = 0.0; + + if ((distance < NBH_DIST_FACTOR*particle[i].radius)&&(k != i)) return(1); +// if ((distance < NBH_DIST_FACTOR*particle[i].radius)) return(1); + else return(0); +} + + +int compute_repelling_force(int i, int j, double force[2], double *torque, t_particle* particle, double krepel) +/* compute repelling force of neighbour #j on particle #i */ +/* returns 1 if distance between particles is smaller than NBH_DIST_FACTOR*MU */ +{ + double distance, ca, sa, cj, sj, ca_rel, sa_rel, f[2], ff[2], torque1, ck, sk, ck_rel, sk_rel; + static double distmin = 10.0*((XMAX - XMIN)/HASHX + (YMAX - YMIN)/HASHY); + int interact, k; + + k = particle[i].hashneighbour[j]; + + distance = module2(particle[i].deltax[j], particle[i].deltay[j]); + + /* for monitoring purposes */ + if (distance > distmin) + { + printf("i = %i, hashcell %i, j = %i, hashcell %i\n", i, particle[i].hashcell, k, particle[k].hashcell); + printf("X = (%.3lg, %.3lg)\n", particle[i].xc, particle[i].yc); + printf("Y = (%.3lg, %.3lg) d = %.3lg\n", particle[k].xc, particle[k].yc, distance); + } + + if ((distance == 0.0)||(i == k)) + { + force[0] = 0.0; + force[1] = 0.0; + *torque = 0.0; + return(1); + } + else if (distance > REPEL_RADIUS*particle[i].radius) + { + force[0] = 0.0; + force[1] = 0.0; + *torque = 0.0; + return(0); + } + else + { + /* to avoid numerical problems, assign minimal value to distance */ + if (distance < 0.1*particle[i].radius) distance = 0.1*particle[i].radius; + + ca = (particle[i].deltax[j])/distance; + sa = (particle[i].deltay[j])/distance; + + /* compute relative angle in case particles can rotate */ + if (ROTATION) + { + cj = cos(particle[j].angle); + sj = sin(particle[j].angle); + ca_rel = ca*cj + sa*sj; + sa_rel = sa*cj - ca*sj; + } + else + { + ca_rel = ca; + sa_rel = sa; + } + + interact = compute_particle_interaction(i, k, f, torque, particle, distance, krepel, ca, sa, ca_rel, sa_rel); + + if (SYMMETRIZE_FORCE) + { + torque1 = *torque; +// compute_particle_interaction(k, i, ff, torque, particle, distance, krepel, ca, sa, ca_rel, sa_rel); + ck = cos(particle[j].angle); + sk = sin(particle[j].angle); + ck_rel = ca*ck + sa*sk; + sk_rel = sa*ck - ca*sk; + compute_particle_interaction(k, i, ff, torque, particle, distance, krepel, -ca, -sa, -ck_rel, -sk_rel); + force[0] = 0.5*(f[0] - ff[0]); + force[1] = 0.5*(f[1] - ff[1]); + *torque = 0.5*(torque1 - *torque); +// *torque = 0.5*(*torque + torque1); + } + else + { + force[0] = f[0]; + force[1] = f[1]; + } + +// printf("force = (%.3lg, %.3lg), torque = %.3lg\n", f[0], f[1], *torque); + return(interact); + } +} + + +int resample_particle(int n, int maxtrials, t_particle particle[NMAXCIRCLES]) +/* resample y coordinate of particle n, returns 1 if no collision is created */ +{ + double x, y, dist, dmin = 10.0; + int i, j, closeby = 0, success = 0, trials = 0; + + while ((!success)&&(trials < maxtrials)) + { + success = 1; + x = particle[n].xc - MU*(double)rand()/RAND_MAX; + y = 0.9*(BCYMIN + (BCYMAX - BCYMIN)*(double)rand()/RAND_MAX); + i = 0; + while ((success)&&(i= NMAXCIRCLES)) + { + printf("Cannot add particle at (%.3lg, %.3lg)\n", x, y); + return(0); + } + else + { + i = ncircles; + + particle[i].type = type; + + particle[i].xc = x; + particle[i].yc = y; + particle[i].radius = MU*sqrt(mass); + particle[i].active = 1; + particle[i].neighb = 0; + particle[i].thermostat = 1; + + particle[i].energy = 0.0; + + if (RANDOM_RADIUS) particle[i].radius = particle[i].radius*(0.75 + 0.5*((double)rand()/RAND_MAX)); + + particle[i].mass_inv = 1.0/mass; + if (particle[i].type == 0) particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT; + else particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT; + + particle[i].vx = vx; + particle[i].vy = vy; + particle[i].energy = (particle[i].vx*particle[i].vx + particle[i].vy*particle[i].vy)*particle[i].mass_inv; + + 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; + } + + ncircles++; + + printf("Added particle at (%.3lg, %.3lg)\n", x, y); + printf("Number of particles: %i\n", ncircles); + + return(1); + } +} + +double neighbour_color(int nnbg) +{ + if (nnbg > 7) nnbg = 7; + switch(nnbg){ + case (7): return(340.0); + case (6): return(310.0); + case (5): return(260.0); + case (4): return(200.0); + case (3): return(140.0); + case (2): return(100.0); + case (1): return(70.0); + default: return(30.0); + } +} + + +void compute_entropy(t_particle particle[NMAXCIRCLES], double entropy[2]) +{ + int i, nleft1 = 0, nleft2 = 0; + double p1, p2; + static int first = 1, ntot1 = 0, ntot2 = 0; + static double log2; + + if (first) + { + log2 = log(2.0); + for (i=0; i 0.0) + { + hue = ENERGY_HUE_MIN + (ENERGY_HUE_MAX - ENERGY_HUE_MIN)*ej/PARTICLE_EMAX; + if (hue > ENERGY_HUE_MIN) hue = ENERGY_HUE_MIN; + if (hue < ENERGY_HUE_MAX) hue = ENERGY_HUE_MAX; + } + radius = particle[j].radius; + width = BOUNDARY_WIDTH; + break; + } + case (P_NEIGHBOURS): + { + hue = neighbour_color(particle[j].neighb); + radius = particle[j].radius; + width = BOUNDARY_WIDTH; + break; + } + case (P_BONDS): + { +// if (particle[j].type == 1) hue = 70.0; /* to make second particle type more visible */ +// if (particle[j].type == 1) hue = neighbour_color(7 - particle[j].neighb); +// else + hue = neighbour_color(particle[j].neighb); + radius = particle[j].radius; + width = 1; + break; + } + case (P_ANGLE): + { + angle = particle[j].angle; + hue = angle*particle[j].spin_freq/DPI; + hue -= (double)((int)hue); + huex = (DPI - angle)*particle[j].spin_freq/DPI; + huex -= (double)((int)huex); + angle = PI - angle; + if (angle < 0.0) angle += DPI; + huey = angle*particle[j].spin_freq/DPI; + huey -= (double)((int)huey); + hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*hue; + huex = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*huex; + huey = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*huey; + radius = particle[j].radius; + width = BOUNDARY_WIDTH; + break; + } + case (P_TYPE): + { +// if (particle[j].type == 0) hue = 310.0; +// else hue = 70.0; + if (particle[j].type == 0) hue = HUE_TYPE0; + else hue = HUE_TYPE1; + radius = particle[j].radius; + width = BOUNDARY_WIDTH; + break; + } + case (P_DIRECTION): + { + hue = argument(particle[j].vx, particle[j].vy); + if (hue > DPI) hue -= DPI; + if (hue < 0.0) hue += DPI; + hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*hue/DPI; + radius = particle[j].radius; + width = BOUNDARY_WIDTH; + break; + } + case (P_ANGULAR_SPEED): + { + hue = 160.0*(1.0 + tanh(SLOPE*particle[j].omega)); +// printf("omega = %.3lg, hue = %.3lg\n", particle[j].omega, hue); + radius = particle[j].radius; + width = BOUNDARY_WIDTH; + break; + } + } + + switch (particle[j].interaction) { + case (I_LJ_DIRECTIONAL): + { + nsides = 4; + break; + } + case (I_LJ_PENTA): + { + nsides = 5; + break; + } + case (I_LJ_QUADRUPOLE): + { + nsides = 4; + break; + } + case (I_LJ_WATER): + { + nsides = NSEG; + radius *= 0.75; + break; + } + default: nsides = NSEG; + } + + switch (plot) { + case (P_KINETIC): + { + hsl_to_rgb_turbo(hue, 0.9, 0.5, rgb); + hsl_to_rgb_turbo(hue, 0.9, 0.5, rgbx); + hsl_to_rgb_turbo(hue, 0.9, 0.5, rgby); + break; + } + case (P_BONDS): + { + hsl_to_rgb_turbo(hue, 0.9, 0.5, rgb); + hsl_to_rgb_turbo(hue, 0.9, 0.5, rgbx); + hsl_to_rgb_turbo(hue, 0.9, 0.5, rgby); + break; + } + case (P_DIRECTION): + { + hsl_to_rgb_twilight(hue, 0.9, 0.5, rgb); + hsl_to_rgb_twilight(hue, 0.9, 0.5, rgbx); + hsl_to_rgb_twilight(hue, 0.9, 0.5, rgby); + break; + } + default: + { + hsl_to_rgb(hue, 0.9, 0.5, rgb); + hsl_to_rgb(hue, 0.9, 0.5, rgbx); + hsl_to_rgb(hue, 0.9, 0.5, rgby); + } + } + angle = particle[j].angle + APOLY*DPI; + + draw_one_particle(particle[j], particle[j].xc, particle[j].yc, radius, angle, nsides, width, rgb); + + /* in case of periodic b.c., draw translates of particles */ + if (PERIODIC_BC) + { + x1 = particle[j].xc; + y1 = particle[j].yc; + + for (i=-2; i<3; i++) + for (k=-1; k<2; k++) + draw_one_particle(particle[j], x1 + (double)i*(BCXMAX - BCXMIN), y1 + (double)k*(BCYMAX - BCYMIN), radius, angle, nsides, width, rgb); + } + else if (BOUNDARY_COND == BC_KLEIN) + { + x1 = particle[j].xc; + y1 = particle[j].yc; + + for (i=-2; i<3; i++) + { + if (vabs(i) == 1) sign = -1.0; + else sign = 1.0; + angle1 = angle*sign; + for (k=-1; k<2; k++) + draw_one_particle(particle[j], x1 + (double)i*(BCXMAX - BCXMIN), sign*(y1 + (double)k*(BCYMAX - BCYMIN)), + radius, angle1, nsides, width, rgb); + } + } + else if (BOUNDARY_COND == BC_BOY) + { + x1 = particle[j].xc; + y1 = particle[j].yc; + + for (i=-1; i<2; i++) for (k=-1; k<2; k++) + { + if (vabs(i) == 1) sign = -1.0; + else sign = 1.0; + if (vabs(k) == 1) signy = -1.0; + else signy = 1.0; + if (signy == 1.0) angle1 = angle*sign; + else angle1 = PI - angle; + if (sign == -1.0) draw_one_particle(particle[j], signy*(x1 + (double)i*(BCXMAX - BCXMIN)), + sign*(y1 + (double)k*(BCYMAX - BCYMIN)), radius, angle1, nsides, width, rgbx); + else if (signy == -1.0) draw_one_particle(particle[j], signy*(x1 + (double)i*(BCXMAX - BCXMIN)), + sign*(y1 + (double)k*(BCYMAX - BCYMIN)), radius, angle1, nsides, width, rgby); + else draw_one_particle(particle[j], signy*(x1 + (double)i*(BCXMAX - BCXMIN)), + sign*(y1 + (double)k*(BCYMAX - BCYMIN)), radius, angle1, nsides, width, rgb); + } + } + } + +// /* draw spin vectors */ + if ((DRAW_SPIN)||(DRAW_SPIN_B)) + { + glLineWidth(width); + for (j=0; j maxpressure) maxpressure = meanpress[j]; + printf("Max pressure = %.5lg\n\n", maxpressure); + } + + y = YMAX - 0.1; + if (INCREASE_BETA) /* print temperature */ + { + logratio = log(beta/BETA)/log(0.5*BETA_FACTOR); + if (logratio > 1.0) logratio = 1.0; + if (BETA_FACTOR > 1.0) hue = PARTICLE_HUE_MAX - (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*logratio; + else hue = PARTICLE_HUE_MIN - (PARTICLE_HUE_MIN - PARTICLE_HUE_MAX)*logratio; + erase_area_hsl_turbo(xbox, y + 0.025, 0.37, 0.05, hue, 0.9, 0.5); +// erase_area_hsl_turbo(xmid + 0.1, y + 0.025, 0.4, 0.05, hue, 0.9, 0.5); + if ((hue < 90)||(hue > 270)) glColor3f(1.0, 1.0, 1.0); + else glColor3f(0.0, 0.0, 0.0); + sprintf(message, "Temperature %.2f", 1.0/beta); + write_text(xtext, y, message); +// write_text(xmidtext, y, message); +// y -= 0.1; + +// erase_area_hsl(xxbox, y + 0.025, 0.37, 0.05, 0.0, 0.9, 0.0); +// glColor3f(1.0, 1.0, 1.0); +// sprintf(message, "Pressure %.3f", meanpressure); +// write_text(xxtext, y, message); + } + if (DECREASE_CONTAINER_SIZE) /* print density */ + { + density = (double)ncircles/((lengthcontainer)*(INITYMAX - INITYMIN)); + erase_area_hsl(xbox, y + 0.025, 0.37, 0.05, 0.0, 0.9, 0.0); + glColor3f(1.0, 1.0, 1.0); + sprintf(message, "Density %.3f", density); + write_text(xtext, y, message); + + erase_area_hsl(xmid, y + 0.025, 0.37, 0.05, 0.0, 0.9, 0.0); + glColor3f(1.0, 1.0, 1.0); + sprintf(message, "Temperature %.2f", temperature); + write_text(xmidtext, y, message); + + erase_area_hsl(xxbox, y + 0.025, 0.37, 0.05, 0.0, 0.9, 0.0); + glColor3f(1.0, 1.0, 1.0); + sprintf(message, "Pressure %.3f", meanpressure); + write_text(xxtext, y, message); + + } + else if (INCREASE_KREPEL) /* print force constant */ + { + erase_area_hsl(xbox, y + 0.025, 0.22, 0.05, 0.0, 0.9, 0.0); + glColor3f(1.0, 1.0, 1.0); + sprintf(message, "Force %.0f", krepel); + write_text(xtext + 0.28, y, message); + } + + if (RECORD_PRESSURES) + { + y = FUNNEL_WIDTH + OBSTACLE_RADIUS; + for (i=0; i PARTICLE_HUE_MIN) hue = PARTICLE_HUE_MIN; + if (hue < PARTICLE_HUE_MAX) hue = PARTICLE_HUE_MAX; + hsl_to_rgb_turbo(hue, 0.9, 0.5, rgb); + + dphi = DPI/(double)N_PRESSURES; +// x = 0.95*OBSTACLE_RADIUS*cos(phi); + sphi = sin(phi); + if (sphi < 0.0) draw_colored_sector(0.0, y, 0.95*OBSTACLE_RADIUS, OBSTACLE_RADIUS, phi, phi + dphi, rgb, 10); + else draw_colored_sector(0.0, -y, 0.95*OBSTACLE_RADIUS, OBSTACLE_RADIUS, phi, phi + dphi, rgb, 10); + } + + glColor3f(1.0, 1.0, 1.0); + for (i=-1; i<2; i++) + { + k = N_PRESSURES/4 + i*N_PRESSURES/9; + phi = DPI*(double)k/(double)N_PRESSURES; + pprint = 0.0; + for (j=-2; j<3; j++) pprint += meanpress[k + j]; + sprintf(message, "p = %.0f", pprint*200.0/MAX_PRESSURE); + write_text(0.85*OBSTACLE_RADIUS*cos(phi) - 0.1, -y + 0.85*OBSTACLE_RADIUS*sin(phi), message); + } + } + + if (PARTIAL_THERMO_COUPLING) + { + printf("Temperature %i in average: %.3lg\n", i_temp, temperature); + temp[i_temp] = temperature; + i_temp++; + if (i_temp >= N_T_AVERAGE) i_temp = 0; + + mean_temp = 0.0; + for (i=0; i 270)) glColor3f(1.0, 1.0, 1.0); + else glColor3f(0.0, 0.0, 0.0); + sprintf(message, "Temperature %.2f", mean_temp); + write_text(xtext, y, message); + + } +} + + +void print_ehrenfest_parameters(t_particle particle[NMAXCIRCLES], double pleft, double pright) +{ + char message[100]; + int i, j, nleft1 = 0, nleft2 = 0, nright1 = 0, nright2 = 0; + double density, hue, rgb[3], logratio, y, shiftx = 0.3; + static double xleftbox, xlefttext, xmidbox, xmidtext, xrightbox, xrighttext, pressures[500][2], meanpressure[2]; + static int first = 1, i_pressure, naverage = 500; + + if (first) + { + xleftbox = -1.0; + xlefttext = xleftbox - 0.5; + xrightbox = 1.0; + xrighttext = xrightbox - 0.45; +// xmid = 0.5*(XMIN + XMAX) - 0.1; +// xmidtext = xmid - 0.24; + + meanpressure[0] = 0.0; + meanpressure[1] = 0.0; + for (i=0; i 1.0 - EHRENFEST_RADIUS) + { + if (particle[i].type == 0) nright1++; + else nright2++; + } + } + + y = YMIN + 0.1; + + erase_area_hsl(xleftbox - shiftx, y + 0.025, 0.22, 0.05, 0.0, 0.9, 0.0); + hsl_to_rgb(310.0, 0.9, 0.5, rgb); + glColor3f(rgb[0], rgb[1], rgb[2]); + sprintf(message, "%i particles", nleft1); + write_text(xlefttext + 0.28 - shiftx, y, message); + + erase_area_hsl(xleftbox + shiftx, y + 0.025, 0.22, 0.05, 0.0, 0.9, 0.0); + hsl_to_rgb(70.0, 0.9, 0.5, rgb); + glColor3f(rgb[0], rgb[1], rgb[2]); + sprintf(message, "%i particles", nleft2); + write_text(xlefttext + 0.28 + shiftx, y, message); + + erase_area_hsl(xrightbox - shiftx, y + 0.025, 0.22, 0.05, 0.0, 0.9, 0.0); + hsl_to_rgb(310.0, 0.9, 0.5, rgb); + glColor3f(rgb[0], rgb[1], rgb[2]); + sprintf(message, "%i particles", nright1); + write_text(xrighttext + 0.28 - shiftx, y, message); + + erase_area_hsl(xrightbox + shiftx, y + 0.025, 0.22, 0.05, 0.0, 0.9, 0.0); + hsl_to_rgb(70.0, 0.9, 0.5, rgb); + glColor3f(rgb[0], rgb[1], rgb[2]); + sprintf(message, "%i particles", nright2); + write_text(xrighttext + 0.28 + shiftx, y, message); + y = YMAX - 0.1; + + erase_area_hsl(xleftbox - 0.1, y + 0.025, 0.22, 0.05, 0.0, 0.9, 0.0); + glColor3f(1.0, 1.0, 1.0); + sprintf(message, "Pressure %.2f", 0.001*meanpressure[0]/(double)ncircles); + write_text(xlefttext + 0.25, y, message); + + erase_area_hsl(xrightbox - 0.1, y + 0.025, 0.22, 0.05, 0.0, 0.9, 0.0); + glColor3f(1.0, 1.0, 1.0); + sprintf(message, "Pressure %.2f", 0.001*meanpressure[1]/(double)ncircles); + write_text(xrighttext + 0.2, y, message); + +} + +void print_particle_number(int npart) +{ + char message[100]; + double y = YMAX - 0.1; + static double xleftbox, xlefttext; + static int first = 1; + + if (first) + { + xleftbox = XMIN + 0.5; + xlefttext = xleftbox - 0.5; + first = 0; + } + + erase_area_hsl(xleftbox, y + 0.025, 0.22, 0.05, 0.0, 0.9, 0.0); + glColor3f(1.0, 1.0, 1.0); + if (npart > 1) sprintf(message, "%i particles", npart); + else sprintf(message, "%i particle", npart); + write_text(xlefttext + 0.28, y, message); +} + +void print_entropy(double entropy[2]) +{ + char message[100]; + double y = YMAX - 0.1, rgb[3]; + static double xleftbox, xlefttext, xrightbox, xrighttext; + static int first = 1; + + if (first) + { + xleftbox = XMIN + 0.5; + xlefttext = xleftbox - 0.55; + xrightbox = XMAX - 0.39; + xrighttext = xrightbox - 0.55; + first = 0; + } + + erase_area_hsl(xleftbox, y + 0.025, 0.35, 0.05, 0.0, 0.9, 0.0); + hsl_to_rgb_turbo(HUE_TYPE1, 0.9, 0.5, rgb); + glColor3f(rgb[0], rgb[1], rgb[2]); + sprintf(message, "Entropy = %.4f", entropy[1]); + write_text(xlefttext + 0.28, y, message); + + erase_area_hsl(xrightbox, y + 0.025, 0.35, 0.05, 0.0, 0.9, 0.0); + hsl_to_rgb_turbo(HUE_TYPE0, 0.9, 0.5, rgb); + glColor3f(rgb[0], rgb[1], rgb[2]); + sprintf(message, "Entropy = %.4f", entropy[0]); + write_text(xrighttext + 0.28, y, message); + +} + + +double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NMAXOBSTACLES], + double xleft, double xright, double *pleft, double *pright, double pressure[N_PRESSURES]) +{ + int i, k; + double xmin, xmax, ymin, ymax, padding, r, rp, r2, cphi, sphi, + f, fperp = 0.0, x, y, xtube, distance, dx, dy, width, ybin, angle, x1, x2, h, ytop, norm, dleft, dplus, dminus; + + /* compute force from fixed obstacles */ + if (ADD_FIXED_OBSTACLES) for (i=0; i BCXMAX) particle[j].fx -= KSPRING_BOUNDARY*(particle[j].xc - BCXMAX); + else if (particle[j].xc < BCXMIN) particle[j].fx += KSPRING_BOUNDARY*(BCXMIN - particle[j].xc); + if (particle[j].yc > BCYMAX) particle[j].fy -= KSPRING_BOUNDARY*(particle[j].yc - BCYMAX); + else if (particle[j].yc < BCYMIN) particle[j].fy += KSPRING_BOUNDARY*(BCYMIN - particle[j].yc); + return(fperp); + } + case (BC_RECTANGLE): + { + /* add harmonic force outside rectangular box */ + padding = MU + 0.01; + xmin = xleft + padding; + xmax = xright - padding; + ymin = BCYMIN + padding; + ymax = BCYMAX - padding; + + if (particle[j].xc > xmax) + { + fperp = KSPRING_BOUNDARY*(particle[j].xc - xmax); + particle[j].fx -= fperp; + } + else if (particle[j].xc < xmin) + { + fperp = KSPRING_BOUNDARY*(xmin - particle[j].xc); + particle[j].fx += fperp; + } + if (particle[j].yc > ymax) + { + fperp = KSPRING_BOUNDARY*(particle[j].yc - ymax); + particle[j].fy -= fperp; + } + else if (particle[j].yc < ymin) + { + fperp = KSPRING_BOUNDARY*(ymin - particle[j].yc); + particle[j].fy += fperp; + } +// if (particle[j].xc > xmax) particle[j].fx -= KSPRING_BOUNDARY*(particle[j].xc - xmax); +// else if (particle[j].xc < xmin) particle[j].fx += KSPRING_BOUNDARY*(xmin - particle[j].xc); +// if (particle[j].yc > ymax) particle[j].fy -= KSPRING_BOUNDARY*(particle[j].yc - ymax); +// else if (particle[j].yc < ymin) particle[j].fy += KSPRING_BOUNDARY*(ymin - particle[j].yc); + + return(fperp); + } + case (BC_CIRCLE): + { + /* add harmonic force outside screen */ + if (particle[j].xc > BCXMAX) particle[j].fx -= KSPRING_BOUNDARY*(particle[j].xc - BCXMAX); + else if (particle[j].xc < BCXMIN) particle[j].fx += KSPRING_BOUNDARY*(BCXMIN - particle[j].xc); + if (particle[j].yc > BCYMAX) particle[j].fy -= KSPRING_BOUNDARY*(particle[j].yc - BCYMAX); + else if (particle[j].yc < BCYMIN) particle[j].fy += KSPRING_BOUNDARY*(BCYMIN - particle[j].yc); + + /* add harmonic force from obstacle */ + for (i=-2; i<2; i++) + { + x = xleft + (double)i*(OBSXMAX - OBSXMIN); + if (vabs(particle[j].xc - x) < 1.1*OBSTACLE_RADIUS) + { + r = module2(particle[j].xc - x, particle[j].yc); + if (r < 1.0e-5) r = 1.0e-05; + cphi = (particle[j].xc - x)/r; + sphi = particle[j].yc/r; + padding = MU + 0.03; + if (r < OBSTACLE_RADIUS + padding) + { + f = KSPRING_OBSTACLE*(OBSTACLE_RADIUS + padding - r); + particle[j].fx += f*cphi; + particle[j].fy += f*sphi; + } + } + } + return(fperp); + } + case (BC_PERIODIC_CIRCLE): + { + x = xleft; + if (vabs(particle[j].xc - x) < 1.1*OBSTACLE_RADIUS) + { + r = module2(particle[j].xc - x, particle[j].yc); + if (r < 1.0e-5) r = 1.0e-05; + cphi = (particle[j].xc - x)/r; + sphi = particle[j].yc/r; + padding = MU + 0.03; + if (r < OBSTACLE_RADIUS + padding) + { + f = KSPRING_OBSTACLE*(OBSTACLE_RADIUS + padding - r); + particle[j].fx += f*cphi; + particle[j].fy += f*sphi; + } + } + return(f); + } + case (BC_PERIODIC_TRIANGLE): + { + x = xleft; + x1 = x + OBSTACLE_RADIUS; + x2 = x - OBSTACLE_RADIUS; + h = 2.0*OBSTACLE_RADIUS*tanh(APOLY*PID); + padding = MU + 0.03; +// ytop = 0.5*h*(1.0 - (particle[j].xc - x)/OBSTACLE_RADIUS); + if ((vabs(particle[j].xc - x) < OBSTACLE_RADIUS + padding)&&(vabs(particle[j].yc < h + padding))) + { + /* signed distances to side of triangle */ + dleft = x2 - particle[j].xc; + norm = module2(h, 2.0*OBSTACLE_RADIUS); + + if (particle[j].yc >= 0.0) + { + dplus = (h*particle[j].xc + 2.0*OBSTACLE_RADIUS*particle[j].yc - h*(x+OBSTACLE_RADIUS))/norm; + if ((dleft < padding)&&(dleft > dplus)) /* left side is closer */ + { + f = KSPRING_OBSTACLE*(padding - dleft); + particle[j].fx -= f; + } + else if (dplus < padding) /* top side is closer */ + { + f = KSPRING_OBSTACLE*(padding - dplus); + particle[j].fx += f*h/norm; + particle[j].fy += 2.0*f*OBSTACLE_RADIUS/norm; + } + } + else + { + dminus = (h*particle[j].xc - 2.0*OBSTACLE_RADIUS*particle[j].yc - h*(x+OBSTACLE_RADIUS))/norm; + if ((dleft < padding)&&(dleft > dminus)) /* left side is closer */ + { + f = KSPRING_OBSTACLE*(padding - dleft); + particle[j].fx -= f; + } + else if (dminus < padding) /* bottom side is closer */ + { + f = KSPRING_OBSTACLE*(padding - dminus); + particle[j].fx += f*h/norm; + particle[j].fy += -2.0*f*OBSTACLE_RADIUS/norm; + } + } + /* force from tip of triangle */ + r = module2(particle[j].xc - x1, particle[j].yc); + if (r < padding) + { + if (r < 1.0e-5) r = 1.0e-05; + cphi = (particle[j].xc - x1)/r; + sphi = particle[j].yc/r; + f = KSPRING_OBSTACLE*(padding - r); + particle[j].fx += f*cphi; + particle[j].fy += f*sphi; + } + } + return(f); + } + case (BC_PERIODIC_FUNNEL): + { + x = xleft; + padding = MU + 0.02; + if (vabs(particle[j].yc) > FUNNEL_WIDTH - padding) for (i=-1; i<2; i+=2) + { + r = module2(particle[j].xc - x, particle[j].yc - (double)i*(FUNNEL_WIDTH + OBSTACLE_RADIUS)); + if (r < 1.0e-5) r = 1.0e-05; + cphi = (particle[j].xc - x)/r; + sphi = (particle[j].yc - (double)i*(FUNNEL_WIDTH + OBSTACLE_RADIUS))/r; + if (r < OBSTACLE_RADIUS + padding) + { + f = KSPRING_OBSTACLE*(OBSTACLE_RADIUS + padding - r); + particle[j].fx += f*cphi; + particle[j].fy += f*sphi; + if (RECORD_PRESSURES) + { + angle = argument(cphi, sphi); + if (angle < 0.0) angle += DPI; + k = (int)((double)N_PRESSURES*angle/DPI); + if (k >= N_PRESSURES) k = N_PRESSURES - 1; + pressure[k] += f; + } + } + } + return(f); + } + case (BC_RECTANGLE_LID): + { + r = particle[j].radius; + + if (particle[j].yc < BCYMIN + r) particle[j].fy += KSPRING_BOUNDARY*(BCYMIN + r - particle[j].yc); + else if (particle[j].yc > ylid - r) + { + fperp = KSPRING_BOUNDARY*(particle[j].yc - ylid + r); + particle[j].fy -= fperp; + } + if (particle[j].yc < BCYMAX + r) + { + if (particle[j].xc > BCXMAX - r) particle[j].fx -= KSPRING_BOUNDARY*(particle[j].xc - BCXMAX + r); + else if (particle[j].xc < BCXMIN + r) particle[j].fx += KSPRING_BOUNDARY*(BCXMIN + r - particle[j].xc); + } + return(fperp); + } + case (BC_EHRENFEST): + { + rp = particle[j].radius; + xtube = 1.0 - sqrt(EHRENFEST_RADIUS*EHRENFEST_RADIUS - EHRENFEST_WIDTH*EHRENFEST_WIDTH); + distance = 0.0; + /* middle tube */ + if (vabs(particle[j].xc) <= xtube) + { + if (particle[j].yc > EHRENFEST_WIDTH - rp) + { + distance = particle[j].yc - EHRENFEST_WIDTH; + particle[j].fy -= KSPRING_BOUNDARY*(distance + rp); + } + else if (particle[j].yc < -EHRENFEST_WIDTH + rp) + { + distance = - EHRENFEST_WIDTH - particle[j].yc; + particle[j].fy += KSPRING_BOUNDARY*(distance + rp); + } + } + /* right container */ + else if (particle[j].xc > 0.0) + { + r = module2(particle[j].xc - 1.0, particle[j].yc); + if ((r > EHRENFEST_RADIUS - rp)&&((particle[j].xc > 1.0)||(vabs(particle[j].yc) > EHRENFEST_WIDTH))) + { + cphi = (particle[j].xc - 1.0)/r; + sphi = particle[j].yc/r; + f = KSPRING_BOUNDARY*(EHRENFEST_RADIUS - r - rp); + particle[j].fx += f*cphi; + particle[j].fy += f*sphi; + *pright -= f; + } + } + /* left container */ + else + { + r = module2(particle[j].xc + 1.0, particle[j].yc); + if ((r > EHRENFEST_RADIUS - rp)&&((particle[j].xc < -1.0)||(vabs(particle[j].yc) > EHRENFEST_WIDTH))) + { + cphi = (particle[j].xc + 1.0)/r; + sphi = particle[j].yc/r; + f = KSPRING_BOUNDARY*(EHRENFEST_RADIUS - r - rp); + particle[j].fx += f*cphi; + particle[j].fy += f*sphi; + *pleft -= f; + } + } + + /* add force from "corners" */ + if ((vabs(particle[j].xc) - xtube < rp)&&(vabs(particle[j].yc) - EHRENFEST_WIDTH < rp)) + { + for (i=-1; i<=1; i+=2) + for (k=-1; k<=1; k+=2) + { + distance = module2(particle[j].xc - (double)i*xtube, particle[j].yc - (double)k*EHRENFEST_WIDTH); + if (distance < rp) + { + cphi = (particle[j].xc - (double)i*xtube)/distance; + sphi = (particle[j].yc - (double)k*EHRENFEST_WIDTH)/distance; + f = KSPRING_BOUNDARY*(rp - distance); + particle[j].fx += f*cphi; + particle[j].fy += f*sphi; + } + } + } + return(fperp); + } + case (BC_SCREEN_BINS): + { + /* add harmonic force outside screen */ + if (particle[j].xc > XMAX) particle[j].fx -= KSPRING_BOUNDARY*(particle[j].xc - XMAX); + else if (particle[j].xc < XMIN) particle[j].fx += KSPRING_BOUNDARY*(XMIN - particle[j].xc); + if (particle[j].yc > YMAX + 10.0*MU) particle[j].fy -= KSPRING_BOUNDARY*(particle[j].yc - YMAX - 10.0*MU); + else if (particle[j].yc < YMIN) particle[j].fy += KSPRING_BOUNDARY*(YMIN - particle[j].yc); + + /* force from the bins */ + dy = (YMAX - YMIN)/((double)NGRIDX + 3); + dx = dy/cos(PI/6.0); + rp = particle[j].radius; + width = rp + 0.05*dx; + ybin = 2.75*dy; + + if (particle[j].yc < YMIN + ybin) for (i=-1; i<=NGRIDX; i++) + { + x = ((double)i - 0.5*(double)NGRIDX + 0.5)*dx; + distance = vabs(particle[j].xc - x); + if (distance < width) + { + if (particle[j].xc > x) particle[j].fx += KSPRING_BOUNDARY*(width - distance); + else particle[j].fx -= KSPRING_BOUNDARY*(width - distance); + } + } + else if (particle[j].yc < YMIN + ybin + particle[j].radius) for (i=-1; i<=NGRIDX; i++) + { + x = ((double)i - 0.5*(double)NGRIDX + 0.5)*dx; + distance = module2(particle[j].xc - x, particle[j].yc - YMIN - ybin); + if (distance < rp) + { + if (distance < 1.0e-8) distance = 1.0e-8; + cphi = (particle[j].xc - x)/distance; + sphi = (particle[j].yc - YMIN - ybin)/distance; + f = KSPRING_BOUNDARY*(rp - distance); + particle[j].fx += f*cphi; + particle[j].fy += f*sphi; + } + } + return(fperp); + } + } +} + + +void compute_particle_force(int j, double krepel, t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HASHX*HASHY]) +/* compute force from other particles on particle j */ +{ + int i0, j0, m0, k, m, q, close; + double fx = 0.0, fy = 0.0, force[2], torque = 0.0, torque_ij, x, y; + + particle[j].neighb = 0; + + for (k=0; k TPYE_PROPORTION)) + { + particle[i].type = 1; + particle[i].radius = MU_B; + } + + particle[i].neighb = 0; + particle[i].thermostat = 1; + +// particle[i].energy = 0.0; +// y = particle[i].yc; +// if (y >= YMAX) y -= particle[i].radius; +// if (y <= YMIN) y += particle[i].radius; + + if (RANDOM_RADIUS) particle[i].radius = particle[i].radius*(0.75 + 0.5*((double)rand()/RAND_MAX)); + + if (particle[i].type == 0) + { + particle[i].interaction = INTERACTION; + particle[i].eq_dist = EQUILIBRIUM_DIST; + particle[i].spin_range = SPIN_RANGE; + particle[i].spin_freq = SPIN_INTER_FREQUENCY; + particle[i].mass_inv = 1.0/PARTICLE_MASS; + particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT; + } + else + { + 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; + particle[i].mass_inv = 1.0/PARTICLE_MASS_B; + particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT_B; + } + + particle[i].vx = V_INITIAL*gaussian(); + particle[i].vy = V_INITIAL*gaussian(); + particle[i].energy = (particle[i].vx*particle[i].vx + particle[i].vy*particle[i].vy)*particle[i].mass_inv; + + px[i] = particle[i].vx; + py[i] = particle[i].vy; + + if (ROTATION) + { + particle[i].angle = DPI*(double)rand()/RAND_MAX; + particle[i].omega = OMEGA_INITIAL*gaussian(); + if (COUPLE_ANGLE_TO_THERMOSTAT) + particle[i].energy += particle[i].omega*particle[i].omega*particle[i].inertia_moment_inv; + } + else + { + particle[i].angle = 0.0; + particle[i].omega = 0.0; + } + pangle[i] = particle[i].omega; + } + /* initialize dummy values in case particles are added */ + for (i=ncircles; i < NMAXCIRCLES; i++) + { + particle[i].type = 0; + particle[i].active = 0; + particle[i].neighb = 0; + particle[i].thermostat = 0; + particle[i].energy = 0.0; + particle[i].mass_inv = 1.0/PARTICLE_MASS; + particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT; + particle[i].vx = 0.0; + particle[i].vy = 0.0; + px[i] = 0.0; + py[i] = 0.0; + particle[i].angle = DPI*(double)rand()/RAND_MAX; + particle[i].omega = 0.0; + pangle[i] = 0.0; + particle[i].interaction = INTERACTION; + particle[i].eq_dist = EQUILIBRIUM_DIST; + particle[i].spin_range = SPIN_RANGE; + particle[i].spin_freq = SPIN_INTER_FREQUENCY; + } + + /* add particles at the bottom as seed */ + if (PART_AT_BOTTOM) for (i=0; i<=NPART_BOTTOM; i++) + { + x = XMIN + (double)i*(XMAX - XMIN)/(double)NPART_BOTTOM; + y = YMIN + 2.0*MU; + add_particle(x, y, 0.0, 0.0, MASS_PART_BOTTOM, 0, particle); + } + if (PART_AT_BOTTOM) for (i=0; i<=NPART_BOTTOM; i++) + { + x = XMIN + (double)i*(XMAX - XMIN)/(double)NPART_BOTTOM; + y = YMIN + 4.0*MU; + add_particle(x, y, 0.0, 0.0, MASS_PART_BOTTOM, 0, particle); + } + + /* add larger copies of particles (for Ehrenfest model)*/ + if (EHRENFEST_COPY) + { + for (i=0; i < ncircles; i++) + { + n = ncircles + i; + particle[n].xc = -particle[i].xc; + particle[n].yc = particle[i].yc; + particle[n].vx = -0.5*particle[i].vx; + particle[n].vy = 0.5*particle[i].vy; + px[n] = -0.5*px[i]; + py[n] = 0.5*py[i]; + particle[n].energy = particle[i].energy; + particle[n].radius = 2.0*particle[i].radius; + particle[n].type = 2; + particle[n].mass_inv = 1.25*particle[i].mass_inv; + particle[n].thermostat = 1; + particle[n].interaction = particle[i].interaction; + particle[n].eq_dist = 0.45*particle[i].eq_dist; + + if ((double)rand()/RAND_MAX > 0.6) particle[n].active = 1; + } + ncircles *= 2; + } + + /* change type of tracer particle */ + if (TRACER_PARTICLE) + { + i = 0; + while ((!particle[i].active)||(module2(particle[i].xc, particle[i].yc) > 0.5)) i++; + tracer_n = i; + particle[tracer_n].type = 2; + particle[tracer_n].radius *= 1.5; + particle[tracer_n].mass_inv *= 1.0/TRACER_PARTICLE_MASS; + particle[tracer_n].vx *= 0.1; + particle[tracer_n].vy *= 0.1; + particle[tracer_n].thermostat = 0; + px[tracer_n] *= 0.1; + py[tracer_n] *= 0.1; + } + + /* position-dependent particle type */ + if (POSITION_DEPENDENT_TYPE) for (i=0; i EHRENFEST_RADIUS) + particle[i].active = 0; + } + if (ADD_FIXED_OBSTACLES) + { + for (i=0; i< ncircles; i++) for (j=0; j < nobstacles; j++) + if (module2(particle[i].xc - obstacle[j].xc, particle[i].yc - obstacle[j].yc) < OBSTACLE_RADIUS + particle[i].radius) + particle[i].active = 0; + } + + /* count number of active particles */ + for (i=0; i< ncircles; i++) nactive += particle[i].active; + printf("%i active particles\n", nactive); + + return(nactive); +} + + +int add_particles(t_particle particle[NMAXCIRCLES], double px[NMAXCIRCLES], double py[NMAXCIRCLES], int nadd_particle) +/* add several particles to the system */ +{ +// add_particle(XMIN + 0.1, 0.0, 50.0, 0.0, 3.0, 0, particle); +// px[ncircles - 1] = particle[ncircles - 1].vx; +// py[ncircles - 1] = particle[ncircles - 1].vy; +// particle[ncircles - 1].radius = 1.5*MU; +// j = 0; +// while (module2(particle[j].xc,particle[j].yc) > 0.7) j = rand()%ncircles; +// x = particle[j].xc + 2.5*MU; +// y = particle[j].yc; + +// x = XMIN + (XMAX - XMIN)*rand()/RAND_MAX; +// y = YMAX + 0.01*rand()/RAND_MAX; +// add_particle(x, y, 0.0, 0.0, 1.0, 0, particle); + +// x = XMIN + 0.25*(XMAX - XMIN); +// y = YMAX + 0.01; +// prop = 1.0 - (double)nadd_particle/5.0; +// vx = 100.0*prop; +// add_particle(x, y, vx, -10.0, 5.0*prop, 0, particle); +// particle[ncircles - 1].radius = 10.0*MU*prop; +// particle[ncircles - 1].eq_dist = 2.0; +// particle[ncircles - 1].thermostat = 0; +// px[ncircles - 1] = particle[ncircles - 1].vx; +// py[ncircles - 1] = particle[ncircles - 1].vy; + + printf("Adding a particle\n\n"); + add_particle(MU*(2.0*rand()/RAND_MAX - 1.0), YMAX + 2.0*MU, 0.0, 0.0, PARTICLE_MASS, 0, particle); + + particle[ncircles - 1].radius = MU; + particle[ncircles - 1].eq_dist = EQUILIBRIUM_DIST; + particle[ncircles - 1].thermostat = 0; + px[ncircles - 1] = particle[ncircles - 1].vx; + py[ncircles - 1] = particle[ncircles - 1].vy; + return (nadd_particle + 1); +} + + +void center_momentum(double p[NMAXCIRCLES]) +{ + int i; + double ptot = 0.0, pmean; + + for (i=0; i PMAX) + { + p[i] = PMAX; + floor = 1; + } + else if (p[i] < -PMAX) + { + p[i] = -PMAX; + floor = 1; + } + } + if (floor) printf("Flooring momentum\n"); + return (floor); +} + +int partial_thermostat_coupling(t_particle particle[NMAXCIRCLES], double xmin) +/* only couple particles with x > xmin to thermostat */ +{ + int i, nthermo = 0; + + for (i=0; i xmin) + { + particle[i].thermostat = 1; + nthermo++; + } + else particle[i].thermostat = 0; + } + return(nthermo); +} + +double compute_mean_energy(t_particle particle[NMAXCIRCLES]) +{ + int i, nactive = 0; + double total_energy = 0.0; + + for (i=0; i 0.0) x = -MU; + else x = MU; + polyline[0].x = x; polyline[0].y = -ymax; - xy_to_pos(-MU, -ymax, pos); + xy_to_pos(x, -ymax, pos); polyline[0].posi = pos[0]; polyline[0].posj = pos[1]; @@ -1362,6 +1364,7 @@ int compute_fresnel_coordinates(t_vertex polyline[NMAXPOLY]) // x = sqrt(LAMBDA*LAMBDA - y*y) - LAMBDA*LAMBDA; while (x <= 0.0) x+= MU; + if (LAMBDA < 0.0) x = -x; polyline[i].x = x; polyline[i].y = y; @@ -1371,9 +1374,11 @@ int compute_fresnel_coordinates(t_vertex polyline[NMAXPOLY]) polyline[i].posj = pos[1]; } - polyline[NSEG].x = -MU; + if (LAMBDA > 0.0) x = -MU; + else x = MU; + polyline[NSEG].x = x; polyline[NSEG].y = ymax; - xy_to_pos(-MU, ymax, pos); + xy_to_pos(x, ymax, pos); polyline[NSEG].posi = pos[0]; polyline[NSEG].posj = pos[1]; @@ -1407,7 +1412,7 @@ int compute_double_fresnel_coordinates(t_vertex polyline[NMAXPOLY], double xshif int compute_noisepanel_coordinates(t_vertex polyline[NMAXPOLY]) -/* compute positions of vertices approximating Fresnel lens */ +/* compute positions of vertices of noise panel */ { int i, n, even; double ymax, dy, x, y, x1, pos[2]; @@ -1452,6 +1457,46 @@ int compute_noisepanel_coordinates(t_vertex polyline[NMAXPOLY]) return(2*n); } +int compute_qrd_coordinates(t_vertex polyline[NMAXPOLY]) +/* compute positions of quadratic noise diffuser */ +{ + int n = 0, b, k, k1, kmin, kmax; + double x, y, x1, y1 = YMIN, pos[2]; + + kmin = (int)(XMIN/LAMBDA) - 2; + kmax = (int)(XMAX/LAMBDA) + 2; + + for (b = -1; b <= 1; b+= 2) + { + if (b == 1) y1 = YMAX; + for (k = kmin; k < kmax; k++) + { + x = LAMBDA*((double)(k) - 0.5); + k1 = (k*k) % 13; + if (b == -1) y = YMIN + (MU/13.0)*(14.0 - (double)k1); + else y = YMAX - (MU/13.0)*(14.0 - (double)k1); + + polyline[n].x = x; + polyline[n].y = y1; + xy_to_pos(x, y1, pos); + polyline[n].posi = pos[0]; + polyline[n].posj = pos[1]; + n++; + + polyline[n].x = x; + polyline[n].y = y; + xy_to_pos(x, y, pos); + polyline[n].posi = pos[0]; + polyline[n].posj = pos[1]; + n++; + + y1 = y; + } + } + + return(n); +} + int init_polyline(int depth, t_vertex polyline[NMAXPOLY]) /* initialise variable polyline, for certain polygonal domain shapes */ { @@ -1500,6 +1545,10 @@ int init_polyline(int depth, t_vertex polyline[NMAXPOLY]) { return(compute_noisepanel_coordinates(polyline)); } + case (D_QRD): + { + return(compute_qrd_coordinates(polyline)); + } default: { return(0); @@ -1940,13 +1989,22 @@ int xy_in_billiard(double x, double y) } case (D_FRESNEL): { - if (vabs(y) > 0.9*LAMBDA) return(1); + if (vabs(y) > 0.9*vabs(LAMBDA)) return(1); if (vabs(x) > MU) return(1); - x1 = sqrt(LAMBDA*LAMBDA - y*y) - LAMBDA; + x1 = sqrt(LAMBDA*LAMBDA - y*y) - vabs(LAMBDA); while (x1 <= 0.0) x1 += MU; - if (x < x1) return(0); - else return(1); + if (LAMBDA > 0.0) + { + if (x < x1) return(0); + else return(1); + } + else + { + x1 = -x1; + if (x > x1) return(0); + else return(1); + } } case (D_DOUBLE_FRESNEL): { @@ -1979,11 +2037,35 @@ int xy_in_billiard(double x, double y) if (x1 <= LAMBDA) y1 = 0.1 + MU*x1/LAMBDA; else y1 = 0.1 + 2.0*MU - MU*x1/LAMBDA; return((y > YMIN + y1)&&(y < YMAX - y1)); -// x1 = vabs(x); -// while (x1 > 2.0*LAMBDA) x1 -= 2.0*LAMBDA; -// if (x1 <= LAMBDA) y1 = MU + x1; -// else y1 = 3.0*MU - x1; -// return((y > YMIN + y1)&&(y < YMAX - y1)); + } + case (D_QRD): + { + x1 = vabs(x)/LAMBDA; + k = (int)(x1 + 0.5); + k1 = (k*k) % 13; + y1 = (MU/13.0)*(14.0 - (double)k1); + return ((y > YMIN + y1)&&(y < YMAX - y1)); + } + case (D_CIRCLE_SEGMENT): + { + if (vabs(y) > 0.9*vabs(LAMBDA)) return(1); + + y1 = 0.9*LAMBDA; + x1 = sqrt(LAMBDA*LAMBDA - y1*y1) - vabs(LAMBDA) + MU; + if ((LAMBDA > 0.0)&&(x < x1)) return(1); + else if ((LAMBDA < 0.0)&&(x > -x1)) return(1); + + x1 = sqrt(LAMBDA*LAMBDA - y*y) - vabs(LAMBDA) + MU; + if (LAMBDA > 0.0) + { + if (x < x1) return(0); + else return(1); + } + else + { + if (x > -x1) return(0); + else return(1); + } } case (D_MENGER): { @@ -2950,6 +3032,37 @@ void draw_billiard() /* draws the billiard boundary */ glEnd(); break; } + case (D_QRD): + { + glLineWidth(BOUNDARY_WIDTH); + glBegin(GL_LINE_STRIP); + for (i=0; i 0.0) x = sqrt(LAMBDA*LAMBDA - y*y) - LAMBDA + MU; + else x = -sqrt(LAMBDA*LAMBDA - y*y) - LAMBDA - MU; + xy_to_pos(x, y, pos); + glVertex2d(pos[0], pos[1]); + } + y = 0.9*LAMBDA; + if (LAMBDA > 0.0) x = sqrt(LAMBDA*LAMBDA - y*y) - LAMBDA + MU; + else x = -sqrt(LAMBDA*LAMBDA - y*y) - LAMBDA - MU; + xy_to_pos(x, y, pos); + glVertex2d(pos[0], pos[1]); + glEnd(); + break; + } case (D_CIRCLES): { glLineWidth(BOUNDARY_WIDTH); @@ -3316,4 +3429,110 @@ void draw_color_scheme(double x1, double y1, double x2, double y2, int plot, dou draw_rectangle(x1, y1, x2, y2); } +void draw_color_scheme_palette(double x1, double y1, double x2, double y2, int plot, double min, double max, int palette) +{ + int j, k, ij_botleft[2], ij_topright[2], imin, imax, jmin, jmax; + double y, dy, dy_e, rgb[3], value, lum, amp; + + xy_to_ij(x1, y1, ij_botleft); + xy_to_ij(x2, y2, ij_topright); + + rgb[0] = 0.0; rgb[1] = 0.0; rgb[2] = 0.0; + erase_area_rgb(0.5*(x1 + x2), x2 - x1, 0.5*(y1 + y2), y2 - y1, rgb); + + if (ROTATE_COLOR_SCHEME) + { + jmin = ij_botleft[0]; + jmax = ij_topright[0]; + imin = ij_botleft[1]; + imax = ij_topright[1]; + } + else + { + imin = ij_botleft[0]; + imax = ij_topright[0]; + jmin = ij_botleft[1]; + jmax = ij_topright[1]; + } + + + glBegin(GL_QUADS); + dy = (max - min)/((double)(jmax - jmin)); + dy_e = max/((double)(jmax - jmin)); + + for (j = jmin; j < jmax; j++) + { + switch (plot) { + case (P_AMPLITUDE): + { + value = min + 1.0*dy*(double)(j - jmin); + color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + break; + } + case (P_ENERGY): + { + value = dy_e*(double)(j - jmin)*100.0/E_SCALE; + if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + break; + } + case (P_MEAN_ENERGY): + { + value = dy_e*(double)(j - jmin)*100.0/E_SCALE; + if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + break; + } + case (P_LOG_ENERGY): + { + value = LOG_SHIFT + 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; + } + case (P_LOG_MEAN_ENERGY): + { + value = LOG_SHIFT + 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; + } + case (P_PHASE): + { + value = min + 1.0*dy*(double)(j - jmin); +// lum = (color_amplitude(value, 1.0, 1))*0.5; +// if (lum < 0.0) lum = 0.0; +// hsl_to_rgb(value*360.0, 0.9, 0.5, rgb); +// color_scheme(COLOR_SCHEME, value, 1.0, 1, rgb); +// amp = color_amplitude_linear(value, 1.0, 1); + amp = 0.5*color_amplitude_linear(value, 1.0, 1); + while (amp > 1.0) amp -= 2.0; + while (amp < -1.0) amp += 2.0; + amp_to_rgb(0.5*(1.0 + amp), rgb); + break; + } + } + glColor3f(rgb[0], rgb[1], rgb[2]); + if (ROTATE_COLOR_SCHEME) + { + glVertex2i(j, imin); + glVertex2i(j, imax); + glVertex2i(j+1, imax); + glVertex2i(j+1, imin); + } + else + { + glVertex2i(imin, j); + glVertex2i(imax, j); + glVertex2i(imax, j+1); + glVertex2i(imin, j+1); + } + } + glEnd (); + + glColor3f(1.0, 1.0, 1.0); + glLineWidth(BOUNDARY_WIDTH); + draw_rectangle(x1, y1, x2, y2); +} + diff --git a/sub_wave_comp.c b/sub_wave_comp.c index 4f86f6d..22d46df 100644 --- a/sub_wave_comp.c +++ b/sub_wave_comp.c @@ -527,7 +527,46 @@ int xy_in_billiard_half(double x, double y, int domain, int pattern, int top) } return(1); } - cdefault: + case (D_FRESNEL): + { + if (vabs(y) > 0.9*vabs(LAMBDA)) return(1); + if (vabs(x) > MU) return(1); + + x1 = sqrt(LAMBDA*LAMBDA - y*y) - vabs(LAMBDA); + while (x1 <= 0.0) x1 += MU; + if (LAMBDA > 0.0) + { + if (x < x1) return(0); + else return(1); + } + else + { + if (x > -x1) return(0); + else return(1); + } + } + case (D_CIRCLE_SEGMENT): + { + if (vabs(y) > 0.9*vabs(LAMBDA)) return(1); + + y1 = 0.9*LAMBDA; + x1 = sqrt(LAMBDA*LAMBDA - y1*y1) - vabs(LAMBDA) + MU; + if ((LAMBDA > 0.0)&&(x < x1)) return(1); + else if ((LAMBDA < 0.0)&&(x > -x1)) return(1); + + x1 = sqrt(LAMBDA*LAMBDA - y*y) - vabs(LAMBDA) + MU; + if (LAMBDA > 0.0) + { + if (x < x1) return(0); + else return(1); + } + else + { + if (x > -x1) return(0); + else return(1); + } + } + default: { printf("Function ij_in_billiard not defined for this billiard \n"); return(0); @@ -563,6 +602,7 @@ void draw_billiard_half(int domain, int pattern, int top) /* draws the bill glEnable(GL_SCISSOR_TEST); if (top) glScissor(0.0, YMID, NX, YMID); +// if (top) glScissor(0.0, YMID, NX, YMAX); else glScissor(0.0, 0.0, NX, YMID); if (BLACK) glColor3f(1.0, 1.0, 1.0); @@ -661,6 +701,34 @@ void draw_billiard_half(int domain, int pattern, int top) /* draws the bill } break; } + case (D_FRESNEL): + { + glLineWidth(BOUNDARY_WIDTH); + glBegin(GL_LINE_LOOP); + for (i=0; i 0.0) x = sqrt(LAMBDA*LAMBDA - y*y) - LAMBDA + MU; + else x = -sqrt(LAMBDA*LAMBDA - y*y) - LAMBDA - MU; + xy_to_pos(x, y, pos); + glVertex2d(pos[0], pos[1]); + } + y = 0.9*LAMBDA; + if (LAMBDA > 0.0) x = sqrt(LAMBDA*LAMBDA - y*y) - LAMBDA + MU; + else x = -sqrt(LAMBDA*LAMBDA - y*y) - LAMBDA - MU; + xy_to_pos(x, y, pos); + glVertex2d(pos[0], pos[1]); + glEnd(); + break; + } case (D_MANDELBROT): { /* Do nothing */ @@ -731,7 +799,7 @@ void init_wave_flat_comp( double *phi[NX], double *psi[NX], short int * xy_in[NX } -void draw_wave_comp(double *phi[NX], double *psi[NX], short int *xy_in[NX], double scale, int time) +void draw_wave_comp(double *phi[NX], double *psi[NX], short int *xy_in[NX], double scale, int time, int plot) /* draw the field */ { int i, j, iplus, iminus, jplus, jminus; @@ -747,7 +815,7 @@ void draw_wave_comp(double *phi[NX], double *psi[NX], short int *xy_in[NX], doub { if (((TWOSPEEDS)&&(xy_in[i][j] != 2))||(xy_in[i][j] == 1)) { - switch (PLOT) { + switch (plot) { case (P_AMPLITUDE): { color_scheme(COLOR_SCHEME, phi[i][j], scale, time, rgb); diff --git a/wave_common.c b/wave_common.c index 6970586..4129c19 100644 --- a/wave_common.c +++ b/wave_common.c @@ -154,14 +154,17 @@ void init_wave_flat( double *phi[NX], double *psi[NX], short int * xy_in[NX]) int i, j; double xy[2], dist2; - for (i=0; i= COL_TURBO) color_scheme_asym(COLOR_SCHEME, energy, scale, time, rgb); + else color_scheme(COLOR_SCHEME, energy, scale, time, rgb); + break; + } + case (P_MIXED): + { + if (j > NY/2) color_scheme(COLOR_SCHEME, phi[i][j], scale, time, rgb); + else color_scheme(COLOR_SCHEME, compute_energy(phi, psi, xy_in, i, j), scale, time, rgb); + break; + } + case (P_MEAN_ENERGY): + { + energy = compute_energy(phi, psi, xy_in, i, j); + total_energy[i][j] += energy; + if (COLOR_PALETTE >= COL_TURBO) + color_scheme_asym(COLOR_SCHEME, total_energy[i][j]/(double)(time+1), scale, time, rgb); + else color_scheme(COLOR_SCHEME, total_energy[i][j]/(double)(time+1), scale, time, rgb); + break; + } + case (P_LOG_ENERGY): + { + energy = compute_energy(phi, psi, xy_in, i, j); + color_scheme(COLOR_SCHEME, LOG_SHIFT + LOG_SCALE*log(energy), scale, time, rgb); + break; + } + case (P_LOG_MEAN_ENERGY): + { + energy = compute_energy(phi, psi, xy_in, i, j); + if (energy == 0.0) energy = 1.0e-20; + total_energy[i][j] += energy; + energy = LOG_SHIFT + LOG_SCALE*log(total_energy[i][j]/(double)(time+1)); + color_scheme(COLOR_SCHEME, energy, scale, time, rgb); + break; + } + } + glColor3f(rgb[0], rgb[1], rgb[2]); + + glVertex2i(i, j); + glVertex2i(i+1, j); + glVertex2i(i+1, j+1); + glVertex2i(i, j+1); + } + } + + glEnd (); +} + + +void draw_wave_highres(int size, double *phi[NX], double *psi[NX], double *total_energy[NX], short int *xy_in[NX], double scale, int time, int plot) +/* draw the field, new version with total energy option */ +{ + int i, j, iplus, iminus, jplus, jminus; + double rgb[3], xy[2], x1, y1, x2, y2, velocity, energy, gradientx2, gradienty2; + static double dtinverse = ((double)NX)/(COURANT*(XMAX-XMIN)), dx = (XMAX-XMIN)/((double)NX); + + glBegin(GL_QUADS); + +// printf("dtinverse = %.5lg\n", dtinverse); + + for (i=0; i= COL_TURBO) color_scheme_asym(COLOR_SCHEME, energy, scale, time, rgb); + else color_scheme(COLOR_SCHEME, energy, scale, time, rgb); + break; + } + case (P_MIXED): + { + if (j > NY/2) color_scheme(COLOR_SCHEME, phi[i][j], scale, time, rgb); + else color_scheme(COLOR_SCHEME, compute_energy(phi, psi, xy_in, i, j), scale, time, rgb); + break; + } + case (P_MEAN_ENERGY): + { + energy = compute_energy(phi, psi, xy_in, i, j); + total_energy[i][j] += energy; + if (COLOR_PALETTE >= COL_TURBO) + color_scheme_asym(COLOR_SCHEME, total_energy[i][j]/(double)(time+1), scale, time, rgb); + else color_scheme(COLOR_SCHEME, total_energy[i][j]/(double)(time+1), scale, time, rgb); + break; + } + case (P_LOG_ENERGY): + { + energy = compute_energy(phi, psi, xy_in, i, j); +// energy = LOG_SHIFT + LOG_SCALE*log(energy); +// if (energy < 0.0) energy = 0.0; + color_scheme(COLOR_SCHEME, LOG_SHIFT + LOG_SCALE*log(energy), scale, time, rgb); + break; + } + case (P_LOG_MEAN_ENERGY): + { + energy = compute_energy(phi, psi, xy_in, i, j); + if (energy == 0.0) energy = 1.0e-20; + total_energy[i][j] += energy; + color_scheme(COLOR_SCHEME, LOG_SHIFT + LOG_SCALE*log(total_energy[i][j]/(double)(time+1)), scale, time, rgb); + break; + } + + } + glColor3f(rgb[0], rgb[1], rgb[2]); + + glVertex2i(i, j); + glVertex2i(i+size, j); + glVertex2i(i+size, j+size); + glVertex2i(i, j+size); + } + } + + glEnd (); +} + + + +void draw_wave_ediss(double *phi[NX], double *psi[NX], double *total_energy[NX], double *color_scale[NX], short int *xy_in[NX], + double *gamma[NX], double gammamax, double scale, int time, int plot) +/* draw the field with luminosity depending on damping */ +{ + int i, j, k, iplus, iminus, jplus, jminus; + double rgb[3], xy[2], x1, y1, x2, y2, velocity, field_value, energy, gradientx2, gradienty2, r2; + static double dtinverse = ((double)NX)/(COURANT*(XMAX-XMIN)), dx = (XMAX-XMIN)/((double)NX); + + glBegin(GL_QUADS); + // printf("dtinverse = %.5lg\n", dtinverse); for (i=0; i= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, energy, scale, time, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, energy, scale, time, rgb); + break; + } + case (P_MIXED): + { + if (j > NY/2) color_scheme_palette(COLOR_SCHEME, palette, phi[i][j], scale, time, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, compute_energy(phi, psi, xy_in, i, j), scale, time, rgb); + break; + } + case (P_MEAN_ENERGY): + { + energy = compute_energy(phi, psi, xy_in, i, j); + total_energy[i][j] += energy; + if (COLOR_PALETTE >= COL_TURBO) + color_scheme_asym_palette(COLOR_SCHEME, palette, total_energy[i][j]/(double)(time+1), scale, time, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, total_energy[i][j]/(double)(time+1), scale, time, rgb); + break; + } + case (P_LOG_ENERGY): + { + energy = compute_energy(phi, psi, xy_in, i, j); + color_scheme_palette(COLOR_SCHEME, palette, LOG_SHIFT + LOG_SCALE*log(energy), scale, time, rgb); + break; + } + case (P_LOG_MEAN_ENERGY): + { + energy = compute_energy(phi, psi, xy_in, i, j); + if (energy == 0.0) energy = 1.0e-20; + total_energy[i][j] += energy; + energy = LOG_SHIFT + LOG_SCALE*log(total_energy[i][j]/(double)(time+1)); + color_scheme_palette(COLOR_SCHEME, palette, energy, scale, time, rgb); + break; + } + } + glColor3f(rgb[0], rgb[1], rgb[2]); + + glVertex2i(i, j); + glVertex2i(i+1, j); + glVertex2i(i+1, j+1); + glVertex2i(i, j+1); + } + } + + glEnd (); +} + + +void draw_wave_highres_palette(int size, double *phi[NX], double *psi[NX], double *total_energy[NX], short int *xy_in[NX], double scale, int time, int plot, int palette) +/* same as draw_wave_highres, but with color scheme option */ +{ + int i, j, iplus, iminus, jplus, jminus; + double rgb[3], xy[2], x1, y1, x2, y2, velocity, energy, gradientx2, gradienty2; + static double dtinverse = ((double)NX)/(COURANT*(XMAX-XMIN)), dx = (XMAX-XMIN)/((double)NX); + + glBegin(GL_QUADS); + +// printf("dtinverse = %.5lg\n", dtinverse); + + for (i=0; i= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, energy, scale, time, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, energy, scale, time, rgb); + break; + } + case (P_MIXED): + { + if (j > NY/2) color_scheme_palette(COLOR_SCHEME, palette, phi[i][j], scale, time, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, compute_energy(phi, psi, xy_in, i, j), scale, time, rgb); + break; + } + case (P_MEAN_ENERGY): + { + energy = compute_energy(phi, psi, xy_in, i, j); + total_energy[i][j] += energy; + if (COLOR_PALETTE >= COL_TURBO) + color_scheme_asym_palette(COLOR_SCHEME, palette, total_energy[i][j]/(double)(time+1), scale, time, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, total_energy[i][j]/(double)(time+1), scale, time, rgb); + break; + } + case (P_LOG_ENERGY): + { + energy = compute_energy(phi, psi, xy_in, i, j); +// energy = LOG_SHIFT + LOG_SCALE*log(energy); +// if (energy < 0.0) energy = 0.0; + color_scheme_palette(COLOR_SCHEME, palette, LOG_SHIFT + LOG_SCALE*log(energy), scale, time, rgb); + break; + } + case (P_LOG_MEAN_ENERGY): + { + energy = compute_energy(phi, psi, xy_in, i, j); + if (energy == 0.0) energy = 1.0e-20; + total_energy[i][j] += energy; + color_scheme_palette(COLOR_SCHEME, palette, LOG_SHIFT + LOG_SCALE*log(total_energy[i][j]/(double)(time+1)), scale, time, rgb); + break; + } } glColor3f(rgb[0], rgb[1], rgb[2]); @@ -494,8 +845,3 @@ void draw_wave_highres(int size, double *phi[NX], double *psi[NX], double *total } - - - - - diff --git a/wave_comparison.c b/wave_comparison.c index 46560fc..c877a51 100644 --- a/wave_comparison.c +++ b/wave_comparison.c @@ -43,29 +43,28 @@ #include #define MOVIE 0 /* set to 1 to generate movie */ +#define DOUBLE_MOVIE 0 /* set to 1 to produce movies for wave height and energy simultaneously */ -#define TIME_LAPSE 1 /* set to 1 to add a time-lapse movie at the end */ +#define TIME_LAPSE 0 /* set to 1 to add a time-lapse movie at the end */ #define TIME_LAPSE_FACTOR 4 /* factor of time-lapse movie */ -/* uncomment for higher resolution version */ // #define WINWIDTH 1920 /* window width */ // #define WINHEIGHT 1000 /* window height */ // #define NX 1920 /* number of grid points on x axis */ // #define NY 1000 /* number of grid points on y axis */ // #define YMID 500 /* mid point of display */ -// #define XMIN -2.0 -// #define XMAX 2.0 /* x interval */ +// #define XMIN -0.5 +// #define XMAX 3.5 /* x interval */ // #define YMIN -1.041666667 // #define YMAX 1.041666667 /* y interval for 9/16 aspect ratio */ -/* comment out for higher resolution version */ #define WINWIDTH 1280 /* window width */ #define WINHEIGHT 720 /* window height */ #define NX 1280 /* number of grid points on x axis */ #define NY 720 /* number of grid points on y axis */ #define YMID 360 /* mid point of display */ -#define XMIN -2.0 -#define XMAX 2.0 /* x interval */ +#define XMIN -0.5 +#define XMAX 3.5 /* x interval */ #define YMIN -1.125 #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ @@ -73,8 +72,8 @@ /* Choice of the billiard table */ -#define B_DOMAIN 40 /* choice of domain shape, see list in global_pdes.c */ -#define B_DOMAIN_B 40 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN 43 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN_B 47 /* choice of domain shape, see list in global_pdes.c */ #define CIRCLE_PATTERN 13 /* pattern of circles, see list in global_pdes.c */ #define CIRCLE_PATTERN_B 13 /* pattern of circles, see list in global_pdes.c */ @@ -90,9 +89,9 @@ #define HEX_NONUNIF_COMPRESSSION 0.15 /* compression factor for HEX_NONUNIF pattern */ #define HEX_NONUNIF_COMPRESSSION_B -0.15 /* compression factor for HEX_NONUNIF pattern */ -#define LAMBDA 0.8 /* parameter controlling the dimensions of domain */ -#define MU 0.035 /* parameter controlling the dimensions of domain */ -#define MUB 0.035 /* parameter controlling the dimensions of domain */ +#define LAMBDA -1.1 /* parameter controlling the dimensions of domain */ +#define MU 0.1 /* parameter controlling the dimensions of domain */ +#define MUB 0.1 /* parameter controlling the dimensions of domain */ // #define MU 0.04665361 /* parameter controlling the dimensions of domain */ // #define MUB 0.04665361 /* parameter controlling the dimensions of domain */ #define NPOLY 3 /* number of sides of polygon */ @@ -122,19 +121,19 @@ /* Physical parameters of wave equation */ -#define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */ -#define OSCILLATE_LEFT 0 /* set to 1 to add oscilating boundary condition on the left */ +#define TWOSPEEDS 1 /* set to 1 to replace hardcore boundary by medium with different speed */ +#define OSCILLATE_LEFT 1 /* set to 1 to add oscilating boundary condition on the left */ #define OSCILLATE_TOPBOT 0 /* set to 1 to enforce a planar wave on top and bottom boundary */ -#define OMEGA 0.0 /* frequency of periodic excitation */ -#define AMPLITUDE 0.025 /* amplitude of periodic excitation */ +#define OMEGA 0.004 /* frequency of periodic excitation */ +#define AMPLITUDE 1.0 /* amplitude of periodic excitation */ #define COURANT 0.02 /* Courant number */ -#define COURANTB 0.004 /* Courant number in medium B */ +#define COURANTB 0.01154 /* Courant number in medium B */ // #define COURANTB 0.005 /* Courant number in medium B */ // #define COURANTB 0.008 /* Courant number in medium B */ #define GAMMA 0.0 /* damping factor in wave equation */ // #define GAMMA 1.0e-8 /* damping factor in wave equation */ -#define GAMMAB 1.0e-8 /* damping factor in wave equation */ +#define GAMMAB 0.0 /* damping factor in wave equation */ // #define GAMMAB 1.0e-6 /* damping factor in wave equation */ // #define GAMMAB 2.0e-4 /* damping factor in wave equation */ // #define GAMMAB 2.5e-4 /* damping factor in wave equation */ @@ -150,22 +149,23 @@ /* Boundary conditions, see list in global_pdes.c */ -#define B_COND 3 +#define B_COND 0 /* Parameters for length and speed of simulation */ -#define NSTEPS 3700 /* number of frames of movie */ +#define NSTEPS 2500 /* number of frames of movie */ // #define NSTEPS 3300 /* number of frames of movie */ #define NVID 25 /* number of iterations between images displayed on screen */ #define NSEG 100 /* number of segments of boundary */ -#define INITIAL_TIME 150 /* time after which to start saving frames */ -#define COMPUTE_ENERGIES 1 /* set to 1 to compute and print energies */ +#define INITIAL_TIME 20 /* time after which to start saving frames */ +#define COMPUTE_ENERGIES 0 /* set to 1 to compute and print energies */ #define BOUNDARY_WIDTH 2 /* width of billiard boundary */ #define PAUSE 100 /* number of frames after which to pause */ #define PSLEEP 1 /* sleep time during pause */ #define SLEEP1 1 /* initial sleeping time */ #define SLEEP2 1 /* final sleeping time */ +#define MID_FRAMES 20 /* number of still frames between movies */ #define END_FRAMES 100 /* number of still frames at end of movie */ /* Parameters of initial condition */ @@ -178,11 +178,13 @@ /* Plot type, see list in global_pdes.c */ -#define PLOT 4 +#define PLOT 0 + +#define PLOT_B 1 /* Color schemes */ -#define COLOR_PALETTE 12 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE 18 /* Color palette, see list in global_pdes.c */ #define BLACK 1 /* background */ #define BLACK_TEXT 1 /* set to 1 to write text in black instead of white */ @@ -193,7 +195,7 @@ #define SLOPE 1.0 /* sensitivity of color on wave amplitude */ // #define SLOPE 0.75 /* sensitivity of color on wave amplitude */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ -#define E_SCALE 2000.0 /* scaling factor for energy representation */ +#define E_SCALE 200.0 /* scaling factor for energy representation */ #define LOG_SCALE 1.5 /* scaling factor for energy log representation */ #define LOG_SHIFT 1.0 /* shift of colors on log scale */ #define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */ @@ -206,8 +208,8 @@ #define HUEAMP -220.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 30.0 /* scale of color scheme bar */ -#define COLORBAR_RANGE_B 30.0 /* scale of color scheme bar for 2nd part */ +#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 ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ @@ -329,6 +331,12 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX phi_out[0][j] = x - tc[0][j]*(x - phi_in[1][j]) - KAPPA_SIDES*x - GAMMA_SIDES*(x-y); break; } + case (BC_ABS_REFLECT): + { + delta = phi_in[1][j] + phi_in[0][j+1] + phi_in[0][j-1] - 3.0*x; + phi_out[0][j] = x - tc[0][j]*(x - phi_in[1][j]) - KAPPA_SIDES*x - GAMMA_SIDES*(x-y); + break; + } } psi_out[0][j] = x; } @@ -365,6 +373,12 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX phi_out[NX-1][j] = x - tc[NX-1][j]*(x - phi_in[NX-2][j]) - KAPPA_SIDES*x - GAMMA_SIDES*(x-y); break; } + case (BC_ABS_REFLECT): + { + delta = phi_in[NX-2][j] + phi_in[NX-1][j+1] + phi_in[NX-1][j-1] - 3.0*x; + phi_out[NX-1][j] = x - tc[NX-1][j]*(x - phi_in[NX-2][j]) - KAPPA_SIDES*x - GAMMA_SIDES*(x-y); + break; + } } psi_out[NX-1][j] = x; } @@ -414,6 +428,15 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX else phi_out[i][jmid-1] = -y + 2*x + tcc[i][jmid-1]*delta - KAPPA*x - tgamma[i][jmid-1]*(x-y); break; } + case (BC_ABS_REFLECT): + { + iplus = i+1; if (iplus == NX) iplus = NX-1; + iminus = i-1; if (iminus == -1) iminus = 0; + + delta = phi_in[iplus][jmid-1] + phi_in[iminus][jmid-1] + phi_in[i][jmid-2] - 3.0*x; + phi_out[i][jmid-1] = -y + 2*x + tcc[i][jmid-1]*delta - KAPPA*x - tgamma[i][jmid-1]*(x-y); + break; + } } psi_out[i][jmid-1] = x; } @@ -463,6 +486,15 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX else phi_out[i][0] = -y + 2*x + tcc[i][0]*delta - KAPPA*x - tgamma[i][0]*(x-y); break; } + case (BC_ABS_REFLECT): + { + iplus = (i+1); if (iplus == NX) iplus = NX-1; + iminus = (i-1); if (iminus == -1) iminus = 0; + + delta = phi_in[iplus][0] + phi_in[iminus][0] + phi_in[i][1] - 3.0*x; + phi_out[i][0] = x - tc[i][0]*(x - phi_in[i][1]) - KAPPA_TOPBOT*x - GAMMA_TOPBOT*(x-y); + break; + } } psi_out[i][0] = x; } @@ -512,6 +544,15 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX else phi_out[i][NY-1] = -y + 2*x + tcc[i][NY-1]*delta - KAPPA*x - tgamma[i][NY-1]*(x-y); break; } + case (BC_ABS_REFLECT): + { + iplus = (i+1); if (iplus == NX) iplus = NX-1; + iminus = (i-1); if (iminus == -1) iminus = 0; + + delta = phi_in[iplus][NY-1] + phi_in[iminus][NY-1] + phi_in[i][NY-2] - 3.0*x; + phi_out[i][NY-1] = x - tc[i][NY-1]*(x - phi_in[i][NY-2]) - KAPPA_TOPBOT*x - GAMMA_TOPBOT*(x-y); + break; + } } psi_out[i][NY-1] = x; } @@ -529,7 +570,7 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX iplus = i+1; if (iplus == NX) iplus = NX-1; iminus = i-1; if (iminus == -1) iminus = 0; - delta = phi_in[iplus][jmid] + phi_in[iminus][jmid] + phi_in[i][1] - 3.0*x; + delta = phi_in[iplus][jmid] + phi_in[iminus][jmid] + phi_in[i][jmid+1] - 3.0*x; phi_out[i][jmid] = -y + 2*x + tcc[i][jmid]*delta - KAPPA*x - tgamma[i][jmid]*(x-y); break; } @@ -561,6 +602,15 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX else phi_out[i][jmid] = -y + 2*x + tcc[i][jmid]*delta - KAPPA*x - tgamma[i][jmid]*(x-y); break; } + case (BC_ABS_REFLECT): + { + iplus = i+1; if (iplus == NX) iplus = NX-1; + iminus = i-1; if (iminus == -1) iminus = 0; + + delta = phi_in[iplus][jmid] + phi_in[iminus][jmid] + phi_in[i][jmid+1] - 3.0*x; + phi_out[i][jmid] = -y + 2*x + tcc[i][jmid]*delta - KAPPA*x - tgamma[i][jmid]*(x-y); + break; + } } psi_out[i][jmid] = x; } @@ -603,7 +653,7 @@ void evolve_wave(double *phi[NX], double *psi[NX], double *phi_tmp[NX], double * void draw_color_bar(int plot, double range) { if (ROTATE_COLOR_SCHEME) draw_color_scheme(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range); - else draw_color_scheme(1.7, YMIN + 0.1, 1.9, YMAX - 0.1, plot, -range, range); + else draw_color_scheme(XMAX - 0.3, YMIN + 0.1, XMAX - 0.1, YMAX - 0.1, plot, -range, range); } @@ -613,7 +663,7 @@ void animation() double time, scale, energies[6], top_energy, bottom_energy; double *phi[NX], *psi[NX], *phi_tmp[NX], *psi_tmp[NX]; short int *xy_in[NX]; - int i, j, s; + int i, j, s, counter = 0; /* Since NX and NY are big, it seemed wiser to use some memory allocation here */ for (i=0; i