Add files via upload

billiard_phase_space.c

@@ -30,7 +30,7 @@
 #include <tiffio.h>     /* Sam Leffler's libtiff library. */
 #include <time.h>
 
-#define MOVIE 1         /* set to 1 to generate movie */
+#define MOVIE 0         /* set to 1 to generate movie */
 #define SAVE_MEMORY 1       /* set to 1 to save memory when writing tiff images */
 
 #define WINWIDTH 	1280  /* window width */
@@ -52,7 +52,7 @@
 
 /* Choice of the billiard table, see global_particles.c */
 
-#define B_DOMAIN 11      /* choice of domain shape */
+#define B_DOMAIN 1      /* choice of domain shape */
 
 #define CIRCLE_PATTERN 6   /* pattern of circles */
 #define POLYLINE_PATTERN 4  /* pattern of polyline */
@@ -67,8 +67,8 @@
 #define NGOLDENSPIRAL 2000  /* max number of points for C_GOLDEN_SPIRAL arrandement */
 #define SDEPTH 2            /* Sierpinski gastket depth */
 
-#define LAMBDAMIN 6.0	/* parameter controlling shape of domain */
+#define LAMBDAMIN 0.0	/* parameter controlling shape of domain */
-#define LAMBDA 20.0	/* parameter controlling shape of domain */
+#define LAMBDA 1.5	/* parameter controlling shape of domain */
 #define MU 1.0          /* second parameter controlling shape of billiard */
 #define FOCI 1          /* set to 1 to draw focal points of ellipse */
 #define NPOLY 6             /* number of sides of polygon */
@@ -96,7 +96,7 @@
 #define PRINT_TRAJECTORY_LENGTH 0   /* set to 1 to print length of trajectory 0 */
 #define PRINT_TRAJECTORY_PERIOD 0   /* set to 1 to print period of trajectory 0 */
 #define DRAW_LENGTHS_PLOT 0         /* set to 1 to plot trajectory lengths */
-#define LENGTHS_LOG_SCALE 1         /* set to 1 to use log scale for plot of lengths */
+#define LENGTHS_LOG_SCALE 0         /* set to 1 to use log scale for plot of lengths */
 #define MAX_ANGLE 90.0         /* range of angles of trajectory */
 
 #define NSTEPS 4000      /* number of frames of movie */
@@ -141,7 +141,7 @@
 #define PLOT_LYAPUNOV 1     /* set to 1 to add plot of Lyapunov exponents */
 #define LOGSCALEX_LYAP 0    /* set to 1 to use log scale on parameter axis of Lyapunov exponents */
 #define LYAP_MAX 1.0        /* maximal Lyapunov exponent */
-#define ADAPT_TO_SYMMETRY 1 /* set to 1 to show only one symmetric part of phase space */
+#define ADAPT_TO_SYMMETRY 0 /* set to 1 to show only one symmetric part of phase space */
 #define SYMMETRY_FACTOR 3   /* proportion of phase space to be shown */
 
 #define PAUSE 1000       /* number of frames after which to pause */

global_3d.c

@@ -65,6 +65,14 @@
 
 #define VP_HORIZONTAL 0     /* rotate in a horizontal plane */
 #define VP_ORBIT 1          /* rotate in a plane containing the origin */
+#define VP_ORBIT2 11        /* rotate in a plane specified by max latitude */
+#define VP_POLAR 2          /* polar orbit */
+
+/* Type of digital elevation model */
+
+#define DEM_EARTH 0         /* DEM of Earth */
+#define DEM_MARS 1          /* DEM of Mars */
+#define DEM_MOON 2          /* DEM of the Moon */
 
 /* macros to avoid unnecessary computations in 3D plots */
 
@@ -85,6 +93,9 @@
 
 #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 PLANET ((B_DOMAIN == D_SPHERE_EARTH)||(B_DOMAIN == D_SPHERE_MARS)||(B_DOMAIN == D_SPHERE_MOON))
+#define OTHER_PLANET ((B_DOMAIN == D_SPHERE_MARS)||(B_DOMAIN == D_SPHERE_MOON))
+
 #define NMAXCIRC_SPHERE 100     /* max number of circles on sphere */
 
 int global_time = 0;
@@ -154,6 +165,9 @@ typedef struct
     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 */
+    double altitude;            /* altitude in case of Earth with digital elevation model */
+    double cos_angle;           /* cosine of light angle */
+    double cos_angle_sphere;    /* cosing of light angle for perfect sphere */
 } t_wave_sphere;
 
 

global_ljones.c

@@ -75,8 +75,10 @@
 #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_MIXER 23              /* exterior of a blender made of rectangles */
 #define S_AIRFOIL 24            /* exterior of an air foil */
+#define S_COANDA 25             /* wall for Coanda effect */
+#define S_COANDA_SHORT 26       /* shorter wall for Coanda effect */
 
 /* particle interaction */
 
@@ -92,6 +94,7 @@
 #define I_VICSEK_REPULSIVE 9  /* Vicsek-type interaction with harmonic repulsion */
 #define I_VICSEK_SPEED 10   /* Vicsek-type interaction with speed adjustment */
 #define I_VICSEK_SHARK 11   /* Vicsek-type interaction with speed adjustment, and one shark */
+#define I_COULOMB_LJ 12     /* Coulomb force regularised by Lennard-Jones repulsion */
 
 /* Boundary conditions */
 
@@ -178,6 +181,17 @@
 #define IC_TWOROCKETS 8    /* type 1 or 2 depending on rocket position */
 #define IC_TWOROCKETS_TWOFUELS 9    /* type 1 and 2 or 1 and 3 depending on rocket */
 
+/* Initial conditions for option TWO_TYPES */
+
+#define TTC_RANDOM 0        /* assign types randomly */
+#define TTC_CHESSBOARD 1    /* assign types according to chessboard, works with hex initial config */
+#define TTC_COANDA 2        /* type 1 in a band of width LAMBDA */
+
+/* Initial speed distribution */
+
+#define VI_RANDOM 0         /* random (Gaussian) initial speed distribution */
+#define VI_COANDA 1         /* nonzero speed in a band of width LAMBDA */
+
 /* Plot types */
 
 #define P_KINETIC 0       /* colors represent kinetic energy of particles */
@@ -197,6 +211,11 @@
 #define P_DIRECT_EMEAN 14 /* averaged version of P_DIRECT_ENERGY */
 #define P_NOPARTICLE 15   /* particles are not drawn (only the links between them) */
 
+/* Rotation schedules */
+
+#define ROT_SPEEDUP_SLOWDOWN 0  /* rotation speeds up and then slows down to zero */
+#define ROT_BACK_FORTH 1        /* rotation goes in one direction and then back */
+
 /* Initial position dependence types */
 
 #define IP_X 0      /* color depends on x coordinate of initial position */
@@ -241,6 +260,7 @@ typedef struct
     double vy;                  /* y velocity of particle */
     double omega;               /* angular velocity of particle */
     double mass_inv;            /* inverse of particle mass */
+    double charge;              /* electric charge */
     double inertia_moment_inv;  /* inverse of moment of inertia */
     double fx;                  /* x component of force on particle */
     double fy;                  /* y component of force on particle */
@@ -313,6 +333,8 @@ typedef struct
     double angle01, angle02;    /* initial values of angles in which concave corners repel */
     double fx, fy;              /* x and y-components of force on segment */
     double torque;              /* torque on segment with respect to its center */
+    double pressure;            /* pressure acting on segement */
+    double avrg_pressure;       /* time-averaged pressure */
     short int inactivate;       /* set to 1 for segment to become inactive at time SEGMENT_DEACTIVATION_TIME */
 } t_segment;
 
@@ -361,6 +383,7 @@ typedef struct
     double angle;               /* orientation of obstacle */
     double omega;               /* angular speed of obstacle */
     double bdry_fx, bdry_fy;    /* components of boundary force */
+    double efield, bfield;      /* electric and magnetic field */
 } t_lj_parameters;
 
 int ncircles, nobstacles, nsegments, ngroups = 1, counter = 0;

global_pdes.c

@@ -77,6 +77,7 @@
 #define D_ONE_FUNNEL 581        /* one funnel */
 #define D_LENSES_RING 59        /* several lenses forming a ring */
 #define D_MAZE_CIRCULAR 60      /* circular maze */
+#define D_LENS 61               /* symmetric lens made of circular faces */
 
 #define D_WING 70               /* complement of wing-shaped domain */
 #define D_TESLA 71              /* Tesla valve */
@@ -89,6 +90,9 @@
 #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 D_SPHERE_JULIA_CUBIC 85 /* Julia set for cubic polynomial on Riemann sphere */
+#define D_SPHERE_MARS 86        /* map of Mars */
+#define D_SPHERE_MOON 87        /* map of the Moon */
 
 #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) */
@@ -149,6 +153,8 @@
 #define IOR_POISSON_WELLS 8     /* wells located on a random Poisson disc process */
 #define IOR_PPP_WELLS 9         /* wells located on a Poisson point process */
 
+#define IOR_EARTH_DEM 10        /* digital elevation model (for waves on sphere) */
+
 /* Boundary conditions */
 
 #define BC_DIRICHLET 0   /* Dirichlet boundary conditions */
@@ -167,6 +173,18 @@
 #define OSC_WAVE_PACKET 2   /* Gaussian wave packet */
 #define OSC_CHIRP 3     /* chirp (linearly accelerating frequency) */
 
+/* Wave packet types */
+
+#define WP_RANDOM1 0    /* random, variant 1 */
+#define WP_RANDOM2 1    /* random, variant 2 */
+#define WP_PAIR 2       /* 2 sources */
+
+/* Wave packet envelope types */
+
+#define WE_SINE 0       /* sine envelope */
+#define WE_CUTOFF 1     /* sharp cut-off */
+#define WE_CONSTANT 2   /* constant envelope */
+
 /* For debugging purposes only */
 // #define FLOOR 0         /* set to 1 to limit wave amplitude to VMAX */
 // #define VMAX 10.0       /* max value of wave amplitude */

heat.c

@@ -212,6 +212,7 @@
 #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 OSCILLATING_SOURCE_PERIOD 20    /* period of oscillating source */
 /* end of constants only used by sub_wave and sub_maze */
 
 #include "global_pdes.c"

lennardjones.c

@@ -22,7 +22,7 @@
 /*  It may be possible to increase parameter PAUSE                               */
 /*                                                                               */
 /*  create movie using                                                           */
-/*  ffmpeg -i lj.%05d.tif -vcodec libx264 lj.mp4                             */
+/*  ffmpeg -i lj.%05d.tif -vcodec libx264 lj.mp4                                 */
 /*                                                                               */
 /*********************************************************************************/
 
@@ -50,8 +50,8 @@
 
 /* General geometrical parameters */
 
-#define WINWIDTH 	1280  /* window width */
+#define WINWIDTH 	1600  /* window width */
-#define WINHEIGHT 	720   /* window height */
+#define WINHEIGHT 	900   /* window height */
 
 #define XMIN -2.0
 #define XMAX 2.0	/* x interval */
@@ -59,7 +59,7 @@
 #define YMAX 1.125	/* y interval for 9/16 aspect ratio */
 
 #define INITXMIN -2.0
-#define INITXMAX 2.0	/* x interval for initial condition */
+#define INITXMAX 2.03	/* x interval for initial condition */
 #define INITYMIN -1.1
 #define INITYMAX 1.1	/* y interval for initial condition */
 
@@ -86,14 +86,15 @@
 #define OBSTACLE_PATTERN 4  /* 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 24      /* pattern of repelling segments, see list in global_ljones.c */
+#define SEGMENT_PATTERN 25      /* 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 TWO_TYPES 1         /* set to 1 to have two types of particles */
 #define TYPE_PROPORTION 0.5 /* proportion of particles of first type */
+#define TWOTYPE_CONFIG 2    /* choice of types, see TTC_ list in global_ljones.c */
 #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 */
@@ -110,9 +111,9 @@
 #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.75	    /* parameter controlling the dimensions of domain */
+#define LAMBDA 0.3	    /* parameter controlling the dimensions of domain */
-#define MU 0.012 	    /* parameter controlling radius of particles */
+#define MU 0.014 	    /* parameter controlling radius of particles */
-#define MU_B 0.010          /* parameter controlling radius of particles of second type */
+#define MU_B 0.01           /* 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 */ 
@@ -121,10 +122,8 @@
 #define MANDELLEVEL 1000    /* iteration level for Mandelbrot set */
 #define MANDELLIMIT 10.0    /* limit value for approximation of Mandelbrot set */
 #define FOCI 1              /* set to 1 to draw focal points of ellipse */
-// #define NGRIDX 140           /* number of grid point for grid of disks */
+#define NGRIDX 101          /* number of grid point for grid of disks */
-// #define NGRIDY 70           /* number of grid point for grid of disks */
+#define NGRIDY 47           /* 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 */
@@ -137,9 +136,8 @@
 
 /* Parameters for length and speed of simulation */
 
-#define NSTEPS 4525      /* number of frames of movie */
+#define NSTEPS 3200      /* number of frames of movie */
-// #define NSTEPS 1000      /* number of frames of movie */
+#define NVID 15         /* number of iterations between images displayed on screen */
-#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 */
@@ -152,7 +150,7 @@
 #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 END_FRAMES 250   /* number of still frames at end of movie */
 
 /* Boundary conditions, see list in global_ljones.c */
 
@@ -160,18 +158,18 @@
 
 /* Plot type, see list in global_ljones.c  */
 
-#define PLOT 13
+#define PLOT 5
 #define PLOT_B 11        /* plot type for second movie */
 
 #define DRAW_BONDS 1    /* set to 1 to draw bonds between neighbours */
 #define COLOR_BONDS 1   /* set to 1 to color bonds according to length */
 #define FILL_TRIANGLES 0    /* set to 1 to fill triangles between neighbours */
 #define ALTITUDE_LINES 0    /* set to 1 to add horizontal lines to show altitude */
-#define COLOR_SEG_GROUPS 1  /* set to 1 to collor segment groups differently */
+#define COLOR_SEG_GROUPS 0  /* 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 ERATIO 0.995          /* ratio for time-averaging in P_EMEAN color scheme */
-#define DRATIO 0.995          /* ratio for time-averagin in P_DIRECT_EMEAN color scheme */
+#define DRATIO 0.995          /* ratio for time-averaging in P_DIRECT_EMEAN color scheme */
 
 /* Color schemes */
 
@@ -195,15 +193,15 @@
 #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 COLOR_HUESHIFT 0.0     /* shift in color hue (for some cyclic palettes) */
 
-#define PRINT_PARAMETERS 1  /* set to 1 to print certain parameters */
+#define PRINT_PARAMETERS 0  /* set to 1 to print certain parameters */
 #define PRINT_TEMPERATURE 0 /* set to 1 to print current temperature */
-#define PRINT_ANGLE 1               /* set to 1 to print obstacle orientation */
+#define PRINT_ANGLE 0               /* set to 1 to print obstacle orientation */
-#define PRINT_OMEGA 1               /* set to 1 to print angular 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 0     /* set to 1 to print velocity of moving segments */
-#define PRINT_SEGMENTS_FORCE 1      /* set to 1 to print force on segments */
+#define PRINT_SEGMENTS_FORCE 0      /* set to 1 to print force on segments */
 #define FORCE_FACTOR 0.1            /* factor controlling length of force vector */
 
 /* particle properties */
@@ -212,7 +210,7 @@
 #define ENERGY_HUE_MAX 50.0         /* color of saturated particle */
 #define PARTICLE_HUE_MIN 359.0      /* color of original particle */
 #define PARTICLE_HUE_MAX 0.0        /* color of saturated particle */
-#define PARTICLE_EMAX 5000.0        /* energy of particle with hottest color */
+#define PARTICLE_EMAX 1000.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 */
@@ -221,35 +219,35 @@
 #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 7.0   /* Lennard-Jones equilibrium distance */
-#define EQUILIBRIUM_DIST_B 5.0  /* Lennard-Jones equilibrium distance for second type of particle */
+#define EQUILIBRIUM_DIST_B 6.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 100.0          /* damping coefficient of particles */
+#define DAMPING 0.02          /* 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 1.0    /* mass of particle of radius MU */
-#define PARTICLE_MASS_B 1.0   /* mass of particle of radius MU */
+#define PARTICLE_MASS_B 6.5   /* 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 V_INITIAL 1000.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 V_INITIAL_TYPE 1    /* type of initial speed distribution (see VI_ in global_ljones.c) */
+
 #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 BETA 0.0002          /* 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.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 0.0             /* gravity acting on all particles */
+#define GRAVITY_X 0.0          /* horizontal gravity acting on all particles */
-#define GRAVITY_X 5000.0          /* horizontal gravity acting on all particles */
+#define CIRCULAR_GRAVITY 0     /* set to 1 to have gravity directed to center */
 #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 */
@@ -258,6 +256,15 @@
 #define KSPRING_VICSEK 0.2   /* spring constant for I_VICSEK_SPEED interaction */
 #define VICSEK_REPULSION 10.0    /* repulsion between particles in Vicsek model */
 
+#define ADD_EFIELD 1     /* set to 1 to add an electric field */
+#define EFIELD 10000.0  /* value of electric field */
+#define ADD_BFIELD 0     /* set to 1 to add a magnetic field */
+#define BFIELD 1000.0    /* value of magnetic field */
+#define CHARGE 0.0       /* charge of particles of first type */
+#define CHARGE_B 1.0     /* charge of particles of second type */
+#define INCREASE_E 0     /* set to 1 to increase electric field */
+#define EFIELD_FACTOR 1000.0    /* factor by which to increase electric field */
+
 #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 */  
@@ -267,7 +274,7 @@
 #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 DRAW_CROSS 0          /* 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 */ 
@@ -290,7 +297,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.25  /* radius of obstacle for circle boundary conditions */
+#define OBSTACLE_RADIUS 0.15  /* radius of obstacle for circle boundary conditions */
 #define FUNNEL_WIDTH  0.25   /* funnel width for funnel boundary conditions */
 #define OBSTACLE_XMIN 0.0    /* initial position of obstacle */
 #define OBSTACLE_XMAX 3.0    /* final position of obstacle */
@@ -325,24 +332,28 @@
 #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 ROTATE_BOUNDARY 0           /* 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 ROTATION_SCHEDULE 0         /* time-dependence of rotation angle, see ROT_* in global_ljones.c */
 #define PERIOD_ROTATE_BOUNDARY 1000  /* period of rotating boundary */
 #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 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 */
 #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_X0 0.0        /* initial position of segments */
-#define SEGMENTS_Y0 0.0        /* initial position of segments */
+#define SEGMENTS_Y0 -0.6      /* 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 SHOW_SEGMENTS_PRESSURE 0    /* set to 1 to show (averaged) pressure acting on segments */
+#define SEGMENT_PMAX 7.5e7        /* pressure of segment with hottest color */
+#define P_AVRG_FACTOR 0.02      /* factor in computation of mean pressure */
 
 #define MOVE_SEGMENT_GROUPS 0       /* set to 1 to group segments into moving units */
 #define SEGMENT_GROUP_MASS 500.0   /* mass of segment group */
@@ -407,8 +418,8 @@
 #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 HASHX 60    /* size of hashgrid in x direction */
-#define HASHY 60    /* size of hashgrid in y direction */
+#define HASHY 30    /* 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 */
 
@@ -473,6 +484,23 @@ double repel_schedule(int i)
 }
 
 
+double efield_schedule(int i)
+{
+    static double efactor;
+    static int first = 1;
+    double efield;
+    
+    if (first) 
+    {
+        efactor = EFIELD_FACTOR/(double)(INITIAL_TIME + NSTEPS);
+        first = 0;
+    }
+    efield = EFIELD*(double)i*efactor;
+    printf("E = %.3lg\n", efield);
+    return(efield);
+}
+
+
 double temperature_schedule(int i)
 {
     static double bexponent, omega, bexp2;
@@ -614,20 +642,26 @@ double rotation_angle(double phase)
 //     return(0.5*OMEGAMAX*(phase - sin(phase)));
     
     /* case of centrifuge remaining at constant speed for a while */
-    if (phase < ROTATE_CHANGE_TIME)
+    switch (ROTATION_SCHEDULE) {
-    {
+        case (ROT_SPEEDUP_SLOWDOWN):
-//         angular_speed = 0.5*OMEGAMAX*(1.0 - cos(phase*PI/ROTATE_CHANGE_TIME));
+        {
-        return(0.5*OMEGAMAX*(phase - (ROTATE_CHANGE_TIME/PI)*sin(phase*PI/ROTATE_CHANGE_TIME)));
+            if (phase < ROTATE_CHANGE_TIME)
-    }
+            {
-    else if (phase < 1.0 - ROTATE_CHANGE_TIME)
+                return(0.5*OMEGAMAX*(phase - (ROTATE_CHANGE_TIME/PI)*sin(phase*PI/ROTATE_CHANGE_TIME)));
-    {
+            }
-//         angular_speed = OMEGAMAX;
+            else if (phase < 1.0 - ROTATE_CHANGE_TIME)
-        return(0.5*OMEGAMAX*(2.0*phase - ROTATE_CHANGE_TIME));
+            {
-    }
+                return(0.5*OMEGAMAX*(2.0*phase - ROTATE_CHANGE_TIME));
-    else
+            }
-    {
+            else
-//         angular_speed = 0.5*OMEGAMAX*(1.0 + cos((phase - 1.0 + ROTATE_CHANGE_TIME)*PI/ROTATE_CHANGE_TIME));
+            {
-        return(0.5*OMEGAMAX*(2.0 - 2.0*ROTATE_CHANGE_TIME + phase - 1.0 + (ROTATE_CHANGE_TIME/PI)*sin((1.0-phase)*PI/ROTATE_CHANGE_TIME)));
+                return(0.5*OMEGAMAX*(2.0 - 2.0*ROTATE_CHANGE_TIME + phase - 1.0 + (ROTATE_CHANGE_TIME/PI)*sin((1.0-phase)*PI/ROTATE_CHANGE_TIME)));
+            }
+        }
+        case (ROT_BACK_FORTH):
+        {
+            return(OMEGAMAX*(1.0 - cos(DPI*phase))/DPI);
+        }
     }
 }
 
@@ -1252,7 +1286,8 @@ void animation()
         printf("Computing frame %d\n",i);
         
         if (INCREASE_KREPEL) params.krepel = repel_schedule(i);
-        if (INCREASE_BETA) params.beta = temperature_schedule(i);        
+        if (INCREASE_BETA) params.beta = temperature_schedule(i);  
+        if (INCREASE_E) params.efield = efield_schedule(i);
         if (DECREASE_CONTAINER_SIZE) 
         {
             params.xmincontainer = container_size_schedule(i);
@@ -1358,12 +1393,33 @@ void animation()
 
                 /* add gravity */
                 if (INCREASE_GRAVITY) particle[j].fy -= params.gravity/particle[j].mass_inv;
+                else if (CIRCULAR_GRAVITY)
+                {
+                    particle[j].fx -= GRAVITY*particle[j].xc/particle[j].mass_inv;
+                    particle[j].fy -= GRAVITY*particle[j].yc/particle[j].mass_inv;
+                }
                 else 
                 {
                     particle[j].fy -= GRAVITY/particle[j].mass_inv;
                     particle[j].fx += GRAVITY_X/particle[j].mass_inv;
                 }
                 
+                /* add electric force */
+                if (ADD_EFIELD) 
+                {
+                    if (INCREASE_E) particle[j].fx += params.efield*particle[j].charge;
+                    else particle[j].fx += EFIELD*particle[j].charge;
+                }
+                
+                /* add magnetic force */
+                if (ADD_BFIELD)
+                {
+//                     particle[j].fx += BFIELD*particle[j].charge*particle[j].vy;
+//                     particle[j].fy -= BFIELD*particle[j].charge*particle[j].vx;
+                    particle[j].fx += BFIELD*particle[j].charge*particle[j].vy*particle[j].mass_inv;
+                    particle[j].fy -= BFIELD*particle[j].charge*particle[j].vx*particle[j].mass_inv;
+                }
+                
                 if (FLOOR_FORCE)
                 {
                     if (particle[j].fx > FMAX) particle[j].fx = FMAX;

makefile

@@ -1,4 +1,4 @@
-LIBS = -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lXmu -lglut -fopenmp
+LIBS = -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lglut -fopenmp
 
 .c:
 	gcc -o $@ $< -O3 $(LIBS)

mangrove.c

@@ -251,16 +251,11 @@
 #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 */
+#define OSCILLATING_SOURCE_PERIOD 20    /* period of oscillating source */
 /* end of constants only used by sub_wave and sub_maze */
 
 #include "global_pdes.c"

schrodinger.c

@@ -192,6 +192,7 @@
 #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 OSCILLATING_SOURCE_PERIOD 20    /* period of oscillating source */
 /* end of constants only used by sub_wave and sub_maze */
 
 

sub_lj.c

@@ -34,8 +34,8 @@ int writetiff(char *filename, char *description, int x, int y, int width, int he
   glPixelStorei(GL_PACK_ALIGNMENT, 1);
 
   glReadPixels(x, y, width, height, GL_RGB, GL_UNSIGNED_BYTE, image);
-  TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32) width);
+  TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32_t) width);
-  TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32) height);
+  TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32_t) height);
   TIFFSetField(file, TIFFTAG_BITSPERSAMPLE, 8);
   TIFFSetField(file, TIFFTAG_COMPRESSION, compression);
   TIFFSetField(file, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB);
@@ -622,15 +622,18 @@ void init_particle_config(t_particle particles[NMAXCIRCLES])
         case (C_SQUARE):
         {
             ncircles = NGRIDX*NGRIDY;
-            dy = (YMAX - YMIN)/((double)NGRIDY);
+            dx = (INITXMAX - INITXMIN)/((double)NGRIDX);
+            dy = (INITYMAX - INITYMIN)/((double)NGRIDY);
             for (i = 0; i < NGRIDX; i++)
                 for (j = 0; j < NGRIDY; j++)
                 {
                     n = NGRIDY*i + j;
-                    particles[n].xc = ((double)(i-NGRIDX/2) + 0.5)*dy;
+                    particles[n].xc = INITXMIN + ((double)i - 0.5)*dx;   
-                    particles[n].yc = YMIN + ((double)j + 0.5)*dy;
+                    particles[n].yc = INITYMIN + ((double)j - 0.5)*dy;
                     particles[n].radius = MU;
-                    particles[n].active = 1;
+                    /* activate only circles that intersect the domain */
+                    if ((particles[n].yc < INITYMAX + MU)&&(particles[n].yc > INITYMIN - MU)&&(particles[n].xc < INITXMAX + MU)&&(particles[n].xc > INITXMIN - MU)) particles[n].active = 1;
+                    else particles[n].active = 0;
                 }
             break;
         }
@@ -1616,7 +1619,7 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
 /* initialise linear obstacle configuration */
 {
     int i, j, cycle = 0, iminus, iplus, nsides, n, concave = 1; 
-    double angle, angle2, dangle, dx, width, height, a, b, length, xmid = 0.5*(BCXMIN + BCXMAX), lpocket, r, x, x1, y1, x2, y2, nozx, nozy, y, dy, ca, sa;
+    double angle, angle2, dangle, dx, width, height, a, b, length, xmid = 0.5*(BCXMIN + BCXMAX), lpocket, r, x, x1, y1, x2, y2, nozx, nozy, y, yy, dy, ca, sa, padding;
     
     switch (SEGMENT_PATTERN) {
         case (S_RECTANGLE):
@@ -2399,6 +2402,82 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
             
             break;
         }
+        case (S_COANDA):
+        {
+            y1 = SEGMENTS_Y0;
+            width = 0.05;
+            padding = 0.01;
+            
+            height = 0.25;
+            dx = (XMAX - XMIN + 2.0*padding)/(double)(NPOLY-1);
+            
+            for (i=0; i<NPOLY; i++)
+            {
+                x = XMIN - padding + (double)i*dx;
+                y = y1 + height*cos(PI*x/XMAX);
+                
+                segment[i].x1 = x;
+                segment[i].y1 = y + width;
+                segment[i].x2 = x + dx;
+                segment[i].x2 = y1 + height*cos(PI*(x+dx)/XMAX) + width;
+                
+                segment[i].concave = 1;
+            }
+            
+            for (i=0; i<NPOLY; i++)
+            {
+                x = XMAX + padding - (double)i*dx;
+                y = y1 + height*cos(PI*x/XMAX);
+                
+                segment[NPOLY+i].x1 = x;
+                segment[NPOLY+i].y1 = y - width;
+                segment[NPOLY+i].x2 = x - dx;
+                segment[NPOLY+i].x2 = y1 + height*cos(PI*(x-dx)/XMAX) - width;
+                
+                segment[NPOLY+i].concave = 1;
+            }
+            
+            nsegments = 2*NPOLY;
+            
+            cycle = 1;
+            concave = 1;    
+            
+            for (i=0; i<nsegments; i++) 
+            {
+                segment[i].inactivate = 0;
+            }
+            
+            break;
+        }
+        case (S_COANDA_SHORT):
+        {
+            y1 = SEGMENTS_Y0;
+            width = 0.05;
+            padding = 0.1;
+            x1 = XMIN + padding;
+            
+            height = 0.2;
+            length = 0.85*(XMAX - XMIN - 2.0*padding);
+            dx = length/(double)(NPOLY-1);
+            
+            for (i=0; i<NPOLY; i++)
+            {
+                x = x1 + (double)i*dx;
+                y = y1 - height*cos(DPI*(x-x1)/length);
+                yy = y1 - height*cos(DPI*(x+dx-x1)/length);
+                add_rotated_angle_to_segments(x, y, x+dx, yy, width, 1, segment, 0);
+            }
+            
+            cycle = 0;
+            concave = 1;    
+            
+            for (i=0; i<nsegments; i++) 
+            {
+                segment[i].inactivate = 0;
+            }
+            
+            break;
+        }
         default: 
         {
             printf("Function init_segment_config not defined for this pattern \n");
@@ -2491,6 +2570,10 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
         segment[i].angle02 = segment[i].angle2;
     }
     
+    /* initialize averaged pressure */
+    if (SHOW_SEGMENTS_PRESSURE) for (i=0; i<nsegments; i++) 
+        segment[i].avrg_pressure = 0.0;
+    
 //     for (i=0; i<nsegments; i++) 
 //     {
 //         printf("Segment %i: (x1, y1) = (%.3lg,%.3lg), (x2, y2) = (%.3lg,%.3lg)\n (nx, ny) = (%.3lg,%.3lg), c = %.3lg, length = %.3lg\n", i, segment[i].x1, segment[i].y1, segment[i].x2, segment[i].y2, segment[i].nx, segment[i].ny, segment[i].c, segment[i].length);
@@ -2550,7 +2633,7 @@ int in_segment_region(double x, double y, t_segment segment[NMAXSEGMENTS])
 /* returns 1 if (x,y) is inside region delimited by obstacle segments */
 {
     int i;
-    double angle, dx, height, width, theta, lx, ly, x1, y1, x2, y2, padding, ca, sa, r;
+    double angle, dx, height, width, length, theta, lx, ly, x1, y1, x2, y2, padding, ca, sa, r;
     
     if (x >= BCXMAX) return(0);
     if (x <= BCXMIN) return(0);
@@ -2750,6 +2833,17 @@ int in_segment_region(double x, double y, t_segment segment[NMAXSEGMENTS])
             y1 += 0.5*x1*x1;
             return(x1*x1 + 25.0*y1*y1 > LAMBDA*LAMBDA*(1.0 + padding));
         }
+        case (S_COANDA):
+        {
+            return(vabs(y - SEGMENTS_Y0 - 0.25*cos(PI*x/XMAX)) > 0.1);
+        }
+        case (S_COANDA_SHORT):
+        {
+            x1 = XMIN + 0.1;
+            length = 0.85*(XMAX - XMIN - 0.2);
+            if (x > x1 + length + 0.1) return(1);
+            return(vabs(y - SEGMENTS_Y0 + 0.2*cos(DPI*(x-x1)/length)) > 0.05);
+        }
         case (S_MAZE):
         {
             for (i=0; i<nsegments; i++)
@@ -3384,7 +3478,7 @@ int compute_particle_interaction(int i, int k, double force[2], double *torque,
 /* compute repelling force and torque of particle #k on particle #i */
 /* returns 1 if distance between particles is smaller than NBH_DIST_FACTOR*MU */
 {
-    double x1, y1, x2, y2, r, f, angle, aniso, fx, fy, ff[2], dist_scaled, spin_f, ck, sk, ck_rel, sk_rel, alpha, amp;
+    double x1, y1, x2, y2, r, f, angle, aniso, fx, fy, ff[2], dist_scaled, spin_f, ck, sk, ck_rel, sk_rel, alpha, amp, charge;
     static double dxhalf = 0.5*(BCXMAX - BCXMIN), dyhalf = 0.5*(BCYMAX - BCYMIN);
     int wwrapx, wwrapy, twrapx, twrapy;
     
@@ -3514,6 +3608,17 @@ int compute_particle_interaction(int i, int k, double force[2], double *torque,
             }
             break;
         }
+        case (I_COULOMB_LJ):
+        {
+            charge = particle[i].charge*particle[k].charge;
+            f = -100.0*krepel*charge/(1.0e-12 + distance*distance);
+            if (charge < 0.0) 
+                f += 0.01*krepel*lennard_jones_force(distance, particle[k]);
+            force[0] = f*ca;
+            force[1] = f*sa;   
+            break;
+        }
+        
     }
 
     if (ROTATION) 
@@ -3998,7 +4103,7 @@ void compute_particle_colors(t_particle particle, int plot, double rgb[3], doubl
             if (hue > DPI) hue -= DPI;
             hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*(hue)/DPI;
             ej = particle.emean;
-            if (ej < 0.1*PARTICLE_EMAX) lum = 10.0*ej/PARTICLE_EMAX;
+            if (ej < 0.5*PARTICLE_EMAX) lum = 2.0*ej/PARTICLE_EMAX;
             else lum = 1.0;
             *radius = particle.radius;
             *width = BOUNDARY_WIDTH;
@@ -4127,6 +4232,21 @@ void set_segment_group_color(int group, double lum, double rgb[3])
     glColor3f(lum*rgb[0], lum*rgb[1], lum*rgb[2]);
 }
 
+void set_segment_pressure_color(double pressure, double lum, double rgb[3])
+{
+    double hue;
+    
+//     if (pressure < 0.0) pressure = 0.0;
+    hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*pressure/SEGMENT_PMAX;
+    if (hue > PARTICLE_HUE_MIN) hue = PARTICLE_HUE_MIN;
+    else if (hue < PARTICLE_HUE_MAX) hue = PARTICLE_HUE_MAX;
+        
+    hsl_to_rgb_turbo(hue, 0.9, lum, rgb);
+    
+    glColor3f(rgb[0], rgb[1], rgb[2]);
+}
+
+
 void draw_altitude_lines()
 {
     int i, i1, i2;
@@ -4625,8 +4745,18 @@ int draw_special_particle(t_particle particle, double xc1, double yc1, double ra
 void draw_one_particle(t_particle particle, double xc, double yc, double radius, double angle, int nsides, double width, double rgb[3])
 /* draw one of the particles */ 
 {
-    double ca, sa, x1, x2, y1, y2, xc1, yc1, wangle, newradius = radius;
+    double x1, x2, y1, y2, xc1, yc1, wangle, newradius = radius;
     int wsign, cont = 1, draw = 1;
+    static double ca, sa;
+    static int first = 1;
+    
+    if (first)
+    {
+        angle = APOLY*PID;
+        ca = cos(angle);
+        sa = sin(angle);
+        first = 0;
+    }
     
     if (CENTER_VIEW_ON_OBSTACLE) xc1 = xc - xshift;
     else xc1 = xc;
@@ -4651,14 +4781,16 @@ void draw_one_particle(t_particle particle, double xc, double yc, double radius,
     
     /* draw crosses on particles of second type */
     if ((TWO_TYPES)&&(DRAW_CROSS))
-        if (particle.type == 1)
+        if (particle.type == 2)
         {
-            if (ROTATION) angle = angle + APOLY*PID;
+            if (ROTATION) 
-            else angle = APOLY*PID;
+            {
-            ca = cos(angle);
+                angle = angle + APOLY*PID;
-            sa = sin(angle);
+                ca = cos(angle);
-            glLineWidth(3);
+                sa = sin(angle);
-            glColor3f(0.0, 0.0, 0.0);
+            }
+            glLineWidth(2);
+            glColor3f(1.0, 1.0, 1.0);
             x1 = xc1 - MU_B*ca;
             y1 = yc1 - MU_B*sa;
             x2 = xc1 + MU_B*ca;
@@ -5002,6 +5134,7 @@ void draw_container(double xmin, double xmax, t_obstacle obstacle[NMAXOBSTACLES]
         for (i = 0; i < nsegments; i++) if (segment[i].active)
         {
             if (COLOR_SEG_GROUPS) set_segment_group_color(segment[i].group, 1.0, rgb);
+            else if (SHOW_SEGMENTS_PRESSURE) set_segment_pressure_color(segment[i].avrg_pressure, 1.0, rgb);
             draw_line(segment[i].x1 - xtrack, segment[i].y1 - ytrack, segment[i].x2 - xtrack, segment[i].y2 - ytrack);
         }
     }
@@ -5248,6 +5381,12 @@ void print_omega(double angle, double angular_speed, double fx, double fy)
         glColor3f(1.0, 1.0, 1.0);
         sprintf(message, "Fy = %.4f", fy);
         write_text(xrighttext + 0.1, y1, message);
+        
+        y1 -= 0.1;
+        erase_area_hsl(xrightbox, y1 + 0.025, 0.42, 0.05, 0.0, 0.9, 0.0);
+        glColor3f(1.0, 1.0, 1.0);
+        sprintf(message, "Fy/Fx = %.4f", fy/fx);
+        write_text(xrighttext + 0.1, y1, message);
     }
 }
 
@@ -5389,10 +5528,17 @@ void print_parameters(t_lj_parameters params, short int left, double pressure[N_
     }   
     else if (INCREASE_KREPEL)  /* print force constant */
     {
-        erase_area_hsl(xbox, y + 0.025*scale, 0.22*scale, 0.05*scale, 0.0, 0.9, 0.0);
+        erase_area_hsl(xbox, y + 0.025*scale, 0.35*scale, 0.05*scale, 0.0, 0.9, 0.0);
         glColor3f(1.0, 1.0, 1.0);
-        sprintf(message, "Force %.0f", params.krepel);
+        sprintf(message, "Force constant %.0f", params.krepel);
-        write_text(xtext + 0.28, y, message);
+        write_text(xtext + 0.03, y, message);
+    }   
+    else if (INCREASE_E)  /* print electric field */
+    {
+        erase_area_hsl(xbox, y + 0.025*scale, 0.27*scale, 0.05*scale, 0.0, 0.9, 0.0);
+        glColor3f(1.0, 1.0, 1.0);
+        sprintf(message, "E field = %.2f", 25.0*NVID*DT_PARTICLE*params.efield);
+        write_text(xtext + 0.08, y, message);
     }   
     
     if (RECORD_PRESSURES) 
@@ -5869,9 +6015,7 @@ void print_particles_speeds(t_particle particle[NMAXCIRCLES])
     write_text(xrighttext + 0.1, y, message);
 }
 
-double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NMAXOBSTACLES], 
+double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NMAXOBSTACLES], t_segment segment[NMAXSEGMENTS], double xleft, double xright, double *pleft, double *pright, double pressure[N_PRESSURES], int wall)
-                              t_segment segment[NMAXSEGMENTS], 
-                              double xleft, double xright, double *pleft, double *pright, double pressure[N_PRESSURES], int wall)
 {
     int i, k;
     double xmin, xmax, ymin, ymax, padding, r, rp, r2, cphi, sphi, f, fperp = 0.0, x, y, xtube, distance, dx, dy, width, ybin, angle, x1, x2, h, ytop, norm, dleft, dplus, dminus, tmp_pleft = 0.0, tmp_pright = 0.0, proj, pscal, pvect, pvmin;
@@ -5914,8 +6058,6 @@ double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacl
                 particle[j].fx += f*segment[i].nx;
                 particle[j].fy += f*segment[i].ny;
                 
-                
-            
                 if ((MOVE_BOUNDARY)||(MOVE_SEGMENT_GROUPS)||(PRINT_SEGMENTS_FORCE))
                 {
                     segment[i].fx -= f*segment[i].nx;
@@ -5923,6 +6065,14 @@ double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacl
                     segment[i].torque -= (x - segment[i].xc)*f*segment[i].ny - (y - segment[i].yc)*f*segment[i].nx;
 //                     printf("Segment %i: f = (%.3lg, %.3lg)\n", i, segment[i].fx, segment[i].fy);
                 }
+                
+                if (SHOW_SEGMENTS_PRESSURE)
+                {
+                    segment[i].pressure = f/segment[i].length;
+                    segment[i].avrg_pressure *= (1.0 - P_AVRG_FACTOR);
+                    segment[i].avrg_pressure += P_AVRG_FACTOR*segment[i].pressure;
+//                      printf("Segment %i: pressure = %.3lg\n", i, segment[i].avrg_pressure);
+                }
             }
             if ((VICSEK_INT)&&(vabs(distance) < 1.5*r))
             {
@@ -6480,21 +6630,51 @@ int reorder_particles(t_particle particle[NMAXCIRCLES], double py[NMAXCIRCLES],
 }
 
 
+int twotype_config(int i, t_particle particle[NMAXCIRCLES])
+/* assign different particle types */
+{
+    switch (TWOTYPE_CONFIG) {
+        case (TTC_RANDOM): return((double)rand()/RAND_MAX > TYPE_PROPORTION);
+        case (TTC_CHESSBOARD): 
+        {
+            switch (CIRCLE_PATTERN) {
+                case (C_SQUARE): return((i%NGRIDY + i/NGRIDY)%2 == 0);
+                case (C_HEX): return(i%2 == 0);
+                default: return(i%2 == 0);
+            }
+        }
+        case (TTC_COANDA): 
+        {
+            if (vabs(particle[i].yc) < LAMBDA) return(1);
+            if (vabs(particle[i].yc) < LAMBDA + MU) return(-1);
+            return(0);
+        }
+        default: return(0);
+    }
+}
+
+
 int initialize_configuration(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HASHX*HASHY], t_obstacle obstacle[NMAXOBSTACLES], double px[NMAXCIRCLES], double py[NMAXCIRCLES], double pangle[NMAXCIRCLES], int tracer_n[N_TRACER_PARTICLES],
 t_segment segment[NMAXSEGMENTS])
 /* initialize all particles, obstacles, and the hashgrid */
 {
-    int i, j, k, n, nactive = 0;
+    int i, j, k, n, type, nactive = 0;
     double x, y, h, xx, yy, rnd;
     
     for (i=0; i < ncircles; i++) 
     {
         /* set particle type */
         particle[i].type = 0;
-        if ((TWO_TYPES)&&((double)rand()/RAND_MAX > TYPE_PROPORTION)) 
+//         if ((TWO_TYPES)&&((double)rand()/RAND_MAX > TYPE_PROPORTION)) 
+        if ((TWO_TYPES)&&(twotype_config(i, particle))) 
         {
-            particle[i].type = 2;
+            type = twotype_config(i, particle);
-            particle[i].radius = MU_B;
+            if (type == 1)
+            {
+                particle[i].type = 2;
+                particle[i].radius = MU_B;
+            }
+            else if (type == -1) particle[i].active = 0;
         }
         if ((INTERACTION == I_VICSEK_SHARK)&&(i==1)) 
         {
@@ -6524,6 +6704,7 @@ t_segment segment[NMAXSEGMENTS])
             particle[i].spin_freq = SPIN_INTER_FREQUENCY;
             particle[i].mass_inv = 1.0/PARTICLE_MASS;
             particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT;
+            particle[i].charge = CHARGE;
         }
         else 
         {
@@ -6533,10 +6714,25 @@ t_segment segment[NMAXSEGMENTS])
             particle[i].spin_freq = SPIN_INTER_FREQUENCY_B;
             particle[i].mass_inv = 1.0/PARTICLE_MASS_B;
             particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT_B;
+            particle[i].charge = CHARGE_B;
+        }
+               
+        switch (V_INITIAL_TYPE)
+        {
+            case (VI_RANDOM):
+            {
+                particle[i].vx = V_INITIAL*gaussian();
+                particle[i].vy = V_INITIAL*gaussian();
+                break;
+            }
+            case (VI_COANDA):
+            {
+                if (vabs(particle[i].yc) < LAMBDA) particle[i].vx = V_INITIAL;
+                else particle[i].vx = 0.0;
+                particle[i].vy = 0.0;
+                break;
+            }            
         }
-                
-        particle[i].vx = V_INITIAL*gaussian();
-        particle[i].vy = V_INITIAL*gaussian();
         
         if ((INTERACTION == I_VICSEK_SHARK)&&(i==1)) 
         {
@@ -6596,6 +6792,7 @@ t_segment segment[NMAXSEGMENTS])
         particle[i].dirmean = 0.0;
         particle[i].mass_inv = 1.0/PARTICLE_MASS;
         particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT;
+        particle[i].charge = 0.0;
         particle[i].vx = 0.0;
         particle[i].vy = 0.0;
         px[i] = 0.0;

sub_part_billiard.c

@@ -110,8 +110,8 @@ int writetiff(char *filename, char *description, int x, int y, int width, int he
   glPixelStorei(GL_PACK_ALIGNMENT, 1);
 
   glReadPixels(x, y, width, height, GL_RGB, GL_UNSIGNED_BYTE, image);
-  TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32) width);
+  TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32_t) width);
-  TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32) height);
+  TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32_t) height);
   TIFFSetField(file, TIFFTAG_BITSPERSAMPLE, 8);
   TIFFSetField(file, TIFFTAG_COMPRESSION, compression);
   TIFFSetField(file, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB);

sub_part_billiard_phasespace.c

@@ -112,8 +112,8 @@ int writetiff(char *filename, char *description, int x, int y, int width, int he
   glPixelStorei(GL_PACK_ALIGNMENT, 1);
 
   glReadPixels(x, y, width, height, GL_RGB, GL_UNSIGNED_BYTE, image);
-  TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32) width);
+  TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32_t) width);
-  TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32) height);
+  TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32_t) height);
   TIFFSetField(file, TIFFTAG_BITSPERSAMPLE, 8);
   TIFFSetField(file, TIFFTAG_COMPRESSION, compression);
   TIFFSetField(file, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB);

sub_perco.c

@@ -50,8 +50,8 @@ int writetiff(char *filename, char *description, int x, int y, int width, int he
   glPixelStorei(GL_PACK_ALIGNMENT, 1);
 
   glReadPixels(x, y, width, height, GL_RGB, GL_UNSIGNED_BYTE, image);
-  TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32) width);
+  TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32_t) width);
-  TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32) height);
+  TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32_t) height);
   TIFFSetField(file, TIFFTAG_BITSPERSAMPLE, 8);
   TIFFSetField(file, TIFFTAG_COMPRESSION, compression);
   TIFFSetField(file, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB);

sub_wave.c

@@ -30,8 +30,8 @@ int writetiff_new(char *filename, char *description, int x, int y, int width, in
   glPixelStorei(GL_PACK_ALIGNMENT, 1);
 
   glReadPixels(x, y, width, height, GL_RGB, GL_UNSIGNED_BYTE, image);
-  TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32) width);
+  TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32_t) width);
-  TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32) height);
+  TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32_t) height);
   TIFFSetField(file, TIFFTAG_BITSPERSAMPLE, 8);
   TIFFSetField(file, TIFFTAG_COMPRESSION, compression);
   TIFFSetField(file, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB);
@@ -85,8 +85,8 @@ int writetiff(char *filename, char *description, int x, int y, int width, int he
   glPixelStorei(GL_PACK_ALIGNMENT, 1);
 
   glReadPixels(x, y, width, height, GL_RGB, GL_UNSIGNED_BYTE, image);
-  TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32) width);
+  TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32_t) width);
-  TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32) height);
+  TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32_t) height);
   TIFFSetField(file, TIFFTAG_BITSPERSAMPLE, 8);
   TIFFSetField(file, TIFFTAG_COMPRESSION, compression);
   TIFFSetField(file, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB);
@@ -3130,6 +3130,13 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_
             }
             return(1);
         }
+        case (D_LENS):
+        {
+            if (vabs(y) > MU) return(1);
+            a = LAMBDA - 0.25;
+            if (x > 0) return ((x+a)*(x+a) + y*y > LAMBDA*LAMBDA);
+            return ((x-a)*(x-a) + y*y > LAMBDA*LAMBDA);
+        }
         case (D_MENGER):       
         {
             x1 = 0.5*(x+1.0);
@@ -4537,6 +4544,22 @@ void draw_billiard(int fade, double fade_value)      /* draws the billiard bound
             }
             break;
         }
+        case (D_LENS):
+        {
+            a = LAMBDA - 0.25;
+            if (first)
+            {
+                arcangle = asin(MU/LAMBDA);
+                ll = LAMBDA*cos(arcangle) - a;
+                first = 0;
+            }
+            
+            draw_circle_arc(-a, 0.0, LAMBDA, -arcangle, 2.0*arcangle, NSEG);
+            draw_line(ll, MU, -ll, MU);
+            draw_circle_arc(a, 0.0, LAMBDA, PI-arcangle, 2.0*arcangle, NSEG);
+            draw_line(-ll, -MU, ll, -MU);
+            break;
+        }
         case (D_MENGER):
         {
             glLineWidth(3);
@@ -6330,7 +6353,7 @@ void init_wave_packets(t_wave_packet *packet, int radius)
     
     printf("Initialising wave packets\n");
     switch (WAVE_PACKET_SOURCE_TYPE) {
-        case (0): 
+        case (WP_RANDOM1): 
         {
             nx = (int)sqrt((double)N_WAVE_PACKETS);
             ny = N_WAVE_PACKETS/nx;
@@ -6356,15 +6379,12 @@ void init_wave_packets(t_wave_packet *packet, int radius)
             
             break;
         }
-        case (1): 
+        case (WP_RANDOM2): 
         {
             for (i=0; i<N_WAVE_PACKETS; i++)
             {
-//                 j = i/nx;
-//                 k = i%nx;
                 packet[i].xc = XMIN + 0.15*(XMAX - XMIN)*(double)rand()/RAND_MAX;
                 packet[i].yc = 0.4*(YMAX - YMIN)*((double)rand()/RAND_MAX - 0.5);
-//                 packet[i].period = 30.0*(1.0 + 0.5*(double)rand()/RAND_MAX);
                 packet[i].period = 50.0*(1.0 + 0.5*(double)rand()/RAND_MAX);
                 packet[i].amp = INITIAL_AMP;
                 packet[i].phase = DPI*(double)rand()/RAND_MAX;
@@ -6377,6 +6397,27 @@ void init_wave_packets(t_wave_packet *packet, int radius)
                 packet[i].iy = ij[1]; 
             }
             
+            break;
+        }
+        case (WP_PAIR):
+        {
+            for (i=0; i<2; i++)
+            {
+                packet[i].xc = -1.25;
+                packet[i].yc = 0.06;
+                if (i==1) packet[i].yc *= -1.0;
+                packet[i].period = OSCILLATING_SOURCE_PERIOD;
+                packet[i].amp = INITIAL_AMP;
+                packet[i].phase = 0.0;
+                packet[i].var_envelope = 100000.0;
+                packet[i].time_shift = 0;
+                
+                xy_to_ij(packet[i].xc, packet[i].yc, ij);
+                if(ij[0] <= radius) ij[0] = radius+1;
+                packet[i].ix = ij[0];
+                packet[i].iy = ij[1]; 
+            }
+            
             break;
         }
     }
@@ -6390,18 +6431,23 @@ double wave_packet_height(double t, t_wave_packet packet, int type_envelope)
     cwave = packet.amp*cos(DPI*t/packet.period + packet.phase);
     
     switch (type_envelope) {
-        case (0): 
+        case (WE_SINE): 
         {
             envelope = 0.1 + sin(DPI*t/packet.var_envelope);
             envelope = envelope*envelope;
             break;
         }
-        case (1): 
+        case (WE_CUTOFF): 
         {
             if (t < (double)packet.var_envelope) envelope = 1.0;
             else envelope = 0.0;
             break;
         }
+        case (WE_CONSTANT):
+        {
+            envelope = 1.0;
+            break;
+        }
     }
     
     return(cwave*envelope);

wave_energy.c

@@ -219,16 +219,11 @@
 #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 */
+#define OSCILLATING_SOURCE_PERIOD 20    /* period of oscillating source */
 /* end of constants only used by sub_wave and sub_maze */
 
 #include "global_pdes.c"        /* constants and global variables */

wave_sphere.c

@@ -53,6 +53,7 @@
 
 #define DPOLE 20         /* safety distance to poles */
 #define SMOOTHPOLE 0.1     /* smoothing coefficient at poles */
+#define ZERO_MERIDIAN 0.0     /* choice of zero meridian (will be at left/right boundary of 2d plot) */
 
 #define XMIN -2.0
 #define XMAX 2.0	/* x interval */
@@ -61,10 +62,10 @@
 
 #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_SCALE 0.6  /* scaling for Julia sets */
-#define JULIA_ROT 90.0       /* rotation of Julia set, in degrees */
+#define JULIA_ROT -20.0       /* rotation of Julia set, in degrees */
-#define JULIA_RE -0.77145    
+#define JULIA_RE 0.5    
-#define JULIA_IM -0.10295    /* parameters for Julia sets */
+#define JULIA_IM 0.462    /* parameters for Julia sets */
 
 /* Choice of the billiard table */
 
@@ -77,7 +78,7 @@
 #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 10              /* 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 */
 
@@ -114,7 +115,7 @@
 
 /* Physical parameters of wave equation */
 
-#define TWOSPEEDS 0          /* set to 1 to replace hardcore boundary by medium with different speed */
+#define TWOSPEEDS 1          /* set to 1 to replace hardcore boundary by medium with different speed */
 #define OSCILLATE_LEFT 0     /* set to 1 to add oscilating boundary condition on the left */
 #define OSCILLATE_TOPBOT 0   /* set to 1 to enforce a planar wave on top and bottom boundary */
 #define OSCILLATION_SCHEDULE 3  /* oscillation schedule, see list in global_pdes.c */
@@ -124,7 +125,7 @@
 #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 COURANTB 0.002    /* 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 */
@@ -155,7 +156,7 @@
 
 /* Parameters for length and speed of simulation */
 
-#define NSTEPS 3600       /* number of frames of movie */
+#define NSTEPS 2500       /* 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 */
@@ -167,15 +168,14 @@
 #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 END_FRAMES 500   /* 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_AMP 0.5            /* amplitude of initial condition */
-#define INITIAL_VARIANCE 0.0005  /* variance of initial condition */
+#define INITIAL_VARIANCE 0.001  /* variance of initial condition */
-#define INITIAL_WAVELENGTH  0.025  /* wavelength of initial condition */
+#define INITIAL_WAVELENGTH  0.05  /* wavelength of initial condition */
 
 /* Plot type, see list in global_pdes.c  */
 
@@ -188,15 +188,15 @@
 #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_3D 0              /* 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_2D 1              /* 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 FADE_IN_OBSTACLE 0      /* 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 */
@@ -208,16 +208,19 @@
 /* 3D representation */
 
 #define REPRESENTATION_3D 1     /* choice of 3D representation */ 
-#define PLOT_2D 0               /* switch to 2D representation, equirectangular projection */
+#define PLOT_2D 1               /* switch to 2D representation, equirectangular projection */
 #define PHISHIFT 0.0            /* shift of phi in 2D plot (in degrees) */
+#define FLOODING 1              /* set to 1 to draw waves above altitude (for Earth representations) */
 
 #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 ROTATE_ANGLE 75.0   /* total angle of rotation during simulation */
+#define ROTATE_VIEW_WHILE_FADE 1    /* set to 1 to keep rotating viewpoint during fade */
 
-#define VIEWPOINT_TRAJ 1    /* type of viewpoint trajectory */
+#define VIEWPOINT_TRAJ 2    /* type of viewpoint trajectory */
+#define MAX_LATITUDE 45.0   /* maximal latitude for viewpoint trajectory VP_ORBIT2 */
 
 /* Color schemes */
 
@@ -233,18 +236,18 @@
 
 #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_AMPLITUDE 1.5   /* additional scaling factor for color scheme P_3D_AMPLITUDE */
-#define VSCALE_ENERGY 10.0     /* additional scaling factor for color scheme P_3D_ENERGY */
+#define VSCALE_ENERGY 4.0     /* additional scaling factor for color scheme P_3D_ENERGY */
 #define PHASE_FACTOR 20.0      /* factor in computation of phase in color scheme P_3D_PHASE */
 #define PHASE_SHIFT 0.0      /* shift of phase in color scheme P_3D_PHASE */
 #define ATTENUATION 0.0    /* exponential attenuation coefficient of contrast with time */
-#define E_SCALE 100.0       /* scaling factor for energy representation */
+#define E_SCALE 150.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_SCALE 1200.0    /* scaling factor for energy flux representation */
-#define FLUX_CSCALE 5.0      /* scaling factor for color in energy flux representation */
+#define FLUX_CSCALE 2.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 */
@@ -261,9 +264,9 @@
 #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 DRAW_COLOR_SCHEME 0       /* set to 1 to plot the color scheme */
-#define COLORBAR_RANGE 3.0      /* scale of color scheme bar */
+#define COLORBAR_RANGE 2.0      /* scale of color scheme bar */
-#define COLORBAR_RANGE_B 2.0    /* scale of color scheme bar for 2nd part */
+#define COLORBAR_RANGE_B 3.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 */
@@ -276,7 +279,6 @@
 #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 */
@@ -286,18 +288,27 @@
 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 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 */ 
+double observer[3] = {-5.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 ADD_DEM 1               /* add DEM (digital elevation model) */
+#define ADD_NEGATIVE_DEM 0      /* add DEM with bathymetric data */
+#define RSCALE_DEM 0.1          /* scaling factor of radial component for DEM */
+#define SMOOTH_DEM 0            /* set to 1 to smoothen DEM (to make altitude less constant) */
+#define DEM_SMOOTH_STEPS 10     /* number of smoothening steps */
+#define DEM_SMOOTH_HEIGHT 0.5   /* relative height below which to smoothen */
+#define DEM_MAXHEIGHT 8000      /* max height of DEM (estimated from Everest) */
+#define PLANET_SEALEVEL 2500.0      /* sea level for flooded planet */
+
 #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 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.75        /* limit on cosine of normal to shown facets */
+#define COS_VISIBLE -0.5        /* limit on cosine of normal to shown facets */
 
 #include "global_pdes.c"        /* constants and global variables */
 #include "global_3d.c"          /* constants and global variables */
@@ -498,25 +509,36 @@ void evolve_wave(double phi[NX*NY], double psi[NX*NY], double tmp[NX*NY], short
 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.12;
 //     double width = 0.2;
     
     width *= (double)NX/(double)WINWIDTH;
     
-    if (ROTATE_COLOR_SCHEME) 
+    if (PLOT_2D)
-        draw_color_scheme_palette_3d(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range, palette, fade, fade_value);
+    {
-    else if (circular)
+        if (ROTATE_COLOR_SCHEME) 
-        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);
+            draw_color_scheme_palette_fade(XMIN + 0.6, YMIN + 0.1, XMAX - 0.6, YMIN + 0.1 + width, plot, -range, range, palette, fade, fade_value);
-    else 
+        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);
+            draw_color_scheme_palette_fade(XMAX - 1.5*width, YMIN + 0.1, XMAX - 0.5*width, YMAX - 0.1, plot, -range, range, palette, fade, fade_value);
+    }
+    else
+    {
+        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;
+    double angle, ca, sa, r1, interpolate, rho;
-    static double observer_initial[3], r, ratio;
+    static double observer_initial[3], r, ratio, rho0, zmax;
     static int first = 1;
     
     if (first)
@@ -525,10 +547,14 @@ void viewpoint_schedule(int i)
         r1 = observer[0]*observer[0] + observer[1]*observer[1];
         r = sqrt(r1 + observer[2]*observer[2]);
         ratio = r/sqrt(r1);
+        rho0 = module2(observer[0], observer[1]);
+        if (vabs(rho0) < 0.001) rho0 = 0.001; 
+        zmax = r*sin(MAX_LATITUDE*PI/180.0);
         first = 0;
     }
     
-    angle = (ROTATE_ANGLE*DPI/360.0)*(double)i/(double)NSTEPS;
+    interpolate = (double)i/(double)NSTEPS;
+    angle = (ROTATE_ANGLE*DPI/360.0)*interpolate;
     ca = cos(angle);
     sa = sin(angle);
     switch (VIEWPOINT_TRAJ)
@@ -546,6 +572,21 @@ void viewpoint_schedule(int i)
             observer[2] = ca*observer_initial[2];
             break;
         }
+        case (VP_ORBIT2):
+        {
+            observer[0] = ca*observer_initial[0] - sa*observer_initial[1]*ratio;
+            observer[1] = ca*observer_initial[1] + sa*observer_initial[0]*ratio;
+            observer[2] = sa*zmax;
+            break;
+        }
+        case (VP_POLAR):
+        {
+            rho = -sa*observer_initial[2] + ca*rho0;
+            observer[0] = observer_initial[0]*rho/rho0;
+            observer[1] = observer_initial[1]*rho/rho0;
+            observer[2] = ca*observer_initial[2] + sa*rho0;
+            break; 
+        }
     }
     
     printf("Angle %.3lg, Observer position (%.3lg, %.3lg, %.3lg)\n", angle, observer[0], observer[1], observer[2]);
@@ -630,10 +671,7 @@ void animation()
     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;
     
@@ -674,16 +712,18 @@ void animation()
 //     for (j=1; j<NPOLY; j++)
 //         add_circular_wave_mod(1.0, lambda1*cos(((double)j+0.5)*angle), lambda1*sin(((double)j+0.5)*angle), phi, psi, xy_in);
         
-    init_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);
+    init_circular_wave_sphere(123.3*PI/180.0, -48.8*PI/180.0, phi, psi, xy_in, wsphere);
-//     add_circular_wave_sphere(-1.0, 0.25 + 3.0*PID, 0.0, phi, psi, xy_in, wsphere);
+    
-
+//     add_circular_wave_sphere(1.0, 1.0 - PI/9.0 + DPI/3.0, 0.15, phi, psi, xy_in, wsphere);
+//     add_circular_wave_sphere(1.0, 1.0 - PI/9.0 + 2.0*DPI/3.0, 0.15, phi, psi, xy_in, wsphere);
     
 //     printf("Wave initialized\n");
 
     /* initialize table of wave speeds/dissipation */
-    init_speed_dissipation(xy_in, tc, tcc, tgamma);
+    init_speed_dissipation_sphere(xy_in, tc, tcc, tgamma, wsphere);
+    
+    
     
     /* initialze potential to add to z coordinate */ 
     if (ADD_POTENTIAL)