diff --git a/global_3d.c b/global_3d.c index 30c003c..44c265b 100644 --- a/global_3d.c +++ b/global_3d.c @@ -1,39 +1,106 @@ -/* global variables and definition used by sub_wave_3d.c */ +/* global variables and definitions used by sub_wave_3d.c */ /* plot types used by wave_3d */ -#define P_3D_AMPLITUDE 101 /* color depends on amplitude */ -#define P_3D_ANGLE 102 /* color depends on angle with fixed direction */ -#define P_3D_AMP_ANGLE 103 /* color depends on amplitude, luminosity depends on angle */ -#define P_3D_ENERGY 104 /* color depends on energy, luminosity depends on angle */ -#define P_3D_LOG_ENERGY 105 /* color depends on logarithm of energy, luminosity depends on angle */ +#define P_3D_AMPLITUDE 101 /* height/color depends on amplitude - DEPRECATED, instead use set SHADE_3D to 0 */ +#define P_3D_ANGLE 102 /* height/color depends on angle with fixed direction - TO DO */ +#define P_3D_AMP_ANGLE 103 /* height/color depends on amplitude, luminosity depends on angle */ +#define P_3D_ENERGY 104 /* height/color depends on energy, luminosity depends on angle */ +#define P_3D_LOG_ENERGY 105 /* height/color depends on logarithm of energy, luminosity depends on angle */ +#define P_3D_TOTAL_ENERGY 106 /* height/color depends on total energy over time, luminosity depends on angle */ +#define P_3D_LOG_TOTAL_ENERGY 107 /* height/color depends on log on total energy over time, luminosity depends on angle */ +#define P_3D_MEAN_ENERGY 108 /* height/color depends on energy averaged over time, luminosity depends on angle */ +#define P_3D_LOG_MEAN_ENERGY 109 /* height/color depends on log on energy averaged over time, luminosity depends on angle */ #define P_3D_PHASE 111 /* phase of wave */ +/* Choice of simulated reaction-diffusion equation in rde.c */ + +#define E_HEAT 0 /* heat equation */ +#define E_ALLEN_CAHN 1 /* Allen-Cahn equation */ +#define E_CAHN_HILLIARD 2 /* Cahn-Hilliard equation */ +#define E_FHN 3 /* FitzHugh-Nagumo equation */ +#define E_RPS 4 /* rock-paper-scissors equation */ +#define E_RPSLZ 41 /* rock-paper-scissors-lizard-Spock equation */ +#define E_SCHRODINGER 5 /* Schrodinger equation */ + +/* Choice of potential */ + +#define POT_HARMONIC 1 /* harmonic oscillator */ + /* plot types used by rde */ #define Z_AMPLITUDE 0 /* amplitude of first field */ -#define Z_NORM_GRADIENT 22 /* gradient of polar angle */ -#define Z_NORM_GRADIENTX 23 /* direction of gradient of u */ +#define Z_RGB 20 /* RGB plot */ +#define Z_POLAR 21 /* polar angle associated with RBG plot */ +#define Z_NORM_GRADIENT 22 /* gradient of polar angle */ +#define Z_ANGLE_GRADIENT 221 /* direction of polar angle */ +#define Z_NORM_GRADIENTX 23 /* norm of gradient of u */ +#define Z_ANGLE_GRADIENTX 231 /* direction of gradient of u */ #define Z_NORM_GRADIENT_INTENSITY 24 /* gradient and intensity of polar angle */ #define Z_VORTICITY 25 /* curl of polar angle */ +/* for Schrodinger equation */ +#define Z_MODULE 30 /* module squared of first two fields */ +#define Z_ARGUMENT 31 /* argument of first two fields, with luminosity depending on module */ -#define COMPUTE_THETA ((plot == P_POLAR)||(plot == P_GRADIENT)||(plot == P_GRADIENTX)||(plot == P_GRADIENT_INTENSITY)||(plot == P_VORTICITY)) -#define COMPUTE_THETAZ ((zplot == P_POLAR)||(zplot == P_GRADIENT)||(zplot == P_GRADIENTX)||(zplot == P_GRADIENT_INTENSITY)||(zplot == P_VORTICITY)) +/* for RPSLZ equation */ +#define Z_MAXTYPE_RPSLZ 40 /* color of type with maximal density */ +#define Z_THETA_RPSLZ 41 /* polar angle */ +#define Z_NORM_GRADIENT_RPSLZ 42 /* gradient of polar angle */ +/* macros to avoid unnecessary computations in 3D plots */ + +#define COMPUTE_THETA ((cplot == Z_POLAR)||(cplot == Z_NORM_GRADIENT)||(cplot == Z_ANGLE_GRADIENT)||(cplot == Z_NORM_GRADIENT_INTENSITY)||(cplot == Z_VORTICITY)) +#define COMPUTE_THETAZ ((zplot == Z_POLAR)||(zplot == Z_NORM_GRADIENT)||(zplot == Z_ANGLE_GRADIENT)||(zplot == Z_NORM_GRADIENT_INTENSITY)||(zplot == Z_VORTICITY)) + +#define COMPUTE_ENERGY ((zplot == P_3D_ENERGY)||(cplot == P_3D_ENERGY)||(zplot == P_3D_LOG_ENERGY)||(cplot == P_3D_LOG_ENERGY)||(zplot == P_3D_TOTAL_ENERGY)||(cplot == P_3D_TOTAL_ENERGY)||(zplot == P_3D_LOG_TOTAL_ENERGY)||(cplot == P_3D_LOG_TOTAL_ENERGY)||(zplot == P_3D_MEAN_ENERGY)||(cplot == P_3D_MEAN_ENERGY)||(zplot == P_3D_LOG_MEAN_ENERGY)||(cplot == P_3D_LOG_MEAN_ENERGY)) + +#define COMPUTE_LOG_TOTAL_ENERGY ((ZPLOT == P_3D_LOG_TOTAL_ENERGY)||(CPLOT == P_3D_LOG_TOTAL_ENERGY)||(ZPLOT_B == P_3D_LOG_TOTAL_ENERGY)||(CPLOT_B == P_3D_LOG_TOTAL_ENERGY)) + +#define COMPUTE_LOG_MEAN_ENERGY ((ZPLOT == P_3D_LOG_MEAN_ENERGY)||(CPLOT == P_3D_LOG_MEAN_ENERGY)||(ZPLOT_B == P_3D_LOG_MEAN_ENERGY)||(CPLOT_B == P_3D_LOG_MEAN_ENERGY)) + +#define COMPUTE_LOG_ENERGY ((ZPLOT == P_3D_LOG_TOTAL_ENERGY)||(CPLOT == P_3D_LOG_TOTAL_ENERGY)||(ZPLOT_B == P_3D_LOG_TOTAL_ENERGY)||(CPLOT_B == P_3D_LOG_TOTAL_ENERGY)||(ZPLOT == P_3D_LOG_MEAN_ENERGY)||(CPLOT == P_3D_LOG_MEAN_ENERGY)||(ZPLOT_B == P_3D_LOG_MEAN_ENERGY)||(CPLOT_B == P_3D_LOG_MEAN_ENERGY)) + +#define COMPUTE_MEAN_ENERGY ((ZPLOT == P_3D_MEAN_ENERGY)||(CPLOT == P_3D_MEAN_ENERGY)||(ZPLOT_B == P_3D_MEAN_ENERGY)||(CPLOT_B == P_3D_MEAN_ENERGY)||(ZPLOT == P_3D_LOG_MEAN_ENERGY)||(CPLOT == P_3D_LOG_MEAN_ENERGY)||(ZPLOT_B == P_3D_LOG_MEAN_ENERGY)||(CPLOT_B == P_3D_LOG_MEAN_ENERGY)) + +#define COMPUTE_TOTAL_ENERGY ((ZPLOT == P_3D_TOTAL_ENERGY)||(CPLOT == P_3D_TOTAL_ENERGY)||(ZPLOT == P_3D_LOG_TOTAL_ENERGY)||(CPLOT == P_3D_LOG_TOTAL_ENERGY)||(ZPLOT == P_3D_MEAN_ENERGY)||(CPLOT == P_3D_MEAN_ENERGY)||(ZPLOT == P_3D_LOG_MEAN_ENERGY)||(CPLOT == P_3D_LOG_MEAN_ENERGY)||(ZPLOT_B == P_3D_TOTAL_ENERGY)||(CPLOT_B == P_3D_TOTAL_ENERGY)||(ZPLOT_B == P_3D_LOG_TOTAL_ENERGY)||(CPLOT_B == P_3D_LOG_TOTAL_ENERGY)||(ZPLOT_B == P_3D_MEAN_ENERGY)||(CPLOT_B == P_3D_MEAN_ENERGY)||(ZPLOT_B == P_3D_LOG_MEAN_ENERGY)||(CPLOT_B == P_3D_LOG_MEAN_ENERGY)) + + +int global_time = 0; /* structure used for color and height representations */ /* possible extra fields: zfield, cfield, interpolated coordinates */ typedef struct { - double energy; /* wave energy */ - double phase; /* wave phase */ - double log_energy; /* log of wave energy */ - double total_energy; /* total energy since beginning of simulation */ - double cos_angle; /* cos of angle between normal vector and direction of light */ - double rgb[3]; /* RGB color code */ - double *p_zfield; /* pointer to z field */ - double *p_cfield; /* pointer to color field */ + double energy; /* wave energy */ + double phase; /* wave phase */ + double log_energy; /* log of wave energy */ + double total_energy; /* total energy since beginning of simulation */ + double log_total_energy; /* log of total energy since beginning of simulation */ + double mean_energy; /* energy averaged since beginning of simulation */ + double log_mean_energy; /* log of energy averaged since beginning of simulation */ + double cos_angle; /* cos of angle between normal vector and direction of light */ + double rgb[3]; /* RGB color code */ + double *p_zfield[2]; /* pointers to z field (second pointer for option DOUBLE_MOVIE) */ + double *p_cfield[2]; /* pointers to color field (second pointer for option DOUBLE_MOVIE) */ } t_wave; + + +typedef struct +{ + double theta; /* angle for Rock-Paper-Scissors equation */ + double nablax; /* gradient of first field */ + double nablay; /* gradient of first field */ + double field_norm; /* norm of field or gradient */ + double field_arg; /* argument of field or gradient */ + double curl; /* curl of field */ + double cos_angle; /* cos of angle between normal vector and direction of light */ + double rgb[3]; /* RGB color code */ + double *p_zfield[2]; /* pointers to z field (second pointer for option DOUBLE_MOVIE) */ + double *p_cfield[2]; /* pointers to color field (second pointer for option DOUBLE_MOVIE) */ +} t_rde; + + + diff --git a/global_ljones.c b/global_ljones.c index efccf51..fa6ab35 100644 --- a/global_ljones.c +++ b/global_ljones.c @@ -14,6 +14,7 @@ #define NMAXCIRCLES 20000 /* total number of circles/polygons (must be at least NCX*NCY for square grid) */ #define MAXNEIGH 20 /* max number of neighbours kept in memory */ #define NMAXOBSTACLES 100 /* max number of obstacles */ +#define NMAXSEGMENTS 100 /* max number of repelling segments */ #define C_SQUARE 0 /* square grid of circles */ #define C_HEX 1 /* hexagonal/triangular grid of circles */ @@ -38,6 +39,14 @@ #define O_GALTON_BOARD 1 /* Galton board pattern */ #define O_GENUS_TWO 2 /* obstacles in corners of L-shape domeain (for genus 2 b.c.) */ +/* pattern of additional repelling segments */ +#define S_RECTANGLE 0 /* segments forming a rectangle */ +#define S_CUP 1 /* segments forming a cup (for increasing gravity) */ +#define S_HOURGLASS 2 /* segments forming an hour glass */ +#define S_PENTA 3 /* segments forming a pentagon with 3 angles of 120° and 2 right angles */ +#define S_CENTRIFUGE 4 /* segments forming "centrifuge" (polygon with radial segments) */ +#define S_POLY_ELLIPSE 5 /* segments forming a polygonal approximation of an ellipse */ + /* particle interaction */ #define I_COULOMB 0 /* Coulomb force */ @@ -66,6 +75,11 @@ #define BC_BOY 13 /* Boy surface/projective plane (periodic with twisted horizontal and vertical parts) */ #define BC_GENUS_TWO 14 /* surface of genus 2, obtained by identifying opposite sides of an L shape */ +/* Regions for partial thermostat couplings */ + +#define TH_VERTICAL 0 /* only particles at the right of x = PARTIAL_THERMO_SHIFT are coupled */ +#define TH_INSEGMENT 1 /* only particles in region defined by segments are coupled */ + /* Plot types */ #define P_KINETIC 0 /* colors represent kinetic energy of particles */ @@ -76,13 +90,15 @@ #define P_TYPE 5 /* colors represent type of particle */ #define P_DIRECTION 6 /* colors represent direction of velocity */ #define P_ANGULAR_SPEED 7 /* colors represent angular speed */ +#define P_DIRECT_ENERGY 8 /* hues represent direction, luminosity represents energy */ +#define P_DIFF_NEIGHB 9 /* colors represent number of neighbours of different type */ /* Color schemes */ #define C_LUM 0 /* color scheme modifies luminosity (with slow drift of hue) */ #define C_HUE 1 /* color scheme modifies hue */ #define C_PHASE 2 /* color scheme shows phase */ -#define C_ONEDIM 3 /* use preset 1d color scheme (for Turbo, Viridis, Magma, Inferno, Plasma) */ +#define C_ONEDIM 3 /* use preset 1d color scheme (for Turbo, Viridis, Magma, Inferno, Plasma, Twilight) */ #define C_ONEDIM_LINEAR 4 /* use preset 1d color scheme with linear scale */ /* Color palettes */ @@ -121,6 +137,7 @@ typedef struct short int thermostat; /* whether particle is coupled to thermostat */ int hashcell; /* hash cell in which particle is located */ int neighb; /* number of neighbours within given distance */ + int diff_neighb; /* number of neighbours of different type */ int hash_nneighb; /* number of neighbours in hashgrid */ int hashneighbour[9*HASHMAX]; /* particle numbers of neighbours in hashgrid */ double deltax[9*HASHMAX]; /* relative position of neighbours */ @@ -164,11 +181,26 @@ typedef struct short int active; /* circle is active */ } t_obstacle; +typedef struct +{ + double x1, x2, y1, y2; /* extremities of segment */ + double nx, ny; /* normal vector */ + double c; /* constant term in cartesian eq nx*x + ny*y = c */ + double length; /* length of segment */ + short int concave; /* corner is concave, to add extra repelling force */ + double angle1, angle2; /* angles in which concave corners repel */ + short int active; /* segment is active */ + double x01, x02, y01, y02; /* initial values of extremities, in case of rotation/translation */ + double nx0, ny0; /* initial normal vector */ + double angle01, angle02; /* initial values of angles in which concave corners repel */ + double fx, fy; /* x and y-components of force on segment */ +} t_segment; + typedef struct { double xc, yc; /* center of circle */ } t_tracer; -int ncircles, nobstacles, counter = 0; +int ncircles, nobstacles, nsegments, counter = 0; diff --git a/global_pdes.c b/global_pdes.c index 7742fb9..07738bd 100644 --- a/global_pdes.c +++ b/global_pdes.c @@ -11,6 +11,7 @@ /* shape of domain */ +#define D_NOTHING 999 /* no boundaries */ #define D_RECTANGLE 0 /* rectangular domain */ #define D_ELLIPSE 1 /* elliptical domain */ #define D_STADIUM 2 /* stadium-shaped domain */ @@ -94,6 +95,9 @@ #define D_MANDELBROT_CIRCLE 26 /* Mandelbrot set with circular conductor */ #define D_VONKOCH_HEATED 27 /* von Koch snowflake in larger circle */ +/* Variable index of refraction */ +#define IOR_MANDELBROT 1 /* index of refraction depends on escape time in Mandelbrot set */ + /* Boundary conditions */ #define BC_DIRICHLET 0 /* Dirichlet boundary conditions */ @@ -114,7 +118,7 @@ #define P_ENERGY 1 /* plot energy of wave */ #define P_MIXED 2 /* plot amplitude in upper half, energy in lower half */ #define P_MEAN_ENERGY 3 /* energy averaged over time */ -#define P_LOG_ENERGY 4 /* log of energy averaged over time */ +#define P_LOG_ENERGY 4 /* log of energy averaged over time */ #define P_LOG_MEAN_ENERGY 5 /* log of energy averaged over time */ /* For Schrodinger equation */ diff --git a/lennardjones.c b/lennardjones.c index d985c31..97c0510 100644 --- a/lennardjones.c +++ b/lennardjones.c @@ -34,11 +34,12 @@ #include #include /* Sam Leffler's libtiff library. */ #include +#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 */ /* so far incompatible with double movie */ #define TIME_LAPSE_FACTOR 3 /* factor of time-lapse movie */ @@ -53,24 +54,27 @@ #define YMIN -1.125 #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ -#define INITXMIN -1.9 -#define INITXMAX 1.9 /* x interval for initial condition */ -#define INITYMIN -1.0 -#define INITYMAX 1.0 /* y interval for initial condition */ +#define INITXMIN -0.7 +#define INITXMAX 0.7 /* x interval for initial condition */ +#define INITYMIN 0.1 +#define INITYMAX 0.6 /* y interval for initial condition */ -#define BCXMIN -2.0 -#define BCXMAX 2.0 /* x interval for boundary condition */ -#define BCYMIN -1.125 -#define BCYMAX 1.125 /* y interval for boundary condition */ +#define BCXMIN -0.7 +#define BCXMAX 0.7 /* x interval for boundary condition */ +#define BCYMIN -0.85 +#define BCYMAX 0.85 /* y interval for boundary condition */ #define OBSXMIN -2.0 #define OBSXMAX 2.0 /* x interval for motion of obstacle */ #define CIRCLE_PATTERN 8 /* pattern of circles, see list in global_ljones.c */ -#define ADD_FIXED_OBSTACLES 1 /* set to 1 do add fixed circular obstacles */ +#define ADD_FIXED_OBSTACLES 0 /* set to 1 do add fixed circular obstacles */ #define OBSTACLE_PATTERN 2 /* pattern of obstacles, see list in global_ljones.c */ +#define ADD_FIXED_SEGMENTS 1 /* set to 1 to add fixed segments as obstacles */ +#define SEGMENT_PATTERN 2 /* pattern of repelling segments, see list in global_ljones.c */ + #define TWO_TYPES 0 /* set to 1 to have two types of particles */ #define TPYE_PROPORTION 0.8 /* proportion of particles of first type */ #define SYMMETRIZE_FORCE 1 /* set to 1 to symmetrize two-particle interaction, only needed if particles are not all the same */ @@ -85,23 +89,22 @@ #define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */ #define NPOISSON 100 /* number of points for Poisson C_RAND_POISSON arrangement */ -#define PDISC_DISTANCE 5.75 /* minimal distance in Poisson disc process, controls density of particles */ -// #define PDISC_DISTANCE 2.25 /* minimal distance in Poisson disc process, controls density of particles */ +#define PDISC_DISTANCE 2.75 /* minimal distance in Poisson disc process, controls density of particles */ #define PDISC_CANDIDATES 100 /* number of candidates in construction of Poisson disc process */ #define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */ #define LAMBDA 2.0 /* parameter controlling the dimensions of domain */ -#define MU 0.015 /* parameter controlling radius of particles */ +#define MU 0.013 /* parameter controlling radius of particles */ #define MU_B 0.0254 /* parameter controlling radius of particles of second type */ #define NPOLY 3 /* number of sides of polygon */ -#define APOLY 0.125 /* angle by which to turn polygon, in units of Pi/2 */ +#define APOLY 0.666666666 /* angle by which to turn polygon, in units of Pi/2 */ #define MDEPTH 4 /* depth of computation of Menger gasket */ #define MRATIO 3 /* ratio defining Menger gasket */ #define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ #define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */ #define FOCI 1 /* set to 1 to draw focal points of ellipse */ -#define NGRIDX 30 /* number of grid point for grid of disks */ -#define NGRIDY 20 /* number of grid point for grid of disks */ +#define NGRIDX 46 /* number of grid point for grid of disks */ +#define NGRIDY 24 /* number of grid point for grid of disks */ #define EHRENFEST_RADIUS 0.9 /* radius of container for Ehrenfest urn configuration */ #define EHRENFEST_WIDTH 0.035 /* width of tube for Ehrenfest urn configuration */ @@ -112,11 +115,10 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 5100 /* number of frames of movie */ -// #define NSTEPS 2750 /* number of frames of movie */ -#define NVID 200 /* number of iterations between images displayed on screen */ +#define NSTEPS 5500 /* number of frames of movie */ +#define NVID 650 /* number of iterations between images displayed on screen */ #define NSEG 250 /* number of segments of boundary */ -#define INITIAL_TIME 20 /* time after which to start saving frames */ +#define INITIAL_TIME 10 /* time after which to start saving frames */ #define BOUNDARY_WIDTH 1 /* width of particle boundary */ #define LINK_WIDTH 2 /* width of links between particles */ #define CONTAINER_WIDTH 4 /* width of container boundary */ @@ -130,12 +132,12 @@ /* Boundary conditions, see list in global_ljones.c */ -#define BOUNDARY_COND 14 +#define BOUNDARY_COND 0 /* Plot type, see list in global_ljones.c */ -#define PLOT 5 -#define PLOT_B 4 /* plot type for second movie */ +#define PLOT 0 +#define PLOT_B 8 /* plot type for second movie */ #define COLOR_BONDS 1 /* set to 1 to color bonds according to length */ @@ -164,7 +166,7 @@ #define ENERGY_HUE_MAX 50.0 /* color of saturated particle */ #define PARTICLE_HUE_MIN 359.0 /* color of original particle */ #define PARTICLE_HUE_MAX 0.0 /* color of saturated particle */ -#define PARTICLE_EMAX 1.0e3 /* energy of particle with hottest color */ +#define PARTICLE_EMAX 2.0e2 /* energy of particle with hottest color */ #define HUE_TYPE0 280.0 /* hue of particles of type 0 */ #define HUE_TYPE1 135.0 /* hue of particles of type 1 */ #define HUE_TYPE2 70.0 /* hue of particles of type 1 */ @@ -173,7 +175,7 @@ #define RANDOM_RADIUS 0 /* set to 1 for random circle radius */ #define DT_PARTICLE 5.0e-7 /* time step for particle displacement */ #define KREPEL 12.0 /* constant in repelling force between particles */ -#define EQUILIBRIUM_DIST 5.0 /* Lennard-Jones equilibrium distance */ +#define EQUILIBRIUM_DIST 4.0 /* Lennard-Jones equilibrium distance */ #define EQUILIBRIUM_DIST_B 5.0 /* Lennard-Jones equilibrium distance for second type of particle */ #define REPEL_RADIUS 20.0 /* radius in which repelling force acts (in units of particle radius) */ #define DAMPING 0.0 /* damping coefficient of particles */ @@ -185,16 +187,18 @@ #define OMEGA_INITIAL 10.0 /* initial angular velocity range */ #define THERMOSTAT 1 /* set to 1 to switch on thermostat */ +#define VARY_THERMOSTAT 0 /* set to 1 for time-dependent thermostat schedule */ #define SIGMA 5.0 /* noise intensity in thermostat */ -#define BETA 0.0001 /* initial inverse temperature */ +#define BETA 0.01 /* initial inverse temperature */ #define MU_XI 0.01 /* friction constant in thermostat */ #define KSPRING_BOUNDARY 5.0e9 /* confining harmonic potential outside simulation region */ -#define KSPRING_OBSTACLE 5.0e8 /* harmonic potential of obstacles */ -#define NBH_DIST_FACTOR 4.0 /* radius in which to count neighbours */ -#define GRAVITY 0.0 /* gravity acting on all particles */ +#define KSPRING_OBSTACLE 5.0e9 /* harmonic potential of obstacles */ +#define NBH_DIST_FACTOR 4.5 /* radius in which to count neighbours */ +#define GRAVITY 10000.0 /* gravity acting on all particles */ #define INCREASE_GRAVITY 0 /* set to 1 to increase gravity during the simulation */ #define GRAVITY_FACTOR 100.0 /* factor by which to increase gravity */ -#define GRAVITY_RESTORE_TIME 750 /* time at end of simulation with gravity restored to initial value */ +#define GRAVITY_INITIAL_TIME 500 /* time at start of simulation with constant gravity */ +#define GRAVITY_RESTORE_TIME 1000 /* time at end of simulation with gravity restored to initial value */ #define ROTATION 0 /* set to 1 to include rotation of particles */ #define COUPLE_ANGLE_TO_THERMOSTAT 0 /* set to 1 to couple angular degrees of freedom to thermostat */ @@ -212,8 +216,8 @@ #define INCREASE_BETA 0 /* set to 1 to increase BETA during simulation */ #define BETA_FACTOR 20.0 /* factor by which to change BETA during simulation */ #define N_TOSCILLATIONS 1.5 /* number of temperature oscillations in BETA schedule */ -#define NO_OSCILLATION 0 /* set to 1 to have exponential BETA change only */ -#define FINAL_CONSTANT_PHASE 0 /* final phase in which temperature is constant */ +#define NO_OSCILLATION 1 /* set to 1 to have exponential BETA change only */ +#define FINAL_CONSTANT_PHASE 2000 /* final phase in which temperature is constant */ #define DECREASE_CONTAINER_SIZE 0 /* set to 1 to decrease size of container */ #define SYMMETRIC_DECREASE 0 /* set tp 1 to decrease container symmetrically */ @@ -235,6 +239,7 @@ #define N_T_AVERAGE 50 /* size of temperature averaging window */ #define MAX_PRESSURE 3.0e10 /* pressure shown in "hottest" color */ #define PARTIAL_THERMO_COUPLING 0 /* set to 1 to couple only particles to the right of obstacle to thermostat */ +#define PARTIAL_THERMO_REGION 1 /* region for partial thermostat coupling (see list in global_ljones.c) */ #define PARTIAL_THERMO_SHIFT 0.5 /* distance from obstacle at the right of which particles are coupled to thermostat */ #define INCREASE_KREPEL 0 /* set to 1 to increase KREPEL during simulation */ @@ -245,22 +250,40 @@ #define NPART_BOTTOM 100 /* number of particles at the bottom */ #define ADD_PARTICLES 0 /* set to 1 to add particles */ -#define ADD_TIME 25 /* time at which to add first particle */ -#define ADD_PERIOD 20 /* time interval between adding further particles */ -#define FINAL_NOADD_PERIOD 250 /* final period where no particles are added */ +#define ADD_TIME 500 /* time at which to add first particle */ +#define ADD_PERIOD 250 /* time interval between adding further particles */ +#define FINAL_NOADD_PERIOD 200 /* final period where no particles are added */ #define SAFETY_FACTOR 2.0 /* no particles are added at distance less than MU*SAFETY_FACTOR of other particles */ -#define TRACER_PARTICLE 1 /* set to 1 to have a tracer particle */ +#define TRACER_PARTICLE 0 /* set to 1 to have a tracer particle */ #define N_TRACER_PARTICLES 3 /* number of tracer particles */ #define TRAJECTORY_LENGTH 8000 /* length of recorded trajectory */ #define TRACER_PARTICLE_MASS 4.0 /* relative mass of tracer particle */ #define TRAJECTORY_WIDTH 3 /* width of tracer particle trajectory */ +#define ROTATE_BOUNDARY 1 /* set to 1 to rotate the repelling segments */ +#define SMOOTH_ROTATION 1 /* set to 1 to update segments at each time step (rather than at each movie frame) */ +#define PERIOD_ROTATE_BOUNDARY 2500 /* period of rotating boundary */ +#define ROTATE_INITIAL_TIME 0 /* initial time without rotation */ +#define ROTATE_FINAL_TIME 500 /* final time without rotation */ +#define PRINT_OMEGA 0 /* set to 1 to print angular speed */ +#define PRINT_PARTICLE_SPEEDS 0 /* set to 1 to print average speeds/momenta of particles */ + +#define MOVE_BOUNDARY 0 /* set to 1 to move repelling segments, due to force from particles */ +#define SEGMENTS_MASS 100.0 /* mass of collection of segments */ +#define DEACTIVATE_SEGMENT 1 /* set to 1 to deactivate last segment after a certain time */ +#define SEGMENT_DEACTIVATION_TIME 1000 /* time at which to deactivate last segment */ + #define POSITION_DEPENDENT_TYPE 0 /* set to 1 to make particle type depend on initial position */ #define POSITION_Y_DEPENDENCE 0 /* set to 1 for the separation between particles to be vertical */ #define PRINT_ENTROPY 0 /* set to 1 to compute entropy */ +#define REACTION_DIFFUSION 0 /* set to 1 to simulate a chemical reaction (particles may change type) */ +#define RD_TYPES 3 /* number of types in reaction-diffusion equation */ +#define REACTION_PROB 0.03 /* probability controlling reaction term */ + #define PRINT_PARTICLE_NUMBER 0 /* set to 1 to print total number of particles */ +#define PRINT_LEFT 0 /* set to 1 to print certain parameters at the top left instead of right */ #define EHRENFEST_COPY 0 /* set to 1 to add equal number of larger particles (for Ehrenfest model) */ @@ -277,8 +300,8 @@ #define FLOOR_OMEGA 1 /* set to 1 to limit particle momentum to PMAX */ #define PMAX 1000.0 /* maximal force */ -#define HASHX 34 /* size of hashgrid in x direction */ -#define HASHY 18 /* size of hashgrid in y direction */ +#define HASHX 40 /* size of hashgrid in x direction */ +#define HASHY 20 /* size of hashgrid in y direction */ #define HASHMAX 100 /* maximal number of particles per hashgrid cell */ #define HASHGRID_PADDING 0.1 /* padding of hashgrid outside simulation window */ @@ -296,6 +319,14 @@ double ylid = 0.9; /* y coordinate of lid (for BC_RECTANGLE_LID b.c.) */ double vylid = 0.0; /* y speed coordinate of lid (for BC_RECTANGLE_LID b.c.) */ double xwall = 0.0; /* x coordinate of wall (for BC_RECTANGLE_WALL b.c.) */ double vxwall = 0.0; /* x speed of wall (for BC_RECTANGLE_WALL b.c.) */ +double angular_speed = 0.0; /* angular speed of rotating segments */ +double xsegments = 0.0; /* x coordinate of segments (for option MOVE_BOUNDARY) */ +double ysegments = 0.0; /* y coordinate of segments (for option MOVE_BOUNDARY) */ +double vxsegments = 0.0; /* vx coordinate of segments (for option MOVE_BOUNDARY) */ +double vysegments = 0.0; /* vy coordinate of segments (for option MOVE_BOUNDARY) */ +int thermostat_on = 1; /* thermostat switch used when VARY_THERMOSTAT is on */ + +#define THERMOSTAT_ON ((THERMOSTAT)&&((!VARY_THERMOSTAT)||(thermostat_on))) #include "global_ljones.c" #include "sub_lj.c" @@ -332,8 +363,8 @@ double temperature_schedule(int i) if (first) { - bexponent = log(BETA_FACTOR)/(double)(INITIAL_TIME + NSTEPS - FINAL_CONSTANT_PHASE); - omega = N_TOSCILLATIONS*DPI/(double)(INITIAL_TIME + NSTEPS - FINAL_CONSTANT_PHASE); + bexponent = log(BETA_FACTOR)/(double)(NSTEPS - FINAL_CONSTANT_PHASE); + omega = N_TOSCILLATIONS*DPI/(double)(NSTEPS - FINAL_CONSTANT_PHASE); first = 0; } if (i < INITIAL_TIME) beta = BETA; @@ -411,16 +442,73 @@ double gravity_schedule(int i, int j) { double time, gravity; - if ((i < INITIAL_TIME)||(i > NSTEPS + INITIAL_TIME - GRAVITY_RESTORE_TIME)) return(GRAVITY); + if ((i < INITIAL_TIME + GRAVITY_INITIAL_TIME)||(i > NSTEPS + INITIAL_TIME - GRAVITY_RESTORE_TIME)) return(GRAVITY); else { - time = ((double)(i - INITIAL_TIME) + (double)j/(double)NVID)/(double)(NSTEPS - GRAVITY_RESTORE_TIME); + time = ((double)(i - INITIAL_TIME - GRAVITY_INITIAL_TIME) + + (double)j/(double)NVID)/(double)(NSTEPS - GRAVITY_RESTORE_TIME - GRAVITY_INITIAL_TIME); gravity = GRAVITY*(1.0 + time*(GRAVITY_FACTOR - 1.0)); // printf("i = %i, time = %.3lg, Gravity = %.3lg\n", i, time, gravity); return(gravity); } } +double rotation_angle(double phase) +{ + double omegamax = 15.0; + + /* case of rotating hourglass */ + while (phase > DPI) phase -= DPI; + return(phase - 0.5*sin(2.0*phase)); + + /* case of centrifuge */ +// while (phase > DPI) phase -= DPI; +// angular_speed = 0.5*omegamax*(1.0 - cos(phase)); +// return(0.5*omegamax*(phase - sin(phase))); +} + +double rotation_schedule(int i) +{ + double phase; + static int imin = INITIAL_TIME + ROTATE_INITIAL_TIME, imax = INITIAL_TIME + NSTEPS - ROTATE_FINAL_TIME; + + if (i < imin) + { + angular_speed = 0.0; + return(0.0); + } + else + { + if (i > imax) i = imax; + phase = (DPI/(double)PERIOD_ROTATE_BOUNDARY)*(double)(i - imin); + return(rotation_angle(phase)); + } +} + +double rotation_schedule_smooth(int i, int j) +{ + double phase; + static int imin = INITIAL_TIME + ROTATE_INITIAL_TIME, imax = INITIAL_TIME + NSTEPS - ROTATE_FINAL_TIME; + + if (i < imin) + { + angular_speed = 0.0; + return(0.0); + } + else + { + if (i > imax) i = imax; + phase = (DPI/(double)PERIOD_ROTATE_BOUNDARY)*((double)(i - imin) + (double)j/(double)NVID); + return(rotation_angle(phase)); + } +} + +int thermostat_schedule(int i) +{ + if (i < INITIAL_TIME) return(1); + else return(0); +} + double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HASHX*HASHY], double qx[NMAXCIRCLES], double qy[NMAXCIRCLES], double qangle[NMAXCIRCLES], double px[NMAXCIRCLES], double py[NMAXCIRCLES], double pangle[NMAXCIRCLES], @@ -449,7 +537,7 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA qy[j] = particle[j].yc + 0.5*DT_PARTICLE*py[j]*particle[j].mass_inv; qangle[j] = particle[j].angle + 0.5*DT_PARTICLE*pangle[j]*particle[j].inertia_moment_inv; - if ((THERMOSTAT)&&(particle[j].thermostat)) + if ((THERMOSTAT_ON)&&(particle[j].thermostat)) { px[j] *= exp(- 0.5*DT_PARTICLE*xi); py[j] *= exp(- 0.5*DT_PARTICLE*xi); @@ -467,7 +555,7 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA // *nactive++; } totalenergy *= DIMENSION_FACTOR; /* normalize energy to take number of degrees of freedom into account */ - if (THERMOSTAT) + if (THERMOSTAT_ON) { a = DT_PARTICLE*(totalenergy - (double)*nactive/beta)/MU_XI; a += SIGMA*sqrt(DT_PARTICLE)*gaussian(); @@ -477,7 +565,7 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA move = 0; for (j=0; j 0) + { + fx = fx/(double)nactive; + fy = fy/(double)nactive; + } + if (FLOOR_FORCE) + { + if (fx > FMAX) fx = FMAX; + else if (fx < -FMAX) fx = -FMAX; + if (fy > FMAX) fy = FMAX; + else if (fy < -FMAX) fy = -FMAX; + } + vxsegments += fx*DT_PARTICLE/(SEGMENTS_MASS); + vysegments += fy*DT_PARTICLE/(SEGMENTS_MASS); + xsegments += vxsegments*DT_PARTICLE; + ysegments += vysegments*DT_PARTICLE; + + /* to avoid numerical instabilities */ + if (xsegments + 1.0 > BCXMAX) + { + xsegments = BCXMAX - 1.0; + vxsegments = 0.0; + } + + translate_segments(segment, xsegments, ysegments); +} + + void animation() { double time, scale, diss, rgb[3], dissip, gradient[2], x, y, dx, dy, dt, xleft, xright, a, b, @@ -584,13 +711,15 @@ void animation() static short int first = 1; t_particle *particle; t_obstacle *obstacle; + t_segment *segment; t_tracer *trajectory; t_hashgrid *hashgrid; char message[100]; particle = (t_particle *)malloc(NMAXCIRCLES*sizeof(t_particle)); /* particles */ - if (ADD_FIXED_OBSTACLES) obstacle = (t_obstacle *)malloc(NMAXOBSTACLES*sizeof(t_obstacle)); /* obstacles */ + if (ADD_FIXED_OBSTACLES) obstacle = (t_obstacle *)malloc(NMAXOBSTACLES*sizeof(t_obstacle)); /* circular obstacles */ + if (ADD_FIXED_SEGMENTS) segment = (t_segment *)malloc(NMAXSEGMENTS*sizeof(t_segment)); /* linear obstacles */ if (TRACER_PARTICLE) trajectory = (t_tracer *)malloc(TRAJECTORY_LENGTH*N_TRACER_PARTICLES*sizeof(t_tracer)); @@ -611,10 +740,11 @@ void animation() xshift = OBSTACLE_XMIN; if (ADD_FIXED_OBSTACLES) init_obstacle_config(obstacle); + if (ADD_FIXED_SEGMENTS) init_segment_config(segment); if (RECORD_PRESSURES) for (i=0; i INITIAL_TIME + SEGMENT_DEACTIVATION_TIME)) + segment[nsegments-1].active = 0; blank(); @@ -662,10 +800,18 @@ void animation() xmincontainer = obstacle_schedule_smooth(i, n); xshift = xmincontainer; } + if ((ROTATE_BOUNDARY)&&(SMOOTH_ROTATION)) rotate_segments(segment, rotation_schedule_smooth(i,n)); + if (INCREASE_GRAVITY) gravity = gravity_schedule(i,n); if ((BOUNDARY_COND == BC_RECTANGLE_WALL)&&(i < INITIAL_TIME + WALL_TIME)) wall = 1; else wall = 0; + if (MOVE_BOUNDARY) for (j=0; jN_T_AVERAGE)) { @@ -770,16 +919,19 @@ void animation() if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length); draw_particles(particle, PLOT); - draw_container(xmincontainer, xmaxcontainer, obstacle, wall); + draw_container(xmincontainer, xmaxcontainer, obstacle, segment, wall); /* add a particle */ - if ((ADD_PARTICLES)&&((i - INITIAL_TIME - ADD_TIME + 1)%ADD_PERIOD == 0)&&(i < NSTEPS - FINAL_NOADD_PERIOD)) + if ((ADD_PARTICLES)&&(i > ADD_TIME)&&((i - INITIAL_TIME - ADD_TIME + 1)%ADD_PERIOD == 0)&&(i < NSTEPS - FINAL_NOADD_PERIOD)) nadd_particle = add_particles(particle, px, py, nadd_particle); + + /* case of reaction-diffusion equation */ + if (REACTION_DIFFUSION) update_types(particle); update_hashgrid(particle, hashgrid, 1); print_parameters(beta, mean_energy, krepel, xmaxcontainer - xmincontainer, - fboundary/(double)(ncircles*NVID), 0, pressure, gravity); + fboundary/(double)(ncircles*NVID), PRINT_LEFT, pressure, gravity); if ((BOUNDARY_COND == BC_EHRENFEST)||(BOUNDARY_COND == BC_RECTANGLE_WALL)) print_ehrenfest_parameters(particle, pleft, pright); else if (PRINT_PARTICLE_NUMBER) print_particle_number(ncircles); @@ -790,6 +942,9 @@ void animation() printf("Entropy 1 = %.5lg, Entropy 2 = %.5lg\n", entropy[0], entropy[1]); print_entropy(entropy); } + + if (PRINT_OMEGA) print_omega(angular_speed); + else if (PRINT_PARTICLE_SPEEDS) print_particles_speeds(particle); glutSwapBuffers(); @@ -811,11 +966,13 @@ void animation() { if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length); draw_particles(particle, PLOT_B); - draw_container(xmincontainer, xmaxcontainer, obstacle, wall); + draw_container(xmincontainer, xmaxcontainer, obstacle, segment, wall); print_parameters(beta, mean_energy, krepel, xmaxcontainer - xmincontainer, - fboundary/(double)(ncircles*NVID), 0, pressure, gravity); + fboundary/(double)(ncircles*NVID), PRINT_LEFT, pressure, gravity); if (BOUNDARY_COND == BC_EHRENFEST) print_ehrenfest_parameters(particle, pleft, pright); else if (PRINT_PARTICLE_NUMBER) print_particle_number(ncircles); + if (PRINT_OMEGA) print_omega(angular_speed); + else if (PRINT_PARTICLE_SPEEDS) print_particles_speeds(particle); glutSwapBuffers(); save_frame_lj_counter(NSTEPS + MID_FRAMES + 1 + counter); counter++; @@ -839,11 +996,13 @@ void animation() { if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length); draw_particles(particle, PLOT); - draw_container(xmincontainer, xmaxcontainer, obstacle, wall); + draw_container(xmincontainer, xmaxcontainer, obstacle, segment, wall); print_parameters(beta, mean_energy, krepel, xmaxcontainer - xmincontainer, - fboundary/(double)(ncircles*NVID), 0, pressure, gravity); + fboundary/(double)(ncircles*NVID), PRINT_LEFT, pressure, gravity); if (BOUNDARY_COND == BC_EHRENFEST) print_ehrenfest_parameters(particle, pleft, pright); else if (PRINT_PARTICLE_NUMBER) print_particle_number(ncircles); + if (PRINT_OMEGA) print_omega(angular_speed); + else if (PRINT_PARTICLE_SPEEDS) print_particles_speeds(particle); glutSwapBuffers(); } for (i=0; i +#include +#include +#include +#include +#include +#include /* Sam Leffler's libtiff library. */ +#include +#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 */ + +/* General geometrical parameters */ + +// #define WINWIDTH 1920 /* window width */ +// #define WINHEIGHT 1000 /* window height */ +// #define NX 480 /* number of grid points on x axis */ +// #define NY 250 /* number of grid points on y axis */ +// // #define NX 1920 /* number of grid points on x axis */ +// // #define NY 1000 /* number of grid points on y axis */ +// +// #define XMIN -2.0 +// #define XMAX 2.0 /* x interval */ +// #define YMIN -1.041666667 +// #define YMAX 1.041666667 /* y interval for 9/16 aspect ratio */ + +#define WINWIDTH 1280 /* window width */ +#define WINHEIGHT 720 /* window height */ + +// #define NX 160 /* number of grid points on x axis */ +// #define NY 90 /* number of grid points on y axis */ +#define NX 320 /* number of grid points on x axis */ +#define NY 180 /* number of grid points on y axis */ + +// #define NX 640 /* number of grid points on x axis */ +// #define NY 360 /* number of grid points on y axis */ + +// #define NX 1280 /* number of grid points on x axis */ +// #define NX 720 /* number of grid points on x axis */ +// #define NY 720 /* number of grid points on y axis */ + +#define XMIN -2.0 +#define XMAX 2.0 /* x interval */ +#define YMIN -1.125 +#define YMAX 1.125 /* y interval for 9/16 aspect ratio */ + +/* Choice of simulated equation */ + +#define RDE_EQUATION 5 /* choice of reaction term, see list in global_3d.c */ +#define NFIELDS 2 /* number of fields in reaction-diffusion equation */ +#define NLAPLACIANS 2 /* number of fields for which to compute Laplacian */ + +#define ADD_POTENTIAL 0 /* set to 1 to add a potential (for Schrodiner equation) */ +#define POTENTIAL 1 /* type of potential, see list in global_3d.c */ + +#define JULIA_SCALE 0.5 /* scaling for Julia sets */ + +/* Choice of the billiard table */ + +#define B_DOMAIN 9 /* choice of domain shape, see list in global_pdes.c */ + +#define CIRCLE_PATTERN 99 /* pattern of circles, see list in global_pdes.c */ + +#define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */ +#define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */ +#define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */ + +#define LAMBDA 0.4 /* parameter controlling the dimensions of domain */ +#define MU 0.06 /* parameter controlling the dimensions of domain */ +#define NPOLY 6 /* number of sides of polygon */ +#define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */ +#define MDEPTH 5 /* depth of computation of Menger gasket */ +#define MRATIO 5 /* ratio defining Menger gasket */ +#define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ +#define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */ +#define FOCI 1 /* set to 1 to draw focal points of ellipse */ +#define NGRIDX 15 /* number of grid point for grid of disks */ +#define NGRIDY 20 /* number of grid point for grid of disks */ + +#define X_SHOOTER -0.2 +#define Y_SHOOTER -0.6 +#define X_TARGET 0.4 +#define Y_TARGET 0.7 /* shooter and target positions in laser fight */ + +#define ISO_XSHIFT_LEFT -1.65 +#define ISO_XSHIFT_RIGHT 0.4 +#define ISO_YSHIFT_LEFT -0.05 +#define ISO_YSHIFT_RIGHT -0.05 +#define ISO_SCALE 0.85 /* coordinates for isospectral billiards */ + +/* You can add more billiard tables by adapting the functions */ +/* xy_in_billiard and draw_billiard in sub_wave.c */ + +/* Physical patameters of wave equation */ + +#define DT 0.00000001 + +#define VISCOSITY 0.1 + +#define RPSA 0.75 /* parameter in Rock-Paper-Scissors-type interaction */ +#define RPSLZB 0.75 /* second parameter in Rock-Paper-Scissors-Lizard-Spock type interaction */ + +#define EPSILON 0.8 /* time scale separation */ +#define DELTA 0.1 /* time scale separation */ +#define FHNA 1.0 /* parameter in FHN equation */ +#define FHNC -0.01 /* parameter in FHN equation */ +#define K_HARMONIC 1.5 /* spring constant of harmonic potential */ +#define BZQ 0.0008 /* parameter in BZ equation */ +#define BZF 1.2 /* parameter in BZ equation */ + +#define T_OUT 2.0 /* outside temperature */ +#define T_IN 0.0 /* inside temperature */ +#define SPEED 0.0 /* speed of drift to the right */ + +#define ADD_NOISE 0 /* set to 1 to add noise, set to 2 to add noise in right half */ +#define NOISE_INTENSITY 0.1 /* noise intensity */ + +#define CHANGE_VISCOSITY 0 /* set to 1 to change the viscosity in the course of the simulation */ +#define ADJUST_INTSTEP 1 /* set to 1 to decrease integration step when viscosity increases */ +#define VISCOSITY_INITIAL_TIME 10 /* initial time during which viscosity remains constant */ +#define VISCOSITY_FACTOR 100.0 /* factor by which to change viscosity */ +#define VISCOSITY_MAX 2.0 /* max value of viscosity beyond which NVID is increased and integration step is decrase, + for numerical stability */ +#define CHANGE_RPSLZB 0 /* set to 1 to change second parameter in Rock-Paper-Scissors-Lizard-Spock equation */ +#define RPSLZB_CHANGE 0.75 /* factor by which to rpslzb parameter */ +#define RPSLZB_INITIAL_TIME 0 /* initial time during which rpslzb remains constant */ +#define RPSLZB_FINAL_TIME 500 /* final time during which rpslzb remains constant */ + +/* Boundary conditions, see list in global_pdes.c */ + +#define B_COND 1 + +/* Parameters for length and speed of simulation */ + +#define NSTEPS 1500 /* number of frames of movie */ +#define NVID 900 /* number of iterations between images displayed on screen */ +#define ACCELERATION_FACTOR 1.0 /* factor by which to increase NVID in course of simulation */ +#define DT_ACCELERATION_FACTOR 1.0 /* factor by which to increase time step in course of simulation */ +#define MAX_DT 0.024 /* maximal value of integration step */ +#define NSEG 100 /* number of segments of boundary */ +#define BOUNDARY_WIDTH 1 /* width of billiard boundary */ + +#define PAUSE 100 /* number of frames after which to pause */ +#define PSLEEP 2 /* sleep time during pause */ +#define SLEEP1 2 /* initial sleeping time */ +#define SLEEP2 1 /* final sleeping time */ +#define INITIAL_TIME 0 /* initial still time */ +#define MID_FRAMES 50 /* number of still frames between parts of two-part movie */ +#define END_FRAMES 100 /* number of still frames at end of movie */ +#define FADE 1 /* set to 1 to fade at end of movie */ + +#define PLOT_3D 1 /* controls whether plot is 2D or 3D */ + +/* Plot type - color scheme */ + +#define CPLOT 30 +#define CPLOT_B 31 + +/* Plot type - height of 3D plot */ + +#define ZPLOT 30 /* z coordinate in 3D plot */ +#define ZPLOT_B 30 /* z coordinate in second 3D plot */ + +#define AMPLITUDE_HIGH_RES 1 /* set to 1 to increase resolution of P_3D_AMPLITUDE plot */ +#define SHADE_3D 1 /* set to 1 to change luminosity according to normal vector */ +#define NON_DIRICHLET_BC 0 /* set to 1 to draw only facets in domain, if field is not zero on boundary */ +#define WRAP_ANGLE 1 /* experimental: wrap angle to [0, 2Pi) for interpolation in angle schemes */ +#define FADE_IN_OBSTACLE 0 /* set to 1 to fade color inside obstacles */ +#define DRAW_OUTSIDE_GRAY 1 /* experimental - draw outside of billiard in gray */ + +#define PLOT_SCALE_ENERGY 0.05 /* vertical scaling in energy plot */ + +#define PRINT_TIME 0 /* set to 1 to print running time */ +#define PRINT_VISCOSITY 0 /* set to 1 to print viscosity */ +#define PRINT_RPSLZB 0 /* set to 1 to print rpslzb parameter */ + +#define DRAW_FIELD_LINES 0 /* set to 1 to draw field lines */ +#define FIELD_LINE_WIDTH 1 /* width of field lines */ +#define N_FIELD_LINES 120 /* number of field lines */ +#define FIELD_LINE_FACTOR 120 /* factor controlling precision when computing origin of field lines */ +#define DRAW_BILLIARD 1 /* set to 1 to draw boundary */ +#define DRAW_BILLIARD_FRONT 1 /* set to 1 to draw boundary */ +#define FILL_BILLIARD_COMPLEMENT 1 /* set to 1 to fill complement of billiard (for certain shapes only) */ + +/* 3D representation */ + +#define REPRESENTATION_3D 1 /* choice of 3D representation */ + +#define REP_AXO_3D 0 /* linear projection (axonometry) */ +#define REP_PROJ_3D 1 /* projection on plane orthogonal to observer line of sight */ + +/* Color schemes, see list in global_pdes.c */ + +#define COLOR_PALETTE 11 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE_B 17 /* Color palette, see list in global_pdes.c */ + +#define BLACK 1 /* black background */ + +#define COLOR_SCHEME 3 /* choice of color scheme */ + +#define COLOR_PHASE_SHIFT 1.0 /* phase shift of color scheme, in units of Pi */ + +#define SCALE 0 /* set to 1 to adjust color scheme to variance of field */ +#define SLOPE 1.0 /* sensitivity of color on wave amplitude */ +#define VSCALE_AMPLITUDE 20.0 /* additional scaling factor for color scheme P_3D_AMPLITUDE */ +#define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ +#define CURL_SCALE 0.000015 /* scaling factor for curl representation */ +#define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */ +#define SLOPE_SCHROD_LUM 50.0 /* sensitivity of luminosity on module, for color scheme Z_ARGUMENT */ + +#define COLORHUE 260 /* initial hue of water color for scheme C_LUM */ +#define COLORDRIFT 0.0 /* how much the color hue drifts during the whole simulation */ +#define LUMMEAN 0.5 /* amplitude of luminosity variation for scheme C_LUM */ +#define LUMAMP 0.3 /* amplitude of luminosity variation for scheme C_LUM */ +#define HUEMEAN 359.0 /* mean value of hue for color scheme C_HUE */ +#define HUEAMP -359.0 /* amplitude of variation of hue for color scheme C_HUE */ +#define E_SCALE 100.0 /* scaling factor for energy representation */ +#define LOG_SCALE 1.0 /* scaling factor for energy log representation */ +#define LOG_SHIFT 0.0 + +#define DRAW_COLOR_SCHEME 0 /* set to 1 to plot the color scheme */ +#define COLORBAR_RANGE 0.55 /* scale of color scheme bar */ +#define COLORBAR_RANGE_B 12.0 /* scale of color scheme bar for 2nd part */ +#define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ + +/* only for compatibility with wave_common.c */ +#define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */ +#define OMEGA 0.005 /* frequency of periodic excitation */ +#define COURANT 0.08 /* Courant number */ +#define COURANTB 0.03 /* Courant number in medium B */ +#define INITIAL_AMP 0.5 /* amplitude of initial condition */ +#define INITIAL_VARIANCE 0.0002 /* variance of initial condition */ +#define INITIAL_WAVELENGTH 0.1 /* wavelength of initial condition */ +#define VSCALE_ENERGY 200.0 /* additional scaling factor for color scheme P_3D_ENERGY */ +#define PHASE_FACTOR 20.0 /* factor in computation of phase in color scheme P_3D_PHASE */ +#define PHASE_SHIFT 0.0 /* shift of phase in color scheme P_3D_PHASE */ +/* end of constants added only for compatibility with wave_common.c */ + + +double u_3d[2] = {0.75, -0.45}; /* projections of basis vectors for REP_AXO_3D representation */ +double v_3d[2] = {-0.75, -0.45}; +double w_3d[2] = {0.0, 0.015}; +double light[3] = {0.816496581, -0.40824829, 0.40824829}; /* vector of "light" direction for P_3D_ANGLE color scheme */ +double observer[3] = {8.0, 8.0, 8.0}; /* location of observer for REP_PROJ_3D representation */ + +#define Z_SCALING_FACTOR 0.75 /* overall scaling factor of z axis for REP_PROJ_3D representation */ +#define XY_SCALING_FACTOR 2.0 /* overall scaling factor for on-screen (x,y) coordinates after projection */ +#define ZMAX_FACTOR 1.0 /* max value of z coordinate for REP_PROJ_3D representation */ +#define XSHIFT_3D 0.0 /* overall x shift for REP_PROJ_3D representation */ +#define YSHIFT_3D 0.05 /* overall y shift for REP_PROJ_3D representation */ + +/* For debugging purposes only */ +#define FLOOR 1 /* set to 1 to limit wave amplitude to VMAX */ +#define VMAX 10.0 /* max value of wave amplitude */ + +#define REFRESH_B (ZPLOT_B != ZPLOT)||(CPLOT_B != CPLOT) /* to save computing time, to be improved */ +#define COMPUTE_WRAP_ANGLE ((WRAP_ANGLE)&&((cplot == Z_ANGLE_GRADIENT)||(Z_ANGLE_GRADIENTX)||(cplot == Z_ARGUMENT)||(Z_ANGLE_GRADIENTX))) + +#include "global_pdes.c" +#include "sub_wave.c" +#include "wave_common.c" /* common functions for wave_billiard, wave_comparison, etc */ + +#include "global_3d.c" /* constants and global variables */ +#include "sub_wave_3d_rde.c" /* should be later replaced by sub_wave_rde.c */ + +#include "sub_rde.c" + + +double potential(int i, int j) +/* compute potential (e.g. for Schrödinger equation) */ +{ + double x, y, xy[2]; + + ij_to_xy(i, j, xy); + x = xy[0]; + y = xy[1]; + + switch (POTENTIAL) { + case (POT_HARMONIC): + { + return (K_HARMONIC*(x*x + y*y)); + } + default: + { + return(0.0); + } + } +} + + +void initialize_potential(double potential_field[NX*NY]) +/* initialize the potential field, e.f. for the Schrödinger equation */ +{ + int i, j; + + #pragma omp parallel for private(i,j) + for (i=0; i VMAX) phi_out[k][i*NY+j] = VMAX; + if (phi_out[k][i*NY+j] < -VMAX) phi_out[k][i*NY+j] = -VMAX; + } + } + } + + for (i=0; i 1280) + { + boxheight = 0.035; + boxwidth = 0.21; + if (left) + { + xbox = XMIN + 0.4; + xtext = XMIN + 0.2; + } + else + { + xbox = XMAX - 0.39; + xtext = XMAX - 0.55; + } + } + else + { + boxwidth = 0.3; + boxheight = 0.05; + if (left) + { + xbox = XMIN + 0.4; + xtext = XMIN + 0.1; + } + else + { + xbox = XMAX - 0.39; + xtext = XMAX - 0.61; + } + } + + first = 0; + } + + y = YMAX - 0.1; + erase_area_hsl(xbox, y + 0.02, boxwidth, boxheight, 0.0, 0.9, 0.0); + glColor3f(1.0, 1.0, 1.0); + if (PRINT_TIME) sprintf(message, "Time %.3f", time); + else if (PRINT_VISCOSITY) sprintf(message, "Viscosity %.3f", viscosity); + else if (PRINT_RPSLZB) sprintf(message, "b = %.3f", rpslzb); + if (PLOT_3D) write_text(xtext, y, message); + else + { + xy_to_pos(xtext, y, pos); + write_text(pos[0], pos[1], message); + } +} + +void draw_color_bar_palette(int plot, double range, int palette, int fade, double fade_value) +{ + double width = 0.14; +// double width = 0.2; + + if (ROTATE_COLOR_SCHEME) + draw_color_scheme_palette_3d(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range, palette, fade, fade_value); + else + draw_color_scheme_palette_3d(XMAX - 1.5*width, YMIN + 0.1, XMAX - 0.5*width, YMAX - 0.1, plot, -range, range, palette, fade, fade_value); +} + + +double viscosity_schedule(int i) +{ + double ratio; + + if (i < VISCOSITY_INITIAL_TIME) return (VISCOSITY); + else + { + ratio = (double)(i - VISCOSITY_INITIAL_TIME)/(double)(NSTEPS - VISCOSITY_INITIAL_TIME); + return (VISCOSITY*(1.0 + ratio*(VISCOSITY_FACTOR - 1.0))); + } +} + +double rpslzb_schedule(int i) +{ + double ratio; + + if (i < RPSLZB_INITIAL_TIME) return (RPSLZB); + else if (i > NSTEPS - RPSLZB_FINAL_TIME) return(RPSLZB - RPSLZB_CHANGE); + else + { + ratio = (double)(i - RPSLZB_INITIAL_TIME)/(double)(NSTEPS - RPSLZB_INITIAL_TIME - RPSLZB_FINAL_TIME); + return (RPSLZB - ratio*RPSLZB_CHANGE); + } +} + + + +void animation() +{ + double time = 0.0, scale, dx, var, jangle, cosj, sinj, sqrintstep, intstep0, viscosity_printed, fade_value; + double *phi[NFIELDS], *phi_tmp[NFIELDS], *potential_field; + short int *xy_in; + int i, j, k, s, nvid, field; + static int counter = 0; + t_rde *rde; + + /* Since NX and NY are big, it seemed wiser to use some memory allocation here */ + for (i=0; i VISCOSITY_MAX)) + { + nvid = (int)((double)NVID*viscosity/VISCOSITY_MAX); +// viscosity = VISCOSITY_MAX; + intstep = intstep0*VISCOSITY_MAX/viscosity; + printf("Nvid = %i, intstep = %.3lg\n", nvid, intstep); + } + } + if (CHANGE_RPSLZB) rpslzb = rpslzb_schedule(i); + + printf("Drawing wave %i\n", i); + draw_wave_rde(0, phi, xy_in, rde, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1); + +// nvid = (int)((double)NVID*(1.0 + (ACCELERATION_FACTOR - 1.0)*(double)i/(double)NSTEPS)); + /* increase integration step */ +// intstep = intstep0*exp(log(DT_ACCELERATION_FACTOR)*(double)i/(double)NSTEPS); +// if (intstep > MAX_DT) +// { +// nvid *= intstep/MAX_DT; +// intstep = MAX_DT; +// } +// printf("Steps per frame: %i\n", nvid); +// printf("Integration step %.5lg\n", intstep); + + printf("Evolving wave\n"); + for (j=0; j= INITIAL_TIME)&&(DOUBLE_MOVIE)) + { + draw_wave_rde(1, phi, xy_in, rde, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 0, 1.0, REFRESH_B); +// draw_billiard(); + if ((PRINT_TIME)||(PRINT_VISCOSITY)||(PRINT_RPSLZB)) print_parameters(time, 0, viscosity_printed); + if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0); + glutSwapBuffers(); + save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter); + counter++; + } + + /* it seems that saving too many files too fast can cause trouble with the file system */ + /* so this is to make a pause from time to time - parameter PAUSE may need adjusting */ + if (i % PAUSE == PAUSE - 1) + { + printf("Making a short pause\n"); + sleep(PSLEEP); + s = system("mv wave*.tif tif_bz/"); + } + } + else printf("Computing frame %i\n", i); + + } + + if (MOVIE) + { + if (DOUBLE_MOVIE) + { + draw_wave_rde(0, phi, xy_in, rde, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1); +// draw_billiard(); + if ((PRINT_TIME)||(PRINT_VISCOSITY)||(PRINT_RPSLZB)) print_parameters(time, 0, viscosity_printed); + if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 0, 1.0); + glutSwapBuffers(); + + if (!FADE) for (i=0; i 2.0) - printf("Grid cell (%i, %i) - neighbour %i = %i = (%i, %i)\n", i, j, k, m, p, q); - } -// sleep(1); - } -// sleep(1); - } +// for (i=0; i 2.0) +// printf("Grid cell (%i, %i) - neighbour %i = %i = (%i, %i)\n", i, j, k, m, p, q); +// } +// // sleep(1); +// } +// // sleep(1); +// } sleep(1); } diff --git a/sub_lj.c b/sub_lj.c index aeb3c92..1cb2eb5 100644 --- a/sub_lj.c +++ b/sub_lj.c @@ -837,7 +837,7 @@ void init_people_config(t_person people[NMAXCIRCLES]) } void init_obstacle_config(t_obstacle obstacle[NMAXOBSTACLES]) -/* initialise particle configuration */ +/* initialise circular obstacle configuration */ { int i, j, n; double x, y, dx, dy; @@ -898,6 +898,354 @@ void init_obstacle_config(t_obstacle obstacle[NMAXOBSTACLES]) } +void init_segment_config(t_segment segment[NMAXSEGMENTS]) +/* initialise linear obstacle configuration */ +{ + int i, cycle = 0, iminus, iplus; + double angle, angle2, dx, width, height, a, b, length; + + switch (SEGMENT_PATTERN) { + case (S_RECTANGLE): + { + segment[0].x1 = BCXMIN; + segment[0].y1 = BCYMAX; + + segment[1].x1 = BCXMIN; + segment[1].y1 = BCYMIN; + + segment[2].x1 = BCXMAX; + segment[2].y1 = BCYMIN; + + segment[3].x1 = BCXMAX; + segment[3].y1 = BCYMAX; + + cycle = 1; + nsegments = 4; + + for (i=0; i= BCXMAX) return(0); + if (x <= BCXMIN) return(0); + if (y >= BCYMAX) return(0); + if (y <= BCYMIN) return(0); + + switch (SEGMENT_PATTERN) { + case (S_CUP): + { + angle = APOLY*PID; + dx = (BCYMAX - BCYMIN)/tan(angle); + + if (y < BCYMAX - (BCYMAX - BCYMIN)*(x - BCXMIN)/dx) return(0); + if (y < BCYMAX - (BCYMAX - BCYMIN)*(BCXMAX - x)/dx) return(0); + } + case (S_HOURGLASS): + { + angle = APOLY*PID; + width = 2.5*MU; + height = 2.5*MU; + + x = vabs(x); + y = vabs(y); + + if ((x >= width)&&(x - width >= (y - height)*(BCXMAX - width)/(BCYMAX - height))) return(0); + return(1); + } + case (S_PENTA): + { + height = 0.5*(BCYMAX - BCYMIN); + width = height/sqrt(3.0); + + if (y < BCYMIN + height*(1.0 - (x - BCXMIN)/width)) return(0); + if (y > BCYMAX - height*(1.0 - (x - BCXMIN)/width)) return(0); + + return(1); + } + case (S_CENTRIFUGE): + { + angle = argument(x,y); + theta = DPI/(double)NPOLY; + while (angle > theta) angle -= theta; + while (angle < 0.0) angle += theta; + if (angle < 0.1) return(0); + if (angle > 0.9) return(0); + return(1); + } + case (S_POLY_ELLIPSE): + { + if (x*x + y*y/(LAMBDA*LAMBDA) < 0.95) return(1); + else return(0); + } + default: return(1); + } +} + + +void rotate_segments(t_segment segment[NMAXSEGMENTS], double angle) +/* rotates the repelling segments by given angle */ +{ + int i; + double ca, sa; + + ca = cos(angle); + sa = sin(angle); + + for (i=0; i DPI) + { + segment[i].angle1 -= DPI; + segment[i].angle2 -= DPI; + } + } + } + +} + +void translate_segments(t_segment segment[NMAXSEGMENTS], double deltax, double deltay) +/* rotates the repelling segments by given angle */ +{ + int i; + + for (i=0; i= 0.0)) -// for (i=-1; i<2; i++) -// for (j=-1; j<2; j++) -// { -// xt1 = x1 + (double)i*periodx; -// yt1 = y1 + (double)j*periody; -// xt2 = x2 + (double)i*periodx; -// yt2 = y2 + (double)j*periody; -// draw_line(xt1, yt1, xt2, yt2); -// } -// else if ((x1 >= 0.0)&&(y1 < 0.0)) -// for (i=-1; i<2; i++) -// for (j=-1; j<2; j++) -// { -// xt1 = x1 + (double)i*periodx; -// yt1 = y1 + (double)j*periody; -// xt2 = x2 + (double)i*periodx; -// yt2 = y2 + (double)j*periody; -// draw_line(xt1, yt1, xt2, yt2); -// } -// } -// } -// } - -// sprintf(message, "%i - %i", particle[j].hash_nneighb, particle[j].hashcell); -// write_text(particle[j].xc, particle[j].yc, message); -// } } void draw_one_particle(t_particle particle, double xc, double yc, double radius, double angle, int nsides, double width, double rgb[3]) @@ -1969,7 +2227,7 @@ void draw_trajectory(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTICLES], void draw_particles(t_particle particle[NMAXCIRCLES], int plot) { int i, j, k, m, width, nnbg, nsides; - double ej, hue, huex, huey, rgb[3], rgbx[3], rgby[3], radius, x1, y1, x2, y2, angle, ca, sa, length, linkcolor, sign = 1.0, angle1, signy = 1.0, periodx, periody, x, y; + double ej, hue, huex, huey, rgb[3], rgbx[3], rgby[3], radius, x1, y1, x2, y2, angle, ca, sa, length, linkcolor, sign = 1.0, angle1, signy = 1.0, periodx, periody, x, y, lum; char message[100]; if (!TRACER_PARTICLE) blank(); @@ -1980,32 +2238,6 @@ void draw_particles(t_particle particle[NMAXCIRCLES], int plot) { glLineWidth(LINK_WIDTH); for (j=0; j DPI) hue -= DPI; + if (hue < 0.0) hue += DPI; + hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*hue/DPI; + if (particle[j].energy < 0.1*PARTICLE_EMAX) lum = 10.0*particle[j].energy/PARTICLE_EMAX; + else lum = 1.0; + radius = particle[j].radius; + width = BOUNDARY_WIDTH; + break; + } case (P_ANGULAR_SPEED): { hue = 160.0*(1.0 + tanh(SLOPE*particle[j].omega)); @@ -2090,6 +2334,15 @@ void draw_particles(t_particle particle[NMAXCIRCLES], int plot) width = BOUNDARY_WIDTH; break; } + case (P_DIFF_NEIGHB): + { + hue = (double)(particle[j].diff_neighb+1)/(double)(particle[j].neighb+1); +// hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*hue; + hue = 180.0*(1.0 + hue); + radius = particle[j].radius; + width = BOUNDARY_WIDTH; + break; + } } switch (particle[j].interaction) { @@ -2139,6 +2392,26 @@ void draw_particles(t_particle particle[NMAXCIRCLES], int plot) hsl_to_rgb_twilight(hue, 0.9, 0.5, rgby); break; } + case (P_DIRECT_ENERGY): + { + hsl_to_rgb_twilight(hue, 0.9, 0.5, rgb); + hsl_to_rgb_twilight(hue, 0.9, 0.5, rgbx); + hsl_to_rgb_twilight(hue, 0.9, 0.5, rgby); + for (i=0; i<3; i++) + { + rgb[i] *= lum; + rgbx[i] *= lum; + rgby[i] *= lum; + } + break; + } + case (P_DIFF_NEIGHB): + { + hsl_to_rgb_twilight(hue, 0.9, 0.5, rgb); + hsl_to_rgb_twilight(hue, 0.9, 0.5, rgbx); + hsl_to_rgb_twilight(hue, 0.9, 0.5, rgby); + break; + } default: { hsl_to_rgb(hue, 0.9, 0.5, rgb); @@ -2257,7 +2530,7 @@ void draw_particles(t_particle particle[NMAXCIRCLES], int plot) } -void draw_container(double xmin, double xmax, t_obstacle obstacle[NMAXOBSTACLES], int wall) +void draw_container(double xmin, double xmax, t_obstacle obstacle[NMAXOBSTACLES], t_segment segment[NMAXSEGMENTS], int wall) /* draw the container, for certain boundary conditions */ { int i, j; @@ -2275,6 +2548,13 @@ void draw_container(double xmin, double xmax, t_obstacle obstacle[NMAXOBSTACLES] for (i = 0; i < nobstacles; i++) draw_circle(obstacle[i].xc, obstacle[i].yc, obstacle[i].radius, NSEG); } + if (ADD_FIXED_SEGMENTS) + { + glLineWidth(CONTAINER_WIDTH); + glColor3f(1.0, 1.0, 1.0); + for (i = 0; i < nsegments; i++) if (segment[i].active) + draw_line(segment[i].x1, segment[i].y1, segment[i].x2, segment[i].y2); + } switch (BOUNDARY_COND) { case (BC_SCREEN): @@ -2462,20 +2742,8 @@ void draw_container(double xmin, double xmax, t_obstacle obstacle[NMAXOBSTACLES] hsl_to_rgb(300.0, 0.1, 0.5, rgb); draw_colored_rectangle(0.0, 0.0, BCXMAX, BCYMAX, rgb); } - } -// /* draw fixed obstacles */ -// if (ADD_FIXED_OBSTACLES) -// { -// glLineWidth(CONTAINER_WIDTH); -// hsl_to_rgb(300.0, 0.1, 0.5, rgb); -// for (i = 0; i < nobstacles; i++) -// draw_colored_circle(obstacle[i].xc, obstacle[i].yc, obstacle[i].radius, NSEG, rgb); -// glColor3f(1.0, 1.0, 1.0); -// for (i = 0; i < nobstacles; i++) -// draw_circle(obstacle[i].xc, obstacle[i].yc, obstacle[i].radius, NSEG); -// } } void print_parameters(double beta, double temperature, double krepel, double lengthcontainer, double boundary_force, @@ -2628,7 +2896,7 @@ void print_parameters(double beta, double temperature, double krepel, double len } } - if (PARTIAL_THERMO_COUPLING) + if ((PARTIAL_THERMO_COUPLING)&&(!INCREASE_BETA)) { printf("Temperature %i in average: %.3lg\n", i_temp, temperature); temp[i_temp] = temperature; @@ -2828,15 +3096,77 @@ void print_entropy(double entropy[2]) } +void print_omega(double angular_speed) +{ + char message[100]; + double rgb[3]; + static double xleftbox, xlefttext, xrightbox, xrighttext, y = YMAX - 0.1, ymin = YMIN + 0.05; + static int first = 1; + + if (first) + { +// xleftbox = XMIN + 0.4; +// xlefttext = xleftbox - 0.55; + xrightbox = XMAX - 0.39; + xrighttext = xrightbox - 0.55; + first = 0; + } + + + erase_area_hsl(xrightbox, y + 0.025, 0.35, 0.05, 0.0, 0.9, 0.0); + glColor3f(1.0, 1.0, 1.0); + sprintf(message, "Angular speed = %.4f", DPI*angular_speed*25.0/(double)(PERIOD_ROTATE_BOUNDARY)); + write_text(xrighttext + 0.1, y, message); + +} + +void print_particles_speeds(t_particle particle[NMAXCIRCLES]) +{ + char message[100]; + double y = YMAX - 0.1, vx = 0.0, vy = 0.0; + int i, nactive = 0; + static double xleftbox, xlefttext, xrightbox, xrighttext, ymin = YMIN + 0.05; + static int first = 1; + + if (first) + { + xrightbox = XMAX - 0.39; + xrighttext = xrightbox - 0.55; + first = 0; + } + + for (i=0; i 0.0)&&(proj < 1.0)) + { + distance = segment[i].nx*x + segment[i].ny*y - segment[i].c; + r = 1.5*particle[j].radius; + if (vabs(distance) < r) + { + f = KSPRING_OBSTACLE*(r - distance); + particle[j].fx += f*segment[i].nx; + particle[j].fy += f*segment[i].ny; + if (MOVE_BOUNDARY) + { + segment[i].fx -= f*segment[i].nx; + segment[i].fy -= f*segment[i].ny; +// segment[i].fx += f*segment[i].nx; +// segment[i].fy += f*segment[i].ny; + } + } + } + + /* compute force from concave corners */ + if (segment[i].concave) + { + distance = module2(x - segment[i].x1, y - segment[i].y1); + angle = argument(x - segment[i].x1, y - segment[i].y1); + if (angle < segment[i].angle1) angle += DPI; + r = 1.5*particle[j].radius; + + if ((distance < r)&&(angle > segment[i].angle1)&&(angle < segment[i].angle2)) + { + f = KSPRING_OBSTACLE*(r - distance); + particle[j].fx += f*cos(angle); + particle[j].fy += f*sin(angle); + if (MOVE_BOUNDARY) + { + segment[i].fx -= f*cos(angle); + segment[i].fy -= f*sin(angle); +// segment[i].fx += f*cos(angle); +// segment[i].fy += f*sin(angle); + } + } + } + } switch(BOUNDARY_COND){ case (BC_SCREEN): { /* add harmonic force outside screen */ - if (particle[j].xc > BCXMAX) particle[j].fx -= KSPRING_BOUNDARY*(particle[j].xc - BCXMAX); - else if (particle[j].xc < BCXMIN) particle[j].fx += KSPRING_BOUNDARY*(BCXMIN - particle[j].xc); - if (particle[j].yc > BCYMAX) particle[j].fy -= KSPRING_BOUNDARY*(particle[j].yc - BCYMAX); - else if (particle[j].yc < BCYMIN) particle[j].fy += KSPRING_BOUNDARY*(BCYMIN - particle[j].yc); + if (particle[j].xc > XMAX) particle[j].fx -= KSPRING_BOUNDARY*(particle[j].xc - XMAX); + else if (particle[j].xc < XMIN) particle[j].fx += KSPRING_BOUNDARY*(XMIN - particle[j].xc); + if (particle[j].yc > YMAX) particle[j].fy -= KSPRING_BOUNDARY*(particle[j].yc - YMAX); + else if (particle[j].yc < YMIN) particle[j].fy += KSPRING_BOUNDARY*(YMIN - particle[j].yc); +// if (particle[j].xc > BCXMAX) particle[j].fx -= KSPRING_BOUNDARY*(particle[j].xc - BCXMAX); +// else if (particle[j].xc < BCXMIN) particle[j].fx += KSPRING_BOUNDARY*(BCXMIN - particle[j].xc); +// if (particle[j].yc > BCYMAX) particle[j].fy -= KSPRING_BOUNDARY*(particle[j].yc - BCYMAX); +// else if (particle[j].yc < BCYMIN) particle[j].fy += KSPRING_BOUNDARY*(BCYMIN - particle[j].yc); return(fperp); } case (BC_RECTANGLE): @@ -3235,6 +3617,7 @@ void compute_particle_force(int j, double krepel, t_particle particle[NMAXCIRCLE double fx = 0.0, fy = 0.0, force[2], torque = 0.0, torque_ij, x, y; particle[j].neighb = 0; + if (REACTION_DIFFUSION) particle[j].diff_neighb = 0; for (k=0; k TPYE_PROPORTION)) { - particle[i].type = 1; + particle[i].type = 2; particle[i].radius = MU_B; } particle[i].neighb = 0; + particle[i].diff_neighb = 0; particle[i].thermostat = 1; // particle[i].energy = 0.0; @@ -3462,17 +3851,29 @@ int initialize_configuration(t_particle particle[NMAXCIRCLES], t_hashgrid hashgr } - if (ADD_FIXED_OBSTACLES) { for (i=0; i< ncircles; i++) for (j=0; j < nobstacles; j++) if (module2(particle[i].xc - obstacle[j].xc, particle[i].yc - obstacle[j].yc) < OBSTACLE_RADIUS + particle[i].radius) particle[i].active = 0; } + + /* case of segment obstacles */ + if (ADD_FIXED_SEGMENTS) for (i=0; i< ncircles; i++) + if (!in_segment_region(particle[i].xc, particle[i].yc)) + particle[i].active = 0; + + /* case of reaction-diffusion equation */ + if (REACTION_DIFFUSION) for (i=0; i< ncircles; i++) + { + particle[i].type = 1 + (int)(RD_TYPES*(double)rand()/(double)RAND_MAX); + } /* count number of active particles */ for (i=0; i< ncircles; i++) nactive += particle[i].active; - printf("%i active particles\n", nactive); + printf("%i active particles\n", nactive); + + for (i=0; i xmin to thermostat */ { - int i, nthermo = 0; + int condition, i, nthermo = 0; for (i=0; i xmin) + switch (PARTIAL_THERMO_REGION) { + case (TH_VERTICAL): + { + condition = (particle[i].xc > xmin); + break; + } + case (TH_INSEGMENT): + { + condition = (in_segment_region(particle[i].xc - xsegments, particle[i].yc - ysegments)); + break; + } + default: condition = 1; + } + if (condition) { particle[i].thermostat = 1; nthermo++; @@ -3578,4 +3996,24 @@ double compute_mean_energy(t_particle particle[NMAXCIRCLES]) nactive++; } return(total_energy/(double)nactive); -} \ No newline at end of file +} + +void update_types(t_particle particle[NMAXCIRCLES]) +/* update the types in case of reaction-diffusion equation */ +{ + int i, j, k; + double distance; + + for (i=0; i= 1 || S == 0); + X = V1 * sqrt(-2 * log(S) / S); + } + else X = V2 * sqrt(-2 * log(S) / S); + + phase = 1 - phase; + + return X; +} + +void init_random(double mean, double amplitude, double *phi[NFIELDS], short int xy_in[NX*NY]) +/* initialise field with gaussian at position (x,y) */ +{ + int i, j, k, in; + double xy[2], dist2, module, phase, scale2; + + printf("Initialising xy_in\n"); + #pragma omp parallel for private(i,j) + for (i=0; i= 2) phi[i][j] = T_IN*pow(0.75, (double)(in-2)); +// else if (in >= 2) phi[i][j] = T_IN*pow(1.0 - 0.5*(double)(in-2), (double)(in-2)); +// else if (in >= 2) phi[i][j] = T_IN*(1.0 - (double)(in-2)/((double)MDEPTH))*(1.0 - (double)(in-2)/((double)MDEPTH)); + else phi[i][j] = T_OUT; + } +} + +void init_coherent_state(double x, double y, double px, double py, double scalex, double *phi[NFIELDS], short int xy_in[NX*NY]) +/* initialise field with coherent state of position (x,y) and momentum (px, py) */ +/* phi[0] is real part, phi[1] is imaginary part */ +{ + int i, j; + double xy[2], dist2, module, phase, scale2; + + scale2 = scalex*scalex; + for (i=0; i= DPI) return (angle - DPI); + else return(angle); +} + +double unwrap_angle(double angle1, double angle2) +{ + if (angle2 - angle1 < -PI) return(angle2 + DPI); + else if (angle2 - angle1 >= PI) return (angle2 - DPI); + else return(angle2); +} + + +void compute_theta(double *phi[NFIELDS], short int xy_in[NX*NY], t_rde rde[NX*NY]) +/* compute angle for rock-paper-scissors equation */ +{ + int i, j; + double u, v, w, rho, xc, yc, angle; + + #pragma omp parallel for private(i,j,u, v, w, rho, xc, yc, angle) + for (i=0; i= DPI) angle -= DPI; + rde[i*NY+j].theta = angle; + } + else rde[i*NY+j].theta = 0.0; +} + +double colors_rpslz(int type) +{ + switch (type) { + case (0): return(0.0); + case (1): return(60.0); + case (2): return(120.0); + case (3): return(200.0); + case (4): return(270.0); + } +} + +void compute_theta_rpslz(double *phi[NFIELDS], short int xy_in[NX*NY], t_rde rde[NX*NY], int cplot) +/* compute angle for rock-paper -scissors-lizard-Spock equation */ +{ + int i, j, k, kmax; + double rho, xc, yc, angle, max; + static double sa[5], ca[5], shift; + static int first = 1; + + if (first) + { +// shift = 0.0; + for (i = 0; i < 5; i++) + { + ca[i] = cos(colors_rpslz(i)*DPI/360.0); + sa[i] = sin(colors_rpslz(i)*DPI/360.0); +// ca[i] = cos(shift + 0.2*DPI*(double)i); +// sa[i] = sin(shift + 0.2*DPI*(double)i); + } + first = 0; + } + + switch (cplot) { + case (Z_MAXTYPE_RPSLZ): + { + #pragma omp parallel for private(i,j,max,kmax) + for (i=0; i max) + { + max = phi[k][i*NY+j]; + kmax = k; + } + } + angle = colors_rpslz(kmax); + rde[i*NY+j].theta = angle; + } + else rde[i*NY+j].theta = 0.0; + } + break; + } +// case (Z_THETA_RPSLZ): + default: + { + #pragma omp parallel for private(i,j,rho,xc,yc,angle) + for (i=0; i= DPI) angle -= DPI; + rde[i*NY+j].theta = angle; + } + else rde[i*NY+j].theta = 0.0; + } + break; + } + } +} + + +void compute_gradient_rde(double phi[NX*NY], t_rde rde[NX*NY]) +/* compute the gradient of the field */ +{ + int i, j, iplus, iminus, jplus, jminus; + double deltaphi; + double dx = (XMAX-XMIN)/((double)NX); + + dx = (XMAX-XMIN)/((double)NX); + + #pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,deltaphi) + for (i=0; i PI) deltaphi -= DPI; + if (vabs(deltaphi) < 1.0e9) rde[i*NY+j].nablax = (deltaphi)/dx; + else rde[i*NY+j].nablax = 0.0; + + deltaphi = phi[i*NY+jplus] - phi[i*NY+jminus]; + if (deltaphi < -PI) deltaphi += DPI; + if (deltaphi > PI) deltaphi -= DPI; + if (vabs(deltaphi) < 1.0e9) rde[i*NY+j].nablay = (deltaphi)/dx; + else rde[i*NY+j].nablay = 0.0; + + /* TO DO: improve this computation */ + rde[i*NY+j].field_norm = module2(rde[i*NY+j].nablax,rde[i*NY+j].nablay); + rde[i*NY+j].field_arg = argument(rde[i*NY+j].nablax,rde[i*NY+j].nablay); + } +} + +void compute_gradient_theta(t_rde rde[NX*NY]) +/* compute the gradient of the theta field */ +{ + int i, j, iplus, iminus, jplus, jminus; + double deltaphi; + double dx = (XMAX-XMIN)/((double)NX); + + dx = (XMAX-XMIN)/((double)NX); + + #pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,deltaphi) + for (i=0; i PI) deltaphi -= DPI; + if (vabs(deltaphi) < 1.0e9) rde[i*NY+j].nablax = (deltaphi)/dx; + else rde[i*NY+j].nablax = 0.0; + + deltaphi = rde[i*NY+jplus].theta - rde[i*NY+jminus].theta; + if (deltaphi < -PI) deltaphi += DPI; + if (deltaphi > PI) deltaphi -= DPI; + if (vabs(deltaphi) < 1.0e9) rde[i*NY+j].nablay = (deltaphi)/dx; + else rde[i*NY+j].nablay = 0.0; + } +} + +void compute_curl(t_rde rde[NX*NY]) +/* compute the curl of the field */ +{ + int i, j, iplus, iminus, jplus, jminus; + double delta; + double dx = (XMAX-XMIN)/((double)NX); + + dx = (XMAX-XMIN)/((double)NX); + + #pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta) + for (i=0; i= DPI) angle -= DPI; + rde[i*NY+j].field_arg = angle; + } +} + +void compute_field_module(double *phi[NFIELDS], t_rde rde[NX*NY], double factor) +/* compute the norm squared of first two fields */ +{ + int i, j; + + #pragma omp parallel for private(i,j) + for (i=0; i= DPI) arg -= DPI; + rde[i*NY+j].field_arg = arg; + } +} + + + +// void compute_field_norm(double phi_x[NX*NY], double phi_y[NX*NY], double phi_norm[NX*NY], double factor) +// /* compute the norm of (phi_x, phi_y) */ +// { +// int i, j; +// +// #pragma omp parallel for private(i,j) +// for (i=0; i= DPI) angle -= DPI; +// phi_arg[i*NY+j] = 360.0*angle/DPI; +// // phi_arg[i*NY+j] = 180.0*angle/DPI; +// } +// } + + + +void compute_laplacian(double phi_in[NX*NY], double phi_out[NX*NY], short int xy_in[NX*NY]) +/* computes the Laplacian of phi_in and stores it in phi_out */ +{ + int i, j, iplus, iminus, jplus, jminus; + + #pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus) + for (i=1; i NX-1) ij[0] = NX-1; + if (ij[1] < 0) ij[1] = 0; + if (ij[1] > NY-1) ij[1] = NY-1; + + nabx = nablax[ij[0]][ij[1]]; + naby = nablay[ij[0]][ij[1]]; + + norm2 = nabx*nabx + naby*naby; + + if (norm2 > 1.0e-14) + { + /* avoid too large step size */ + if (norm2 < 1.0e-9) norm2 = 1.0e-9; + norm = sqrt(norm2); + x1 = x1 + delta*nabx/norm; + y1 = y1 + delta*naby/norm; + } + else cont = 0; + + if (!xy_in[ij[0]*NY + ij[1]]) cont = 0; + + /* stop if the boundary is hit */ +// if (xy_in[ij[0]][ij[1]] != 1) cont = 0; + +// printf("x1 = %.3lg \t y1 = %.3lg \n", x1, y1); + + xy_to_pos(x1, y1, pos); + glVertex2d(pos[0], pos[1]); + + i++; + } + glEnd(); +} + + + +void compute_field_color_rde(double value, int cplot, int palette, double rgb[3]) +/* compute the color depending on the field value and color palette */ +{ + switch (cplot) { + case (Z_AMPLITUDE): + { + color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 0, rgb); + break; + } + case (Z_RGB): + { + /* TO DO */ + break; + } + case (Z_POLAR): + { +// hsl_to_rgb_palette(360.0*value/DPI, 0.9, 0.5, rgb, palette); + color_scheme_palette(C_ONEDIM_LINEAR, palette, value/DPI, 1.0, 1, rgb); + break; + } + case (Z_NORM_GRADIENT): + { + color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 0, rgb); + break; + } + case (Z_ANGLE_GRADIENT): + { + color_scheme_palette(C_ONEDIM_LINEAR, palette, value/DPI, 1.0, 1, rgb); + break; + } + case (Z_NORM_GRADIENTX): + { + hsl_to_rgb_palette(360.0*value/DPI, 0.9, 0.5, rgb, palette); + break; + } + case (Z_ANGLE_GRADIENTX): + { + color_scheme_palette(C_ONEDIM_LINEAR, palette, value/DPI, 1.0, 1, rgb); + break; + } + case (Z_NORM_GRADIENT_INTENSITY): + { + hsl_to_rgb_palette(360.0*value/DPI, 0.9, 0.5, rgb, palette); + break; + } + case (Z_VORTICITY): + { + hsl_to_rgb_palette(180.0*(1.0 - color_amplitude(value, 1.0, 0)), 0.9, 0.5, rgb, palette); + break; + } + case (Z_MAXTYPE_RPSLZ): + { + hsl_to_rgb_palette(value, 0.9, 0.5, rgb, palette); + break; + } + case (Z_THETA_RPSLZ): + { + color_scheme_palette(C_ONEDIM_LINEAR, palette, value/DPI, 1.0, 1, rgb); + break; + } + case (Z_NORM_GRADIENT_RPSLZ): + { + color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 0, rgb); + break; + } + case (Z_MODULE): + { + color_scheme_asym_palette(COLOR_SCHEME, palette, VSCALE_AMPLITUDE*value, 1.0, 0, rgb); + break; + } + case (Z_ARGUMENT): + { + color_scheme_palette(C_ONEDIM_LINEAR, palette, value/DPI, 1.0, 1, rgb); + break; + } + } +} + + +double compute_interpolated_colors_rde(int i, int j, t_rde rde[NX*NY], double palette, int cplot, + double rgb_e[3], double rgb_w[3], double rgb_n[3], double rgb_s[3], + int fade, double fade_value, int movie) +{ + int k; + double cw, ce, cn, cs, c_sw, c_se, c_nw, c_ne, c_mid, ca, z_mid; + double *z_sw, *z_se, *z_nw, *z_ne, lum; + + z_sw = rde[i*NY+j].p_zfield[movie]; + z_se = rde[(i+1)*NY+j].p_zfield[movie]; + z_nw = rde[i*NY+j+1].p_zfield[movie]; + z_ne = rde[(i+1)*NY+j+1].p_zfield[movie]; + + z_mid = 0.25*(*z_sw + *z_se + *z_nw + *z_ne); + + c_sw = *rde[i*NY+j].p_cfield[movie]; + c_se = *rde[(i+1)*NY+j].p_cfield[movie]; + c_nw = *rde[i*NY+j+1].p_cfield[movie]; + c_ne = *rde[(i+1)*NY+j+1].p_cfield[movie]; + + if (COMPUTE_WRAP_ANGLE) + { + c_se = unwrap_angle(c_sw,c_se); + c_nw = unwrap_angle(c_sw,c_nw); + c_ne = unwrap_angle(c_sw,c_ne); + } + + c_mid = 0.25*(c_sw + c_se + c_nw + c_ne); + + cw = (c_sw + c_nw + c_mid)/3.0; + ce = (c_se + c_ne + c_mid)/3.0; + cs = (c_sw + c_se + c_mid)/3.0; + cn = (c_nw + c_ne + c_mid)/3.0; + + if (COMPUTE_WRAP_ANGLE) + { + cw = wrap_angle(cw); + ce = wrap_angle(ce); + cs = wrap_angle(cs); + cn = wrap_angle(cn); + } + + compute_field_color_rde(ce, cplot, palette, rgb_e); + compute_field_color_rde(cw, cplot, palette, rgb_w); + compute_field_color_rde(cn, cplot, palette, rgb_n); + compute_field_color_rde(cs, cplot, palette, rgb_s); + + if (cplot == Z_ARGUMENT) + { + lum = tanh(SLOPE_SCHROD_LUM*rde[i*NY+j].field_norm); + for (k=0; k<3; k++) + { + rgb_e[k] *= lum; + rgb_w[k] *= lum; + rgb_n[k] *= lum; + rgb_s[k] *= lum; + } + } + if (SHADE_3D) + { + ca = rde[i*NY+j].cos_angle; + ca = (ca + 1.0)*0.4 + 0.2; + for (k=0; k<3; k++) + { + rgb_e[k] *= ca; + rgb_w[k] *= ca; + rgb_n[k] *= ca; + rgb_s[k] *= ca; + } + } + if (fade) + for (k=0; k<3; k++) + { + rgb_e[k] *= fade_value; + rgb_w[k] *= fade_value; + rgb_n[k] *= fade_value; + rgb_s[k] *= fade_value; + } + + return(z_mid); +} + + +void compute_rde_fields(double *phi[NFIELDS], short int xy_in[NX*NY], int zplot, int cplot, t_rde rde[NX*NY]) +/* compute the necessary auxiliary fields */ +{ + int i, j; + + switch (RDE_EQUATION) { + case (E_RPS): + { + if ((COMPUTE_THETA)||(COMPUTE_THETAZ)) + compute_theta(phi, xy_in, rde); + + if ((zplot == Z_NORM_GRADIENT)||(cplot == Z_NORM_GRADIENT)) + { + compute_gradient_theta(rde); +// compute_gradient_polar(rde, 1.0); + compute_gradient_polar(rde, 0.003); + } + + if ((zplot == Z_NORM_GRADIENTX)||(cplot == Z_NORM_GRADIENTX)||(zplot == Z_ANGLE_GRADIENTX)||(cplot == Z_ANGLE_GRADIENTX)) + { + compute_gradient_rde(phi[0], rde); + compute_gradient_polar(rde, 0.03); + } + + if ((zplot == Z_VORTICITY)||(cplot == Z_VORTICITY)) + { + compute_gradient_theta(rde); + compute_curl(rde); + } + break; + } + case (E_RPSLZ): + { + compute_theta_rpslz(phi, xy_in, rde, cplot); + if ((zplot == Z_NORM_GRADIENT_RPSLZ)||(cplot == Z_NORM_GRADIENT_RPSLZ)) + { + compute_gradient_theta(rde); + compute_gradient_polar(rde, 0.005); + } + break; + } + case (E_SCHRODINGER): + { + compute_field_module(phi, rde, 1.0); + if ((zplot == Z_ARGUMENT)||(cplot == Z_ARGUMENT)) + { + compute_field_argument(phi, rde); + } + break; + } + default : break; + } +} + +void init_zfield_rde(double *phi[NFIELDS], short int xy_in[NX*NY], int zplot, t_rde rde[NX*NY], int movie) +/* compute the necessary fields for the z coordinate */ +{ + int i, j; + + switch(zplot) { + case (Z_AMPLITUDE): + { + #pragma omp parallel for private(i,j) + for (i=0; i 0) + { + z_mid = compute_interpolated_colors_rde(i, j, rde, palette, cplot, + rgb_e, rgb_w, rgb_n, rgb_s, fade, fade_value, movie); + + ij_to_xy(i, j, xy_sw); + ij_to_xy(i+1, j, xy_se); + ij_to_xy(i, j+1, xy_nw); + ij_to_xy(i+1, j+1, xy_ne); + + for (k=0; k<2; k++) xy_mid[k] = 0.25*(xy_sw[k] + xy_se[k] + xy_nw[k] + xy_ne[k]); + + if (AMPLITUDE_HIGH_RES == 1) + { + glBegin(GL_TRIANGLE_FAN); + glColor3f(rgb_w[0], rgb_w[1], rgb_w[2]); + draw_vertex_xyz(xy_mid, z_mid); + draw_vertex_xyz(xy_nw, *rde[i*NY+j+1].p_zfield[movie]); + draw_vertex_xyz(xy_sw, *rde[i*NY+j].p_zfield[movie]); + + glColor3f(rgb_s[0], rgb_s[1], rgb_s[2]); + draw_vertex_xyz(xy_se, *rde[(i+1)*NY+j].p_zfield[movie]); + + glColor3f(rgb_e[0], rgb_e[1], rgb_e[2]); + draw_vertex_xyz(xy_ne, *rde[(i+1)*NY+j+1].p_zfield[movie]); + + glColor3f(rgb_n[0], rgb_n[1], rgb_n[2]); + draw_vertex_xyz(xy_nw, *rde[i*NY+j+1].p_zfield[movie]); + glEnd (); + } + } + else + { + glColor3f(rde[i*NY+j].rgb[0], rde[i*NY+j].rgb[1], rde[i*NY+j].rgb[2]); + + glBegin(GL_TRIANGLE_FAN); + ij_to_xy(i, j, xy); + draw_vertex_xyz(xy, *rde[i*NY+j].p_zfield[movie]); + ij_to_xy(i+1, j, xy); + draw_vertex_xyz(xy, *rde[(i+1)*NY+j].p_zfield[movie]); + ij_to_xy(i+1, j+1, xy); + draw_vertex_xyz(xy, *rde[(i+1)*NY+j+1].p_zfield[movie]); + ij_to_xy(i, j+1, xy); + draw_vertex_xyz(xy, *rde[i*NY+j+1].p_zfield[movie]); + glEnd (); + } + } + } + + if (DRAW_BILLIARD_FRONT) draw_billiard_3d_front(fade, fade_value); +} + + +void draw_wave_rde(int movie, double *phi[NFIELDS], short int xy_in[NX*NY], t_rde rde[NX*NY], + int zplot, int cplot, int palette, int fade, double fade_value, int refresh) +{ + int i, j, k, l, draw = 1; + double xy[2], xy_screen[2], rgb[3], pos[2], ca, rgb_e[3], rgb_w[3], rgb_n[3], rgb_s[3]; + double z_sw, z_se, z_nw, z_ne, z_mid, zw, ze, zn, zs, min = 1000.0, max = 0.0; + double xy_sw[2], xy_se[2], xy_nw[2], xy_ne[2], xy_mid[2]; + double energy; + + if (refresh) + { +// printf("Computing fields\n"); + compute_rde_fields(phi, xy_in, zplot, cplot, rde); +// printf("Computed fields\n"); + if ((PLOT_3D)&&(SHADE_3D)) compute_light_angle_rde(xy_in, rde, movie); + } + compute_cfield_rde(xy_in, cplot, palette, rde, fade, fade_value, movie); + + if (PLOT_3D) draw_wave_3d_rde(movie, phi, xy_in, rde, zplot, cplot, palette, fade, fade_value); + else draw_wave_2d_rde(xy_in, rde); +} + + diff --git a/sub_wave.c b/sub_wave.c index 0766ff6..8e428ac 100644 --- a/sub_wave.c +++ b/sub_wave.c @@ -1434,10 +1434,14 @@ int compute_noisepanel_coordinates(t_vertex polyline[NMAXPOLY]) if (even) y = YMIN + 0.1; else y = YMIN + 0.1 + MU; - polyline[i].x = x; + x1 = x; + if (x1 > XMAX) x1 = XMAX; + else if (x1 < XMIN) x1 = XMIN; + + polyline[i].x = x1; polyline[i].y = y; - xy_to_pos(x, y, pos); + xy_to_pos(x1, y, pos); polyline[i].posi = pos[0]; polyline[i].posj = pos[1]; @@ -1612,6 +1616,11 @@ int xy_in_billiard(double x, double y) static int first = 1, nsides; switch (B_DOMAIN) { + case (D_NOTHING): + { + return(1); + break; + } case (D_RECTANGLE): { if ((vabs(x) margin*length*olength)&&(x1*observer[0] + y1*observer[1] > margin*length1*olength)) + { + glBegin(GL_LINE_STRIP); + tvertex_lineto_3d(polyline[j]); + tvertex_lineto_3d(polyline[k]); + glEnd(); + } +} + + void draw_billiard_3d_front(int fade, double fade_value) /* hack to draw the billiard boundary in front of the wave */ /* only parts of the boundary having a small enough angle with the observer vector are drawn */ @@ -672,27 +781,51 @@ void draw_billiard_3d_front(int fade, double fade_value) glEnable(GL_LINE_SMOOTH); switch (B_DOMAIN) { + case (D_YOUNG): + { + glBegin(GL_LINE_STRIP); + draw_vertex_x_y_z(-MU, YMIN, 0.0); + draw_vertex_x_y_z(-MU, -LAMBDA-MU, 0.0); +// draw_vertex_x_y_z(MU, -LAMBDA-MU, 0.0); +// draw_vertex_x_y_z(MU, YMIN, 0.0); + glEnd(); + + glBegin(GL_LINE_STRIP); + draw_vertex_x_y_z(-MU, YMAX, 0.0); + draw_vertex_x_y_z(-MU, LAMBDA+MU, 0.0); +// draw_vertex_x_y_z(MU, LAMBDA+MU, 0.0); +// draw_vertex_x_y_z(MU, YMAX, 0.0); + glEnd(); + + glBegin(GL_LINE_LOOP); + draw_vertex_x_y_z(-MU, -LAMBDA+MU, 0.0); + draw_vertex_x_y_z(-MU, LAMBDA-MU, 0.0); +// draw_vertex_x_y_z(MU, LAMBDA-MU, 0.0); +// draw_vertex_x_y_z(MU, -LAMBDA+MU, 0.0); + glEnd(); + break; + } + case (D_TOKA_PRIME): + { + draw_polyline_visible(0, 4, 0.2); + for (i=4; i<41; i++) draw_polyline_visible(i, i+1, -0.1); +// draw_polyline_visible(42, 3, 0.2); + for (i=84; i>46; i--) draw_polyline_visible(i, i-1, -0.1); + draw_polyline_visible(46, 0, 0.2); + break; + } case (D_STAR): { glLineWidth(BOUNDARY_WIDTH); - for (i=0; i 0.2*length*olength)&&(x1*observer[0] + y1*observer[1] > 0.2*length1*olength)) - { - glBegin(GL_LINE_STRIP); - tvertex_lineto_3d(polyline[j]); - tvertex_lineto_3d(polyline[k]); - glEnd(); - } - } + for (i=0; i LOG_ENERGY_FLOOR) wave[i*NY+j].log_energy = LOG_SHIFT + PLOT_SCALE_LOG_ENERGY*logenergy; + else wave[i*NY+j].log_energy = LOG_SHIFT + PLOT_SCALE_LOG_ENERGY*LOG_ENERGY_FLOOR; + } + + if (COMPUTE_LOG_TOTAL_ENERGY) + { + logenergy = log(wave[i*NY+j].total_energy); + if (logenergy > LOG_ENERGY_FLOOR) wave[i*NY+j].log_total_energy = LOG_SHIFT + PLOT_SCALE_LOG_ENERGY*logenergy; + else wave[i*NY+j].log_total_energy = LOG_SHIFT + PLOT_SCALE_LOG_ENERGY*LOG_ENERGY_FLOOR; + } + + if (COMPUTE_MEAN_ENERGY) + { + wave[i*NY+j].mean_energy = wave[i*NY+j].total_energy/((double)(global_time+1)); + } + + if (COMPUTE_LOG_MEAN_ENERGY) + { + logenergy = log(wave[i*NY+j].mean_energy) + LOG_MEAN_ENERGY_SHIFT; + if (logenergy > LOG_ENERGY_FLOOR) wave[i*NY+j].log_mean_energy = LOG_SHIFT + PLOT_SCALE_LOG_ENERGY*logenergy; + else wave[i*NY+j].mean_energy = LOG_SHIFT + PLOT_SCALE_LOG_ENERGY*LOG_ENERGY_FLOOR; + } + } + else if (first) + { + wave[i*NY+j].energy = 0.0; + wave[i*NY+j].total_energy = 0.0; + wave[i*NY+j].log_total_energy = LOG_ENERGY_FLOOR; + wave[i*NY+j].mean_energy = 0.0; + wave[i*NY+j].log_mean_energy = LOG_ENERGY_FLOOR; + } } + + first = 0; } void compute_log_energy_field(double phi[NX*NY], double psi[NX*NY], short int xy_in[NX*NY], t_wave wave[NX*NY]) /* computes cosine of angle between normal vector and vector light */ +/* TO DO: include into compute_energy_field */ { int i, j; double value; @@ -729,15 +913,16 @@ void compute_log_energy_field(double phi[NX*NY], double psi[NX*NY], short int xy { if ((TWOSPEEDS)||(xy_in[i*NY+j])) { - value = compute_energy_mod(phi, psi, xy_in, i, j); - if (value > 1.0e-100) wave[i*NY+j].log_energy = LOG_SHIFT + PLOT_SCALE_LOG_ENERGY*log(value); - else wave[i*NY+j].log_energy = LOG_SHIFT + log(1.0e-100)*PLOT_SCALE_LOG_ENERGY; + value = log(compute_energy_mod(phi, psi, xy_in, i, j)); + if (value > LOG_ENERGY_FLOOR) wave[i*NY+j].log_energy = LOG_SHIFT + PLOT_SCALE_LOG_ENERGY*value; + else wave[i*NY+j].log_energy = LOG_SHIFT + PLOT_SCALE_LOG_ENERGY*LOG_ENERGY_FLOOR; } else wave[i*NY+j].log_energy = 0.0; } } + void compute_phase_field(double phi[NX*NY], double psi[NX*NY], short int xy_in[NX*NY], t_wave wave[NX*NY]) /* computes cosine of angle between normal vector and vector light */ { @@ -759,7 +944,7 @@ void compute_phase_field(double phi[NX*NY], double psi[NX*NY], short int xy_in[N } } -void compute_light_angle(short int xy_in[NX*NY], t_wave wave[NX*NY]) +void compute_light_angle(short int xy_in[NX*NY], t_wave wave[NX*NY], int movie) /* computes cosine of angle between normal vector and vector light */ { int i, j; @@ -780,8 +965,8 @@ void compute_light_angle(short int xy_in[NX*NY], t_wave wave[NX*NY]) { if ((TWOSPEEDS)||(xy_in[i*NY+j])) { - gradx = (*wave[(i+1)*NY+j].p_zfield - *wave[(i-1)*NY+j].p_zfield)/dx; - grady = (*wave[i*NY+j+1].p_zfield - *wave[i*NY+j-1].p_zfield)/dy; + gradx = (*wave[(i+1)*NY+j].p_zfield[movie] - *wave[(i-1)*NY+j].p_zfield[movie])/dx; + grady = (*wave[i*NY+j+1].p_zfield[movie] - *wave[i*NY+j-1].p_zfield[movie])/dy; norm = sqrt(1.0 + gradx*gradx + grady*grady); pscal = -gradx*light[0] - grady*light[1] + 1.0; @@ -791,25 +976,6 @@ void compute_light_angle(short int xy_in[NX*NY], t_wave wave[NX*NY]) } -// void energy_color_scheme(int palette, double energy, double rgb[3]) -// { -// if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, energy, 1.0, 0, rgb); -// else color_scheme_palette(COLOR_SCHEME, palette, energy, 1.0, 0, rgb); -// } -// -// void log_energy_color_scheme(int palette, double log_energy, double rgb[3]) -// { -// color_scheme_palette(COLOR_SCHEME, palette, LOG_SHIFT + LOG_SCALE*log_energy, 1.0, 0, rgb); -// // if (energy > 1.0e-8) printf("energy = %.3lg, log energy = %.3lg\n", energy, LOG_SHIFT + LOG_SCALE*log(energy)); -// } -// -// void phase_color_scheme(int palette, double phase, double rgb[3]) -// { -// // color_scheme_palette(C_ONEDIM_LINEAR, palette, phase/DPI, 1.0, 0, rgb); -// amp_to_rgb_palette(phase, rgb, palette); -// } - - void compute_field_color(double value, int cplot, int palette, double rgb[3]) /* compute the color depending on the field value and color palette */ { @@ -821,20 +987,39 @@ void compute_field_color(double value, int cplot, int palette, double rgb[3]) } case (P_3D_ENERGY): { -// energy_color_scheme(palette, VSCALE_ENERGY*value, rgb); if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 0, rgb); else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 0, rgb); break; } case (P_3D_LOG_ENERGY): { -// log_energy_color_scheme(palette, value, rgb); + color_scheme_palette(COLOR_SCHEME, palette, LOG_SHIFT + LOG_SCALE*value, 1.0, 0, rgb); + break; + } + case (P_3D_TOTAL_ENERGY): + { + if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 0, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 0, rgb); + break; + } + case (P_3D_LOG_TOTAL_ENERGY): + { + color_scheme_palette(COLOR_SCHEME, palette, LOG_SHIFT + LOG_SCALE*value, 1.0, 0, rgb); + break; + } + case (P_3D_MEAN_ENERGY): + { + if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 0, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 0, rgb); + break; + } + case (P_3D_LOG_MEAN_ENERGY): + { color_scheme_palette(COLOR_SCHEME, palette, LOG_SHIFT + LOG_SCALE*value, 1.0, 0, rgb); break; } case (P_3D_PHASE): { -// phase_color_scheme(palette, value, rgb); amp_to_rgb_palette(value, rgb, palette); break; } @@ -842,25 +1027,25 @@ void compute_field_color(double value, int cplot, int palette, double rgb[3]) } -double compute_interpolated_colors_wave(int i, int j, t_wave wave[NX*NY], double palette, int cplot, +double compute_interpolated_colors_wave(int i, int j, short int xy_in[NX*NY], t_wave wave[NX*NY], double palette, int cplot, double rgb_e[3], double rgb_w[3], double rgb_n[3], double rgb_s[3], - int fade, double fade_value) + int fade, double fade_value, int movie) { int k; double cw, ce, cn, cs, c_sw, c_se, c_nw, c_ne, c_mid, ca, z_mid; double *z_sw, *z_se, *z_nw, *z_ne; - z_sw = wave[i*NY+j].p_zfield; - z_se = wave[(i+1)*NY+j].p_zfield; - z_nw = wave[i*NY+j+1].p_zfield; - z_ne = wave[(i+1)*NY+j+1].p_zfield; + z_sw = wave[i*NY+j].p_zfield[movie]; + z_se = wave[(i+1)*NY+j].p_zfield[movie]; + z_nw = wave[i*NY+j+1].p_zfield[movie]; + z_ne = wave[(i+1)*NY+j+1].p_zfield[movie]; z_mid = 0.25*(*z_sw + *z_se + *z_nw + *z_ne); - c_sw = *wave[i*NY+j].p_cfield; - c_se = *wave[(i+1)*NY+j].p_cfield; - c_nw = *wave[i*NY+j+1].p_cfield; - c_ne = *wave[(i+1)*NY+j+1].p_cfield; + c_sw = *wave[i*NY+j].p_cfield[movie]; + c_se = *wave[(i+1)*NY+j].p_cfield[movie]; + c_nw = *wave[i*NY+j+1].p_cfield[movie]; + c_ne = *wave[(i+1)*NY+j+1].p_cfield[movie]; c_mid = 0.25*(c_sw + c_se + c_nw + c_ne); @@ -878,6 +1063,8 @@ double compute_interpolated_colors_wave(int i, int j, t_wave wave[NX*NY], double { ca = wave[i*NY+j].cos_angle; ca = (ca + 1.0)*0.4 + 0.2; +// if ((FADE_IN_OBSTACLE)&&(!xy_in[i*NY+j])) ca *= 1.6; + if ((FADE_IN_OBSTACLE)&&(!xy_in[i*NY+j])) ca = (ca + 0.1)*1.6; for (k=0; k<3; k++) { rgb_e[k] *= ca; @@ -904,11 +1091,11 @@ void compute_wave_fields(double phi[NX*NY], double psi[NX*NY], short int xy_in[N { int i, j; - if ((zplot == P_3D_ENERGY)||(cplot == P_3D_ENERGY)) + if (COMPUTE_ENERGY) compute_energy_field(phi, psi, xy_in, wave); - if ((zplot == P_3D_LOG_ENERGY)||(cplot == P_3D_LOG_ENERGY)) - compute_log_energy_field(phi, psi, xy_in, wave); +// if ((zplot == P_3D_LOG_ENERGY)||(cplot == P_3D_LOG_ENERGY)) +// compute_log_energy_field(phi, psi, xy_in, wave); if ((zplot == P_3D_PHASE)||(cplot == P_3D_PHASE)) compute_phase_field(phi, psi, xy_in, wave); @@ -916,7 +1103,90 @@ void compute_wave_fields(double phi[NX*NY], double psi[NX*NY], short int xy_in[N } -void compute_zfield(double phi[NX*NY], double psi[NX*NY], short int xy_in[NX*NY], int zplot, t_wave wave[NX*NY]) +void init_speed_dissipation(short int xy_in[NX*NY], double tc[NX*NY], double tcc[NX*NY], double tgamma[NX*NY]) +/* initialise fields for wave speed and dissipation */ +{ + int i, j, k; + double courant2 = COURANT*COURANT, courantb2 = COURANTB*COURANTB; + double u, v, u1, x, y, xy[2], norm2, speed; + + if (VARIABLE_IOR) + { + switch (IOR) { + case (IOR_MANDELBROT): + { + #pragma omp parallel for private(i,j) + for (i=0; i zfloor)&&(*wave[(i+1)*NY+j].p_zfield[movie] > zfloor)&&(*wave[i*NY+j+1].p_zfield[movie] > zfloor)&&(*wave[(i+1)*NY+j+1].p_zfield[movie] > zfloor); + if (draw) { if (AMPLITUDE_HIGH_RES > 0) { - z_mid = compute_interpolated_colors_wave(i, j, wave, palette, cplot, - rgb_e, rgb_w, rgb_n, rgb_s, fade, fade_value); + z_mid = compute_interpolated_colors_wave(i, j, xy_in, wave, palette, cplot, + rgb_e, rgb_w, rgb_n, rgb_s, fade, fade_value, movie); ij_to_xy(i, j, xy_sw); ij_to_xy(i+1, j, xy_se); @@ -1041,17 +1372,17 @@ void draw_wave_3d(double phi[NX*NY], double psi[NX*NY], short int xy_in[NX*NY], glBegin(GL_TRIANGLE_FAN); glColor3f(rgb_w[0], rgb_w[1], rgb_w[2]); draw_vertex_xyz(xy_mid, z_mid); - draw_vertex_xyz(xy_nw, *wave[i*NY+j+1].p_zfield); - draw_vertex_xyz(xy_sw, *wave[i*NY+j].p_zfield); + draw_vertex_xyz(xy_nw, *wave[i*NY+j+1].p_zfield[movie]); + draw_vertex_xyz(xy_sw, *wave[i*NY+j].p_zfield[movie]); glColor3f(rgb_s[0], rgb_s[1], rgb_s[2]); - draw_vertex_xyz(xy_se, *wave[(i+1)*NY+j].p_zfield); + draw_vertex_xyz(xy_se, *wave[(i+1)*NY+j].p_zfield[movie]); glColor3f(rgb_e[0], rgb_e[1], rgb_e[2]); - draw_vertex_xyz(xy_ne, *wave[(i+1)*NY+j+1].p_zfield); + draw_vertex_xyz(xy_ne, *wave[(i+1)*NY+j+1].p_zfield[movie]); glColor3f(rgb_n[0], rgb_n[1], rgb_n[2]); - draw_vertex_xyz(xy_nw, *wave[i*NY+j+1].p_zfield); + draw_vertex_xyz(xy_nw, *wave[i*NY+j+1].p_zfield[movie]); glEnd (); } else /* experimental */ @@ -1059,29 +1390,29 @@ void draw_wave_3d(double phi[NX*NY], double psi[NX*NY], short int xy_in[NX*NY], glColor3f(rgb_w[0], rgb_w[1], rgb_w[2]); glBegin(GL_TRIANGLE_STRIP); draw_vertex_xyz(xy_mid, z_mid); - draw_vertex_xyz(xy_nw, *wave[i*NY+j+1].p_zfield); - draw_vertex_xyz(xy_sw, *wave[i*NY+j].p_zfield); + draw_vertex_xyz(xy_nw, *wave[i*NY+j+1].p_zfield[movie]); + draw_vertex_xyz(xy_sw, *wave[i*NY+j].p_zfield[movie]); glEnd (); glColor3f(rgb_s[0], rgb_s[1], rgb_s[2]); glBegin(GL_TRIANGLE_STRIP); draw_vertex_xyz(xy_mid, z_mid); - draw_vertex_xyz(xy_sw, *wave[i*NY+j].p_zfield); - draw_vertex_xyz(xy_se, *wave[(i+1)*NY+j].p_zfield); + draw_vertex_xyz(xy_sw, *wave[i*NY+j].p_zfield[movie]); + draw_vertex_xyz(xy_se, *wave[(i+1)*NY+j].p_zfield[movie]); glEnd (); glColor3f(rgb_e[0], rgb_e[1], rgb_e[2]); glBegin(GL_TRIANGLE_STRIP); draw_vertex_xyz(xy_mid, z_mid); - draw_vertex_xyz(xy_se, *wave[(i+1)*NY+j].p_zfield); - draw_vertex_xyz(xy_ne, *wave[(i+1)*NY+j+1].p_zfield); + draw_vertex_xyz(xy_se, *wave[(i+1)*NY+j].p_zfield[movie]); + draw_vertex_xyz(xy_ne, *wave[(i+1)*NY+j+1].p_zfield[movie]); glEnd (); glColor3f(rgb_n[0], rgb_n[1], rgb_n[2]); glBegin(GL_TRIANGLE_STRIP); draw_vertex_xyz(xy_mid, z_mid); - draw_vertex_xyz(xy_nw, *wave[i*NY+j+1].p_zfield); - draw_vertex_xyz(xy_ne, *wave[(i+1)*NY+j+1].p_zfield); + draw_vertex_xyz(xy_nw, *wave[i*NY+j+1].p_zfield[movie]); + draw_vertex_xyz(xy_ne, *wave[(i+1)*NY+j+1].p_zfield[movie]); glEnd (); } } @@ -1091,22 +1422,38 @@ void draw_wave_3d(double phi[NX*NY], double psi[NX*NY], short int xy_in[NX*NY], glBegin(GL_TRIANGLE_FAN); ij_to_xy(i, j, xy); - draw_vertex_xyz(xy, *wave[i*NY+j].p_zfield); + draw_vertex_xyz(xy, *wave[i*NY+j].p_zfield[movie]); ij_to_xy(i+1, j, xy); - draw_vertex_xyz(xy, *wave[(i+1)*NY+j].p_zfield); + draw_vertex_xyz(xy, *wave[(i+1)*NY+j].p_zfield[movie]); ij_to_xy(i+1, j+1, xy); - draw_vertex_xyz(xy, *wave[(i+1)*NY+j+1].p_zfield); + draw_vertex_xyz(xy, *wave[(i+1)*NY+j+1].p_zfield[movie]); ij_to_xy(i, j+1, xy); - draw_vertex_xyz(xy, *wave[i*NY+j+1].p_zfield); + draw_vertex_xyz(xy, *wave[i*NY+j+1].p_zfield[movie]); glEnd (); } } + + if ((DRAW_OUTSIDE_GRAY)&&((!xy_in[i*NY+j]))) + { + glColor3f(0.5, 0.5, 0.5); + glBegin(GL_TRIANGLE_FAN); + ij_to_xy(i, j, xy); + draw_vertex_xyz(xy, 0.0); + ij_to_xy(i+1, j, xy); + draw_vertex_xyz(xy, 0.0); + ij_to_xy(i+1, j+1, xy); + draw_vertex_xyz(xy, 0.0); + ij_to_xy(i, j+1, xy); + draw_vertex_xyz(xy, 0.0); + glEnd (); + } } if (DRAW_BILLIARD_FRONT) draw_billiard_3d_front(fade, fade_value); } -void draw_color_scheme_palette_3d(double x1, double y1, double x2, double y2, int plot, double min, double max, int palette) +void draw_color_scheme_palette_3d(double x1, double y1, double x2, double y2, int plot, double min, double max, + int palette, int fade, double fade_value) { int j, k, ij_botleft[2], ij_topright[2], imin, imax, jmin, jmax; double y, dy, dy_e, dy_phase, rgb[3], value, lum, amp; @@ -1168,19 +1515,44 @@ void draw_color_scheme_palette_3d(double x1, double y1, double x2, double y2, in } case (P_3D_LOG_ENERGY): { -// value = LOG_SHIFT + LOG_SCALE*dy_e*(double)(j - jmin)*100.0/E_SCALE; + value = LOG_SCALE*dy_e*(double)(j - jmin)*100.0/E_SCALE; + color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + break; + } + case (P_3D_TOTAL_ENERGY): + { + value = dy_e*(double)(j - 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_3D_LOG_TOTAL_ENERGY): + { + value = LOG_SCALE*dy_e*(double)(j - jmin)*100.0/E_SCALE; + color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + break; + } + case (P_3D_MEAN_ENERGY): + { + value = dy_e*(double)(j - 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_3D_LOG_MEAN_ENERGY): + { value = LOG_SCALE*dy_e*(double)(j - jmin)*100.0/E_SCALE; color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); break; } case (P_3D_PHASE): { -// value = min + 1.0*dy*(double)(j - jmin); value = dy_phase*(double)(j - jmin); color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb); break; } } + if (fade) for (k=0; k<3; k++) rgb[k] *= fade_value; glColor3f(rgb[0], rgb[1], rgb[2]); if (ROTATE_COLOR_SCHEME) { @@ -1199,7 +1571,8 @@ void draw_color_scheme_palette_3d(double x1, double y1, double x2, double y2, in } glEnd (); - glColor3f(1.0, 1.0, 1.0); + if (fade) glColor3f(fade_value, fade_value, fade_value); + else glColor3f(1.0, 1.0, 1.0); glLineWidth(BOUNDARY_WIDTH); draw_rectangle_noscale(x1, y1, x2, y2); } diff --git a/sub_wave_3d_rde.c b/sub_wave_3d_rde.c new file mode 100644 index 0000000..2158c0d --- /dev/null +++ b/sub_wave_3d_rde.c @@ -0,0 +1,1430 @@ +/*********************/ +/* animation part */ +/*********************/ + + + +void init_3d() /* initialisation of window */ +{ + glLineWidth(3); + + glClearColor(0.0, 0.0, 0.0, 1.0); + glClear(GL_COLOR_BUFFER_BIT); + + glOrtho(XMIN, XMAX, YMIN, YMAX , -1.0, 1.0); +} + +void xyz_to_xy(double x, double y, double z, double xy_out[2]) +{ + int i; + double s, t, xinter[3]; + static double n2, m2, d, sm2, sn2, v[3], h[2], plane_ratio = 0.5; + static int first = 1; + + if ((first)&&(REPRESENTATION_3D == REP_PROJ_3D)) + { + m2 = observer[0]*observer[0] + observer[1]*observer[1]; + n2 = m2 + observer[2]*observer[2]; + d = plane_ratio*n2; + sm2 = sqrt(m2); + sn2 = sqrt(n2); + h[0] = observer[1]/sm2; + h[1] = -observer[0]/sm2; + v[0] = -observer[0]*observer[2]/(sn2*sm2); + v[1] = -observer[1]*observer[2]/(sn2*sm2); + v[2] = m2/(sn2*sm2); + first = 0; + printf("h = (%.3lg, %.3lg)\n", h[0], h[1]); + printf("v = (%.3lg, %.3lg, %.3lg)\n", v[0], v[1], v[2]); + } + + switch (REPRESENTATION_3D) { + case (REP_AXO_3D): + { + for (i=0; i<2; i++) + xy_out[i] = x*u_3d[i] + y*v_3d[i] + z*w_3d[i]; + break; + } + case (REP_PROJ_3D): + { + if (z > ZMAX_FACTOR*n2) z = ZMAX_FACTOR*n2; + z *= Z_SCALING_FACTOR; + s = observer[0]*x + observer[1]*y + observer[2]*z; + t = (d - s)/(n2 - s); + xinter[0] = t*observer[0] + (1.0-t)*x; + xinter[1] = t*observer[1] + (1.0-t)*y; + xinter[2] = t*observer[2] + (1.0-t)*z; + + xy_out[0] = XSHIFT_3D + XY_SCALING_FACTOR*(xinter[0]*h[0] + xinter[1]*h[1]); + xy_out[1] = YSHIFT_3D + XY_SCALING_FACTOR*(xinter[0]*v[0] + xinter[1]*v[1] + xinter[2]*v[2]); + break; + } + } +} + + +void draw_vertex_ij(int i, int j) +{ + double xy[2]; + + ij_to_xy(i, j, xy); +// if (xy[1] > 0.0) printf("(i,j) = (%i,%i), (x,y) = (%.2lg,%.2lg)\n", i, j, xy[0], xy[1]); + glVertex2d(xy[0], xy[1]); +} + + +void draw_vertex_xyz(double xy[2], double z) +{ + double xy_screen[2]; + + xyz_to_xy(xy[0], xy[1], z, xy_screen); + glVertex2d(xy_screen[0], xy_screen[1]); +} + +void draw_vertex_x_y_z(double x, double y, double z) +{ + double xy_screen[2]; + + xyz_to_xy(x, y, z, xy_screen); + glVertex2d(xy_screen[0], xy_screen[1]); +} + +void draw_rectangle_3d(double x1, double y1, double x2, double y2) +{ + glBegin(GL_LINE_LOOP); + draw_vertex_x_y_z(x1, y1, 0.0); + draw_vertex_x_y_z(x2, y1, 0.0); + draw_vertex_x_y_z(x2, y2, 0.0); + draw_vertex_x_y_z(x1, y2, 0.0); + glEnd(); +} + + +void draw_rectangle_noscale(double x1, double y1, double x2, double y2) +{ + glBegin(GL_LINE_LOOP); + glVertex2d(x1, y1); + glVertex2d(x2, y1); + glVertex2d(x2, y2); + glVertex2d(x1, y2); + glEnd(); +} + +void draw_circle_3d(double x, double y, double r, int nseg) +{ + int i; + double pos[2], alpha, dalpha; + + dalpha = DPI/(double)nseg; + + glBegin(GL_LINE_LOOP); + for (i=0; i<=nseg; i++) + { + alpha = (double)i*dalpha; + draw_vertex_x_y_z(x + r*cos(alpha), y + r*sin(alpha), 0.0); + } + glEnd(); +} + +void tvertex_lineto_3d(t_vertex z) +/* draws boundary segments of isospectral billiard */ +{ + draw_vertex_x_y_z(z.x, z.y, 0.0); +} + + +void draw_billiard_3d(int fade, double fade_value) /* draws the billiard boundary */ +{ + double x0, x, y, x1, y1, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z, l, width, a, b, c, ymax; + int i, j, k, k1, k2, mr2; + static int first = 1, nsides; + + if (BLACK) + { + if (fade) glColor3f(fade_value, fade_value, fade_value); + else glColor3f(1.0, 1.0, 1.0); + } + else + { + if (fade) glColor3f(1.0 - fade_value, 1.0 - fade_value, 1.0 - fade_value); + else glColor3f(0.0, 0.0, 0.0); + } + glLineWidth(BOUNDARY_WIDTH); + + glEnable(GL_LINE_SMOOTH); + + switch (B_DOMAIN) { + case (D_ELLIPSE): + { + glBegin(GL_LINE_LOOP); + for (i=0; i<=NSEG; i++) + { + phi = (double)i*DPI/(double)NSEG; + x = LAMBDA*cos(phi); + y = sin(phi); + draw_vertex_x_y_z(x, y, 0.0); + } + glEnd (); + + /* draw foci */ +// if (FOCI) +// { +// glColor3f(0.3, 0.3, 0.3); +// x0 = sqrt(LAMBDA*LAMBDA-1.0); +// +// glLineWidth(2); +// glEnable(GL_LINE_SMOOTH); +// +// draw_circle(x0, 0.0, r, NSEG); +// draw_circle(-x0, 0.0, r, NSEG); +// } + break; + } + case (D_STADIUM): + { + glBegin(GL_LINE_LOOP); + for (i=0; i<=NSEG; i++) + { + phi = -PID + (double)i*PI/(double)NSEG; + x = 0.5*LAMBDA + cos(phi); + y = sin(phi); + draw_vertex_x_y_z(x, y, 0.0); + } + for (i=0; i<=NSEG; i++) + { + phi = PID + (double)i*PI/(double)NSEG; + x = -0.5*LAMBDA + cos(phi); + y = sin(phi); + draw_vertex_x_y_z(x, y, 0.0); + } + glEnd(); + break; + } + case (D_SINAI): + { + draw_circle_3d(0.0, 0.0, LAMBDA, NSEG); + break; + } + case (D_YOUNG): + { + if (FILL_BILLIARD_COMPLEMENT) + { + if (fade) glColor3f(0.75*fade_value, 0.75*fade_value, 0.75*fade_value); + else glColor3f(0.75, 0.75, 0.75); + + glBegin(GL_TRIANGLE_FAN); + draw_vertex_x_y_z(-MU, YMIN, 0.0); + draw_vertex_x_y_z(-MU, -LAMBDA-MU, 0.0); + draw_vertex_x_y_z(MU, -LAMBDA-MU, 0.0); + draw_vertex_x_y_z(MU, YMIN, 0.0); + glEnd(); + + glBegin(GL_TRIANGLE_FAN); + draw_vertex_x_y_z(-MU, YMAX, 0.0); + draw_vertex_x_y_z(-MU, LAMBDA+MU, 0.0); + draw_vertex_x_y_z(MU, LAMBDA+MU, 0.0); + draw_vertex_x_y_z(MU, YMAX, 0.0); + glEnd(); + + glBegin(GL_TRIANGLE_FAN); + draw_vertex_x_y_z(-MU, -LAMBDA+MU, 0.0); + draw_vertex_x_y_z(-MU, LAMBDA-MU, 0.0); + draw_vertex_x_y_z(MU, LAMBDA-MU, 0.0); + draw_vertex_x_y_z(MU, -LAMBDA+MU, 0.0); + glEnd(); + } + else + { + glBegin(GL_LINE_STRIP); + draw_vertex_x_y_z(-MU, YMIN, 0.0); + draw_vertex_x_y_z(-MU, -LAMBDA-MU, 0.0); + draw_vertex_x_y_z(MU, -LAMBDA-MU, 0.0); + draw_vertex_x_y_z(MU, YMIN, 0.0); + glEnd(); + + glBegin(GL_LINE_STRIP); + draw_vertex_x_y_z(-MU, YMAX, 0.0); + draw_vertex_x_y_z(-MU, LAMBDA+MU, 0.0); + draw_vertex_x_y_z(MU, LAMBDA+MU, 0.0); + draw_vertex_x_y_z(MU, YMAX, 0.0); + glEnd(); + + glBegin(GL_LINE_LOOP); + draw_vertex_x_y_z(-MU, -LAMBDA+MU, 0.0); + draw_vertex_x_y_z(-MU, LAMBDA-MU, 0.0); + draw_vertex_x_y_z(MU, LAMBDA-MU, 0.0); + draw_vertex_x_y_z(MU, -LAMBDA+MU, 0.0); + glEnd(); + } + } + case (D_GRATING): + { + k1 = -(int)(-YMIN/LAMBDA); + k2 = (int)(YMAX/LAMBDA); + for (i=k1; i<= k2; i++) + { + z = (double)i*LAMBDA; + draw_circle_3d(0.0, z, MU, NSEG); + } + break; + } + case (D_EHRENFEST): + { + alpha = asin(MU/LAMBDA); + x0 = 1.0 - sqrt(LAMBDA*LAMBDA - MU*MU); + dphi = 2.0*(PI-alpha)/((double)NSEG); + glBegin(GL_LINE_LOOP); + for (i=0; i<=NSEG; i++) + { + phi = -PI + alpha + (double)i*dphi; + x = 1.0 + (LAMBDA + 0.01)*cos(phi); + y = (LAMBDA + 0.01)*sin(phi); + draw_vertex_x_y_z(x, y, 0.0); + } + for (i=0; i<=NSEG; i++) + { + phi = alpha + (double)i*dphi; + x = -1.0 + (LAMBDA + 0.01)*cos(phi); + y = (LAMBDA + 0.01)*sin(phi); + draw_vertex_x_y_z(x, y, 0.0); + } + glEnd (); + break; + } + case (D_TWO_PARABOLAS): + { + dy = 3.0*MU/(double)NSEG; + width = 0.25*MU; + if (width > 0.2) width = 0.2; + glBegin(GL_LINE_LOOP); + for (i = 0; i < NSEG+1; i++) + { + y = -1.5*MU + dy*(double)i; + x = 0.25*y*y/MU - MU - LAMBDA; + draw_vertex_x_y_z(x, y, 0.0); + } + for (i = 0; i < NSEG+1; i++) + { + y = 1.5*MU - dy*(double)i; + x = 0.25*y*y/MU - (MU + width) - LAMBDA; + draw_vertex_x_y_z(x, y, 0.0); + } + glEnd (); + + glBegin(GL_LINE_LOOP); + for (i = 0; i < NSEG+1; i++) + { + y = -1.5*MU + dy*(double)i; + x = LAMBDA + MU - 0.25*y*y/MU; + draw_vertex_x_y_z(x, y, 0.0); + } + for (i = 0; i < NSEG+1; i++) + { + y = 1.5*MU - dy*(double)i; + x = LAMBDA + (MU + width) - 0.25*y*y/MU; + draw_vertex_x_y_z(x, y, 0.0); + } + glEnd (); + +// if (FOCI) +// { +// glColor3f(0.3, 0.3, 0.3); +// draw_circle(-LAMBDA, 0.0, r, NSEG); +// draw_circle(LAMBDA, 0.0, r, NSEG); +// } + + break; + } + case (D_POLY_PARABOLAS): + { + omega = PI/((double)NPOLY); + a = 0.25/MU; + b = 1.0/tan(omega); + c = LAMBDA + MU; + ymax = (-b + sqrt(b*b + 4.0*a*c))/(2.0*a); + dy = 2.0*ymax/(double)NSEG; + + glBegin(GL_LINE_LOOP); + for (k=0; k 0) + { + glLineWidth(2); + x = 1.0/((double)MRATIO); + draw_rectangle_3d(x, x, -x, -x); + } + + /* level 2 */ + if (MDEPTH > 1) + { + glLineWidth(1); + mr2 = MRATIO*MRATIO; + l = 2.0/((double)mr2); + + for (i=0; i 2) + { + glLineWidth(1); + l = 2.0/((double)(mr2*MRATIO)); + + for (i=0; i45; i--) tvertex_lineto_3d(polyline[i]); + glEnd(); + + /* inner lines */ +// glLineWidth(BOUNDARY_WIDTH/2); + glLineWidth(1); + glColor3f(0.75, 0.75, 0.75); + glBegin(GL_LINE_STRIP); + tvertex_lineto_3d(polyline[0]); + tvertex_lineto_3d(polyline[1]); + tvertex_lineto_3d(polyline[2]); + tvertex_lineto_3d(polyline[0]); + tvertex_lineto_3d(polyline[3]); + tvertex_lineto_3d(polyline[4]); + glEnd(); + + glBegin(GL_LINE_STRIP); + tvertex_lineto_3d(polyline[0]); + tvertex_lineto_3d(polyline[44]); + tvertex_lineto_3d(polyline[45]); + tvertex_lineto_3d(polyline[0]); + tvertex_lineto_3d(polyline[46]); + tvertex_lineto_3d(polyline[45]); + glEnd(); + + for (i=3; i<43; i++) + { + glBegin(GL_LINE_STRIP); + tvertex_lineto_3d(polyline[i]); + tvertex_lineto_3d(polyline[43]); + glEnd(); + glBegin(GL_LINE_STRIP); + tvertex_lineto_3d(polyline[i+42]); + tvertex_lineto_3d(polyline[85]); + glEnd(); + } + + break; + } + case (D_ISOSPECTRAL): + { + /* 1st triangle */ + glBegin(GL_LINE_LOOP); + tvertex_lineto_3d(polyline[0]); + tvertex_lineto_3d(polyline[4]); + tvertex_lineto_3d(polyline[7]); + tvertex_lineto_3d(polyline[1]); + tvertex_lineto_3d(polyline[5]); + tvertex_lineto_3d(polyline[8]); + tvertex_lineto_3d(polyline[2]); + tvertex_lineto_3d(polyline[3]); + tvertex_lineto_3d(polyline[6]); + glEnd(); + + /* inner lines */ + glBegin(GL_LINE_LOOP); + tvertex_lineto_3d(polyline[0]); + tvertex_lineto_3d(polyline[1]); + tvertex_lineto_3d(polyline[2]); + tvertex_lineto_3d(polyline[0]); + tvertex_lineto_3d(polyline[3]); + tvertex_lineto_3d(polyline[2]); + tvertex_lineto_3d(polyline[5]); + tvertex_lineto_3d(polyline[1]); + tvertex_lineto_3d(polyline[4]); + glEnd(); + + /* 2nd triangle */ + glBegin(GL_LINE_LOOP); + tvertex_lineto_3d( polyline[9]); + tvertex_lineto_3d(polyline[16]); + tvertex_lineto_3d(polyline[13]); + tvertex_lineto_3d(polyline[10]); + tvertex_lineto_3d(polyline[17]); + tvertex_lineto_3d(polyline[14]); + tvertex_lineto_3d(polyline[11]); + tvertex_lineto_3d(polyline[15]); + tvertex_lineto_3d(polyline[12]); + glEnd(); + + /* inner lines */ + glBegin(GL_LINE_LOOP); + tvertex_lineto_3d( polyline[9]); + tvertex_lineto_3d(polyline[10]); + tvertex_lineto_3d(polyline[11]); + tvertex_lineto_3d( polyline[9]); + tvertex_lineto_3d(polyline[13]); + tvertex_lineto_3d(polyline[10]); + tvertex_lineto_3d(polyline[14]); + tvertex_lineto_3d(polyline[11]); + tvertex_lineto_3d(polyline[12]); + glEnd(); + break; + } + case (D_HOMOPHONIC): + { + /* 1st triangle */ + glBegin(GL_LINE_LOOP); + tvertex_lineto_3d(polyline[1]); + tvertex_lineto_3d(polyline[3]); + tvertex_lineto_3d(polyline[4]); + tvertex_lineto_3d(polyline[5]); + tvertex_lineto_3d(polyline[6]); + tvertex_lineto_3d(polyline[8]); + tvertex_lineto_3d(polyline[9]); + tvertex_lineto_3d(polyline[10]); + tvertex_lineto_3d(polyline[12]); + tvertex_lineto_3d(polyline[13]); + tvertex_lineto_3d(polyline[15]); + tvertex_lineto_3d(polyline[16]); + tvertex_lineto_3d(polyline[17]); + tvertex_lineto_3d(polyline[18]); + tvertex_lineto_3d(polyline[20]); + glEnd(); + + /* inner lines */ + glLineWidth(BOUNDARY_WIDTH/2); + glBegin(GL_LINE_STRIP); + tvertex_lineto_3d(polyline[9]); + tvertex_lineto_3d(polyline[1]); + tvertex_lineto_3d(polyline[2]); + tvertex_lineto_3d(polyline[5]); + tvertex_lineto_3d(polyline[7]); + tvertex_lineto_3d(polyline[2]); + tvertex_lineto_3d(polyline[8]); + tvertex_lineto_3d(polyline[21]); + tvertex_lineto_3d(polyline[10]); + tvertex_lineto_3d(polyline[2]); + tvertex_lineto_3d(polyline[21]); + tvertex_lineto_3d(polyline[11]); + tvertex_lineto_3d(polyline[13]); + tvertex_lineto_3d(polyline[21]); + tvertex_lineto_3d(polyline[14]); + tvertex_lineto_3d(polyline[20]); + tvertex_lineto_3d(polyline[15]); + tvertex_lineto_3d(polyline[19]); + tvertex_lineto_3d(polyline[16]); + tvertex_lineto_3d(polyline[18]); + glEnd(); + + /* 2nd triangle */ + glLineWidth(BOUNDARY_WIDTH); + glBegin(GL_LINE_LOOP); + tvertex_lineto_3d(polyline[22+10]); + tvertex_lineto_3d(polyline[22+16]); + tvertex_lineto_3d(polyline[22+17]); + tvertex_lineto_3d(polyline[22+18]); + tvertex_lineto_3d(polyline[22+12]); + tvertex_lineto_3d(polyline[22+13]); + tvertex_lineto_3d(polyline[22+15]); + tvertex_lineto_3d(polyline[22+19]); + tvertex_lineto_3d(polyline[22+20]); + tvertex_lineto_3d(polyline[22+1]); + tvertex_lineto_3d(polyline[22+4]); + tvertex_lineto_3d(polyline[22+5]); + tvertex_lineto_3d(polyline[22+7]); + tvertex_lineto_3d(polyline[22+8]); + tvertex_lineto_3d(polyline[22+9]); + glEnd(); + + /* inner lines */ + glLineWidth(BOUNDARY_WIDTH/2); + glBegin(GL_LINE_STRIP); + tvertex_lineto_3d(polyline[22+2]); + tvertex_lineto_3d(polyline[22+6]); + tvertex_lineto_3d(polyline[22+8]); + tvertex_lineto_3d(polyline[22+2]); + tvertex_lineto_3d(polyline[22+5]); + tvertex_lineto_3d(polyline[22+3]); + tvertex_lineto_3d(polyline[22+2]); + tvertex_lineto_3d(polyline[22+1]); + tvertex_lineto_3d(polyline[22+0]); + tvertex_lineto_3d(polyline[22+21]); + tvertex_lineto_3d(polyline[22+18]); + tvertex_lineto_3d(polyline[22+16]); + tvertex_lineto_3d(polyline[22+13]); + tvertex_lineto_3d(polyline[22+21]); + tvertex_lineto_3d(polyline[22+10]); + tvertex_lineto_3d(polyline[22+12]); + tvertex_lineto_3d(polyline[22+21]); + tvertex_lineto_3d(polyline[22+14]); + tvertex_lineto_3d(polyline[22+20]); + tvertex_lineto_3d(polyline[22+15]); + glEnd(); + break; + } + default: + { + break; + } + } +} + +void draw_polyline_visible(int j, int k, double margin) +/* hack to draw the billiard boundary in front of the wave */ +/* only parts of the boundary having a small enough angle with the observer vector are drawn */ +{ + double x, y, x1, y1, length, length1; + static int first = 1; + static double olength; + + if (first) + { + olength = module2(observer[0], observer[1]); + first = 0; + } + + x = polyline[j].x; + y = polyline[j].y; + x1 = polyline[k].x; + y1 = polyline[k].y; + length = module2(x,y); + length1 = module2(x1,y1); + if ((x*observer[0] + y*observer[1] > margin*length*olength)&&(x1*observer[0] + y1*observer[1] > margin*length1*olength)) + { + glBegin(GL_LINE_STRIP); + tvertex_lineto_3d(polyline[j]); + tvertex_lineto_3d(polyline[k]); + glEnd(); + } +} + +void draw_vertex_visible(double x, double y, double margin) +/* hack to draw the billiard boundary in front of the wave */ +/* only parts of the boundary having a small enough angle with the observer vector are drawn */ +{ + static int first = 1; + static double olength; + + if (first) + { + olength = module2(observer[0], observer[1]); + first = 0; + } + + if (x*observer[0] + y*observer[1] > margin*olength*module2(x,y)) + draw_vertex_x_y_z(x, y, 0.0); + else + { + glEnd(); + glBegin(GL_LINE_STRIP); + } +} + +void draw_billiard_3d_front(int fade, double fade_value) +/* hack to draw the billiard boundary in front of the wave */ +/* only parts of the boundary having a small enough angle with the observer vector are drawn */ +{ + double x0, x, y, x1, y1, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z, l, width, a, b, c, ymax, length, length1; + int i, j, k, k1, k2, mr2; + static int first = 1, nsides; + static double olength; + + if (first) + { + olength = module2(observer[0], observer[1]); + first = 0; + } + + if (BLACK) + { + if (fade) glColor3f(fade_value, fade_value, fade_value); + else glColor3f(1.0, 1.0, 1.0); + } + else + { + if (fade) glColor3f(1.0 - fade_value, 1.0 - fade_value, 1.0 - fade_value); + else glColor3f(0.0, 0.0, 0.0); + } + glLineWidth(BOUNDARY_WIDTH); + + glEnable(GL_LINE_SMOOTH); + + switch (B_DOMAIN) { + case (D_STADIUM): + { + glBegin(GL_LINE_STRIP); + for (i=0; i<=NSEG; i++) + { + phi = -PID + (double)i*PI/(double)NSEG; + x = 0.5*LAMBDA + cos(phi); + y = sin(phi); + draw_vertex_visible(x, y, 0.0); + } + for (i=0; i<=NSEG; i++) + { + phi = PID + (double)i*PI/(double)NSEG; + x = -0.5*LAMBDA + cos(phi); + y = sin(phi); + draw_vertex_visible(x, y, 0.0); + } + glEnd(); + break; + } + case (D_YOUNG): + { + glBegin(GL_LINE_STRIP); + draw_vertex_x_y_z(-MU, YMIN, 0.0); + draw_vertex_x_y_z(-MU, -LAMBDA-MU, 0.0); + glEnd(); + + glBegin(GL_LINE_STRIP); + draw_vertex_x_y_z(-MU, YMAX, 0.0); + draw_vertex_x_y_z(-MU, LAMBDA+MU, 0.0); + glEnd(); + + glBegin(GL_LINE_LOOP); + draw_vertex_x_y_z(-MU, -LAMBDA+MU, 0.0); + draw_vertex_x_y_z(-MU, LAMBDA-MU, 0.0); + glEnd(); + break; + } + case (D_TOKA_PRIME): + { + draw_polyline_visible(0, 4, 0.2); + for (i=4; i<41; i++) draw_polyline_visible(i, i+1, -0.1); +// draw_polyline_visible(42, 3, 0.2); + for (i=84; i>46; i--) draw_polyline_visible(i, i-1, -0.1); + draw_polyline_visible(46, 0, 0.2); + break; + } + case (D_STAR): + { + glLineWidth(BOUNDARY_WIDTH); + for (i=0; i= 1.0) angle -= 1.0; + else if (angle < 0.0) angle += 1.0; + phase[i*NY+j] = angle; + } + else phase[i*NY+j] = 0.0; + } +} + + +// void compute_log_energy_field(double *phi[NX], double *psi[NX], short int xy_in[NX*NY], double *energies[NX]) +// /* computes cosine of angle between normal vector and vector light */ +// { +// int i, j; +// +// for (i=0; i= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, energy, 1.0, 0, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, energy, 1.0, 0, rgb); +} + +void log_energy_color_scheme(int palette, double energy, double rgb[3]) +{ + color_scheme_palette(COLOR_SCHEME, palette, LOG_SHIFT + LOG_SCALE*log(energy), 1.0, 0, rgb); +// if (energy > 1.0e-8) printf("energy = %.3lg, log energy = %.3lg\n", energy, LOG_SHIFT + LOG_SCALE*log(energy)); +} + +void phase_color_scheme(int palette, double phase, double rgb[3]) +{ +// color_scheme_palette(C_ONEDIM_LINEAR, palette, phase/DPI, 1.0, 0, rgb); + amp_to_rgb_palette(phase, rgb, palette); +} + + +void compute_interpolated_colors(int i, int j, double field[NX*NY], double palette, int plot, + double *z_sw, double *z_se, double *z_nw, double *z_ne, double *z_mid, + double rgb_e[3], double rgb_w[3], double rgb_n[3], double rgb_s[3]) +{ + double zw, ze, zn, zs; + + *z_sw = field[i*NY+j]; + *z_se = field[(i+1)*NY+j]; + *z_nw = field[i*NY+j+1]; + *z_ne = field[(i+1)*NY+j+1]; + + *z_mid = 0.25*(*z_sw + *z_se + *z_nw + *z_ne); + + zw = (*z_sw + *z_nw + *z_mid)/3.0; + ze = (*z_se + *z_ne + *z_mid)/3.0; + zs = (*z_sw + *z_se + *z_mid)/3.0; + zn = (*z_nw + *z_ne + *z_mid)/3.0; + + if (plot == P_3D_ENERGY) + { + energy_color_scheme(palette, VSCALE_ENERGY*ze, rgb_e); + energy_color_scheme(palette, VSCALE_ENERGY*zw, rgb_w); + energy_color_scheme(palette, VSCALE_ENERGY*zn, rgb_n); + energy_color_scheme(palette, VSCALE_ENERGY*zs, rgb_s); + } + else if (plot == P_3D_LOG_ENERGY) + { + log_energy_color_scheme(palette, ze, rgb_e); + log_energy_color_scheme(palette, zw, rgb_w); + log_energy_color_scheme(palette, zn, rgb_n); + log_energy_color_scheme(palette, zs, rgb_s); + } + else + { + color_scheme_palette(COLOR_SCHEME, palette, VSCALE_AMPLITUDE*ze, 1.0, 0, rgb_e); + color_scheme_palette(COLOR_SCHEME, palette, VSCALE_AMPLITUDE*zw, 1.0, 0, rgb_w); + color_scheme_palette(COLOR_SCHEME, palette, VSCALE_AMPLITUDE*zn, 1.0, 0, rgb_n); + color_scheme_palette(COLOR_SCHEME, palette, VSCALE_AMPLITUDE*zs, 1.0, 0, rgb_s); + } +} + +void compute_interpolated_colors_new(int i, int j, double zfield[NX*NY], double colorfield[NX*NY], double palette, int cplot, + double *z_sw, double *z_se, double *z_nw, double *z_ne, double *z_mid, + double rgb_e[3], double rgb_w[3], double rgb_n[3], double rgb_s[3]) +{ + double zw, ze, zn, zs, c_sw, c_se, c_nw, c_ne, c_mid; + + *z_sw = zfield[i*NY+j]; + *z_se = zfield[(i+1)*NY+j]; + *z_nw = zfield[i*NY+j+1]; + *z_ne = zfield[(i+1)*NY+j+1]; + + *z_mid = 0.25*(*z_sw + *z_se + *z_nw + *z_ne); + + c_sw = colorfield[i*NY+j]; + c_se = colorfield[(i+1)*NY+j]; + c_nw = colorfield[i*NY+j+1]; + c_ne = colorfield[(i+1)*NY+j+1]; + + c_mid = 0.25*(c_sw + c_se + c_nw + c_ne); + + zw = (c_sw + c_nw + c_mid)/3.0; + ze = (c_se + c_ne + c_mid)/3.0; + zs = (c_sw + c_se + c_mid)/3.0; + zn = (c_nw + c_ne + c_mid)/3.0; + + switch (cplot){ + case (P_3D_ENERGY): + { + energy_color_scheme(palette, VSCALE_ENERGY*ze, rgb_e); + energy_color_scheme(palette, VSCALE_ENERGY*zw, rgb_w); + energy_color_scheme(palette, VSCALE_ENERGY*zn, rgb_n); + energy_color_scheme(palette, VSCALE_ENERGY*zs, rgb_s); + break; + } + case (P_3D_LOG_ENERGY): + { + log_energy_color_scheme(palette, ze, rgb_e); + log_energy_color_scheme(palette, zw, rgb_w); + log_energy_color_scheme(palette, zn, rgb_n); + log_energy_color_scheme(palette, zs, rgb_s); + break; + } + case (P_3D_PHASE): + { + phase_color_scheme(palette, ze, rgb_e); + phase_color_scheme(palette, zw, rgb_w); + phase_color_scheme(palette, zn, rgb_n); + phase_color_scheme(palette, zs, rgb_s); + break; + } + default: + { +// hsl_to_rgb_palette(VSCALE_AMPLITUDE*ze, 0.9, 0.5, rgb_e, palette); +// hsl_to_rgb_palette(VSCALE_AMPLITUDE*zw, 0.9, 0.5, rgb_w, palette); +// hsl_to_rgb_palette(VSCALE_AMPLITUDE*zn, 0.9, 0.5, rgb_n, palette); +// hsl_to_rgb_palette(VSCALE_AMPLITUDE*zs, 0.9, 0.5, rgb_s, palette); + color_scheme_palette(COLOR_SCHEME, palette, VSCALE_AMPLITUDE*ze, 1.0, 0, rgb_e); + color_scheme_palette(COLOR_SCHEME, palette, VSCALE_AMPLITUDE*zw, 1.0, 0, rgb_w); + color_scheme_palette(COLOR_SCHEME, palette, VSCALE_AMPLITUDE*zn, 1.0, 0, rgb_n); + color_scheme_palette(COLOR_SCHEME, palette, VSCALE_AMPLITUDE*zs, 1.0, 0, rgb_s); + } + } +} + +void compute_interpolated_colors_rde_temp(int i, int j, double zfield[NX*NY], double colorfield[NX*NY], double palette, int cplot, + double *z_sw, double *z_se, double *z_nw, double *z_ne, double *z_mid, + double rgb_e[3], double rgb_w[3], double rgb_n[3], double rgb_s[3]) +{ + double zw, ze, zn, zs, c_sw, c_se, c_nw, c_ne, c_mid; + + *z_sw = zfield[i*NY+j]; + *z_se = zfield[(i+1)*NY+j]; + *z_nw = zfield[i*NY+j+1]; + *z_ne = zfield[(i+1)*NY+j+1]; + + *z_mid = 0.25*(*z_sw + *z_se + *z_nw + *z_ne); + + c_sw = colorfield[i*NY+j]; + c_se = colorfield[(i+1)*NY+j]; + c_nw = colorfield[i*NY+j+1]; + c_ne = colorfield[(i+1)*NY+j+1]; + + c_mid = 0.25*(c_sw + c_se + c_nw + c_ne); + + zw = (c_sw + c_nw + c_mid)/3.0; + ze = (c_se + c_ne + c_mid)/3.0; + zs = (c_sw + c_se + c_mid)/3.0; + zn = (c_nw + c_ne + c_mid)/3.0; + + switch (cplot){ + case (P_3D_ENERGY): + { + energy_color_scheme(palette, VSCALE_ENERGY*ze, rgb_e); + energy_color_scheme(palette, VSCALE_ENERGY*zw, rgb_w); + energy_color_scheme(palette, VSCALE_ENERGY*zn, rgb_n); + energy_color_scheme(palette, VSCALE_ENERGY*zs, rgb_s); + break; + } + case (P_3D_LOG_ENERGY): + { + log_energy_color_scheme(palette, ze, rgb_e); + log_energy_color_scheme(palette, zw, rgb_w); + log_energy_color_scheme(palette, zn, rgb_n); + log_energy_color_scheme(palette, zs, rgb_s); + break; + } + case (P_3D_PHASE): + { + phase_color_scheme(palette, ze, rgb_e); + phase_color_scheme(palette, zw, rgb_w); + phase_color_scheme(palette, zn, rgb_n); + phase_color_scheme(palette, zs, rgb_s); + break; + } + default: + { + hsl_to_rgb_palette(VSCALE_AMPLITUDE*ze, 0.9, 0.5, rgb_e, palette); + hsl_to_rgb_palette(VSCALE_AMPLITUDE*zw, 0.9, 0.5, rgb_w, palette); + hsl_to_rgb_palette(VSCALE_AMPLITUDE*zn, 0.9, 0.5, rgb_n, palette); + hsl_to_rgb_palette(VSCALE_AMPLITUDE*zs, 0.9, 0.5, rgb_s, palette); +// color_scheme_palette(COLOR_SCHEME, palette, VSCALE_AMPLITUDE*ze, 1.0, 0, rgb_e); +// color_scheme_palette(COLOR_SCHEME, palette, VSCALE_AMPLITUDE*zw, 1.0, 0, rgb_w); +// color_scheme_palette(COLOR_SCHEME, palette, VSCALE_AMPLITUDE*zn, 1.0, 0, rgb_n); +// color_scheme_palette(COLOR_SCHEME, palette, VSCALE_AMPLITUDE*zs, 1.0, 0, rgb_s); + } + } +} + +void draw_wave_3d(double phi[NX*NY], double psi[NX*NY], short int xy_in[NX*NY], int plot, int cplot, int palette, int fade, double fade_value) +{ + int i, j, k, l, draw = 1; + double xy[2], xy_screen[2], rgb[3], pos[2], ca, rgb_e[3], rgb_w[3], rgb_n[3], rgb_s[3]; + double z_sw, z_se, z_nw, z_ne, z_mid, zw, ze, zn, zs, min = 1000.0, max = 0.0; + double xy_sw[2], xy_se[2], xy_nw[2], xy_ne[2], xy_mid[2]; + double energy; + double *cos_angle, *energies, *phases; + + blank(); + draw_billiard_3d(fade, fade_value); + cos_angle = (double *)malloc(NX*NY*sizeof(double)); + energies = (double *)malloc(NX*NY*sizeof(double)); + phases = (double *)malloc(NX*NY*sizeof(double)); + + + if ((plot == P_3D_ENERGY)||(plot == P_3D_LOG_ENERGY)) +// if (plot == P_3D_ENERGY) + { + compute_energy_field(phi, psi, xy_in, energies); + compute_light_angle(energies, xy_in, cos_angle); + } +// else if (plot == P_3D_LOG_ENERGY) +// { +// compute_log_energy_field(phi, psi, xy_in, energies); +// compute_light_angle(energies, xy_in, cos_angle); +// } + else if (plot != P_3D_AMPLITUDE) compute_light_angle(phi, xy_in, cos_angle); + + if (cplot == P_3D_PHASE) + { + compute_phase_field(phi, psi, xy_in, phases); +// compute_energy_field(phi, psi, xy_in, energies); + } + +// if ((plot == P_3D_ANGLE)||(plot == P_3D_AMP_ANGLE)) compute_light_angle(phi, xy_in, cos_angle); + + for (i=0; i= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + break; + } + case (P_3D_LOG_ENERGY): + { +// value = LOG_SHIFT + LOG_SCALE*dy_e*(double)(j - jmin)*100.0/E_SCALE; + value = LOG_SCALE*dy_e*(double)(j - jmin)*100.0/E_SCALE; + color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + break; + } + case (P_3D_PHASE): + { + value = dy_phase*(double)(j - jmin); + color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb); + break; + } + case (Z_MODULE): + { + value = min + 1.0*dy*(double)(j - jmin); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + case (Z_ARGUMENT): + { + value = dy_phase*(double)(j - jmin); + color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb); + break; + } + } + if (fade) for (k=0; k<3; k++) rgb[k] *= fade_value; + glColor3f(rgb[0], rgb[1], rgb[2]); + if (ROTATE_COLOR_SCHEME) + { + draw_vertex_ij(j, imin); + draw_vertex_ij(j, imax); + draw_vertex_ij(j+1, imax); + draw_vertex_ij(j+1, imin); + } + else + { + draw_vertex_ij(imin, j); + draw_vertex_ij(imax, j); + draw_vertex_ij(imax, j+1); + draw_vertex_ij(imin, j+1); + } + } + glEnd (); + + if (fade) glColor3f(fade_value, fade_value, fade_value); + else glColor3f(1.0, 1.0, 1.0); + glLineWidth(BOUNDARY_WIDTH); + draw_rectangle_noscale(x1, y1, x2, y2); +} + diff --git a/wave_3d.c b/wave_3d.c index 17b9a07..f6a7057 100644 --- a/wave_3d.c +++ b/wave_3d.c @@ -41,34 +41,38 @@ #include #include /* Sam Leffler's libtiff library. */ #include +#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 */ /* General geometrical parameters */ -/* uncomment for higher resolution */ // #define WINWIDTH 1920 /* window width */ // #define WINHEIGHT 1000 /* window height */ -// #define NX 1920 /* number of grid points on x axis */ -// #define NY 1000 /* number of grid points on y axis */ -// // #define NX 3840 /* number of grid points on x axis */ -// // #define NY 2000 /* number of grid points on y axis */ +// // #define NX 900 /* number of grid points on x axis */ +// // #define NY 450 /* number of grid points on y axis */ +// #define NX 1800 /* number of grid points on x axis */ +// #define NY 900 /* number of grid points on y axis */ +// // #define NX 1920 /* number of grid points on x axis */ +// // #define NY 1000 /* number of grid points on y axis */ // -// #define XMIN -2.0 -// #define XMAX 2.0 /* x interval */ +// // #define XMIN -2.0 +// // #define XMAX 2.0 /* x interval */ // #define YMIN -1.041666667 // #define YMAX 1.041666667 /* y interval for 9/16 aspect ratio */ #define HIGHRES 0 /* set to 1 if resolution of grid is double that of displayed image */ -/* comment out for higher resolution */ #define WINWIDTH 1280 /* window width */ #define WINHEIGHT 720 /* window height */ #define NX 1280 /* number of grid points on x axis */ #define NY 720 /* number of grid points on y axis */ +// #define NX 640 /* number of grid points on x axis */ +// #define NY 360 /* number of grid points on y axis */ + #define XMIN -2.0 #define XMAX 2.0 /* x interval */ #define YMIN -1.125 @@ -78,25 +82,29 @@ /* Choice of the billiard table */ -#define B_DOMAIN 16 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN 44 /* choice of domain shape, see list in global_pdes.c */ -#define CIRCLE_PATTERN 201 /* pattern of circles or polygons, see list in global_pdes.c */ +#define CIRCLE_PATTERN 202 /* pattern of circles or polygons, see list in global_pdes.c */ + +#define VARIABLE_IOR 0 /* set to 1 for a variable index of refraction */ +#define IOR 1 /* choice of index of refraction, see list in global_pdes.c */ +#define MANDEL_IOR_SCALE -0.05 /* parameter controlling dependece of IoR on Mandelbrot escape speed */ #define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */ #define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */ -#define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */ +#define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */ -#define LAMBDA 0.6 /* parameter controlling the dimensions of domain */ -#define MU 0.6 /* parameter controlling the dimensions of domain */ -#define NPOLY 6 /* number of sides of polygon */ -#define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */ +#define LAMBDA 0.2 /* parameter controlling the dimensions of domain */ +#define MU 0.5 /* parameter controlling the dimensions of domain */ +#define NPOLY 3 /* number of sides of polygon */ +#define APOLY 2.0 /* angle by which to turn polygon, in units of Pi/2 */ #define MDEPTH 3 /* depth of computation of Menger gasket */ #define MRATIO 3 /* ratio defining Menger gasket */ #define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ #define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */ #define FOCI 1 /* set to 1 to draw focal points of ellipse */ #define NGRIDX 36 /* number of grid point for grid of disks */ -#define NGRIDY 6 /* number of grid point for grid of disks */ +#define NGRIDY 6 /* number of grid point for grid of disks */ #define X_SHOOTER -0.2 #define Y_SHOOTER -0.6 @@ -115,18 +123,16 @@ /* Physical parameters of wave equation */ -// #define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */ -#define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */ +#define TWOSPEEDS 1 /* set to 1 to replace hardcore boundary by medium with different speed */ #define OSCILLATE_LEFT 0 /* set to 1 to add oscilating boundary condition on the left */ #define OSCILLATE_TOPBOT 0 /* set to 1 to enforce a planar wave on top and bottom boundary */ -#define OMEGA 0.005 /* frequency of periodic excitation */ -#define AMPLITUDE 0.8 /* amplitude of periodic excitation */ -#define COURANT 0.06 /* Courant number */ -#define COURANTB 0.03 /* Courant number in medium B */ -// #define COURANTB 0.016363636 /* Courant number in medium B */ +#define OMEGA 0.017 /* frequency of periodic excitation */ +#define AMPLITUDE 0.9 /* amplitude of periodic excitation */ +#define COURANT 0.1 /* Courant number */ +#define COURANTB 0.05 /* Courant number in medium B */ #define GAMMA 0.0 /* damping factor in wave equation */ -#define GAMMAB 1.0e-7 /* damping factor in wave equation */ +#define GAMMAB 2.0e-5 /* damping factor in wave equation */ #define GAMMA_SIDES 1.0e-4 /* damping factor on boundary */ #define GAMMA_TOPBOT 1.0e-7 /* damping factor on boundary */ #define KAPPA 0.0 /* "elasticity" term enforcing oscillations */ @@ -139,53 +145,54 @@ #define ADD_OSCILLATING_SOURCE 0 /* set to 1 to add an oscillating wave source */ #define OSCILLATING_SOURCE_PERIOD 30 /* period of oscillating source */ -// #define OSCILLATING_SOURCE_PERIOD 14 /* period of oscillating source */ /* Boundary conditions, see list in global_pdes.c */ #define B_COND 2 -// #define B_COND 2 /* Parameters for length and speed of simulation */ -#define NSTEPS 2500 /* number of frames of movie */ -#define NVID 10 /* number of iterations between images displayed on screen */ +#define NSTEPS 1500 /* number of frames of movie */ +#define NVID 8 /* number of iterations between images displayed on screen */ #define NSEG 1000 /* number of segments of boundary */ #define INITIAL_TIME 0 /* time after which to start saving frames */ -#define BOUNDARY_WIDTH 3 /* width of billiard boundary */ +#define BOUNDARY_WIDTH 1 /* width of billiard boundary */ -#define PAUSE 200 /* number of frames after which to pause */ +#define PAUSE 100 /* number of frames after which to pause */ #define PSLEEP 2 /* sleep time during pause */ #define SLEEP1 1 /* initial sleeping time */ #define SLEEP2 1 /* final sleeping time */ -#define MID_FRAMES 200 /* number of still frames between parts of two-part movie */ +#define MID_FRAMES 100 /* number of still frames between parts of two-part movie */ #define END_FRAMES 100 /* number of still frames at end of movie */ #define FADE 1 /* set to 1 to fade at end of movie */ /* Parameters of initial condition */ -#define INITIAL_AMP 0.5 /* amplitude of initial condition */ -#define INITIAL_VARIANCE 0.0005 /* variance of initial condition */ -#define INITIAL_WAVELENGTH 0.1 /* wavelength of initial condition */ +#define INITIAL_AMP 0.75 /* amplitude of initial condition */ +#define INITIAL_VARIANCE 0.0003 /* variance of initial condition */ +#define INITIAL_WAVELENGTH 0.02 /* wavelength of initial condition */ /* Plot type, see list in global_pdes.c */ -#define ZPLOT 103 /* wave height */ -#define CPLOT 103 /* color scheme */ +#define ZPLOT 108 /* wave height */ +#define CPLOT 108 /* color scheme */ -#define ZPLOT_B 104 -#define CPLOT_B 104 /* plot type for second movie */ +#define ZPLOT_B 109 +#define CPLOT_B 109 /* plot type for second movie */ #define AMPLITUDE_HIGH_RES 1 /* set to 1 to increase resolution of plot */ #define SHADE_3D 1 /* set to 1 to change luminosity according to normal vector */ #define NON_DIRICHLET_BC 0 /* set to 1 to draw only facets in domain, if field is not zero on boundary */ +#define FLOOR_ZCOORD 1 /* set to 1 to draw only facets with z not too negative */ #define DRAW_BILLIARD 1 /* set to 1 to draw boundary */ -#define DRAW_BILLIARD_FRONT 1 /* set to 1 to draw front of boundary after drawing wave */ +#define DRAW_BILLIARD_FRONT 0 /* set to 1 to draw front of boundary after drawing wave */ +#define DRAW_CONSTRUCTION_LINES 0 /* set to 1 to draw construction lines of certain domains */ #define FADE_IN_OBSTACLE 1 /* set to 1 to fade color inside obstacles */ +#define DRAW_OUTSIDE_GRAY 0 /* experimental, draw outside of billiard in gray */ -#define PLOT_SCALE_ENERGY 0.05 /* vertical scaling in energy plot */ -#define PLOT_SCALE_LOG_ENERGY 0.6 /* vertical scaling in log energy plot */ +#define PLOT_SCALE_ENERGY 0.02 /* vertical scaling in energy plot */ +#define PLOT_SCALE_LOG_ENERGY 1.0 /* vertical scaling in log energy plot */ /* 3D representation */ @@ -194,10 +201,9 @@ #define REP_AXO_3D 0 /* linear projection (axonometry) */ #define REP_PROJ_3D 1 /* projection on plane orthogonal to observer line of sight */ - /* Color schemes */ -#define COLOR_PALETTE 14 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE 14 /* Color palette, see list in global_pdes.c */ #define COLOR_PALETTE_B 11 /* Color palette, see list in global_pdes.c */ #define BLACK 1 /* background */ @@ -205,15 +211,17 @@ #define COLOR_SCHEME 3 /* choice of color scheme, see list in global_pdes.c */ #define SCALE 0 /* set to 1 to adjust color scheme to variance of field */ -#define SLOPE 1.0 /* sensitivity of color on wave amplitude */ -#define VSCALE_AMPLITUDE 0.2 /* additional scaling factor for color scheme P_3D_AMPLITUDE */ -#define VSCALE_ENERGY 0.35 /* additional scaling factor for color scheme P_3D_ENERGY */ +#define SLOPE 1.0 /* sensitivity of color on wave amplitude */ +#define VSCALE_AMPLITUDE 1.0 /* additional scaling factor for color scheme P_3D_AMPLITUDE */ +#define VSCALE_ENERGY 10.0 /* additional scaling factor for color scheme P_3D_ENERGY */ #define PHASE_FACTOR 20.0 /* factor in computation of phase in color scheme P_3D_PHASE */ #define PHASE_SHIFT 0.0 /* shift of phase in color scheme P_3D_PHASE */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ -#define E_SCALE 200.0 /* scaling factor for energy representation */ -#define LOG_SCALE 1.0 /* scaling factor for energy log representation */ -#define LOG_SHIFT 1.0 /* shift of colors on log scale */ +#define E_SCALE 1000.0 /* scaling factor for energy representation */ +#define LOG_SCALE 0.25 /* scaling factor for energy log representation */ +#define LOG_SHIFT 0.0 /* shift of colors on log scale */ +#define LOG_ENERGY_FLOOR -1000.0 /* floor value for log of (total) energy */ +#define LOG_MEAN_ENERGY_SHIFT 5.0 /* additional shift for log of mean energy */ #define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */ #define COLORHUE 260 /* initial hue of water color for scheme C_LUM */ @@ -223,9 +231,9 @@ #define HUEMEAN 240.0 /* mean value of hue for color scheme C_HUE */ #define HUEAMP -200.0 /* amplitude of variation of hue for color scheme C_HUE */ -#define DRAW_COLOR_SCHEME 0 /* set to 1 to plot the color scheme */ -#define COLORBAR_RANGE 3.0 /* scale of color scheme bar */ -#define COLORBAR_RANGE_B 5.0 /* scale of color scheme bar for 2nd part */ +#define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */ +#define COLORBAR_RANGE 20.0 /* scale of color scheme bar */ +#define COLORBAR_RANGE_B 100.0 /* scale of color scheme bar for 2nd part */ #define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ #define SAVE_TIME_SERIES 0 /* set to 1 to save wave time series at a point */ @@ -240,13 +248,13 @@ double u_3d[2] = {0.75, -0.45}; /* projections of basis vectors for REP_AXO_ double v_3d[2] = {-0.75, -0.45}; double w_3d[2] = {0.0, 0.015}; double light[3] = {0.816496581, -0.40824829, 0.40824829}; /* vector of "light" direction for P_3D_ANGLE color scheme */ -double observer[3] = {10.0, 6.0, 8.5}; /* location of observer for REP_PROJ_3D representation */ +double observer[3] = {6.0, 8.0, 6.0}; /* location of observer for REP_PROJ_3D representation */ -#define Z_SCALING_FACTOR 0.018 /* overall scaling factor of z axis for REP_PROJ_3D representation */ -#define XY_SCALING_FACTOR 3.75 /* overall scaling factor for on-screen (x,y) coordinates after projection */ +#define Z_SCALING_FACTOR 0.05 /* overall scaling factor of z axis for REP_PROJ_3D representation */ +#define XY_SCALING_FACTOR 2.0 /* overall scaling factor for on-screen (x,y) coordinates after projection */ #define ZMAX_FACTOR 1.0 /* max value of z coordinate for REP_PROJ_3D representation */ -#define XSHIFT_3D 0.0 /* overall x shift for REP_PROJ_3D representation */ -#define YSHIFT_3D 0.0 /* overall y shift for REP_PROJ_3D representation */ +#define XSHIFT_3D -0.1 /* overall x shift for REP_PROJ_3D representation */ +#define YSHIFT_3D 0.15 /* overall y shift for REP_PROJ_3D representation */ #include "global_pdes.c" /* constants and global variables */ @@ -500,10 +508,15 @@ void evolve_wave(double phi[NX*NY], double psi[NX*NY], double phi_tmp[NX*NY], do } -void draw_color_bar_palette(int plot, double range, int palette) +void draw_color_bar_palette(int plot, double range, int palette, int fade, double fade_value) { - if (ROTATE_COLOR_SCHEME) draw_color_scheme_palette_3d(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range, palette); - else draw_color_scheme_palette_3d(XMAX - 0.3, YMIN + 0.1, XMAX - 0.1, YMAX - 0.1, plot, -range, range, palette); + double width = 0.14; +// double width = 0.2; + + if (ROTATE_COLOR_SCHEME) + draw_color_scheme_palette_3d(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range, palette, fade, fade_value); + else + draw_color_scheme_palette_3d(XMAX - 1.5*width, YMIN + 0.1, XMAX - 0.5*width, YMAX - 0.1, plot, -range, range, palette, fade, fade_value); } @@ -550,9 +563,6 @@ void animation() courant2 = COURANT*COURANT; courantb2 = COURANTB*COURANTB; - - - /* initialize color scale, for option RESCALE_COLOR_IN_CENTER */ if (RESCALE_COLOR_IN_CENTER) { @@ -576,84 +586,33 @@ void animation() ratio = (XMAX - XMIN)/8.4; /* for Tokarsky billiard */ -// isospectral_initial_point(0.2, 0.0, startleft, startright); /* for isospectral billiards */ -// homophonic_initial_point(0.5, -0.25, 1.5, -0.25, startleft, startright); -// homophonic_initial_point(0.5, -0.25, 1.5, -0.25, startleft, startright); -// printf("xleft = (%.3f, %.3f) xright = (%.3f, %.3f)\n", startleft[0], startleft[1], startright[0], startright[1]); - -// xy_to_ij(startleft[0], startleft[1], sample_left); -// xy_to_ij(startright[0], startright[1], sample_right); -// printf("xleft = (%.3f, %.3f) xright = (%.3f, %.3f)\n", xin_left, yin_left, xin_right, yin_right); - -// init_wave_flat(phi, psi, xy_in); - -// init_wave_plus(LAMBDA - 0.3*MU, 0.5*MU, phi, psi, xy_in); -// init_wave(LAMBDA - 0.3*MU, 0.5*MU, phi, psi, xy_in); -// init_circular_wave(X_SHOOTER, Y_SHOOTER, phi, psi, xy_in); -// printf("Initializing wave\n"); // init_circular_wave_mod(polyline[85].x, polyline[85].y, phi, psi, xy_in); -// init_circular_wave_mod(0.0, 0.0, phi, psi, xy_in); - - init_circular_wave_mod(0.2, 0.4, phi, psi, xy_in); - add_circular_wave_mod(-1.0, -0.2, -0.4, phi, psi, xy_in); - - + init_circular_wave_mod(0.0, 0.0, phi, psi, xy_in); // add_circular_wave(-1.0, -0.2, -0.4, phi, psi, xy_in); // printf("Wave initialized\n"); - -// init_circular_wave(0.6*cos((double)(period)*DPI/3.0), 0.6*sin((double)(period)*DPI/3.0), phi, psi, xy_in); -// period++; -// for (i=0; i<3; i++) -// { -// add_circular_wave(-1.0, 0.6*cos(PID + (double)(i)*DPI/3.0), 0.6*sin(PID + (double)(i)*DPI/3.0), phi, psi, xy_in); -// } -// add_circular_wave(1.0, -LAMBDA, 0.0, phi, psi, xy_in); -// add_circular_wave(-1.0, 0.0, -LAMBDA, phi, psi, xy_in); - -// init_circular_wave_xplusminus(startleft[0], startleft[1], startright[0], startright[1], phi, psi, xy_in); - -// init_circular_wave_xplusminus(-0.9, 0.0, 0.81, 0.0, phi, psi, xy_in); -// init_circular_wave(-2.0*ratio, 0.0, phi, psi, xy_in); -// init_planar_wave(XMIN + 0.015, 0.0, phi, psi, xy_in); -// init_planar_wave(XMIN + 0.02, 0.0, phi, psi, xy_in); -// init_planar_wave(XMIN + 0.5, 0.0, phi, psi, xy_in); -// init_wave(-1.5, 0.0, phi, psi, xy_in); -// init_wave(0.0, 0.0, phi, psi, xy_in); - - /* add a drop at another point */ -// add_drop_to_wave(1.0, 0.7, 0.0, phi, psi); -// add_drop_to_wave(1.0, -0.7, 0.0, phi, psi); -// add_drop_to_wave(1.0, 0.0, -0.7, phi, psi); /* initialize table of wave speeds/dissipation */ - for (i=0; i= INITIAL_TIME)&&(DOUBLE_MOVIE)) { - draw_wave_3d(phi, psi, xy_in, wave, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 0, 1.0); - if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B); + draw_wave_3d(1, phi, psi, xy_in, wave, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 0, 1.0, 1); + if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0); glutSwapBuffers(); save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter); -// save_frame_counter(NSTEPS + 21 + counter); counter++; } @@ -735,27 +695,29 @@ void animation() { if (DOUBLE_MOVIE) { - draw_wave_3d(phi, psi, xy_in, wave, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0); - if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE); + draw_wave_3d(0, phi, psi, xy_in, wave, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1); + if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 0, 1.0); glutSwapBuffers(); if (!FADE) for (i=0; i