Update README.md

Parameters.md

@@ -10,367 +10,393 @@ updated gradually.
 
 
 
-### 2 August 23 - Faraday waves with 14-fold symmetry ###
+### 2 September 23 - A rotating triangle in a granular medium ###
 
-**Program:** `rde.c` 
-
-**Initial condition in function `animation()`:** `init_circular_vibration(0.0, 0.0, 0.015, LAMBDA, 0.001, 0.1, 14.0, SWATER_MIN_HEIGHT, phi, xy_in);`
+**Program:** `lennardjones.c` 
 
 ```
-#define MOVIE 1         /* set to 1 to generate movie */
-#define DOUBLE_MOVIE 1  /* set to 1 to produce movies for wave height and energy simultaneously */
-#define SAVE_MEMORY 1           /* set to 1 to save memory when writing tiff images */
+#define MOVIE 1        /* set to 1 to generate movie */
+#define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */
+#define SAVE_MEMORY 1   /* set to 1 to save memory while saving frames */
 #define NO_EXTRA_BUFFER_SWAP 1    /* some OS require one less buffer swap when recording images */
 
+#define TIME_LAPSE 0     /* set to 1 to add a time-lapse movie at the end */
+                         /* so far incompatible with double movie */
+#define TIME_LAPSE_FACTOR 3    /* factor of time-lapse movie */
+#define TIME_LAPSE_FIRST 1  /* set to 1 to show time-lapse version first */
+
+#define SAVE_TIME_SERIES 0  /* set to 1 to save time series of particle positions */
+
 /* General geometrical parameters */
 
-#define WINWIDTH 	1920  /* window width */
-#define WINHEIGHT 	1150  /* window height */
-#define NX 960          /* number of grid points on x axis */
-#define NY 575          /* number of grid points on y axis */
+#define WINWIDTH 	1280  /* window width */
+#define WINHEIGHT 	720   /* window height */
 
 #define XMIN -2.0
 #define XMAX 2.0	/* x interval */
-#define YMIN -1.197916667
-#define YMAX 1.197916667	/* y interval for 9/16 aspect ratio */
+#define YMIN -1.125
+#define YMAX 1.125	/* y interval for 9/16 aspect ratio */
 
-/* Choice of simulated equation */
+#define INITXMIN -2.0
+#define INITXMAX 2.0	/* x interval for initial condition */
+#define INITYMIN -1.125
+#define INITYMAX 1.125	/* y interval for initial condition */
 
-#define RDE_EQUATION 8  /* choice of reaction term, see list in global_3d.c */
-#define NFIELDS 3       /* number of fields in reaction-diffusion equation */
-#define NLAPLACIANS 0   /* number of fields for which to compute Laplacian */
+#define ADDXMIN -1.97
+#define ADDXMAX -0.8	/* x interval for adding particles */
+#define ADDYMIN -1.125
+#define ADDYMAX 1.125	/* y interval for adding particles */
 
-#define ADD_POTENTIAL 0 /* set to 1 to add a potential (for Schrodinger equation) */
-#define ADD_MAGNETIC_FIELD 0    /* set to 1 to add a magnetic field (for Schrodinger equation) - then set POTENTIAL 1 */
-#define ADD_FORCE_FIELD 1   /* set to 1 to add a foce field (for compressible Euler equation) */
-#define POTENTIAL 7         /* type of potential or vector potential, see list in global_3d.c  */
-#define FORCE_FIELD 5       /* type of force field, see list in global_3d.c  */
-#define ADD_CORIOLIS_FORCE 0    /* set to 1 to add Coriolis force (quasigeostrophic Euler equations) */
-#define VARIABLE_DEPTH 0    /* set to 1 for variable depth in shallow water equation */
-#define SWATER_DEPTH 4      /* variable depth in shallow water equation */
 
-#define ANTISYMMETRIZE_WAVE_FCT 0   /* set tot 1 to make wave function antisymmetric */
-#define ADAPT_STATE_TO_BC 1     /* to smoothly adapt initial state to obstacles */
-#define OBSTACLE_GEOMETRY 1     /* geometry of obstacles, as in B_DOMAIN */
-#define BC_STIFFNESS 50.0        /* controls region of boundary condition control */
-#define CHECK_INTEGRAL 1     /* set to 1 to check integral of first field */
+#define 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 JULIA_SCALE 0.5 /* scaling for Julia sets */
+#define OBSXMIN -2.0
+#define OBSXMAX 2.0     /* x interval for motion of obstacle */
 
-/* Choice of the billiard table */
+#define CIRCLE_PATTERN 1  /* pattern of circles, see list in global_ljones.c */
 
-#define B_DOMAIN 999          /* choice of domain shape, see list in global_pdes.c  */
+#define ADD_INITIAL_PARTICLES 0 /* set to 1 to add a second type of particles */
+#define CIRCLE_PATTERN_B 1  /* pattern of circles for additional particles */
 
-#define CIRCLE_PATTERN 8     /* pattern of circles, see list in global_pdes.c */
+#define ADD_FIXED_OBSTACLES 0   /* set to 1 do add fixed circular obstacles */
+#define OBSTACLE_PATTERN 5  /* 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 21      /* pattern of repelling segments, see list in global_ljones.c */
+#define ROCKET_SHAPE 3        /* shape of rocket combustion chamber, see list in global_ljones.c */
+#define ROCKET_SHAPE_B 3      /* shape of second rocket */
+#define NOZZLE_SHAPE 6        /* shape of nozzle, see list in global_ljones.c */
+#define NOZZLE_SHAPE_B 6      /* shape of nozzle for second rocket, see list in global_ljones.c */
+
+#define TWO_TYPES 0         /* set to 1 to have two types of particles */
+#define TYPE_PROPORTION 0.5 /* proportion of particles of first type */
+#define SYMMETRIZE_FORCE 1  /* set to 1 to symmetrize two-particle interaction, only needed if particles are not all the same */
+#define CENTER_PX 0         /* set to 1 to center horizontal momentum */
+#define CENTER_PY 0         /* set to 1 to center vertical momentum */
+#define CENTER_PANGLE 0     /* set to 1 to center angular momentum */
+
+#define INTERACTION 1       /* particle interaction, see list in global_ljones.c */
+#define INTERACTION_B 1     /* particle interaction for second type of particle, see list in global_ljones.c */
+#define SPIN_INTER_FREQUENCY 4.0 /* angular frequency of spin-spin interaction */
+#define SPIN_INTER_FREQUENCY_B 4.0 /* angular frequency of spin-spin interaction for second particle type */
 
 #define P_PERCOL 0.25       /* probability of having a circle in C_RAND_PERCOL arrangement */
-#define NPOISSON 300        /* number of points for Poisson C_RAND_POISSON arrangement */
+#define NPOISSON 100        /* number of points for Poisson C_RAND_POISSON arrangement */
+#define PDISC_DISTANCE 2.8  /* minimal distance in Poisson disc process, controls density of particles */
+#define PDISC_CANDIDATES 100 /* number of candidates in construction of Poisson disc process */
 #define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */
 
-#define LAMBDA 1.0	    /* parameter controlling the dimensions of domain */
-#define MU 0.06	            /* parameter controlling the dimensions of domain */
-#define NPOLY 5             /* number of sides of polygon */
-#define APOLY 2.0          /* angle by which to turn polygon, in units of Pi/2 */
-#define MDEPTH 7            /* depth of computation of Menger gasket */
-#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 LAMBDA 0.75	    /* parameter controlling the dimensions of domain */
+#define MU 0.015 	    /* parameter controlling radius of particles */
+#define MU_B 0.015          /* parameter controlling radius of particles of second type */
+#define NPOLY 3             /* number of sides of polygon */
+#define APOLY 0.0           /* angle by which to turn polygon, in units of Pi/2 */ 
+#define AWEDGE 0.5          /* opening angle of wedge, 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 6            /* number of grid point for grid of disks */
-#define NGRIDY 8            /* number of grid point for grid of disks */
-#define REVERSE_TESLA_VALVE 1   /* set to 1 to orient Tesla valve in blocking configuration */
+#define NGRIDX 128           /* number of grid point for grid of disks */
+#define NGRIDY 64           /* number of grid point for grid of disks */
+#define EHRENFEST_RADIUS 0.9    /* radius of container for Ehrenfest urn configuration */
+#define EHRENFEST_WIDTH 0.035     /* width of tube for Ehrenfest urn configuration */
+#define TWO_CIRCLES_RADIUS_RATIO 0.8    /* ratio of radii for S_TWO_CIRCLES_EXT segment configuration */
+#define DAM_WIDTH 0.05       /* width of dam for S_DAM segment configuration */
 
 #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 parameters of wave equation */
-
-#define DT 0.00000025
-
-#define VISCOSITY 1.5e-5
-#define DISSIPATION 1.0e-8
-
-#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.0  /* spring constant of harmonic potential */
-#define K_COULOMB 0.5   /* constant in Coulomb potential */
-#define V_MAZE 0.4      /* potential in walls of maze */
-#define BZQ 0.0008      /* parameter in BZ equation */
-#define BZF 1.2         /* parameter in BZ equation */
-#define B_FIELD 10.0    /* magnetic field */
-#define G_FIELD 0.75     /* gravity */
-#define BC_FIELD 1.0e-5     /* constant in repulsive field from obstacles */
-#define AB_RADIUS 0.2   /* radius of region with magnetic field for Aharonov-Bohm effect */
-#define K_EULER 50.0    /* constant in stream function integration of Euler equation */
-#define K_EULER_INC 0.5    /* constant in incompressible Euler equation */
-
-#define SMOOTHEN_VORTICITY 0    /* set to 1 to smoothen vorticity field in Euler equation */
-#define SMOOTHEN_VELOCITY 1     /* set to 1 to smoothen velocity field in Euler equation */
-#define SMOOTHEN_PERIOD 10      /* period between smoothenings */
-#define SMOOTH_FACTOR 0.15       /* factor by which to smoothen */
-
-#define ADD_TRACERS 1    /* set to 1 to add tracer particles (for Euler equations) */
-#define N_TRACERS 1000    /* number of tracer particles */
-#define TRACERS_STEP 0.005  /* step size in tracer evolution */
-
-#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.005      /* noise intensity */
-#define CHANGE_NOISE 1      /* set to 1 to increase noise intensity */
-#define NOISE_FACTOR 40.0   /* factor by which to increase noise intensity */
-#define NOISE_INITIAL_TIME 100  /* initial time during which noise remains constant */
-
-#define CHANGE_VISCOSITY 0      /* set to 1 to change the viscosity in the course of the simulation */
-#define ADJUST_INTSTEP 0       /* 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 */
-                                      
-#define CHANGE_FLOW_SPEED 0     /* set to 1 to change speed of laminar flow */
-#define IN_OUT_FLOW_BC 0          /* type of in-flow/out-flow boundary conditions for Euler equation, 0 for no b.c. */
-#define IN_OUT_BC_FACTOR 0.001    /* factor of convex combination between old and new flow */
-#define BC_FLOW_TYPE 1            /* type of initial condition */
-                                  /* see list in global_pdes.c */
-#define IN_OUT_FLOW_MIN_AMP 0.25   /* amplitude of in-flow/out-flow boundary conditions (for Euler equation) - min value */
-#define IN_OUT_FLOW_AMP 0.25       /* amplitude of in-flow/out-flow boundary conditions (for Euler equation) - max value */
-#define LAMINAR_FLOW_MODULATION 0.01     /* asymmetry of laminar flow */
-#define LAMINAR_FLOW_YPERIOD 1.0    /* period of laminar flow in y direction */
-#define PRESSURE_GRADIENT 0.2       /* amplitude of pressure gradient for Euler equation */
-
-#define SWATER_MIN_HEIGHT 0.5      /* min height of initial condition for shallow water equation */
-#define DEPTH_FACTOR 0.015        /* proportion of min height in variable depth */
-#define TANH_FACTOR 1.0         /* steepness of variable depth */
-
-#define EULER_GRADIENT_YSHIFT 0.0    /* y-shift in computation of gradient in Euler equation */
-
-/* Boundary conditions, see list in global_pdes.c  */
-
-#define B_COND 0
-
-#define B_COND_LEFT 0
-#define B_COND_RIGHT 0
-#define B_COND_TOP 0
-#define B_COND_BOTTOM 0
-
-
 /* Parameters for length and speed of simulation */
 
-#define NSTEPS 1400         /* number of frames of movie */
-#define NVID 50            /* 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 999         /* number of segments of boundary */
-#define BOUNDARY_WIDTH 2    /* width of billiard boundary */
+#define NSTEPS 3600      /* number of frames of movie */
+#define NVID 30         /* number of iterations between images displayed on screen */
+#define NSEG 250         /* number of segments of boundary */
+#define INITIAL_TIME 0     /* time after which to start saving frames */
+#define OBSTACLE_INITIAL_TIME 150     /* time after which to start moving obstacle */
+#define BOUNDARY_WIDTH 1    /* width of particle boundary */
+#define LINK_WIDTH 2        /* width of links between particles */
+#define CONTAINER_WIDTH 4   /* width of container boundary */
 
-#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 50    /* number of still frames at end of movie */
-#define FADE 1           /* set to 1 to fade at end of movie */
+#define PAUSE 1000         /* number of frames after which to pause */
+#define PSLEEP 1         /* sleep time during pause */
+#define SLEEP1  1        /* initial sleeping time */
+#define SLEEP2  1   /* final sleeping time */
+#define MID_FRAMES 20    /* number of still frames between parts of two-part movie */
+#define END_FRAMES 100   /* number of still frames at end of movie */
 
-/* Visualisation */
+/* Boundary conditions, see list in global_ljones.c */
 
-#define PLOT_3D 1    /* controls whether plot is 2D or 3D */
+#define BOUNDARY_COND 3
 
-#define ROTATE_VIEW 0       /* set to 1 to rotate position of observer */
-#define ROTATE_ANGLE 360.0  /* total angle of rotation during simulation */
+/* Plot type, see list in global_ljones.c  */
 
-#define DRAW_PERIODICISED 0     /* set to 1 to repeat wave periodically in x and y directions */
+#define PLOT 11
+#define PLOT_B 15        /* plot type for second movie */
 
-/* Plot type - color scheme */
+#define DRAW_BONDS 1    /* set to 1 to draw bonds between neighbours */
+#define COLOR_BONDS 1   /* set to 1 to color bonds according to length */
+#define FILL_TRIANGLES 0    /* set to 1 to fill triangles between neighbours */
+#define ALTITUDE_LINES 0    /* set to 1 to add horizontal lines to show altitude */
+#define COLOR_SEG_GROUPS 1  /* set to 1 to collor segment groups differently */
+#define N_PARTICLE_COLORS 200   /* number of colors for P_NUMBER color scheme */
+#define INITIAL_POS_TYPE 1      /* type of initial position dependence */
+#define ERATIO 0.995          /* ratio for time-averagin in P_EMEAN color scheme */
+#define DRATIO 0.995          /* ratio for time-averagin in P_DIRECT_EMEAN color scheme */
 
-#define CPLOT 70
-#define CPLOT_B 74
+/* Color schemes */
 
-/* Plot type - height of 3D plot */
+#define COLOR_PALETTE 17     /* Color palette, see list in global_ljones.c  */
+#define COLOR_PALETTE_EKIN 10   /* Color palette for kinetic energy */
+#define COLOR_PALETTE_ANGLE 10  /* Color palette for angle representation */
+#define COLOR_PALETTE_DIRECTION 17  /* Color palette for direction representation */
+#define COLOR_PALETTE_INITIAL_POS 0  /* Color palette for initial position representation */
 
-#define ZPLOT 70     /* z coordinate in 3D plot */
-#define ZPLOT_B 71    /* z coordinate in second 3D plot */
+#define BLACK 1          /* background */
 
-#define AMPLITUDE_HIGH_RES 1    /* set to 1 to increase resolution of P_3D_AMPLITUDE plot */
-#define SHADE_3D 1              /* set to 1 to change luminosity according to normal vector */
-#define NON_DIRICHLET_BC 0      /* set to 1 to draw only facets in domain, if field is not zero on boundary */
-#define WRAP_ANGLE 1            /* experimental: wrap angle to [0, 2Pi) for interpolation in angle schemes */
-#define FADE_IN_OBSTACLE 0      /* set to 1 to fade color inside obstacles */
-#define FADE_WATER_DEPTH 0      /* set to 1 to make wave color depth-dependent */
-#define DRAW_OUTSIDE_GRAY 0     /* experimental - draw outside of billiard in gray */
-#define ADD_POTENTIAL_TO_Z 0    /* set to 1 to add the external potential to z-coordinate of plot */
-#define ADD_POT_CONSTANT 0.35   /* constant added to wave height */
-#define ADD_HEIGHT_CONSTANT -0.5   /* constant added to wave height */
-#define DRAW_DEPTH 0            /* set to 1 to draw water depth */
-#define DEPTH_SCALE 0.75         /* vertical scaling of depth plot */
-#define DEPTH_SHIFT -0.015      /* vertical shift of depth plot */
-
-#define PLOT_SCALE_ENERGY 0.05      /* vertical scaling in energy plot */
-
-#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 PRINT_PROBABILITIES 0   /* set to 1 to print probabilities (for Ehrenfest urn configuration) */
-#define PRINT_NOISE 0       /* set to 1 to print noise intensity */
-#define PRINT_FLOW_SPEED 0      /* set to 1 to print speed of flow */
-#define PRINT_AVERAGE_SPEED 0   /* set to 1 to print average speed of flow */
-#define PRINT_LEFT 1            /* set to 1 to print parameters at left side */
-
-#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 13       /* Color palette, see list in global_pdes.c  */
-#define COLOR_PALETTE_B 0      /* Color palette, see list in global_pdes.c  */
-
-#define BLACK 1          /* black background */
-
-#define COLOR_SCHEME 3   /* choice of color scheme */
-
-#define COLOR_PHASE_SHIFT 0.5   /* phase shift of color scheme, in units of Pi */
+#define COLOR_SCHEME 3   /* choice of color scheme, see list in global_ljones.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 10.0      /* additional scaling factor for color scheme P_3D_AMPLITUDE */
+#define SLOPE 0.5        /* sensitivity of color on wave 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 200.0       /* sensitivity of luminosity on module, for color scheme Z_ARGUMENT */
-#define MIN_SCHROD_LUM 0.2       /* minimal luminosity in color scheme Z_ARGUMENT*/
-#define VSCALE_PRESSURE 2.0      /* additional scaling factor for color scheme Z_EULER_PRESSURE */
-#define PRESSURE_SHIFT 10.0        /* shift for color scheme Z_EULER_PRESSURE */
-#define PRESSURE_LOG_SHIFT -2.5     /* shift for color scheme Z_EULER_PRESSURE */
-#define VSCALE_WATER_HEIGHT 0.4     /* vertical scaling of water height */
 
 #define COLORHUE 260     /* initial hue of water color for scheme C_LUM */
 #define COLORDRIFT 0.0   /* how much the color hue drifts during the whole simulation */
 #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 FLUX_SCALE 100.0 /* scaling factor for energy representation */
-#define LOG_SCALE 0.5    /* scaling factor for energy log representation */
-#define LOG_SHIFT 1.0   
-#define LOG_MIN 1.0e-3   /* floor value for log vorticity plot */
-#define VSCALE_SPEED 200.0      /* additional scaling factor for color scheme Z_EULER_SPEED */
-#define VMEAN_SPEED 0.0       /* mean value around which to scale for color scheme Z_EULER_SPEED */
-#define SHIFT_DENSITY 8.5        /* shift for color scheme Z_EULER_DENSITY */
-#define VSCALE_DENSITY 3.0        /* additional scaling factor for color scheme Z_EULER_DENSITY */
-#define VSCALE_VORTICITY 20.0     /* additional scaling factor for color scheme Z_EULERC_VORTICITY */
-#define VORTICITY_SHIFT 0.0     /* vertical shift of vorticity */
-#define ZSCALE_SPEED 0.3        /* additional scaling factor for z-coord Z_EULER_SPEED and Z_SWATER_SPEED */
-#define VSCALE_SWATER 300.0        /* additional scaling factor for color scheme Z_EULER_DENSITY */
+#define HUEMEAN 220.0    /* mean value of hue for color scheme C_HUE */
+#define HUEAMP -50.0      /* amplitude of variation of hue for color scheme C_HUE */
+#define COLOR_HUESHIFT 0.5     /* shift in color hue (for some cyclic palettes) */
 
-#define NXMAZE 7      /* width of maze */
-#define NYMAZE 7      /* height of maze */
+#define PRINT_PARAMETERS 1  /* set to 1 to print certain parameters */
+#define PRINT_TEMPERATURE 0 /* set to 1 to print current temperature */
+
+/* particle properties */
+
+#define ENERGY_HUE_MIN 330.0        /* color of original particle */
+#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 10000.0        /* energy of particle with hottest color */
+#define HUE_TYPE0 60.0     /* hue of particles of type 0 */
+#define HUE_TYPE1 280.0      /* hue of particles of type 1 */
+#define HUE_TYPE2 140.0      /* hue of particles of type 2 */
+#define HUE_TYPE3 200.0     /* hue of particles of type 3 */
+
+#define RANDOM_RADIUS 0     /* set to 1 for random circle radius */
+#define DT_PARTICLE 3.0e-6    /* time step for particle displacement */
+#define KREPEL 1.0          /* constant in repelling force between particles */
+#define EQUILIBRIUM_DIST 5.0    /* Lennard-Jones equilibrium distance */
+#define EQUILIBRIUM_DIST_B 5.0  /* Lennard-Jones equilibrium distance for second type of particle */
+#define REPEL_RADIUS 15.0    /* radius in which repelling force acts (in units of particle radius) */
+#define DAMPING 100.0          /* damping coefficient of particles */
+#define OMEGA_INITIAL 10.0        /* initial angular velocity range */
+#define INITIAL_DAMPING 5.0  /* damping coefficient of particles during initial phase */
+#define DAMPING_ROT 100.0      /* dampint coefficient for rotation of particles */
+#define PARTICLE_MASS 2.5    /* mass of particle of radius MU */
+#define PARTICLE_MASS_B 1.0   /* mass of particle of radius MU */
+#define PARTICLE_INERTIA_MOMENT 0.2     /* moment of inertia of particle */
+#define PARTICLE_INERTIA_MOMENT_B 0.02     /* moment of inertia of second type of particle */
+#define V_INITIAL 0.0        /* initial velocity range */
+#define OMEGA_INITIAL 10.0        /* initial angular velocity range */
+#define VICSEK_VMIN 1.0    /* minimal speed of particles in Vicsek model */
+#define VICSEK_VMAX 40.0    /* minimal speed of particles in Vicsek model */
+
+#define THERMOSTAT 0        /* set to 1 to switch on thermostat */
+#define VARY_THERMOSTAT 0   /* set to 1 for time-dependent thermostat schedule */
+#define SIGMA 5.0           /* noise intensity in thermostat */
+#define BETA 0.02          /* initial inverse temperature */
+#define MU_XI 0.01           /* friction constant in thermostat */
+#define KSPRING_BOUNDARY 1.0e7    /* confining harmonic potential outside simulation region */
+#define KSPRING_OBSTACLE 1.0e11    /* harmonic potential of obstacles */
+#define NBH_DIST_FACTOR 2.7       /* radius in which to count neighbours */
+#define GRAVITY 0.0             /* gravity acting on all particles */
+#define GRAVITY_X 0.0          /* horizontal gravity acting on all particles */
+#define INCREASE_GRAVITY 0     /* set to 1 to increase gravity during the simulation */
+#define GRAVITY_SCHEDULE 2     /* type of gravity schedule, see list in global_ljones.c */
+#define GRAVITY_FACTOR 100.0    /* factor by which to increase gravity */
+#define GRAVITY_INITIAL_TIME 200    /* time at start of simulation with constant gravity */
+#define GRAVITY_RESTORE_TIME 700    /* time at end of simulation with gravity restored to initial value */
+#define KSPRING_VICSEK 0.2   /* spring constant for I_VICSEK_SPEED interaction */
+#define VICSEK_REPULSION 10.0    /* repulsion between particles in Vicsek model */
+
+#define ROTATION 0           /* set to 1 to include rotation of particles */
+#define COUPLE_ANGLE_TO_THERMOSTAT 1    /* set to 1 to couple angular degrees of freedom to thermostat */
+#define DIMENSION_FACTOR 1.0  /* scaling factor taking into account number of degrees of freedom */  
+#define KTORQUE 50.0         /* force constant in angular dynamics */
+#define KTORQUE_BOUNDARY 1.0e6  /* constant in torque from the boundary */
+#define KTORQUE_B 10.0        /* force constant in angular dynamics */
+#define KTORQUE_DIFF 150.0    /* force constant in angular dynamics for different particles */
+#define DRAW_SPIN 0           /* set to 1 to draw spin vectors of particles */
+#define DRAW_SPIN_B 0         /* set to 1 to draw spin vectors of particles */
+#define DRAW_CROSS 1          /* set to 1 to draw cross on particles of second type */
+#define SPIN_RANGE 10.0       /* range of spin-spin interaction */
+#define SPIN_RANGE_B 5.0     /* range of spin-spin interaction for second type of particle */
+#define QUADRUPOLE_RATIO 0.6  /* anisotropy in quadrupole potential */ 
+
+#define INCREASE_BETA 0  /* set to 1 to increase BETA during simulation */
+#define BETA_FACTOR 0.02   /* factor by which to change BETA during simulation */
+#define N_TOSCILLATIONS 1.5   /* number of temperature oscillations in BETA schedule */
+#define NO_OSCILLATION 1        /* set to 1 to have exponential BETA change only */
+#define MIDDLE_CONSTANT_PHASE 2000   /* final phase in which temperature is constant */
+#define FINAL_DECREASE_PHASE 1300    /* final phase in which temperature decreases */ 
+#define FINAL_CONSTANT_PHASE -1     /* 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 */
+#define COMPRESSION_RATIO 0.3       /* final size of container */
+#define RESTORE_CONTAINER_SIZE 1    /* set to 1 to restore container to initial size at end of simulation */
+#define RESTORE_TIME 700            /* time before end of sim at which to restore size */
+
+#define MOVE_OBSTACLE 0     /* set to 1 to have a moving obstacle */
+#define CENTER_VIEW_ON_OBSTACLE 0   /* set to 1 to center display on moving obstacle */
+#define RESAMPLE_Y 0         /* set to 1 to resample y coordinate of moved particles (for shock waves) */
+#define NTRIALS 2000         /* number of trials when resampling */
+#define OBSTACLE_RADIUS 0.25  /* radius of obstacle for circle boundary conditions */
+#define FUNNEL_WIDTH  0.25   /* funnel width for funnel boundary conditions */
+#define OBSTACLE_XMIN 0.0    /* initial position of obstacle */
+#define OBSTACLE_XMAX 3.0    /* final position of obstacle */
+#define RECORD_PRESSURES 0   /* set to 1 to record pressures on obstacle */
+#define N_PRESSURES 100      /* number of intervals to record pressure */
+#define N_P_AVERAGE 100      /* size of pressure averaging window */
+#define N_T_AVERAGE 1        /* size of temperature averaging window */
+#define MAX_PRESSURE 3.0e10  /* pressure shown in "hottest" color */
+#define PARTIAL_THERMO_COUPLING 0   /* set to 1 to couple only some particles to thermostat */
+#define PARTIAL_THERMO_REGION 1     /* region for partial thermostat coupling (see list in global_ljones.c) */
+#define PARTIAL_THERMO_SHIFT 0.2    /* distance from obstacle at the right of which particles are coupled to thermostat */
+#define PARTIAL_THERMO_WIDTH 0.5    /* vertical size of partial thermostat coupling */
+#define PARTIAL_THERMO_HEIGHT 0.25   /* vertical size of partial thermostat coupling */
+
+#define INCREASE_KREPEL 0   /* set to 1 to increase KREPEL during simulation */
+#define KREPEL_FACTOR 1000.0   /* factor by which to change KREPEL during simulation */
+
+#define PART_AT_BOTTOM 0     /* set to 1 to include "seed" particles at bottom */
+#define MASS_PART_BOTTOM 10000.0 /* mass of particles at bottom */
+#define NPART_BOTTOM 100     /* number of particles at the bottom */
+
+#define ADD_PARTICLES 0   /* set to 1 to add particles */
+#define ADD_TIME 0       /* time at which to add first particle */
+#define ADD_PERIOD 10       /* time interval between adding further particles */
+#define N_ADD_PARTICLES 2  /* number of particles to add */
+#define FINAL_NOADD_PERIOD 0  /* final period where no particles are added */
+#define SAFETY_FACTOR 2.0  /* no particles are added at distance less than MU*SAFETY_FACTOR of other particles */
+
+#define TRACER_PARTICLE 0   /* set to 1 to have a tracer particle */
+#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 1000  /* period of rotating boundary */
+#define ROTATE_INITIAL_TIME 300       /* initial time without rotation */
+#define ROTATE_FINAL_TIME 600       /* final time without rotation */
+#define ROTATE_CHANGE_TIME 0.5     /* relative duration of acceleration/deceleration phases */
+#define OMEGAMAX 8.0*PI              /* maximal rotation speed */
+#define PRINT_ANGLE 1               /* set to 1 to print obstacle orientation */
+#define PRINT_OMEGA 1               /* set to 1 to print angular speed */
+#define PRINT_PARTICLE_SPEEDS 0     /* set to 1 to print average speeds/momenta of particles */
+#define PRINT_SEGMENTS_SPEEDS 0     /* set to 1 to print velocity of moving segments */
+
+#define MOVE_BOUNDARY 0        /* set to 1 to move repelling segments, due to force from particles */
+#define SEGMENTS_MASS 40.0     /* mass of collection of segments */
+#define DEACTIVATE_SEGMENT 1    /* set to 1 to deactivate last segment after a certain time */
+#define SEGMENT_DEACTIVATION_TIME 20   /* time at which to deactivate last segment */
+#define RELEASE_ROCKET_AT_DEACTIVATION 0    /* set to 1 to limit segments velocity before segment release */
+#define SEGMENTS_X0 1.5        /* initial position of segments */
+#define SEGMENTS_Y0 0.0        /* initial position of segments */
+#define SEGMENTS_VX0 0.0       /* initial velocity of segments */
+#define SEGMENTS_VY0 0.0      /* initial velocity of segments */
+#define DAMP_SEGS_AT_NEGATIVE_Y 0   /* set to 1 to dampen segments when y coordinate is negative */
+
+#define MOVE_SEGMENT_GROUPS 0       /* set to 1 to group segments into moving units */
+#define SEGMENT_GROUP_MASS 500.0   /* mass of segment group */
+#define SEGMENT_GROUP_I 1000.0      /* moment of inertia of segment group */
+#define SEGMENT_GROUP_DAMPING 0.0   /* damping of segment groups */
+#define GROUP_REPULSION 0           /* set to 1 for groups of segments to repel each other */
+#define KSPRING_GROUPS 5.0e11       /* harmonic potential between segment groups */
+#define GROUP_WIDTH 0.05            /* interaction width of groups */
+#define GROUP_G_REPEL 0             /* set to 1 to add repulsion between centers of mass of groups */
+#define GROUP_G_REPEL_RADIUS 1.2    /* radius within which centers of mass of groups repel each other */
+#define TRACK_SEGMENT_GROUPS 0      /* set to 1 for view to track group of segments */
+#define TRACK_X_PADDING 2.0         /* distance from x boundary where tracking starts */
+
+#define POSITION_DEPENDENT_TYPE 0   /* set to 1 to make particle type depend on initial position */
+#define POSITION_Y_DEPENDENCE 0     /* set to 1 for the separation between particles to be horizontal */
+#define POSITION_DEP_SIGN -1.0      /* sign in position dependence condition */
+#define POSITION_DEP_X -0.625         /* threshold value for position-dependent type */
+#define PRINT_ENTROPY 0     /* set to 1 to compute entropy */
+
+#define REACTION_DIFFUSION 0    /* set to 1 to simulate a chemical reaction (particles may change type) */
+#define RD_REACTION 16          /* type of reaction, see list in global_ljones.c */
+#define RD_TYPES 5              /* number of types in reaction-diffusion equation */
+#define RD_INITIAL_COND 9       /* initial condition of particles */
+#define REACTION_DIST 4.0       /* maximal distance for reaction to occur */
+#define REACTION_PROB 1.0        /* probability controlling reaction term */ 
+#define DISSOCIATION_PROB 0.002  /* probability controlling dissociation reaction */ 
+#define CENTER_COLLIDED_PARTICLES 0  /* set to 1 to recenter particles upon reaction (may interfere with thermostat) */
+#define EXOTHERMIC 1            /* set to 1 to make reaction exo/endothermic */
+#define DELTA_EKIN 2000.0        /* change of kinetic energy in reaction */
+#define COLLISION_TIME 15       /* time during which collisions are shown */
+
+#define CHANGE_RADIUS 0         /* set to 1 to change particle radius during simulation */
+#define MU_RATIO 0.666666667    /* ratio by which to increase radius */
+
+#define PRINT_PARTICLE_NUMBER 0     /* set to 1 to print total number of particles */
+#define PLOT_PARTICLE_NUMBER 0      /* set to 1 to make of plot of particle number over time */
+#define PARTICLE_NB_PLOT_FACTOR 0.5 /* expected final number of particles over initial number */
+#define PRINT_LEFT 0        /* set to 1 to print certain parameters at the top left instead of right */
+#define PLOT_SPEEDS 0       /* set to 1 to add a plot of obstacle speeds (e.g. for rockets) */
+#define PLOT_TRAJECTORIES 0     /* set to 1 to add a plot of obstacle trajectories (e.g. for rockets) */
+#define VMAX_PLOT_SPEEDS 0.25    /* vertical scale of plot of obstacle speeds */
+
+#define EHRENFEST_COPY 0    /* set to 1 to add equal number of larger particles (for Ehrenfest model) */
+
+#define LID_MASS 1000.0     /* mass of lid for BC_RECTANGLE_LID b.c. */
+#define LID_WIDTH 0.1       /* width of lid for BC_RECTANGLE_LID b.c. */
+#define WALL_MASS 2000.0    /* mass of wall for BC_RECTANGLE_WALL b.c. */
+#define WALL_FRICTION 0.0   /* friction on wall for BC_RECTANGLE_WALL b.c. */
+#define WALL_WIDTH 0.1      /* width of wall for BC_RECTANGLE_WALL b.c. */
+#define WALL_VMAX 100.0     /* max speed of wall */
+#define WALL_TIME 0         /* time during which to keep wall */
+
+#define NXMAZE 12      /* width of maze */
+#define NYMAZE 12      /* height of maze */
 #define MAZE_MAX_NGBH 4     /* max number of neighbours of maze cell */
-#define RAND_SHIFT 3        /* seed of random number generator */
-#define MAZE_XSHIFT 0.0     /* horizontal shift of maze */
-#define MAZE_WIDTH 0.04     /* half width of maze walls */
+#define RAND_SHIFT 4        /* seed of random number generator */
+#define MAZE_XSHIFT 0.5     /* horizontal shift of maze */
+#define MAZE_WIDTH 0.01     /* width of maze walls */
 
-#define DRAW_COLOR_SCHEME 1     /* set to 1 to plot the color scheme */
-#define COLORBAR_RANGE 2.5      /* scale of color scheme bar */
-#define COLORBAR_RANGE_B 2.5    /* scale of color scheme bar for 2nd part */
+#define FLOOR_FORCE 1      /* set to 1 to limit force on particle to FMAX */
+#define FMAX 1.0e9         /* maximal force */
+#define FLOOR_OMEGA 0      /* set to 1 to limit particle momentum to PMAX */
+#define PMAX 1000.0        /* maximal force */
+
+#define HASHX 120   /* size of hashgrid in x direction */
+#define HASHY 60    /* 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 */
+
+#define DRAW_COLOR_SCHEME 0     /* set to 1 to plot the color scheme */
+#define COLORBAR_RANGE 8.0    /* 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 */
-#define CIRC_COLORBAR 0         /* set to 1 to draw circular color scheme */
-#define CIRC_COLORBAR_B 1       /* set to 1 to draw circular color scheme */
 
-/* only for compatibility with wave_common.c */
-#define TWOSPEEDS 0          /* set to 1 to replace hardcore boundary by medium with different speed */
-#define VARIABLE_IOR 0      /* set to 1 for a variable index of refraction */
-#define IOR 4               /* choice of index of refraction, see list in global_pdes.c */
-#define IOR_TOTAL_TURNS 1.5 /* total angle of rotation for IOR_PERIODIC_WELLS_ROTATING */
-#define MANDEL_IOR_SCALE -0.05   /* parameter controlling dependence of IoR on Mandelbrot escape speed */
-#define 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 */
-#define OSCILLATION_SCHEDULE 0  /* oscillation schedule, see list in global_pdes.c */
-#define AMPLITUDE 0.8      /* amplitude of periodic excitation */ 
-#define ACHIRP 0.2        /* acceleration coefficient in chirp */
-#define DAMPING 0.0        /* damping of periodic excitation */
-#define COMPARISON 0        /* set to 1 to compare two different patterns (beta) */
-#define B_DOMAIN_B 20       /* second domain shape, for comparisons */
-#define CIRCLE_PATTERN_B 0  /* second pattern of circles or polygons */
-#define FLUX_WINDOW 20      /* averaging window for energy flux */
-#define ADD_WAVE_PACKET_SOURCES 1       /* set to 1 to add several sources emitting wave packets */
-#define WAVE_PACKET_SOURCE_TYPE 1       /* type of wave packet sources */
-#define N_WAVE_PACKETS 15               /* number of wave packets */
-#define WAVE_PACKET_RADIUS 20            /* radius of wave packets */
-#define OSCIL_LEFT_YSHIFT 25.0   /* y-dependence of left oscillation (for non-horizontal waves) */
-/* 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, 7.0};    /* location of observer for REP_PROJ_3D representation */ 
-int reset_view = 0;         /* switch to reset 3D view parameters (for option ROTATE_VIEW) */
-
-#define Z_SCALING_FACTOR 75.0   /* overall scaling factor of z axis for REP_PROJ_3D representation */
-#define XY_SCALING_FACTOR 2.4  /* overall scaling factor for on-screen (x,y) coordinates after projection */
-#define ZMAX_FACTOR 1.0        /* max value of z coordinate for REP_PROJ_3D representation */
-#define XSHIFT_3D 0.0          /* overall x shift for REP_PROJ_3D representation */
-#define YSHIFT_3D 0.1         /* overall y shift for REP_PROJ_3D representation */
-#define BORDER_PADDING 0       /* distance from boundary at which to plot points, to avoid boundary effects due to gradient */
-
-/* For debugging purposes only */
-#define FLOOR 1         /* set to 1 to limit wave amplitude to VMAX */
-#define VMAX 100.0        /* max value of wave amplitude */
-#define TEST_GRADIENT 0 /* print norm squared of gradient */
+#define LIMIT_ENERGY 0      /* set to 1 to limit energy, when there is no thermostat */
 
 ```
 
-### 1 August 23 - Bragg diffraction ###
+### 1 September 23 - Simulation of a meteor impact in the South Pacific Ocean, viewpoint moving west ###
 
-**Program:** `wave_billiard.c` 
+**Program:** `wave_sphere.c` 
 
-**Initial condition in function `animation()`:** `init_wave_flat(phi, psi, xy_in);`
+**Initial condition in function `animation()`:** `init_circular_wave_sphere(2.6, -0.5, phi, psi, xy_in, wsphere);`
 
 ```
 #define MOVIE 1         /* set to 1 to generate movie */
@@ -378,53 +404,57 @@ int reset_view = 0;         /* switch to reset 3D view parameters (for option RO
 #define SAVE_MEMORY 1   /* set to 1 to save memory when writing tiff images */
 #define NO_EXTRA_BUFFER_SWAP 1    /* some OS require one less buffer swap when recording images */
 
-#define VARIABLE_IOR 0      /* set to 1 for a variable index of refraction */
-#define IOR 9               /* choice of index of refraction, see list in global_pdes.c */
-#define IOR_TOTAL_TURNS 1.5 /* total angle of rotation for IOR_PERIODIC_WELLS_ROTATING */
-#define MANDEL_IOR_SCALE -0.05   /* parameter controlling dependence of IoR on Mandelbrot escape speed */
-
-
 /* General geometrical parameters */
 
 #define WINWIDTH 	1920  /* window width */
 #define WINHEIGHT 	1150  /* window height */
-#define NX 3840          /* number of grid points on x axis */
-#define NY 2300          /* number of grid points on y axis */
+#define NX 2560          /* number of grid points on x axis */
+#define NY 1280           /* number of grid points on y axis */
+
+#define DPOLE 20         /* safety distance to poles */
+#define SMOOTHPOLE 0.1     /* smoothing coefficient at poles */
 
 #define XMIN -2.0
 #define XMAX 2.0	/* x interval */
-#define YMIN -0.797916667
-#define YMAX 1.597916667	/* y interval for 9/16 aspect ratio */
+#define YMIN -1.041666667
+#define YMAX 1.041666667	/* y interval for 9/16 aspect ratio */
 
-#define HIGHRES 1       /* set to 1 if resolution of grid is double that of displayed image */
+#define HIGHRES 0        /* set to 1 if resolution of grid is double that of displayed image */
 
-#define JULIA_SCALE 1.0 /* scaling for Julia sets */
+#define JULIA_SCALE 0.25  /* scaling for Julia sets */
+#define JULIA_RE -0.77145    
+#define JULIA_IM -0.10295    /* parameters for Julia sets */
 
 /* Choice of the billiard table */
 
-#define B_DOMAIN 20        /* choice of domain shape, see list in global_pdes.c */
+#define B_DOMAIN 84         /* choice of domain shape, see list in global_pdes.c */
 
-#define CIRCLE_PATTERN 101   /* pattern of circles or polygons, see list in global_pdes.c */
+#define CIRCLE_PATTERN 33   /* pattern of circles or polygons, see list in global_pdes.c */
 
-#define COMPARISON 0        /* set to 1 to compare two different patterns (beta) */
+#define COMPARISON 0        /* set to 1 to compare two different patterns */
 #define B_DOMAIN_B 20       /* second domain shape, for comparisons */
 #define CIRCLE_PATTERN_B 0  /* second pattern of circles or polygons */
 
+#define VARIABLE_IOR 0      /* set to 1 for a variable index of refraction */
+#define IOR 9               /* choice of index of refraction, see list in global_pdes.c */
+#define IOR_TOTAL_TURNS 1.0 /* total angle of rotation for IOR_PERIODIC_WELLS_ROTATING */
+#define MANDEL_IOR_SCALE -0.05   /* parameter controlling dependence of IoR on Mandelbrot escape speed */
+
 #define P_PERCOL 0.25       /* probability of having a circle in C_RAND_PERCOL arrangement */
 #define NPOISSON 1000        /* number of points for Poisson C_RAND_POISSON arrangement */
 #define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */
 
-#define LAMBDA 0.1197916667	    /* parameter controlling the dimensions of domain */
-#define MU 0.01             /* parameter controlling the dimensions of domain */
+#define LAMBDA 0.75	    /* parameter controlling the dimensions of domain */
+#define MU 0.1             /* parameter controlling the dimensions of domain */
 #define NPOLY 6             /* number of sides of polygon */
 #define APOLY 0.0           /* angle by which to turn polygon, in units of Pi/2 */ 
-#define MDEPTH 6            /* depth of computation of Menger gasket */
+#define MDEPTH 7            /* 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 MANDELLEVEL 2000    /* iteration level for Mandelbrot set */
+#define MANDELLIMIT 20.0    /* limit value for approximation of Mandelbrot set */
 #define FOCI 1              /* set to 1 to draw focal points of ellipse */
-#define NGRIDX 40            /* number of grid point for grid of disks */
-#define NGRIDY 8            /* number of grid point for grid of disks */
+#define NGRIDX 30            /* number of grid point for grid of disks */
+#define NGRIDY 18            /* number of grid point for grid of disks */
 
 #define X_SHOOTER -0.2
 #define Y_SHOOTER -0.6
@@ -437,30 +467,31 @@ int reset_view = 0;         /* switch to reset 3D view parameters (for option RO
 #define ISO_YSHIFT_RIGHT -0.15 
 #define ISO_SCALE 0.5           /* coordinates for isospectral billiards */
 
+
 /* You can add more billiard tables by adapting the functions */
 /* xy_in_billiard and draw_billiard below */
 
 /* Physical parameters of wave equation */
 
 #define TWOSPEEDS 0          /* set to 1 to replace hardcore boundary by medium with different speed */
-#define OSCILLATE_LEFT 1     /* set to 1 to add oscilating boundary condition on the left */
+#define OSCILLATE_LEFT 0     /* set to 1 to add oscilating boundary condition on the left */
 #define OSCILLATE_TOPBOT 0   /* set to 1 to enforce a planar wave on top and bottom boundary */
-#define OSCILLATION_SCHEDULE 0  /* oscillation schedule, see list in global_pdes.c */
+#define OSCILLATION_SCHEDULE 3  /* oscillation schedule, see list in global_pdes.c */
 
-#define OMEGA 0.012        /* frequency of periodic excitation */
-#define AMPLITUDE 1.0      /* amplitude of periodic excitation */ 
-#define ACHIRP 0.25        /* acceleration coefficient in chirp */
-#define DAMPING 0.0        /* damping of periodic excitation */
-#define COURANT 0.1        /* Courant number */
+#define OMEGA 0.001       /* frequency of periodic excitation */
+#define AMPLITUDE 0.8     /* amplitude of periodic excitation */ 
+#define ACHIRP 0.2        /* acceleration coefficient in chirp */
+#define DAMPING 0.0       /* damping of periodic excitation */
+#define COURANT 0.05       /* Courant number */
 #define COURANTB 0.01      /* Courant number in medium B */
 #define GAMMA 0.0          /* damping factor in wave equation */
-#define GAMMAB 0.0         /* damping factor in wave equation */
+#define GAMMAB 1.0e-6         /* 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 */
 #define KAPPA_SIDES 5.0e-4  /* "elasticity" term on absorbing boundary */
 #define KAPPA_TOPBOT 0.0    /* "elasticity" term on absorbing boundary */
-#define OSCIL_LEFT_YSHIFT -200.0   /* y-dependence of left oscillation (for non-horizontal waves) */
+#define OSCIL_LEFT_YSHIFT 0.0   /* y-dependence of left oscillation (for non-horizontal waves) */
 /* The Courant number is given by c*DT/DX, where DT is the time step and DX the lattice spacing */
 /* The physical damping coefficient is given by GAMMA/(DT)^2 */
 /* Increasing COURANT speeds up the simulation, but decreases accuracy */
@@ -479,88 +510,146 @@ int reset_view = 0;         /* switch to reset 3D view parameters (for option RO
 
 #define B_COND 2
 
+#define PRECOMPUTE_BC 0     /* set to 1 to compute neighbours for Laplacian in advance */
+
 /* Parameters for length and speed of simulation */
 
-#define NSTEPS 2800       /* number of frames of movie */
-#define NVID 7            /* number of iterations between images displayed on screen */
-#define NSEG 1000         /* number of segments of boundary */
+#define NSTEPS 3000       /* number of frames of movie */
+#define NVID 4            /* number of iterations between images displayed on screen */
+#define NSEG 1000          /* number of segments of boundary */
 #define INITIAL_TIME 0      /* time after which to start saving frames */
 #define BOUNDARY_WIDTH 2    /* width of billiard boundary */
-#define PRINT_SPEED 0       /* print speed of moving source */
-#define PRINT_FREQUENCY 0       /* print frequency (for phased array) */
+#define PRINT_SPEED 0       /* set to 1 to print speed of moving source */
 
-#define PAUSE 200       /* number of frames after which to pause */
-#define PSLEEP 1         /* sleep time during pause */
+#define PAUSE 100       /* number of frames after which to pause */
+#define PSLEEP 3         /* sleep time during pause */
 #define SLEEP1  1        /* initial sleeping time */
 #define SLEEP2  1        /* final sleeping time */
-#define MID_FRAMES 20    /* number of still frames between parts of two-part movie */
-#define END_FRAMES 100    /* number of still frames at end of 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 */
+#define ROTATE_VIEW_WHILE_FADE 1    /* set to 1 to keep rotating viewpoint during fade */
 
 /* Parameters of initial condition */
 
-#define INITIAL_AMP 2.0            /* amplitude of initial condition */
-#define INITIAL_VARIANCE 0.000025    /* variance of initial condition */
-#define INITIAL_WAVELENGTH  0.05  /* wavelength of initial condition */
+#define INITIAL_AMP 1.0            /* amplitude of initial condition */
+#define INITIAL_VARIANCE 0.0005  /* variance of initial condition */
+#define INITIAL_WAVELENGTH  0.025  /* wavelength of initial condition */
 
 /* Plot type, see list in global_pdes.c  */
 
-#define PLOT 0
+#define ZPLOT 103     /* wave height */
+#define CPLOT 103     /* color scheme */
 
-#define PLOT_B 5        /* plot type for second movie */
+#define ZPLOT_B 108 
+#define CPLOT_B 108        /* plot type for second movie */
+
+#define CHANGE_LUMINOSITY 1     /* set to 1 to let luminosity depend on energy flux intensity */
+#define FLUX_WINDOW 30          /* size of averaging window of flux intensity */
+#define AMPLITUDE_HIGH_RES 1    /* set to 1 to increase resolution of plot */
+#define SHADE_3D 1              /* set to 1 to change luminosity according to normal vector */
+#define SHADE_WAVE 1            /* set to 1 to have luminosity depend on wave height */
+#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 0         /* set to 1 to draw boundary */
+#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.4          /* vertical scaling in energy plot */
+#define PLOT_SCALE_LOG_ENERGY 0.5       /* vertical scaling in log energy plot */
+
+/* 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 */
+
+#define ROTATE_VIEW 1       /* set to 1 to rotate position of observer */
+#define ROTATE_ANGLE 360.0   /* total angle of rotation during simulation */
 
 /* Color schemes */
 
-#define COLOR_PALETTE 17      /* Color palette, see list in global_pdes.c  */
-#define COLOR_PALETTE_B 13     /* Color palette, see list in global_pdes.c  */
+#define COLOR_PALETTE 11      /* Color palette, see list in global_pdes.c  */
+#define COLOR_PALETTE_B 16     /* Color palette, see list in global_pdes.c  */
 
 #define BLACK 1          /* background */
+#define COLOR_OUT_R 1.0    /* color outside domain */
+#define COLOR_OUT_G 1.0    
+#define COLOR_OUT_B 1.0    
 
 #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 PHASE_FACTOR 1.0       /* factor in computation of phase in color scheme P_3D_PHASE */
+#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 100.0     /* scaling factor for energy representation */
-#define LOG_SCALE 1.0     /* scaling factor for energy log representation */
-#define LOG_SHIFT 1.5     /* shift of colors on log scale */
-#define FLUX_SCALE 5.0e3    /* scaling factor for enegy flux represtnation */
+#define ATTENUATION 0.0    /* exponential attenuation coefficient of contrast with time */
+#define E_SCALE 100.0       /* scaling factor for energy representation */
+#define LOG_SCALE 0.75     /* scaling factor for energy log representation */
+#define LOG_SHIFT 0.5      /* shift of colors on log scale */
+#define LOG_ENERGY_FLOOR -10.0    /* floor value for log of (total) energy */
+#define LOG_MEAN_ENERGY_SHIFT 1.0   /* additional shift for log of mean energy */
+#define FLUX_SCALE 600.0    /* scaling factor for energy flux representation */
+#define FLUX_CSCALE 5.0      /* scaling factor for color in energy flux representation */
 #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 */
 #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 180.0    /* mean value of hue for color scheme C_HUE */
-#define HUEAMP -180.0      /* amplitude of variation of hue for color scheme C_HUE */
-
-#define DRAW_COLOR_SCHEME 1    /* set to 1 to plot the color scheme */
-#define COLORBAR_RANGE 2.0     /* scale of color scheme bar */
-#define COLORBAR_RANGE_B 0.8   /* scale of color scheme bar for 2nd part */
-#define ROTATE_COLOR_SCHEME 0   /* set to 1 to draw color scheme horizontally */
-#define CIRC_COLORBAR 0         /* set to 1 to draw circular color scheme */
-#define CIRC_COLORBAR_B 0       /* set to 1 to draw circular color scheme */
-
-#define SAVE_TIME_SERIES 0      /* set to 1 to save wave time series at a point */
+#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 NXMAZE 8      /* width of maze */
 #define NYMAZE 32      /* height of maze */
 #define MAZE_MAX_NGBH 5     /* max number of neighbours of maze cell */
-#define RAND_SHIFT 0        /* seed of random number generator */
+#define RAND_SHIFT 5        /* seed of random number generator */
 #define MAZE_XSHIFT 0.0     /* horizontal shift of maze */
 #define MAZE_WIDTH 0.02     /* half width of maze walls */
 
-/* for compatibility with sub_wave and sub_maze */
-#define ADD_POTENTIAL 0
-#define POT_MAZE 7
-#define POTENTIAL 0
+#define DRAW_COLOR_SCHEME 1       /* set to 1 to plot the color scheme */
+#define COLORBAR_RANGE 3.0      /* scale of color scheme bar */
+#define COLORBAR_RANGE_B 2.0    /* scale of color scheme bar for 2nd part */
+#define ROTATE_COLOR_SCHEME 0     /* set to 1 to draw color scheme horizontally */
+#define CIRC_COLORBAR 0         /* set to 1 to draw circular color scheme */
+#define CIRC_COLORBAR_B 0       /* set to 1 to draw circular color scheme */
+
+#define DRAW_WAVE_PROFILE 0     /* set to 1 to draw a profile of the wave */
+#define SAVE_TIME_SERIES 0      /* set to 1 to save wave time series at a point */
+
+#define ADD_POTENTIAL 0         /* set to 1 to add potential to z coordinate */
+#define POTENTIAL 10
+#define POT_FACT 20.0
 /* end of constants only used by sub_wave and sub_maze */
 
+
 /* For debugging purposes only */
-#define FLOOR 0         /* set to 1 to limit wave amplitude to VMAX */
+#define FLOOR 1         /* set to 1 to limit wave amplitude to VMAX */
 #define VMAX 10.0       /* max value of wave amplitude */
 
+/* Parameters controlling 3D projection */
+
+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};
+P_3D_ANGLE color scheme */
+double light[3] = {-0.816496581, -0.40824829, 0.40824829};      /* vector of "light" direction for P_3D_ANGLE color scheme */
+double observer[3] = {-10.0, 4.0, -3.5};    /* location of observer for REP_PROJ_3D representation */ 
+int reset_view = 0;         /* switch to reset 3D view parameters (for option ROTATE_VIEW) */
+
+#define RSCALE 0.01             /* scaling factor of radial component */
+#define RMAX 10.0               /* max value of radial component */
+#define Z_SCALING_FACTOR 0.85    /* 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 COS_VISIBLE -0.5        /* limit on cosine of normal to shown facets */
+
 ```
 

Parameters_July23.md

@@ -97,8 +97,6 @@ updated gradually.
 #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 140           /* number of grid point for grid of disks */
-// #define NGRIDY 70           /* number of grid point for grid of disks */
 #define NGRIDX 128           /* number of grid point for grid of disks */
 #define NGRIDY 64           /* number of grid point for grid of disks */
 #define EHRENFEST_RADIUS 0.9    /* radius of container for Ehrenfest urn configuration */
@@ -114,7 +112,6 @@ updated gradually.
 /* Parameters for length and speed of simulation */
 
 #define NSTEPS 3000      /* number of frames of movie */
-// #define NSTEPS 800      /* number of frames of movie */
 #define NVID 80         /* number of iterations between images displayed on screen */
 #define NSEG 250         /* number of segments of boundary */
 #define INITIAL_TIME 0     /* time after which to start saving frames */
@@ -175,8 +172,6 @@ updated gradually.
 
 /* particle properties */
 
-// #define ENERGY_HUE_MIN 330.0        /* color of original particle */
-// #define ENERGY_HUE_MAX 50.0         /* color of saturated particle */
 #define ENERGY_HUE_MIN 360.0        /* color of original particle */
 #define ENERGY_HUE_MAX 0.0         /* color of saturated particle */
 #define PARTICLE_HUE_MIN 359.0      /* color of original particle */
@@ -193,8 +188,6 @@ updated gradually.
 #define EQUILIBRIUM_DIST 3.5    /* Lennard-Jones equilibrium distance */
 #define EQUILIBRIUM_DIST_B 3.5  /* Lennard-Jones equilibrium distance for second type of particle */
 #define REPEL_RADIUS 15.0    /* radius in which repelling force acts (in units of particle radius) */
-// #define DAMPING 45.0          /* damping coefficient of particles */
-// #define DAMPING 35.0          /* damping coefficient of particles */
 #define DAMPING 18.0          /* damping coefficient of particles */
 #define OMEGA_INITIAL 10.0        /* initial angular velocity range */
 #define INITIAL_DAMPING 5.0  /* damping coefficient of particles during initial phase */
@@ -215,9 +208,7 @@ updated gradually.
 #define MU_XI 0.01           /* friction constant in thermostat */
 #define KSPRING_BOUNDARY 1.0e7    /* confining harmonic potential outside simulation region */
 #define KSPRING_OBSTACLE 1.0e11    /* harmonic potential of obstacles */
-// #define NBH_DIST_FACTOR 2.3       /* radius in which to count neighbours */
 #define NBH_DIST_FACTOR 2.7       /* radius in which to count neighbours */
-// #define NBH_DIST_FACTOR 4.0       /* radius in which to count neighbours */
 #define GRAVITY 5000.0             /* gravity acting on all particles */
 #define GRAVITY_X 0.0          /* horizontal gravity acting on all particles */
 #define INCREASE_GRAVITY 0     /* set to 1 to increase gravity during the simulation */

README.md

@@ -6,7 +6,7 @@ C code for videos on YouTube Channel https://www.youtube.com/c/NilsBerglund
 
 Parameter values used in specific simulations will be gradually added to file `Parameters.md`, `Parameters_June21.md` and so on.
 
-There are four groups of 6 files, 17 files, 5 files and 4 files. 
+There are four groups of 6 files, 19 files, 5 files and 4 files. 
 In addition the following files handling color schemes have been included:
 
 1. `hsluv.c`and `hsluv.h` from https://github.com/adammaj1/hsluv-color-gradient 
@@ -17,6 +17,12 @@ The following file (beta version) provides support for creating mazes:
 
 4. `sub_maze.c`
 
+The file 
+
+5. `Earth_Map_Blue_Marble_2002_large.ppm.gz` 
+
+is required by `wave_sphere.c` and should be unzipped before compiling. 
+
 ### Simulations of classical particles in billiards.
 
 1. *particle_billiard.c*:   simulation of a collection of non-interacting particles in a billiard
@@ -29,7 +35,7 @@ The following file (beta version) provides support for creating mazes:
 
 - Create subfolders `tif_part`, `tif_drop`
 - Customize constants at beginning of .c file
-- Compile with 
+- Compile with `make particle_billiard`, `make_drop_billiard`, etc, or
 
 `gcc -o particle_billiard particle_billiard.c-O3 -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lXmu -lglut`
 
@@ -49,26 +55,28 @@ in the shell before running the program
 
 1. *wave_billiard.c*:    simulation of the (linear) wave equation
 2. *wave_3d.c*:          3d rendering of wave equation
-3. *wave_comparison.c*: comparison of the wave equation in two different domains
-4. *wave_energy.c*:     a version of `wave_billiard` plotting the energy profile of the wave
-5. *mangrove.c*:        a version of `wave_billiard` with additional features to animate mangroves
-6. *heat.c*:            simulation of the heat equation, with optional drawing of gradient field lines
-7. *rde.c*:             simulation of reaction-diffusion equations, plots in 2D and 3D (including Schrödinger equation, 
-                        Euler equation, and shallow water equation)
-8. *schrodinger.c*:     simulation of the Schrodinger equation in 2D (old version)
-9. *global_pdes.c*:      global variables and parameters
-10. *global_3d.c*:        additional global variables for 3d version
-11. *sub_wave.c*:         drawing/computation routines common to `wave_billiard`, `heat` and `schrodinger`
-12. *sub_wave_comp.c*:    some modified functions needed by `wave_comparison`
-13. *sub_wave_3d.c*:      additional functions for 3d version
-14. *common_wave.c*:      common functions of `wave_billiard` and `wave_comparison`
-15. *colors_waves.c*:     colormaps used by wave simulations
-16. *sub_rde.c*:          additional routines for rde.c
-17. *sub_wave_rde_3d.c*:  additional 3d drawing routines for rde.c
+3. *wave_sphere.c*:      wave equation on a sphere (3D and 2D render)
+4. *wave_comparison.c*:  comparison of the wave equation in two different domains
+5. *wave_energy.c*:      a version of `wave_billiard` plotting the energy profile of the wave
+6. *mangrove.c*:         a version of `wave_billiard` with additional features to animate mangroves
+7. *heat.c*:             simulation of the heat equation, with optional drawing of gradient field lines
+8. *rde.c*:              simulation of reaction-diffusion equations, plots in 2D and 3D (including Schrödinger equation, 
+                         Euler equation, and shallow water equation)
+9. *schrodinger.c*:      simulation of the Schrodinger equation in 2D (old version)
+10. *global_pdes.c*:      global variables and parameters
+11. *global_3d.c*:        additional global variables for 3d version
+12. *sub_wave.c*:         drawing/computation routines common to `wave_billiard`, `heat` and `schrodinger`
+13. *sub_wave_comp.c*:    some modified functions needed by `wave_comparison`
+14. *sub_wave_3d.c*:      additional functions for 3d version
+15. *common_wave.c*:      common functions of `wave_billiard` and `wave_comparison`
+16. *colors_waves.c*:     colormaps used by wave simulations
+17. *sub_rde.c*:          additional routines for rde.c
+18. *sub_wave_rde_3d.c*:  additional 3d drawing routines for rde.c
+19. *sub_sphere.c*:       additional routines for wave_sphere.c
 
 - Create subfolders `tif_wave`, `tif_heat`, `tif_bz`, `tif_schrod`
 - Customize constants at beginning of .c file
-- Compile with 
+- Compile with `make wave_billiard`, etc, or
 
 `gcc -o wave_billiard wave_billiard.c -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lXmu -lglut -O3 -fopenmp`
 
@@ -99,7 +107,7 @@ in the shell before running the program
 
 - Create subfolder `tif_ljones`
 - Customize constants at beginning of .c file
-- Compile with 
+- Compile with `make lennardjones` or
 
 `gcc -o lennardjones lennardjones.c -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lXmu -lglut -O3 -fopenmp`
 
@@ -116,7 +124,7 @@ in the shell before running the program
 
 - Create subfolder `tif_perc`
 - Customize constants at beginning of .c file
-- Compile with 
+- Compile with `make percolation` or 
 
 `gcc -o percolation percolation.c -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lXmu -lglut -O3 -fopenmp`
 

global_3d.c

@@ -61,6 +61,11 @@
 #define SH_COAST 3      /* depth varying with x-coordinate */
 #define SH_COAST_MONOTONE 4      /* depth decreasing with x-coordinate */
 
+/* Type of rotating viewpoint */
+
+#define VP_HORIZONTAL 0     /* rotate in a horizontal plane */
+#define VP_ORBIT 1          /* rotate in a plane containing the origin */
+
 /* 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)||(cplot == Z_VORTICITY_ABS))
@@ -80,6 +85,7 @@
 
 #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))
 
+#define NMAXCIRC_SPHERE 100     /* max number of circles on sphere */
 
 int global_time = 0;
 double max_depth = 1.0;
@@ -137,3 +143,26 @@ typedef struct
     double gradx, grady;        /* gradient of water depth */
 } t_swater_depth;
 
+
+typedef struct
+{
+    double phi, theta;          /* phi, theta angles */
+    double cphi, sphi;          /* cos and sin of phi */
+    double ctheta, stheta, cottheta;   /* cos, sin and cotangent of theta */
+    double x, y, z;             /* x, y, z coordinates of point on sphere */
+    double radius;              /* radius with wave height */
+    double r, g, b;             /* RGB values for image */
+    short int indomain;         /* has value 1 if lattice point is in domain */
+    double x2d, y2d;            /* x and y coordinates for 2D representation */
+} t_wave_sphere;
+
+
+typedef struct 
+{
+    double phi, theta;          /* longitude, latitude */
+    double radius;              /* radius */
+    double x, y, z;             /* x, y, z coordinates of point on sphere */
+} t_circles_sphere;
+
+
+t_circles_sphere circ_sphere[NMAXCIRC_SPHERE];      /* circular scatterers on sphere */

global_ljones.c

@@ -45,6 +45,7 @@
 #define O_POOL_TABLE 3      /* obstacles around pockets of pool table */
 #define O_HLINE_HOLE_SPOKES 181    /* tips of spokes for S_HLINE_HOLE_SPOKES segment pattern */
 #define O_CIRCLE 4          /* one circle at the origin */
+#define O_FOUR_CIRCLES 5    /* four circles  */
 
 /* pattern of additional repelling segments */
 #define S_RECTANGLE 0       /* segments forming a rectangle */
@@ -72,6 +73,10 @@
 #define S_HLINE_HOLE_SPOKES 181    /* horizontal line with a hole in the bottom and extra spokes */
 #define S_EXT_CIRCLE_RECT 19    /* particles outside a circle and a rectangle */
 #define S_BIN_OPENING 20        /* bin containing particles opening at deactivation time */
+#define S_POLYGON_EXT 21        /* exterior of a regular polygon */
+#define S_WEDGE_EXT 22          /* exterior of a wedge */ 
+#define S_MIXER 23              /* exterior of a set of fan of rectangles */
+#define S_AIRFOIL 24            /* exterior of an air foil */
 
 /* particle interaction */
 
@@ -183,11 +188,19 @@
 #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_DIRECT_ENERGY 8 /* hue represents direction, luminosity represents energy */
 #define P_DIFF_NEIGHB 9   /* colors represent number of neighbours of different type */
 #define P_THERMOSTAT 10   /* colors show which particles are coupled to the thermostat */
 #define P_INITIAL_POS 11  /* colors depend on initial position of particle */
 #define P_NUMBER 12       /* colors depend on particle number */
+#define P_EMEAN 13        /* averaged kinetic energy (with exponential damping) */
+#define P_DIRECT_EMEAN 14 /* averaged version of P_DIRECT_ENERGY */
+#define P_NOPARTICLE 15   /* particles are not drawn (only the links between them) */
+
+/* Initial position dependence types */
+
+#define IP_X 0      /* color depends on x coordinate of initial position */
+#define IP_Y 1      /* color depends on y coordinate of initial position */
 
 /* Color schemes */
 
@@ -223,6 +236,7 @@ typedef struct
     double angle;               /* angle of particle's "spin" */
     short int active;           /* circle is active */
     double energy;              /* dissipated energy */
+    double emean;               /* mean energy */
     double vx;                  /* x velocity of particle */
     double vy;                  /* y velocity of particle */
     double omega;               /* angular velocity of particle */
@@ -231,6 +245,7 @@ typedef struct
     double fx;                  /* x component of force on particle */
     double fy;                  /* y component of force on particle */
     double torque;              /* torque on particle */
+    double dirmean;             /* time averaged direction */
     int close_to_boundary;      /* has value 1 if particle is close to a boundary */
     short int thermostat;       /* whether particle is coupled to thermostat */
     int hashcell;               /* hash cell in which particle is located */
@@ -331,5 +346,22 @@ typedef struct
     int color;                  /* color hue in case of different collisions */
 } t_collision;
 
+
+typedef struct
+{
+    int nactive;                /* number of active particles */
+    double beta;                /* inverse temperature */
+    double mean_energy;         /* mean energy */
+    double krepel;              /* force constant */
+    double xmincontainer, xmaxcontainer;      /* container size */
+    double fboundary;           /* boundary force */
+    double pressure;            /* pressure */
+    double gravity;             /* gravity */
+    double radius;              /* particle radius */
+    double angle;               /* orientation of obstacle */
+    double omega;               /* angular speed of obstacle */
+    double bdry_fx, bdry_fy;    /* components of boundary force */
+} t_lj_parameters;
+
 int ncircles, nobstacles, nsegments, ngroups = 1, counter = 0;
 

global_pdes.c

@@ -82,6 +82,14 @@
 #define D_TESLA 71              /* Tesla valve */
 #define D_TESLA_FOUR 72         /* four Tesla valves */
 
+/* for wave_sphere.c */
+
+#define D_LATITUDE 80           /* strip between two latitudes */
+#define D_SPHERE_CIRCLES 81     /* circles on the sphere */
+#define D_SPHERE_JULIA 82       /* Julia set on Riemann sphere */
+#define D_SPHERE_JULIA_INV 83   /* inverted Julia set on Riemann sphere */
+#define D_SPHERE_EARTH 84       /* map of the Earth */
+
 #define NMAXCIRCLES 10000       /* total number of circles/polygons (must be at least NCX*NCY for square grid) */
 #define NMAXPOLY 50000          /* maximal number of vertices of polygonal lines (for von Koch et al) */
 
@@ -104,6 +112,15 @@
 #define C_RINGS_T 201       /* obstacles arranged in concentric rings, triangular lattice */
 #define C_RINGS_SPIRAL 202  /* obstacles arranged on a "subflower" spiral, similar to C_GOLDEN_SPIRAL */
 
+#define C_HEX_BOTTOM 101    /* hex/triangular lattice in lower half */
+#define C_HEX_BOTTOM2 102   /* smaller hex/triangular lattice in lower half */
+#define C_SQUARE_BOTTOM 103 /* square lattice in lower half */
+
+#define C_SPH_DODECA 30     /* dodecahedron (on sphere) */
+#define C_SPH_ICOSA 31      /* icosahedron (on sphere) */
+#define C_SPH_OCTA 32       /* octahedron (on sphere) */
+#define C_SPH_CUBE 33       /* cube (on sphere) */
+
 #define C_ONE 97            /* one single circle, as for Sinai */
 #define C_TWO 98            /* two concentric circles of different type */
 #define C_NOTHING 99        /* no circle at all, for comparisons */

heat.c

@@ -211,7 +211,7 @@
 #define TWOSPEEDS 0          /* set to 1 to replace hardcore boundary by medium with different speed */
 #define WAVE_PACKET_SOURCE_TYPE 1       /* type of wave packet sources */
 #define N_WAVE_PACKETS 15               /* number of wave packets */
-
+#define OSCIL_LEFT_YSHIFT 0.0   /* y-dependence of left oscillation (for non-horizontal waves) */
 /* end of constants only used by sub_wave and sub_maze */
 
 #include "global_pdes.c"

lennardjones.c

@@ -62,8 +62,6 @@
 #define INITXMAX 2.0	/* x interval for initial condition */
 #define INITYMIN -1.1
 #define INITYMAX 1.1	/* y interval for initial condition */
-// #define INITYMAX 7.0	/* y interval for initial condition */
-// #define INITYMAX 9.0	/* y interval for initial condition */
 
 #define ADDXMIN -1.97
 #define ADDXMAX -0.8	/* x interval for adding particles */
@@ -84,11 +82,11 @@
 #define ADD_INITIAL_PARTICLES 0 /* set to 1 to add a second type of particles */
 #define CIRCLE_PATTERN_B 1  /* pattern of circles for additional particles */
 
-#define ADD_FIXED_OBSTACLES 1   /* set to 1 do add fixed circular obstacles */
+#define ADD_FIXED_OBSTACLES 0   /* set to 1 do add fixed circular obstacles */
 #define OBSTACLE_PATTERN 4  /* pattern of obstacles, see list in global_ljones.c */
 
-#define ADD_FIXED_SEGMENTS 0    /* set to 1 to add fixed segments as obstacles */
-#define SEGMENT_PATTERN 151     /* pattern of repelling segments, see list in global_ljones.c */
+#define ADD_FIXED_SEGMENTS 1    /* set to 1 to add fixed segments as obstacles */
+#define SEGMENT_PATTERN 24      /* pattern of repelling segments, see list in global_ljones.c */
 #define ROCKET_SHAPE 3        /* shape of rocket combustion chamber, see list in global_ljones.c */
 #define ROCKET_SHAPE_B 3      /* shape of second rocket */
 #define NOZZLE_SHAPE 6        /* shape of nozzle, see list in global_ljones.c */
@@ -101,8 +99,8 @@
 #define CENTER_PY 0         /* set to 1 to center vertical momentum */
 #define CENTER_PANGLE 0     /* set to 1 to center angular momentum */
 
-#define INTERACTION 2       /* particle interaction, see list in global_ljones.c */
-#define INTERACTION_B 2     /* particle interaction for second type of particle, see list in global_ljones.c */
+#define INTERACTION 1       /* particle interaction, see list in global_ljones.c */
+#define INTERACTION_B 1     /* particle interaction for second type of particle, see list in global_ljones.c */
 #define SPIN_INTER_FREQUENCY 4.0 /* angular frequency of spin-spin interaction */
 #define SPIN_INTER_FREQUENCY_B 4.0 /* angular frequency of spin-spin interaction for second particle type */
 
@@ -112,18 +110,21 @@
 #define PDISC_CANDIDATES 100 /* number of candidates in construction of Poisson disc process */
 #define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */
 
-#define LAMBDA 0.8	    /* parameter controlling the dimensions of domain */
-#define MU 0.015 	    /* parameter controlling radius of particles */
-#define MU_B 0.015          /* parameter controlling radius of particles of second type */
-#define NPOLY 25            /* number of sides of polygon */
+#define LAMBDA 0.75	    /* parameter controlling the dimensions of domain */
+#define MU 0.012 	    /* parameter controlling radius of particles */
+#define MU_B 0.010          /* parameter controlling radius of particles of second type */
+#define NPOLY 40            /* number of sides of polygon */
 #define APOLY 0.0           /* angle by which to turn polygon, in units of Pi/2 */ 
+#define AWEDGE 0.5          /* opening angle of wedge, 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 140           /* number of grid point for grid of disks */
-#define NGRIDY 70           /* number of grid point for grid of disks */
+// #define NGRIDX 140           /* number of grid point for grid of disks */
+// #define NGRIDY 70           /* number of grid point for grid of disks */
+#define NGRIDX 130           /* number of grid point for grid of disks */
+#define NGRIDY 65           /* number of grid point for grid of disks */
 #define EHRENFEST_RADIUS 0.9    /* radius of container for Ehrenfest urn configuration */
 #define EHRENFEST_WIDTH 0.035     /* width of tube for Ehrenfest urn configuration */
 #define TWO_CIRCLES_RADIUS_RATIO 0.8    /* ratio of radii for S_TWO_CIRCLES_EXT segment configuration */
@@ -136,10 +137,10 @@
 
 /* Parameters for length and speed of simulation */
 
-#define NSTEPS 2800      /* number of frames of movie */
-// #define NSTEPS 200      /* number of frames of movie */
-#define NVID 60         /* number of iterations between images displayed on screen */
-#define NSEG 250         /* number of segments of boundary */
+#define NSTEPS 4525      /* number of frames of movie */
+// #define NSTEPS 1000      /* number of frames of movie */
+#define NVID 30         /* number of iterations between images displayed on screen */
+#define NSEG 25         /* number of segments of boundary of circles */
 #define INITIAL_TIME 0     /* time after which to start saving frames */
 #define OBSTACLE_INITIAL_TIME 150     /* time after which to start moving obstacle */
 #define BOUNDARY_WIDTH 1    /* width of particle boundary */
@@ -159,7 +160,7 @@
 
 /* Plot type, see list in global_ljones.c  */
 
-#define PLOT 4
+#define PLOT 13
 #define PLOT_B 11        /* plot type for second movie */
 
 #define DRAW_BONDS 1    /* set to 1 to draw bonds between neighbours */
@@ -168,12 +169,17 @@
 #define ALTITUDE_LINES 0    /* set to 1 to add horizontal lines to show altitude */
 #define COLOR_SEG_GROUPS 1  /* set to 1 to collor segment groups differently */
 #define N_PARTICLE_COLORS 200   /* number of colors for P_NUMBER color scheme */
+#define INITIAL_POS_TYPE 0      /* type of initial position dependence */
+#define ERATIO 0.995          /* ratio for time-averagin in P_EMEAN color scheme */
+#define DRATIO 0.995          /* ratio for time-averagin in P_DIRECT_EMEAN color scheme */
 
 /* Color schemes */
 
-#define COLOR_PALETTE 18     /* Color palette, see list in global_ljones.c  */
+#define COLOR_PALETTE 10     /* Color palette, see list in global_ljones.c  */
+#define COLOR_PALETTE_EKIN 10   /* Color palette for kinetic energy */
 #define COLOR_PALETTE_ANGLE 10  /* Color palette for angle representation */
-#define COLOR_PALETTE_INITIAL_POS 10  /* Color palette for initial position representation */
+#define COLOR_PALETTE_DIRECTION 17  /* Color palette for direction representation */
+#define COLOR_PALETTE_INITIAL_POS 0  /* Color palette for initial position representation */
 
 #define BLACK 1          /* background */
 
@@ -189,9 +195,16 @@
 #define LUMAMP 0.3       /* amplitude of luminosity variation for scheme C_LUM */
 #define HUEMEAN 220.0    /* mean value of hue for color scheme C_HUE */
 #define HUEAMP -50.0      /* amplitude of variation of hue for color scheme C_HUE */
+#define COLOR_HUESHIFT 0.5     /* shift in color hue (for some cyclic palettes) */
 
-#define PRINT_PARAMETERS 0  /* set to 1 to print certain parameters */
+#define PRINT_PARAMETERS 1  /* set to 1 to print certain parameters */
 #define PRINT_TEMPERATURE 0 /* set to 1 to print current temperature */
+#define PRINT_ANGLE 1               /* set to 1 to print obstacle orientation */
+#define PRINT_OMEGA 1               /* set to 1 to print angular speed */
+#define PRINT_PARTICLE_SPEEDS 0     /* set to 1 to print average speeds/momenta of particles */
+#define PRINT_SEGMENTS_SPEEDS 0     /* set to 1 to print velocity of moving segments */
+#define PRINT_SEGMENTS_FORCE 1      /* set to 1 to print force on segments */
+#define FORCE_FACTOR 0.1            /* factor controlling length of force vector */
 
 /* particle properties */
 
@@ -199,9 +212,9 @@
 #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 6000.0         /* energy of particle with hottest color */
-#define HUE_TYPE0 60.0     /* hue of particles of type 0 */
-#define HUE_TYPE1 280.0      /* hue of particles of type 1 */
+#define PARTICLE_EMAX 5000.0        /* energy of particle with hottest color */
+#define HUE_TYPE0 280.0     /* hue of particles of type 0 */
+#define HUE_TYPE1 60.0      /* hue of particles of type 1 */
 #define HUE_TYPE2 140.0      /* hue of particles of type 2 */
 #define HUE_TYPE3 200.0     /* hue of particles of type 3 */
 
@@ -215,7 +228,7 @@
 #define OMEGA_INITIAL 10.0        /* initial angular velocity range */
 #define INITIAL_DAMPING 5.0  /* damping coefficient of particles during initial phase */
 #define DAMPING_ROT 100.0      /* dampint coefficient for rotation of particles */
-#define PARTICLE_MASS 2.5    /* mass of particle of radius MU */
+#define PARTICLE_MASS 1.0    /* mass of particle of radius MU */
 #define PARTICLE_MASS_B 1.0   /* mass of particle of radius MU */
 #define PARTICLE_INERTIA_MOMENT 0.2     /* moment of inertia of particle */
 #define PARTICLE_INERTIA_MOMENT_B 0.02     /* moment of inertia of second type of particle */
@@ -231,7 +244,9 @@
 #define MU_XI 0.01           /* friction constant in thermostat */
 #define KSPRING_BOUNDARY 1.0e7    /* confining harmonic potential outside simulation region */
 #define KSPRING_OBSTACLE 1.0e11    /* harmonic potential of obstacles */
-#define NBH_DIST_FACTOR 2.5       /* radius in which to count neighbours */
+// #define NBH_DIST_FACTOR 2.3       /* radius in which to count neighbours */
+// #define NBH_DIST_FACTOR 2.7       /* radius in which to count neighbours */
+#define NBH_DIST_FACTOR 3.8       /* radius in which to count neighbours */
 // #define NBH_DIST_FACTOR 4.0       /* radius in which to count neighbours */
 #define GRAVITY 0.0             /* gravity acting on all particles */
 #define GRAVITY_X 5000.0          /* horizontal gravity acting on all particles */
@@ -243,7 +258,7 @@
 #define KSPRING_VICSEK 0.2   /* spring constant for I_VICSEK_SPEED interaction */
 #define VICSEK_REPULSION 10.0    /* repulsion between particles in Vicsek model */
 
-#define ROTATION 1           /* set to 1 to include rotation of particles */
+#define ROTATION 0           /* set to 1 to include rotation of particles */
 #define COUPLE_ANGLE_TO_THERMOSTAT 1    /* set to 1 to couple angular degrees of freedom to thermostat */
 #define DIMENSION_FACTOR 1.0  /* scaling factor taking into account number of degrees of freedom */  
 #define KTORQUE 50.0         /* force constant in angular dynamics */
@@ -275,7 +290,7 @@
 #define CENTER_VIEW_ON_OBSTACLE 0   /* set to 1 to center display on moving obstacle */
 #define RESAMPLE_Y 0         /* set to 1 to resample y coordinate of moved particles (for shock waves) */
 #define NTRIALS 2000         /* number of trials when resampling */
-#define OBSTACLE_RADIUS 0.35  /* radius of obstacle for circle boundary conditions */
+#define OBSTACLE_RADIUS 0.25  /* radius of obstacle for circle boundary conditions */
 #define FUNNEL_WIDTH  0.25   /* funnel width for funnel boundary conditions */
 #define OBSTACLE_XMIN 0.0    /* initial position of obstacle */
 #define OBSTACLE_XMAX 3.0    /* final position of obstacle */
@@ -310,16 +325,13 @@
 #define TRACER_PARTICLE_MASS 4.0    /* relative mass of tracer particle */
 #define TRAJECTORY_WIDTH 3      /* width of tracer particle trajectory */
 
-#define ROTATE_BOUNDARY 0           /* set to 1 to rotate the repelling segments */
+#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 1000  /* period of rotating boundary */
-#define ROTATE_INITIAL_TIME 0       /* initial time without rotation */
-#define ROTATE_FINAL_TIME 100       /* final time without rotation */
-#define ROTATE_CHANGE_TIME 0.33     /* relative duration of acceleration/deceleration phases */
-#define OMEGAMAX 100.0              /* maximal rotation speed */
-#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 PRINT_SEGMENTS_SPEEDS 1     /* set to 1 to print velocity of moving segments */
+#define ROTATE_INITIAL_TIME 300       /* initial time without rotation */
+#define ROTATE_FINAL_TIME 300       /* final time without rotation */
+#define ROTATE_CHANGE_TIME 0.5     /* relative duration of acceleration/deceleration phases */
+#define OMEGAMAX  -2.0*PI              /* maximal rotation speed */
 
 #define MOVE_BOUNDARY 0        /* set to 1 to move repelling segments, due to force from particles */
 #define SEGMENTS_MASS 40.0     /* mass of collection of segments */
@@ -362,6 +374,9 @@
 #define DELTA_EKIN 2000.0        /* change of kinetic energy in reaction */
 #define COLLISION_TIME 15       /* time during which collisions are shown */
 
+#define CHANGE_RADIUS 0         /* set to 1 to change particle radius during simulation */
+#define MU_RATIO 0.666666667    /* ratio by which to increase radius */
+
 #define PRINT_PARTICLE_NUMBER 0     /* set to 1 to print total number of particles */
 #define PLOT_PARTICLE_NUMBER 0      /* set to 1 to make of plot of particle number over time */
 #define PARTICLE_NB_PLOT_FACTOR 0.5 /* expected final number of particles over initial number */
@@ -392,8 +407,8 @@
 #define FLOOR_OMEGA 0      /* set to 1 to limit particle momentum to PMAX */
 #define PMAX 1000.0        /* maximal force */
 
-#define HASHX 100   /* size of hashgrid in x direction */
-#define HASHY 50    /* size of hashgrid in y direction */
+#define HASHX 120   /* size of hashgrid in x direction */
+#define HASHY 60    /* 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 */
 
@@ -407,6 +422,8 @@
 #define NO_WRAP_BC ((BOUNDARY_COND != BC_PERIODIC)&&(BOUNDARY_COND != BC_PERIODIC_CIRCLE)&&(BOUNDARY_COND != BC_PERIODIC_TRIANGLE)&&(BOUNDARY_COND != BC_KLEIN)&&(BOUNDARY_COND != BC_PERIODIC_FUNNEL)&&(BOUNDARY_COND != BC_BOY)&&(BOUNDARY_COND != BC_GENUS_TWO))
 #define PERIODIC_BC ((BOUNDARY_COND == BC_PERIODIC)||(BOUNDARY_COND == BC_PERIODIC_CIRCLE)||(BOUNDARY_COND == BC_PERIODIC_FUNNEL)||(BOUNDARY_COND == BC_PERIODIC_TRIANGLE))
 #define TWO_OBSTACLES ((SEGMENT_PATTERN == S_TWO_CIRCLES_EXT)||(SEGMENT_PATTERN == S_TWO_ROCKETS))
+#define COMPUTE_EMEAN ((PLOT == P_EMEAN)||(PLOT_B == P_EMEAN)||(PLOT == P_DIRECT_EMEAN)||(PLOT_B == P_DIRECT_EMEAN))
+#define COMPUTE_DIRMEAN ((PLOT == P_DIRECT_EMEAN)||(PLOT_B == P_DIRECT_EMEAN))
 
 double xshift = 0.0;      /* x shift of shown window */
 double xspeed = 0.0;      /* x speed of obstacle */
@@ -422,6 +439,8 @@ double ysegments[2] = {SEGMENTS_Y0, SEGMENTS_Y0};  /* y coordinate of segments (
 double vxsegments[2] = {SEGMENTS_VX0, SEGMENTS_VX0};       /* vx coordinate of segments (for option MOVE_BOUNDARY) */
 double vysegments[2] = {SEGMENTS_VY0, SEGMENTS_VY0};       /* vy coordinate of segments (for option MOVE_BOUNDARY) */
 int thermostat_on = 1;    /* thermostat switch used when VARY_THERMOSTAT is on */
+double cosangle[NSEG+1];      
+double sinangle[NSEG+1];      /* precomputed trig functions of angles to draw circles faster */
 
 #define THERMOSTAT_ON ((THERMOSTAT)&&((!VARY_THERMOSTAT)||(thermostat_on)))
 
@@ -667,12 +686,17 @@ int thermostat_schedule(int i)
     else return(0);
 }
 
+double radius_schedule(int i)
+{
+    return(1.0 + (MU_RATIO - 1.0)*(double)i/(double)NSTEPS);
+}
+
 double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HASHX*HASHY], 
                         double qx[NMAXCIRCLES], double qy[NMAXCIRCLES], double qangle[NMAXCIRCLES],
                         double px[NMAXCIRCLES], double py[NMAXCIRCLES], double pangle[NMAXCIRCLES], 
                         double beta, int *nactive, int *nsuccess, int *nmove, int initial_phase)
 {
-    double a, totalenergy = 0.0, damping;  
+    double a, totalenergy = 0.0, damping, direction, dmean;  
     static double b = 0.25*SIGMA*SIGMA*DT_PARTICLE/MU_XI, xi = 0.0;
     int j, move;
     
@@ -691,6 +715,21 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA
         pangle[j] = particle[j].omega + 0.5*DT_PARTICLE*particle[j].torque;
         
         particle[j].energy = (px[j]*px[j] + py[j]*py[j])*particle[j].mass_inv;
+        
+        if (COMPUTE_EMEAN)
+            particle[j].emean = ERATIO*particle[j].emean + (1.0-ERATIO)*particle[j].energy;
+        
+        if (COMPUTE_DIRMEAN)
+        {
+            direction = argument(particle[j].vx, particle[j].vy);
+            dmean = particle[j].dirmean;
+            if (dmean < direction - PI) dmean += DPI;
+            else if (dmean > direction + PI) dmean -= DPI;
+            particle[j].dirmean = DRATIO*dmean + (1.0-DRATIO)*direction;
+            if (particle[j].dirmean < 0.0) particle[j].dirmean += DPI;
+            else if (particle[j].dirmean > DPI) particle[j].dirmean -= DPI;
+        }
+        
         if ((COUPLE_ANGLE_TO_THERMOSTAT)&&(particle[j].thermostat))
             particle[j].energy += pangle[j]*pangle[j]*particle[j].inertia_moment_inv;
         
@@ -1091,10 +1130,10 @@ void evolve_segment_groups(t_segment segment[NMAXSEGMENTS], int time, t_group_se
 void animation()
 {
     double time, scale, diss, rgb[3], dissip, gradient[2], x, y, dx, dy, dt, xleft, xright, 
-            a, b, length, fx, fy, force[2], totalenergy = 0.0, krepel = KREPEL, pos[2], prop, vx, beta = BETA, xi = 0.0, xmincontainer = BCXMIN, xmaxcontainer = BCXMAX, torque, torque_ij, fboundary = 0.0, pleft = 0.0, pright = 0.0, entropy[2], mean_energy, gravity = GRAVITY, speed_ratio, ymin, ymax, delta_energy, speed, ratio = 1.0, ratioc, cum_etot = 0.0, emean = 0.0;
+            a, b, length, fx, fy, force[2], totalenergy = 0.0, pos[2], prop, vx, xi = 0.0, torque, torque_ij, pleft = 0.0, pright = 0.0, entropy[2], speed_ratio, ymin, ymax, delta_energy, speed, ratio = 1.0, ratioc, cum_etot = 0.0, emean = 0.0, radius_ratio;
     double *qx, *qy, *px, *py, *qangle, *pangle, *pressure, *obstacle_speeds;
     int i, j, k, n, m, s, ij[2], i0, iplus, iminus, j0, jplus, jminus, p, q, p1, q1, p2, q2, total_neighbours = 0, 
-        min_nb, max_nb, close, wrapx = 0, wrapy = 0, nactive = 0, nadd_particle = 0, nmove = 0, nsuccess = 0, 
+        min_nb, max_nb, close, wrapx = 0, wrapy = 0, nadd_particle = 0, nmove = 0, nsuccess = 0, 
         tracer_n[N_TRACER_PARTICLES], traj_position = 0, traj_length = 0, move = 0, old, m0, floor, nthermo, wall = 0,
         group, gshift, n_total_active = 0, ncollisions = 0;
     int *particle_numbers;
@@ -1108,10 +1147,20 @@ void animation()
     t_group_data *group_speeds;
     t_collision *collisions;
     t_hashgrid *hashgrid;
+    t_lj_parameters params;
     char message[100];
 
     ratioc = 1.0 - ratio;
     
+    /* parameter values, grouped in a structure to simplify parameter printing */
+    params.beta = BETA;
+    params.krepel = KREPEL; 
+    params.xmincontainer = BCXMIN;
+    params.xmaxcontainer = BCXMAX; 
+    params.fboundary = 0.0;
+    params.gravity = GRAVITY;
+    params.radius = MU;
+    
     particle = (t_particle *)malloc(NMAXCIRCLES*sizeof(t_particle));    /* particles */  
     if (ADD_FIXED_OBSTACLES) obstacle = (t_obstacle *)malloc(NMAXOBSTACLES*sizeof(t_obstacle));    /* circular obstacles */  
     if (ADD_FIXED_SEGMENTS) 
@@ -1154,6 +1203,9 @@ void animation()
         group_speeds = (t_group_data *)malloc(ngroups*(INITIAL_TIME + NSTEPS)*sizeof(t_group_data));
     }
     
+    /* initialise array of trig functions to speed up drawing particles */
+    init_angles();
+    
     /* initialise positions and radii of circles */
     init_particle_config(particle);
     
@@ -1174,7 +1226,7 @@ void animation()
     
     printf("Initializing configuration\n");
 
-    nactive = initialize_configuration(particle, hashgrid, obstacle, px, py, pangle, tracer_n, segment);
+    params.nactive = initialize_configuration(particle, hashgrid, obstacle, px, py, pangle, tracer_n, segment);
      
 //     xi = 0.0;
     
@@ -1199,12 +1251,12 @@ void animation()
     {
         printf("Computing frame %d\n",i);
         
-        if (INCREASE_KREPEL) krepel = repel_schedule(i);
-        if (INCREASE_BETA) beta = temperature_schedule(i);        
+        if (INCREASE_KREPEL) params.krepel = repel_schedule(i);
+        if (INCREASE_BETA) params.beta = temperature_schedule(i);        
         if (DECREASE_CONTAINER_SIZE) 
         {
-            xmincontainer = container_size_schedule(i);
-            if (SYMMETRIC_DECREASE) xmaxcontainer = -container_size_schedule(i);
+            params.xmincontainer = container_size_schedule(i);
+            if (SYMMETRIC_DECREASE) params.xmaxcontainer = -container_size_schedule(i);
         }
         if ((ROTATE_BOUNDARY)&&(!SMOOTH_ROTATION)) rotate_segments(segment, rotation_schedule(i));
         if (VARY_THERMOSTAT) 
@@ -1235,7 +1287,7 @@ void animation()
 	
 	blank();
         
-        fboundary = 0.0;
+        params.fboundary = 0.0;
         pleft = 0.0;
         pright = 0.0;
         if (RECORD_PRESSURES) for (j=0; j<N_PRESSURES; j++) pressure[j] = 0.0;
@@ -1245,16 +1297,22 @@ void animation()
         {
             if (MOVE_OBSTACLE) 
             {
-                xmincontainer = obstacle_schedule_smooth(i, n);
-                xshift = xmincontainer;
+                params.xmincontainer = obstacle_schedule_smooth(i, n);
+                xshift = params.xmincontainer;
+            }
+            if ((ROTATE_BOUNDARY)&&(SMOOTH_ROTATION)) 
+                rotate_segments(segment, rotation_schedule_smooth(i,n));
+            if (ROTATE_BOUNDARY) 
+            {
+                params.omega = angular_speed;
+                params.angle = rotation_schedule_smooth(i,n);
             }
-            if ((ROTATE_BOUNDARY)&&(SMOOTH_ROTATION)) rotate_segments(segment, rotation_schedule_smooth(i,n));
             
-            if (INCREASE_GRAVITY) gravity = gravity_schedule(i,n); 
+            if (INCREASE_GRAVITY) params.gravity = gravity_schedule(i,n); 
             if ((BOUNDARY_COND == BC_RECTANGLE_WALL)&&(i < INITIAL_TIME + WALL_TIME)) wall = 1;
             else wall = 0;
             
-            if ((MOVE_BOUNDARY)||(MOVE_SEGMENT_GROUPS)) for (j=0; j<nsegments; j++)
+            if ((MOVE_BOUNDARY)||(MOVE_SEGMENT_GROUPS)||(PRINT_SEGMENTS_FORCE)) for (j=0; j<nsegments; j++)
             {
                 segment[j].fx = 0.0;
                 segment[j].fy = 0.0;
@@ -1272,10 +1330,10 @@ void animation()
                 particle[j].torque = 0.0;
                 
                 /* compute force from other particles */
-                compute_particle_force(j, krepel, particle, hashgrid);
+                compute_particle_force(j, params.krepel, particle, hashgrid);
                 
                 /* take care of boundary conditions */
-                fboundary += compute_boundary_force(j, particle, obstacle, segment, xmincontainer, xmaxcontainer, &pleft, &pright, pressure, wall);
+                params.fboundary += compute_boundary_force(j, particle, obstacle, segment, params.xmincontainer, params.xmaxcontainer, &pleft, &pright, pressure, wall);
 
                 /* align velocities in case of Vicsek models */
 //                 if (VICSEK_INT) 
@@ -1299,7 +1357,7 @@ void animation()
                 }
 
                 /* add gravity */
-                if (INCREASE_GRAVITY) particle[j].fy -= gravity/particle[j].mass_inv;
+                if (INCREASE_GRAVITY) particle[j].fy -= params.gravity/particle[j].mass_inv;
                 else 
                 {
                     particle[j].fy -= GRAVITY/particle[j].mass_inv;
@@ -1318,14 +1376,14 @@ void animation()
             }
             
             /* timestep of thermostat algorithm */
-            totalenergy = evolve_particles(particle, hashgrid, qx, qy, qangle, px, py, pangle, beta, &nactive, &nsuccess, &nmove, i < INITIAL_TIME);
+            totalenergy = evolve_particles(particle, hashgrid, qx, qy, qangle, px, py, pangle, params.beta, &params.nactive, &nsuccess, &nmove, i < INITIAL_TIME);
             
             
             /* evolution of lid coordinate */
-            if (BOUNDARY_COND == BC_RECTANGLE_LID) evolve_lid(fboundary);
+            if (BOUNDARY_COND == BC_RECTANGLE_LID) evolve_lid(params.fboundary);
             if (BOUNDARY_COND == BC_RECTANGLE_WALL) 
             {
-                if (i < INITIAL_TIME + WALL_TIME) evolve_wall(fboundary);
+                if (i < INITIAL_TIME + WALL_TIME) evolve_wall(params.fboundary);
                 else xwall = 0.0;
             }
             if ((MOVE_BOUNDARY)&&(i > OBSTACLE_INITIAL_TIME)) evolve_segments(segment, i);
@@ -1378,10 +1436,10 @@ void animation()
         if ((PARTIAL_THERMO_COUPLING)&&(i>N_T_AVERAGE)) 
         {
             nthermo = partial_thermostat_coupling(particle, xshift + PARTIAL_THERMO_SHIFT, segment);
-            printf("%i particles coupled to thermostat out of %i active\n", nthermo, nactive);
-            mean_energy = compute_mean_energy(particle);
+            printf("%i particles coupled to thermostat out of %i active\n", nthermo, params.nactive);
+            params.mean_energy = compute_mean_energy(particle);
         }
-        else mean_energy = totalenergy/(double)ncircles;
+        else params.mean_energy = totalenergy/(double)ncircles;
         
         if (CENTER_PX) center_momentum(px);
         if (CENTER_PY) center_momentum(py);
@@ -1434,9 +1492,9 @@ void animation()
             cum_etot += totalenergy;
         }
         
-        printf("Boundary force: %.3f\n", fboundary/(double)(ncircles*NVID)); 
+        printf("Boundary force: %.3f\n", params.fboundary/(double)(ncircles*NVID)); 
         if (RESAMPLE_Y) printf("%i succesful moves out of %i trials\n", nsuccess, nmove);
-        if (INCREASE_GRAVITY) printf("Gravity: %.3f\n", gravity);
+        if (INCREASE_GRAVITY) printf("Gravity: %.3f\n", params.gravity);
                 
         total_neighbours = 0;
         min_nb = 100;
@@ -1456,11 +1514,11 @@ void animation()
         if ((i > INITIAL_TIME)&&(REACTION_DIFFUSION)) 
         {
             ncollisions = update_types(particle, collisions, ncollisions, particle_numbers, i - INITIAL_TIME - 1,  &delta_energy);
-            if (EXOTHERMIC) beta *= 1.0/(1.0 + delta_energy/totalenergy);
-            nactive = 0;
+            if (EXOTHERMIC) params.beta *= 1.0/(1.0 + delta_energy/totalenergy);
+            params.nactive = 0;
             for (j=0; j<ncircles; j++) if (particle[j].active)
             {
-                nactive++;
+                params.nactive++;
                 qx[j] = particle[j].xc;
                 qy[j] = particle[j].yc;
                 px[j] = particle[j].vx;
@@ -1470,8 +1528,8 @@ void animation()
         }
  
         if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length);
-        draw_particles(particle, PLOT, beta, collisions, ncollisions);
-        draw_container(xmincontainer, xmaxcontainer, obstacle, segment, wall);
+        draw_particles(particle, PLOT, params.beta, collisions, ncollisions);
+        draw_container(params.xmincontainer, params.xmaxcontainer, obstacle, segment, wall);
 
         /* add a particle */
         if ((ADD_PARTICLES)&&(i > ADD_TIME)&&((i - INITIAL_TIME - ADD_TIME)%ADD_PERIOD == 1)&&(i < NSTEPS - FINAL_NOADD_PERIOD))
@@ -1480,11 +1538,22 @@ void animation()
             nadd_particle = add_particles(particle, px, py, nadd_particle);
         }
         
+        /* change particle radius */
+        if (CHANGE_RADIUS)
+        {
+            radius_ratio = radius_schedule(i+1)/radius_schedule(i);
+            printf("Particle radius factor %.5lg\t", radius_schedule(i+1));
+            for (j=0; j<ncircles; j++) particle[j].radius *= radius_ratio;
+            printf("Particle 0 radius %.5lg\n", particle[0].radius);
+            params.radius *= radius_ratio; 
+        }
+        
+        /* compute force on segments */
+        if (PRINT_SEGMENTS_FORCE) compute_segments_force(&params, segment);
+        
         update_hashgrid(particle, hashgrid, 1);
         
-        if (PRINT_PARAMETERS)
-            print_parameters(beta, mean_energy, krepel, xmaxcontainer - xmincontainer, 
-                         fboundary/(double)(ncircles*NVID), PRINT_LEFT, pressure, gravity);
+        if (PRINT_PARAMETERS) print_parameters(params, PRINT_LEFT, pressure, 1);
         if ((BOUNDARY_COND == BC_EHRENFEST)||(BOUNDARY_COND == BC_RECTANGLE_WALL)) 
             print_ehrenfest_parameters(particle, pleft, pright);
         else if (PRINT_PARTICLE_NUMBER) 
@@ -1504,7 +1573,6 @@ void animation()
         if (PLOT_SPEEDS) draw_speed_plot(group_speeds, i);
         if (PLOT_TRAJECTORIES) draw_trajectory_plot(group_speeds, i);
         
-        if (PRINT_OMEGA) print_omega(angular_speed); 
         else if (PRINT_PARTICLE_SPEEDS) print_particles_speeds(particle);
         else if (PRINT_SEGMENTS_SPEEDS)
         {
@@ -1546,11 +1614,9 @@ void animation()
             if ((i >= INITIAL_TIME)&&(DOUBLE_MOVIE))
             {
                 if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length);
-                draw_particles(particle, PLOT_B, beta, collisions, ncollisions);
-                draw_container(xmincontainer, xmaxcontainer, obstacle, segment, wall);
-                if (PRINT_PARAMETERS)
-                    print_parameters(beta, mean_energy, krepel, xmaxcontainer - xmincontainer, 
-                         fboundary/(double)(ncircles*NVID), PRINT_LEFT, pressure, gravity);
+                draw_particles(particle, PLOT_B, params.beta, collisions, ncollisions);
+                draw_container(params.xmincontainer, params.xmaxcontainer, obstacle, segment, wall);
+                if (PRINT_PARAMETERS) print_parameters(params, PRINT_LEFT, pressure, 0);
                 if (PLOT_SPEEDS) draw_speed_plot(group_speeds, i);
                 if (PLOT_TRAJECTORIES) draw_trajectory_plot(group_speeds, i);
                 if (BOUNDARY_COND == BC_EHRENFEST) print_ehrenfest_parameters(particle, pleft, pright);
@@ -1560,7 +1626,6 @@ void animation()
                         print_particle_types_number(particle, RD_TYPES);
                     else print_particle_number(ncircles);
                 }
-                if (PRINT_OMEGA) print_omega(angular_speed); 
                 else if (PRINT_PARTICLE_SPEEDS) print_particles_speeds(particle);
                 else if (PRINT_SEGMENTS_SPEEDS) print_segment_group_speeds(segment_group);
 //                     print_segments_speeds(vxsegments, vysegments);
@@ -1601,11 +1666,9 @@ void animation()
         {
             blank();
             if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length);
-            draw_particles(particle, PLOT, beta, collisions, ncollisions); 
-            draw_container(xmincontainer, xmaxcontainer, obstacle, segment, wall);
-            if (PRINT_PARAMETERS)
-                print_parameters(beta, mean_energy, krepel, xmaxcontainer - xmincontainer, 
-                    fboundary/(double)(ncircles*NVID), PRINT_LEFT, pressure, gravity);
+            draw_particles(particle, PLOT, params.beta, collisions, ncollisions); 
+            draw_container(params.xmincontainer, params.xmaxcontainer, obstacle, segment, wall);
+            if (PRINT_PARAMETERS) print_parameters(params, PRINT_LEFT, pressure, 1);
             if (PLOT_SPEEDS) draw_speed_plot(group_speeds, i);
             if (PLOT_TRAJECTORIES) draw_trajectory_plot(group_speeds, i);
             if (BOUNDARY_COND == BC_EHRENFEST) print_ehrenfest_parameters(particle, pleft, pright);
@@ -1615,7 +1678,6 @@ void animation()
                     print_particle_types_number(particle, RD_TYPES);
                 else print_particle_number(ncircles);
             }
-            if (PRINT_OMEGA) print_omega(angular_speed); 
             else if (PRINT_PARTICLE_SPEEDS) print_particles_speeds(particle);
             else if (PRINT_SEGMENTS_SPEEDS) print_segment_group_speeds(segment_group);
 //                 print_segments_speeds(vxsegments, vysegments);
@@ -1630,11 +1692,9 @@ void animation()
         if (DOUBLE_MOVIE) 
         {
             if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length);
-            draw_particles(particle, PLOT_B, beta, collisions, ncollisions);
-            draw_container(xmincontainer, xmaxcontainer, obstacle, segment, wall);
-            if (PRINT_PARAMETERS)
-                print_parameters(beta, mean_energy, krepel, xmaxcontainer - xmincontainer, 
-                    fboundary/(double)(ncircles*NVID), PRINT_LEFT, pressure, gravity);
+            draw_particles(particle, PLOT_B, params.beta, collisions, ncollisions);
+            draw_container(params.xmincontainer, params.xmaxcontainer, obstacle, segment, wall);
+            if (PRINT_PARAMETERS) print_parameters(params, PRINT_LEFT, pressure, 0);
             if (PLOT_SPEEDS) draw_speed_plot(group_speeds, i);
             if (PLOT_TRAJECTORIES) draw_trajectory_plot(group_speeds, i);
             if (BOUNDARY_COND == BC_EHRENFEST) print_ehrenfest_parameters(particle, pleft, pright);
@@ -1644,7 +1704,6 @@ void animation()
                     print_particle_types_number(particle, RD_TYPES);
                 else print_particle_number(ncircles);
             }
-            if (PRINT_OMEGA) print_omega(angular_speed); 
             else if (PRINT_PARTICLE_SPEEDS) print_particles_speeds(particle);
             else if (PRINT_SEGMENTS_SPEEDS) print_segment_group_speeds(segment_group);
 //                 print_segments_speeds(vxsegments, vysegments);
@@ -1660,45 +1719,45 @@ void animation()
         s = system("mv lj*.tif tif_ljones/");
     }
     
-    nactive = 0;
-    for (j=0; j<ncircles; j++) if (particle[j].active) nactive++;
-    printf("%i active particles\n", nactive);  
+    params.nactive = 0;
+    for (j=0; j<ncircles; j++) if (particle[j].active) params.nactive++;
+    printf("%i active particles\n", params.nactive);  
     
-    printf("1\n");
+//     printf("1\n");
     free(particle);
-    printf("2\n");
+//     printf("2\n");
     if (ADD_FIXED_OBSTACLES) free(obstacle);
-    printf("3\n");
+//     printf("3\n");
     if (ADD_FIXED_SEGMENTS) 
     {
         free(segment);
         free(segment_group);
     }
-    printf("4\n");
+//     printf("4\n");
     if (MOVE_SEGMENT_GROUPS) free(group_speeds);
-    printf("5\n");
+//     printf("5\n");
     if (TRACER_PARTICLE) free(trajectory);
-    printf("6\n");
+//     printf("6\n");
     if (PLOT_SPEEDS) free(obstacle_speeds);
-    printf("7\n");
+//     printf("7\n");
     if (PLOT_PARTICLE_NUMBER) free(particle_numbers);
-    printf("8\n");
+//     printf("8\n");
     free(hashgrid);
-    printf("9\n");
+//     printf("9\n");
     free(qx);
-    printf("10\n");
+//     printf("10\n");
     free(qy);
-    printf("11\n");
+//     printf("11\n");
     free(px);
-    printf("12\n");
+//     printf("12\n");
     free(py);
-    printf("13\n");
+//     printf("13\n");
     free(qangle);
-    printf("14\n");
+//     printf("14\n");
     free(pangle);
-    printf("15\n");
+//     printf("15\n");
     free(pressure);
-    printf("16\n");
+//     printf("16\n");
     if (REACTION_DIFFUSION) free(collisions);
     
     if (SAVE_TIME_SERIES) 

mangrove.c

@@ -259,6 +259,8 @@
 #define TWOSPEEDS 0          /* set to 1 to replace hardcore boundary by medium with different speed */
 #define WAVE_PACKET_SOURCE_TYPE 1       /* type of wave packet sources */
 #define N_WAVE_PACKETS 15               /* number of wave packets */
+#define OSCIL_LEFT_YSHIFT 0.0   /* y-dependence of left oscillation (for non-horizontal waves) */
+#define DRAW_WAVE_PROFILE 0     /* set to 1 to draw a profile of the wave */
 /* end of constants only used by sub_wave and sub_maze */
 
 #include "global_pdes.c"

rde.c

@@ -49,13 +49,13 @@
 #define WINWIDTH 	1920  /* window width */
 #define WINHEIGHT 	1150  /* window height */
 // #define NY 575          /* number of grid points on y axis */
-// #define NX 960          /* number of grid points on x axis */
-// #define NY 575          /* number of grid points on y axis */
+#define NX 960          /* number of grid points on x axis */
+#define NY 575          /* number of grid points on y axis */
 
 // #define WINWIDTH 	1280  /* window width */
 // #define WINHEIGHT 	720   /* window height */
-#define NX 640          /* number of grid points on x axis */
-#define NY 360          /* number of grid points on y axis */
+// #define NX 640          /* number of grid points on x axis */
+// #define NY 360          /* number of grid points on y axis */
 // #define NX 320          /* number of grid points on x axis */
 // #define NY 180          /* number of grid points on y axis */
 
@@ -76,16 +76,16 @@
 
 #define ADD_POTENTIAL 0 /* set to 1 to add a potential (for Schrodinger equation) */
 #define ADD_MAGNETIC_FIELD 0    /* set to 1 to add a magnetic field (for Schrodinger equation) - then set POTENTIAL 1 */
-#define ADD_FORCE_FIELD 0   /* set to 1 to add a foce field (for compressible Euler equation) */
+#define ADD_FORCE_FIELD 1   /* set to 1 to add a foce field (for compressible Euler equation) */
 #define POTENTIAL 7         /* type of potential or vector potential, see list in global_3d.c  */
 #define FORCE_FIELD 5       /* type of force field, see list in global_3d.c  */
 #define ADD_CORIOLIS_FORCE 0    /* set to 1 to add Coriolis force (quasigeostrophic Euler equations) */
-#define VARIABLE_DEPTH 1    /* set to 1 for variable depth in shallow water equation */
+#define VARIABLE_DEPTH 0    /* set to 1 for variable depth in shallow water equation */
 #define SWATER_DEPTH 4      /* variable depth in shallow water equation */
 
 #define ANTISYMMETRIZE_WAVE_FCT 0   /* set tot 1 to make wave function antisymmetric */
-#define ADAPT_STATE_TO_BC 0     /* to smoothly adapt initial state to obstacles */
-#define OBSTACLE_GEOMETRY 7     /* geometry of obstacles, as in B_DOMAIN */
+#define ADAPT_STATE_TO_BC 1     /* to smoothly adapt initial state to obstacles */
+#define OBSTACLE_GEOMETRY 1     /* geometry of obstacles, as in B_DOMAIN */
 // #define BC_STIFFNESS 100.0        /* controls region of boundary condition control */
 #define BC_STIFFNESS 50.0        /* controls region of boundary condition control */
 #define CHECK_INTEGRAL 1     /* set to 1 to check integral of first field */
@@ -102,7 +102,7 @@
 #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.5	    /* parameter controlling the dimensions of domain */
+#define LAMBDA 1.2	    /* parameter controlling the dimensions of domain */
 #define MU 0.06	            /* parameter controlling the dimensions of domain */
 #define NPOLY 5             /* number of sides of polygon */
 #define APOLY 2.0          /* angle by which to turn polygon, in units of Pi/2 */
@@ -129,7 +129,7 @@
 /* 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 */
+/* Physical parameters of wave equation */
 
 #define DT 0.00000025
 
@@ -164,6 +164,10 @@
 #define SMOOTHEN_PERIOD 10      /* period between smoothenings */
 #define SMOOTH_FACTOR 0.15       /* factor by which to smoothen */
 
+#define ADD_OSCILLATING_SOURCE 1        /* set to 1 to add an oscillating wave source */
+#define OSCILLATING_SOURCE_PERIOD 1     /* period of oscillating source */
+#define OSCILLATING_SOURCE_OMEGA 0.2    /* frequency of oscillating source */
+
 #define ADD_TRACERS 1    /* set to 1 to add tracer particles (for Euler equations) */
 #define N_TRACERS 1000    /* number of tracer particles */
 #define TRACERS_STEP 0.005  /* step size in tracer evolution */
@@ -201,7 +205,7 @@
 #define LAMINAR_FLOW_YPERIOD 1.0    /* period of laminar flow in y direction */
 #define PRESSURE_GRADIENT 0.2       /* amplitude of pressure gradient for Euler equation */
 
-#define SWATER_MIN_HEIGHT 2.0      /* min height of initial condition for shallow water equation */
+#define SWATER_MIN_HEIGHT 0.5      /* min height of initial condition for shallow water equation */
 // #define DEPTH_FACTOR 0.0125        /* proportion of min height in variable depth */
 // #define DEPTH_FACTOR 0.001        /* proportion of min height in variable depth */
 // #define DEPTH_FACTOR 0.0012        /* proportion of min height in variable depth */
@@ -216,20 +220,20 @@
 
 /* Boundary conditions, see list in global_pdes.c  */
 
-#define B_COND 1
+#define B_COND 0
 
 #define B_COND_LEFT 0
 #define B_COND_RIGHT 0
-#define B_COND_TOP 1
-#define B_COND_BOTTOM 1
+#define B_COND_TOP 0
+#define B_COND_BOTTOM 0
 
 
 /* Parameters for length and speed of simulation */
 
-#define NSTEPS 2100        /* number of frames of movie */
+#define NSTEPS 1300        /* number of frames of movie */
 // #define NSTEPS 500         /* number of frames of movie */
 // #define NVID 100           /* number of iterations between images displayed on screen */
-#define NVID 30            /* number of iterations between images displayed on screen */
+#define NVID 45            /* 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 */
@@ -249,7 +253,7 @@
 
 #define PLOT_3D 1    /* controls whether plot is 2D or 3D */
 
-#define ROTATE_VIEW 0       /* set to 1 to rotate position of observer */
+#define ROTATE_VIEW 1       /* set to 1 to rotate position of observer */
 #define ROTATE_ANGLE 360.0  /* total angle of rotation during simulation */
 
 #define DRAW_PERIODICISED 0     /* set to 1 to repeat wave periodically in x and y directions */
@@ -273,8 +277,8 @@
 #define DRAW_OUTSIDE_GRAY 0     /* experimental - draw outside of billiard in gray */
 #define ADD_POTENTIAL_TO_Z 0    /* set to 1 to add the external potential to z-coordinate of plot */
 #define ADD_POT_CONSTANT 0.35   /* constant added to wave height */
-#define ADD_HEIGHT_CONSTANT -2.0   /* constant added to wave height */
-#define DRAW_DEPTH 1            /* set to 1 to draw water depth */
+#define ADD_HEIGHT_CONSTANT -0.5   /* constant added to wave height */
+#define DRAW_DEPTH 0            /* set to 1 to draw water depth */
 #define DEPTH_SCALE 0.75         /* vertical scaling of depth plot */
 #define DEPTH_SHIFT -0.015      /* vertical shift of depth plot */
 
@@ -306,9 +310,9 @@
 
 /* Color schemes, see list in global_pdes.c  */
 
-#define COLOR_PALETTE 11       /* Color palette, see list in global_pdes.c  */
-#define COLOR_PALETTE_B 10      /* Color palette, see list in global_pdes.c  */
-// #define COLOR_PALETTE_B 0      /* Color palette, see list in global_pdes.c  */
+#define COLOR_PALETTE 14       /* Color palette, see list in global_pdes.c  */
+// #define COLOR_PALETTE_B 17      /* Color palette, see list in global_pdes.c  */
+#define COLOR_PALETTE_B 0      /* Color palette, see list in global_pdes.c  */
 
 #define BLACK 1          /* black background */
 
@@ -318,11 +322,11 @@
 
 #define SCALE 0          /* set to 1 to adjust color scheme to variance of field */
 #define SLOPE 1.0        /* sensitivity of color on wave amplitude */
-#define VSCALE_AMPLITUDE 10.0      /* additional scaling factor for color scheme P_3D_AMPLITUDE */
+#define VSCALE_AMPLITUDE 15.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 200.0       /* sensitivity of luminosity on module, for color scheme Z_ARGUMENT */
+#define SLOPE_SCHROD_LUM 400.0       /* sensitivity of luminosity on module, for color scheme Z_ARGUMENT */
 #define MIN_SCHROD_LUM 0.2       /* minimal luminosity in color scheme Z_ARGUMENT*/
 #define VSCALE_PRESSURE 2.0      /* additional scaling factor for color scheme Z_EULER_PRESSURE */
 #define PRESSURE_SHIFT 10.0        /* shift for color scheme Z_EULER_PRESSURE */
@@ -347,7 +351,7 @@
 #define VSCALE_VORTICITY 20.0     /* additional scaling factor for color scheme Z_EULERC_VORTICITY */
 #define VORTICITY_SHIFT 0.0     /* vertical shift of vorticity */
 #define ZSCALE_SPEED 0.3        /* additional scaling factor for z-coord Z_EULER_SPEED and Z_SWATER_SPEED */
-#define VSCALE_SWATER 300.0        /* additional scaling factor for color scheme Z_EULER_DENSITY */
+#define VSCALE_SWATER 250.0        /* additional scaling factor for color scheme Z_EULER_DENSITY */
 
 #define NXMAZE 7      /* width of maze */
 #define NYMAZE 7      /* height of maze */
@@ -390,6 +394,8 @@
 #define WAVE_PACKET_SOURCE_TYPE 1       /* type of wave packet sources */
 #define N_WAVE_PACKETS 15               /* number of wave packets */
 #define WAVE_PACKET_RADIUS 20            /* radius of wave packets */
+#define OSCIL_LEFT_YSHIFT 25.0   /* y-dependence of left oscillation (for non-horizontal waves) */
+#define DRAW_WAVE_PROFILE 1     /* set to 1 to draw a profile of the wave */
 /* end of constants added only for compatibility with wave_common.c */
 
 
@@ -400,8 +406,8 @@ double light[3] = {0.816496581, 0.40824829, 0.40824829};      /* vector of "ligh
 double observer[3] = {8.0, 8.0, 7.0};    /* location of observer for REP_PROJ_3D representation */ 
 int reset_view = 0;         /* switch to reset 3D view parameters (for option ROTATE_VIEW) */
 
-#define Z_SCALING_FACTOR 35.0   /* overall scaling factor of z axis for REP_PROJ_3D representation */
-#define XY_SCALING_FACTOR 1.7  /* overall scaling factor for on-screen (x,y) coordinates after projection */
+#define Z_SCALING_FACTOR 150.0   /* overall scaling factor of z axis for REP_PROJ_3D representation */
+#define XY_SCALING_FACTOR 2.5  /* overall scaling factor for on-screen (x,y) coordinates after projection */
 #define ZMAX_FACTOR 1.0        /* max value of z coordinate for REP_PROJ_3D representation */
 #define XSHIFT_3D 0.0          /* overall x shift for REP_PROJ_3D representation */
 #define YSHIFT_3D 0.1         /* overall y shift for REP_PROJ_3D representation */
@@ -1565,7 +1571,7 @@ void viewpoint_schedule(int i)
 void animation()
 {
     double time = 0.0, scale, dx, var, jangle, cosj, sinj, sqrintstep, 
-        intstep0, viscosity_printed, fade_value, noise = NOISE_INTENSITY, x, y, sign;
+        intstep0, viscosity_printed, fade_value, noise = NOISE_INTENSITY, x, y, sign, phase;
     double *phi[NFIELDS], *phi_tmp[NFIELDS], *potential_field, *vector_potential_field, *tracers, *gfield, *bc_field, *bc_field2;
     short int *xy_in;
     int i, j, k, s, nvid, field;
@@ -1672,9 +1678,11 @@ void animation()
 
 //     init_gaussian_wave(-1.0, 0.0, 0.005, 0.25, SWATER_MIN_HEIGHT, phi, xy_in);
     
-    init_linear_wave(-1.0, 0.01, 5.0e-8, 0.25, SWATER_MIN_HEIGHT, phi, xy_in);
+//     init_linear_wave(-1.0, 0.01, 5.0e-8, 0.25, SWATER_MIN_HEIGHT, phi, xy_in);
     
 //     init_linear_blob(0.0, -0.75, 0.025, 0.0, 0.001, 0.25, 0.5, SWATER_MIN_HEIGHT, phi, xy_in);
+        
+    init_elliptical_vibration(0.0, 0.0, 0.0075, LAMBDA, 1.0, 0.0003, 0.1, 1.0, SWATER_MIN_HEIGHT, phi, xy_in);
     
 //     add_gaussian_wave(-1.6, -0.5, 0.015, 0.25, SWATER_MIN_HEIGHT, phi, xy_in);
     
@@ -1759,6 +1767,13 @@ void animation()
 
         if (ADAPT_STATE_TO_BC) adapt_state_to_bc(phi, bc_field, xy_in);
         
+        if ((ADD_OSCILLATING_SOURCE)&&(i%OSCILLATING_SOURCE_PERIOD == 0))
+        {
+            phase = (double)i*OSCILLATING_SOURCE_OMEGA;
+            printf("Adding vibration with phase %.3lg\n", phase); 
+            add_elliptical_vibration(0.0, 0.0, 0.00007*cos(phase), LAMBDA, 1.0, 0.00004*cos(phase), 0.1, 1.0, SWATER_MIN_HEIGHT, phi, xy_in);
+        }
+                
         if (ADD_TRACERS)
         {
             printf("Evolving tracer particles\n");

schrodinger.c

@@ -191,6 +191,7 @@
 #define TWOSPEEDS 0          /* set to 1 to replace hardcore boundary by medium with different speed */
 #define WAVE_PACKET_SOURCE_TYPE 1       /* type of wave packet sources */
 #define N_WAVE_PACKETS 15               /* number of wave packets */
+#define OSCIL_LEFT_YSHIFT 0.0   /* y-dependence of left oscillation (for non-horizontal waves) */
 /* end of constants only used by sub_wave and sub_maze */
 
 

sub_rde.c

@@ -794,7 +794,93 @@ double init_linear_blob(double x, double y, double vx, double vy, double amp, do
 }
 
 
+double init_circular_vibration(double x, double y, double v, double radius, double amp, double var, double omega, double min, double *phi[NFIELDS], short int xy_in[NX*NY])
+/* initialise gaussian wave height for shallow water equation */
+{
+    int i, j;
+    double xy[2], r, angle, a, height;
+        
+    for (i=0; i<NX; i++)
+        for (j=0; j<NY; j++)
+        {
+            ij_to_xy(i, j, xy);
+            xy_in[i*NY+j] = xy_in_billiard(xy[0],xy[1]);
+            
+            if (xy_in[i*NY+j])
+            {
+                r = module2(xy[0]-x, xy[1]-y);
+                angle = argument(xy[0]-x, xy[1]-y);
+                height = exp(-(r - radius)*(r - radius)/var)*cos(omega*angle);
+                phi[0][i*NY+j] = min + amp*height;
+                phi[1][i*NY+j] = v*cos(angle)*height;
+                phi[2][i*NY+j] = v*sin(angle)*height;
+            }
+            else
+            {
+                phi[0][i*NY+j] = 0.0;
+                phi[1][i*NY+j] = 0.0;
+                phi[2][i*NY+j] = 0.0;
+            }
+        }
+}
 
+
+double init_elliptical_vibration(double x, double y, double v, double radiusx, double radiusy, double amp, double var, double omega, double min, double *phi[NFIELDS], short int xy_in[NX*NY])
+/* initialise gaussian wave height for shallow water equation */
+{
+    int i, j;
+    double xy[2], r, angle, a, height;
+        
+    for (i=0; i<NX; i++)
+        for (j=0; j<NY; j++)
+        {
+            ij_to_xy(i, j, xy);
+            xy_in[i*NY+j] = xy_in_billiard(xy[0],xy[1]);
+            
+            if (xy_in[i*NY+j])
+            {
+                r = module2((xy[0]-x)/radiusx, (xy[1]-y)/radiusy);
+                angle = argument((xy[0]-x)/radiusx, (xy[1]-y)/radiusy);
+                height = exp(-(r - 1.0)*(r - 1.0)/var)*cos(omega*angle);
+                phi[0][i*NY+j] = min + amp*height;
+                phi[1][i*NY+j] = v*cos(angle)*height*radiusx;
+                phi[2][i*NY+j] = v*sin(angle)*height*radiusy;
+            }
+            else
+            {
+                phi[0][i*NY+j] = 0.0;
+                phi[1][i*NY+j] = 0.0;
+                phi[2][i*NY+j] = 0.0;
+            }
+        }
+}
+
+double add_elliptical_vibration(double x, double y, double v, double radiusx, double radiusy, double amp, double var, double omega, double min, double *phi[NFIELDS], short int xy_in[NX*NY])
+/* initialise gaussian wave height for shallow water equation */
+{
+    int i, j;
+    double xy[2], r, angle, a, height;
+        
+    for (i=0; i<NX; i++)
+        for (j=0; j<NY; j++)
+        {
+            ij_to_xy(i, j, xy);
+            xy_in[i*NY+j] = xy_in_billiard(xy[0],xy[1]);
+            
+            if (xy_in[i*NY+j])
+            {
+                r = module2((xy[0]-x)/radiusx, (xy[1]-y)/radiusy);
+                if (r < 1.0 + var)
+                {
+                    angle = argument((xy[0]-x)/radiusx, (xy[1]-y)/radiusy);
+                    height = exp(-(r - 1.0)*(r - 1.0)/var)*cos(omega*angle);
+                    phi[0][i*NY+j] += amp*height;
+                    phi[1][i*NY+j] += v*cos(angle)*height*radiusx;
+                    phi[2][i*NY+j] += v*sin(angle)*height*radiusy;
+                }
+            }
+        }
+}
 double add_gaussian_wave(double x, double y, double amp, double radius, double min, double *phi[NFIELDS], short int xy_in[NX*NY])
 /* initialise gaussian wave height for shallow water equation */
 {

sub_sphere.c

@@ -0,0 +1,895 @@
+void init_sphere_2D()		/* initialisation of window */
+{
+    glLineWidth(3);
+
+    glClearColor(0.0, 0.0, 0.0, 1.0);
+    glClear(GL_COLOR_BUFFER_BIT);
+
+    glOrtho(0.0, NX, DPOLE, NY-DPOLE, -1.0, 1.0);
+}
+
+void draw_vertex_sphere(double xyz[3])
+{
+    double xy_screen[2];
+    
+    xyz_to_xy(xyz[0], xyz[1], xyz[2], xy_screen);
+    glVertex2d(xy_screen[0], xy_screen[1]);
+}
+
+void init_wave_sphere(t_wave_sphere wsphere[NX*NY])
+/* initialize sphere data */
+{
+    int i, j;
+    double dphi, dtheta, theta0, xy[2], phishift;
+    
+    printf("Initializing wsphere\n");
+    
+    dphi = DPI/(double)NX;
+//     dtheta = PI/(double)NY;
+    dtheta = PI/(double)(NY-2*(DPOLE));
+    theta0 = (double)(DPOLE)*dtheta;
+    phishift = PHISHIFT*(XMAX-XMIN)/360.0;
+    
+    #pragma omp parallel for private(i,j)
+    for (i=0; i<NX; i++)
+    {
+        for (j=DPOLE; j<NY-DPOLE; j++)
+            wsphere[i*NY+j].theta = (double)j*dtheta - theta0;
+        for (j=0; j<DPOLE; j++) wsphere[i*NY+j].theta = 0.0;
+        for (j=NY-DPOLE; j<NY; j++) wsphere[i*NY+j].theta = PI;
+        
+        for (j=0; j<NY; j++)
+        {
+            wsphere[i*NY+j].phi = (double)i*dphi;
+//             wsphere[i*NY+j].theta = (double)j*dtheta;
+            
+            wsphere[i*NY+j].cphi = cos(wsphere[i*NY+j].phi);
+            wsphere[i*NY+j].sphi = sin(wsphere[i*NY+j].phi);
+            
+            wsphere[i*NY+j].ctheta = cos(wsphere[i*NY+j].theta);
+            wsphere[i*NY+j].stheta = sin(wsphere[i*NY+j].theta);
+            
+            wsphere[i*NY+j].x = wsphere[i*NY+j].cphi*wsphere[i*NY+j].stheta;
+            wsphere[i*NY+j].y = wsphere[i*NY+j].sphi*wsphere[i*NY+j].stheta;
+            wsphere[i*NY+j].z = -wsphere[i*NY+j].ctheta;
+            
+            ij_to_xy(NX-1-i,j,xy);
+//             xy[0] = XMIN + ((double)(NX-i-1))*(XMAX-XMIN)/((double)NX);
+//             xy[1] = YMIN + ((double)(j-DPOLE))*(YMAX-YMIN)/((double)(NY-2*DPOLE));
+            
+            xy[0] += phishift;
+            if (xy[0] > XMAX) xy[0] += XMIN - XMAX;
+            
+            xy[1] *= (double)NY/(double)(NY-2*DPOLE);
+            
+            wsphere[i*NY+j].x2d = xy[0]; 
+            wsphere[i*NY+j].y2d = xy[1]; 
+        }
+        
+        /* cotangent, taking care of not dividing by zero */
+        for (j=1; j<NY-1; j++) wsphere[i*NY+j].cottheta = wsphere[i*NY+j].ctheta/wsphere[i*NY+j].stheta;
+        wsphere[i*NY].cottheta = wsphere[i*NY+1].cottheta;
+        wsphere[i*NY + NY-1].cottheta = wsphere[i*NY+1].cottheta;
+    }
+}
+
+
+void init_earth_map(t_wave_sphere wsphere[NX*NY])
+/* init file from earth map */
+{
+    int i, j, ii, jj, k, nx, ny, maxrgb, nmaxpixels = 4915200, scan, rgbval, diff, sshift, nshift;
+    int *rgb_values;
+    double cratio, rx, ry, cy;
+    FILE *image_file;
+    
+    printf("Reading Earth map\n");
+    
+    rgb_values = (int *)malloc(3*nmaxpixels*sizeof(int));
+    
+    image_file = fopen("Earth_Map_Blue_Marble_2002_large.ppm", "r");
+    scan = fscanf(image_file,"%i %i\n", &nx, &ny);
+    scan = fscanf(image_file,"%i\n", &maxrgb);
+
+    if (nx*ny > nmaxpixels)
+    {
+        printf("Image too large, increase nmaxpixels in choose_colors()\n");
+        exit(0);
+    }
+    
+    /* shift due to min/max latitudes of image */
+    sshift = 0 + DPOLE;
+    nshift = 0 + DPOLE; 
+    
+    /* read rgb values */
+    for (j=0; j<ny; j++)
+        for (i=0; i<nx; i++)
+            for (k=0; k<3; k++)
+                {
+                    scan = fscanf(image_file,"%i\n", &rgbval);
+                    rgb_values[3*(j*nx+i)+k] = rgbval;
+                }
+                    
+    cratio = 1.0/(double)maxrgb;
+    rx = (double)nx/(double)NX;
+    ry = (double)ny/(double)(NY - sshift - nshift);
+//     cy = rx*(double)(NY - nshift);
+    
+    /* build wave table */
+    for (i=0; i<NX; i++)
+        for (j=0; j<NY; j++)
+        {
+            ii = (int)(rx*(double)(NX-1 - i)) + nx/2;
+            if (ii > nx-1) ii -= nx;
+//             jj = (int)(-ry*(double)j + cy);
+//             jj = (int)(ry*(double)(NY-nshift - j)) + sshift;
+            jj = (int)(ry*(double)(NY-nshift - j));
+            if (jj > ny-1) jj = ny-1;
+            if (jj < 0) jj = 0;
+            wsphere[i*NY+j].r = (double)rgb_values[3*(jj*nx+ii)]*cratio;
+            wsphere[i*NY+j].g = (double)rgb_values[3*(jj*nx+ii)+1]*cratio;
+            wsphere[i*NY+j].b = (double)rgb_values[3*(jj*nx+ii)+2]*cratio;
+            
+            /* decide which points are in the Sea */
+            diff = iabs(rgb_values[3*(jj*nx+ii)] - 10);
+            diff += iabs(rgb_values[3*(jj*nx+ii)+1] - 10);
+            diff += iabs(rgb_values[3*(jj*nx+ii)+2] - 51);
+            wsphere[i*NY+j].indomain = (diff < 10);
+        }
+    
+    free(rgb_values);
+    fclose(image_file);
+}
+
+
+int ij_to_sphere(int i, int j, double r, t_wave_sphere wsphere[NX*NY], double xyz[3])
+/* convert spherical to rectangular coordinates */
+{
+    double pscal, newr;
+    static double norm_observer;
+    static int first = 1;
+    
+    if (first)
+    {
+        norm_observer = sqrt(observer[0]*observer[0] + observer[1]*observer[1] + observer[2]*observer[2]);
+        first = 0;
+    }
+    
+    newr = wsphere[i*NY+j].radius;
+    
+    xyz[0] = wsphere[i*NY+j].x;
+    xyz[1] = wsphere[i*NY+j].y;
+    xyz[2] = wsphere[i*NY+j].z;
+    
+    pscal = xyz[0]*observer[0] + xyz[1]*observer[1] + xyz[2]*observer[2];
+    
+    xyz[0] *= newr;
+    xyz[1] *= newr;
+    xyz[2] *= newr;
+    
+    return(pscal/norm_observer > COS_VISIBLE);
+}
+
+int init_circle_sphere(t_circles_sphere *circles, int circle_pattern)
+/* initialise circles on sphere */
+/* for billiard shape D_SPHERE_CIRCLES */
+{
+    int ncircles, k;
+    double alpha, beta, gamma;
+    
+    switch (circle_pattern) {
+        case (C_SPH_DODECA):
+        {
+            alpha = acos(sqrt(5.0)/3.0);
+            beta = acos(1.0/3.0);
+            gamma = asin(sqrt(3.0/8.0));
+            
+            circles[0].theta = 0.0;
+            circles[0].phi = 0.0;
+            
+            for (k=0; k<3; k++)
+            {
+                circles[1+k].theta = alpha;
+                circles[1+k].phi = (double)k*DPI/3.0;
+            }
+
+            for (k=0; k<3; k++)
+            {
+                circles[4+k].theta = beta;
+                circles[4+k].phi = (double)k*DPI/3.0 + gamma;
+                circles[7+k].theta = beta;
+                circles[7+k].phi = (double)k*DPI/3.0 - gamma;
+            }
+            
+            for (k=0; k<3; k++)
+            {
+                circles[10+k].theta = PI - beta;
+                circles[10+k].phi = (double)k*DPI/3.0 + PI/3.0 + gamma;
+                circles[13+k].theta = PI - beta;
+                circles[13+k].phi = (double)k*DPI/3.0 + PI/3.0 - gamma;
+            }
+
+            for (k=0; k<3; k++)
+            {
+                circles[16+k].theta = PI - alpha;
+                circles[16+k].phi = (double)k*DPI/3.0 + PI/3.0;
+            }
+            
+            circles[19].theta = PI;
+            circles[19].phi = 0.0;
+
+            for (k=0; k<20; k++)
+            {
+                circles[k].radius = MU;
+                circles[k].x = sin(circles[k].theta)*cos(circles[k].phi);
+                circles[k].y = sin(circles[k].theta)*sin(circles[k].phi);
+                circles[k].z = cos(circles[k].theta);
+            }
+            
+            ncircles = 20;
+            break;
+        }
+        case (C_SPH_ICOSA):
+        {
+            alpha = acos(1.0/sqrt(5.0));
+            
+            circles[0].theta = 0.0;
+            circles[0].phi = 0.0;
+       
+            for (k=0; k<5; k++)
+            {
+                circles[1+k].theta = alpha;
+                circles[1+k].phi = (double)k*DPI/5.0;
+            }
+            
+            for (k=0; k<5; k++)
+            {
+                circles[6+k].theta = PI - alpha;
+                circles[6+k].phi = (double)k*DPI/5.0 + PI/5.0;
+            }
+
+            circles[11].theta = PI;
+            circles[11].phi = 0.0;
+
+            for (k=0; k<12; k++) 
+            {
+                circles[k].radius = MU;
+                circles[k].x = sin(circles[k].theta)*cos(circles[k].phi);
+                circles[k].y = sin(circles[k].theta)*sin(circles[k].phi);
+                circles[k].z = cos(circles[k].theta);
+            }
+            
+            ncircles = 12;
+            break;
+        }
+        case (C_SPH_OCTA):
+        {
+            circles[0].theta = 0.0;
+            circles[0].phi = 0.0;
+       
+            for (k=0; k<4; k++)
+            {
+                circles[1+k].theta = PID;
+                circles[1+k].phi = (double)k*PID;
+            }
+            
+            circles[5].theta = PI;
+            circles[5].phi = 0.0;
+
+            for (k=0; k<6; k++) 
+            {
+                circles[k].radius = MU;
+                circles[k].x = sin(circles[k].theta)*cos(circles[k].phi);
+                circles[k].y = sin(circles[k].theta)*sin(circles[k].phi);
+                circles[k].z = cos(circles[k].theta);
+            }
+            
+            ncircles = 6;
+            break;
+        }
+        case (C_SPH_CUBE):
+        {
+            alpha = acos(1.0/3.0);
+            
+            circles[0].theta = 0.0;
+            circles[0].phi = 0.0;
+       
+            for (k=0; k<3; k++)
+            {
+                circles[1+k].theta = alpha;
+                circles[1+k].phi = (double)k*DPI/3.0;
+                circles[4+k].theta = PI - alpha;
+                circles[4+k].phi = (double)k*DPI/3.0 + PI/3.0;
+            }
+            
+            circles[7].theta = PI;
+            circles[7].phi = 0.0;
+
+            for (k=0; k<8; k++) 
+            {
+                circles[k].radius = MU;
+                circles[k].x = sin(circles[k].theta)*cos(circles[k].phi);
+                circles[k].y = sin(circles[k].theta)*sin(circles[k].phi);
+                circles[k].z = cos(circles[k].theta);
+            }
+            
+            ncircles = 8;
+            break;
+        }
+        default: 
+        {
+            printf("Function init_circle_sphere not defined for this pattern \n");
+        }
+    }
+    return(ncircles);
+}
+
+int xy_in_billiard_sphere(int i, int j, t_wave_sphere wsphere[NX*NY])
+/* returns 1 if (x,y) represents a point in the billiard */
+{
+    int k;
+    double pscal, dist, r, u, v, u1, v1;
+    static double cos_rot, sin_rot;
+    static int first = 1;
+    
+    if (first)
+    {
+        if (B_DOMAIN == D_SPHERE_EARTH) init_earth_map(wsphere);
+        else if ((B_DOMAIN == D_SPHERE_JULIA)||(B_DOMAIN == D_SPHERE_JULIA_INV))
+        {
+            cos_rot = cos(JULIA_ROT*DPI/360.0);
+            sin_rot = sin(JULIA_ROT*DPI/360.0);
+        }
+        first = 0;
+    }
+    
+    switch (B_DOMAIN) {
+        case (D_NOTHING):
+        {
+            return(1);
+        }
+        case (D_LATITUDE):
+        {
+            return(vabs(wsphere[i*NY+j].theta - PID) < LAMBDA*PID);
+        }
+        case (D_SPHERE_CIRCLES):
+        {
+            for (k=0; k<ncircles; k++)
+            {
+                pscal = cos(wsphere[i*NY+j].phi - circ_sphere[k].phi)*wsphere[i*NY+j].stheta*sin(circ_sphere[k].theta);
+                pscal += wsphere[i*NY+j].ctheta*cos(circ_sphere[k].theta);
+            
+                dist = acos(pscal);
+                if (dist < circ_sphere[k].radius) return(0);
+            }
+            return(1);
+        }
+        case (D_SPHERE_JULIA):
+        {
+            if (wsphere[i*NY+j].z == 1.0) return(1);
+            if (wsphere[i*NY+j].z == -1.0) return(0);
+            r = (1.0 + wsphere[i*NY+j].ctheta)/wsphere[i*NY+j].stheta;
+            u1 = r*wsphere[i*NY+j].cphi*JULIA_SCALE;
+            v1 = r*wsphere[i*NY+j].sphi*JULIA_SCALE;
+            u = u1*cos_rot + v1*sin_rot;
+            v = -u1*sin_rot + v1*cos_rot;
+            i = 0;
+            while ((i<MANDELLEVEL)&&(u*u+v*v < 1000.0*MANDELLIMIT))
+            {
+                u1 = u*u - v*v + JULIA_RE;
+                v = 2.0*u*v + JULIA_IM;
+                u = u1;
+                i++;
+            }
+            if (u*u + v*v < MANDELLIMIT) return(1);
+            return(0);
+        }
+        case (D_SPHERE_JULIA_INV):
+        {
+            if (wsphere[i*NY+j].z == 1.0) return(1);
+            if (wsphere[i*NY+j].z == -1.0) return(0);
+            r = (1.0 - wsphere[i*NY+j].ctheta)/wsphere[i*NY+j].stheta;
+            u1 = r*wsphere[i*NY+j].cphi*JULIA_SCALE;
+            v1 = r*wsphere[i*NY+j].sphi*JULIA_SCALE;
+            u = u1*cos_rot + v1*sin_rot;
+            v = -u1*sin_rot + v1*cos_rot;
+            i = 0;
+            while ((i<MANDELLEVEL)&&(u*u+v*v < 1000.0*MANDELLIMIT))
+            {
+                u1 = u*u - v*v + JULIA_RE;
+                v = 2.0*u*v + JULIA_IM;
+                u = u1;
+                i++;
+            }
+            if (u*u + v*v < MANDELLIMIT) return(0);
+            return(1);
+        }
+        case (D_SPHERE_EARTH):
+        {
+            return(wsphere[i*NY+j].indomain);
+        }
+        default:
+        {
+            printf("Function xy_in_billiard_sphere not defined for this billiard \n");
+            return(0);
+        }
+    }
+}
+
+void init_wave_flat_sphere(double phi[NX*NY], double psi[NX*NY], short int xy_in[NX*NY], t_wave_sphere wsphere[NX*NY])
+/* initialise field with drop - phi is wave height, psi is phi at time t-1 */
+/* phi0 is longidude, theta0 is latitude */
+{
+    int i, j;
+    
+    printf("Initializing wave\n"); 
+//     #pragma omp parallel for private(i,j,xy,dist2)
+    for (i=0; i<NX; i++)
+    {
+        if (i%100 == 0) printf("Initializing column %i of %i\n", i, NX);
+        for (j=0; j<NY; j++)
+        {
+            xy_in[i*NY+j] = xy_in_billiard_sphere(i, j, wsphere);
+            
+	    phi[i*NY+j] = 0.0;
+            psi[i*NY+j] = 0.0;
+        }
+    }
+}
+
+
+void init_circular_wave_sphere(double phi0, double theta0, double phi[NX*NY], double psi[NX*NY], short int xy_in[NX*NY], t_wave_sphere wsphere[NX*NY])
+/* initialise field with drop - phi is wave height, psi is phi at time t-1 */
+/* phi0 is longidude, theta0 is latitude */
+{
+    int i, j;
+    double xy[2], pscal, dist, dist2, stheta, ctheta;
+
+    ctheta = cos(theta0 + PID);
+    stheta = sin(theta0 + PID);
+    
+    printf("Initializing wave\n"); 
+//     #pragma omp parallel for private(i,j,xy,dist2)
+    for (i=0; i<NX; i++)
+    {
+        if (i%100 == 0) printf("Initializing column %i of %i\n", i, NX);
+        for (j=0; j<NY; j++)
+        {
+            pscal = cos(wsphere[i*NY+j].phi - phi0)*wsphere[i*NY+j].stheta*stheta;
+            pscal += wsphere[i*NY+j].ctheta*ctheta;
+            
+            dist = acos(pscal);
+            dist2 = dist*dist;
+            
+	    xy_in[i*NY+j] = xy_in_billiard_sphere(i, j, wsphere);
+            
+	    if ((xy_in[i*NY+j])||(TWOSPEEDS)) 
+                phi[i*NY+j] = INITIAL_AMP*exp(-dist2/INITIAL_VARIANCE)*cos(-sqrt(dist2)/INITIAL_WAVELENGTH);
+            else phi[i*NY+j] = 0.0;
+            psi[i*NY+j] = 0.0;
+        }
+    }
+}
+
+void add_circular_wave_sphere(double amp, double phi0, double theta0, double phi[NX*NY], double psi[NX*NY], short int xy_in[NX*NY], t_wave_sphere wsphere[NX*NY])
+/* initialise field with drop - phi is wave height, psi is phi at time t-1 */
+/* phi0 is longidude, theta0 is latitude */
+{
+    int i, j;
+    double xy[2], pscal, dist, dist2, stheta, ctheta;
+
+    ctheta = cos(theta0 + PID);
+    stheta = sin(theta0 + PID);
+    
+    printf("Initializing wave\n"); 
+//     #pragma omp parallel for private(i,j,xy,dist2)
+    for (i=0; i<NX; i++)
+    {
+        if (i%100 == 0) printf("Initializing column %i of %i\n", i, NX);
+        for (j=0; j<NY; j++)
+        {
+            pscal = cos(wsphere[i*NY+j].phi - phi0)*wsphere[i*NY+j].stheta*stheta;
+            pscal += wsphere[i*NY+j].ctheta*ctheta;
+            
+            dist = acos(pscal);
+            dist2 = dist*dist;
+            
+	    if ((xy_in[i*NY+j])||(TWOSPEEDS)) 
+                phi[i*NY+j] += amp*INITIAL_AMP*exp(-dist2/INITIAL_VARIANCE)*cos(-sqrt(dist2)/INITIAL_WAVELENGTH);
+            else phi[i*NY+j] = 0.0;
+            psi[i*NY+j] = 0.0;
+        }
+    }
+}
+
+
+void compute_light_angle_sphere(short int xy_in[NX*NY], t_wave wave[NX*NY], t_wave_sphere wsphere[NX*NY], int movie)
+/* computes cosine of angle between normal vector and vector light */
+{
+    int i, j;
+    double x, y, z, norm, pscal, deltai[3], deltaj[3], deltar, n[3], r;
+    static double dphi, dtheta;
+    static int first = 1;
+    
+    if (first)
+    {
+        dphi = DPI/(double)NX;
+        dtheta = PI/(double)NY;
+        first = 0;
+    }
+    
+    #pragma omp parallel for private(i,j)
+    for (i=0; i<NX; i++)
+        for (j=0; j<NY; j++)
+        {
+            r = 1.0 + RSCALE*(*wave[i*NY+j].p_zfield[movie]);
+            if (r > RMAX) r = RMAX;
+            wsphere[i*NY+j].radius = r;
+        }
+    
+    
+    if (SHADE_WAVE)
+    {
+        #pragma omp parallel for private(i,j,norm,pscal,deltar,deltai,deltaj,n)
+        for (i=0; i<NX-1; i++)
+            for (j=0; j<NY-1; j++)
+            {
+                if ((TWOSPEEDS)||(xy_in[i*NY+j]))
+                {
+                    /* computation of tangent vectors */
+                    deltar = (wsphere[(i+1)*NY+j].radius - wsphere[i*NY+j].radius)/dphi;
+                    
+                    deltai[0] = -wsphere[i*NY+j].radius*wsphere[i*NY+j].sphi;
+                    deltai[0] += deltar*wsphere[i*NY+j].cphi;
+                    
+                    deltai[1] = wsphere[i*NY+j].radius*wsphere[i*NY+j].cphi;
+                    deltai[1] += deltar*wsphere[i*NY+j].sphi;
+                    
+                    deltai[2] = -deltar*wsphere[i*NY+j].cottheta;
+                    
+                    deltar = (wsphere[i*NY+j+1].radius - wsphere[i*NY+j].radius)/dtheta;
+                    
+                    deltaj[0] = wsphere[i*NY+j].radius*wsphere[i*NY+j].cphi*wsphere[i*NY+j].ctheta;
+                    deltaj[0] += deltar*wsphere[i*NY+j].cphi*wsphere[i*NY+j].stheta;
+                    
+                    deltaj[1] = wsphere[i*NY+j].radius*wsphere[i*NY+j].sphi*wsphere[i*NY+j].ctheta;
+                    deltaj[1] += deltar*wsphere[i*NY+j].sphi*wsphere[i*NY+j].stheta;
+                    
+                    deltaj[2] = wsphere[i*NY+j].radius*wsphere[i*NY+j].stheta;
+                    deltaj[2] += -deltar*wsphere[i*NY+j].ctheta;
+                    
+                    /* computation of normal vector */
+                    n[0] = deltai[1]*deltaj[2] - deltai[2]*deltaj[1];
+                    n[1] = deltai[2]*deltaj[0] - deltai[0]*deltaj[2];
+                    n[2] = deltai[0]*deltaj[1] - deltai[1]*deltaj[0];
+                                        
+                    norm = sqrt(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]);
+                    
+                    pscal = n[0]*light[0] + n[1]*light[1] + n[2]*light[2];
+                
+                    wave[i*NY+j].cos_angle = pscal/norm;
+                }
+                else 
+                {
+                    pscal = wsphere[i*NY+j].x*light[0] + wsphere[i*NY+j].y*light[1] + wsphere[i*NY+j].z*light[2];
+                
+                    wave[i*NY+j].cos_angle = pscal;
+                }
+            }
+            
+        for (i=0; i<NX-1; i++) wave[i*NY+NY-1].cos_angle = wave[i*NY+NY-2].cos_angle;
+        for (j=0; j<NY-1; j++) wave[(NX-1)*NY+j].cos_angle = wave[(NX-2)*NY+j].cos_angle;
+        wave[(NX-1)*NY+NY-1].cos_angle = wave[(NX-1)*NY+NY-2].cos_angle;
+    }
+    else
+    {
+        #pragma omp parallel for private(i,j,pscal)
+        for (i=0; i<NX; i++)
+            for (j=0; j<NY; j++)
+            {
+                if ((TWOSPEEDS)||(xy_in[i*NY+j]))
+                {
+                    pscal = wsphere[i*NY+j].x*light[0] + wsphere[i*NY+j].y*light[1] + wsphere[i*NY+j].z*light[2];
+                
+                    wave[i*NY+j].cos_angle = pscal;
+                }
+            }
+    }
+}
+
+void compute_light_angle_sphere_2d(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;
+    double gradx, grady, norm, pscal;
+    static double dx, dy, vscale2;
+    static int first = 1;
+    
+    if (first)
+    {
+        dx = 2.0*(XMAX - XMIN)/(double)NX;
+        dy = 2.0*(YMAX - YMIN)/(double)NY;
+        vscale2 = SHADE_SCALE_2D*SHADE_SCALE_2D;
+        first = 0;
+    }
+    
+    #pragma omp parallel for private(i,j,gradx, grady, norm, pscal)
+    for (i=1; i<NX-1; i++)
+        for (j=1; j<NY-1; j++)
+        {
+            if ((TWOSPEEDS)||(xy_in[i*NY+j]))
+            {
+                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(vscale2 + gradx*gradx + grady*grady);
+                pscal = -gradx*light[0] - grady*light[1] + SHADE_SCALE_2D;
+                
+                wave[i*NY+j].cos_angle = pscal/norm;
+            }
+        }
+    
+    /* i=0 */
+    for (j=1; j<NY-1; j++)
+    {
+        if ((TWOSPEEDS)||(xy_in[j]))
+        {
+            gradx = (*wave[NY+j].p_zfield[movie] - *wave[(NY-1)*NY+j].p_zfield[movie])/dx;
+            grady = (*wave[j+1].p_zfield[movie] - *wave[j-1].p_zfield[movie])/dy;
+                
+            norm = sqrt(vscale2 + gradx*gradx + grady*grady);
+            pscal = -gradx*light[0] - grady*light[1] + SHADE_SCALE_2D;
+                
+            wave[j].cos_angle = pscal/norm;
+        }
+    }
+    
+    /* i=NX-1 */
+    for (j=1; j<NY-1; j++)
+    {
+        if ((TWOSPEEDS)||(xy_in[(NX-1)*NY+j]))
+        {
+            gradx = (*wave[j].p_zfield[movie] - *wave[(NX-2)*NY+j].p_zfield[movie])/dx;
+            grady = (*wave[(NX-1)*NY+j+1].p_zfield[movie] - *wave[(NX-1)*NY+j-1].p_zfield[movie])/dy;
+                
+            norm = sqrt(vscale2 + gradx*gradx + grady*grady);
+            pscal = -gradx*light[0] - grady*light[1] + SHADE_SCALE_2D;
+                
+            wave[(NX-1)*NY+j].cos_angle = pscal/norm;
+        }
+    }
+}
+
+
+void compute_cfield_sphere(short int xy_in[NX*NY], int cplot, int palette, t_wave wave[NX*NY], int fade, double fade_value, int movie)
+/* compute the colors */
+{
+    int i, j, k;
+    double ca;
+    static double fact_max;
+    static int first = 1;
+    
+    if (first)
+    {
+        fact_max = 0.5*COS_LIGHT_MAX;
+        first = 0;
+    }
+       
+    #pragma omp parallel for private(i,j,ca)
+    for (i=0; i<NX; i++) for (j=0; j<NY; j++)
+    {
+        compute_field_color(*wave[i*NY+j].p_cfield[movie], *wave[i*NY+j].p_cfield[movie+2], cplot, palette, wave[i*NY+j].rgb);
+        if ((SHADE_3D)||(SHADE_2D))
+        {
+            ca = wave[i*NY+j].cos_angle;
+//             ca = (ca + 1.0)*0.4 + 0.2;
+            ca = (ca + 1.0)*fact_max + COS_LIGHT_MIN;
+            if ((FADE_IN_OBSTACLE)&&(!xy_in[i*NY+j])) ca = (ca + 0.1)*1.6;
+            for (k=0; k<3; k++) wave[i*NY+j].rgb[k] *= ca;
+        }
+        if (fade) for (k=0; k<3; k++) wave[i*NY+j].rgb[k] *= fade_value;
+    }
+}
+
+
+void draw_wave_sphere_2D(int movie, double phi[NX*NY], double psi[NX*NY], short int xy_in[NX*NY], t_wave wave[NX*NY], t_wave_sphere wsphere[NX*NY], int zplot, int cplot, int palette, int fade, double fade_value, int refresh)
+{
+    int i, j, ii;
+    double x, y, ca;
+    static int ishift;
+    static double dx, dy;
+    
+    if (refresh)
+    {
+        compute_wave_fields(phi, psi, xy_in, zplot, cplot, wave);
+        if (SHADE_3D) compute_light_angle_sphere(xy_in, wave, wsphere, movie);
+        else if (SHADE_2D) compute_light_angle_sphere_2d(xy_in, wave, movie);
+        compute_cfield_sphere(xy_in, cplot, palette, wave, fade, fade_value, movie);
+        dx = (XMAX - XMIN)/(double)NX;
+        dy = (YMAX - YMIN)/(double)(NY-2*DPOLE);
+        ishift = (int)((double)NX*PHISHIFT/360.0);
+    }
+    
+    glBegin(GL_QUADS);
+    
+    for (i=0; i<NX; i++)
+        for (j=0; j<NY; j++)
+        {
+            if ((TWOSPEEDS)||(xy_in[i*NY+j]))
+                glColor3f(wave[i*NY+j].rgb[0], wave[i*NY+j].rgb[1], wave[i*NY+j].rgb[2]);
+            else
+            {
+                ca = wave[i*NY+j].cos_angle;
+                ca = (ca + 1.0)*0.4 + 0.2;
+                if (fade) ca *= fade_value;
+                if (B_DOMAIN == D_SPHERE_EARTH)
+                    glColor3f(wsphere[i*NY+j].r*ca, wsphere[i*NY+j].g*ca, wsphere[i*NY+j].b*ca);
+                else glColor3f(COLOR_OUT_R*ca, COLOR_OUT_G*ca, COLOR_OUT_B*ca);
+            }
+
+            x = wsphere[i*NY+j].x2d;
+            y = wsphere[i*NY+j].y2d;
+            
+            ii = NX-i-1+ishift;
+            if (ii > NX) ii -= NX;
+            
+            glVertex2i(ii, j);
+            glVertex2i(ii+1, j);
+            glVertex2i(ii+1, j+1);
+            glVertex2i(ii, j+1);
+        }
+    glEnd ();
+}
+
+void draw_wave_sphere_ij(int i, int iplus, int j, int jplus, int jcolor, int movie, double phi[NX*NY], double psi[NX*NY], short int xy_in[NX*NY], t_wave wave[NX*NY], t_wave_sphere wsphere[NX*NY], int zplot, int cplot, int palette, int fade, double fade_value)
+/* draw wave at simulation grid point (i,j) */
+{
+    int k, l;
+    double xyz[3], ca;
+    
+    if ((TWOSPEEDS)||(xy_in[i*NY+j]))
+        glColor3f(wave[i*NY+jcolor].rgb[0], wave[i*NY+jcolor].rgb[1], wave[i*NY+jcolor].rgb[2]);
+    else
+    {
+        ca = wave[i*NY+j].cos_angle;
+        ca = (ca + 1.0)*0.4 + 0.2;
+        if (fade) ca *= fade_value;
+        if (B_DOMAIN == D_SPHERE_EARTH)
+            glColor3f(wsphere[i*NY+j].r*ca, wsphere[i*NY+j].g*ca, wsphere[i*NY+j].b*ca);
+        else glColor3f(COLOR_OUT_R*ca, COLOR_OUT_G*ca, COLOR_OUT_B*ca);
+    }
+    glBegin(GL_TRIANGLE_FAN);
+    if (ij_to_sphere(i, j, *wave[i*NY+j].p_zfield[movie], wsphere, xyz))
+        draw_vertex_sphere(xyz);
+    if (ij_to_sphere(iplus, j, *wave[iplus*NY+j].p_zfield[movie], wsphere, xyz))
+        draw_vertex_sphere(xyz);
+    if (ij_to_sphere(iplus, jplus, *wave[iplus*NY+j+1].p_zfield[movie], wsphere, xyz))
+        draw_vertex_sphere(xyz);
+    if (ij_to_sphere(i, jplus, *wave[i*NY+j+1].p_zfield[movie], wsphere, xyz))
+        draw_vertex_sphere(xyz);
+    glEnd ();
+}
+
+
+void draw_wave_sphere_3D(int movie, double phi[NX*NY], double psi[NX*NY], short int xy_in[NX*NY], t_wave wave[NX*NY], t_wave_sphere wsphere[NX*NY], int zplot, int cplot, int palette, int fade, double fade_value, int refresh)
+{
+    int i, j, imax, imin;
+    double observer_angle, angle2;
+    
+    blank();
+            
+    if (refresh)
+    {
+        compute_wave_fields(phi, psi, xy_in, zplot, cplot, wave);
+        if (SHADE_3D) compute_light_angle_sphere(xy_in, wave, wsphere, movie);
+        else if (SHADE_2D) compute_light_angle_sphere_2d(xy_in, wave, movie);
+        compute_cfield_sphere(xy_in, cplot, palette, wave, fade, fade_value, movie);
+    }
+    
+    observer_angle = argument(observer[0], observer[1]);
+    if (observer_angle < 0.0) observer_angle += DPI;
+    
+    angle2 = observer_angle + PI;
+    if (angle2 > DPI) angle2 -= DPI;
+    
+    imin = (int)(observer_angle*(double)NX/DPI);
+    imax = (int)(angle2*(double)NX/DPI);
+    if (imin >= NX-1) imin = NX-2;
+    if (imax >= NX-1) imax = NX-2;
+    
+//     printf("Angle = %.5lg, angle2 = %.5lg, imin = %i, imax = %i\n", observer_angle, angle2, imin, imax);
+    
+    if (observer[2] > 0.0)
+    {
+        if (imin < imax)
+        {
+            for (i=imax; i>imin; i--)
+                for (j=0; j<NY-2; j++)
+                    draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
+        
+            for (i=imax+1; i<NX-1; i++)
+                for (j=0; j<NY-2; j++)
+                    draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
+        
+            for (j=0; j<NY-2; j++)
+                draw_wave_sphere_ij(NX-1, 0, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
+        
+            for (i=0; i<=imin; i++)
+                for (j=0; j<NY-2; j++)
+                    draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
+        }
+        else
+        {
+            for (i=imax; i<imin; i++)
+                for (j=0; j<NY-2; j++)
+                    draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
+        
+            for (i=imax-1; i>=0; i--)
+                for (j=0; j<NY-2; j++)
+                    draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
+        
+            for (j=0; j<NY-2; j++)
+                draw_wave_sphere_ij(NX-1, 0, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
+        
+            for (i=NX-2; i>=imin; i--)
+                for (j=0; j<NY-2; j++)
+                    draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
+        }
+    
+        /* North pole */
+        for (i=0; i<NX-1; i++) for (j=NY-2; j<NY-1; j++)
+            draw_wave_sphere_ij(i, i+1, j-1, j, j, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
+        
+        for (j=NY-2; j<NY-1; j++)
+            draw_wave_sphere_ij(NX-1, 0, j-1, j, NY-DPOLE, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
+    }
+    else
+    {
+        if (imin < imax)
+        {
+            for (i=imax; i>imin; i--)
+                for (j=NY-3; j>=0; j--)
+                    draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
+        
+            for (i=imax+1; i<NX-1; i++)
+                for (j=NY-3; j>=0; j--)
+                    draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
+        
+            for (j=NY-3; j>=0; j--)
+                draw_wave_sphere_ij(NX-1, 0, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
+        
+            for (i=0; i<=imin; i++)
+                for (j=NY-3; j>=0; j--)
+                    draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
+        }
+        else
+        {
+            for (i=imax; i<imin; i++)
+                for (j=NY-3; j>=0; j--)
+                    draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
+        
+            for (i=imax-1; i>=0; i--)
+                for (j=NY-3; j>=0; j--)
+                    draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
+        
+            for (j=NY-3; j>=0; j--)
+                draw_wave_sphere_ij(NX-1, 0, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
+        
+            for (i=NX-2; i>=imin; i--)
+                for (j=NY-3; j>=0; j--)
+                    draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
+        }
+    
+        /* South pole */
+        for (i=0; i<NX-1; i++) for (j=2; j>0; j--)
+            draw_wave_sphere_ij(i, i+1, j-1, j, j, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
+        
+        for (j=2; j>0; j--)
+            draw_wave_sphere_ij(NX-1, 0, j-1, j, DPOLE, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
+    }
+}
+
+void draw_wave_sphere(int movie, double phi[NX*NY], double psi[NX*NY], short int xy_in[NX*NY], t_wave wave[NX*NY], t_wave_sphere wsphere[NX*NY], int zplot, int cplot, int palette, int fade, double fade_value, int refresh)
+{
+    if (PLOT_2D) draw_wave_sphere_2D(movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade,fade_value, refresh);
+    else draw_wave_sphere_3D(movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade,fade_value, refresh);
+}

sub_wave.c

@@ -293,6 +293,15 @@ void write_text( double x, double y, char *st)
 	return(alph);
  }
  
+ double iabs(int i)     /* absolute value */
+ {
+	int res;
+
+	if (i<0) res = -i;
+	else res = i;
+	return(i);
+ }
+
  
 double gaussian()
 /* returns standard normal random variable, using Box-Mueller algorithm */
@@ -446,6 +455,16 @@ void erase_area(double x, double y, double dx, double dy)
     glEnd();
 }
 
+void erase_area_ij(int imin, int jmin, int imax, int jmax)
+{
+    glColor3f(0.0, 0.0, 0.0);
+    glBegin(GL_QUADS);
+    glVertex2i(imin, jmin);
+    glVertex2i(imax, jmin);
+    glVertex2i(imax, jmax);
+    glVertex2i(imin, jmax);
+    glEnd();
+}
 
 void erase_area_rgb(double x, double y, double dx, double dy, double rgb[3])
 {
@@ -1004,6 +1023,63 @@ int init_circle_config_pattern(t_circle circles[NMAXCIRCLES], int circle_pattern
             }
             break;
         }
+        case (C_HEX_BOTTOM):
+        {
+            ncircles = NGRIDX*(NGRIDY+1);
+            dy = (- YMIN)/((double)NGRIDY);
+//             dx = dy*0.5*sqrt(3.0);
+            dx = (XMAX - XMIN)/((double)NGRIDX);
+            for (i = 0; i < NGRIDX; i++)
+                for (j = 0; j < NGRIDY+1; j++)
+                {
+                    n = (NGRIDY+1)*i + j;
+                    circles[n].xc = ((double)(i-NGRIDX/2) + 0.5)*dx;   /* is +0.5 needed? */
+                    circles[n].yc = YMIN + ((double)j - 0.5)*dy;
+                    if ((i+NGRIDX)%2 == 1) circles[n].yc += 0.5*dy;
+                    circles[n].radius = MU;
+                    /* activate only circles that intersect the domain */
+                    if ((circles[n].yc < YMAX + MU)&&(circles[n].yc > YMIN - MU)) circles[n].active = 1;
+                    else circles[n].active = 0;
+                }
+            break;
+        }
+        case (C_HEX_BOTTOM2):
+        {
+            ncircles = NGRIDX*(NGRIDY+1);
+            dy = (- YMIN)/((double)NGRIDY);
+            dx = 2.0*LAMBDA/((double)NGRIDX);
+            for (i = 0; i < NGRIDX; i++)
+                for (j = 0; j < NGRIDY+1; j++)
+                {
+                    n = (NGRIDY+1)*i + j;
+                    circles[n].xc = ((double)(i-NGRIDX/2) + 0.5)*dx;   /* is +0.5 needed? */
+                    circles[n].yc = YMIN + ((double)j - 0.5)*dy;
+                    if ((i+NGRIDX)%2 == 1) circles[n].yc += 0.5*dy;
+                    circles[n].radius = MU;
+                    /* activate only circles that intersect the domain */
+                    if ((circles[n].yc < YMAX + MU)&&(circles[n].yc > YMIN - MU)) circles[n].active = 1;
+                    else circles[n].active = 0;
+                }
+            break;
+        }
+        case (C_SQUARE_BOTTOM):
+        {
+            ncircles = NGRIDX*(NGRIDY+1);
+            dy = (- YMIN)/((double)NGRIDY);
+            dx = 2.0*LAMBDA/((double)NGRIDX);
+            for (i = 0; i < NGRIDX; i++)
+                for (j = 0; j < NGRIDY+1; j++)
+                {
+                    n = (NGRIDY+1)*i + j;
+                    circles[n].xc = ((double)(i-NGRIDX/2) + 0.5)*dx;  
+                    circles[n].yc = YMIN + ((double)j - 0.5)*dy;
+                    circles[n].radius = MU;
+                    /* activate only circles that intersect the domain */
+                    if ((circles[n].yc < YMAX + MU)&&(circles[n].yc > YMIN - MU)) circles[n].active = 1;
+                    else circles[n].active = 0;
+                }
+            break;
+        }
         case (C_ONE):
         {
             circles[ncircles].xc = 0.0;
@@ -5752,7 +5828,7 @@ void compute_laplacian(double phi[NX*NY], t_laplacian laplace[NX*NY], double del
     }
 }
 
-double oscillating_bc(int time)
+double oscillating_bc(int time, int j)
 {
     double t, phase, a, envelope, omega;
     
@@ -5760,7 +5836,12 @@ double oscillating_bc(int time)
     {
         case (OSC_PERIODIC): 
         {
-            return(AMPLITUDE*cos((double)time*OMEGA)*exp(-(double)time*DAMPING));
+            if (OSCIL_LEFT_YSHIFT > 0.0)
+                t = (double)time*OMEGA - (double)j*OSCIL_LEFT_YSHIFT/(double)NY;
+            else 
+                t = (double)time*OMEGA + (double)(NY-j)*OSCIL_LEFT_YSHIFT/(double)NY;
+            if (t < 0.0) return(0.0);
+            else return(AMPLITUDE*cos(t)*exp(-(double)t*DAMPING));
         }
         case (OSC_SLOWING):
         {

sub_wave_3d.c

@@ -1850,6 +1850,7 @@ void draw_wave_3d_ij(int i, int j, int movie, double phi[NX*NY], double psi[NX*N
         {
             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);
+            /* TODO: incorporate these in wave structure */
             ij_to_xy(i, j, xy_sw);
             ij_to_xy(i+1, j, xy_se);
             ij_to_xy(i, j+1, xy_nw);
@@ -1863,6 +1864,7 @@ void draw_wave_3d_ij(int i, int j, int movie, double phi[NX*NY], double psi[NX*N
 //                 z_mid += *wave[(i+1)*NY+j+1].potential*POT_FACT*0.25;
 //             }
                     
+            /* TODO: incorporate these in wave structure */
             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)
@@ -2198,6 +2200,162 @@ void draw_color_scheme_palette_3d(double x1, double y1, double x2, double y2, in
     draw_rectangle_noscale(x1, y1, x2, y2);
 }
 
+void draw_circular_color_scheme_palette_3d(double x1, double y1, double radiusx, double radiusy, int plot, double min, double max, int palette, int fade, double fade_value)
+{
+    int j, k, ij[2], jmin=0;
+    double x, y, dy, dy_e, dy_phase, rgb[3], value, lum, amp, dphi, pos[2], phi, xy[2];
+    
+    printf("Drawing circular color scheme\n");
+    
+    glBegin(GL_TRIANGLE_FAN);
+//     xy_to_pos(x1, y1, xy);
+//     glVertex2d(xy[0], xy[1]);
+    xy_to_ij(x1, y1, ij);
+    draw_vertex_ij(ij[0], ij[1]);
+    dy = (max - min)/360.0;
+    dy_e = max/360.0;
+    dy_phase = 1.0/360.0;
+    dphi = DPI/360.0;
+    
+    for (j = 0; j <= 360; j++)
+    {
+        switch (plot) {
+            case (P_3D_AMPLITUDE):
+            {
+                value = min + 1.0*dy*(double)(j - jmin);
+                color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb);
+                break;
+            }
+            case (P_3D_ANGLE):
+            {
+                value = 1.0*dy*(double)(j - jmin);
+                color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 0, rgb);
+                break;
+            }
+            case (P_3D_AMP_ANGLE):
+            {
+                value = min + 1.0*dy*(double)(j - jmin);
+                color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb);
+                break;
+            }
+            case (P_3D_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_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_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 = dy_phase*(double)(j - jmin);
+                color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb);
+                break;
+            }
+            case (P_3D_FLUX_INTENSITY):
+            {
+                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_FLUX_DIRECTION):
+            {
+                value = dy_phase*(double)(j - jmin);
+                color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb);
+                break;
+            }
+            case (Z_EULER_VORTICITY):
+            {
+                value = min + 1.0*dy*(double)(j - jmin);
+                color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb);
+                break;
+            }
+            case (Z_EULER_LOG_VORTICITY):
+            {
+                value = min + 1.0*dy*(double)(j - jmin);
+                color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb);
+                break;
+            }
+            case (Z_EULER_VORTICITY_ASYM):
+            {
+                value = min + 1.0*dy*(double)(j - jmin);
+                color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb);
+                break;
+            }
+            case (Z_EULER_LPRESSURE):
+            {
+                value = min + 1.0*dy*(double)(j - jmin);
+                color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb);
+                break;
+            }
+            case (Z_EULERC_VORTICITY):
+             {
+                value = min + 1.0*dy*(double)(j - jmin);
+                printf("Palette value %.3lg\n", value);
+                color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb);
+                break;
+            }
+               
+        }
+        if (fade) for (k=0; k<3; k++) rgb[k] *= fade_value;
+        glColor3f(rgb[0], rgb[1], rgb[2]);
+        
+        xy_to_ij(x1 + radiusx*cos(dphi*(double)j), y1 + radiusy*sin(dphi*(double)j), ij);
+        draw_vertex_ij(ij[0], ij[1]);
+    }
+    
+    xy_to_ij(x1 + radiusx*cos(dphi), y1 + radiusy*sin(dphi), ij);
+    draw_vertex_ij(ij[0], ij[1]);
+    glEnd ();
+    
+    if (fade) glColor3f(fade_value, fade_value, fade_value);
+    else glColor3f(1.0, 1.0, 1.0);
+    glLineWidth(BOUNDARY_WIDTH*3/2);
+    glEnable(GL_LINE_SMOOTH);
+    
+    dphi = DPI/(double)NSEG;
+    glBegin(GL_LINE_LOOP);
+    for (j = 0; j < NSEG; j++)
+    {               
+        xy_to_ij(x1 + radiusx*cos(dphi*(double)j), y1 + radiusy*sin(dphi*(double)j), ij);
+        draw_vertex_ij(ij[0], ij[1]);
+    }
+    glEnd ();
+}
+
+
 void print_speed_3d(double speed, int fade, double fade_value)
 {
     char message[100];

wave_3d.c

@@ -156,6 +156,7 @@
 #define KAPPA 0.0           /* "elasticity" term enforcing oscillations */
 #define KAPPA_SIDES 5.0e-4  /* "elasticity" term on absorbing boundary */
 #define KAPPA_TOPBOT 0.0    /* "elasticity" term on absorbing boundary */
+#define OSCIL_LEFT_YSHIFT 0.0   /* y-dependence of left oscillation (for non-horizontal waves) */
 /* The Courant number is given by c*DT/DX, where DT is the time step and DX the lattice spacing */
 /* The physical damping coefficient is given by GAMMA/(DT)^2 */
 /* Increasing COURANT speeds up the simulation, but decreases accuracy */
@@ -289,6 +290,7 @@
 // #define POT_MAZE 7
 #define POTENTIAL 10
 #define POT_FACT 20.0
+#define DRAW_WAVE_PROFILE 0     /* set to 1 to draw a profile of the wave */
 /* end of constants only used by sub_wave and sub_maze */
 
 
@@ -364,8 +366,8 @@ void evolve_wave_half_new(double phi_in[NX*NY], double psi_in[NX*NY], double phi
 //     if (OSCILLATE_LEFT) for (j=1; j<NY-1; j++) phi_out[j] = AMPLITUDE*cos((double)time*OMEGA);
     if (OSCILLATE_LEFT) 
     {
-        for (j=1; j<NY-1; j++) phi_out[j] = oscillating_bc(time);
-        printf("Boundary condition %.3lg\n", oscillating_bc(time));
+        for (j=1; j<NY-1; j++) phi_out[j] = oscillating_bc(time, j);
+        printf("Boundary condition %.3lg\n", oscillating_bc(time, 0));
     }
     else for (j=1; j<NY-1; j++){
         if ((TWOSPEEDS)||(xy_in[j] != 0)){
@@ -574,7 +576,7 @@ void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out
 //     if (OSCILLATE_LEFT) for (j=1; j<NY-1; j++) phi_out[j] = AMPLITUDE*cos((double)time*OMEGA);
     if (OSCILLATE_LEFT) 
     {
-        for (j=1; j<NY-1; j++) phi_out[j] = oscillating_bc(time);
+        for (j=1; j<NY-1; j++) phi_out[j] = oscillating_bc(time, j);
 //         printf("Boundary condition %.3lg\n", oscillating_bc(time));
     }
     else for (j=1; j<NY-1; j++){
@@ -865,8 +867,8 @@ void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out
     /* add oscillating boundary condition on the left corners */
     if (OSCILLATE_LEFT)
     {
-        phi_out[0] = oscillating_bc(time);
-        phi_out[NY-1] = oscillating_bc(time);
+        phi_out[0] = oscillating_bc(time, 0);
+        phi_out[NY-1] = oscillating_bc(time, NY-1);
     }
     
     

wave_billiard.c

@@ -24,7 +24,7 @@
 /*  create movie using                                                           */
 /*  ffmpeg -i wave.%05d.tif -vcodec libx264 wave.mp4                             */
 /*                                                                               */
-/******************************F***************************************************/
+/**********************************************************************************/
 
 /*********************************************************************************/
 /*                                                                               */
@@ -61,10 +61,12 @@
 #define NX 3840          /* number of grid points on x axis */
 #define NY 2300          /* number of grid points on y axis */
 
-#define XMIN -1.0
-#define XMAX 3.0	/* x interval */
-#define YMIN -1.197916667
-#define YMAX 1.197916667	/* y interval for 9/16 aspect ratio */
+#define XMIN -2.0
+#define XMAX 2.0	/* x interval */
+#define YMIN -0.397916667
+#define YMAX 1.997916667	/* y interval for 9/16 aspect ratio */
+// #define YMIN -1.197916667
+// #define YMAX 1.197916667	/* y interval for 9/16 aspect ratio */
 
 #define HIGHRES 1       /* set to 1 if resolution of grid is double that of displayed image */
 
@@ -87,9 +89,9 @@
 
 /* Choice of the billiard table */
 
-#define B_DOMAIN 10        /* choice of domain shape, see list in global_pdes.c */
+#define B_DOMAIN 20        /* choice of domain shape, see list in global_pdes.c */
 
-#define CIRCLE_PATTERN 1   /* pattern of circles or polygons, see list in global_pdes.c */
+#define CIRCLE_PATTERN 103   /* pattern of circles or polygons, see list in global_pdes.c */
 
 #define COMPARISON 0        /* set to 1 to compare two different patterns (beta) */
 #define B_DOMAIN_B 20       /* second domain shape, for comparisons */
@@ -99,8 +101,8 @@
 #define NPOISSON 1000        /* number of points for Poisson C_RAND_POISSON arrangement */
 #define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */
 
-#define LAMBDA 0.1197916667	    /* parameter controlling the dimensions of domain */
-#define MU 0.035             /* parameter controlling the dimensions of domain */
+#define LAMBDA 1.0	    /* parameter controlling the dimensions of domain */
+#define MU 0.005            /* parameter controlling the dimensions of domain */
 #define NPOLY 6             /* number of sides of polygon */
 #define APOLY 0.0           /* angle by which to turn polygon, in units of Pi/2 */ 
 #define MDEPTH 6            /* depth of computation of Menger gasket */
@@ -108,8 +110,8 @@
 #define MANDELLEVEL 1000    /* iteration level for Mandelbrot set */
 #define MANDELLIMIT 10.0    /* limit value for approximation of Mandelbrot set */
 #define FOCI 1              /* set to 1 to draw focal points of ellipse */
-#define NGRIDX 30            /* number of grid point for grid of disks */
-#define NGRIDY 30            /* number of grid point for grid of disks */
+#define NGRIDX 60           /* number of grid point for grid of disks */
+#define NGRIDY 10           /* number of grid point for grid of disks */
 // #define NGRIDY 18            /* number of grid point for grid of disks */
 
 #define X_SHOOTER -0.2
@@ -129,11 +131,11 @@
 /* Physical parameters of wave equation */
 
 #define TWOSPEEDS 0          /* set to 1 to replace hardcore boundary by medium with different speed */
-#define OSCILLATE_LEFT 0     /* set to 1 to add oscilating boundary condition on the left */
+#define OSCILLATE_LEFT 1     /* set to 1 to add oscilating boundary condition on the left */
 #define OSCILLATE_TOPBOT 0   /* set to 1 to enforce a planar wave on top and bottom boundary */
 #define OSCILLATION_SCHEDULE 0  /* oscillation schedule, see list in global_pdes.c */
 
-#define OMEGA 0.004        /* frequency of periodic excitation */
+#define OMEGA 0.024        /* frequency of periodic excitation */
 #define AMPLITUDE 1.0      /* amplitude of periodic excitation */ 
 #define ACHIRP 0.25        /* acceleration coefficient in chirp */
 #define DAMPING 0.0        /* damping of periodic excitation */
@@ -146,12 +148,13 @@
 #define KAPPA 0.0           /* "elasticity" term enforcing oscillations */
 #define KAPPA_SIDES 5.0e-4  /* "elasticity" term on absorbing boundary */
 #define KAPPA_TOPBOT 0.0    /* "elasticity" term on absorbing boundary */
+#define OSCIL_LEFT_YSHIFT -400.0   /* y-dependence of left oscillation (for non-horizontal waves) */
 /* The Courant number is given by c*DT/DX, where DT is the time step and DX the lattice spacing */
 /* The physical damping coefficient is given by GAMMA/(DT)^2 */
 /* Increasing COURANT speeds up the simulation, but decreases accuracy */
 /* For similar wave forms, COURANT^2*GAMMA should be kept constant */
 
-#define ADD_OSCILLATING_SOURCE 1        /* set to 1 to add an oscillating wave source */
+#define ADD_OSCILLATING_SOURCE 0        /* set to 1 to add an oscillating wave source */
 #define OSCILLATING_SOURCE_PERIOD 30    /* period of oscillating source */
 #define ALTERNATE_OSCILLATING_SOURCE 1  /* set to 1 to alternate sign of oscillating source */
 
@@ -166,11 +169,11 @@
 
 /* Parameters for length and speed of simulation */
 
-#define NSTEPS 2400       /* number of frames of movie */
-// #define NSTEPS 1300       /* number of frames of movie */
-#define NVID 6            /* number of iterations between images displayed on screen */
+#define NSTEPS 3600       /* number of frames of movie */
+// #define NSTEPS 500       /* number of frames of movie */
+#define NVID 7            /* number of iterations between images displayed on screen */
 #define NSEG 1000         /* number of segments of boundary */
-#define INITIAL_TIME 0      /* time after which to start saving frames */
+#define INITIAL_TIME 700      /* time after which to start saving frames */
 #define BOUNDARY_WIDTH 2    /* width of billiard boundary */
 #define PRINT_SPEED 0       /* print speed of moving source */
 #define PRINT_FREQUENCY 0       /* print frequency (for phased array) */
@@ -199,9 +202,9 @@
 
 /* Color schemes */
 
-// #define COLOR_PALETTE 15     /* Color palette, see list in global_pdes.c  */
-#define COLOR_PALETTE 17      /* Color palette, see list in global_pdes.c  */
-#define COLOR_PALETTE_B 13     /* Color palette, see list in global_pdes.c  */
+#define COLOR_PALETTE 18     /* Color palette, see list in global_pdes.c  */
+// #define COLOR_PALETTE 17      /* Color palette, see list in global_pdes.c  */
+#define COLOR_PALETTE_B 12     /* Color palette, see list in global_pdes.c  */
 
 #define BLACK 1          /* background */
 
@@ -212,9 +215,9 @@
 #define PHASE_FACTOR 1.0       /* factor in computation of phase in color scheme P_3D_PHASE */
 #define PHASE_SHIFT 0.0      /* shift of phase in color scheme P_3D_PHASE */
 #define ATTENUATION 0.0  /* exponential attenuation coefficient of contrast with time */
-#define E_SCALE 100.0     /* scaling factor for energy representation */
+#define E_SCALE 50.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 LOG_SHIFT 1.5     /* shift of colors on log scale */
 #define FLUX_SCALE 5.0e3    /* scaling factor for enegy flux represtnation */
 #define RESCALE_COLOR_IN_CENTER 0   /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */
 
@@ -226,12 +229,14 @@
 #define HUEAMP -180.0      /* amplitude of variation of hue for color scheme C_HUE */
 
 #define DRAW_COLOR_SCHEME 1    /* set to 1 to plot the color scheme */
-#define COLORBAR_RANGE 2.0     /* scale of color scheme bar */
-#define COLORBAR_RANGE_B 2.0   /* scale of color scheme bar for 2nd part */
+#define COLORBAR_RANGE 1.0     /* scale of color scheme bar */
+#define COLORBAR_RANGE_B 0.4   /* scale of color scheme bar for 2nd part */
 #define ROTATE_COLOR_SCHEME 0   /* set to 1 to draw color scheme horizontally */
 #define CIRC_COLORBAR 0         /* set to 1 to draw circular color scheme */
 #define CIRC_COLORBAR_B 0       /* set to 1 to draw circular color scheme */
 
+#define DRAW_WAVE_PROFILE 1     /* set to 1 to draw a profile of the wave */
+
 #define SAVE_TIME_SERIES 0      /* set to 1 to save wave time series at a point */
 
 #define NXMAZE 8      /* width of maze */
@@ -328,7 +333,7 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
     }
     
     /* left boundary */
-    if (OSCILLATE_LEFT) for (j=1; j<NY-1; j++) phi_out[0][j] = oscillating_bc(time);
+    if (OSCILLATE_LEFT) for (j=1; j<NY-1; j++) phi_out[0][j] = oscillating_bc(time, j);
     else for (j=1; j<NY-1; j++){
         if ((TWOSPEEDS)||(xy_in[0][j] != 0)){
             x = phi_in[0][j];
@@ -519,8 +524,8 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
     /* add oscillating boundary condition on the left corners */
     if (OSCILLATE_LEFT)
     {
-        phi_out[0][0] = oscillating_bc(time);
-        phi_out[0][NY-1] = oscillating_bc(time);
+        phi_out[0][0] = oscillating_bc(time, 0);
+        phi_out[0][NY-1] = oscillating_bc(time, NY-1);
     }
     
     /* for debugging purposes/if there is a risk of blow-up */
@@ -676,9 +681,9 @@ void animation()
 //     xy_to_ij(startright[0], startright[1], sample_right);
 //     printf("xleft = (%.3f, %.3f) xright = (%.3f, %.3f)\n", xin_left, yin_left, xin_right, yin_right);
     
-//     init_wave_flat(phi, psi, xy_in);
+    init_wave_flat(phi, psi, xy_in);
 
-    init_circular_wave(-0.5, 0.0, phi, psi, xy_in);
+//     init_circular_wave(-0.5, 0.0, phi, psi, xy_in);
 //     x = XMIN + (XMAX - XMIN)*rand()/RAND_MAX;
 //     y = YMIN + (YMAX - YMIN)*rand()/RAND_MAX;
 //     init_circular_wave(0.0, -0.8, phi, psi, xy_in);

wave_common.c

@@ -853,12 +853,171 @@ void draw_wave_epalette(double *phi[NX], double *psi[NX], double *total_energy[N
 }
 
 
+double wave_value(int i, int j, double *phi[NX], double *psi[NX], double *total_energy[NX], double *total_flux, short int *xy_in[NX], double scale, int time, int plot, int palette, double rgb[3])
+/* compute value of wave and color */
+{
+    int k;
+    double value, velocity, energy, gradientx2, gradienty2, arg, mod, flux_factor, gx, gy, mgx, mgy;
+    
+    switch (plot) {
+        case (P_AMPLITUDE):
+        {
+            value = phi[i][j];
+            color_scheme_palette(COLOR_SCHEME, palette, value, scale, time, rgb);
+            break;
+        }
+        case (P_ENERGY):
+        {
+            value = compute_energy(phi, psi, xy_in, i, j);
+            /* adjust energy to color palette */
+            if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, scale, time, rgb);
+            else color_scheme_palette(COLOR_SCHEME, palette, value, scale, time, rgb);
+            break;
+        }
+        case (P_MIXED):
+        {
+            if (j > NY/2) value = phi[i][j];
+            else value = compute_energy(phi, psi, xy_in, i, j);
+            color_scheme_palette(COLOR_SCHEME, palette, value, scale, time, rgb);
+            break;
+        }
+        case (P_MEAN_ENERGY):
+        {
+            energy = compute_energy(phi, psi, xy_in, i, j);
+            total_energy[i][j] += energy;
+            value = total_energy[i][j]/(double)(time+1);
+            if (COLOR_PALETTE >= COL_TURBO) 
+                color_scheme_asym_palette(COLOR_SCHEME, palette, value, scale, time, rgb);
+            else color_scheme_palette(COLOR_SCHEME, palette, value, scale, time, rgb);
+            break;
+        }
+        case (P_LOG_ENERGY):
+        {
+            energy = compute_energy(phi, psi, xy_in, i, j);
+//                       energy = LOG_SHIFT + LOG_SCALE*log(energy);
+//                         if (energy < 0.0) energy = 0.0;
+            value = LOG_SHIFT + LOG_SCALE*log(energy);
+            color_scheme_palette(COLOR_SCHEME, palette, value, scale, time, rgb);
+            break;
+        }
+        case (P_LOG_MEAN_ENERGY):
+        {
+            energy = compute_energy(phi, psi, xy_in, i, j);
+            if (energy == 0.0) energy = 1.0e-20;
+            total_energy[i][j] += energy;
+            value = LOG_SHIFT + LOG_SCALE*log(total_energy[i][j]/(double)(time+1));
+            color_scheme_palette(COLOR_SCHEME, palette, value, scale, time, rgb);
+            break;
+        }
+        case (P_ENERGY_FLUX):
+        {
+            compute_energy_flux(phi, psi, xy_in, i, j, &gx, &gy, &arg, &mod);
+//                         color_scheme_palette(C_ONEDIM_LINEAR, palette, arg/DPI, 1.0, 1, rgb);
+//                         flux_factor = tanh(mod*E_SCALE);
+//                         for (k=0; k<3; k++) rgb[k] *= flux_factor;
+            value = mod*FLUX_SCALE;
+            color_scheme_asym_palette(COLOR_SCHEME, palette, value, scale, time, rgb);
+            break;
+        }
+        case (P_TOTAL_ENERGY_FLUX):
+        {
+            compute_energy_flux(phi, psi, xy_in, i, j, &gx, &gy, &arg, &mod);
+            total_flux[2*j*NX + 2*i] *= 0.99;
+            total_flux[2*j*NX + 2*i + 1] *= 0.99;
+            total_flux[2*j*NX + 2*i] += gx;
+            total_flux[2*j*NX + 2*i + 1] += gy;
+//                         mgx = total_flux[2*j*NX + 2*i]/(double)(time+1);
+//                         mgy = total_flux[2*j*NX + 2*i + 1]/(double)(time+1);
+            mgx = total_flux[2*j*NX + 2*i];
+            mgy = total_flux[2*j*NX + 2*i + 1];
+//                         mgx = total_flux[2*j*NX + 2*i]/log((double)(time+2));
+//                         mgy = total_flux[2*j*NX + 2*i + 1]/log((double)(time+2));
+            mod = module2(mgx, mgy);
+            arg = argument(mgx, mgy);
+            if (arg < 0.0) arg += DPI;
+            value = arg/DPI;
+            color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb);
+            flux_factor = tanh(mod*E_SCALE);
+            for (k=0; k<3; k++) rgb[k] *= flux_factor;
+            break;
+        }
+    }
+    
+    return(value);
+}
+
+
+void draw_wave_profile(double *values, int size, int fade, double fade_value)
+/* draw a profile of the wave, if option DRAW_WAVE_PROFILE is active */
+{
+    int i, k;
+    double vmin, vmax, deltav, a, b, y;
+    static int imin, imax, jmin, jmax, jmid, d, first = 1;
+    static double deltaj, ymin;
+    
+    if (first)
+    {
+        imin = 100;
+        imax = NX - 250;
+        jmin = NY - 150;
+        jmax = NY - 50;
+        jmid = (jmin + jmax)/2;
+        jmid -= (jmid%size);
+        d = (jmax - jmin)/10 + 1;
+        deltaj = (double)(jmax - jmin - 2*d); 
+        ymin = (double)(jmin + d);
+        first = 0;
+    }
+    
+    vmin = 1.0e10;
+    vmax = -vmin;
+    for (i=imin; i<imax; i+=size)
+    {
+        if (values[i*NY+jmin] > vmax) vmax = values[i*NY+jmin];
+        if (values[i*NY+jmin] < vmin) vmin = values[i*NY+jmin];        
+    }
+    if ((vmin < 0.0)&&(vmax > 0.0))
+    {
+        if (vmax > -vmin) vmin = -vmax;
+        else if (vmax < -vmin) vmax = -vmin;
+    }
+//     printf("vmin = %.3lg, vmax = %.3lg\n", vmin, vmax);
+    deltav = vmax-vmin;
+    if (deltav == 0.0) deltav = 0.01;
+    a = deltaj/deltav;
+    b = ymin - a*vmin;
+
+    erase_area_ij(imin, jmin, imax, jmax);
+    
+    if (fade) glColor3f(fade_value, fade_value, fade_value);
+    else glColor3f(1.0, 1.0, 1.0);
+    glBegin(GL_LINE_STRIP);
+    for (i=imin; i<imax; i+=size)
+    {
+        y = a*values[i*NY+jmin] + b;
+        glVertex2d((double)i, y);
+    }
+    glEnd();
+    
+    glBegin(GL_LINE_LOOP);
+    glVertex2i(imin, jmin);
+    glVertex2i(imax, jmin);
+    glVertex2i(imax, jmax);
+    glVertex2i(imin, jmax);
+    glEnd();
+}
+
+
 void draw_wave_highres_palette(int size, double *phi[NX], double *psi[NX], double *total_energy[NX], double *total_flux, short int *xy_in[NX], double scale, int time, int plot, int palette, int fade, double fade_value)
 /* same as draw_wave_highres, but with color scheme option */
 {
     int i, j, k, iplus, iminus, jplus, jminus;
-    double rgb[3], xy[2], x1, y1, x2, y2, velocity, energy, gradientx2, gradienty2, arg, mod, flux_factor, gx, gy, mgx, mgy;
+    double value, rgb[3], xy[2], x1, y1, x2, y2, velocity, energy, gradientx2, gradienty2, arg, mod, flux_factor, gx, gy, mgx, mgy;
+//     double vmin, vmax, deltav;
     static double dtinverse = ((double)NX)/(COURANT*(XMAX-XMIN)), dx = (XMAX-XMIN)/((double)NX);
+    double *values;
+    
+    if (DRAW_WAVE_PROFILE) values = (double *)malloc(NX*NY*sizeof(double));
 
     glBegin(GL_QUADS);
     
@@ -869,86 +1028,10 @@ void draw_wave_highres_palette(int size, double *phi[NX], double *psi[NX], doubl
         {
             if ((TWOSPEEDS)||(xy_in[i][j]))
             {
-                switch (plot) {
-                    case (P_AMPLITUDE):
-                    {
-//                         /* make wave luminosity larger inside obstacles */
-//                         if (!(xy_in[i][j])) color_scheme_lum(COLOR_SCHEME, phi[i][j], scale, time, 0.7, rgb);
-//                         else 
-                        color_scheme_palette(COLOR_SCHEME, palette, phi[i][j], scale, time, rgb);
-                        break;
-                    }
-                    case (P_ENERGY):
-                    {
-                        energy = compute_energy(phi, psi, xy_in, i, j);
-                        /* adjust energy to color palette */
-                        if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, energy, scale, time, rgb);
-                        else color_scheme_palette(COLOR_SCHEME, palette, energy, scale, time, rgb);
-                        break;
-                    }
-                    case (P_MIXED):
-                    {
-                        if (j > NY/2) color_scheme_palette(COLOR_SCHEME, palette, phi[i][j], scale, time, rgb);
-                        else color_scheme_palette(COLOR_SCHEME, palette, compute_energy(phi, psi, xy_in, i, j), scale, time, rgb);
-                        break;
-                    }
-                    case (P_MEAN_ENERGY):
-                    {
-                        energy = compute_energy(phi, psi, xy_in, i, j);
-                        total_energy[i][j] += energy;
-                        if (COLOR_PALETTE >= COL_TURBO) 
-                            color_scheme_asym_palette(COLOR_SCHEME, palette, total_energy[i][j]/(double)(time+1), scale, time, rgb);
-                        else color_scheme_palette(COLOR_SCHEME, palette, total_energy[i][j]/(double)(time+1), scale, time, rgb);
-                        break;
-                    }
-                    case (P_LOG_ENERGY):
-                    {
-                        energy = compute_energy(phi, psi, xy_in, i, j);
-//                         energy = LOG_SHIFT + LOG_SCALE*log(energy);
-//                         if (energy < 0.0) energy = 0.0;
-                        color_scheme_palette(COLOR_SCHEME, palette, LOG_SHIFT + LOG_SCALE*log(energy), scale, time, rgb);
-                        break;
-                    }
-                    case (P_LOG_MEAN_ENERGY):
-                    {
-                        energy = compute_energy(phi, psi, xy_in, i, j);
-                        if (energy == 0.0) energy = 1.0e-20;
-                        total_energy[i][j] += energy;
-                        color_scheme_palette(COLOR_SCHEME, palette, LOG_SHIFT + LOG_SCALE*log(total_energy[i][j]/(double)(time+1)), scale, time, rgb);
-                        break;
-                    }
-                    case (P_ENERGY_FLUX):
-                    {
-                        compute_energy_flux(phi, psi, xy_in, i, j, &gx, &gy, &arg, &mod);
-//                         color_scheme_palette(C_ONEDIM_LINEAR, palette, arg/DPI, 1.0, 1, rgb);
-//                         flux_factor = tanh(mod*E_SCALE);
-//                         for (k=0; k<3; k++) rgb[k] *= flux_factor;
-                        color_scheme_asym_palette(COLOR_SCHEME, palette, mod*FLUX_SCALE, scale, time, rgb);
-                        break;
-                    }
-                    case (P_TOTAL_ENERGY_FLUX):
-                    {
-                        compute_energy_flux(phi, psi, xy_in, i, j, &gx, &gy, &arg, &mod);
-                        total_flux[2*j*NX + 2*i] *= 0.99;
-                        total_flux[2*j*NX + 2*i + 1] *= 0.99;
-                        total_flux[2*j*NX + 2*i] += gx;
-                        total_flux[2*j*NX + 2*i + 1] += gy;
-//                         mgx = total_flux[2*j*NX + 2*i]/(double)(time+1);
-//                         mgy = total_flux[2*j*NX + 2*i + 1]/(double)(time+1);
-                        mgx = total_flux[2*j*NX + 2*i];
-                        mgy = total_flux[2*j*NX + 2*i + 1];
-//                         mgx = total_flux[2*j*NX + 2*i]/log((double)(time+2));
-//                         mgy = total_flux[2*j*NX + 2*i + 1]/log((double)(time+2));
-                        mod = module2(mgx, mgy);
-                        arg = argument(mgx, mgy);
-                        if (arg < 0.0) arg += DPI;
-                        color_scheme_palette(C_ONEDIM_LINEAR, palette, arg/DPI, 1.0, 1, rgb);
-                        flux_factor = tanh(mod*E_SCALE);
-                        for (k=0; k<3; k++) rgb[k] *= flux_factor;
-                        break;
-                    }
-                    
-                }
+                value = wave_value(i, j, phi, psi, total_energy, total_flux, xy_in, scale, time, plot, palette, rgb);
+                
+                if (DRAW_WAVE_PROFILE) values[i*NY+j] = value;
+                
                 if (fade) for (k=0; k<3; k++) rgb[k] *= fade_value;
                 glColor3f(rgb[0], rgb[1], rgb[2]);
                 
@@ -960,6 +1043,12 @@ void draw_wave_highres_palette(int size, double *phi[NX], double *psi[NX], doubl
         }
 
     glEnd ();
+    
+    if (DRAW_WAVE_PROFILE) 
+    {
+        draw_wave_profile(values, size, fade, fade_value);
+        free(values);
+    }
 }
 
 /* modified function for "flattened" wave tables */

wave_comparison.c

@@ -43,7 +43,7 @@
 #include <omp.h>
 #include <time.h>
 
-#define MOVIE 1         /* set to 1 to generate movie */
+#define MOVIE 0         /* set to 1 to generate movie */
 #define DOUBLE_MOVIE 1  /* set to 1 to produce movies for wave height and energy simultaneously */
 #define SAVE_MEMORY 1   /* set to 1 to save memory when writing tiff images */
 #define NO_EXTRA_BUFFER_SWAP 1    /* some OS require one less buffer swap when recording images */
@@ -263,6 +263,8 @@
 // #define TWOSPEEDS 0          /* set to 1 to replace hardcore boundary by medium with different speed */
 #define WAVE_PACKET_SOURCE_TYPE 1       /* type of wave packet sources */
 #define N_WAVE_PACKETS 15               /* number of wave packets */
+#define OSCIL_LEFT_YSHIFT 0.0   /* y-dependence of left oscillation (for non-horizontal waves) */
+#define DRAW_WAVE_PROFILE 0     /* set to 1 to draw a profile of the wave */
 /* end of constants only used by sub_wave and sub_maze */
 
 

wave_energy.c

@@ -219,8 +219,16 @@
 #define IOR 7               /* choice of index of refraction, see list in global_pdes.c */
 #define IOR_TOTAL_TURNS 1.5 /* total angle of rotation for IOR_PERIODIC_WELLS_ROTATING */
 #define MANDEL_IOR_SCALE -0.05   /* parameter controlling dependence of IoR on Mandelbrot escape speed */
+#define COURANT 0.04       /* Courant number */
+#define COURANTB 0.0       /* Courant number in medium B */
+#define INITIAL_AMP 0.5            /* amplitude of initial condition */
+#define INITIAL_VARIANCE 0.0003    /* variance of initial condition */
+#define INITIAL_WAVELENGTH  0.015  /* wavelength of initial condition */
+#define TWOSPEEDS 0          /* set to 1 to replace hardcore boundary by medium with different speed */
 #define WAVE_PACKET_SOURCE_TYPE 1       /* type of wave packet sources */
 #define N_WAVE_PACKETS 15               /* number of wave packets */
+#define OSCIL_LEFT_YSHIFT 0.0   /* y-dependence of left oscillation (for non-horizontal waves) */
+#define DRAW_WAVE_PROFILE 0     /* set to 1 to draw a profile of the wave */
 /* end of constants only used by sub_wave and sub_maze */
 
 #include "global_pdes.c"        /* constants and global variables */

wave_sphere.c

@@ -0,0 +1,944 @@
+/*********************************************************************************/
+/*                                                                               */
+/*  Animation of wave equation on a sphere                                       */
+/*                                                                               */
+/*  N. Berglund, july 2023                                                       */
+/*                                                                               */
+/*  Feel free to reuse, but if doing so it would be nice to drop a               */
+/*  line to nils.berglund@univ-orleans.fr - Thanks!                              */
+/*                                                                               */
+/*  compile with                                                                 */
+/*  gcc -o wave_sphere wave_sphere.c                                             */
+/* -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lXmu -lglut -O3 -fopenmp        */
+/*                                                                               */
+/*  OMP acceleration may be more effective after executing                       */
+/*  export OMP_NUM_THREADS=2 in the shell before running the program             */
+/*                                                                               */
+/*  To make a video, set MOVIE to 1 and create subfolder tif_wave                */
+/*  It may be possible to increase parameter PAUSE                               */
+/*                                                                               */
+/*  create movie using                                                           */
+/*  ffmpeg -i wave.%05d.tif -vcodec libx264 wave.mp4                             */
+/*                                                                               */
+/*********************************************************************************/
+
+/*********************************************************************************/
+/*                                                                               */
+/* NB: The algorithm used to simulate the wave equation is highly paralellizable */
+/* One could make it much faster by using a GPU                                  */
+/*                                                                               */
+/*********************************************************************************/
+
+#include <math.h>
+#include <string.h>
+#include <GL/glut.h>
+#include <GL/glu.h>
+#include <unistd.h>
+#include <sys/types.h>
+#include <tiffio.h>     /* Sam Leffler's libtiff library. */
+#include <omp.h>
+#include <time.h>
+
+#define MOVIE 0         /* set to 1 to generate movie */
+#define DOUBLE_MOVIE 1  /* set to 1 to produce movies for wave height and energy simultaneously */
+#define SAVE_MEMORY 1   /* set to 1 to save memory when writing tiff images */
+#define NO_EXTRA_BUFFER_SWAP 1    /* some OS require one less buffer swap when recording images */
+
+/* General geometrical parameters */
+
+#define WINWIDTH 	1920  /* window width */
+#define WINHEIGHT 	1150  /* window height */
+#define NX 2560          /* number of grid points on x axis */
+#define NY 1280           /* number of grid points on y axis */
+
+#define DPOLE 20         /* safety distance to poles */
+#define SMOOTHPOLE 0.1     /* smoothing coefficient at poles */
+
+#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 */
+
+#define JULIA_SCALE 0.25     /* scaling for Julia sets */
+#define JULIA_ROT 90.0       /* rotation of Julia set, in degrees */
+#define JULIA_RE -0.77145    
+#define JULIA_IM -0.10295    /* parameters for Julia sets */
+
+/* Choice of the billiard table */
+
+#define B_DOMAIN 84         /* choice of domain shape, see list in global_pdes.c */
+
+#define CIRCLE_PATTERN 33   /* pattern of circles or polygons, see list in global_pdes.c */
+
+#define COMPARISON 0        /* set to 1 to compare two different patterns */
+#define B_DOMAIN_B 20       /* second domain shape, for comparisons */
+#define CIRCLE_PATTERN_B 0  /* second pattern of circles or polygons */
+
+#define VARIABLE_IOR 0      /* set to 1 for a variable index of refraction */
+#define IOR 9               /* choice of index of refraction, see list in global_pdes.c */
+#define IOR_TOTAL_TURNS 1.0 /* total angle of rotation for IOR_PERIODIC_WELLS_ROTATING */
+#define MANDEL_IOR_SCALE -0.05   /* parameter controlling dependence of IoR on Mandelbrot escape speed */
+
+#define P_PERCOL 0.25       /* probability of having a circle in C_RAND_PERCOL arrangement */
+#define NPOISSON 1000        /* number of points for Poisson C_RAND_POISSON arrangement */
+#define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */
+
+#define LAMBDA 0.75	    /* parameter controlling the dimensions of domain */
+#define MU 0.1             /* parameter controlling the dimensions of domain */
+#define NPOLY 6             /* number of sides of polygon */
+#define APOLY 0.0           /* angle by which to turn polygon, in units of Pi/2 */ 
+#define MDEPTH 7            /* depth of computation of Menger gasket */
+#define MRATIO 3            /* ratio defining Menger gasket */
+#define MANDELLEVEL 2000    /* iteration level for Mandelbrot set */
+#define MANDELLIMIT 20.0    /* limit value for approximation of Mandelbrot set */
+#define FOCI 1              /* set to 1 to draw focal points of ellipse */
+#define NGRIDX 30            /* number of grid point for grid of disks */
+#define NGRIDY 18            /* number of grid point for grid of disks */
+
+#define X_SHOOTER -0.2
+#define Y_SHOOTER -0.6
+#define X_TARGET 0.4
+#define Y_TARGET 0.7        /* shooter and target positions in laser fight */
+
+#define ISO_XSHIFT_LEFT -2.9
+#define ISO_XSHIFT_RIGHT 1.4
+#define ISO_YSHIFT_LEFT -0.15
+#define ISO_YSHIFT_RIGHT -0.15 
+#define ISO_SCALE 0.5           /* coordinates for isospectral billiards */
+
+
+/* You can add more billiard tables by adapting the functions */
+/* xy_in_billiard and draw_billiard below */
+
+/* Physical parameters of wave equation */
+
+#define TWOSPEEDS 0          /* set to 1 to replace hardcore boundary by medium with different speed */
+#define OSCILLATE_LEFT 0     /* set to 1 to add oscilating boundary condition on the left */
+#define OSCILLATE_TOPBOT 0   /* set to 1 to enforce a planar wave on top and bottom boundary */
+#define OSCILLATION_SCHEDULE 3  /* oscillation schedule, see list in global_pdes.c */
+
+#define OMEGA 0.001       /* frequency of periodic excitation */
+#define AMPLITUDE 0.8     /* amplitude of periodic excitation */ 
+#define ACHIRP 0.2        /* acceleration coefficient in chirp */
+#define DAMPING 0.0       /* damping of periodic excitation */
+#define COURANT 0.05       /* Courant number */
+#define COURANTB 0.01      /* Courant number in medium B */
+#define GAMMA 0.0          /* damping factor in wave equation */
+#define GAMMAB 1.0e-6         /* damping factor in wave equation */
+#define GAMMA_SIDES 1.0e-4      /* damping factor on boundary */
+#define GAMMA_TOPBOT 1.0e-7     /* damping factor on boundary */
+#define KAPPA 0.0           /* "elasticity" term enforcing oscillations */
+#define KAPPA_SIDES 5.0e-4  /* "elasticity" term on absorbing boundary */
+#define KAPPA_TOPBOT 0.0    /* "elasticity" term on absorbing boundary */
+#define OSCIL_LEFT_YSHIFT 0.0   /* y-dependence of left oscillation (for non-horizontal waves) */
+/* The Courant number is given by c*DT/DX, where DT is the time step and DX the lattice spacing */
+/* The physical damping coefficient is given by GAMMA/(DT)^2 */
+/* Increasing COURANT speeds up the simulation, but decreases accuracy */
+/* For similar wave forms, COURANT^2*GAMMA should be kept constant */
+
+#define ADD_OSCILLATING_SOURCE 0        /* set to 1 to add an oscillating wave source */
+#define OSCILLATING_SOURCE_PERIOD 30    /* period of oscillating source */
+#define ALTERNATE_OSCILLATING_SOURCE 1  /* set to 1 to alternate sign of oscillating source */
+
+#define ADD_WAVE_PACKET_SOURCES 0       /* set to 1 to add several sources emitting wave packets */
+#define WAVE_PACKET_SOURCE_TYPE 1       /* type of wave packet sources */
+#define N_WAVE_PACKETS 15               /* number of wave packets */
+#define WAVE_PACKET_RADIUS 20            /* radius of wave packets */
+
+/* Boundary conditions, see list in global_pdes.c  */
+
+#define B_COND 2
+
+#define PRECOMPUTE_BC 0     /* set to 1 to compute neighbours for Laplacian in advance */
+
+/* Parameters for length and speed of simulation */
+
+#define NSTEPS 3600       /* number of frames of movie */
+#define NVID 4            /* number of iterations between images displayed on screen */
+#define NSEG 1000          /* number of segments of boundary */
+#define INITIAL_TIME 0      /* time after which to start saving frames */
+#define BOUNDARY_WIDTH 2    /* width of billiard boundary */
+#define PRINT_SPEED 0       /* set to 1 to print speed of moving source */
+
+#define PAUSE 100       /* number of frames after which to pause */
+#define PSLEEP 3         /* sleep time during pause */
+#define SLEEP1  1        /* initial sleeping time */
+#define SLEEP2  1        /* final sleeping time */
+#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 */
+#define ROTATE_VIEW_WHILE_FADE 1    /* set to 1 to keep rotating viewpoint during fade */
+
+/* Parameters of initial condition */
+
+#define INITIAL_AMP 0.75            /* amplitude of initial condition */
+#define INITIAL_VARIANCE 0.0005  /* variance of initial condition */
+#define INITIAL_WAVELENGTH  0.025  /* wavelength of initial condition */
+
+/* Plot type, see list in global_pdes.c  */
+
+#define ZPLOT 103     /* wave height */
+#define CPLOT 103     /* color scheme */
+
+#define ZPLOT_B 108 
+#define CPLOT_B 108        /* plot type for second movie */
+
+#define CHANGE_LUMINOSITY 1     /* set to 1 to let luminosity depend on energy flux intensity */
+#define FLUX_WINDOW 30          /* size of averaging window of flux intensity */
+#define AMPLITUDE_HIGH_RES 1    /* set to 1 to increase resolution of plot */
+#define SHADE_3D 1              /* set to 1 to change luminosity according to normal vector */
+#define SHADE_2D 0              /* set to 1 to change luminosity according to normal vector to plane */
+#define SHADE_WAVE 1            /* set to 1 to have luminosity depend on wave height */
+#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 0         /* set to 1 to draw boundary */
+#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 SHADE_SCALE_2D 10.0     /* controls "depth" of 2D shading */
+#define COS_LIGHT_MIN 0.0       /* controls angle-dependence of 2D shading */
+#define COS_LIGHT_MAX 0.8       /* controls angle-dependence of 2D shading */
+
+#define PLOT_SCALE_ENERGY 0.4          /* vertical scaling in energy plot */
+#define PLOT_SCALE_LOG_ENERGY 0.5       /* vertical scaling in log energy plot */
+
+/* 3D representation */
+
+#define REPRESENTATION_3D 1     /* choice of 3D representation */ 
+#define PLOT_2D 0               /* switch to 2D representation, equirectangular projection */
+#define PHISHIFT 0.0            /* shift of phi in 2D plot (in degrees) */
+
+#define REP_AXO_3D 0        /* linear projection (axonometry) */
+#define REP_PROJ_3D 1       /* projection on plane orthogonal to observer line of sight */
+
+#define ROTATE_VIEW 1       /* set to 1 to rotate position of observer */
+#define ROTATE_ANGLE -360.0   /* total angle of rotation during simulation */
+
+#define VIEWPOINT_TRAJ 1    /* type of viewpoint trajectory */
+
+/* Color schemes */
+
+#define COLOR_PALETTE 11      /* Color palette, see list in global_pdes.c  */
+#define COLOR_PALETTE_B 16     /* Color palette, see list in global_pdes.c  */
+
+#define BLACK 1          /* background */
+#define COLOR_OUT_R 1.0    /* color outside domain */
+#define COLOR_OUT_G 1.0    
+#define COLOR_OUT_B 1.0    
+
+#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 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 100.0       /* scaling factor for energy representation */
+#define LOG_SCALE 0.75     /* scaling factor for energy log representation */
+#define LOG_SHIFT 0.5      /* shift of colors on log scale */
+#define LOG_ENERGY_FLOOR -10.0    /* floor value for log of (total) energy */
+#define LOG_MEAN_ENERGY_SHIFT 1.0   /* additional shift for log of mean energy */
+#define FLUX_SCALE 600.0    /* scaling factor for energy flux representation */
+#define FLUX_CSCALE 5.0      /* scaling factor for color in energy flux representation */
+#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 */
+#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 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 NXMAZE 8      /* width of maze */
+#define NYMAZE 32      /* height of maze */
+#define MAZE_MAX_NGBH 5     /* max number of neighbours of maze cell */
+#define RAND_SHIFT 5        /* seed of random number generator */
+#define MAZE_XSHIFT 0.0     /* horizontal shift of maze */
+#define MAZE_WIDTH 0.02     /* half width of maze walls */
+
+#define DRAW_COLOR_SCHEME 1       /* set to 1 to plot the color scheme */
+#define COLORBAR_RANGE 3.0      /* scale of color scheme bar */
+#define COLORBAR_RANGE_B 2.0    /* scale of color scheme bar for 2nd part */
+#define ROTATE_COLOR_SCHEME 0     /* set to 1 to draw color scheme horizontally */
+#define CIRC_COLORBAR 0         /* set to 1 to draw circular color scheme */
+#define CIRC_COLORBAR_B 0       /* set to 1 to draw circular color scheme */
+
+#define DRAW_WAVE_PROFILE 0     /* set to 1 to draw a profile of the wave */
+#define SAVE_TIME_SERIES 0      /* set to 1 to save wave time series at a point */
+
+#define ADD_POTENTIAL 0         /* set to 1 to add potential to z coordinate */
+#define POTENTIAL 10
+#define POT_FACT 20.0
+/* end of constants only used by sub_wave and sub_maze */
+
+
+/* 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 */
+
+/* Parameters controlling 3D projection */
+
+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] = {6.0, 8.0, 2.5};    /* location of observer for REP_PROJ_3D representation */ 
+int reset_view = 0;         /* switch to reset 3D view parameters (for option ROTATE_VIEW) */
+
+#define RSCALE 0.01             /* scaling factor of radial component */
+#define RMAX 10.0               /* max value of radial component */
+#define Z_SCALING_FACTOR 0.8    /* 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 COS_VISIBLE -0.75        /* limit on cosine of normal to shown facets */
+
+#include "global_pdes.c"        /* constants and global variables */
+#include "global_3d.c"          /* constants and global variables */
+#include "sub_maze.c"           /* support for generating mazes */
+#include "sub_wave.c"           /* common functions for wave_billiard, heat and schrodinger */
+#include "wave_common.c"        /* common functions for wave_billiard, wave_comparison, etc */
+
+#include "sub_wave_3d.c"        /* graphical functions specific to wave_3d */
+#include "sub_sphere.c"         /* graphical functions specific to wave_sphere */
+
+FILE *time_series_left, *time_series_right, *image_file;
+
+double courant2, courantb2;  /* Courant parameters squared */
+
+
+void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out[NX*NY], 
+                      short int xy_in[NX*NY], double tc[NX*NY], double tcc[NX*NY], double tgamma[NX*NY])
+/* time step of field evolution */
+/* phi is value of field at time t, psi at time t-1 */
+/* this version of the function has been rewritten in order to minimize the number of if-branches */
+{
+    int i, j, iplus, iminus, jplus, jminus, jtop, jbot;
+    double delta, x, y, c, cc, gamma, sintheta, cottheta, invstheta, sum, avrg;
+    static long time = 0;
+    static short int first = 1;
+    static double dphi, dtheta, cphiphi, ctheta;
+    
+    if (first)
+    {
+        dphi = DPI/(double)NX;
+        dtheta = PI/(double)NY;
+        cphiphi = dphi*dphi/(dtheta*dtheta);
+        ctheta = dphi*dphi/(2.0*dtheta);
+        
+        printf("dphi = %.5lg, dtheta = %.5lg, cphiphi = %.5lg, ctheta = %.5lg\n", dphi, dtheta, cphiphi, ctheta); 
+        
+        first = 0;
+    }
+    
+    time++;
+    
+    #pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta,x,y)
+    /* evolution in the bulk */
+    for (j=DPOLE; j<NY-DPOLE; j++){
+        sintheta = sin(j*dtheta);
+//         invstheta = 1.0/(sintheta*sintheta);
+        invstheta = 1.0/(sintheta*sintheta + SMOOTHPOLE*SMOOTHPOLE);
+//         cottheta = ctheta*cos(j*dtheta)/sintheta;
+        cottheta = ctheta*cos(j*dtheta)/(sintheta + SMOOTHPOLE);
+        
+        for (i=1; i<NX-1; i++){
+            if ((TWOSPEEDS)||(xy_in[i*NY+j] != 0)){
+                x = phi_in[i*NY+j];
+		y = psi_in[i*NY+j];
+                
+                /* discretized Laplacian */
+                /* 2nd phi derivative */
+                delta = invstheta*(phi_in[(i+1)*NY+j] + phi_in[(i-1)*NY+j] - 2.0*x);
+                
+                /* 2nd theta derivative */
+                delta += cphiphi*(phi_in[i*NY+j+1] + phi_in[i*NY+j-1] - 2.0*x);
+                
+                /* first theta derivative */
+                delta += cottheta*(phi_in[i*NY+j+1] - phi_in[i*NY+j-1]);
+
+                /* evolve phi */
+                phi_out[i*NY+j] = -y + 2*x + tcc[i*NY+j]*delta - KAPPA*x - tgamma[i*NY+j]*(x-y);
+            }
+        }
+    }
+    
+    /* evolution at longitude zero */
+    for (j=DPOLE; j<NY-DPOLE; j++){
+        sintheta = sin(j*dtheta);
+        invstheta = 1.0/(sintheta*sintheta);
+        cottheta = ctheta*cos(j*dtheta)/sintheta;
+        
+        /* i = 0 */
+        if ((TWOSPEEDS)||(xy_in[j] != 0)){
+            x = phi_in[j];
+            y = psi_in[j];
+                
+            /* discretized Laplacian */
+            /* 2nd phi derivative */
+            delta = invstheta*(phi_in[NY+j] + phi_in[(NX-1)*NY+j] - 2.0*x);
+                
+            /* 2nd theta derivative */
+            delta += cphiphi*(phi_in[j+1] + phi_in[j-1] - 2.0*x);
+                
+            /* first theta derivative */
+            delta += cottheta*(phi_in[j+1] - phi_in[j-1]);
+
+            /* evolve phi */
+            phi_out[j] = -y + 2*x + tcc[j]*delta - KAPPA*x - tgamma[j]*(x-y);
+        }
+        
+        /* i = NX-1 */
+        if ((TWOSPEEDS)||(xy_in[(NX-1)*NY+j] != 0)){
+            x = phi_in[(NX-1)*NY+j];
+            y = psi_in[(NX-1)*NY+j];
+                
+            /* discretized Laplacian */
+            /* 2nd phi derivative */
+            delta = invstheta*(phi_in[j] + phi_in[(NX-2)*NY+j] - 2.0*x);
+                
+            /* 2nd theta derivative */
+            delta += cphiphi*(phi_in[(NX-1)*NY+j+1] + phi_in[(NX-1)*NY+j-1] - 2.0*x);
+                
+            /* first theta derivative */
+            delta += cottheta*(phi_in[(NX-1)*NY+j+1] - phi_in[(NX-1)*NY+j-1]);
+
+            /* evolve phi */
+            phi_out[(NX-1)*NY+j] = -y + 2*x + tcc[(NX-1)*NY+j]*delta - KAPPA*x - tgamma[(NX-1)*NY+j]*(x-y);
+        }
+    }
+    
+    /* compute average at north pole */
+    sum = 0.0;
+    for (i=0; i<NX; i++) sum += phi_out[i*NY + DPOLE];
+    avrg = sum/(double)NX;
+    for (i=0; i<NX; i++) for (j=0; j<DPOLE; j++)
+        phi_out[i*NY + j] = avrg;
+//     {
+//         x = phi_in[(NX-1)*NY+j];
+//         y = psi_in[(NX-1)*NY+j];
+//         phi_out[(NX-1)*NY+j] = -y + 2*x + tcc[(NX-1)*NY+j]*avrg - KAPPA*x - tgamma[(NX-1)*NY+j]*(x-y);
+// //         
+//     }
+    
+    /* compute average at south pole */
+    sum = 0.0;
+    for (i=0; i<NX; i++) sum += phi_out[i*NY + NY-1-DPOLE];
+    avrg = sum/(double)NX;
+    for (i=0; i<NX; i++) for (j=NY-DPOLE; j<NY; j++) 
+        phi_out[i*NY + j] = avrg;
+//     {
+//         x = phi_in[(NX-1)*NY+j];
+//         y = psi_in[(NX-1)*NY+j];
+//         phi_out[(NX-1)*NY+j] = -y + 2*x + tcc[(NX-1)*NY+j]*avrg - KAPPA*x - tgamma[(NX-1)*NY+j]*(x-y);
+// //         
+//     }
+    
+//     /* north pole, j = 0 */
+//     if ((TWOSPEEDS)||(xy_in[0] != 0)){
+//         x = phi_in[0];
+//         y = psi_in[0];
+//         
+//         /* discretized Laplacian */
+//         delta = cphiphi*(phi_in[1] + phi_in[(NX/4)*NY+1] + phi_in[(NX/2)*NY+1] + phi_in[(3*NX/4)*NY+1] - 4.0*x);
+//         
+//         /* evolve phi */
+//         phi_out[0] = -y + 2*x + tcc[0]*delta - KAPPA*x - tgamma[0]*(x-y);
+//         
+//         /* set same values for all i */
+//         for (i=1; i<NX; i++) phi_out[i*NY] = phi_out[0];
+//     }
+//     
+//     /* south pole, j = NY-1 */
+//     if ((TWOSPEEDS)||(xy_in[NY-1] != 0)){
+//         x = phi_in[NY-1];
+//         y = psi_in[NY-1];
+//         
+//         /* discretized Laplacian */
+//         delta = cphiphi*(phi_in[NY-2] + phi_in[(NX/4)*NY+NY-2] + phi_in[(NX/2)*NY+NY-2] + phi_in[(3*NX/4)*NY+NY-2] - 4.0*x);
+//         
+//         /* evolve phi */
+//         phi_out[NY-1] = -y + 2*x + tcc[0]*delta - KAPPA*x - tgamma[0]*(x-y);
+//         
+//         /* set same values for all i */
+//         for (i=1; i<NX; i++) phi_out[i*NY+NY-1] = phi_out[NY-1];
+//     }
+    
+    /* for debugging purposes/if there is a risk of blow-up */
+    if (FLOOR) for (i=0; i<NX; i++){
+        for (j=0; j<NY; j++){
+            if (xy_in[i*NY+j] != 0) 
+            {
+                if (phi_out[i*NY+j] > VMAX) phi_out[i*NY+j] = VMAX;
+                if (phi_out[i*NY+j] < -VMAX) phi_out[i*NY+j] = -VMAX;
+            }
+        }
+    }
+}
+
+
+void evolve_wave(double phi[NX*NY], double psi[NX*NY], double tmp[NX*NY], short int xy_in[NX*NY],
+    double tc[NX*NY], double tcc[NX*NY], double tgamma[NX*NY], t_laplacian laplace[NX*NY], t_laplacian laplace1[NX*NY], 
+    t_laplacian laplace2[NX*NY])
+/* time step of field evolution */
+/* phi is value of field at time t, psi at time t-1 */
+{
+    evolve_wave_half(phi, psi, tmp, xy_in, tc, tcc, tgamma);
+    evolve_wave_half(tmp, phi, psi, xy_in, tc, tcc, tgamma);
+    evolve_wave_half(psi, tmp, phi, xy_in, tc, tcc, tgamma);
+}
+
+
+void draw_color_bar_palette(int plot, double range, int palette, int circular, int fade, double fade_value)
+{
+//     double width = 0.2;
+    double width = 0.14;
+//     double width = 0.2;
+    
+    width *= (double)NX/(double)WINWIDTH;
+    
+    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 if (circular)
+        draw_circular_color_scheme_palette_3d(XMAX - 2.0*width, YMIN + 2.0*width, 1.5*width, 1.3*width, plot, -range, range, palette, fade, fade_value);
+    else 
+        draw_color_scheme_palette_3d(XMAX - 1.5*width, YMIN + 0.1, XMAX - 0.5*width, YMAX - 0.1, plot, -range, range, palette, fade, fade_value);
+}
+
+void viewpoint_schedule(int i)
+/* change position of observer */
+{
+    int j;
+    double angle, ca, sa, r1;
+    static double observer_initial[3], r, ratio;
+    static int first = 1;
+    
+    if (first)
+    {
+        for (j=0; j<3; j++) observer_initial[j] = observer[j];
+        r1 = observer[0]*observer[0] + observer[1]*observer[1];
+        r = sqrt(r1 + observer[2]*observer[2]);
+        ratio = r/sqrt(r1);
+        first = 0;
+    }
+    
+    angle = (ROTATE_ANGLE*DPI/360.0)*(double)i/(double)NSTEPS;
+    ca = cos(angle);
+    sa = sin(angle);
+    switch (VIEWPOINT_TRAJ)
+    {
+        case (VP_HORIZONTAL):
+        {
+            observer[0] = ca*observer_initial[0] - sa*observer_initial[1];
+            observer[1] = sa*observer_initial[0] + ca*observer_initial[1];
+            break;
+        }
+        case (VP_ORBIT):
+        {
+            observer[0] = ca*observer_initial[0] - sa*observer_initial[1]*ratio;
+            observer[1] = ca*observer_initial[1] + sa*observer_initial[0]*ratio;
+            observer[2] = ca*observer_initial[2];
+            break;
+        }
+    }
+    
+    printf("Angle %.3lg, Observer position (%.3lg, %.3lg, %.3lg)\n", angle, observer[0], observer[1], observer[2]);
+}
+
+
+void animation()
+{
+    double time, scale, ratio, startleft[2], startright[2], sign = 1.0, r2, xy[2], fade_value, yshift, speed = 0.0, a, b, c, angle = 0.0, lambda1, y, x1, sign1, omega, phase_shift, theta, amp; 
+    double *phi, *psi, *tmp, *color_scale, *tc, *tcc, *tgamma;
+//     double *total_energy;
+    short int *xy_in;
+    int i, j, s, sample_left[2], sample_right[2], period = 0, fade, source_counter = 0, k, p, q;
+    static int counter = 0, first_source = 1;
+    long int wave_value;
+    t_wave *wave;
+    t_laplacian *laplace, *laplace1, *laplace2;
+    t_wave_source wave_source[25];
+    t_wave_sphere *wsphere;
+    
+    if (SAVE_TIME_SERIES)
+    {
+        time_series_left = fopen("wave_left.dat", "w");
+        time_series_right = fopen("wave_right.dat", "w");
+    }
+
+    /* Since NX and NY are big, it seemed wiser to use some memory allocation here */
+    xy_in = (short int *)malloc(NX*NY*sizeof(short int));
+    phi = (double *)malloc(NX*NY*sizeof(double));
+    psi = (double *)malloc(NX*NY*sizeof(double));
+    tmp = (double *)malloc(NX*NY*sizeof(double));
+//     total_energy = (double *)malloc(NX*NY*sizeof(double));
+    color_scale = (double *)malloc(NX*NY*sizeof(double));
+    tc = (double *)malloc(NX*NY*sizeof(double));
+    tcc = (double *)malloc(NX*NY*sizeof(double));
+    tgamma = (double *)malloc(NX*NY*sizeof(double));
+    
+    wave = (t_wave *)malloc(NX*NY*sizeof(t_wave));
+    wsphere = (t_wave_sphere *)malloc(NX*NY*sizeof(t_wave_sphere));
+    
+    laplace = (t_laplacian *)malloc(NX*NY*sizeof(t_laplacian));
+    laplace1 = (t_laplacian *)malloc(NX*NY*sizeof(t_laplacian));
+    laplace2 = (t_laplacian *)malloc(NX*NY*sizeof(t_laplacian));
+    
+    /* initialise positions and radii of circles */
+    if (COMPARISON)
+    {
+        if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN == D_CIRCLES_IN_RECT))
+            ncircles = init_circle_config_pattern(circles, CIRCLE_PATTERN);
+        else if (B_DOMAIN == D_POLYGONS) ncircles = init_polygon_config_pattern(polygons, CIRCLE_PATTERN);
+        
+        if ((B_DOMAIN_B == D_CIRCLES)||(B_DOMAIN_B == D_CIRCLES_IN_RECT))
+            ncircles_b = init_circle_config_pattern(circles_b, CIRCLE_PATTERN_B);
+        else if (B_DOMAIN_B == D_POLYGONS) ncircles_b = init_polygon_config_pattern(polygons_b, CIRCLE_PATTERN_B);
+        /* TO DO: adapt to different polygon patterns */
+    }
+    else
+    {
+        if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN == D_CIRCLES_IN_RECT)) ncircles = init_circle_config(circles);
+        else if (B_DOMAIN == D_POLYGONS) ncircles = init_polygon_config(polygons);
+    }
+    
+    printf("Polygons initialized\n");
+    
+    /* initialise polyline for von Koch and similar domains */
+    npolyline = init_polyline(MDEPTH, polyline);
+    for (i=0; i<npolyline; i++) printf("vertex %i: (%.3f, %.3f)\n", i, polyline[i].x, polyline[i].y);
+    if (COMPARISON) npolyline_b = init_polyline(MDEPTH, polyline);
+    
+    npolyrect = init_polyrect(polyrect);
+    for (i=0; i<npolyrect; i++) printf("polyrect vertex %i: (%.3f, %.3f) - (%.3f, %.3f)\n", i, polyrect[i].x1, polyrect[i].y1, polyrect[i].x2, polyrect[i].y2);
+
+    init_polyrect_arc(polyrectrot, polyarc, &npolyrect_rot, &npolyarc);
+    printf("Rotated rectangles and arcs initialized\n");
+    printf("%i rotated rectangles, %i arcs\n", npolyrect_rot, npolyarc);
+    
+    if ((B_DOMAIN == D_SPHERE_CIRCLES)||(B_DOMAIN_B == D_SPHERE_CIRCLES))
+    {
+        ncircles = init_circle_sphere(circ_sphere, CIRCLE_PATTERN); 
+    }
+    
+    courant2 = COURANT*COURANT;
+    courantb2 = COURANTB*COURANTB;
+    c = COURANT*(XMAX - XMIN)/(double)NX;
+//     a = 0.015;
+//     b = 0.0003;
+//     a = 0.04;
+//     b = 0.0018;
+    a = 0.05;
+    b = 0.0016;
+    
+    /* initialize color scale, for option RESCALE_COLOR_IN_CENTER */
+    if (RESCALE_COLOR_IN_CENTER)
+    {
+        for (i=0; i<NX; i++)
+            for (j=0; j<NY; j++)
+            {
+                ij_to_xy(i, j, xy);
+                r2 = xy[0]*xy[0] + xy[1]*xy[1];
+                color_scale[i*NY+j] = 1.0 - exp(-4.0*r2/LAMBDA*LAMBDA);
+            }
+    }
+
+    /* initialize wave with a drop at one point, zero elsewhere */
+//     init_circular_wave(0.0, -LAMBDA, phi, psi, xy_in);
+
+    
+    init_wave_fields(wave);
+    
+    init_wave_sphere(wsphere); 
+    
+    /* initialize total energy table - no longer needed */
+//     if ((ZPLOT == P_MEAN_ENERGY)||(ZPLOT_B == P_MEAN_ENERGY)||(ZPLOT == P_LOG_MEAN_ENERGY)||(ZPLOT_B == P_LOG_MEAN_ENERGY))
+//         for (i=0; i<NX; i++)
+//             for (j=0; j<NY; j++) 
+//                 total_energy[i*NY+j] = 0.0;
+    
+    ratio = (XMAX - XMIN)/8.4;  /* for Tokarsky billiard */
+    
+
+//     init_circular_wave_mod(polyline[85].x, polyline[85].y, phi, psi, xy_in);
+//     init_circular_wave_mod(LAMBDA*cos(APOLY*PID), LAMBDA*sin(APOLY*PID), phi, psi, xy_in);
+    lambda1 = LAMBDA;
+    angle = DPI/(double)NPOLY;
+//     init_circular_wave_mod(lambda1*cos(0.5*angle), lambda1*sin(0.5*angle), phi, psi, xy_in);
+//     for (j=1; j<NPOLY; j++)
+//         add_circular_wave_mod(1.0, lambda1*cos(((double)j+0.5)*angle), lambda1*sin(((double)j+0.5)*angle), phi, psi, xy_in);
+        
+    init_circular_wave_sphere(0.7, 0.5, phi, psi, xy_in, wsphere);
+//     init_wave_flat_sphere(phi, psi, xy_in, wsphere);
+//     init_circular_wave_sphere(0.25 + PID, 0.0, phi, psi, xy_in, wsphere);
+//     add_circular_wave_sphere(-1.0, 0.25 + 3.0*PID, 0.0, phi, psi, xy_in, wsphere);
+
+    
+//     printf("Wave initialized\n");
+
+    /* initialize table of wave speeds/dissipation */
+    init_speed_dissipation(xy_in, tc, tcc, tgamma);
+    
+    /* initialze potential to add to z coordinate */ 
+    if (ADD_POTENTIAL)
+    {
+        if (POTENTIAL == POT_IOR) 
+            for (i=0; i<NX*NY; i++)
+                wave[i].potential = &tcc[i];
+    }
+    
+    init_zfield(phi, psi, xy_in, ZPLOT, wave, 0);
+    init_cfield(phi, psi, xy_in, CPLOT, wave, 0);
+    
+    if (DOUBLE_MOVIE)
+    {
+        init_zfield(phi, psi, xy_in, ZPLOT_B, wave, 1);
+        init_cfield(phi, psi, xy_in, CPLOT_B, wave, 1);
+    }
+
+    blank();
+    glColor3f(0.0, 0.0, 0.0);
+    draw_wave_sphere(0, phi, psi, xy_in, wave, wsphere, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1);
+//     draw_billiard();
+    
+    
+    if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, 0, 1.0);
+
+    glutSwapBuffers();
+
+
+
+    sleep(SLEEP1);
+
+    for (i=0; i<=INITIAL_TIME + NSTEPS; i++)
+    {
+        global_time++;
+        
+	//printf("%d\n",i);
+        /* compute the variance of the field to adjust color scheme */
+        /* the color depends on the field divided by sqrt(1 + variance) */
+        if (SCALE)
+        {
+            scale = sqrt(1.0 + compute_variance_mod(phi,psi, xy_in));
+//             printf("Scaling factor: %5lg\n", scale);
+        }
+        else scale = 1.0;
+        
+        if (ROTATE_VIEW) 
+        {
+            viewpoint_schedule(i - INITIAL_TIME);
+            reset_view = 1;
+        }
+        
+        draw_wave_sphere(0, phi, psi, xy_in, wave, wsphere, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1);
+        for (j=0; j<NVID; j++) 
+        {
+            evolve_wave(phi, psi, tmp, xy_in, tc, tcc, tgamma, laplace, laplace1, laplace2);
+            if (SAVE_TIME_SERIES)
+            {
+                wave_value = (long int)(phi[sample_left[0]*NY+sample_left[1]]*1.0e16);
+                fprintf(time_series_left, "%019ld\n", wave_value);
+                wave_value = (long int)(phi[sample_right[0]*NY+sample_right[1]]*1.0e16);
+                fprintf(time_series_right, "%019ld\n", wave_value);
+                if ((j == 0)&&(i%10 == 0)) printf("Frame %i of %i\n", i, NSTEPS);
+//                 fprintf(time_series_right, "%.15f\n", phi[sample_right[0]][sample_right[1]]);
+            }
+//             if (i % 10 == 9) oscillate_linear_wave(0.2*scale, 0.15*(double)(i*NVID + j), -1.5, YMIN, -1.5, YMAX, phi, psi);
+        }
+        
+//         draw_billiard();
+        
+        if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, fade, fade_value); 
+        
+        /* add oscillating waves */
+//         if ((ADD_OSCILLATING_SOURCE)&&(i%OSCILLATING_SOURCE_PERIOD == OSCILLATING_SOURCE_PERIOD - 1))
+        if ((ADD_OSCILLATING_SOURCE)&&(i%OSCILLATING_SOURCE_PERIOD == 1))
+        {
+            if (ALTERNATE_OSCILLATING_SOURCE) sign = -sign;
+            add_circular_wave_mod(sign, -0.5, 0.0, phi, psi, xy_in);
+        }
+        if (PRINT_SPEED) print_speed_3d(speed, 0, 1.0);
+
+	if (!((NO_EXTRA_BUFFER_SWAP)&&(MOVIE))) glutSwapBuffers();
+
+	if (MOVIE)
+        {
+            if (i >= INITIAL_TIME) save_frame();
+//             if (i >= INITIAL_TIME) save_frame_counter(i);
+//             if (i >= INITIAL_TIME) save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter);
+            else printf("Initial phase time %i of %i\n", i, INITIAL_TIME);
+            
+            if ((i >= INITIAL_TIME)&&(DOUBLE_MOVIE))
+            {
+                draw_wave_sphere(1, phi, psi, xy_in, wave, wsphere, 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, CIRC_COLORBAR_B, 0, 1.0);  
+                if (PRINT_SPEED) print_speed_3d(speed, 0, 1.0);
+                glutSwapBuffers();
+                save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter);
+//                 save_frame_counter(i);
+                counter++;
+            }
+            else if (NO_EXTRA_BUFFER_SWAP) glutSwapBuffers();
+
+            /* 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_wave/");
+            }
+        }
+        else printf("Computing frame %i\n", i);
+
+    }
+
+    if (MOVIE) 
+    {
+        if (DOUBLE_MOVIE) 
+        {
+            draw_wave_sphere(0, phi, psi, xy_in, wave, wsphere, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1);
+            if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, 0, 1.0);   
+            if (PRINT_SPEED) print_speed_3d(speed, 0, 1.0);
+            glutSwapBuffers();
+            
+            if (!FADE) for (i=0; i<MID_FRAMES; i++) save_frame();
+            else for (i=0; i<MID_FRAMES; i++) 
+            {
+                if ((ROTATE_VIEW)&&(ROTATE_VIEW_WHILE_FADE))
+                {
+                    viewpoint_schedule(NSTEPS - INITIAL_TIME + i);
+                    reset_view = 1;
+                }
+                fade_value = 1.0 - (double)i/(double)MID_FRAMES;
+                draw_wave_sphere(0, phi, psi, xy_in, wave, wsphere, ZPLOT, CPLOT, COLOR_PALETTE, 1, fade_value, 1);
+                if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, 1, fade_value);   
+                if (PRINT_SPEED) print_speed_3d(speed, 1, fade_value);
+                if (!NO_EXTRA_BUFFER_SWAP) glutSwapBuffers();
+                save_frame_counter(NSTEPS + i + 1);
+            }
+            
+            if ((ROTATE_VIEW)&&(ROTATE_VIEW_WHILE_FADE)) 
+            {
+                viewpoint_schedule(NSTEPS - INITIAL_TIME);
+                reset_view = 1;
+            }
+            draw_wave_sphere(1, phi, psi, xy_in, wave, wsphere, 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, CIRC_COLORBAR_B, 0, 1.0); 
+            if (PRINT_SPEED) print_speed_3d(speed, 0, 1.0);
+            glutSwapBuffers();
+            
+            if (!FADE) for (i=0; i<END_FRAMES; i++) save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter + i);
+            else for (i=0; i<END_FRAMES; i++) 
+            {
+                if ((ROTATE_VIEW)&&(ROTATE_VIEW_WHILE_FADE)) 
+                {
+                    viewpoint_schedule(NSTEPS - INITIAL_TIME + i);
+                    reset_view = 1;
+                }
+                fade_value = 1.0 - (double)i/(double)END_FRAMES;
+                draw_wave_sphere(1, phi, psi, xy_in, wave, wsphere, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 1, fade_value, 1);
+                if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, CIRC_COLORBAR_B, 1, fade_value);   
+                if (PRINT_SPEED) print_speed_3d(speed, 1, fade_value);
+                glutSwapBuffers();
+                save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter + i);
+            }
+        }
+        else
+        {
+            if (!FADE) for (i=0; i<END_FRAMES; i++) save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter + i);
+            else for (i=0; i<END_FRAMES; i++) 
+            {
+                if ((ROTATE_VIEW)&&(ROTATE_VIEW_WHILE_FADE))
+                {
+                    viewpoint_schedule(NSTEPS - INITIAL_TIME + i);
+                    reset_view = 1;
+                }
+                fade_value = 1.0 - (double)i/(double)END_FRAMES;
+                draw_wave_sphere(0, phi, psi, xy_in, wave, wsphere, ZPLOT, CPLOT, COLOR_PALETTE, 1, fade_value, 1);
+                if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, 1, fade_value); 
+                if (PRINT_SPEED) print_speed_3d(speed, 1, fade_value);
+                glutSwapBuffers();
+                save_frame_counter(NSTEPS + 1 + counter + i);
+            }
+        }
+        
+        s = system("mv wave*.tif tif_wave/");
+    }
+    
+    free(xy_in);
+    free(phi);
+    free(psi);
+    free(tmp);
+//     free(total_energy);
+    free(color_scale);
+    free(tc);
+    free(tcc);
+    free(tgamma);
+    
+    free(wave);
+    free(laplace); 
+    free(laplace1); 
+    free(laplace2); 
+
+    
+    if (SAVE_TIME_SERIES)
+    {
+        fclose(time_series_left);
+        fclose(time_series_right);
+    }
+}
+
+
+void display(void)
+{
+    time_t rawtime;
+    struct tm * timeinfo;
+
+    time(&rawtime);
+    timeinfo = localtime(&rawtime);
+    
+    glPushMatrix();
+
+    blank();
+    glutSwapBuffers();
+    blank();
+    glutSwapBuffers();
+
+    animation();
+    sleep(SLEEP2);
+
+    glPopMatrix();
+
+    glutDestroyWindow(glutGetWindow());
+    
+    printf("Start local time and date: %s", asctime(timeinfo));
+    time(&rawtime);
+    timeinfo = localtime(&rawtime);
+    printf("Current local time and date: %s", asctime(timeinfo));
+}
+
+
+int main(int argc, char** argv)
+{
+    glutInit(&argc, argv);
+    glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH);
+    glutInitWindowSize(WINWIDTH,WINHEIGHT);
+    glutCreateWindow("Wave equation in a planar domain");
+
+    if (PLOT_2D) init_sphere_2D(); 
+    else init_3d();
+
+    glutDisplayFunc(display);
+
+    glutMainLoop();
+
+    return 0;
+}
+