daviddoji/YouTube-simulations
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Add files via upload

Browse Source 
Main changes: lennardjones, wave_billiard, rde
This commit is contained in:
Nils Berglund 2025-11-02 19:13:40 +01:00 committed by GitHub
parent cade2054f3
commit ab63c4d200
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
22 changed files with 4706 additions and 1179 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
colors_waves.c
View File

@ -3,7 +3,16 @@
double color_amplitude(double value, double scale, int time) double color_amplitude(double value, double scale, int time)
/* transforms the wave amplitude into a double in [-1,1] to feed into color scheme */ /* transforms the wave amplitude into a double in [-1,1] to feed into color scheme */
{ {
return(tanh(SLOPE*value/scale)*exp(-((double)time*ATTENUATION))); if (ATTENUATION == 0.0)
{
if (COLOR_RANGE == 1.0) return(tanh(SLOPE*value/scale));
else return(COLOR_RANGE*tanh(SLOPE*value/scale));
}
else
{
if (COLOR_RANGE == 1.0) return(tanh(SLOPE*value/scale)*exp(-((double)time*ATTENUATION)));
else return(COLOR_RANGE*tanh(SLOPE*value/scale)*exp(-((double)time*ATTENUATION)));
}
} }
double color_amplitude_asym(double value, double scale, int time) double color_amplitude_asym(double value, double scale, int time)

8
global_3d.c
View File

@ -13,8 +13,9 @@
#define E_EULER_COMP 7 /* compressible Euler equation */ #define E_EULER_COMP 7 /* compressible Euler equation */
#define E_SHALLOW_WATER 8 /* shallow water equation */ #define E_SHALLOW_WATER 8 /* shallow water equation */
#define E_KURAMOTO 9 /* Kuramoto-type equation (nearest neighbor coupling) */ #define E_KURAMOTO 9 /* Kuramoto-type equation (nearest neighbor coupling) */
#define E_KURAMOTO_SPHERE 91 /* Kuramoto on the sphere, with Poisson grid */ #define E_KURAMOTO_SPHERE 91 /* Kuramoto on the sphere, with cubic grid */
#define E_KELLER_SEGEL 10 /* Keller-Segel model */ #define E_KELLER_SEGEL 10 /* Keller-Segel model */
#define E_KELLER_SEGEL_SPHERE 11 /* Keller-Segel model on the sphere, cubic grid */
/* Choice of potential */ /* Choice of potential */
@ -109,6 +110,7 @@ int global_time = 0;
double max_depth = 1.0; double max_depth = 1.0;
int moon_position; int moon_position;
int ngridpoints; int ngridpoints;
int nfloorchi = 0;
/* structure used for color and height representations */ /* structure used for color and height representations */
/* possible extra fields: zfield, cfield, interpolated coordinates */ /* possible extra fields: zfield, cfield, interpolated coordinates */
@ -185,14 +187,14 @@ typedef struct
double x2d, y2d; /* x and y coordinates for 2D representation */ double x2d, y2d; /* x and y coordinates for 2D representation */
double altitude; /* altitude in case of Earth with digital elevation model */ double altitude; /* altitude in case of Earth with digital elevation model */
double cos_angle; /* cosine of light angle */ double cos_angle; /* cosine of light angle */
double cos_angle_sphere; /* cosing of light angle for perfect sphere */ double cos_angle_sphere; /* cosine of light angle for perfect sphere */
double force; /* external forcing */ double force; /* external forcing */
double phigrid, thetagrid; /* phi, theta angles for alt simulation grid on sphere */ double phigrid, thetagrid; /* phi, theta angles for alt simulation grid on sphere */
short int nneighb; /* number of neighbours, for Kuramoto model on sphere */ short int nneighb; /* number of neighbours, for Kuramoto model on sphere */
int neighbor[NMAX_SPHERE_NEIGHB]; /* list of neighbours */ int neighbor[NMAX_SPHERE_NEIGHB]; /* list of neighbours */
int convert_grid; /* convert field from simulation grid to longitude-latitude */ int convert_grid; /* convert field from simulation grid to longitude-latitude */
short int edge; /* has value 1 on edges of cubic simulation grid */ short int edge; /* has value 1 on edges of cubic simulation grid */
// short int jcoord; /* j coordinate of grid */ double cos_grid, sin_grid; /* cosine and sine of grid angle */
} t_wave_sphere; } t_wave_sphere;

53
global_ljones.c
View File

@ -1,5 +1,7 @@
/* Global variables and parameters for lennardjones */ /* Global variables and parameters for lennardjones */
#define HASHMAXNEIGH 15 /* max number of neighbouring hash grid cells */
/* Basic math */ /* Basic math */
#define PI 3.141592654 #define PI 3.141592654
@ -80,6 +82,7 @@
/* pattern of additional repelling segments */ /* pattern of additional repelling segments */
#define S_RECTANGLE 0 /* segments forming a rectangle */ #define S_RECTANGLE 0 /* segments forming a rectangle */
#define S_CYLINDRICAL 100 /* two lines at top and bottom */
#define S_CUP 1 /* segments forming a cup (for increasing gravity) */ #define S_CUP 1 /* segments forming a cup (for increasing gravity) */
#define S_HOURGLASS 2 /* segments forming an hour glass */ #define S_HOURGLASS 2 /* segments forming an hour glass */
#define S_PENTA 3 /* segments forming a pentagon with 3 angles of 120° and 2 right angles */ #define S_PENTA 3 /* segments forming a pentagon with 3 angles of 120° and 2 right angles */
@ -154,6 +157,8 @@
#define I_SEGMENT_CHARGED 161 /* harmonic interaction between segments and Coulomb interaction between ends*/ #define I_SEGMENT_CHARGED 161 /* harmonic interaction between segments and Coulomb interaction between ends*/
#define I_POLYGON 17 /* harmonic interaction between regular polygons */ #define I_POLYGON 17 /* harmonic interaction between regular polygons */
#define I_POLYGON_ALIGN 171 /* harmonic interaction between polygons with an aligning torque */ #define I_POLYGON_ALIGN 171 /* harmonic interaction between polygons with an aligning torque */
#define I_LJ_SPHERE 20 /* Lennard-Jones interaction on a sphere, redundant with 1 */
#define I_LJ_DIPOLE_SPHERE 25 /* Lennard-Jones with a dipolar angle dependence on a sphere */
/* Boundary conditions */ /* Boundary conditions */
@ -175,6 +180,7 @@
#define BC_REFLECT_ABS 21 /* reflecting on lower boundary, and "no-return" boundary conditions outside BC area */ #define BC_REFLECT_ABS 21 /* reflecting on lower boundary, and "no-return" boundary conditions outside BC area */
#define BC_REFLECT_ABS_BOTTOM 22 /* absorbing on lower boundary, and reflecting elsewhere */ #define BC_REFLECT_ABS_BOTTOM 22 /* absorbing on lower boundary, and reflecting elsewhere */
#define BC_REFLECT_ABS_RIGHT 23 /* absorbing on right boundary, and reflecting elsewhere */ #define BC_REFLECT_ABS_RIGHT 23 /* absorbing on right boundary, and reflecting elsewhere */
#define BC_SPHERE 30 /* particles on a sphere */
/* Regions for partial thermostat couplings */ /* Regions for partial thermostat couplings */
@ -259,6 +265,9 @@
#define CHEM_POLYGON_AGGREGATION 26 /* aggregation of polygons */ #define CHEM_POLYGON_AGGREGATION 26 /* aggregation of polygons */
#define CHEM_POLYGON_CLUSTER 261 /* clustering of polygons into new clusters */ #define CHEM_POLYGON_CLUSTER 261 /* clustering of polygons into new clusters */
#define CHEM_POLYGON_ONECLUSTER 262 /* clustering of polygons, with only one cluster allowed */ #define CHEM_POLYGON_ONECLUSTER 262 /* clustering of polygons, with only one cluster allowed */
#define CHEM_FISSION 27 /* nuclear fission reaction, U235 + n -> Sr95 + Xe139 + 2n */
#define CHEM_FISSION_CD 271 /* nuclear fission 27 and Cd112 + n -> Cd113 */
#define CHEM_FISSION_ABS 272 /* nuclear fission 27 and some extra moderation */
/* Initial conditions for chemical reactions */ /* Initial conditions for chemical reactions */
@ -320,6 +329,7 @@
#define P_COLLISION 22 /* colors depend on number of collision/reaction */ #define P_COLLISION 22 /* colors depend on number of collision/reaction */
#define P_RADIUS 23 /* colors depend on particle radius */ #define P_RADIUS 23 /* colors depend on particle radius */
#define P_MOL_ANG_MOMENTUM 24 /* colors depend on angular momentum of cluster */ #define P_MOL_ANG_MOMENTUM 24 /* colors depend on angular momentum of cluster */
#define P_DENSITY 25 /* colors represent local density (weighted average of neighbours */
/* Rotation schedules */ /* Rotation schedules */
@ -369,14 +379,17 @@
#define BG_NONE 0 /* no background color */ #define BG_NONE 0 /* no background color */
#define BG_DENSITY 1 /* background color depends on number of particles */ #define BG_DENSITY 1 /* background color depends on number of particles */
#define BG_NEIGHBOURS 11 /* background color depends on number of neighbours */
#define BG_CHARGE 2 /* background color depends on charge density */ #define BG_CHARGE 2 /* background color depends on charge density */
#define BG_EKIN 3 /* background color depends on kinetic energy */ #define BG_EKIN 3 /* background color depends on kinetic energy */
#define BG_LOG_EKIN 31 /* background color depends on log of kinetic energy */
#define BG_FORCE 4 /* background color depends on total force */ #define BG_FORCE 4 /* background color depends on total force */
#define BG_EOBSTACLES 5 /* background color depends on obstacle energy */ #define BG_EOBSTACLES 5 /* background color depends on obstacle energy */
#define BG_EKIN_OBSTACLES 6 /* background color depends on kinetic energy plus obstacle energy */ #define BG_EKIN_OBSTACLES 6 /* background color depends on kinetic energy plus obstacle energy */
#define BG_DIR_OBSTACLES 7 /* background color depends on direction of velocity of obstacles */ #define BG_DIR_OBSTACLES 7 /* background color depends on direction of velocity of obstacles */
#define BG_POS_OBSTACLES 8 /* background color depends on displacement of obstacles */ #define BG_POS_OBSTACLES 8 /* background color depends on displacement of obstacles */
#define BG_CURRENTX 9 /* background color depends on x-component of current */ #define BG_CURRENTX 9 /* background color depends on x-component of current */
#define BG_DIRECTION 10 /* background color depends on particle orientation */
/* Obstacle color schemes */ /* Obstacle color schemes */
@ -389,6 +402,13 @@
#define ADD_RECTANGLE 0 /* rectangular region, defined by ADDXMIN, etc */ #define ADD_RECTANGLE 0 /* rectangular region, defined by ADDXMIN, etc */
#define ADD_RING 1 /* ring_shaped region, defined by ADDRMIN, ADDRMAX */ #define ADD_RING 1 /* ring_shaped region, defined by ADDRMIN, ADDRMAX */
/* Type of rotating viewpoint */
#define VP_HORIZONTAL 0 /* rotate in a horizontal plane (constant latitude) */
#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 */
/* Color schemes */ /* Color schemes */
#define C_LUM 0 /* color scheme modifies luminosity (with slow drift of hue) */ #define C_LUM 0 /* color scheme modifies luminosity (with slow drift of hue) */
@ -452,6 +472,7 @@ typedef struct
double spin_freq; /* angular frequency of spin-spin interaction */ double spin_freq; /* angular frequency of spin-spin interaction */
double color_hue; /* color hue of particle, for P_INITIAL_POS plot type */ double color_hue; /* color hue of particle, for P_INITIAL_POS plot type */
int color_rgb[3]; /* RGB colors code of particle, for use in ljones_movie.c */ int color_rgb[3]; /* RGB colors code of particle, for use in ljones_movie.c */
double rgb[3], rgbx[3], rgby[3]; /* particle colors */
int partner[NMAXPARTNERS]; /* partner particles for option PAIR_PARTICLES */ int partner[NMAXPARTNERS]; /* partner particles for option PAIR_PARTICLES */
short int npartners; /* number of partner particles */ short int npartners; /* number of partner particles */
double partner_eqd[NMAXPARTNERS]; /* equilibrium distances between partners */ double partner_eqd[NMAXPARTNERS]; /* equilibrium distances between partners */
@ -499,10 +520,12 @@ typedef struct
int nobs; /* number of obstacles in cell */ int nobs; /* number of obstacles in cell */
int obstacle; /* obstacle number */ int obstacle; /* obstacle number */
int nneighb; /* number of neighbouring cells */ int nneighb; /* number of neighbouring cells */
int neighbour[9]; /* numbers of neighbouring cells */ int neighbour[HASHMAXNEIGH]; /* numbers of neighbouring cells */
double x1, y1, x2, y2; /* coordinates of hashcell corners */ double x1, y1, x2, y2; /* coordinates of hashcell corners */
double hue1, hue2; /* color hues */ double hue1, hue2; /* color hues */
double r, g, b; /* backgroud color RGB code */
double charge; /* charge of fixed obstacles */ double charge; /* charge of fixed obstacles */
double area; /* hash cell area (for sphere) */
} t_hashgrid; } t_hashgrid;
typedef struct typedef struct
@ -599,6 +622,7 @@ typedef struct
typedef struct typedef struct
{ {
int active; /* set to 1 for tracer to be active */
double xc, yc; /* center of circle */ double xc, yc; /* center of circle */
} t_tracer; } t_tracer;
@ -624,7 +648,7 @@ typedef struct
{ {
double x, y; /* location of collision */ double x, y; /* location of collision */
int time; /* time since collision */ int time; /* time since collision */
int color; /* color hue in case of different collisions */ double color; /* color hue in case of different collisions */
} t_collision; } t_collision;
typedef struct typedef struct
@ -672,6 +696,31 @@ typedef struct
} t_lj_parameters; } t_lj_parameters;
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 reg_cottheta; /* regularized cotangent of theta */
double x, y, z; /* x, y, z coordinates of point on sphere */
double radius; /* radius with wave height */
// double radius_dem; /* radius with digital elevation model */
double r, g, b; /* RGB values for image */
// short int indomain; /* has value 1 if lattice point is in domain */
// short int draw_wave; /* has value 1 if wave instead of DEM is drawn */
// short int evolve_wave; /* has value 1 where there is wave evolution */
// 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; /* cosine of light angle for perfect sphere */
// double force; /* external forcing */
// double phigrid, thetagrid; /* phi, theta angles for alt simulation grid on sphere */
// short int nneighb; /* number of neighbours, for Kuramoto model on sphere */
// int neighbor[NMAX_SPHERE_NEIGHB]; /* list of neighbours */
// int convert_grid; /* convert field from simulation grid to longitude-latitude */
// short int edge; /* has value 1 on edges of cubic simulation grid */
// double cos_grid, sin_grid; /* cosine and sine of grid angle */
} t_lj_sphere;
int frame_time = 0, ncircles, nobstacles, nsegments, ngroups = 1, counter = 0, nmolecules = 0, nbelts = 0, n_tracers = 0, n_otriangles = 0, n_ofacets = 0; int frame_time = 0, ncircles, nobstacles, nsegments, ngroups = 1, counter = 0, nmolecules = 0, nbelts = 0, n_tracers = 0, n_otriangles = 0, n_ofacets = 0;
FILE *lj_log; FILE *lj_log;

13
global_pdes.c
View File

@ -14,6 +14,7 @@
#define D_NOTHING 999 /* no boundaries */ #define D_NOTHING 999 /* no boundaries */
#define D_RECTANGLE 0 /* rectangular domain */ #define D_RECTANGLE 0 /* rectangular domain */
#define D_ELLIPSE 1 /* elliptical domain */ #define D_ELLIPSE 1 /* elliptical domain */
#define D_CIRCLE_LAYERS 196 /* concentric circles */
#define D_EXT_ELLIPSE 199 /* exterior of elliptical domain */ #define D_EXT_ELLIPSE 199 /* exterior of elliptical domain */
#define D_EXT_ELLIPSE_CURVED 198 /* exterior of curved elliptical domain */ #define D_EXT_ELLIPSE_CURVED 198 /* exterior of curved elliptical domain */
#define D_EXT_ELLIPSE_CURVED_BDRY 197 /* exterior of curved elliptical domain, with horizontal boundaries */ #define D_EXT_ELLIPSE_CURVED_BDRY 197 /* exterior of curved elliptical domain, with horizontal boundaries */
@ -38,6 +39,8 @@
#define D_MENGER_ROTATED 17 /* rotated Menger-Sierpinski carpet */ #define D_MENGER_ROTATED 17 /* rotated Menger-Sierpinski carpet */
#define D_PARABOLA 18 /* parabolic domain */ #define D_PARABOLA 18 /* parabolic domain */
#define D_TWO_PARABOLAS 19 /* two facing parabolic antennas */ #define D_TWO_PARABOLAS 19 /* two facing parabolic antennas */
#define D_TWO_PARABOLAS_ASYM 191 /* confocal facing parabolas of different size */
#define D_TWO_PARABOLAS_ASYM_GUIDE 192 /* confocal facing parabolas of different size with a wave guide */
#define D_CIRCLES 20 /* several circles */ #define D_CIRCLES 20 /* several circles */
#define D_CIRCLES_IN_RECT 201 /* several circles in a rectangle */ #define D_CIRCLES_IN_RECT 201 /* several circles in a rectangle */
@ -131,6 +134,9 @@
#define D_CIRCLE_LATTICE_HONEY 95 /* honeycomb lattice of connected circles */ #define D_CIRCLE_LATTICE_HONEY 95 /* honeycomb lattice of connected circles */
#define D_CIRCLE_LATTICE_POISSON 96 /* Poisson disc process of connected circles */ #define D_CIRCLE_LATTICE_POISSON 96 /* Poisson disc process of connected circles */
#define D_CIRCLE_LATTICE_RHOMBUS 97 /* rhombus-based lattice of connected circles */ #define D_CIRCLE_LATTICE_RHOMBUS 97 /* rhombus-based lattice of connected circles */
#define D_DISC_WAVEGUIDE 98 /* a disc with an attached waveguide */
#define D_DISC_WAVEGUIDE_SHIFTED 981 /* a disc with a shifted attached waveguide */
#define D_DISC_WAVEGUIDE_ELLIPSE 982 /* a disc with a shifted attached waveguide */
/* for wave_sphere.c */ /* for wave_sphere.c */
@ -231,7 +237,9 @@
#define IOR_MICHELSON 17 /* Michelson interferometer, to use with D_MICHELSON */ #define IOR_MICHELSON 17 /* Michelson interferometer, to use with D_MICHELSON */
#define IOR_GRADIENT_INDEX_LENS 18 /* gradient index lens (parabolic c(r)^2) */ #define IOR_GRADIENT_INDEX_LENS 18 /* gradient index lens (parabolic c(r)^2) */
#define IOR_GRADIENT_INDEX_LENS_B 181 /* gradient index lens (parabolic c(r)) */ #define IOR_GRADIENT_INDEX_LENS_B 181 /* gradient index lens (parabolic c(r)) */
#define IOR_LINEAR_X_A 19 /* IoR depending linearly on x */ #define IOR_LUNEBURG_LENS 182 /* Luneburg lens */
#define IOR_LUNEBURG_LAYERS 183 /* Luneburg lens with layers */
#define IOR_LINEAR_X_A 19 /* IoR depending linearly on x */
#define IOR_LINEAR_X_B 191 /* IoR depending linearly on x, with correct boundary values */ #define IOR_LINEAR_X_B 191 /* IoR depending linearly on x, with correct boundary values */
#define IOR_EARTH_DEM 20 /* digital elevation model (for waves on sphere) */ #define IOR_EARTH_DEM 20 /* digital elevation model (for waves on sphere) */
@ -257,7 +265,10 @@
#define OSC_BEAM 4 /* periodic oscillation modulated by y cut-off */ #define OSC_BEAM 4 /* periodic oscillation modulated by y cut-off */
#define OSC_BEAM_GAUSSIAN 41 /* periodic oscillation modulated by Gaussian in y */ #define OSC_BEAM_GAUSSIAN 41 /* periodic oscillation modulated by Gaussian in y */
#define OSC_BEAM_SINE 42 /* periodic oscillation modulated by sine in y */ #define OSC_BEAM_SINE 42 /* periodic oscillation modulated by sine in y */
#define OSC_BEAM_SINE_DECREASING 43 /* decreasing oscillation modulated by sine in y */
#define OSC_BEAM_SINE_CHIRP 44 /* decreasing socillation modulated by sine in y */
#define OSC_BEAM_TWOPERIODS 5 /* sum of two periodic oscillations modulated by y cut-off */ #define OSC_BEAM_TWOPERIODS 5 /* sum of two periodic oscillations modulated by y cut-off */
#define OSC_BEAM_SINE_TWOPERIODS 51 /* sum of two periodic oscillations modulated by y cut-off */
#define OSC_TWO_WAVES 6 /* two linear waves at an angle, separate */ #define OSC_TWO_WAVES 6 /* two linear waves at an angle, separate */
#define OSC_TWO_WAVES_ADDED 61 /* two linear waves at an angle, superimposed */ #define OSC_TWO_WAVES_ADDED 61 /* two linear waves at an angle, superimposed */

1
heat.c
View File

@ -216,6 +216,7 @@
#define WAVE_PACKET_SOURCE_TYPE 1 /* type of wave packet sources */ #define WAVE_PACKET_SOURCE_TYPE 1 /* type of wave packet sources */
#define N_WAVE_PACKETS 15 /* number of wave packets */ #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 OSCIL_LEFT_YSHIFT 0.0 /* y-dependence of left oscillation (for non-horizontal waves) */
#define OSCIL_YMID -0.9 /* defines oscilling beam midpoint */
#define OSCILLATING_SOURCE_PERIOD 20 /* period of oscillating source */ #define OSCILLATING_SOURCE_PERIOD 20 /* period of oscillating source */
#define MU_B 1.0 /* parameter controlling the dimensions of domain */ #define MU_B 1.0 /* parameter controlling the dimensions of domain */
#define GAMMA 0.0 /* damping factor in wave equation */ #define GAMMA 0.0 /* damping factor in wave equation */

367
lennardjones.c
View File

@ -53,15 +53,15 @@
#define WINWIDTH 1760 /* window width */ #define WINWIDTH 1760 /* window width */
#define WINHEIGHT 990 /* window height */ #define WINHEIGHT 990 /* window height */
#define XMIN -2.0 #define XMIN 0.0
#define XMAX 2.0 /* x interval */ #define XMAX 6.283185307 /* x interval */
#define YMIN -1.125 #define YMIN 0.0
#define YMAX 1.125 /* y interval for 9/16 aspect ratio */ #define YMAX 3.141592654 /* y interval for 9/16 aspect ratio */
#define INITXMIN -0.1 #define INITXMIN 0.025
#define INITXMAX 0.1 /* x interval for initial condition */ #define INITXMAX 6.255 /* x interval for initial condition */
#define INITYMIN -0.1 #define INITYMIN 0.4
#define INITYMAX 0.1 /* y interval for initial condition */ #define INITYMAX 2.74 /* y interval for initial condition */
#define THERMOXMIN -1.25 #define THERMOXMIN -1.25
#define THERMOXMAX 1.25 /* x interval for initial condition */ #define THERMOXMAX 1.25 /* x interval for initial condition */
@ -72,13 +72,13 @@
#define ADDXMAX 1.9 /* x interval for adding particles */ #define ADDXMAX 1.9 /* x interval for adding particles */
#define ADDYMIN 1.2 #define ADDYMIN 1.2
#define ADDYMAX 1.3 /* y interval for adding particles */ #define ADDYMAX 1.3 /* y interval for adding particles */
#define ADDRMIN 4.75 #define ADDRMIN 2.0
#define ADDRMAX 6.0 /* r interval for adding particles */ #define ADDRMAX 2.1 /* r interval for adding particles */
#define BCXMIN -2.0 #define BCXMIN 0.0
#define BCXMAX 2.0 /* x interval for boundary condition */ #define BCXMAX 6.283185307 /* x interval for boundary condition */
#define BCYMIN -1.125 #define BCYMIN 0.05
#define BCYMAX 1.325 /* y interval for boundary condition */ #define BCYMAX 3.091592654 /* y interval for boundary condition */
#define OBSXMIN -2.0 #define OBSXMIN -2.0
#define OBSXMAX 2.0 /* x interval for motion of obstacle */ #define OBSXMAX 2.0 /* x interval for motion of obstacle */
@ -108,7 +108,7 @@
#define COUPLE_MINLENGTH 0.5 /* length at which bonds decouple */ #define COUPLE_MINLENGTH 0.5 /* length at which bonds decouple */
#define ADD_FIXED_SEGMENTS 0 /* set to 1 to add fixed segments as obstacles */ #define ADD_FIXED_SEGMENTS 0 /* set to 1 to add fixed segments as obstacles */
#define SEGMENT_PATTERN 15 /* pattern of repelling segments, see list in global_ljones.c */ #define SEGMENT_PATTERN 100 /* 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 3 /* shape of rocket combustion chamber, see list in global_ljones.c */
#define ROCKET_SHAPE_B 3 /* shape of second rocket */ #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 6 /* shape of nozzle, see list in global_ljones.c */
@ -119,14 +119,14 @@
#define OBSTACLE_OMEGA 300.0 /* obstacle rotation speed */ #define OBSTACLE_OMEGA 300.0 /* obstacle rotation speed */
#define TWO_TYPES 0 /* set to 1 to have two types of particles */ #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 TYPE_PROPORTION 0.75 /* proportion of particles of first type */
#define TWOTYPE_CONFIG 0 /* choice of types, see TTC_ list in global_ljones.c */ #define TWOTYPE_CONFIG 0 /* 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 SYMMETRIZE_FORCE 0 /* set to 1 to symmetrize two-particle interaction, only needed if particles are not all the same */
#define CENTER_PX 0 /* set to 1 to center horizontal momentum */ #define CENTER_PX 0 /* set to 1 to center horizontal momentum */
#define CENTER_PY 0 /* set to 1 to center vertical momentum */ #define CENTER_PY 0 /* set to 1 to center vertical momentum */
#define CENTER_PANGLE 0 /* set to 1 to center angular momentum */ #define CENTER_PANGLE 0 /* set to 1 to center angular momentum */
#define INTERACTION 12 /* particle interaction, 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 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 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 SPIN_INTER_FREQUENCY_B 4.0 /* angular frequency of spin-spin interaction for second particle type */
@ -139,9 +139,10 @@
#define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */ #define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */
#define LAMBDA 0.2 /* parameter controlling the dimensions of domain */ #define LAMBDA 0.2 /* parameter controlling the dimensions of domain */
#define MU 0.016 /* parameter controlling radius of particles */ #define MU 0.022 /* parameter controlling radius of particles */
#define MU_B 0.012 /* parameter controlling radius of particles of second type */ #define MU_B 0.022 /* parameter controlling radius of particles of second type */
#define MU_ADD 0.016 /* parameter controlling radius of added particles */ #define MU_ADD 0.022 /* parameter controlling radius of added particles */
#define MU_ADD_B 0.022 /* parameter controlling radius of added particles */
#define NPOLY 3 /* number of sides of polygon */ #define NPOLY 3 /* number of sides of polygon */
#define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */ #define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */
#define AWEDGE 0.5 /* opening angle of wedge, in units of Pi/2 */ #define AWEDGE 0.5 /* opening angle of wedge, in units of Pi/2 */
@ -150,8 +151,8 @@
#define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ #define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */
#define MANDELLIMIT 10.0 /* limit value for approximation of 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 FOCI 1 /* set to 1 to draw focal points of ellipse */
#define NGRIDX 1 /* number of grid point for grid of disks */ #define NGRIDX 40 /* number of grid point for grid of disks */
#define NGRIDY 1 /* number of grid point for grid of disks */ #define NGRIDY 20 /* number of grid point for grid of disks */
#define EHRENFEST_RADIUS 0.9 /* radius of container for Ehrenfest urn configuration */ #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 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 TWO_CIRCLES_RADIUS_RATIO 0.8 /* ratio of radii for S_TWO_CIRCLES_EXT segment configuration */
@ -167,8 +168,8 @@
/* Parameters for length and speed of simulation */ /* Parameters for length and speed of simulation */
#define NSTEPS 4800 /* number of frames of movie */ #define NSTEPS 2300 /* number of frames of movie */
#define NVID 100 /* number of iterations between images displayed on screen */ #define NVID 100 /* number of iterations between images displayed on screen */
#define NSEG 25 /* number of segments of boundary of circles */ #define NSEG 25 /* number of segments of boundary of circles */
#define INITIAL_TIME 0 /* time after which to start saving frames */ #define INITIAL_TIME 0 /* time after which to start saving frames */
#define OBSTACLE_INITIAL_TIME 0 /* time after which to start moving obstacle */ #define OBSTACLE_INITIAL_TIME 0 /* time after which to start moving obstacle */
@ -181,22 +182,22 @@
#define SLEEP1 1 /* initial sleeping time */ #define SLEEP1 1 /* initial sleeping time */
#define SLEEP2 1 /* final sleeping time */ #define SLEEP2 1 /* final sleeping time */
#define MID_FRAMES 100 /* number of still frames between parts of two-part movie */ #define MID_FRAMES 100 /* number of still frames between parts of two-part movie */
#define END_FRAMES 100 /* number of still frames at end of movie */ #define END_FRAMES 200 /* number of still frames at end of movie */
/* Boundary conditions, see list in global_ljones.c */ /* Boundary conditions, see list in global_ljones.c */
#define BOUNDARY_COND 1 #define BOUNDARY_COND 30
/* Plot type, see list in global_ljones.c */ /* Plot type, see list in global_ljones.c */
#define PLOT 13 #define PLOT 5
#define PLOT_B 17 /* plot type for second movie */ #define PLOT_B 13 /* plot type for second movie */
/* Background color depending on particle properties */ /* Background color depending on particle properties */
#define COLOR_BACKGROUND 0 /* set to 1 to color background */ #define COLOR_BACKGROUND 0 /* set to 1 to color background */
#define BG_COLOR 6 /* type of background coloring, see list in global_ljones.c */ #define BG_COLOR 10 /* type of background coloring, see list in global_ljones.c */
#define BG_COLOR_B 9 /* type of background coloring, see list in global_ljones.c */ #define BG_COLOR_B 3 /* type of background coloring, see list in global_ljones.c */
#define OBSTACLE_COLOR 0 /* type of obstacle, see OC_ in global_ljones.c */ #define OBSTACLE_COLOR 0 /* type of obstacle, see OC_ in global_ljones.c */
#define DRAW_BONDS 0 /* set to 1 to draw bonds between neighbours */ #define DRAW_BONDS 0 /* set to 1 to draw bonds between neighbours */
@ -249,8 +250,8 @@
#define HUEAMP -50.0 /* amplitude of variation of hue for color scheme C_HUE */ #define HUEAMP -50.0 /* amplitude of variation of hue for color scheme C_HUE */
#define COLOR_HUESHIFT 1.0 /* shift in color hue (for some cyclic palettes) */ #define COLOR_HUESHIFT 1.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_TEMPERATURE 1 /* set to 1 to print current temperature */
#define PRINT_ANGLE 0 /* set to 1 to print obstacle orientation */ #define PRINT_ANGLE 0 /* set to 1 to print obstacle orientation */
#define PRINT_OMEGA 0 /* 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_PARTICLE_SPEEDS 0 /* set to 1 to print average speeds/momenta of particles */
@ -267,7 +268,9 @@
#define PARTICLE_HUE_MIN 359.0 /* color of original particle */ #define PARTICLE_HUE_MIN 359.0 /* color of original particle */
#define PARTICLE_HUE_MAX 0.0 /* color of saturated particle */ #define PARTICLE_HUE_MAX 0.0 /* color of saturated particle */
#define PARTICLE_EMIN 0.0 /* energy of particle with coolest color */ #define PARTICLE_EMIN 0.0 /* energy of particle with coolest color */
#define PARTICLE_EMAX 50000.0 /* energy of particle with hottest color */ #define PARTICLE_EMAX 2000.0 /* energy of particle with hottest color */
#define PARTICLE_DMIN 200.0 /* energy of particle with largest local density */
#define PARTICLE_DMAX 500.0 /* energy of particle with largest local density */
#define SEGMENT_HUE_MIN 275.0 /* color of original segment */ #define SEGMENT_HUE_MIN 275.0 /* color of original segment */
#define SEGMENT_HUE_MAX 30.0 /* color of saturated segment */ #define SEGMENT_HUE_MAX 30.0 /* color of saturated segment */
#define OBSTACLE_EMAX 1000000.0 /* energy of obstacle with hottest color */ #define OBSTACLE_EMAX 1000000.0 /* energy of obstacle with hottest color */
@ -275,12 +278,14 @@
#define HUE_TYPE0 320.0 /* hue of particles of type 0 */ #define HUE_TYPE0 320.0 /* hue of particles of type 0 */
#define HUE_TYPE1 60.0 /* hue of particles of type 1 */ #define HUE_TYPE1 60.0 /* hue of particles of type 1 */
#define HUE_TYPE2 320.0 /* hue of particles of type 2 */ #define HUE_TYPE2 320.0 /* hue of particles of type 2 */
#define HUE_TYPE3 260.0 /* hue of particles of type 3 */ #define HUE_TYPE3 140.0 /* hue of particles of type 3 */
#define HUE_TYPE4 200.0 /* hue of particles of type 4 */ #define HUE_TYPE4 200.0 /* hue of particles of type 4 */
#define HUE_TYPE5 60.0 /* hue of particles of type 5 */ #define HUE_TYPE5 60.0 /* hue of particles of type 5 */
#define HUE_TYPE6 130.0 /* hue of particles of type 6 */ #define HUE_TYPE6 130.0 /* hue of particles of type 6 */
#define HUE_TYPE7 150.0 /* hue of particles of type 7 */ #define HUE_TYPE7 200.0 /* hue of particles of type 7 */
#define BG_FORCE_SLOPE 1.0e-6 /* contant in BG_FORCE backgound color scheme*/ #define BG_LOG_EKIN_SHIFT 1.0 /* constant in BG_LOG_EKIN background color scheme */
#define BG_FORCE_SLOPE 1.0e-6 /* constant in BG_FORCE backgound color scheme */
#define BG_CHARGE_SLOPE 0.75 /* constant in BG_CHARGE backgound color scheme (default: 0.5) */
#define PARTICLE_LMAX 1.5e4 /* angular momentum particle with brightest color */ #define PARTICLE_LMAX 1.5e4 /* angular momentum particle with brightest color */
#define RANDOM_RADIUS 0 /* set to 1 for random particle radius */ #define RANDOM_RADIUS 0 /* set to 1 for random particle radius */
@ -290,39 +295,42 @@
#define ADAPT_DAMPING_TO_RADIUS 0.0 /* set to positive value to for friction prop to power of radius */ #define ADAPT_DAMPING_TO_RADIUS 0.0 /* set to positive value to for friction prop to power of radius */
#define ADAPT_DAMPING_FACTOR 0.0 /* factor by which damping is adapted to radius */ #define ADAPT_DAMPING_FACTOR 0.0 /* factor by which damping is adapted to radius */
#define DT_PARTICLE 1.0e-6 /* time step for particle displacement */ #define DT_PARTICLE 1.0e-6 /* time step for particle displacement */
#define KREPEL 20.0 /* constant in repelling force between particles */ #define KREPEL 50.0 /* constant in repelling force between particles */
#define EQUILIBRIUM_DIST 4.0 /* Lennard-Jones equilibrium distance */ #define EQUILIBRIUM_DIST 3.5 /* Lennard-Jones equilibrium distance */
#define EQUILIBRIUM_DIST_B 4.0 /* Lennard-Jones equilibrium distance for second type of particle */ #define EQUILIBRIUM_DIST_B 3.5 /* Lennard-Jones equilibrium distance for second type of particle */
#define SEGMENT_FORCE_EQR 1.0 /* equilibrium distance factor for force from segments (default 1.5) */ #define SEGMENT_FORCE_EQR 1.0 /* equilibrium distance factor for force from segments (default 1.5) */
#define REPEL_RADIUS 25.0 /* radius in which repelling force acts (in units of particle radius) */ #define REPEL_RADIUS 25.0 /* radius in which repelling force acts (in units of particle radius) */
#define DAMPING 50.0 /* damping coefficient of particles */ #define DAMPING 0.0 /* damping coefficient of particles */
#define INITIAL_DAMPING 1000.0 /* damping coefficient of particles during initial phase */ #define INITIAL_DAMPING 1000.0 /* damping coefficient of particles during initial phase */
#define DAMPING_ROT 5.0 /* damping coefficient for rotation of particles */ #define DAMPING_ROT 0.0 /* damping coefficient for rotation of particles */
#define DAMPING_PAIRS 0.0 /* damping between paired particles */ #define DAMPING_PAIRS 0.0 /* damping between paired particles */
#define PARTICLE_MASS 2.0 /* 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_B */ #define PARTICLE_MASS_B 1.0 /* mass of particle of radius MU_B */
#define PARTICLE_ADD_MASS 2.0 /* mass of added particles */ #define PARTICLE_ADD_MASS 1.0 /* mass of added particles */
#define PARTICLE_ADD_MASS_B 1.0 /* mass of added particles */
#define PARTICLE_INERTIA_MOMENT 0.1 /* moment of inertia of particle */ #define PARTICLE_INERTIA_MOMENT 0.1 /* moment of inertia of particle */
#define PARTICLE_INERTIA_MOMENT_B 0.1 /* moment of inertia of second type of particle */ #define PARTICLE_INERTIA_MOMENT_B 0.1 /* moment of inertia of second type of particle */
#define V_INITIAL 100.0 /* initial velocity range */ #define V_INITIAL 25.0 /* initial velocity range */
#define V_INITIAL_ADD 100.0 /* initial velocity range for added particles */ #define V_INITIAL_ADD 0.0 /* initial velocity range for added particles */
#define OMEGA_INITIAL 100.0 /* initial angular velocity range */ #define OMEGA_INITIAL 100.0 /* initial angular velocity range */
#define VICSEK_VMIN 1.0 /* minimal speed of particles in Vicsek model */ #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 VICSEK_VMAX 40.0 /* minimal speed of particles in Vicsek model */
#define COULOMB_LJ_FACTOR 1.0 /* relative intensity of LJ interaction in I_COULOMB_LJ interaction (default: 0.01) */ #define COULOMB_LJ_FACTOR 1.0 /* relative intensity of LJ interaction in I_COULOMB_LJ interaction (default: 0.01) */
#define KCOULOMB_FACTOR 50.0 /* relative intensity of Coulomb interaction in I_COULOMB_LJ (default: 100.0) */
#define OBSTACLE_DAMPING 0.0 /* damping of oscillating obstacles */ #define OBSTACLE_DAMPING 0.0 /* damping of oscillating obstacles */
#define V_INITIAL_TYPE 0 /* type of initial speed distribution (see VI_ in global_ljones.c) */ #define V_INITIAL_TYPE 0 /* type of initial speed distribution (see VI_ in global_ljones.c) */
#define THERMOSTAT 0 /* set to 1 to switch on thermostat */ #define THERMOSTAT 1 /* set to 1 to switch on thermostat */
#define VARY_THERMOSTAT 0 /* set to 1 for time-dependent thermostat schedule */ #define VARY_THERMOSTAT 0 /* set to 1 for time-dependent thermostat schedule */
#define SIGMA 5.0 /* noise intensity in thermostat */ #define SIGMA 5.0 /* noise intensity in thermostat */
#define BETA 0.004 /* initial inverse temperature */ #define BETA 0.0002 /* initial inverse temperature */
#define MU_XI 0.005 /* friction constant in thermostat */ #define MU_XI 0.005 /* friction constant in thermostat */
#define KSPRING_BOUNDARY 2.0e11 /* confining harmonic potential outside simulation region */ #define KSPRING_BOUNDARY 2.0e11 /* confining harmonic potential outside simulation region */
#define KSPRING_OBSTACLE 5.0e11 /* harmonic potential of obstacles */ #define KSPRING_OBSTACLE 5.0e11 /* harmonic potential of obstacles */
#define NBH_DIST_FACTOR 4.0 /* radius in which to count neighbours */ #define NBH_DIST_FACTOR 6.0 /* radius in which to count neighbours */
#define GRAVITY 10000.0 /* gravity acting on all particles */ #define BOND_DIST_FACTOR 6.0 /* radius in which to draw bonds */
#define GRAVITY 0.0 /* gravity acting on all particles */
#define GRAVITY_X 0.0 /* horizontal gravity acting on all particles */ #define GRAVITY_X 0.0 /* horizontal gravity acting on all particles */
#define CIRCULAR_GRAVITY 0 /* set to 1 to have gravity directed to center */ #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 INCREASE_GRAVITY 0 /* set to 1 to increase gravity during the simulation */
@ -338,9 +346,10 @@
#define EFIELD_Y 0.0 /* value of electric field */ #define EFIELD_Y 0.0 /* value of electric field */
#define ADD_BFIELD 0 /* set to 1 to add a magnetic field */ #define ADD_BFIELD 0 /* set to 1 to add a magnetic field */
#define BFIELD 20000.0 /* value of magnetic field */ #define BFIELD 20000.0 /* value of magnetic field */
#define CHARGE 1.0 /* charge of particles of first type */ #define CHARGE 0.0 /* charge of particles of first type */
#define CHARGE_B -1.0 /* charge of particles of second type */ #define CHARGE_B 0.0 /* charge of particles of second type */
#define CHARGE_ADD 1.0 /* charge of added particles */ #define CHARGE_ADD 0.0 /* charge of added particles */
#define CHARGE_ADD_B 0.0 /* charge of added particles */
#define INCREASE_E 0 /* set to 1 to increase electric field */ #define INCREASE_E 0 /* set to 1 to increase electric field */
#define OSCILLATE_E 0 /* set to 1 for oscillating electric field */ #define OSCILLATE_E 0 /* set to 1 for oscillating electric field */
#define E_PERIOD 1000 /* period of oscillating electric field */ #define E_PERIOD 1000 /* period of oscillating electric field */
@ -367,10 +376,10 @@
#define WIND_FORCE 1.35e6 /* force of wind */ #define WIND_FORCE 1.35e6 /* force of wind */
#define WIND_YMIN -0.6 /* min altitude of region with wind */ #define WIND_YMIN -0.6 /* min altitude of region with wind */
#define ROTATION 0 /* set to 1 to include rotation of particles */ #define ROTATION 1 /* set to 1 to include rotation of particles */
#define COUPLE_ANGLE_TO_THERMOSTAT 0 /* set to 1 to couple angular degrees of freedom to thermostat */ #define COUPLE_ANGLE_TO_THERMOSTAT 0 /* 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 DIMENSION_FACTOR 0.25 /* scaling factor taking into account number of degrees of freedom */
#define KTORQUE 1.0e5 /* force constant in angular dynamics */ #define KTORQUE 2.0e3 /* force constant in angular dynamics */
#define KTORQUE_BOUNDARY 1.0e5 /* constant in torque from the boundary */ #define KTORQUE_BOUNDARY 1.0e5 /* constant in torque from the boundary */
#define KTORQUE_B 10.0 /* force constant in angular dynamics */ #define KTORQUE_B 10.0 /* force constant in angular dynamics */
#define KTORQUE_DIFF 500.0 /* force constant in angular dynamics for different particles */ #define KTORQUE_DIFF 500.0 /* force constant in angular dynamics for different particles */
@ -378,19 +387,20 @@
#define DRAW_SPIN_B 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 1 /* set to 1 to draw cross on particles of second type */
#define DRAW_MINUS 1 /* set to 1 to draw cross on particles of negative charge */ #define DRAW_MINUS 1 /* set to 1 to draw cross on particles of negative charge */
#define SPIN_RANGE 10.0 /* range of spin-spin interaction */ #define SPIN_RANGE 5.0 /* range of spin-spin interaction */
#define SPIN_RANGE_B 10.0 /* range of spin-spin interaction for second type of particle */ #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 QUADRUPOLE_RATIO 0.6 /* anisotropy in quadrupole potential */
#define INCREASE_BETA 0 /* set to 1 to increase BETA during simulation */ #define INCREASE_BETA 0 /* set to 1 to increase BETA during simulation */
#define BETA_SCHEDULE 3 /* type of temperature schedule, see TS_* in global_ljones */ #define BETA_SCHEDULE 5 /* type of temperature schedule, see TS_* in global_ljones */
#define BETA_FACTOR 0.06 /* factor by which to change BETA during simulation */ #define BETA_FACTOR 50.0 /* factor by which to change BETA during simulation */
#define TS_SLOPE 8.5 /* controls speed of change of BETA for TS_TANH schedule (default 1.0) */ #define TS_SLOPE 8.5 /* controls speed of change of BETA for TS_TANH schedule (default 1.0) */
#define N_TOSCILLATIONS 1.0 /* number of temperature oscillations in BETA schedule */ #define N_TOSCILLATIONS 1.5 /* number of temperature oscillations in BETA schedule */
#define NO_OSCILLATION 0 /* set to 1 to have exponential BETA change only */ #define NO_OSCILLATION 0 /* set to 1 to have exponential BETA change only */
#define MIDDLE_CONSTANT_PHASE 0 /* final phase in which temperature is constant */ #define INITIAL_CONSTANT_PHASE 500 /* initial phase in which temperature is constant */
#define FINAL_DECREASE_PHASE 0 /* final phase in which temperature decreases */ #define MIDDLE_CONSTANT_PHASE 0 /* middle phase in which temperature is constant */
#define FINAL_CONSTANT_PHASE -1 /* final phase in which temperature is constant */ #define FINAL_DECREASE_PHASE 1 /* final phase in which temperature decreases */
#define FINAL_CONSTANT_PHASE 200 /* final phase in which temperature is constant */
#define DECREASE_CONTAINER_SIZE 0 /* set to 1 to decrease size of container */ #define DECREASE_CONTAINER_SIZE 0 /* set to 1 to decrease size of container */
#define SMOOTH_CONTAINER_DECREASE 1 /* set to 1 to decrease size smoothly at each simulation step */ #define SMOOTH_CONTAINER_DECREASE 1 /* set to 1 to decrease size smoothly at each simulation step */
@ -427,23 +437,25 @@
#define MASS_PART_BOTTOM 10000.0 /* mass of 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 NPART_BOTTOM 100 /* number of particles at the bottom */
#define ADD_PARTICLES 1 /* set to 1 to add particles */ #define ADD_PARTICLES 0 /* set to 1 to add particles */
#define ADD_REGION 0 /* shape of add regions, cf ADD_* in global_ljones */ #define ADD_REGION 1 /* shape of add regions, cf ADD_* in global_ljones */
#define ADD_TIME 0 /* time at which to add first particle */ #define ADD_TIME 0 /* time at which to add first particle */
#define ADD_PERIOD 25 /* time interval between adding further particles */ #define ADD_PERIOD 6 /* time interval between adding further particles */
#define ADD_TYPE 1 /* type of added particles */
#define N_ADD_PARTICLES 1 /* number of particles to add */ #define N_ADD_PARTICLES 1 /* number of particles to add */
#define FINAL_NOADD_PERIOD 1000 /* final period where no particles are added */ #define FINAL_NOADD_PERIOD 1800 /* final period where no particles are added */
#define SAFETY_FACTOR 4.0 /* no particles are added at distance less than MU*SAFETY_FACTOR of other particles */ #define SAFETY_FACTOR 10.0 /* no particles are added at distance less than MU*SAFETY_FACTOR of other particles */
#define ADD_ALTERNATE_CHARGE 1 /* set to 1 to randomly select sign of added charge */ #define ADD_ALTERNATE_CHARGE 1 /* set to 1 to randomly select sign of added charge */
#define TIME_DEPENDENT_ADD_CHARGE 0 /* set to 1 to have added charge depend on time */
#define ALTERNATE_CHARGE_PROPORTION 0.5 /* proportion of particles of opposite charge */ #define ALTERNATE_CHARGE_PROPORTION 0.5 /* proportion of particles of opposite charge */
#define TRACER_PARTICLE 1 /* set to 1 to have a tracer particle */ #define TRACER_PARTICLE 1 /* set to 1 to have a tracer particle */
#define N_TRACER_PARTICLES 500 /* number of tracer particles */ #define N_TRACER_PARTICLES 3000 /* number of tracer particles */
#define TRACER_STEPS 5 /* number of tracer steps recorded between images */ #define TRACER_STEPS 5 /* number of tracer steps recorded between images */
#define TRAJECTORY_LENGTH 5200 /* length of recorded trajectory */ #define TRAJECTORY_LENGTH 7000 /* length of recorded trajectory */
#define TRAJECTORY_DRAW_LENGTH 1000 /* length of drawn trajectory */ #define TRAJECTORY_DRAW_LENGTH 250 /* length of drawn trajectory */
#define TRACER_LUM_FACTOR 25.0 /* controls luminosity decrease of trajectories with time */ #define TRACER_LUM_FACTOR 100.0 /* controls luminosity decrease of trajectories with time */
#define TRACER_PARTICLE_MASS 2.0 /* relative mass of tracer particle */ #define TRACER_PARTICLE_MASS 1.0 /* relative mass of tracer particle */
#define TRAJECTORY_WIDTH 2 /* width of tracer particle trajectory */ #define TRAJECTORY_WIDTH 2 /* width of tracer particle trajectory */
#define TRACK_PARTICLE 0 /* set to 1 to track a given particle */ #define TRACK_PARTICLE 0 /* set to 1 to track a given particle */
@ -494,26 +506,28 @@
#define POSITION_DEP_X -0.625 /* threshold value for position-dependent type */ #define POSITION_DEP_X -0.625 /* threshold value for position-dependent type */
#define PRINT_ENTROPY 0 /* set to 1 to compute entropy */ #define PRINT_ENTROPY 0 /* set to 1 to compute entropy */
#define SPECIAL_IC 0 /* set to 1 for choosing specaial initial condition RD_INITIAL_COND */ #define SPECIAL_IC 1 /* set to 1 for choosing specaial initial condition RD_INITIAL_COND */
#define REACTION_DIFFUSION 0 /* set to 1 to simulate a chemical reaction (particles may change type) */ #define REACTION_DIFFUSION 1 /* set to 1 to simulate a chemical reaction (particles may change type) */
#define REACTION_MAX_TIME 100 /* time after which no reactions take place */ #define REACTION_MAX_TIME 100000 /* time after which no reactions take place */
#define RD_REACTION 231 /* type of reaction, see list in global_ljones.c */ #define RD_REACTION 11 /* type of reaction, see list in global_ljones.c */
#define RD_TYPES 2 /* number of types in reaction-diffusion equation */ #define RD_TYPES 5 /* number of types in reaction-diffusion equation */
#define RD_INITIAL_COND 51 /* initial condition of particles */ #define RD_PLOT_TYPES 4 /* number of types shown in graph */
#define REACTION_DIST 1.8 /* maximal distance for reaction to occur */ #define RD_INITIAL_COND 5 /* initial condition of particles */
#define REACTION_PROB 1.0 /* probability controlling reaction term */ #define REACTION_DIST 2.0 /* maximal distance for reaction to occur */
#define REACTION_PROB 0.5 /* probability controlling reaction term */
#define DISSOCIATION_PROB 0.0 /* probability controlling dissociation reaction */ #define DISSOCIATION_PROB 0.0 /* probability controlling dissociation reaction */
#define KILLING_PROB 0.0015 /* probability of enzymes being killed */ #define KILLING_PROB 0.0015 /* probability of enzymes being killed */
#define DELTAMAX 0.1 /* max orientation difference for pairing polygons */ #define DELTAMAX 0.1 /* max orientation difference for pairing polygons */
#define CENTER_COLLIDED_PARTICLES 1 /* set to 1 to recenter particles upon reaction (may interfere with thermostat) */ #define CENTER_COLLIDED_PARTICLES 1 /* set to 1 to recenter particles upon reaction (may interfere with thermostat) */
#define EXOTHERMIC 1 /* set to 1 to make reaction exo/endothermic */ #define EXOTHERMIC 1 /* set to 1 to make reaction exo/endothermic */
#define DELTA_EKIN -2.0e10 /* change of kinetic energy in reaction */ #define DELTA_EKIN -2.0e3 /* change of kinetic energy in reaction */
#define CORRECT_EQUILIBRIUM_POSITION 1 /* set to 1 to nudge particle dist towards eq dist */ #define CORRECT_EQUILIBRIUM_POSITION 1 /* set to 1 to nudge particle dist towards eq dist */
#define COLLISION_TIME 25 /* time during which collisions are shown */ #define COLLISION_TIME 25 /* time during which collisions are shown */
#define COLLISION_RADIUS 2.0 /* radius of discs showing collisions, in units of MU */ #define COLLISION_RADIUS 2.0 /* radius of discs showing collisions, in units of MU */
#define DELTAVMAX 200.0 /* maximal deltav allowed for pairing molecules */ #define DELTAVMAX 200.0 /* maximal deltav allowed for pairing molecules */
#define AGREGMAX 1 /* maximal number of partners for CHEM_AGGREGATION reaction */ #define AGREGMAX 4 /* maximal number of partners for CHEM_AGGREGATION reaction */
#define AGREG_DECOUPLE 12 /* minimal number of partners to decouple from thermostat */ #define AGREG_DECOUPLE 12 /* minimal number of partners to decouple from thermostat */
#define NEUTRALIZE_REACTING_PARTICLES 0 /* set to 1 for reacting particles to become neutral */
#define CLUSTER_PARTICLES 0 /* set to 1 for particles to form rigid clusters */ #define CLUSTER_PARTICLES 0 /* set to 1 for particles to form rigid clusters */
#define CLUSTER_MAXSIZE 2 /* max size of clusters */ #define CLUSTER_MAXSIZE 2 /* max size of clusters */
#define SMALL_CLUSTER_MAXSIZE 2 /* size limitation on smaller cluster */ #define SMALL_CLUSTER_MAXSIZE 2 /* size limitation on smaller cluster */
@ -526,7 +540,7 @@
#define MU_RATIO 0.666666667 /* ratio by which to increase radius */ #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 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 PLOT_PARTICLE_NUMBER 1 /* set to 1 to make of plot of particle number over time */
#define PARTICLE_NB_PLOT_FACTOR 1.0 /* expected final number of particles over initial number */ #define PARTICLE_NB_PLOT_FACTOR 1.0 /* expected final number of particles over initial number */
#define PRINT_LEFT 0 /* set to 1 to print certain parameters at the top left instead of right */ #define 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_SPEEDS 0 /* set to 1 to add a plot of obstacle speeds (e.g. for rockets) */
@ -545,10 +559,12 @@
#define WALL_TIME 0 /* time during which to keep wall */ #define WALL_TIME 0 /* time during which to keep wall */
#define CHANGE_TYPES 0 /* set to 1 to change type proportion in course of simulation */ #define CHANGE_TYPES 0 /* set to 1 to change type proportion in course of simulation */
#define PROP_MIN 0.1 /* min proportion of type 1 particles */ #define PROP_MIN 0.0 /* min proportion of type 1 particles */
#define PROP_MAX 0.9 /* max proportion of type 1 particles */ #define PROP_MAX 1.0 /* max proportion of type 1 particles */
#define PROP_TINITIAL 250 /* initial time without change */
#define PROP_TFINAL 250 /* final time without change */
#define PAIR_PARTICLES 1 /* set to 1 to form particles pairs */ #define PAIR_PARTICLES 0 /* set to 1 to form particles pairs */
#define RANDOMIZE_ANGLE 0 /* set to 1 for random orientation */ #define RANDOMIZE_ANGLE 0 /* set to 1 for random orientation */
#define DEACIVATE_CLOSE_PAIRS 0 /* set to 1 to test for closeness to other particles */ #define DEACIVATE_CLOSE_PAIRS 0 /* set to 1 to test for closeness to other particles */
#define PAIR_SAFETY_FACTOR 1.2 /* distance to deactivate divided by sum of radii */ #define PAIR_SAFETY_FACTOR 1.2 /* distance to deactivate divided by sum of radii */
@ -564,7 +580,7 @@
#define PARTNER_ANGLE 104.45 /* angle (in degrees) between ions for POLY_WATER case */ #define PARTNER_ANGLE 104.45 /* angle (in degrees) between ions for POLY_WATER case */
#define PAIR_DRATIO 1.0 /* ratio between equilibrium distance and radius (default: 1.0) */ #define PAIR_DRATIO 1.0 /* ratio between equilibrium distance and radius (default: 1.0) */
#define MU_C 0.035 /* radius of partner particle */ #define MU_C 0.035 /* radius of partner particle */
#define PARTICLE_MASS_C 2.0 /* mass or partner particle */ #define PARTICLE_MASS_C 1.0 /* mass or partner particle */
#define CHARGE_C 1.0 /* charge of partner particle */ #define CHARGE_C 1.0 /* charge of partner particle */
#define CLUSTER_COLOR_FACTOR 40 /* factor for initialization of cluster colors */ #define CLUSTER_COLOR_FACTOR 40 /* factor for initialization of cluster colors */
#define ALTERNATE_POLY_CHARGE 1 /* set to 1 for alternating charges in molecule */ #define ALTERNATE_POLY_CHARGE 1 /* set to 1 for alternating charges in molecule */
@ -576,7 +592,7 @@
#define NARMS_B 1 /* number of "arms" for certain paring types */ #define NARMS_B 1 /* number of "arms" for certain paring types */
#define PAIRING_TYPE_B 81 /* type of pairing, see POLY_ in global_ljones.c */ #define PAIRING_TYPE_B 81 /* type of pairing, see POLY_ in global_ljones.c */
#define MU_D 0.035 /* radius of partner particle */ #define MU_D 0.035 /* radius of partner particle */
#define PARTICLE_MASS_D 2.0 /* mass or partner particle */ #define PARTICLE_MASS_D 1.0 /* mass or partner particle */
#define CHARGE_D 1.0 /* charge of partner particle */ #define CHARGE_D 1.0 /* charge of partner particle */
#define NXMAZE 16 /* width of maze */ #define NXMAZE 16 /* width of maze */
@ -587,12 +603,12 @@
#define MAZE_WIDTH 0.015 /* width of maze walls */ #define MAZE_WIDTH 0.015 /* width of maze walls */
#define FLOOR_FORCE 1 /* set to 1 to limit force on particle to FMAX */ #define FLOOR_FORCE 1 /* set to 1 to limit force on particle to FMAX */
#define FMAX 1.0e9 /* maximal force */ #define FMAX 1.0e8 /* maximal force */
#define FLOOR_OMEGA 0 /* set to 1 to limit particle momentum to PMAX */ #define FLOOR_OMEGA 0 /* set to 1 to limit particle momentum to PMAX */
#define PMAX 1000.0 /* maximal force */ #define PMAX 1000.0 /* maximal force */
#define HASHX 40 /* size of hashgrid in x direction */ #define HASHX 60 /* size of hashgrid in x direction */
#define HASHY 20 /* 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 HASHMAX 100 /* maximal number of particles per hashgrid cell */
#define HASHGRID_PADDING 0.1 /* padding of hashgrid outside simulation window */ #define HASHGRID_PADDING 0.1 /* padding of hashgrid outside simulation window */
@ -603,6 +619,29 @@
#define LIMIT_ENERGY 0 /* set to 1 to limit energy, when there is no thermostat */ #define LIMIT_ENERGY 0 /* set to 1 to limit energy, when there is no thermostat */
/* constants related to evolution on a sphere */
#define SPHERE 1 /* set to 1 to compute evolution in spherical geometry */
#define SIN_THETA_REG 0.05 /* regularization of sin(theta) for motion on sphere */
#define POLAR_PADDING 0.1 /* region around poles that belong to the same hashcell */
#define DRAW_SPHERE 1 /* set to 1 to draw 3D sphere */
#define NX_SPHERE 3000
#define NY_SPHERE 1500 /* number of points on sphere */
#define Z_SCALING_FACTOR 0.75 /* overall scaling factor of z axis for REP_PROJ_3D representation */
#define XY_SCALING_FACTOR 1.9 /* 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.35 /* limit on cosine of normal to shown facets */
#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 */
#define MAX_LATITUDE 45.0 /* maximal latitude for viewpoint trajectory VP_ORBIT2 */
double light[3] = {0.40824829, -0.816496581, 0.40824829}; /* vector of "light" direction for P_3D_ANGLE color scheme */
double observer[3] = {3.0, -3.0, -2.5}; /* location of observer for REP_PROJ_3D representation */
#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 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 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 TWO_OBSTACLES ((SEGMENT_PATTERN == S_TWO_CIRCLES_EXT)||(SEGMENT_PATTERN == S_TWO_ROCKETS))
@ -633,12 +672,14 @@ double vysegments[2] = {SEGMENTS_VY0, SEGMENTS_VY0}; /* vy coordinate of s
int thermostat_on = 1; /* thermostat switch used when VARY_THERMOSTAT is on */ int thermostat_on = 1; /* thermostat switch used when VARY_THERMOSTAT is on */
double cosangle[NSEG+1]; double cosangle[NSEG+1];
double sinangle[NSEG+1]; /* precomputed trig functions of angles to draw circles faster */ double sinangle[NSEG+1]; /* precomputed trig functions of angles to draw circles faster */
short int reset_view = 0; /* for moving observer in 3D plot */
#define THERMOSTAT_ON ((THERMOSTAT)&&((!VARY_THERMOSTAT)||(thermostat_on))) #define THERMOSTAT_ON ((THERMOSTAT)&&((!VARY_THERMOSTAT)||(thermostat_on)))
#include "global_ljones.c" #include "global_ljones.c"
#include "sub_maze.c" #include "sub_maze.c"
#include "sub_hashgrid.c" #include "sub_hashgrid.c"
#include "sub_lj_sphere.c"
#include "sub_lj.c" #include "sub_lj.c"
FILE *lj_time_series, *lj_final_position; FILE *lj_time_series, *lj_final_position;
@ -688,16 +729,17 @@ double bfield_schedule(int i)
double temperature_schedule(int i) double temperature_schedule(int i)
{ {
static double bexponent, omega, bexp2, factor2, logf, ac, bc; static double bexponent, omega, bexp2, factor2, logf, ac, bc;
static int first = 1, t1, t2, t3; static int first = 1, t1, t2, t3, t4;
double beta, t; double beta, t;
if (first) if (first)
{ {
t1 = NSTEPS - MIDDLE_CONSTANT_PHASE - FINAL_DECREASE_PHASE - FINAL_CONSTANT_PHASE; t1 =INITIAL_CONSTANT_PHASE;
t2 = NSTEPS - FINAL_DECREASE_PHASE - FINAL_CONSTANT_PHASE; t2 = NSTEPS - MIDDLE_CONSTANT_PHASE - FINAL_DECREASE_PHASE - FINAL_CONSTANT_PHASE;
t3 = NSTEPS - FINAL_CONSTANT_PHASE; t3 = NSTEPS - FINAL_DECREASE_PHASE - FINAL_CONSTANT_PHASE;
bexponent = log(BETA_FACTOR)/(double)(t1); t4 = NSTEPS - FINAL_CONSTANT_PHASE;
omega = N_TOSCILLATIONS*DPI/(double)(t1); bexponent = log(BETA_FACTOR)/(double)(t2-t1);
omega = N_TOSCILLATIONS*DPI/(double)(t2-t1);
logf = log(BETA_FACTOR); logf = log(BETA_FACTOR);
switch (BETA_SCHEDULE) switch (BETA_SCHEDULE)
@ -709,7 +751,7 @@ double temperature_schedule(int i)
} }
case (TS_CYCLING): case (TS_CYCLING):
{ {
factor2 = BETA_FACTOR*2.0/(1.0 + cos(N_TOSCILLATIONS*DPI)); factor2 = BETA_FACTOR*2.0/(1.0 + 1.0e-10 + cos(N_TOSCILLATIONS*DPI));
break; break;
} }
case (TS_PERIODIC): case (TS_PERIODIC):
@ -751,75 +793,78 @@ double temperature_schedule(int i)
factor2 = 1.0/BETA_FACTOR - 1.0; factor2 = 1.0/BETA_FACTOR - 1.0;
break; break;
} }
} }
bexp2 = -log(factor2)/(double)(FINAL_DECREASE_PHASE); bexp2 = -log(factor2)/(double)(FINAL_DECREASE_PHASE);
first = 0; first = 0;
printf("bexp2 = %.3lg\n", bexp2);
printf("factor2 = %.3lg\n", factor2);
// exit(0);
// printf("t1 = %i, factor2 = %.3lg\n", t1, factor2); // printf("t1 = %i, factor2 = %.3lg\n", t1, factor2);
} }
if (i < INITIAL_TIME) beta = BETA;
else if (i < INITIAL_TIME + t1) //
if (i < INITIAL_TIME + t1) beta = BETA;
else if (i < INITIAL_TIME + t2)
{ {
switch (BETA_SCHEDULE) switch (BETA_SCHEDULE)
{ {
case (TS_EXPONENTIAL): case (TS_EXPONENTIAL):
{ {
beta = BETA*exp(bexponent*(double)(i - INITIAL_TIME)); beta = BETA*exp(bexponent*(double)(i - INITIAL_TIME - t1));
break; break;
} }
case (TS_CYCLING): case (TS_CYCLING):
{ {
beta = BETA*exp(bexponent*(double)(i - INITIAL_TIME)); beta = BETA*exp(bexponent*(double)(i - INITIAL_TIME - t1));
beta = beta*2.0/(1.0 + cos(omega*(double)(i - INITIAL_TIME))); beta = beta*2.0/(1.0 + 1.0e-10 + cos(omega*(double)(i - INITIAL_TIME - t1)));
break; break;
} }
case (TS_PERIODIC): case (TS_PERIODIC):
{ {
beta = BETA*exp(logf*sin(omega*(double)(i - INITIAL_TIME))); beta = BETA*exp(logf*sin(omega*(double)(i - INITIAL_TIME - t1)));
break; break;
} }
case (TS_LINEAR): case (TS_LINEAR):
{ {
beta = BETA/(1.0 + (1.0/BETA_FACTOR - 1.0)*(double)(i - INITIAL_TIME)/(double)(t1)); beta = BETA/(1.0 + (1.0/BETA_FACTOR - 1.0)*(double)(i - INITIAL_TIME - t1)/(double)(t2-t1));
// printf("i = %i, beta = %.3lg\n", i, beta); // printf("i = %i, beta = %.3lg\n", i, beta);
break; break;
} }
case (TS_COSINE): case (TS_COSINE):
{ {
beta = ac/(1.0 + bc*cos(omega*(double)(i - INITIAL_TIME))); beta = ac/(1.0 + bc*cos(omega*(double)(i - INITIAL_TIME - t1)));
printf("i = %i, beta = %.3lg\n", i, beta); printf("i = %i, beta = %.3lg\n", i, beta);
break; break;
} }
case (TS_EXPCOS): case (TS_EXPCOS):
{ {
beta = BETA*exp(bc*(-1.0 + cos(omega*(double)(i - INITIAL_TIME)))); beta = BETA*exp(bc*(-1.0 + cos(omega*(double)(i - INITIAL_TIME - t1))));
// printf("i = %i, beta = %.3lg\n", i, beta); // printf("i = %i, beta = %.3lg\n", i, beta);
break; break;
} }
case (TS_ASYM_EXPCOS): case (TS_ASYM_EXPCOS):
{ {
t = (double)(i - INITIAL_TIME)/(double)(t1); t = (double)(i - INITIAL_TIME - t1)/(double)(t2 - t1);
beta = BETA*exp(bc*(-1.0 + cos(N_TOSCILLATIONS*DPI*(t - 0.5*t*(1.0-t))))); beta = BETA*exp(bc*(-1.0 + cos(N_TOSCILLATIONS*DPI*(t - 0.5*t*(1.0-t)))));
break; break;
} }
case (TS_ATAN): case (TS_ATAN):
{ {
beta = BETA/(1.0 + factor2*atan(2.0*(double)(i - INITIAL_TIME)/(double)(t1))); beta = BETA/(1.0 + factor2*atan(2.0*(double)(i - INITIAL_TIME - t1)/(double)(t2 - t1)));
break; break;
} }
case (TS_TANH): case (TS_TANH):
{ {
beta = BETA/(1.0 + factor2*tanh(TS_SLOPE*(double)(i - INITIAL_TIME)/(double)(t1))); beta = BETA/(1.0 + factor2*tanh(TS_SLOPE*(double)(i - INITIAL_TIME - t1)/(double)(t2 - t1)));
break; break;
} }
} }
} }
else if (i < INITIAL_TIME + t2) beta = BETA*factor2; // else if ((i < INITIAL_TIME + t3)&&(t3 > t2)) beta = BETA*factor2;
else if (i < INITIAL_TIME + t3) // else if ((i < INITIAL_TIME + t4)&&(t4 > t3))
{ // beta = BETA*exp(bexp2*(double)(i - INITIAL_TIME - t3) + 1.0e-10);
beta = BETA*exp(bexp2*(double)(i - INITIAL_TIME - t3)); else beta = BETA*1.0e10;
} // else beta = BETA;
else beta = BETA;
printf("beta = %.3lg\n", beta); printf("beta = %.3lg\n", beta);
return(beta); return(beta);
} }
@ -1035,6 +1080,7 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NM
double beta, int *nactive, int *nsuccess, int *nmove, int *ncoupled, int initial_phase) double beta, int *nactive, int *nsuccess, int *nmove, int *ncoupled, int initial_phase)
{ {
double a, totalenergy = 0.0, damping, damping1, damping_rot1, direction, dmean, dratio; double a, totalenergy = 0.0, damping, damping1, damping_rot1, direction, dmean, dratio;
double ctheta, stheta, stheta_reg, ffx, ffy;
static double b = 0.25*SIGMA*SIGMA*DT_PARTICLE/MU_XI, xi = 0.0; static double b = 0.25*SIGMA*SIGMA*DT_PARTICLE/MU_XI, xi = 0.0;
int j, move, ncoup; int j, move, ncoup;
@ -1047,15 +1093,40 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NM
for (j=0; j<ncircles; j++) if (particle[j].active) for (j=0; j<ncircles; j++) if (particle[j].active)
{ {
// printf("Particle %i\n", j); // printf("Particle %i\n", j);
particle[j].vx = px[j] + 0.5*DT_PARTICLE*particle[j].fx; if (SPHERE) /* add inertial terms for equation on a sphere */
particle[j].vy = py[j] + 0.5*DT_PARTICLE*particle[j].fy; {
particle[j].omega = pangle[j] + 0.5*DT_PARTICLE*particle[j].torque; ctheta = cos(particle[j].yc);
stheta = sin(particle[j].yc);
stheta_reg = stheta + SIN_THETA_REG;
ffx = -2.0*ctheta*px[j]*py[j];
ffy = stheta*ctheta*px[j]*px[j];
particle[j].vx = px[j] + 0.5*DT_PARTICLE*(particle[j].fx + ffx)/stheta_reg;
particle[j].vy = py[j] + 0.5*DT_PARTICLE*(particle[j].fy + ffy);
particle[j].omega = pangle[j] + 0.5*DT_PARTICLE*particle[j].torque;
// particle[j].omega = pangle[j] + 0.5*DT_PARTICLE*particle[j].torque/stheta_reg;
/* TODO: check/fix angular evolution */
ffx = -2.0*ctheta*particle[j].vx*particle[j].vy;
ffy = stheta*ctheta*particle[j].vx*particle[j].vx;
px[j] = particle[j].vx + 0.5*DT_PARTICLE*(particle[j].fx + ffx)/stheta_reg;
py[j] = particle[j].vy + 0.5*DT_PARTICLE*(particle[j].fy + ffy);
pangle[j] = particle[j].omega + 0.5*DT_PARTICLE*particle[j].torque;
// pangle[j] = particle[j].omega + 0.5*DT_PARTICLE*particle[j].torque/stheta_reg;
particle[j].energy = (stheta*stheta*px[j]*px[j] + py[j]*py[j])*particle[j].mass_inv;
}
else
{
particle[j].vx = px[j] + 0.5*DT_PARTICLE*particle[j].fx;
particle[j].vy = py[j] + 0.5*DT_PARTICLE*particle[j].fy;
particle[j].omega = pangle[j] + 0.5*DT_PARTICLE*particle[j].torque;
px[j] = particle[j].vx + 0.5*DT_PARTICLE*particle[j].fx; px[j] = particle[j].vx + 0.5*DT_PARTICLE*particle[j].fx;
py[j] = particle[j].vy + 0.5*DT_PARTICLE*particle[j].fy; py[j] = particle[j].vy + 0.5*DT_PARTICLE*particle[j].fy;
pangle[j] = particle[j].omega + 0.5*DT_PARTICLE*particle[j].torque; 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; particle[j].energy = (px[j]*px[j] + py[j]*py[j])*particle[j].mass_inv;
}
if (COMPUTE_EMEAN) if (COMPUTE_EMEAN)
particle[j].emean = ERATIO*particle[j].emean + (1.0-ERATIO)*particle[j].energy; particle[j].emean = ERATIO*particle[j].emean + (1.0-ERATIO)*particle[j].energy;
@ -1170,7 +1241,7 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NM
if (particle[j].yc > BCYMAX) particle[j].yc += BCYMIN - BCYMAX; if (particle[j].yc > BCYMAX) particle[j].yc += BCYMIN - BCYMAX;
else if (particle[j].yc < BCYMIN) particle[j].yc += BCYMAX - BCYMIN; else if (particle[j].yc < BCYMIN) particle[j].yc += BCYMAX - BCYMIN;
} }
else if (!NO_WRAP_BC) else if ((!NO_WRAP_BC)||(BOUNDARY_COND == BC_SPHERE))
{ {
move += wrap_particle(&particle[j], &px[j], &py[j]); move += wrap_particle(&particle[j], &px[j], &py[j]);
} }
@ -1757,6 +1828,7 @@ void animation()
t_belt *conveyor_belt; t_belt *conveyor_belt;
t_otriangle *otriangle; t_otriangle *otriangle;
t_ofacet *ofacet; t_ofacet *ofacet;
t_lj_sphere *wsphere;
char message[100]; char message[100];
ratioc = 1.0 - ratio; ratioc = 1.0 - ratio;
@ -1772,6 +1844,13 @@ void animation()
particle = (t_particle *)malloc(NMAXCIRCLES*sizeof(t_particle)); /* particles */ particle = (t_particle *)malloc(NMAXCIRCLES*sizeof(t_particle)); /* particles */
/* sphere data */
if (DRAW_SPHERE)
{
wsphere = (t_lj_sphere *)malloc(NX_SPHERE*NY_SPHERE*sizeof(t_lj_sphere));
init_lj_sphere(wsphere);
}
if (CLUSTER_PARTICLES) if (CLUSTER_PARTICLES)
cluster = (t_cluster *)malloc(NMAXCIRCLES*sizeof(t_cluster)); /* particle clusters */ cluster = (t_cluster *)malloc(NMAXCIRCLES*sizeof(t_cluster)); /* particle clusters */
@ -1920,6 +1999,12 @@ void animation()
printf("Termostat: %i\n", thermostat_on); printf("Termostat: %i\n", thermostat_on);
} }
if (ROTATE_VIEW)
{
viewpoint_schedule(i - INITIAL_TIME);
reset_view = 1;
}
/* deactivate some segments */ /* deactivate some segments */
if ((ADD_FIXED_SEGMENTS)&&(DEACTIVATE_SEGMENT)&&(i == INITIAL_TIME + SEGMENT_DEACTIVATION_TIME + 1)) if ((ADD_FIXED_SEGMENTS)&&(DEACTIVATE_SEGMENT)&&(i == INITIAL_TIME + SEGMENT_DEACTIVATION_TIME + 1))
@ -2411,7 +2496,9 @@ void animation()
/* case of varying type proportion */ /* case of varying type proportion */
if (CHANGE_TYPES) if (CHANGE_TYPES)
{ {
t = (double)i/(double)NSTEPS; if (i < PROP_TINITIAL) t = 0.0;
else if (i > NSTEPS - PROP_TFINAL) t = 1.0;
else t = (double)(i - PROP_TINITIAL)/(double)(NSTEPS - PROP_TINITIAL - PROP_TFINAL);
params.prop = PROP_MIN*(1.0-t) + PROP_MAX*t; params.prop = PROP_MIN*(1.0-t) + PROP_MAX*t;
change_type_proportion(particle, params.prop); change_type_proportion(particle, params.prop);
} }
@ -2427,11 +2514,12 @@ void animation()
if ((REACTION_DIFFUSION)&&((RD_REACTION == CHEM_DNA_ENZYME)||(RD_REACTION == CHEM_DNA_ENZYME_REPAIR))) if ((REACTION_DIFFUSION)&&((RD_REACTION == CHEM_DNA_ENZYME)||(RD_REACTION == CHEM_DNA_ENZYME_REPAIR)))
{ {
for (k=0; k<N_ADD_PARTICLES; k++) for (k=0; k<N_ADD_PARTICLES; k++)
nadd_particle = add_particles(particle, px, py, nadd_particle, 7, molecule, tracer_n); nadd_particle = add_particles(particle, px, py, nadd_particle, 7, molecule, tracer_n, i);
} }
else for (k=0; k<N_ADD_PARTICLES; k++) else for (k=0; k<N_ADD_PARTICLES; k++)
nadd_particle = add_particles(particle, px, py, nadd_particle, 3, molecule, tracer_n); nadd_particle = add_particles(particle, px, py, nadd_particle, ADD_TYPE, molecule, tracer_n, i);
// params.nactive = nadd_particle; // params.nactive = nadd_particle;
if (PAIR_PARTICLES) init_molecule_data(particle, molecule);
params.nactive = 0; params.nactive = 0;
for (j=0; j<ncircles; j++) if (particle[j].active) for (j=0; j<ncircles; j++) if (particle[j].active)
{ {
@ -2480,7 +2568,7 @@ void animation()
draw_frame(i, PLOT, BG_COLOR, ncollisions, traj_position, traj_length, draw_frame(i, PLOT, BG_COLOR, ncollisions, traj_position, traj_length,
wall, pressure, pleft, pright, currents, particle_numbers, 1, params, particle, cluster, wall, pressure, pleft, pright, currents, particle_numbers, 1, params, particle, cluster,
collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n, otriangle, ofacet); collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n, otriangle, ofacet, wsphere);
if (!((NO_EXTRA_BUFFER_SWAP)&&(MOVIE))) glutSwapBuffers(); if (!((NO_EXTRA_BUFFER_SWAP)&&(MOVIE))) glutSwapBuffers();
@ -2512,7 +2600,7 @@ void animation()
{ {
draw_frame(i, PLOT_B, BG_COLOR_B, ncollisions, traj_position, traj_length, draw_frame(i, PLOT_B, BG_COLOR_B, ncollisions, traj_position, traj_length,
wall, pressure, pleft, pright, currents, particle_numbers, 0, params, particle, cluster, wall, pressure, pleft, pright, currents, particle_numbers, 0, params, particle, cluster,
collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n, otriangle, ofacet); collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n, otriangle, ofacet, wsphere);
glutSwapBuffers(); glutSwapBuffers();
save_frame_lj_counter(NSTEPS + MID_FRAMES + 1 + counter); save_frame_lj_counter(NSTEPS + MID_FRAMES + 1 + counter);
counter++; counter++;
@ -2553,7 +2641,7 @@ void animation()
blank(); blank();
draw_frame(NSTEPS, PLOT, BG_COLOR, ncollisions, traj_position, traj_length, draw_frame(NSTEPS, PLOT, BG_COLOR, ncollisions, traj_position, traj_length,
wall, pressure, pleft, pright, currents, particle_numbers, 0, params, particle, cluster, wall, pressure, pleft, pright, currents, particle_numbers, 0, params, particle, cluster,
collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n, otriangle, ofacet); collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n, otriangle, ofacet, wsphere);
} }
if (DOUBLE_MOVIE) for (i=0; i<MID_FRAMES; i++) if (DOUBLE_MOVIE) for (i=0; i<MID_FRAMES; i++)
{ {
@ -2565,7 +2653,7 @@ void animation()
{ {
draw_frame(NSTEPS, PLOT_B, BG_COLOR_B, ncollisions, traj_position, traj_length, draw_frame(NSTEPS, PLOT_B, BG_COLOR_B, ncollisions, traj_position, traj_length,
wall, pressure, pleft, pright, currents, particle_numbers, 0, params, particle, cluster, wall, pressure, pleft, pright, currents, particle_numbers, 0, params, particle, cluster,
collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n, otriangle, ofacet); collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n, otriangle, ofacet, wsphere);
if (!((NO_EXTRA_BUFFER_SWAP)&&(MOVIE))) glutSwapBuffers(); if (!((NO_EXTRA_BUFFER_SWAP)&&(MOVIE))) glutSwapBuffers();
} }
if ((TIME_LAPSE)&&(!DOUBLE_MOVIE)) if ((TIME_LAPSE)&&(!DOUBLE_MOVIE))
@ -2651,6 +2739,8 @@ void animation()
fclose(lj_final_position); fclose(lj_final_position);
} }
if (DRAW_SPHERE) free(wsphere);
fclose(lj_log); fclose(lj_log);
} }
@ -2691,7 +2781,8 @@ int main(int argc, char** argv)
glutInitWindowSize(WINWIDTH,WINHEIGHT); glutInitWindowSize(WINWIDTH,WINHEIGHT);
glutCreateWindow("Particles with Lennard-Jones interaction in a planar domain"); glutCreateWindow("Particles with Lennard-Jones interaction in a planar domain");
init(); if (DRAW_SPHERE) init_3d();
else init();
glutDisplayFunc(display); glutDisplayFunc(display);

1
mangrove.c
View File

@ -260,6 +260,7 @@
#define WAVE_PACKET_SOURCE_TYPE 1 /* type of wave packet sources */ #define WAVE_PACKET_SOURCE_TYPE 1 /* type of wave packet sources */
#define N_WAVE_PACKETS 15 /* number of wave packets */ #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 OSCIL_LEFT_YSHIFT 0.0 /* y-dependence of left oscillation (for non-horizontal waves) */
#define OSCIL_YMID -0.9 /* defines oscilling beam midpoint */
#define DRAW_WAVE_PROFILE 0 /* set to 1 to draw a profile of the wave */ #define DRAW_WAVE_PROFILE 0 /* set to 1 to draw a profile of the wave */
#define OSCILLATING_SOURCE_PERIOD 20 /* period of oscillating source */ #define OSCILLATING_SOURCE_PERIOD 20 /* period of oscillating source */
#define MU_B 1.0 /* parameter controlling the dimensions of domain */ #define MU_B 1.0 /* parameter controlling the dimensions of domain */

149
rde.c
View File

@ -46,31 +46,27 @@
/* General geometrical parameters */ /* General geometrical parameters */
// #define WINWIDTH 1920 /* window width */ #define WINWIDTH 1920 /* window width */
#define WINWIDTH 1150 /* window width */
#define WINHEIGHT 1150 /* window height */ #define WINHEIGHT 1150 /* window height */
// #define NX 720 /* number of grid points on x axis */ #define NX 870 /* number of grid points on x axis */
#define NX 432 /* number of grid points on x axis */ #define NY 520 /* number of grid points on y axis */
#define NY 432 /* number of grid points on y axis */
#define HRES 1 /* factor for high resolution plots */ #define HRES 1 /* factor for high resolution plots */
// #define XMIN -2.0 #define XMIN -2.0
// #define XMAX 2.0 /* x interval */ #define XMAX 2.0 /* x interval */
#define XMIN -1.041666667
#define XMAX 1.041666667 /* x interval */
#define YMIN -1.041666667 #define YMIN -1.041666667
#define YMAX 1.041666667 /* y interval for 9/16 aspect ratio */ #define YMAX 1.041666667 /* y interval for 9/16 aspect ratio */
/* Choice of simulated equation */ /* Choice of simulated equation */
#define RDE_EQUATION 10 /* choice of reaction term, see list in global_3d.c */ #define RDE_EQUATION 11 /* choice of reaction term, see list in global_3d.c */
#define NFIELDS 2 /* number of fields in reaction-diffusion equation */ #define NFIELDS 3 /* number of fields in reaction-diffusion equation */
#define NLAPLACIANS 2 /* number of fields for which to compute Laplacian */ #define NLAPLACIANS 0 /* number of fields for which to compute Laplacian */
#define SPHERE 0 /* set to 1 to simulate equation on sphere */ #define SPHERE 1 /* set to 1 to simulate equation on sphere */
#define DPOLE 1 /* safety distance to poles */ #define DPOLE 1 /* safety distance to poles */
#define DSMOOTH 1 /* size of neighbourhood of poles that are smoothed */ #define DSMOOTH 1 /* size of neighbourhood of poles that are smoothed */
#define SMOOTHPOLE 0.00 /* smoothing coefficient at poles */ #define SMOOTHPOLE 0.01 /* smoothing coefficient at poles */
#define SMOOTHCOTPOLE 0.05 /* smoothing coefficient of cotangent at poles */ #define SMOOTHCOTPOLE 0.05 /* smoothing coefficient of cotangent at poles */
#define PHISHIFT 0.0 /* shift of phi in 2D plot (in degrees) */ #define PHISHIFT 0.0 /* shift of phi in 2D plot (in degrees) */
#define SMOOTHBLOCKS 0 /* set to 1 to use blocks of points near the poles */ #define SMOOTHBLOCKS 0 /* set to 1 to use blocks of points near the poles */
@ -145,8 +141,7 @@
/* Physical parameters of wave equation */ /* Physical parameters of wave equation */
#define DT 0.000003 #define DT 0.000002
// #define DT 0.00000025
#define VISCOSITY 0.02 #define VISCOSITY 0.02
#define POISSON_STIFFNESS 1.0 /* stiffness of Poisson equation solver for incompressible Euler */ #define POISSON_STIFFNESS 1.0 /* stiffness of Poisson equation solver for incompressible Euler */
@ -167,11 +162,12 @@
#define NU_KURAMOTO 0.0 /* viscosity in Kuramoto model */ #define NU_KURAMOTO 0.0 /* viscosity in Kuramoto model */
#define OMEGA_KURAMOTO 0.02 /* frequency in Kuramoto model */ #define OMEGA_KURAMOTO 0.02 /* frequency in Kuramoto model */
#define A_KS 0.15 /* coupling constant in Keller-Segel model */ #define A_KS 0.15 /* coupling constant in Keller-Segel model */
#define C_KS 3.0 /* coupling constant in Keller-Segel model */ #define C_KS 3.5 /* coupling constant in Keller-Segel model */
#define D_KSU 1.2 /* diffusion in Keller-Segel model */ #define D_KSU 1.0 /* organism diffusion in Keller-Segel model */
#define D_KSV 0.75 /* diffusion in Keller-Segel model */ #define D_KSV 1.2 /* nutrient diffusion in Keller-Segel model */
#define KS_RSCALE 0.001 /* scaling factor for reaction term of Keller-Segel model */ #define KS_RSCALE 0.07 /* scaling factor for reaction term of Keller-Segel model */
#define KS_UMAX 10.0 /* max value for u component of Keller-Segel model */ #define KS_UMAX 100.0 /* max value for u component of Keller-Segel model */
#define KS_CHIMAX_FACTOR 30.0 /* limiter on the value of chi/D_KSU */
#define V_MAZE 0.4 /* potential in walls of maze */ #define V_MAZE 0.4 /* potential in walls of maze */
#define BZQ 0.0008 /* parameter in BZ equation */ #define BZQ 0.0008 /* parameter in BZ equation */
#define BZF 1.2 /* parameter in BZ equation */ #define BZF 1.2 /* parameter in BZ equation */
@ -254,8 +250,8 @@
/* Parameters for length and speed of simulation */ /* Parameters for length and speed of simulation */
#define NSTEPS 1800 /* number of frames of movie */ #define NSTEPS 1400 /* number of frames of movie */
#define NVID 100 /* number of iterations between images displayed on screen */ #define NVID 12 /* number of iterations between images displayed on screen */
#define ACCELERATION_FACTOR 1.0 /* factor by which to increase NVID in course of simulation */ #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 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 MAX_DT 0.024 /* maximal value of integration step */
@ -274,14 +270,15 @@
/* Visualisation */ /* Visualisation */
#define PLOT_3D 1 /* controls whether plot is 2D or 3D */ #define PLOT_3D 1 /* controls whether plot is 2D or 3D */
#define PLOT_SPHERE 0 /* draws fields on a sphere */ #define PLOT_SPHERE 1 /* draws fields on a sphere */
#define SIMULATION_GRID 1 /* type of simulation grid, for certain equations */ #define SIMULATION_GRID 1 /* type of simulation grid, for certain equations */
#define OPTIMIZE_GRID 1 /* set to 1 to optimize grid by evolving grid points */
#define ROTATE_VIEW 1 /* 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 ROTATE_ANGLE 360.0 /* total angle of rotation during simulation */
#define SHADE_3D 1 /* set to 1 to change luminosity according to normal vector */ #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 */ #define SHADE_2D 0 /* set to 1 to change luminosity according to normal vector */
#define VIEWPOINT_TRAJ 0 /* type of viewpoint trajectory */ #define VIEWPOINT_TRAJ 1 /* type of viewpoint trajectory */
#define MAX_LATITUDE 45.0 /* maximal latitude for viewpoint trajectory VP_ORBIT2 */ #define MAX_LATITUDE 45.0 /* maximal latitude for viewpoint trajectory VP_ORBIT2 */
#define DRAW_PERIODICISED 0 /* set to 1 to repeat wave periodically in x and y directions */ #define DRAW_PERIODICISED 0 /* set to 1 to repeat wave periodically in x and y directions */
@ -340,7 +337,7 @@
/* Color schemes, see list in global_pdes.c */ /* Color schemes, see list in global_pdes.c */
#define COLOR_PALETTE 11 /* Color palette, see list in global_pdes.c */ #define COLOR_PALETTE 15 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 11 /* Color palette, see list in global_pdes.c */ #define COLOR_PALETTE_B 11 /* Color palette, see list in global_pdes.c */
#define BLACK 1 /* black background */ #define BLACK 1 /* black background */
@ -355,8 +352,8 @@
#define SCALE 0 /* set to 1 to adjust color scheme to variance of field */ #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 SLOPE 1.0 /* sensitivity of color on wave amplitude */
#define COLOR_RANGE 1.0 /* max range of color (default: 1.0) */ #define COLOR_RANGE 1.0 /* max range of color (default: 1.0) */
#define VSHIFT_AMPLITUDE -2.7 /* additional shift for wave amplitude */ #define VSHIFT_AMPLITUDE -0.75 /* additional shift for wave amplitude */
#define VSCALE_AMPLITUDE 3.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 ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */
#define CURL_SCALE 1.25 /* scaling factor for curl representation */ #define CURL_SCALE 1.25 /* 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 RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */
@ -398,7 +395,7 @@
#define MAZE_XSHIFT 0.0 /* horizontal shift of maze */ #define MAZE_XSHIFT 0.0 /* horizontal shift of maze */
#define MAZE_WIDTH 0.04 /* half width of maze walls */ #define MAZE_WIDTH 0.04 /* half width of maze walls */
#define DRAW_COLOR_SCHEME 0 /* set to 1 to plot the color scheme */ #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 2.5 /* scale of color scheme bar */
#define COLORBAR_RANGE_B 2.5 /* scale of color scheme bar for 2nd part */ #define COLORBAR_RANGE_B 2.5 /* scale of color scheme bar for 2nd part */
#define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ #define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */
@ -413,6 +410,7 @@
#define MANDEL_IOR_SCALE -0.05 /* parameter controlling dependence of IoR on Mandelbrot escape speed */ #define MANDEL_IOR_SCALE -0.05 /* parameter controlling dependence of IoR on Mandelbrot escape speed */
#define OMEGA 0.005 /* frequency of periodic excitation */ #define OMEGA 0.005 /* frequency of periodic excitation */
#define OSCIL_YMAX 0.2 /* defines oscillation range */ #define OSCIL_YMAX 0.2 /* defines oscillation range */
#define OSCIL_YMID -0.75 /* defines oscilling beam midpoint */
#define COURANT 0.08 /* Courant number */ #define COURANT 0.08 /* Courant number */
#define COURANTB 0.03 /* Courant number in medium B */ #define COURANTB 0.03 /* Courant number in medium B */
#define INITIAL_AMP 0.5 /* amplitude of initial condition */ #define INITIAL_AMP 0.5 /* amplitude of initial condition */
@ -480,12 +478,12 @@ int reset_view = 0; /* switch to reset 3D view parameters (for option RO
#define VENUS_NODATA_FACTOR 0.5 /* altitude to assign to DEM points without data (fraction of mean altitude) */ #define VENUS_NODATA_FACTOR 0.5 /* altitude to assign to DEM points without data (fraction of mean altitude) */
#define Z_SCALING_FACTOR 0.75 /* overall scaling factor of z axis for REP_PROJ_3D representation */ #define Z_SCALING_FACTOR 0.75 /* overall scaling factor of z axis for REP_PROJ_3D representation */
#define XY_SCALING_FACTOR 1.5 /* overall scaling factor for on-screen (x,y) coordinates after projection */ #define XY_SCALING_FACTOR 1.9 /* 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 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 XSHIFT_3D -0.0 /* overall x shift for REP_PROJ_3D representation */
#define YSHIFT_3D -0.5 /* overall y shift for REP_PROJ_3D representation */ #define YSHIFT_3D -0.0 /* 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 */ #define BORDER_PADDING 0 /* distance from boundary at which to plot points, to avoid boundary effects due to gradient */
#define DRAW_ARROW 0 /* set to 1 to draw arrow above sphere */ #define DRAW_ARROW 1 /* set to 1 to draw arrow above sphere */
#define ZMAX 3.0 /* maximal value of z coordinate */ #define ZMAX 3.0 /* maximal value of z coordinate */
#define ZMIN -3.0 /* maximal value of z coordinate */ #define ZMIN -3.0 /* maximal value of z coordinate */
@ -493,7 +491,7 @@ int reset_view = 0; /* switch to reset 3D view parameters (for option RO
#define RSCALE_B 1.00 /* experimental, additional radial scaling factor in ij_to_sphere */ #define RSCALE_B 1.00 /* experimental, additional radial scaling factor in ij_to_sphere */
#define RSHIFT 0.0 /* shift in radial component */ #define RSHIFT 0.0 /* shift in radial component */
#define RMAX 4.0 /* max value of radial component */ #define RMAX 4.0 /* max value of radial component */
#define RMIN 0.5 /* min value of radial component */ #define RMIN 0.0 /* min value of radial component */
#define COS_VISIBLE -0.35 /* limit on cosine of normal to shown facets */ #define COS_VISIBLE -0.35 /* limit on cosine of normal to shown facets */
/* For debugging purposes only */ /* For debugging purposes only */
@ -1251,10 +1249,11 @@ double gfield[2*NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY])
{ {
int i, j, k, iplus, iminus, jplus, jminus, ropening, w, n; int i, j, k, iplus, iminus, jplus, jminus, ropening, w, n;
double x, y, z, deltax, deltay, deltaz, rho, rhox, rhoy, pot, u, v, ux, uy, vx, vy, test = 0.0, dx, dy, xy[2], padding, a, eta, etax, etay, sum, phi0, chi; double x, y, z, deltax, deltay, deltaz, rho, rhox, rhoy, pot, u, v, ux, uy, vx, vy, test = 0.0, dx, dy, xy[2], padding, a, eta, etax, etay, sum, phi0, chi;
double *delta_phi[NLAPLACIANS], *nabla_phi, *nabla_psi, *nabla_omega, *delta_vorticity, *delta_pressure, *delta_p, *delta_u, *delta_v, *nabla_rho, *nabla_u, *nabla_v, *nabla_eta, *kuramoto_int, *keller_segel_drift, *keller_segel_div; double *delta_phi[NLAPLACIANS], *nabla_phi, *nabla_psi, *nabla_omega, *delta_vorticity, *delta_pressure, *delta_p, *delta_u, *delta_v, *nabla_rho, *nabla_u, *nabla_v, *nabla_eta, *kuramoto_int, *chi_u, *prod_gradients;
// double u_bc[NY], v_bc[NY]; // double u_bc[NY], v_bc[NY];
static double invsqr3 = 0.577350269; /* 1/sqrt(3) */ static double invsqr3 = 0.577350269; /* 1/sqrt(3) */
static int smooth = 0, y_channels, y_channels1, imin, imax, first = 1; static int smooth = 0, y_channels, y_channels1, imin, imax, first = 1;
int floor_chi, floor_u;
if (first) /* for D_MAZE_CHANNELS boundary conditions in Euler equation */ if (first) /* for D_MAZE_CHANNELS boundary conditions in Euler equation */
{ {
@ -1424,10 +1423,7 @@ double gfield[2*NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY])
case (E_KELLER_SEGEL): case (E_KELLER_SEGEL):
{ {
// keller_segel_drift = (double *)malloc(NX*NY*sizeof(double));
// keller_segel_div = (double *)malloc(NX*NY*sizeof(double));
nabla_psi = (double *)malloc(2*NX*NY*sizeof(double)); nabla_psi = (double *)malloc(2*NX*NY*sizeof(double));
// compute_gradient_euler_test(phi_in[1], nabla_psi, xy_in, wsphere);
compute_gradient_euler_plane(phi_in[1], nabla_psi, rde, 0.0); compute_gradient_euler_plane(phi_in[1], nabla_psi, rde, 0.0);
/* multiply gradient by chi */ /* multiply gradient by chi */
#pragma omp parallel for private(i) #pragma omp parallel for private(i)
@ -1435,15 +1431,51 @@ double gfield[2*NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY])
{ {
x = phi_in[0][i]; x = phi_in[0][i];
chi = C_KS*x*(1.0+x*x); chi = C_KS*x*(1.0+x*x);
// chi = C_KS*x/(1.0+x*x);
nabla_psi[i] *= chi; nabla_psi[i] *= chi;
nabla_psi[NX*NY + i] *= chi; nabla_psi[NX*NY + i] *= chi;
} }
/* compute divergence */ /* compute divergence */
// compute_divergence(nabla_psi, keller_segel_div, xy_in, wsphere);
compute_divergence_rde(nabla_psi, rde, xy_in, wsphere); compute_divergence_rde(nabla_psi, rde, xy_in, wsphere);
break; break;
} }
case (E_KELLER_SEGEL_SPHERE):
{
delta_u = (double *)malloc(NX*NY*sizeof(double));
delta_v = (double *)malloc(NX*NY*sizeof(double));
chi_u = (double *)malloc(NX*NY*sizeof(double));
prod_gradients = (double *)malloc(NX*NY*sizeof(double));
compute_laplacian_grid(phi_in[1], delta_u, wsphere, ngridpoints);
compute_laplacian_grid(phi_in[2], delta_v, wsphere, ngridpoints);
/* multiply gradient by chi */
floor_chi = 0;
#pragma omp parallel for private(i)
for (i=0; i<ngridpoints; i++)
{
x = phi_in[1][i];
chi_u[i] = C_KS*x*(1.0+x*x);
// chi_u[i] = C_KS*x/(1.0+x*x);
if (chi_u[i] > KS_CHIMAX_FACTOR*D_KSU)
{
// printf("chi_i[%i] = %.5lg\n", i, chi_u[i]);
chi_u[i] = KS_CHIMAX_FACTOR*D_KSU;
floor_chi++;
}
}
if (floor_chi > 0)
{
nfloorchi++;
// printf("Floored chi at %i points\n", floor_chi);
}
// printf("Floored chi %i times\n", nfloorchi);
/* compute product of gradients */
compute_grad_product_grid(phi_in[2], chi_u, prod_gradients, wsphere, ngridpoints);
break;
}
default: default:
{ {
/* do nothing */ /* do nothing */
@ -1507,10 +1539,8 @@ double gfield[2*NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY])
y = phi_in[1][i*NY+j]; y = phi_in[1][i*NY+j];
deltax = D_KSU*delta_phi[0][i*NY+j]; deltax = D_KSU*delta_phi[0][i*NY+j];
deltay = D_KSV*delta_phi[1][i*NY+j]; deltay = D_KSV*delta_phi[1][i*NY+j];
// phi_out[0][i*NY+j] = x + intstep*(deltax - keller_segel_div[i*NY+j] + KS_RSCALE*x*(1.0-x));
phi_out[0][i*NY+j] = x + intstep*(deltax - rde[i*NY+j].divergence + KS_RSCALE*x*(1.0-x)); phi_out[0][i*NY+j] = x + intstep*(deltax - rde[i*NY+j].divergence + KS_RSCALE*x*(1.0-x));
if (phi_out[0][i*NY+j] > KS_UMAX) phi_out[0][i*NY+j] = KS_UMAX; if (phi_out[0][i*NY+j] > KS_UMAX) phi_out[0][i*NY+j] = KS_UMAX;
// if (phi_out[0][i*NY+j] < -1.0) phi_out[0][i*NY+j] = -1.0;
phi_out[1][i*NY+j] = y + intstep*(deltay + KS_RSCALE*(x - A_KS*y)); phi_out[1][i*NY+j] = y + intstep*(deltay + KS_RSCALE*(x - A_KS*y));
break; break;
@ -1781,6 +1811,26 @@ double gfield[2*NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY])
/* TEST */ /* TEST */
// for (i=0; i<NX*NY; i++) rde[i].theta = phi_out[0][i]; // for (i=0; i<NX*NY; i++) rde[i].theta = phi_out[0][i];
} }
else if (RDE_EQUATION == E_KELLER_SEGEL_SPHERE)
{
floor_u = 0;
#pragma omp parallel for private(i,x)
for (i=0; i<ngridpoints; i++)
{
x = phi_in[1][i];
phi_out[1][i] = x + intstep*(D_KSU*delta_u[i] - chi_u[i]*delta_v[i] - prod_gradients[i] + KS_RSCALE*x*(1.0-x));
if (phi_out[1][i] > KS_UMAX)
{
phi_out[1][i] = KS_UMAX;
floor_u = 1;
}
phi_out[2][i] = phi_in[2][i] + intstep*(D_KSV*delta_v[i] + KS_RSCALE*(x - A_KS*y));
}
if (floor_u) printf("u floored to KS_UMAX\n");
/* convert field 1 on simulation grid to field 0 on longitude_latitude grid */
convert_fields_from_grids(phi_out[1], phi_out[0], wsphere);
}
/* in-flow/out-flow b.c. for incompressible Euler equation */ /* in-flow/out-flow b.c. for incompressible Euler equation */
if (((RDE_EQUATION == E_EULER_INCOMP)||(RDE_EQUATION == E_EULER_COMP))&&(IN_OUT_FLOW_BC > 0)) if (((RDE_EQUATION == E_EULER_INCOMP)||(RDE_EQUATION == E_EULER_COMP))&&(IN_OUT_FLOW_BC > 0))
@ -1846,11 +1896,17 @@ double gfield[2*NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY])
} }
case (E_KELLER_SEGEL): case (E_KELLER_SEGEL):
{ {
// free(keller_segel_drift);
// free(keller_segel_div);
free(nabla_psi); free(nabla_psi);
break; break;
} }
case (E_KELLER_SEGEL_SPHERE):
{
free(delta_u);
free(delta_v);
free(chi_u);
free(prod_gradients);
break;
}
} }
if (COMPUTE_PRESSURE) if (COMPUTE_PRESSURE)
@ -2339,7 +2395,7 @@ void animation()
// if (ADD_TRACERS) tracers = (double *)malloc(2*NSTEPS*N_TRACERS*sizeof(double)); // if (ADD_TRACERS) tracers = (double *)malloc(2*NSTEPS*N_TRACERS*sizeof(double));
if (ADD_TRACERS) tracers = (double *)malloc(4*NSTEPS*N_TRACERS*sizeof(double)); if (ADD_TRACERS) tracers = (double *)malloc(4*NSTEPS*N_TRACERS*sizeof(double));
if (RDE_EQUATION == E_KURAMOTO_SPHERE) if ((RDE_EQUATION == E_KURAMOTO_SPHERE)||(RDE_EQUATION == E_KELLER_SEGEL_SPHERE))
ngridpoints = initialize_simulation_grid_sphere(wsphere); ngridpoints = initialize_simulation_grid_sphere(wsphere);
dx = (XMAX-XMIN)/((double)NX); dx = (XMAX-XMIN)/((double)NX);
@ -2418,7 +2474,7 @@ void animation()
// init_linear_blob_sphere(0, 1.3*PI, 0.65*PI, 0.0, 0.0, 0.3, 0.1, 0.1, SWATER_MIN_HEIGHT, phi, xy_in, wsphere); // init_linear_blob_sphere(0, 1.3*PI, 0.65*PI, 0.0, 0.0, 0.3, 0.1, 0.1, SWATER_MIN_HEIGHT, phi, xy_in, wsphere);
init_random(0.4, 0.4, phi, xy_in, wsphere); init_random(0.5, 0.4, phi, xy_in, wsphere);
// init_onefield_random(1, 0.2, 0.1, phi, xy_in, wsphere); // init_onefield_random(1, 0.2, 0.1, phi, xy_in, wsphere);
// init_random_smoothed(0.0, 0.4, phi, xy_in, wsphere); // init_random_smoothed(0.0, 0.4, phi, xy_in, wsphere);
@ -2470,6 +2526,7 @@ void animation()
for (i=0; i<=NSTEPS; i++) for (i=0; i<=NSTEPS; i++)
{ {
printf("floored chi %i times\n", nfloorchi);
nvid = NVID; nvid = NVID;
if (CHANGE_VISCOSITY) if (CHANGE_VISCOSITY)
{ {

1
schrodinger.c
View File

@ -196,6 +196,7 @@
#define WAVE_PACKET_SOURCE_TYPE 1 /* type of wave packet sources */ #define WAVE_PACKET_SOURCE_TYPE 1 /* type of wave packet sources */
#define N_WAVE_PACKETS 15 /* number of wave packets */ #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 OSCIL_LEFT_YSHIFT 0.0 /* y-dependence of left oscillation (for non-horizontal waves) */
#define OSCIL_YMID -0.9 /* defines oscilling beam midpoint */
#define OSCILLATING_SOURCE_PERIOD 20 /* period of oscillating source */ #define OSCILLATING_SOURCE_PERIOD 20 /* period of oscillating source */
#define MU_B 1.0 /* parameter controlling the dimensions of domain */ #define MU_B 1.0 /* parameter controlling the dimensions of domain */
#define GAMMA 0.0 /* damping factor in wave equation */ #define GAMMA 0.0 /* damping factor in wave equation */

430
sub_hashgrid.c
View File

@ -2,6 +2,9 @@
/* Note: a possible improvement would be to use pointers instead of integers */ /* Note: a possible improvement would be to use pointers instead of integers */
/* when referencing other hashgrid cells or particles */ /* when referencing other hashgrid cells or particles */
int hashx_sphere[HASHY]; /* number of hash cells in x direction for sphere bc */
double dx_sphere[HASHY]; /* width of hash cells in x direction for sphere bc */
double module2(double x, double y) /* Euclidean norm */ double module2(double x, double y) /* Euclidean norm */
{ {
double m; double m;
@ -32,6 +35,7 @@ int bc_grouped(int bc)
case (BC_REFLECT_ABS): return(0); case (BC_REFLECT_ABS): return(0);
case (BC_REFLECT_ABS_BOTTOM): return(0); case (BC_REFLECT_ABS_BOTTOM): return(0);
case (BC_REFLECT_ABS_RIGHT): return(0); case (BC_REFLECT_ABS_RIGHT): return(0);
case (BC_SPHERE): return(5);
default: default:
{ {
printf("Warning: Hashgrid will not be properly initialised, update bc_grouped()\n\n"); printf("Warning: Hashgrid will not be properly initialised, update bc_grouped()\n\n");
@ -52,7 +56,7 @@ int hash_cell(double x, double y)
/* returns number of hash cell */ /* returns number of hash cell */
{ {
static int first = 1; static int first = 1;
static double lx, ly, padding; static double lx, ly, padding, dy_inverse;
int i, j; int i, j;
if (first) if (first)
@ -61,19 +65,35 @@ int hash_cell(double x, double y)
else padding = HASHGRID_PADDING; else padding = HASHGRID_PADDING;
lx = BCXMAX - BCXMIN + 2.0*padding; lx = BCXMAX - BCXMIN + 2.0*padding;
ly = BCYMAX - (BCYMIN) + 2.0*padding; ly = BCYMAX - (BCYMIN) + 2.0*padding;
dy_inverse = (double)HASHY/(PI - 2.0*POLAR_PADDING);
first = 0; first = 0;
} }
if (CENTER_VIEW_ON_OBSTACLE) x -= xshift; if (bc_grouped(BOUNDARY_COND) == 5) /* spherical bc */
{
j = (int)((y-POLAR_PADDING)*dy_inverse);
if (j<0) j = 0;
else if (j>=HASHY) j = HASHY-1;
i = (int)(x/dx_sphere[j]);
if (i<0) i = 0;
else if (i>=hashx_sphere[j]) i = hashx_sphere[j]-1;
}
else
{
if (CENTER_VIEW_ON_OBSTACLE) x -= xshift;
i = (int)((double)HASHX*(x - BCXMIN + padding)/lx); i = (int)((double)HASHX*(x - BCXMIN + padding)/lx);
j = (int)((double)HASHY*(y - BCYMIN + padding)/ly); j = (int)((double)HASHY*(y - BCYMIN + padding)/ly);
if (i<0) i = 0; if (i<0) i = 0;
else if (i>=HASHX) i = HASHX-1; else if (i>=HASHX) i = HASHX-1;
if (j<0) j = 0; if (j<0) j = 0;
else if (j>=HASHY) j = HASHY-1; else if (j>=HASHY) j = HASHY-1;
}
// printf("Hash_cell (%.3lg, %.3lg): %i\n", x, y, mhash(i,j));
return(mhash(i,j)); return(mhash(i,j));
// printf("Mapped (%.3f,%.3f) to (%i, %i)\n", x, y, ij[0], ij[1]); // printf("Mapped (%.3f,%.3f) to (%i, %i)\n", x, y, ij[0], ij[1]);
} }
@ -85,7 +105,7 @@ void init_hashcell_coordinates(t_hashgrid hashgrid[HASHX*HASHY])
static int first = 1; static int first = 1;
static double lx, ly, padding, dx, dy; static double lx, ly, padding, dx, dy;
int i, j, n; int i, j, n;
double x, y; double x, y, dummy = -100.0;
if (first) if (first)
{ {
@ -94,11 +114,50 @@ void init_hashcell_coordinates(t_hashgrid hashgrid[HASHX*HASHY])
lx = BCXMAX - BCXMIN + 2.0*padding; lx = BCXMAX - BCXMIN + 2.0*padding;
ly = BCYMAX - (BCYMIN) + 2.0*padding; ly = BCYMAX - (BCYMIN) + 2.0*padding;
dx = 1.0*lx/(double)HASHX; dx = 1.0*lx/(double)HASHX;
dy = 1.0*ly/(double)HASHY; if (bc_grouped(BOUNDARY_COND) == 5)
{
dy = (PI - 2.0*POLAR_PADDING)/(double)(HASHY);
}
else dy = 1.0*ly/(double)HASHY;
first = 0; first = 0;
} }
for (i=0; i<HASHX; i++) if (bc_grouped(BOUNDARY_COND) == 5) /* spherical bc */
{
for (j=0; j<HASHY; j++)
{
for (i=0; i<hashx_sphere[j]; i++)
{
n = mhash(i, j);
hashgrid[n].x1 = (double)i*dx_sphere[j];
hashgrid[n].x2 = (double)(i+1)*dx_sphere[j];
hashgrid[n].y1 = POLAR_PADDING + (double)j*dy;
hashgrid[n].y2 = POLAR_PADDING + (double)(j+1)*dy;
/* TODO: correct area */
hashgrid[n].area = dx_sphere[j]*dy*sin(0.5*(hashgrid[n].y1 + hashgrid[n].y2));
printf("Cell %i corners (%.3lg, %.3lg), (%.3lg, %.3lg)\n", n, hashgrid[n].x1, hashgrid[n].y1, hashgrid[n].x2, hashgrid[n].y2);
}
for (i=hashx_sphere[j]; i<HASHX; i++)
{
n = mhash(i, j);
hashgrid[n].x1 = dummy;
hashgrid[n].x2 = dummy;
hashgrid[n].y1 = dummy;
hashgrid[n].y2 = dummy;
hashgrid[n].area = 1.0;
}
}
/* correction at poles */
n = mhash(0,0);
hashgrid[n].y1 = 0.0;
hashgrid[n].area = DPI*(POLAR_PADDING + dy)*sin(0.5*hashgrid[n].y2);
n = mhash(0, HASHY-1);
hashgrid[n].y2 = PI;
hashgrid[n].area = DPI*(POLAR_PADDING + dy)*sin(0.5*(PI + hashgrid[n].y1));
}
else for (i=0; i<HASHX; i++)
for (j=0; j<HASHY; j++) for (j=0; j<HASHY; j++)
{ {
x = BCXMIN - padding + (double)i*dx; x = BCXMIN - padding + (double)i*dx;
@ -110,19 +169,106 @@ void init_hashcell_coordinates(t_hashgrid hashgrid[HASHX*HASHY])
hashgrid[n].y2 = y + dy; hashgrid[n].y2 = y + dy;
hashgrid[n].charge = 0.0; hashgrid[n].charge = 0.0;
} }
} }
void init_hashgrid(t_hashgrid hashgrid[HASHX*HASHY]) void init_hashgrid(t_hashgrid hashgrid[HASHX*HASHY])
/* initialise table of neighbouring cells for each hashgrid cell, depending on boundary condition */ /* initialise table of neighbouring cells for each hashgrid cell, depending on boundary condition */
{ {
int i, j, k, p, q, m, i1, j1; int i, j, k, p, q, m, i1, j1, m1, pmin, pmax;
double dy, x1, x2, xx1, xx2, padding;
short int sym;
printf("Initializing hash grid\n"); printf("Initializing hash grid\n");
/* initialize hashx_sphere[] and dx_sphere[] in case of spherical bc */
if (bc_grouped(BOUNDARY_COND) == 5)
{
dy = PI/((double)HASHY);
for (j=0; j<HASHY/2 + 1; j++)
{
hashx_sphere[j] = (int)((double)HASHX*(sin(dy*(double)j)));
if (hashx_sphere[j] == 0) hashx_sphere[j] = 1;
if (hashx_sphere[j] > HASHX-1) hashx_sphere[j] = HASHX-1;
dx_sphere[j] = DPI/(double)hashx_sphere[j];
printf("hashx_sphere[%i] = %i, dx_sphere[%i] = %.3lg\n", j, hashx_sphere[j], j, dx_sphere[j]);
fprintf(lj_log, "hashx_sphere[%i] = %i, dx_sphere[%i] = %.3lg\n", j, hashx_sphere[j], j, dx_sphere[j]);
}
for (j=HASHY/2 + 1; j<HASHY; j++)
{
hashx_sphere[j] = hashx_sphere[HASHY-1-j];
dx_sphere[j] = dx_sphere[HASHY-1-j];
printf("hashx_sphere[%i] = %i, dx_sphere[%i] = %.3lg\n", j, hashx_sphere[j], j, dx_sphere[j]);
}
}
if ((COLOR_BACKGROUND)||(bc_grouped(BOUNDARY_COND) == 5))
init_hashcell_coordinates(hashgrid);
/* bulk of the table */ /* bulk of the table */
for (i=0; i<HASHX-1; i++) if (bc_grouped(BOUNDARY_COND) == 5) /* spherical bc */
{
/* dummy values for safety */
for (i=0; i<HASHX*HASHY; i++) hashgrid[i].nneighb = 0;
for (j=1; j<HASHY-1; j++)
{
padding = 1.0*dx_sphere[j];
for (i=1; i<hashx_sphere[j]-1; i++)
{
m = mhash(i, j);
// printf("m = %i\n", m);
hashgrid[m].nneighb = 3;
hashgrid[m].neighbour[0] = mhash(i-1,j);
hashgrid[m].neighbour[1] = mhash(i,j);
hashgrid[m].neighbour[2] = mhash(i+1,j);
/* neighbours in layer above and below */
if (j==1) pmin = 1;
else pmin = -1;
if (j==HASHY-2) pmax = 0;
else pmax = 2;
for (p=pmin; p<pmax; p+=2)
for (i1=0; i1<hashx_sphere[j+p]; i1++)
{
m1 = mhash(i1, j+p);
x1 = hashgrid[m].x1 - padding;
x2 = hashgrid[m].x2 + padding;
xx1 = hashgrid[m1].x1;
xx2 = hashgrid[m1].x2;
if (((xx2 >= x1)&&(xx2 <= x2))||((xx1 >= x1)&&(xx1 <= x2)))
{
hashgrid[m].neighbour[hashgrid[m].nneighb] = m1;
hashgrid[m].nneighb++;
}
}
/* extra treatment for pole neighbours */
if (j==1)
{
hashgrid[m].neighbour[hashgrid[m].nneighb] = 0;
hashgrid[m].nneighb++;
}
if (j==HASHY-2)
{
hashgrid[m].neighbour[hashgrid[m].nneighb] = HASHY-1;
hashgrid[m].nneighb++;
}
// printf("hashgrid[%i].nneighb = %i\n", m, hashgrid[m].nneighb);
/* check there are not too many neighbours */
if (hashgrid[m].nneighb >= HASHMAXNEIGH)
{
printf("(i,j) = (%i,%i) Error: HASHMAXNEIGH should be at least %i\n", i, j, hashgrid[m].nneighb + 1);
exit(1);
}
}
}
}
else for (i=0; i<HASHX-1; i++)
for (j=0; j<HASHY-1; j++) for (j=0; j<HASHY-1; j++)
{ {
m = mhash(i, j); m = mhash(i, j);
@ -392,6 +538,201 @@ void init_hashgrid(t_hashgrid hashgrid[HASHX*HASHY])
/* TO DO : add more cells for "corners" ? */ /* TO DO : add more cells for "corners" ? */
break; break;
} }
case (5): /* spherical boundary conditions */
{
printf("Left border\n");
/* left border */
for (j=1; j<HASHY-1; j++)
{
padding = dx_sphere[j];
m = mhash(0, j);
// printf("j = %i, m = %i\n", j, m);
if (j==1)
{
hashgrid[m].nneighb = 4;
hashgrid[m].neighbour[0] = mhash(hashx_sphere[j]-1,j);
hashgrid[m].neighbour[1] = mhash(0,j);
hashgrid[m].neighbour[2] = mhash(1,j);
hashgrid[m].neighbour[3] = mhash(hashx_sphere[j+1]-1,j+1);
}
else if (j==HASHY-2)
{
hashgrid[m].nneighb = 4;
hashgrid[m].neighbour[0] = mhash(hashx_sphere[j]-1,j);
hashgrid[m].neighbour[1] = mhash(0,j);
hashgrid[m].neighbour[2] = mhash(1,j);
hashgrid[m].neighbour[3] = mhash(hashx_sphere[j-1]-1,j-1);
}
else
{
hashgrid[m].nneighb = 5;
hashgrid[m].neighbour[0] = mhash(hashx_sphere[j]-1,j);
hashgrid[m].neighbour[1] = mhash(0,j);
hashgrid[m].neighbour[2] = mhash(1,j);
hashgrid[m].neighbour[3] = mhash(hashx_sphere[j+1]-1,j+1);
hashgrid[m].neighbour[4] = mhash(hashx_sphere[j-1]-1,j-1);
// hashgrid[m].nneighb = 7;
// hashgrid[m].neighbour[5] = mhash(2,j);
// hashgrid[m].neighbour[6] = mhash(hashx_sphere[j]-2,j);
}
/* neighbours in layer above and below */
for (p=-1; p<2; p+=2)
for (i1=0; i1<hashx_sphere[j+p]; i1++)
{
m1 = mhash(i1, j+p);
x1 = hashgrid[m].x1 - padding;
x2 = hashgrid[m].x2 + padding;
xx1 = hashgrid[m1].x1;
xx2 = hashgrid[m1].x2;
if (((xx2 >= x1)&&(xx2 <= x2))||((xx1 >= x1)&&(xx1 <= x2)))
{
hashgrid[m].neighbour[hashgrid[m].nneighb] = m1;
hashgrid[m].nneighb++;
}
}
/* check there are not too many neighbours */
if (hashgrid[m].nneighb >= HASHMAXNEIGH)
{
printf("(i,j) = (0,%i) Error: HASHMAXNEIGH should be at least %i\n", j, hashgrid[m].nneighb + 1);
exit(1);
}
}
printf("Right border\n");
/* right border */
for (j=1; j<HASHY-1; j++)
{
padding = dx_sphere[j];
m = mhash(hashx_sphere[j]-1, j);
if (j==1)
{
hashgrid[m].nneighb = 4;
hashgrid[m].neighbour[0] = mhash(hashx_sphere[j]-2,j);
hashgrid[m].neighbour[1] = mhash(hashx_sphere[j]-1,j);
hashgrid[m].neighbour[2] = mhash(0,j);
hashgrid[m].neighbour[3] = mhash(0,j+1);
}
else if (j==HASHY-2)
{
hashgrid[m].nneighb = 4;
hashgrid[m].neighbour[0] = mhash(hashx_sphere[j]-2,j);
hashgrid[m].neighbour[1] = mhash(hashx_sphere[j]-1,j);
hashgrid[m].neighbour[2] = mhash(0,j);
hashgrid[m].neighbour[3] = mhash(0,j-1);
}
else
{
hashgrid[m].nneighb = 5;
hashgrid[m].neighbour[0] = mhash(hashx_sphere[j]-2,j);
hashgrid[m].neighbour[1] = mhash(hashx_sphere[j]-1,j);
hashgrid[m].neighbour[2] = mhash(0,j);
hashgrid[m].neighbour[3] = mhash(0,j+1);
hashgrid[m].neighbour[4] = mhash(0,j-1);
// hashgrid[m].nneighb = 7;
// hashgrid[m].neighbour[5] = mhash(1,j);
// hashgrid[m].neighbour[6] = mhash(hashx_sphere[j]-3,j);
}
/* neighbours in layer above and below */
for (p=-1; p<2; p+=2)
for (i1=0; i1<hashx_sphere[j+p]; i1++)
{
m1 = mhash(i1, j+p);
x1 = hashgrid[m].x1 - padding;
x2 = hashgrid[m].x2 + padding;
xx1 = hashgrid[m1].x1;
xx2 = hashgrid[m1].x2;
if (((xx2 >= x1)&&(xx2 <= x2))||((xx1 >= x1)&&(xx1 <= x2)))
{
hashgrid[m].neighbour[hashgrid[m].nneighb] = m1;
hashgrid[m].nneighb++;
}
}
/* check there are not too many neighbours */
if (hashgrid[m].nneighb >= HASHMAXNEIGH)
{
printf("(i,j) = (%i,%i) Error: HASHMAXNEIGH should be at least %i\n", hashx_sphere[j]-1, j, hashgrid[m].nneighb + 1);
exit(1);
}
}
printf("Top border\n");
/* top border/North pole */
m = mhash(0, HASHY-1);
hashgrid[m].nneighb = hashx_sphere[HASHY-2] + 1;
hashgrid[m].neighbour[0] = m;
fprintf(lj_log, "Top border: m = %i, %i neighbours\n", m, hashgrid[m].nneighb);
for (i1=0; i1<hashx_sphere[HASHY-2]; i1++)
hashgrid[m].neighbour[i1+1] = mhash(i1,HASHY-2);
/* check there are not too many neighbours */
if (hashgrid[m].nneighb >= HASHMAXNEIGH)
{
printf("Error at North Pole: HASHMAXNEIGH should be at least %i\n", hashgrid[m].nneighb + 1);
exit(1);
}
printf("Bottom border\n");
/* bottom border/South pole */
m = mhash(0, 0);
hashgrid[m].nneighb = hashx_sphere[1] + 1;
hashgrid[m].neighbour[0] = m;
for (i1=0; i1<hashx_sphere[1]; i1++)
hashgrid[m].neighbour[i1+1] = mhash(i1,1);
/* check there are not too many neighbours */
if (hashgrid[m].nneighb >= HASHMAXNEIGH)
{
printf("Error at South Pole: HASHMAXNEIGH should be at least %i\n", hashgrid[m].nneighb + 1);
exit(1);
}
/* symmetrize hashgrid */
for (i=0; i<HASHX*HASHY; i++)
{
for (j=0; j<hashgrid[i].nneighb; j++)
{
m = hashgrid[i].neighbour[j];
sym = 0;
for (k=0; k<hashgrid[m].nneighb; k++)
if (hashgrid[m].neighbour[k] == i) sym = 1;
if (!sym)
{
fprintf(lj_log, "Symmetrising hashcells %i and %i\n", i, m);
hashgrid[m].neighbour[hashgrid[m].nneighb] = i;
hashgrid[m].nneighb++;
}
}
}
for (i=0; i<HASHX*HASHY; i++)
{
p = hashgrid[i].nneighb;
printf("Hash cell %i has %i neighbours: ", i, p);
fprintf(lj_log, "Hash cell %i has %i neighbours: ", i, p);
for (j=0; j<p; j++)
{
printf("%i ", hashgrid[i].neighbour[j]);
fprintf(lj_log, "%i ", hashgrid[i].neighbour[j]);
}
printf("\n");
fprintf(lj_log, "\n");
}
// sleep(5);
printf("Hashgrid initialised\n");
break;
}
default: /* do nothing */; default: /* do nothing */;
} }
@ -413,14 +754,14 @@ void init_hashgrid(t_hashgrid hashgrid[HASHX*HASHY])
// // sleep(1); // // sleep(1);
// } // }
if (COLOR_BACKGROUND) init_hashcell_coordinates(hashgrid); // if (COLOR_BACKGROUND) init_hashcell_coordinates(hashgrid);
// sleep(1); // sleep(1);
} }
void update_hashgrid(t_particle* particle, t_obstacle *obstacle, t_hashgrid* hashgrid, int verbose) void update_hashgrid(t_particle* particle, t_obstacle *obstacle, t_hashgrid* hashgrid, int verbose)
{ {
int i, j, k, n, m, max = 0, hashcell; int i, j, k, n, m, max = 0, hashcell, maxcell = -1;
// printf("Updating hashgrid_number\n"); // printf("Updating hashgrid_number\n");
for (i=0; i<HASHX*HASHY; i++) hashgrid[i].number = 0; for (i=0; i<HASHX*HASHY; i++) hashgrid[i].number = 0;
@ -438,7 +779,11 @@ void update_hashgrid(t_particle* particle, t_obstacle *obstacle, t_hashgrid* has
hashgrid[hashcell].number++; hashgrid[hashcell].number++;
particle[k].hashcell = hashcell; particle[k].hashcell = hashcell;
if (n > max) max = n; if (n > max)
{
max = n;
maxcell = hashcell;
}
} }
if (OSCILLATE_OBSTACLES) if (OSCILLATE_OBSTACLES)
@ -454,7 +799,7 @@ void update_hashgrid(t_particle* particle, t_obstacle *obstacle, t_hashgrid* has
} }
} }
if(verbose) printf("Maximal number of particles per hash cell: %i\n", max); if(verbose) printf("Maximal number of particles per hash cell: %i, in cell %i\n", max, maxcell);
} }
@ -669,6 +1014,43 @@ int wrap_particle(t_particle* particle, double *px, double *py)
return(move); return(move);
break; break;
} }
case (5): /* spherical bc */
{
if (x < 0.0)
{
x1 += DPI;
// printf("Moving particle by 2Pi\n");
move++;
}
else if (x > DPI)
{
x1 -= DPI;
// printf("Moving particle by -2Pi\n");
move++;
}
if (y > PI)
{
x1 += PI;
if (x1 > DPI) x1 -= DPI;
y1 = PI;
particle->vy *= -1.0;
*py *= -1.0;
move++;
}
else if (y < 0.0)
{
x1 += PI;
if (x1 > DPI) x1 -= DPI;
y1 = 0.0;
particle->vy *= -1.0;
*py *= -1.0;
move++;
}
particle->xc = x1;
particle->yc = y1;
return(move);
break;
}
default: default:
{ {
/* do nothing */ /* do nothing */
@ -809,7 +1191,17 @@ int verbose = 0;
// printf("(x2,y2) = (%.3lg, %.3lg)\n", *x2, *y2); // printf("(x2,y2) = (%.3lg, %.3lg)\n", *x2, *y2);
break; break;
} }
} case (5): /* spherical b.c. */
{
if (dx > dxhalf) *x2 -= BCXMAX - BCXMIN;
else if (-dx > dxhalf) *x2 += BCXMAX - BCXMIN;
// if (dy > dyhalf) *y2 -= BCYMAX - BCYMIN;
// else if (-dy > dyhalf) *y2 += BCYMAX - BCYMIN;
/* TODO: poles? */
// printf("(x2,y2) = (%.3lg, %.3lg)\n", *x2, *y2);
break;
}
}
} }

1426
sub_lj.c
View File

File diff suppressed because it is too large Load Diff

2171
sub_lj_sphere.c Normal file
View File

File diff suppressed because it is too large Load Diff

276
sub_rde.c
View File

@ -5941,6 +5941,132 @@ void compute_kuramoto_interaction_grid(double phi_in[NX*NY], double phi_out[NX*N
} }
} }
void compute_laplacian_grid(double phi_in[NX*NY], double phi_out[NX*NY], t_wave_sphere wsphere[NX*NY], int ngridpoints)
/* computes laplacian on adapted grid */
/* assumes each grid point has 4 neighbours, result is not divided by dx */
{
int i, k, n;
#pragma omp parallel for private(i,k)
for (i=0; i<ngridpoints; i++)
{
phi_out[i] = -4.0*phi_in[i];
for (k = 0; k < 4; k++)
{
n = wsphere[i].neighbor[k];
phi_out[i] += phi_in[n];
}
}
}
void compute_gradient_grid(double phi_in[NX*NY], double gradient[2*NX*NY], t_wave_sphere wsphere[NX*NY], int ngridpoints)
/* computes gradient on adapted grid */
/* assumes each grid point has 4 neighbours, result is not divided by dx */
{
int i, n;
double gradx, grady;
#pragma omp parallel for private(i)
for (i=0; i<ngridpoints; i++)
{
n = wsphere[i].neighbor[1];
gradx = phi_in[n];
n = wsphere[i].neighbor[0];
gradx -= phi_in[n];
n = wsphere[i].neighbor[3];
grady = phi_in[n];
n = wsphere[i].neighbor[2];
grady -= phi_in[n];
gradient[i] = wsphere[i].cos_grid*gradx - wsphere[i].sin_grid*grady;
gradient[NX*NY+i] = wsphere[i].sin_grid*gradx + wsphere[i].cos_grid*grady;
// n = wsphere[i].neighbor[1];
// gradient[i] = phi_in[n];
// n = wsphere[i].neighbor[0];
// gradient[i] -= phi_in[n];
//
//
// n = wsphere[i].neighbor[3];
// gradient[NX*NY+i] = phi_in[n];
// n = wsphere[i].neighbor[2];
// gradient[NX*NY+i] -= phi_in[n];
}
}
void compute_divergence_grid(double phi_in[NX*NY], double psi_in[NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY], int ngridpoints)
/* computes divergence on adapted grid */
/* assumes each grid point has 4 neighbours, result is not divided by dx */
{
int i, n;
double gradx, grady;
#pragma omp parallel for private(i)
for (i=0; i<ngridpoints; i++)
{
n = wsphere[i].neighbor[1];
gradx = phi_in[n];
n = wsphere[i].neighbor[0];
gradx -= phi_in[n];
n = wsphere[i].neighbor[3];
grady = phi_in[n];
n = wsphere[i].neighbor[2];
grady -= phi_in[n];
rde[i].divergence = wsphere[i].cos_grid*gradx - wsphere[i].sin_grid*grady;
n = wsphere[i].neighbor[1];
gradx = psi_in[n];
n = wsphere[i].neighbor[0];
gradx -= psi_in[n];
n = wsphere[i].neighbor[3];
grady = psi_in[n];
n = wsphere[i].neighbor[2];
grady -= psi_in[n];
rde[i].divergence += wsphere[i].sin_grid*gradx + wsphere[i].cos_grid*grady;
// n = wsphere[i].neighbor[1];
// rde[i].divergence = phi_in[n];
// n = wsphere[i].neighbor[0];
// rde[i].divergence -= phi_in[n];
// n = wsphere[i].neighbor[3];
// rde[i].divergence += psi_in[n];
// n = wsphere[i].neighbor[2];
// rde[i].divergence -= psi_in[n];
}
}
void compute_grad_product_grid(double phi[NX*NY], double psi[NX*NY], double prod_gradient[NX*NY], t_wave_sphere wsphere[NX*NY], int ngridpoints)
/* computes gradient on adapted grid */
/* assumes each grid point has 4 neighbours, result is not divided by dx */
{
int i, n;
double gradx1, grady1, gradx2, grady2;
#pragma omp parallel for private(i)
for (i=0; i<ngridpoints; i++)
{
n = wsphere[i].neighbor[1];
gradx1 = phi[n];
gradx2 = psi[n];
n = wsphere[i].neighbor[0];
gradx1 -= phi[n];
gradx2 -= psi[n];
n = wsphere[i].neighbor[3];
grady1 = phi[n];
grady2 = psi[n];
n = wsphere[i].neighbor[2];
grady1 -= phi[n];
grady2 -= psi[n];
prod_gradient[i] = 0.25*(gradx1*gradx2 + grady1*grady2);
}
}
void compute_light_angle_sphere_rde(short int xy_in[NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY], double potential[NX*NY], int movie, int transparent) void compute_light_angle_sphere_rde(short int xy_in[NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY], double potential[NX*NY], int movie, int transparent)
/* computes cosine of angle between normal vector and vector light */ /* computes cosine of angle between normal vector and vector light */
{ {
@ -6196,7 +6322,7 @@ void compute_light_angle_sphere_rde_2d(short int xy_in[NX*NY], t_rde rde[NX*NY],
first = 0; first = 0;
} }
printf("[compute_light_angle_sphere_rde_2d] computing gradient\n"); // printf("[compute_light_angle_sphere_rde_2d] computing gradient\n");
#pragma omp parallel for private(i,j,gradx, grady, norm, pscal) #pragma omp parallel for private(i,j,gradx, grady, norm, pscal)
for (i=1; i<NX-1; i++) for (i=1; i<NX-1; i++)
@ -6300,7 +6426,7 @@ void compute_light_angle_sphere_rde_2d(short int xy_in[NX*NY], t_rde rde[NX*NY],
else rde[(NX-1)*NY+j].cos_angle = wsphere[(NX-1)*NY+j].cos_angle; else rde[(NX-1)*NY+j].cos_angle = wsphere[(NX-1)*NY+j].cos_angle;
} }
printf("[compute_light_angle_sphere_rde_2d] computing gradient done\n"); // printf("[compute_light_angle_sphere_rde_2d] computing gradient done\n");
} }
@ -8964,7 +9090,7 @@ int sphere_to_cube(t_wave_sphere *wsphere, double coord[2])
int max; int max;
theta = wsphere->theta; theta = wsphere->theta;
phi = wsphere->phi + PI; phi = wsphere->phi;
x = sin(theta)*cos(phi); x = sin(theta)*cos(phi);
y = sin(theta)*sin(phi); y = sin(theta)*sin(phi);
z = cos(theta); z = cos(theta);
@ -9017,7 +9143,7 @@ int generate_cubic_grid_sphere(t_wave_sphere wsphere[NX*NY])
{ {
int i, j, k, l, n, m, face, nerrors; int i, j, k, l, n, m, face, nerrors;
double dx, dx2, x, y, z, coord[2]; double dx, dx2, x, y, z, coord[2];
double phi0, theta0, phi1, theta1; double phi0, theta0, phi1, theta1, angle;
/* number of sites per face is l */ /* number of sites per face is l */
l = (int)(sqrt((double)(NX*NY/6))); l = (int)(sqrt((double)(NX*NY/6)));
@ -9172,7 +9298,8 @@ int generate_cubic_grid_sphere(t_wave_sphere wsphere[NX*NY])
for (i=0; i<l; i++) for (i=0; i<l; i++)
{ {
wsphere[3*l*l + i].neighbor[2] = 5*l*l + i*l; wsphere[3*l*l + i].neighbor[2] = 5*l*l + i*l;
wsphere[5*l*l + i*l].neighbor[3] = 3*l*l + i; // wsphere[5*l*l + i*l].neighbor[3] = 3*l*l + i;
wsphere[5*l*l + i*l].neighbor[0] = 3*l*l + i;
wsphere[3*l*l + i].edge = 1; wsphere[3*l*l + i].edge = 1;
wsphere[5*l*l + i*l].edge = 1; wsphere[5*l*l + i*l].edge = 1;
} }
@ -9251,26 +9378,134 @@ int generate_cubic_grid_sphere(t_wave_sphere wsphere[NX*NY])
} }
/* test */ /* test */
// nerrors = 0; nerrors = 0;
// for (n=0; n<NX*NY; n++) for (n=0; n<NX*NY; n++)
// { {
// phi0 = wsphere[n].phi; phi0 = wsphere[n].phi;
// theta0 = wsphere[n].theta; theta0 = wsphere[n].theta;
// m = wsphere[n].convert_grid; m = wsphere[n].convert_grid;
// phi1 = wsphere[m].phigrid; phi1 = wsphere[m].phigrid;
// theta1 = wsphere[m].thetagrid; theta1 = wsphere[m].thetagrid;
// if (dist_sphere(phi0, theta0, phi1, theta1) > 0.2) if (dist_sphere(phi0, theta0, phi1, theta1) > 0.01)
// { {
// printf("Warning: Face %i, grid points %i (%.0f, %.0f) and its grid image (%.0f, %.0f) seem too far apart\n", m/(l*l), n, phi0*180.0/PI, theta0*180.0/PI, phi1*180.0/PI, theta1*180.0/PI); printf("Warning: Face %i, grid points %i (%.0f, %.0f) and its grid image (%.0f, %.0f) seem too far apart\n", m/(l*l), n, phi0*180.0/PI, theta0*180.0/PI, phi1*180.0/PI, theta1*180.0/PI);
// nerrors++; nerrors++;
// } }
// } }
// printf("%i errors among %i grid points\n", nerrors, 6*l*l); printf("%i errors among %i grid points\n", nerrors, 6*l*l);
/* compute orientation of grid wrt sphere */
// for (i=0; i<6*l*l; i++)
// {
// phi0 = wsphere[i].phigrid;
// theta0 = wsphere[i].thetagrid;
// n = wsphere[i].neighbor[0];
// phi1 = wsphere[n].phigrid;
// theta1 = wsphere[n].thetagrid;
// angle = argument(phi1 - phi0, theta1 - theta0);
// // angle*= cos(PID*wsphere[i].z);
// wsphere[i].cos_grid = cos(angle);
// wsphere[i].sin_grid = sin(angle);
// // printf("Grid angle %i = %.3lg, edge = %i\n", i, angle*180.0/PI, wsphere[i].edge);
// }
return(6*l*l); return(6*l*l);
} }
void optimize_cubic_grid_sphere(t_wave_sphere wsphere[NX*NY], int npoints, int nsteps)
/* optimize the cubic grid by evolving the grid points */
{
int i, j, k, n, m, p, nmoved;
double phi0, theta0, phi1, theta1, f, *fphi, *ftheta, dist, dist2, dt = 2.0e-6, dmin, dmax;
double phimax0, phimax1, thetamax0, thetamax1, dphi, sinthetamin = 0.05;
int imax, jmax, kmax, edge;
printf("Optimizing cubic grid\n");
fphi = (double *)malloc(npoints*sizeof(double));
ftheta = (double *)malloc(npoints*sizeof(double));
for (n=0; n<nsteps; n++)
{
dmin = 10.0;
dmax = 0.0;
printf("Step %i of %i\n", n, nsteps);
// #pragma omp parallel for private(i)
for (i=0; i<npoints; i++) if (sin(wsphere[i].thetagrid) > sinthetamin)
{
fphi[i] = 0.0;
ftheta[i] = 0.0;
phi0 = wsphere[i].phigrid;
theta0 = wsphere[i].thetagrid;
for (j=0; j<4; j++)
{
k = wsphere[i].neighbor[j];
phi1 = wsphere[k].phigrid;
theta1 = wsphere[k].thetagrid;
dphi = phi1 - phi0;
if (dphi > PI) dphi -= DPI;
else if (dphi < -PI) dphi += DPI;
// if (phi1 > phi0 + PI) phi1 -= DPI;
// if (phi1 < phi0 - PI) phi1 += DPI;
dist = dist_sphere(phi0, theta0, phi1, theta1);
if (dist < dmin) dmin = dist;
if (dist > dmax)
{
dmax = dist;
phimax0 = phi0;
thetamax0 = theta0;
phimax1 = phi1;
thetamax1 = theta1;
imax = i;
jmax = j;
kmax = k;
}
// dist2 = module2(dphi, theta1 - theta0);
f = 1.0/(dist + 1.0);
fphi[i] -= dphi*f/dist;
ftheta[i] += (theta0 - theta1)*f/dist;
}
}
// #pragma omp parallel for private(i)
for (i=0; i<npoints; i++) if (sin(wsphere[i].thetagrid) > sinthetamin)
{
wsphere[i].phigrid += fphi[i]*dt;
if (wsphere[i].phigrid > DPI) wsphere[i].phigrid -= DPI;
if (wsphere[i].phigrid < 0.0) wsphere[i].phigrid += DPI;
wsphere[i].thetagrid += ftheta[i]*dt;
// printf("Moved point %i by (%.5lg, %.5lg)\n", i, fphi*dt, ftheta*dt);
}
printf("dmin = %.5lg, dmax = %.5lg\n", dmin, dmax);
// printf("dmax between pts %i (%.5lg,%.5lg) and %i (%.5lg,%.5lg) neighbor %i, edge = %i\n", imax, phimax0*180.0/PI, thetamax0*180.0/PI, kmax, phimax1*180.0/PI, thetamax1*180.0/PI, jmax, wsphere[imax].edge);
if (n%10 == 0)
{
/* recompute corresponding points in latitude-longitude grid */
printf("Recomputing grid correspondence\n");
nmoved = 0;
for (i=0; i<NX*NY; i++)
{
m = wsphere[i].convert_grid;
dist = dist_sphere(wsphere[i].phi, wsphere[i].theta, wsphere[m].phigrid, wsphere[m].thetagrid);
for (k=0; k<4; k++)
{
p = wsphere[m].neighbor[k];
dist2 = dist_sphere(wsphere[i].phi, wsphere[i].theta, wsphere[p].phigrid, wsphere[p].thetagrid);
if (dist2 < dist)
{
dist = dist2;
wsphere[i].convert_grid = p;
nmoved++;
}
}
}
printf("Moved %i points, distance = %.5lg\n", nmoved, dist);
}
}
free(fphi);
free(ftheta);
}
int initialize_simulation_grid_sphere(t_wave_sphere wsphere[NX*NY]) int initialize_simulation_grid_sphere(t_wave_sphere wsphere[NX*NY])
/* initialize Poisson random simulation grid on a sphere */ /* initialize Poisson random simulation grid on a sphere */
{ {
@ -9291,6 +9526,7 @@ int initialize_simulation_grid_sphere(t_wave_sphere wsphere[NX*NY])
case (GRID_CUBIC): case (GRID_CUBIC):
{ {
npoints = generate_cubic_grid_sphere(wsphere); npoints = generate_cubic_grid_sphere(wsphere);
if (OPTIMIZE_GRID) optimize_cubic_grid_sphere(wsphere, npoints, 600);
break; break;
} }
} }

362
sub_wave.c
View File

@ -4938,6 +4938,21 @@ int init_poly(int depth, t_vertex polyline[NMAXPOLY], t_rectangle polyrect[NMAXP
{ {
return(compute_disc_hex_lattice_strip_coordinates(polyline, polyarc, circles, npolyarc, ncircles, top)); return(compute_disc_hex_lattice_strip_coordinates(polyline, polyarc, circles, npolyarc, ncircles, top));
} }
case (D_TWO_PARABOLAS_ASYM_GUIDE):
{
polyrect[0].x1 = XMIN - WALL_WIDTH;
polyrect[0].y1 = OSCIL_YMID - OSCIL_YMAX - WALL_WIDTH_B;
polyrect[0].x2 = -0.5*MU;
polyrect[0].y2 = OSCIL_YMID - OSCIL_YMAX;
polyrect[1].x1 = XMIN - WALL_WIDTH;
polyrect[1].y1 = OSCIL_YMID + OSCIL_YMAX;
polyrect[1].x2 = -0.5*MU;
polyrect[1].y2 = OSCIL_YMID + OSCIL_YMAX + WALL_WIDTH_B;
*npolyrect = 2;
return(0);
}
default: default:
{ {
if ((ADD_POTENTIAL)&&(POTENTIAL == POT_MAZE)) if ((ADD_POTENTIAL)&&(POTENTIAL == POT_MAZE))
@ -5372,6 +5387,12 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_
else return(0); else return(0);
break; break;
} }
case (D_CIRCLE_LAYERS):
{
if (x*x + y*y < LAMBDA*LAMBDA) return(1);
else return(0);
break;
}
case (D_EXT_ELLIPSE): case (D_EXT_ELLIPSE):
{ {
if (x*x/(LAMBDA*LAMBDA) + y*y/(MU*MU) > 1.0) return(1); if (x*x/(LAMBDA*LAMBDA) + y*y/(MU*MU) > 1.0) return(1);
@ -5515,6 +5536,40 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_
else if ((x > x1)&&(x < x2)) return(1); else if ((x > x1)&&(x < x2)) return(1);
else return(0); else return(0);
} }
case (D_TWO_PARABOLAS_ASYM):
{
if (x > 0.0)
{
if (vabs(y) > 0.5*LAMBDA) return(1);
x1 = -y*y/LAMBDA + 0.25*LAMBDA;
return((x < x1)||(x > x1 + WALL_WIDTH));
}
else
{
if (vabs(y) > 0.5*MU) return(1);
x1 = y*y/MU - 0.25*MU;
return((x > x1)||(x < x1 - WALL_WIDTH));
}
}
case (D_TWO_PARABOLAS_ASYM_GUIDE):
{
if (x > 0.0)
{
if (vabs(y) > 0.5*LAMBDA) return(1);
x1 = -y*y/LAMBDA + 0.25*LAMBDA;
return((x < x1)||(x > x1 + WALL_WIDTH));
}
else
{
for (i=0; i<npolyrect; i++)
if ((x > polyrect[i].x1)&&(x < polyrect[i].x2)&&(y > polyrect[i].y1)&&(y < polyrect[i].y2)) return(0);
if (vabs(y) > 0.5*MU) return(1);
x1 = y*y/MU - 0.25*MU;
return((x > x1)||(x < x1 - WALL_WIDTH));
}
}
case (D_FOUR_PARABOLAS): case (D_FOUR_PARABOLAS):
{ {
x1 = MU + LAMBDA - 0.25*y*y/MU; x1 = MU + LAMBDA - 0.25*y*y/MU;
@ -6767,6 +6822,34 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_
} }
return(0); return(0);
} }
case (D_DISC_WAVEGUIDE):
{
if (x > 0.0) return(x*x + y*y < LAMBDA*LAMBDA);
y1 = y - 0.5*WALL_WIDTH;
r = LAMBDA - 0.5*WALL_WIDTH;
if (x*x + y1*y1 < r*r) return(1);
if ((y > -LAMBDA)&&(y < -LAMBDA + WALL_WIDTH)) return(1);
if ((y > -LAMBDA - 0.5*WALL_WIDTH)&&(y < -LAMBDA + 1.5*WALL_WIDTH)) return(2);
return(0);
}
case (D_DISC_WAVEGUIDE_SHIFTED):
{
if (x*x + y*y < LAMBDA*LAMBDA) return(1);
if (x > 0.0) return(0.0);
y1 = y + LAMBDA - MU - 0.5*WALL_WIDTH;
if (vabs(y1) < 0.5*WALL_WIDTH) return(1);
if (vabs(y1) < 0.75*WALL_WIDTH) return(2);
return(0);
}
case (D_DISC_WAVEGUIDE_ELLIPSE):
{
if (x*x/(LAMBDA*LAMBDA) + y*y < 1.0) return(1);
if (x > 0.0) return(0.0);
y1 = y + 1.0 - MU - 0.5*WALL_WIDTH;
if (vabs(y1) < 0.5*WALL_WIDTH) return(1);
if (vabs(y1) < 0.75*WALL_WIDTH) return(2);
return(0);
}
case (D_MENGER): case (D_MENGER):
{ {
x1 = 0.5*(x+1.0); x1 = 0.5*(x+1.0);
@ -7119,6 +7202,26 @@ void draw_billiard(int fade, double fade_value) /* draws the billiard bound
} }
break; break;
} }
case (D_CIRCLE_LAYERS):
{
for (j=1; j<=MDEPTH; j++)
{
r = (double)j*LAMBDA/(double)MDEPTH;
glBegin(GL_LINE_LOOP);
if (j<MDEPTH) glLineWidth(BOUNDARY_WIDTH/2);
else glLineWidth(BOUNDARY_WIDTH);
for (i=0; i<=NSEG; i++)
{
phi = (double)i*DPI/(double)NSEG;
x = r*cos(phi);
y = r*sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
}
break;
}
case (D_EXT_ELLIPSE): case (D_EXT_ELLIPSE):
{ {
glBegin(GL_LINE_LOOP); glBegin(GL_LINE_LOOP);
@ -7459,6 +7562,99 @@ void draw_billiard(int fade, double fade_value) /* draws the billiard bound
break; break;
} }
case (D_TWO_PARABOLAS_ASYM):
{
dy = LAMBDA/(double)NSEG;
glBegin(GL_LINE_LOOP);
for (i = 0; i < NSEG+1; i++)
{
y = -0.5*LAMBDA + dy*(double)i;
x = 0.25*LAMBDA - y*y/LAMBDA;
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
for (i = 0; i < NSEG+1; i++)
{
y = 0.5*LAMBDA - dy*(double)i;
x = 0.25*LAMBDA - y*y/LAMBDA + WALL_WIDTH;
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
dy = MU/(double)NSEG;
glBegin(GL_LINE_LOOP);
for (i = 0; i < NSEG+1; i++)
{
y = -0.5*MU + dy*(double)i;
x = -0.25*MU + y*y/MU;
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
for (i = 0; i < NSEG+1; i++)
{
y = 0.5*MU - dy*(double)i;
x = -0.25*MU + y*y/MU - WALL_WIDTH;
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
if (FOCI)
{
glColor3f(0.3, 0.3, 0.3);
draw_circle(0.0, 0.0, r, NSEG);
}
break;
}
case (D_TWO_PARABOLAS_ASYM_GUIDE):
{
dy = LAMBDA/(double)NSEG;
glBegin(GL_LINE_LOOP);
for (i = 0; i < NSEG+1; i++)
{
y = -0.5*LAMBDA + dy*(double)i;
x = 0.25*LAMBDA - y*y/LAMBDA;
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
for (i = 0; i < NSEG+1; i++)
{
y = 0.5*LAMBDA - dy*(double)i;
x = 0.25*LAMBDA - y*y/LAMBDA + WALL_WIDTH;
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
dy = MU/(double)NSEG;
glBegin(GL_LINE_LOOP);
for (i = 0; i < NSEG+1; i++)
{
y = -0.5*MU + dy*(double)i;
x = -0.25*MU + y*y/MU;
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
for (i = 0; i < NSEG+1; i++)
{
y = 0.5*MU - dy*(double)i;
x = -0.25*MU + y*y/MU - WALL_WIDTH;
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
for (i=0; i<npolyrect; i++)
draw_rectangle(polyrect[i].x1, polyrect[i].y1, polyrect[i].x2, polyrect[i].y2);
if (FOCI)
{
glColor3f(0.3, 0.3, 0.3);
draw_circle(0.0, 0.0, r, NSEG);
}
break;
}
case (D_FOUR_PARABOLAS): case (D_FOUR_PARABOLAS):
{ {
x1 = 2.0*(sqrt(MU*(2.0*MU + LAMBDA)) - MU); x1 = 2.0*(sqrt(MU*(2.0*MU + LAMBDA)) - MU);
@ -9271,6 +9467,46 @@ void draw_billiard(int fade, double fade_value) /* draws the billiard bound
break; break;
} }
case (D_DISC_WAVEGUIDE):
{
draw_circle_arc(0.0, 0.0, LAMBDA, -PID, PI, NSEG);
draw_circle_arc(0.0, 0.5*WALL_WIDTH, LAMBDA - 0.5*WALL_WIDTH, PID, PI, NSEG);
draw_line(XMIN, -LAMBDA, 0.0, -LAMBDA);
draw_line(XMIN, -LAMBDA+WALL_WIDTH, 0.0, -LAMBDA+WALL_WIDTH);
break;
}
case (D_DISC_WAVEGUIDE_SHIFTED):
{
alpha = acos((LAMBDA - MU)/LAMBDA);
alpha2 = acos((LAMBDA - MU - WALL_WIDTH)/LAMBDA);
draw_circle_arc(0.0, 0.0, LAMBDA, -PID - alpha, DPI - alpha2 + alpha, NSEG);
draw_line(XMIN, -LAMBDA + MU, -LAMBDA*sin(alpha), -LAMBDA + MU);
draw_line(XMIN, -LAMBDA + MU + WALL_WIDTH, -LAMBDA*sin(alpha2), -LAMBDA + MU + WALL_WIDTH);
break;
}
case (D_DISC_WAVEGUIDE_ELLIPSE):
{
alpha = acos(1.0 - MU);
alpha2 = acos(1.0 - MU - WALL_WIDTH);
draw_ellipse_arc(0.0, 0.0, LAMBDA, 1.0, -PID - alpha, DPI - alpha2 + alpha, NSEG);
draw_line(XMIN, -1.0 + MU, -LAMBDA*sin(alpha), -1.0 + MU);
draw_line(XMIN, -1.0 + MU + WALL_WIDTH, -LAMBDA*sin(alpha2), -1.0 + MU + WALL_WIDTH);
/* draw foci */
if (FOCI)
{
if (fade) glColor3f(0.3*fade_value, 0.3*fade_value, 0.3*fade_value);
else glColor3f(0.3, 0.3, 0.3);
x0 = sqrt(LAMBDA*LAMBDA-1.0);
glLineWidth(2);
glEnable(GL_LINE_SMOOTH);
draw_circle(x0, 0.0, r, NSEG);
draw_circle(-x0, 0.0, r, NSEG);
}
break;
}
case (D_MENGER): case (D_MENGER):
{ {
glLineWidth(3); glLineWidth(3);
@ -10621,7 +10857,7 @@ void compute_laplacian(double phi[NX*NY], t_laplacian laplace[NX*NY], double del
double oscillating_bc(int time, int j) double oscillating_bc(int time, int j)
{ {
int ij[2], jmin, jmax; int ij[2], jmin, jmax, jmid;
double t, t1, phase, a, envelope, omega, amp, dist2; double t, t1, phase, a, envelope, omega, amp, dist2;
switch (OSCILLATION_SCHEDULE) switch (OSCILLATION_SCHEDULE)
@ -10692,10 +10928,12 @@ double oscillating_bc(int time, int j)
} }
case (OSC_BEAM_SINE): case (OSC_BEAM_SINE):
{ {
xy_to_ij(0.0, OSCIL_YMAX, ij); xy_to_ij(0.0, OSCIL_YMID + OSCIL_YMAX, ij);
jmax = ij[1]; jmax = ij[1];
jmin = NY - jmax; xy_to_ij(0.0, OSCIL_YMID - OSCIL_YMAX, ij);
dist2 = (double)(j - NY/2)/(double)(jmax - NY/2); jmin = ij[1];
jmid = (jmin + jmax)/2;
dist2 = (double)(j - jmid)/(double)(jmax - jmid);
if (j > jmax) return(0.0); if (j > jmax) return(0.0);
if (j < jmin) return(0.0); if (j < jmin) return(0.0);
t = (double)time*OMEGA; t = (double)time*OMEGA;
@ -10703,6 +10941,22 @@ double oscillating_bc(int time, int j)
amp *= cos(PID*dist2); amp *= cos(PID*dist2);
return(amp); return(amp);
} }
case (OSC_BEAM_SINE_DECREASING):
{
xy_to_ij(0.0, OSCIL_YMID + OSCIL_YMAX, ij);
jmax = ij[1];
xy_to_ij(0.0, OSCIL_YMID - OSCIL_YMAX, ij);
jmin = ij[1];
jmid = (jmin + jmax)/2;
dist2 = (double)(j - jmid)/(double)(jmax - jmid);
if (j > jmax) return(0.0);
if (j < jmin) return(0.0);
t = (double)time*OMEGA;
a = INITIAL_VARIANCE/(OMEGA*OMEGA);
amp = AMPLITUDE*cos(t)*exp(-(double)(t-INITIAL_SHIFT)*(t-INITIAL_SHIFT)/a);
amp *= cos(PID*dist2);
return(amp);
}
case (OSC_BEAM_TWOPERIODS): case (OSC_BEAM_TWOPERIODS):
{ {
xy_to_ij(0.0, OSCIL_YMAX, ij); xy_to_ij(0.0, OSCIL_YMAX, ij);
@ -10717,6 +10971,38 @@ double oscillating_bc(int time, int j)
amp *= exp(-dist2/INITIAL_VARIANCE); amp *= exp(-dist2/INITIAL_VARIANCE);
return(amp); return(amp);
} }
case (OSC_BEAM_SINE_TWOPERIODS):
{
xy_to_ij(0.0, OSCIL_YMID + OSCIL_YMAX, ij);
jmax = ij[1];
xy_to_ij(0.0, OSCIL_YMID - OSCIL_YMAX, ij);
jmin = ij[1];
jmid = (jmin + jmax)/2;
dist2 = (double)(j - jmid)/(double)(jmax - jmid);
if (j > jmax) return(0.0);
if (j < jmin) return(0.0);
t = (double)time*OMEGA;
amp = AMPLITUDE*(0.5*cos(t) + 0.5*cos(sqrt(7.0)*t))*exp(-(double)t*DAMPING);
amp *= cos(PID*dist2);
return(amp);
}
case (OSC_BEAM_SINE_CHIRP):
{
xy_to_ij(0.0, OSCIL_YMID + OSCIL_YMAX, ij);
jmax = ij[1];
xy_to_ij(0.0, OSCIL_YMID - OSCIL_YMAX, ij);
jmin = ij[1];
jmid = (jmin + jmax)/2;
dist2 = (double)(j - jmid)/(double)(jmax - jmid);
if (j > jmax) return(0.0);
if (j < jmin) return(0.0);
t = (double)time*OMEGA;
phase = t + ACHIRP*t*t;
// if (phase > DPI) phase = DPI;
amp = AMPLITUDE*sin(phase)*exp(-phase*DAMPING);
amp *= cos(PID*dist2);
return(amp);
}
case (OSC_TWO_WAVES): case (OSC_TWO_WAVES):
{ {
t = (double)time*OMEGA + (double)(NY-j)*OSCIL_LEFT_YSHIFT/(double)NY; t = (double)time*OMEGA + (double)(NY-j)*OSCIL_LEFT_YSHIFT/(double)NY;
@ -10737,13 +11023,13 @@ double oscillating_bc(int time, int j)
} }
} }
void init_ior_2d(short int *xy_in[NX], double *tcc_table[NX], double *tgamma_table[NX], double ior_angle) void init_ior_2d(short int *xy_in[NX], double *tcc_table[NX], double *tc_table[NX], double *tgamma_table[NX], double ior_angle)
/* compute variable index of refraction */ /* compute variable index of refraction */
/* should be at some point merged with 3D version in suv_wave_3d.c */ /* should be at some point merged with 3D version in suv_wave_3d.c */
{ {
int i, j, k, n, inlens, ncircles; int i, j, k, n, inlens, ncircles;
double courant2 = COURANT*COURANT, courantb2 = COURANTB*COURANTB, lambda1, mu1, courant_med, courant_med2; double courant2 = COURANT*COURANT, courantb2 = COURANTB*COURANTB, lambda1, mu1, courant_med, courant_med2;
double a, u, v, u1, x, y, xy[2], norm2, speed, r2, c, salpha, h, ll, ca, sa, x1, y1, dx, dy, sum, sigma, x0, y0, rgb[3]; double a, u, v, u1, x, y, xy[2], norm2, speed, r, r2, c, salpha, h, ll, ca, sa, x1, y1, dx, dy, sum, sigma, x0, y0, rgb[3];
double xc[NGRIDX*NGRIDY], yc[NGRIDX*NGRIDY], height[NGRIDX*NGRIDY]; double xc[NGRIDX*NGRIDY], yc[NGRIDX*NGRIDY], height[NGRIDX*NGRIDY];
static double xc_stat[NGRIDX*NGRIDY], yc_stat[NGRIDX*NGRIDY], sigma_stat; static double xc_stat[NGRIDX*NGRIDY], yc_stat[NGRIDX*NGRIDY], sigma_stat;
static int first = 1; static int first = 1;
@ -11454,6 +11740,54 @@ void init_ior_2d(short int *xy_in[NX], double *tcc_table[NX], double *tgamma_tab
} }
break; break;
} }
case (IOR_LUNEBURG_LENS):
{
/* n = sqrt(2-r/R) */
printf("Initialising Luneburg lens IOR");
for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
ij_to_xy(i, j, xy);
r = module2(xy[0], xy[1]);
if (r > LAMBDA)
{
tcc_table[i][j] = courant2;
tgamma_table[i][j] = GAMMA;
}
else
{
tcc_table[i][j] = courant2/(2.0 - r/LAMBDA);
tgamma_table[i][j] = GAMMAB;
// printf("tcc[%i][%i] = %.3lg\n", i, j, tcc_table[i][j]);
}
}
}
break;
}
case (IOR_LUNEBURG_LAYERS):
{
/* n = sqrt(2-r/R) discretized in layers */
printf("Initialising Luneburg lens IOR");
for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
ij_to_xy(i, j, xy);
r = module2(xy[0], xy[1])/LAMBDA;
r2 = r*(double)MDEPTH;
r2 = (double)((int)r2)/(double)MDEPTH;
if (r > 1.0)
{
tcc_table[i][j] = courant2;
tgamma_table[i][j] = GAMMA;
}
else
{
tcc_table[i][j] = courant2/(2.0 - r2);
tgamma_table[i][j] = GAMMAB;
// printf("tcc[%i][%i] = %.3lg\n", i, j, tcc_table[i][j]);
}
}
}
break;
}
case (IOR_LINEAR_X_A): case (IOR_LINEAR_X_A):
{ {
/* Warning: Depending on COURANT and COURANTB */ /* Warning: Depending on COURANT and COURANTB */
@ -11495,6 +11829,18 @@ void init_ior_2d(short int *xy_in[NX], double *tcc_table[NX], double *tgamma_tab
} }
} }
} }
for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
if (xy_in[i][j] != 0)
{
tc_table[i][j] = sqrt(tcc_table[i][j]);
}
else if (TWOSPEEDS)
{
tc_table[i][j] = COURANTB;
}
}
}
} }
else else
{ {
@ -11503,14 +11849,14 @@ void init_ior_2d(short int *xy_in[NX], double *tcc_table[NX], double *tgamma_tab
for (j=0; j<NY; j++){ for (j=0; j<NY; j++){
if (xy_in[i][j] != 0) if (xy_in[i][j] != 0)
{ {
// tc[i*NY+j] = COURANT; tc_table[i][j] = COURANT;
tcc_table[i][j] = courant2; tcc_table[i][j] = courant2;
if (xy_in[i][j] == 1) tgamma_table[i][j] = GAMMA; if (xy_in[i][j] == 1) tgamma_table[i][j] = GAMMA;
else tgamma_table[i][j] = GAMMAB; else tgamma_table[i][j] = GAMMAB;
} }
else if (TWOSPEEDS) else if (TWOSPEEDS)
{ {
// tc[i*NY+j] = COURANTB; tc_table[i][j] = COURANTB;
tcc_table[i][j] = courantb2; tcc_table[i][j] = courantb2;
tgamma_table[i][j] = GAMMAB; tgamma_table[i][j] = GAMMAB;
} }

32
sub_wave_3d.c
View File

@ -1242,7 +1242,8 @@ double compute_interpolated_colors_wave(int i, int j, short int xy_in[NX*NY], t_
ca = wave[i*NY+j].cos_angle; ca = wave[i*NY+j].cos_angle;
ca = (ca + 1.0)*0.4 + 0.2; ca = (ca + 1.0)*0.4 + 0.2;
// if ((FADE_IN_OBSTACLE)&&(!xy_in[i*NY+j])) ca *= 1.6; // if ((FADE_IN_OBSTACLE)&&(!xy_in[i*NY+j])) ca *= 1.6;
if ((FADE_IN_OBSTACLE)&&(!xy_in[i*NY+j])) ca = (ca + 0.1)*1.6; if ((FADE_IN_OBSTACLE == 1)&&(!xy_in[i*NY+j])) ca = (ca + 0.1)*1.6;
if ((FADE_IN_OBSTACLE == 2)&&(xy_in[i*NY+j])) ca = (ca + 0.1)*1.6;
for (k=0; k<3; k++) for (k=0; k<3; k++)
{ {
rgb_e[k] *= ca; rgb_e[k] *= ca;
@ -1296,7 +1297,7 @@ void init_speed_dissipation(short int xy_in[NX*NY], double tc[NX*NY], double tcc
{ {
int i, j, k, n, inlens; int i, j, k, n, inlens;
double courant2 = COURANT*COURANT, courantb2 = COURANTB*COURANTB, lambda1, mu1; double courant2 = COURANT*COURANT, courantb2 = COURANTB*COURANTB, lambda1, mu1;
double u, v, u1, x, y, xy[2], norm2, speed, r2, c, salpha, h, ll, ca, sa, x1, y1, dx, dy, sum, sigma, x0, y0, rgb[3]; double u, v, u1, x, y, xy[2], norm2, speed, r2, c, salpha, h, ll, ca, sa, x1, y1, dx, dy, sum, sigma, x0, y0, rgb[3], r;
double xc[NGRIDX*NGRIDY], yc[NGRIDX*NGRIDY], height[NGRIDX*NGRIDY]; double xc[NGRIDX*NGRIDY], yc[NGRIDX*NGRIDY], height[NGRIDX*NGRIDY];
t_circle circles[NMAXCIRCLES]; t_circle circles[NMAXCIRCLES];
@ -1728,6 +1729,30 @@ void init_speed_dissipation(short int xy_in[NX*NY], double tc[NX*NY], double tcc
} }
break; break;
} }
case (IOR_LUNEBURG_LENS):
{
/* n = sqrt(2-r/R) */
printf("Initialising Luneburg lens IOR");
for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
ij_to_xy(i, j, xy);
r = module2(xy[0], xy[1]);
tc[i*NY+j] = COURANT;
if (r > LAMBDA)
{
tcc[i*NY+j] = courant2;
tgamma[i*NY+j] = GAMMA;
}
else
{
tcc[i*NY+j] = courant2/(2.0 - r/LAMBDA);
tgamma[i*NY+j] = GAMMAB;
// printf("tcc[%i][%i] = %.3lg\n", i, j, tcc_table[i][j]);
}
}
}
break;
}
default: default:
{ {
for (i=0; i<NX; i++){ for (i=0; i<NX; i++){
@ -1930,7 +1955,8 @@ void compute_cfield(short int xy_in[NX*NY], int cplot, int palette, t_wave wave[
ca = wave[i*NY+j].cos_angle; ca = wave[i*NY+j].cos_angle;
ca = (ca + 1.0)*0.4 + 0.2; ca = (ca + 1.0)*0.4 + 0.2;
// if ((FADE_IN_OBSTACLE)&&(!xy_in[i*NY+j])) ca *= 1.6; // if ((FADE_IN_OBSTACLE)&&(!xy_in[i*NY+j])) ca *= 1.6;
if ((FADE_IN_OBSTACLE)&&(!xy_in[i*NY+j])) ca = (ca + 0.1)*1.6; if ((FADE_IN_OBSTACLE == 1)&&(!xy_in[i*NY+j])) ca = (ca + 0.1)*1.6;
else if ((FADE_IN_OBSTACLE == 2)&&(xy_in[i*NY+j])) ca = (ca + 0.1)*1.6;
for (k=0; k<3; k++) wave[i*NY+j].rgb[k] *= ca; 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; if (fade) for (k=0; k<3; k++) wave[i*NY+j].rgb[k] *= fade_value;

182
sub_wave_comp.c
View File

@ -489,6 +489,18 @@ int xy_in_billiard_half(double x, double y, int domain, int pattern, int top)
return(1); return(1);
break; break;
} }
case (D_ELLIPSE):
{
if (x*x/(LAMBDA*LAMBDA) + y*y < 1.0) return(1);
else return(0);
break;
}
case (D_CIRCLE_LAYERS):
{
if (x*x + y*y < LAMBDA*LAMBDA) return(1);
else return(0);
break;
}
case (D_MENGER): case (D_MENGER):
{ {
x1 = 0.5*(x+1.0); x1 = 0.5*(x+1.0);
@ -752,7 +764,7 @@ void draw_billiard_half(int domain, int pattern, int top, int fade, double fade_
/* two domain version implemented for D_CIRCLES */ /* two domain version implemented for D_CIRCLES */
{ {
double x0, x, y, x1, y1, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z, l, signtop, width, a, b, arcangle; double x0, x, y, x1, y1, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z, l, signtop, width, a, b, arcangle;
int i, j, k, k1, k2, mr2; int i, j, k, k1, k2, mr2, mdepth;
static int first = 1; static int first = 1;
glEnable(GL_SCISSOR_TEST); glEnable(GL_SCISSOR_TEST);
@ -776,12 +788,68 @@ void draw_billiard_half(int domain, int pattern, int top, int fade, double fade_
} }
glLineWidth(BOUNDARY_WIDTH); glLineWidth(BOUNDARY_WIDTH);
if (top) signtop = 1.0; if (top)
else signtop = -1.0; {
signtop = 1.0;
mdepth = MDEPTH;
}
else
{
signtop = -1.0;
mdepth = MDEPTH_B;
}
glEnable(GL_LINE_SMOOTH); glEnable(GL_LINE_SMOOTH);
switch (domain) { switch (domain) {
case (D_ELLIPSE):
{
glBegin(GL_LINE_LOOP);
for (i=0; i<=NSEG; i++)
{
phi = (double)i*DPI/(double)NSEG;
x = LAMBDA*cos(phi);
y = sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
/* draw foci */
if (FOCI)
{
if (fade) glColor3f(0.3*fade_value, 0.3*fade_value, 0.3*fade_value);
else glColor3f(0.3, 0.3, 0.3);
x0 = sqrt(LAMBDA*LAMBDA-1.0);
glLineWidth(2);
glEnable(GL_LINE_SMOOTH);
draw_circle(x0, 0.0, r, NSEG);
draw_circle(-x0, 0.0, r, NSEG);
}
break;
}
case (D_CIRCLE_LAYERS):
{
for (j=1; j<=mdepth; j++)
{
r = (double)j*LAMBDA/(double)mdepth;
glBegin(GL_LINE_LOOP);
if (j<mdepth) glLineWidth(BOUNDARY_WIDTH/2);
else glLineWidth(BOUNDARY_WIDTH);
for (i=0; i<=NSEG; i++)
{
phi = (double)i*DPI/(double)NSEG;
x = r*cos(phi);
y = r*sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
}
break;
}
case (D_MENGER): case (D_MENGER):
{ {
glLineWidth(3); glLineWidth(3);
@ -1615,13 +1683,13 @@ void print_energies(double energies[6], double top_energy, double bottom_energy)
} }
void init_ior_2d_comp(short int *xy_in[NX], double *tcc_table[NX], double ior_angle) void init_ior_2d_comp_half(short int *xy_in[NX], double *tcc_table[NX], double ior_angle, int ior, int mdepth, int top)
/* compute variable index of refraction */ /* compute variable index of refraction */
/* should be at some point merged with 3D version in suv_wave_3d.c */ /* should be at some point merged with 3D version in suv_wave_3d.c */
{ {
int i, j, k, n, inlens, ncircles; int i, j, k, n, inlens, ncircles, jmin, jmax;
double courant2 = COURANT*COURANT, courantb2 = COURANTB*COURANTB, lambda1, mu1; double courant2 = COURANT*COURANT, courantb2 = COURANTB*COURANTB, lambda1, mu, mu1;
double u, v, u1, x, y, xy[2], norm2, speed, r2, c, salpha, h, ll, ca, sa, x1, y1, dx, dy, sum, sigma, x0, y0, rgb[3]; double u, v, u1, x, y, xy[2], norm2, speed, r2, c, salpha, h, ll, ca, sa, x1, y1, dx, dy, sum, sigma, x0, y0, r, rgb[3];
double xc[NGRIDX*NGRIDY], yc[NGRIDX*NGRIDY], height[NGRIDX*NGRIDY]; double xc[NGRIDX*NGRIDY], yc[NGRIDX*NGRIDY], height[NGRIDX*NGRIDY];
static double xc_stat[NGRIDX*NGRIDY], yc_stat[NGRIDX*NGRIDY], sigma_stat; static double xc_stat[NGRIDX*NGRIDY], yc_stat[NGRIDX*NGRIDY], sigma_stat;
static int first = 1; static int first = 1;
@ -1631,29 +1699,31 @@ void init_ior_2d_comp(short int *xy_in[NX], double *tcc_table[NX], double ior_an
rgb[1] = 1.0; rgb[1] = 1.0;
rgb[2] = 1.0; rgb[2] = 1.0;
if (top)
{
jmin = NY/2;
jmax = NY;
mu = MU;
}
else
{
jmin = 0;
jmax = NY/2;
mu = MUB;
}
if (VARIABLE_IOR) if (VARIABLE_IOR)
{ {
switch (IOR) { switch (ior) {
case (IOR_LENS_WALL): case (IOR_LENS_WALL):
{ {
printf("Initializing IOR_LENS_WALL\n"); printf("Initializing IOR_LENS_WALL\n");
for (i=0; i<NX; i++){ for (i=0; i<NX; i++){
for (j=0; j<NY/2; j++){ for (j=jmin; j<jmax; j++){
ij_to_xy(i, j, xy); ij_to_xy(i, j, xy);
x = xy[0]; x = xy[0];
y = xy[1]; y = xy[1];
if ((vabs(x) < 0.05)&&(vabs(y) > MUB)) tcc_table[i][j] = 0.0; if ((vabs(x) < 0.05)&&(vabs(y) > mu)) tcc_table[i][j] = 0.0;
else
{
if (xy_in[i][j] != 0) tcc_table[i][j] = courant2;
else tcc_table[i][j] = courantb2;
}
}
for (j=NY/2; j<NY; j++){
ij_to_xy(i, j, xy);
x = xy[0];
y = xy[1];
if ((vabs(x) < 0.05)&&(vabs(y) > MU)) tcc_table[i][j] = 0.0;
else else
{ {
if (xy_in[i][j] != 0) tcc_table[i][j] = courant2; if (xy_in[i][j] != 0) tcc_table[i][j] = courant2;
@ -1664,26 +1734,63 @@ void init_ior_2d_comp(short int *xy_in[NX], double *tcc_table[NX], double ior_an
break; break;
} }
case (IOR_LUNEBURG_LENS):
{
/* n = sqrt(2-r/R) */
printf("Initialising Luneburg lens IOR");
for (i=0; i<NX; i++){
for (j=jmin; j<jmax; j++){
ij_to_xy(i, j, xy);
r = module2(xy[0], xy[1]);
if (r > 1.0)
{
tcc_table[i][j] = courant2;
// tgamma_table[i][j] = GAMMA;
}
else
{
tcc_table[i][j] = courant2/(2.0 - r/LAMBDA);
// tgamma_table[i][j] = GAMMAB;
// printf("tcc[%i][%i] = %.3lg\n", i, j, tcc_table[i][j]);
}
}
}
break;
}
case (IOR_LUNEBURG_LAYERS):
{
/* n = sqrt(2-r/R) discretized in layers */
printf("Initialising Luneburg lens IOR");
for (i=0; i<NX; i++){
for (j=jmin; j<jmax; j++){
ij_to_xy(i, j, xy);
r = module2(xy[0], xy[1])/LAMBDA;
r2 = r*(double)mdepth;
r2 = (double)((int)r2)/(double)mdepth;
if (r > 1.0)
{
tcc_table[i][j] = courant2;
// tgamma_table[i][j] = GAMMA;
}
else
{
tcc_table[i][j] = courant2/(2.0 - r2);
// tgamma_table[i][j] = GAMMAB;
// printf("tcc[%i][%i] = %.3lg\n", i, j, tcc_table[i][j]);
}
}
}
break;
}
case (IOR_LENS_OBSTACLE): case (IOR_LENS_OBSTACLE):
{ {
printf("Initializing IOR_LENS_OBSTACLE\n"); printf("Initializing IOR_LENS_OBSTACLE\n");
for (i=0; i<NX; i++){ for (i=0; i<NX; i++){
for (j=0; j<NY/2; j++){ for (j=jmin; j<jmax; j++){
ij_to_xy(i, j, xy); ij_to_xy(i, j, xy);
x = xy[0]; x = xy[0];
y = xy[1]; y = xy[1];
if ((vabs(x) < 0.05)&&(vabs(y) < MUB)) tcc_table[i][j] = 0.0; if ((vabs(x) < 0.05)&&(vabs(y) < mu)) tcc_table[i][j] = 0.0;
else
{
if (xy_in[i][j] != 0) tcc_table[i][j] = courant2;
else tcc_table[i][j] = courantb2;
}
}
for (j=NY/2; j<NY; j++){
ij_to_xy(i, j, xy);
x = xy[0];
y = xy[1];
if ((vabs(x) < 0.05)&&(vabs(y) < MU)) tcc_table[i][j] = 0.0;
else else
{ {
if (xy_in[i][j] != 0) tcc_table[i][j] = courant2; if (xy_in[i][j] != 0) tcc_table[i][j] = courant2;
@ -1721,3 +1828,10 @@ void init_ior_2d_comp(short int *xy_in[NX], double *tcc_table[NX], double ior_an
} }
} }
} }
void init_ior_2d_comp(short int *xy_in[NX], double *tcc_table[NX], double ior_angle)
/* compute variable index of refraction */
{
init_ior_2d_comp_half(xy_in, tcc_table, ior_angle, IOR, MDEPTH, 1);
init_ior_2d_comp_half(xy_in, tcc_table, ior_angle, IOR_B, MDEPTH_B, 0);
}

130
wave_3d.c
View File

@ -44,7 +44,7 @@
#include <time.h> #include <time.h>
#define MOVIE 0 /* set to 1 to generate movie */ #define MOVIE 0 /* set to 1 to generate movie */
#define DOUBLE_MOVIE 0 /* set to 1 to produce movies for wave height and energy simultaneously */ #define 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 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 NO_EXTRA_BUFFER_SWAP 1 /* some OS require one less buffer swap when recording images */
@ -52,8 +52,12 @@
#define WINWIDTH 1920 /* window width */ #define WINWIDTH 1920 /* window width */
#define WINHEIGHT 1150 /* window height */ #define WINHEIGHT 1150 /* window height */
#define NX 1920 /* number of grid points on x axis */ // #define NX 1920 /* number of grid points on x axis */
#define NY 1150 /* number of grid points on y axis */ // #define NY 1150 /* number of grid points on y axis */
#define NX 2880 /* number of grid points on x axis */
#define NY 1725 /* number of grid points on y axis */
// #define NX 2400 /* number of grid points on x axis */
// #define NY 1400 /* number of grid points on y axis */
#define XMIN -2.0 #define XMIN -2.0
#define XMAX 2.0 /* x interval */ #define XMAX 2.0 /* x interval */
@ -66,7 +70,7 @@
/* Choice of the billiard table */ /* Choice of the billiard table */
#define B_DOMAIN 97 /* choice of domain shape, see list in global_pdes.c */ #define B_DOMAIN 192 /* choice of domain shape, see list in global_pdes.c */
#define CIRCLE_PATTERN 2 /* pattern of circles or polygons, see list in global_pdes.c */ #define CIRCLE_PATTERN 2 /* pattern of circles or polygons, see list in global_pdes.c */
@ -76,7 +80,7 @@
#define IMAGE_FILE 5 /* for option D_IMAGE */ #define IMAGE_FILE 5 /* for option D_IMAGE */
#define VARIABLE_IOR 0 /* set to 1 for a variable index of refraction */ #define VARIABLE_IOR 0 /* set to 1 for a variable index of refraction */
#define IOR 181 /* choice of index of refraction, see list in global_pdes.c */ #define IOR 182 /* 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 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 MANDEL_IOR_SCALE -0.05 /* parameter controlling dependence of IoR on Mandelbrot escape speed */
@ -86,11 +90,11 @@
#define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */ #define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */
#define PDISC_CONNECT_FACTOR 1.5 /* controls which discs are connected for D_CIRCLE_LATTICE_POISSON domain */ #define PDISC_CONNECT_FACTOR 1.5 /* controls which discs are connected for D_CIRCLE_LATTICE_POISSON domain */
#define LAMBDA 0.6 /* parameter controlling the dimensions of domain */ #define LAMBDA 2.3 /* parameter controlling the dimensions of domain */
#define MU 0.11 /* parameter controlling the dimensions of domain */ #define MU 1.15 /* parameter controlling the dimensions of domain */
#define MU_B 1.0 /* parameter controlling the dimensions of domain */ #define MU_B 1.0 /* parameter controlling the dimensions of domain */
#define NPOLY 6 /* number of sides of polygon */ #define NPOLY 3 /* number of sides of polygon */
#define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */ #define APOLY -1.0 /* angle by which to turn polygon, in units of Pi/2 */
#define MDEPTH 6 /* depth of computation of Menger gasket */ #define MDEPTH 6 /* depth of computation of Menger gasket */
#define MRATIO 3 /* ratio defining Menger gasket */ #define MRATIO 3 /* ratio defining Menger gasket */
#define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ #define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */
@ -98,8 +102,8 @@
#define FOCI 1 /* set to 1 to draw focal points of ellipse */ #define FOCI 1 /* set to 1 to draw focal points of ellipse */
#define NGRIDX 12 /* number of grid point for grid of disks */ #define NGRIDX 12 /* number of grid point for grid of disks */
#define NGRIDY 8 /* number of grid point for grid of disks */ #define NGRIDY 8 /* number of grid point for grid of disks */
#define WALL_WIDTH 0.022 /* width of wall separating lenses */ #define WALL_WIDTH 0.1 /* width of wall separating lenses */
#define WALL_WIDTH_B 0.01 /* width of wall separating lenses */ #define WALL_WIDTH_B 0.05 /* width of wall separating lenses */
#define WALL_WIDTH_RND 0.0 /* proportion of width of width for random arrangements */ #define WALL_WIDTH_RND 0.0 /* proportion of width of width for random arrangements */
#define RADIUS_FACTOR 0.3 /* controls inner radius for C_RING arrangements */ #define RADIUS_FACTOR 0.3 /* controls inner radius for C_RING arrangements */
#define WALL_WIDTH_ASYM 0.75 /* asymmetry of wall width (D_CIRCLE_LATTICE_NONISO) */ #define WALL_WIDTH_ASYM 0.75 /* asymmetry of wall width (D_CIRCLE_LATTICE_NONISO) */
@ -123,19 +127,20 @@
/* Physical parameters of wave equation */ /* Physical parameters of wave equation */
#define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */ #define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */
#define 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 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 OSCILLATION_SCHEDULE 42 /* oscillation schedule, see list in global_pdes.c */
#define OSCIL_YMAX 0.35 /* defines oscillation range */ #define OSCIL_YMAX 0.1 /* defines oscilling beam range */
#define OSCIL_YMID -0.75 /* defines oscilling beam midpoint */
#define INITIAL_SHIFT 20.0 /* time shift of initial wave packet (in oscillation periods) */ #define INITIAL_SHIFT 20.0 /* time shift of initial wave packet (in oscillation periods) */
#define WAVE_PACKET_SHIFT 200.0 /* time shift between wave packets (in oscillation periods) */ #define WAVE_PACKET_SHIFT 200.0 /* time shift between wave packets (in oscillation periods) */
#define OMEGA 0.025 /* frequency of periodic excitation */ #define OMEGA 0.02 /* frequency of periodic excitation */
#define AMPLITUDE 0.5 /* amplitude of periodic excitation */ #define AMPLITUDE 1.0 /* amplitude of periodic excitation */
#define ACHIRP 0.2 /* acceleration coefficient in chirp */ #define ACHIRP 0.2 /* acceleration coefficient in chirp */
#define DAMPING 0.0 /* damping of periodic excitation */ #define DAMPING 0.0 /* damping of periodic excitation */
#define COURANT 0.1 /* Courant number */ #define COURANT 0.15 /* Courant number in medium B */
#define COURANTB 0.08 /* Courant number in medium B */ #define COURANTB 0.1 /* Courant number */
#define GAMMA 0.0 /* damping factor in wave equation */ #define GAMMA 0.0 /* damping factor in wave equation */
#define GAMMAB 1.0e-6 /* 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_SIDES 1.0e-4 /* damping factor on boundary */
@ -149,7 +154,7 @@
/* Increasing COURANT speeds up the simulation, but decreases accuracy */ /* Increasing COURANT speeds up the simulation, but decreases accuracy */
/* For similar wave forms, COURANT^2*GAMMA should be kept constant */ /* 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 18 /* period of oscillating source */ #define OSCILLATING_SOURCE_PERIOD 18 /* period of oscillating source */
#define ALTERNATE_OSCILLATING_SOURCE 1 /* set to 1 to alternate sign of oscillating source */ #define ALTERNATE_OSCILLATING_SOURCE 1 /* set to 1 to alternate sign of oscillating source */
#define MAX_PULSING_TIME 1500 /* max time for adding pulses */ #define MAX_PULSING_TIME 1500 /* max time for adding pulses */
@ -161,16 +166,16 @@
/* Boundary conditions, see list in global_pdes.c */ /* Boundary conditions, see list in global_pdes.c */
#define B_COND 1 #define B_COND 2
#define PRECOMPUTE_BC 0 /* set to 1 to compute neighbours for Laplacian in advance */ #define PRECOMPUTE_BC 0 /* set to 1 to compute neighbours for Laplacian in advance */
/* Parameters for length and speed of simulation */ /* Parameters for length and speed of simulation */
#define NSTEPS 2800 /* number of frames of movie */ #define NSTEPS 2000 /* number of frames of movie */
#define NVID 3 /* number of iterations between images displayed on screen */ #define NVID 10 /* number of iterations between images displayed on screen */
#define NSEG 1000 /* number of segments of boundary */ #define NSEG 200 /* number of segments of boundary */
#define INITIAL_TIME 0 /* time after which to start saving frames */ #define INITIAL_TIME 100 /* time after which to start saving frames */
#define BOUNDARY_WIDTH 2 /* width of billiard boundary */ #define BOUNDARY_WIDTH 2 /* width of billiard boundary */
#define PRINT_SPEED 0 /* set to 1 to print speed of moving source */ #define PRINT_SPEED 0 /* set to 1 to print speed of moving source */
@ -201,12 +206,12 @@
#define FLUX_WINDOW 30 /* size of averaging window of 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 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 1 /* set to 1 to change luminosity according to normal vector */
#define NON_DIRICHLET_BC 1 /* set to 1 to draw only facets in domain, if field is not zero on boundary */ #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 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 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_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 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 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_ENERGY 0.4 /* vertical scaling in energy plot */
@ -219,28 +224,28 @@
#define REP_AXO_3D 0 /* linear projection (axonometry) */ #define REP_AXO_3D 0 /* linear projection (axonometry) */
#define REP_PROJ_3D 1 /* projection on plane orthogonal to observer line of sight */ #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_VIEW 0 /* set to 1 to rotate position of observer */
#define ROTATE_ANGLE 360.0 /* total angle of rotation during simulation */ #define ROTATE_ANGLE 360.0 /* total angle of rotation during simulation */
/* Color schemes */ /* Color schemes */
#define COLOR_PALETTE 10 /* Color palette, see list in global_pdes.c */ #define COLOR_PALETTE 11 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 17 /* Color palette, see list in global_pdes.c */ #define COLOR_PALETTE_B 13 /* Color palette, see list in global_pdes.c */
#define BLACK 1 /* background */ #define BLACK 1 /* background */
#define COLOR_SCHEME 3 /* choice of color scheme, see list in global_pdes.c */ #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 SCALE 0 /* set to 1 to adjust color scheme to variance of field */
#define SLOPE 0.75 /* sensitivity of color on wave amplitude */ #define SLOPE 1.0 /* sensitivity of color on wave amplitude */
#define COLOR_RANGE 1.0 /* max range of color (default: 1.0) */ #define COLOR_RANGE 1.0 /* max range of color (default: 1.0) */
#define VSCALE_AMPLITUDE 0.5 /* additional scaling factor for color scheme P_3D_AMPLITUDE */ #define VSCALE_AMPLITUDE 1.0 /* additional scaling factor for color scheme P_3D_AMPLITUDE */
#define VSHIFT_AMPLITUDE 0.0 /* additional shift for wave amplitude */ #define VSHIFT_AMPLITUDE 0.0 /* additional shift for wave amplitude */
#define VSCALE_ENERGY 1.0 /* additional scaling factor for color scheme P_3D_ENERGY */ #define VSCALE_ENERGY 2.5 /* 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_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 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 ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */
#define E_SCALE 30.0 /* scaling factor for energy representation */ #define E_SCALE 100.0 /* scaling factor for energy representation */
#define LOG_SCALE 1.5 /* scaling factor for energy log representation */ #define LOG_SCALE 1.5 /* scaling factor for energy log representation */
#define LOG_SHIFT 0.25 /* shift of colors on log scale */ #define LOG_SHIFT 0.25 /* shift of colors on log scale */
#define LOG_ENERGY_FLOOR -10.0 /* floor value for log of (total) energy */ #define LOG_ENERGY_FLOOR -10.0 /* floor value for log of (total) energy */
@ -272,8 +277,8 @@
#define MAZE_WIDTH 0.02 /* half width of maze walls */ #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 1 /* set to 1 to plot the color scheme */
#define COLORBAR_RANGE 2.5 /* scale of color scheme bar */ #define COLORBAR_RANGE 3.0 /* scale of color scheme bar */
#define COLORBAR_RANGE_B 1.0 /* scale of color scheme bar for 2nd part */ #define COLORBAR_RANGE_B 1.5 /* scale of color scheme bar for 2nd part */
#define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ #define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */
#define SAVE_TIME_SERIES 0 /* set to 1 to save wave time series at a point */ #define SAVE_TIME_SERIES 0 /* set to 1 to save wave time series at a point */
@ -310,15 +315,16 @@ double u_3d[2] = {0.75, -0.45}; /* projections of basis vectors for REP_AXO_
double v_3d[2] = {-0.75, -0.45}; double v_3d[2] = {-0.75, -0.45};
double w_3d[2] = {0.0, 0.015}; 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] = {-8.0, -6.0, 6.0}; /* location of observer for REP_PROJ_3D representation */ double observer[3] = {8.0, 6.0, 6.0}; /* location of observer for REP_PROJ_3D representation */
int reset_view = 0; /* switch to reset 3D view parameters (for option ROTATE_VIEW) */ int reset_view = 0; /* switch to reset 3D view parameters (for option ROTATE_VIEW) */
#define Z_SCALING_FACTOR 0.05 /* overall scaling factor of z axis for REP_PROJ_3D representation */ #define Z_SCALING_FACTOR 0.1 /* overall scaling factor of z axis for REP_PROJ_3D representation */
#define XY_SCALING_FACTOR 1.8 /* overall scaling factor for on-screen (x,y) coordinates after projection */ #define XY_SCALING_FACTOR 2.3 /* 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 ZMAX_FACTOR 1.0 /* max value of z coordinate for REP_PROJ_3D representation */
#define XSHIFT_3D -0.1 /* overall x shift for REP_PROJ_3D representation */ #define XSHIFT_3D 0.0 /* overall x shift for REP_PROJ_3D representation */
#define YSHIFT_3D 0.2 /* overall y shift for REP_PROJ_3D representation */ #define YSHIFT_3D -0.25 /* overall y shift for REP_PROJ_3D representation */
#define BEAM_BC ((OSCILLATION_SCHEDULE == OSC_BEAM_SINE)||(OSCILLATION_SCHEDULE == OSC_BEAM_SINE_DECREASING)||(OSCILLATION_SCHEDULE == OSC_BEAM_SINE_TWOPERIODS)||(OSCILLATION_SCHEDULE == OSC_BEAM_SINE_CHIRP))
#include "global_pdes.c" /* constants and global variables */ #include "global_pdes.c" /* constants and global variables */
#include "global_3d.c" /* constants and global variables */ #include "global_3d.c" /* constants and global variables */
@ -340,10 +346,34 @@ void evolve_wave_half_new(double phi_in[NX*NY], double psi_in[NX*NY], double phi
/* phi is value of field at time t, psi at time t-1 */ /* 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 */ /* this version of the function has been rewritten in order to minimize the number of if-branches */
{ {
int i, j, iplus, iminus, jplus, jminus; int i, j, iplus, iminus, jplus, jminus, ij[2];
double x, y, c, cc, gamma, tb_shift; double x, y, c, cc, gamma, tb_shift;
static long time = 0; static long time = 0;
double *delta; double *delta;
static short int first = 1, left_bc[NY];
static int bc_jmin, bc_jmax;
if (first)
{
if ((OSCILLATE_LEFT)&&(BEAM_BC))
{
xy_to_ij(0.0, OSCIL_YMID + OSCIL_YMAX, ij);
bc_jmax = ij[1];
xy_to_ij(0.0, OSCIL_YMID - OSCIL_YMAX, ij);
bc_jmin = ij[1];
}
else
{
bc_jmin = 0;
bc_jmax = NY-1;
}
for (j=0; j<NY-1; j++)
{
if ((OSCILLATE_LEFT)&&(j > bc_jmin)&&(j < bc_jmax)) left_bc[j] = 1;
else left_bc[j] = 0;
}
first = 0;
}
delta = (double *)malloc(NX*NY*sizeof(double)); delta = (double *)malloc(NX*NY*sizeof(double));
@ -370,13 +400,15 @@ void evolve_wave_half_new(double phi_in[NX*NY], double psi_in[NX*NY], double phi
/* left boundary */ /* left boundary */
// 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] = AMPLITUDE*cos((double)time*OMEGA);
if (OSCILLATE_LEFT) // if (OSCILLATE_LEFT)
{ // {
for (j=1; j<NY-1; j++) phi_out[j] = oscillating_bc(time, j); // for (j=1; j<NY-1; j++) phi_out[j] = oscillating_bc(time, j);
printf("Boundary condition %.3lg\n", oscillating_bc(time, 0)); // printf("Boundary condition %.3lg\n", oscillating_bc(time, 0));
} // }
else for (j=1; j<NY-1; j++){ // else
if ((TWOSPEEDS)||(xy_in[j] != 0)){ for (j=1; j<NY-1; j++){
if (left_bc[j]) phi_out[j] = oscillating_bc(time, j);
else if ((TWOSPEEDS)||(xy_in[j] != 0)){
x = phi_in[j]; x = phi_in[j];
y = psi_in[j]; y = psi_in[j];

168
wave_billiard.c
View File

@ -44,12 +44,12 @@
#include <time.h> #include <time.h>
#define MOVIE 0 /* set to 1 to generate movie */ #define MOVIE 0 /* set to 1 to generate movie */
#define DOUBLE_MOVIE 0 /* set to 1 to produce movies for wave height and energy simultaneously */ #define 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 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 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 VARIABLE_IOR 1 /* set to 1 for a variable index of refraction */
#define IOR 191 /* choice of index of refraction, see list in global_pdes.c */ #define IOR 17 /* 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 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 MANDEL_IOR_SCALE -0.05 /* parameter controlling dependence of IoR on Mandelbrot escape speed */
@ -60,8 +60,8 @@
#define NX 3840 /* number of grid points on x axis */ #define NX 3840 /* number of grid points on x axis */
#define NY 2300 /* number of grid points on y axis */ #define NY 2300 /* number of grid points on y axis */
#define XMIN -1.5 #define XMIN -2.0
#define XMAX 2.5 /* x interval */ #define XMAX 2.0 /* x interval */
#define YMIN -1.197916667 #define YMIN -1.197916667
#define YMAX 1.197916667 /* y interval for 9/16 aspect ratio */ #define YMAX 1.197916667 /* y interval for 9/16 aspect ratio */
@ -72,7 +72,7 @@
/* Choice of the billiard table */ /* Choice of the billiard table */
#define B_DOMAIN 947 /* choice of domain shape, see list in global_pdes.c */ #define B_DOMAIN 982 /* choice of domain shape, see list in global_pdes.c */
#define CIRCLE_PATTERN 202 /* pattern of circles or polygons, see list in global_pdes.c */ #define CIRCLE_PATTERN 202 /* pattern of circles or polygons, see list in global_pdes.c */
#define IMAGE_FILE 5 /* for option D_IMAGE */ #define IMAGE_FILE 5 /* for option D_IMAGE */
@ -87,20 +87,20 @@
#define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */ #define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */
#define PDISC_CONNECT_FACTOR 1.5 /* controls which discs are connected for D_CIRCLE_LATTICE_POISSON domain */ #define PDISC_CONNECT_FACTOR 1.5 /* controls which discs are connected for D_CIRCLE_LATTICE_POISSON domain */
#define LAMBDA 0.6 /* parameter controlling the dimensions of domain */ #define LAMBDA 1.4 /* parameter controlling the dimensions of domain */
#define MU 0.075 /* parameter controlling the dimensions of domain */ #define MU 0.05 /* parameter controlling the dimensions of domain */
#define MU_B 1.0 /* parameter controlling the dimensions of domain */ #define MU_B 1.0 /* parameter controlling the dimensions of domain */
#define NPOLY 6 /* number of sides of polygon */ #define NPOLY 3 /* number of sides of polygon */
#define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */ #define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */
#define MDEPTH 6 /* depth of computation of Menger gasket */ #define MDEPTH 6 /* depth of computation of Menger gasket */
#define MRATIO 3 /* ratio defining Menger gasket */ #define MRATIO 3 /* ratio defining Menger gasket */
#define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ #define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */
#define MANDELLIMIT 10.0 /* limit value for approximation of 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 FOCI 1 /* set to 1 to draw focal points of ellipse */
#define NGRIDX 6 /* number of grid point for grid of disks */ #define NGRIDX 6 /* number of grid point for grid of disks */
#define NGRIDY 14 /* number of grid point for grid of disks */ #define NGRIDY 14 /* number of grid point for grid of disks */
#define WALL_WIDTH 0.013 /* width of wall separating lenses */ #define WALL_WIDTH 0.1 /* width of wall separating lenses */
#define WALL_WIDTH_B 0.01 /* width of wall separating lenses */ #define WALL_WIDTH_B 0.05 /* width of wall separating lenses */
#define WALL_WIDTH_RND 0.0 /* proportion of width of width for random arrangements */ #define WALL_WIDTH_RND 0.0 /* proportion of width of width for random arrangements */
#define RADIUS_FACTOR 0.3 /* controls inner radius for C_RING arrangements */ #define RADIUS_FACTOR 0.3 /* controls inner radius for C_RING arrangements */
#define WALL_WIDTH_ASYM 0.75 /* asymmetry of wall width (D_CIRCLE_LATTICE_NONISO) */ #define WALL_WIDTH_ASYM 0.75 /* asymmetry of wall width (D_CIRCLE_LATTICE_NONISO) */
@ -122,38 +122,39 @@
/* Physical parameters of wave equation */ /* Physical parameters of wave equation */
#define TWOSPEEDS 1 /* 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_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 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 42 /* oscillation schedule, see list in global_pdes.c */
#define OSCIL_YMAX 0.35 /* defines oscillation range */ #define OSCIL_YMAX 0.05 /* defines oscilling beam range */
#define OSCIL_YMID -0.9 /* defines oscilling beam midpoint */
#define INITIAL_SHIFT 20.0 /* time shift of initial wave packet (in oscillation periods) */ #define INITIAL_SHIFT 20.0 /* time shift of initial wave packet (in oscillation periods) */
#define WAVE_PACKET_SHIFT 200.0 /* time shift between wave packets (in oscillation periods) */ #define WAVE_PACKET_SHIFT 200.0 /* time shift between wave packets (in oscillation periods) */
#define OMEGA 0.0125 /* frequency of periodic excitation */ #define OMEGA 0.01 /* frequency of periodic excitation */
#define AMPLITUDE 1.0 /* amplitude of periodic excitation */ #define AMPLITUDE 2.0 /* amplitude of periodic excitation */
#define ACHIRP 0.25 /* acceleration coefficient in chirp */ #define ACHIRP 0.25 /* acceleration coefficient in chirp */
#define DAMPING 0.0 /* damping of periodic excitation */ #define DAMPING 0.0 /* damping of periodic excitation */
#define COURANT 0.1 /* Courant number */ #define COURANT 0.25 /* Courant number in medium B */
#define COURANTB 0.035 /* Courant number in medium B */ #define COURANTB 0.1 /* Courant number */
#define GAMMA 0.0 /* damping factor in wave equation */ #define GAMMA 5.0e-6 /* damping factor in wave equation */
#define GAMMAB 0.0 /* damping factor in wave equation */ #define GAMMAB 0.0 /* damping factor in wave equation */
#define GAMMA_SIDES 1.0e-4 /* damping factor on boundary */ #define GAMMA_SIDES 1.0e-4 /* damping factor on boundary */
#define GAMMA_TOPBOT 1.0e-7 /* damping factor on boundary */ #define GAMMA_TOPBOT 1.0e-7 /* damping factor on boundary */
#define KAPPA 0.0 /* "elasticity" term enforcing oscillations */ #define KAPPA 0.0 /* "elasticity" term enforcing oscillations */
#define KAPPA_SIDES 5.0e-4 /* "elasticity" term on absorbing boundary */ #define KAPPA_SIDES 5.0e-4 /* "elasticity" term on absorbing boundary */
#define KAPPA_TOPBOT 0.0 /* "elasticity" term on absorbing boundary */ #define KAPPA_TOPBOT 0.0 /* "elasticity" term on absorbing boundary */
#define OSCIL_LEFT_YSHIFT 40.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 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 */ /* The physical damping coefficient is given by GAMMA/(DT)^2 */
/* Increasing COURANT speeds up the simulation, but decreases accuracy */ /* Increasing COURANT speeds up the simulation, but decreases accuracy */
/* For similar wave forms, COURANT^2*GAMMA should be kept constant */ /* 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 12.5 /* period of oscillating source */ #define OSCILLATING_SOURCE_PERIOD 25.0 /* period of oscillating source */
#define ALTERNATE_OSCILLATING_SOURCE 1 /* set to 1 to alternate sign of oscillating source */ #define ALTERNATE_OSCILLATING_SOURCE 1 /* set to 1 to alternate sign of oscillating source */
#define N_SOURCES 1 /* number of sources, for option draw_sources */ #define N_SOURCES 1 /* number of sources, for option draw_sources */
#define ALTERNATE_SOURCE_PHASES 0 /* set to 1 to alternate initial phases of sources */ #define ALTERNATE_SOURCE_PHASES 0 /* set to 1 to alternate initial phases of sources */
#define MAX_PULSING_TIME 500 /* max time for adding pulses */ #define MAX_PULSING_TIME 750 /* max time for adding pulses */
#define ADD_WAVE_PACKET_SOURCES 0 /* set to 1 to add several sources emitting wave packets */ #define ADD_WAVE_PACKET_SOURCES 0 /* set to 1 to add several sources emitting wave packets */
#define WAVE_PACKET_SOURCE_TYPE 3 /* type of wave packet sources */ #define WAVE_PACKET_SOURCE_TYPE 3 /* type of wave packet sources */
@ -168,10 +169,10 @@
/* Parameters for length and speed of simulation */ /* Parameters for length and speed of simulation */
#define NSTEPS 4800 /* number of frames of movie */ #define NSTEPS 2500 /* number of frames of movie */
#define NVID 8 /* number of iterations between images displayed on screen */ #define NVID 12 /* number of iterations between images displayed on screen */
#define NSEG 1000 /* number of segments of boundary */ #define NSEG 1000 /* number of segments of boundary */
#define INITIAL_TIME 0 /* time after which to start saving frames */ #define INITIAL_TIME 150 /* time after which to start saving frames */
#define BOUNDARY_WIDTH 2 /* width of billiard boundary */ #define BOUNDARY_WIDTH 2 /* width of billiard boundary */
#define PRINT_SPEED 0 /* print speed of moving source */ #define PRINT_SPEED 0 /* print speed of moving source */
#define PRINT_FREQUENCY 0 /* print frequency (for phased array) */ #define PRINT_FREQUENCY 0 /* print frequency (for phased array) */
@ -186,9 +187,9 @@
/* Parameters of initial condition */ /* 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.00025 /* variance of initial condition */
#define INITIAL_WAVELENGTH 0.025 /* wavelength of initial condition */ #define INITIAL_WAVELENGTH 0.015 /* wavelength of initial condition */
/* Plot type, see list in global_pdes.c */ /* Plot type, see list in global_pdes.c */
@ -199,7 +200,7 @@
/* Color schemes */ /* Color schemes */
#define COLOR_PALETTE 11 /* 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 COLOR_PALETTE_B 14 /* Color palette, see list in global_pdes.c */
#define BLACK 1 /* background */ #define BLACK 1 /* background */
@ -213,7 +214,7 @@
#define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */
#define VSHIFT_AMPLITUDE -0.1 /* additional shift for wave amplitude */ #define VSHIFT_AMPLITUDE -0.1 /* additional shift for wave amplitude */
#define VSCALE_AMPLITUDE 1.0 /* additional scaling factor for wave amplitude */ #define VSCALE_AMPLITUDE 1.0 /* additional scaling factor for wave amplitude */
#define E_SCALE 750.0 /* scaling factor for energy representation */ #define E_SCALE 50.0 /* scaling factor for energy representation */
#define LOG_SCALE 0.75 /* scaling factor for energy log representation */ #define LOG_SCALE 0.75 /* scaling factor for energy log representation */
#define LOG_SHIFT 0.75 /* shift of colors on log scale */ #define LOG_SHIFT 0.75 /* shift of colors on log scale */
#define FLUX_SCALE 250.0 /* scaling factor for energy flux represtnation */ #define FLUX_SCALE 250.0 /* scaling factor for energy flux represtnation */
@ -230,14 +231,14 @@
#define HUEMEAN 180.0 /* mean value of hue for color scheme C_HUE */ #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 HUEAMP -180.0 /* amplitude of variation of hue for color scheme C_HUE */
#define DRAW_COLOR_SCHEME 0 /* set to 1 to plot the color scheme */ #define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */
#define COLORBAR_RANGE 1.5 /* scale of color scheme bar */ #define COLORBAR_RANGE 1.7 /* scale of color scheme bar */
#define COLORBAR_RANGE_B 0.12 /* scale of color scheme bar for 2nd part */ #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 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 0 /* set to 1 to draw circular color scheme */
#define CIRC_COLORBAR_B 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 DRAW_WAVE_PROFILE 0 /* set to 1 to draw a profile of the wave */
#define HORIZONTAL_WAVE_PROFILE 0 /* set to 1 to draw wave profile vertically */ #define HORIZONTAL_WAVE_PROFILE 0 /* set to 1 to draw wave profile vertically */
#define VERTICAL_WAVE_PROFILE 1 /* set to 1 to draw wave profile vertically */ #define VERTICAL_WAVE_PROFILE 1 /* set to 1 to draw wave profile vertically */
#define WAVE_PROFILE_X 1.9 /* value of x to sample wave profile */ #define WAVE_PROFILE_X 1.9 /* value of x to sample wave profile */
@ -273,9 +274,11 @@
#define FLOOR 0 /* set to 1 to limit wave amplitude to VMAX */ #define FLOOR 0 /* set to 1 to limit wave amplitude to VMAX */
#define VMAX 10.0 /* max value of wave amplitude */ #define VMAX 10.0 /* max value of wave amplitude */
#define MEAN_FLUX (PLOT == P_TOTAL_ENERGY_FLUX)||(PLOT_B == P_TOTAL_ENERGY_FLUX) #define MEAN_FLUX (PLOT == P_TOTAL_ENERGY_FLUX)||(PLOT_B == P_TOTAL_ENERGY_FLUX)
#define REFRESH_IOR ((IOR == IOR_PERIODIC_WELLS_ROTATING)||(IOR == IOR_PERIODIC_WELLS_ROTATING_LARGE)) #define REFRESH_IOR ((IOR == IOR_PERIODIC_WELLS_ROTATING)||(IOR == IOR_PERIODIC_WELLS_ROTATING_LARGE))
#define XYIN_INITIALISED (B_DOMAIN == D_IMAGE) #define XYIN_INITIALISED (B_DOMAIN == D_IMAGE)
#define BEAM_BC ((OSCILLATION_SCHEDULE == OSC_BEAM_SINE)||(OSCILLATION_SCHEDULE == OSC_BEAM_SINE_DECREASING)||(OSCILLATION_SCHEDULE == OSC_BEAM_SINE_TWOPERIODS)||(OSCILLATION_SCHEDULE == OSC_BEAM_SINE_CHIRP))
double light[2] = {0.40824829, 0.816496581}; /* location of light source for SHADE_2D option*/ double light[2] = {0.40824829, 0.816496581}; /* location of light source for SHADE_2D option*/
@ -293,7 +296,7 @@ FILE *time_series_left, *time_series_right;
// void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX], double *psi_out[NX], // void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX], double *psi_out[NX],
// short int *xy_in[NX]) // short int *xy_in[NX])
void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX], void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX],
short int *xy_in[NX], double *tcc[NX], double *tgamma[NX]) short int *xy_in[NX], double *tcc[NX], double *tc[NX], double *tgamma[NX])
/* time step of field evolution */ /* time step of field evolution */
/* phi is value of field at time t, psi at time t-1 */ /* 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 */ /* this version of the function has been rewritten in order to minimize the number of if-branches */
@ -301,31 +304,55 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
int i, j, iplus, iminus, jplus, jminus, ij[2]; int i, j, iplus, iminus, jplus, jminus, ij[2];
double delta, x, y, c, cc, gamma, tb_shift; double delta, x, y, c, cc, gamma, tb_shift;
static long time = 0; static long time = 0;
static double tc[NX][NY]; // static double tc[NX][NY];
static short int first = 1; static short int first = 1, left_bc[NY];
static int bc_jmin, bc_jmax;
time++; time++;
// if (OSCILLATE_TOPBOT) tb_shift = (int)((X_SHIFT - XMIN)*(double)NX/(XMAX - XMIN)); // if (OSCILLATE_TOPBOT) tb_shift = (int)((X_SHIFT - XMIN)*(double)NX/(XMAX - XMIN));
if (OSCILLATE_TOPBOT) tb_shift = (int)((XMAX - XMIN)*(double)NX/(XMAX - XMIN)); if (OSCILLATE_TOPBOT) tb_shift = (int)((XMAX - XMIN)*(double)NX/(XMAX - XMIN));
/* initialize tables with wave speeds and dissipation */ /* initialize beam range of OSCILLATE_LEFT b.c. */
if (first) if (first)
{ {
for (i=0; i<NX; i++){ if ((OSCILLATE_LEFT)&&(BEAM_BC))
for (j=0; j<NY; j++){ {
if (xy_in[i][j] != 0) xy_to_ij(0.0, OSCIL_YMID + OSCIL_YMAX, ij);
{ bc_jmax = ij[1];
tc[i][j] = sqrt(tcc[i][j]); xy_to_ij(0.0, OSCIL_YMID - OSCIL_YMAX, ij);
} bc_jmin = ij[1];
else if (TWOSPEEDS) }
{ else
tc[i][j] = COURANTB; {
} bc_jmin = 0;
} bc_jmax = NY-1;
}
for (j=0; j<NY-1; j++)
{
if ((OSCILLATE_LEFT)&&(j > bc_jmin)&&(j < bc_jmax)) left_bc[j] = 1;
else left_bc[j] = 0;
} }
first = 0; first = 0;
} }
/* initialize tables with wave speeds and dissipation */
// if (first)
// {
// for (i=0; i<NX; i++){
// for (j=0; j<NY; j++){
// if (xy_in[i][j] != 0)
// {
// tc[i][j] = sqrt(tcc[i][j]);
// }
// else if (TWOSPEEDS)
// {
// tc[i][j] = COURANTB;
// }
// }
// }
// first = 0;
// }
#pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta,x,y) #pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta,x,y)
/* evolution in the bulk */ /* evolution in the bulk */
@ -346,9 +373,11 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
} }
/* left boundary */ /* left boundary */
if (OSCILLATE_LEFT) for (j=1; j<NY-1; j++) phi_out[0][j] = oscillating_bc(time, j); // if (OSCILLATE_LEFT) for (j=bc_jmin+1; j<bc_jmax; j++) phi_out[0][j] = oscillating_bc(time, j);
else for (j=1; j<NY-1; j++){ // else
if ((TWOSPEEDS)||(xy_in[0][j] != 0)){ for (j=1; j<NY-1; j++){
if (left_bc[j]) phi_out[0][j] = oscillating_bc(time, j);
else if ((TWOSPEEDS)||(xy_in[0][j] != 0)){
x = phi_in[0][j]; x = phi_in[0][j];
y = psi_in[0][j]; y = psi_in[0][j];
@ -381,6 +410,7 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
// psi_out[0][j] = x; // psi_out[0][j] = x;
} }
} }
// if (OSCILLATE_LEFT) for (j=bc_jmin+1; j<bc_jmax; j++) phi_out[0][j] = oscillating_bc(time, j);
/* right boundary */ /* right boundary */
for (j=1; j<NY-1; j++){ for (j=1; j<NY-1; j++){
@ -535,10 +565,10 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
} }
/* add oscillating boundary condition on the left corners */ /* add oscillating boundary condition on the left corners */
if (OSCILLATE_LEFT) if ((OSCILLATE_LEFT)&&(!BEAM_BC))
{ {
phi_out[0][0] = oscillating_bc(time, 0); phi_out[0][0] = oscillating_bc(time, 0);
phi_out[0][NY-1] = oscillating_bc(time, NY-1); phi_out[0][NX-1] = oscillating_bc(time, NX-1);
} }
/* for debugging purposes/if there is a risk of blow-up */ /* for debugging purposes/if there is a risk of blow-up */
@ -556,7 +586,7 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
} }
void evolve_wave(double *phi[NX], double *psi[NX], double *tmp[NX], short int *xy_in[NX], double *tcc_table[NX], double *tgamma_table[NX]) void evolve_wave(double *phi[NX], double *psi[NX], double *tmp[NX], short int *xy_in[NX], double *tcc_table[NX], double *tc_table[NX], double *tgamma_table[NX])
/* time step of field evolution */ /* time step of field evolution */
/* phi is value of field at time t, psi at time t-1 */ /* phi is value of field at time t, psi at time t-1 */
{ {
@ -566,11 +596,11 @@ void evolve_wave(double *phi[NX], double *psi[NX], double *tmp[NX], short int *x
// At the beginning w[t] is saved in phi, w[t-1] in psi and tmp is space // At the beginning w[t] is saved in phi, w[t-1] in psi and tmp is space
// for the next wave state w[t+1]. Take w[t] and w[t-1] to calculate the // for the next wave state w[t+1]. Take w[t] and w[t-1] to calculate the
// next wave state. Write this new state in temp // next wave state. Write this new state in temp
evolve_wave_half(phi, psi, tmp, xy_in, tcc_table, tgamma_table); evolve_wave_half(phi, psi, tmp, xy_in, tcc_table, tc_table, tgamma_table);
// now w[t] is saved in tmp, w[t-1] in phi and the result is written to psi // now w[t] is saved in tmp, w[t-1] in phi and the result is written to psi
evolve_wave_half(tmp, phi, psi, xy_in, tcc_table, tgamma_table); evolve_wave_half(tmp, phi, psi, xy_in, tcc_table, tc_table, tgamma_table);
// now w[t] is saved in psi, w[t-1] in tmp and the result is written to phi // now w[t] is saved in psi, w[t-1] in tmp and the result is written to phi
evolve_wave_half(psi, tmp, phi, xy_in, tcc_table, tgamma_table); evolve_wave_half(psi, tmp, phi, xy_in, tcc_table, tc_table, tgamma_table);
// now w[t] is saved in phi, w[t-1] in psi and tmp is free again to take // now w[t] is saved in phi, w[t-1] in psi and tmp is free again to take
// the new wave state w[t+1] in the next call to this function, thus // the new wave state w[t+1] in the next call to this function, thus
// matching the given parameter names again // matching the given parameter names again
@ -600,7 +630,7 @@ void draw_color_bar_palette(int plot, double range, int palette, int circular, i
void animation() void animation()
{ {
double time, scale, ratio, startleft[2], startright[2], sign[N_SOURCES], r2, xy[2], fade_value, yshift, speed = 0.0, a, b, c, x, y, angle = 0.0, x1, ior_angle = 0.0, omega, phase_shift, vshift, dsource, finv, source_amp[N_SOURCES], nx, ny, r, source_periods[N_SOURCES], dperiod; double time, scale, ratio, startleft[2], startright[2], sign[N_SOURCES], r2, xy[2], fade_value, yshift, speed = 0.0, a, b, c, x, y, angle = 0.0, x1, ior_angle = 0.0, omega, phase_shift, vshift, dsource, finv, source_amp[N_SOURCES], nx, ny, r, source_periods[N_SOURCES], dperiod;
double *phi[NX], *psi[NX], *tmp[NX], *total_energy[NX], *average_energy[NX], *color_scale[NX], *total_flux, *tcc_table[NX], *tgamma_table[NX], *fade_table; double *phi[NX], *psi[NX], *tmp[NX], *total_energy[NX], *average_energy[NX], *color_scale[NX], *total_flux, *tcc_table[NX], *tc_table[NX], *tgamma_table[NX], *fade_table;
short int *xy_in[NX]; short int *xy_in[NX];
int i, j, k, s, sample_left[2], sample_right[2], period = 0, fade, source_counter = 0, p, q, first_source = 1, imin, imax, ij[2], source, source_period, source_shift[N_SOURCES], phase; int i, j, k, s, sample_left[2], sample_right[2], period = 0, fade, source_counter = 0, p, q, first_source = 1, imin, imax, ij[2], source, source_period, source_shift[N_SOURCES], phase;
// static int image_counter = 0; // static int image_counter = 0;
@ -626,6 +656,7 @@ void animation()
xy_in[i] = (short int *)malloc(NY*sizeof(short int)); xy_in[i] = (short int *)malloc(NY*sizeof(short int));
color_scale[i] = (double *)malloc(NY*sizeof(double)); color_scale[i] = (double *)malloc(NY*sizeof(double));
tcc_table[i] = (double *)malloc(NX*sizeof(double)); tcc_table[i] = (double *)malloc(NX*sizeof(double));
tc_table[i] = (double *)malloc(NX*sizeof(double));
tgamma_table[i] = (double *)malloc(NX*sizeof(double)); tgamma_table[i] = (double *)malloc(NX*sizeof(double));
} }
@ -737,8 +768,8 @@ void animation()
init_wave_flat(phi, psi, xy_in); init_wave_flat(phi, psi, xy_in);
init_ior_2d(xy_in, tcc_table, tgamma_table, ior_angle); init_ior_2d(xy_in, tcc_table, tc_table, tgamma_table, ior_angle);
if (FADE_IN_OBSTACLE) init_fade_table(tcc_table, fade_table); if (FADE_IN_OBSTACLE) init_fade_table(tcc_table, xy_in, fade_table);
// init_circular_wave(-LAMBDA, 0.0, phi, psi, xy_in); // init_circular_wave(-LAMBDA, 0.0, phi, psi, xy_in);
// x = XMIN + (XMAX - XMIN)*rand()/RAND_MAX; // x = XMIN + (XMAX - XMIN)*rand()/RAND_MAX;
@ -830,7 +861,7 @@ void animation()
else draw_wave_epalette(phi, psi, total_energy, average_energy, total_flux, color_scale, xy_in, scale, i, PLOT, COLOR_PALETTE, 0, 1.0); else draw_wave_epalette(phi, psi, total_energy, average_energy, total_flux, color_scale, xy_in, scale, i, PLOT, COLOR_PALETTE, 0, 1.0);
for (j=0; j<NVID; j++) for (j=0; j<NVID; j++)
{ {
evolve_wave(phi, psi, tmp, xy_in, tcc_table, tgamma_table); evolve_wave(phi, psi, tmp, xy_in, tcc_table, tc_table, tgamma_table);
if (SAVE_TIME_SERIES) if (SAVE_TIME_SERIES)
{ {
wave_value = (long int)(phi[sample_left[0]][sample_left[1]]*1.0e16); wave_value = (long int)(phi[sample_left[0]][sample_left[1]]*1.0e16);
@ -848,7 +879,7 @@ void animation()
/* add oscillating waves */ /* add oscillating waves */
wave_source_x[0] = -1.0; wave_source_x[0] = -1.0;
wave_source_y[0] = 0.0; wave_source_y[0] = -0.5;
source_periods[0] = OSCILLATING_SOURCE_PERIOD; source_periods[0] = OSCILLATING_SOURCE_PERIOD;
source_amp[0] = INITIAL_AMP; source_amp[0] = INITIAL_AMP;
for (source = 0; source < N_SOURCES; source++) for (source = 0; source < N_SOURCES; source++)
@ -871,7 +902,7 @@ void animation()
{ {
ior_angle = ior_angle_schedule(i); ior_angle = ior_angle_schedule(i);
printf("IOR angle = %.5lg\n", ior_angle); printf("IOR angle = %.5lg\n", ior_angle);
init_ior_2d(xy_in, tcc_table, tgamma_table, ior_angle); init_ior_2d(xy_in, tcc_table, tc_table, tgamma_table, ior_angle);
printf("speed = %.5lg\n", tcc_table[3*NX/4][NY/2]); printf("speed = %.5lg\n", tcc_table[3*NX/4][NY/2]);
} }
if (B_DOMAIN == D_MICHELSON_MOVING) if (B_DOMAIN == D_MICHELSON_MOVING)
@ -985,6 +1016,7 @@ void animation()
free(xy_in[i]); free(xy_in[i]);
free(color_scale[i]); free(color_scale[i]);
free(tcc_table[i]); free(tcc_table[i]);
free(tc_table[i]);
free(tgamma_table[i]); free(tgamma_table[i]);
} }

7
wave_common.c
View File

@ -1734,14 +1734,17 @@ void draw_wave_timeseries_old(double *values, int size, int fade, double fade_va
} }
void init_fade_table(double *tcc_table[NX], double *fade_table) void init_fade_table(double *tcc_table[NX], short int *xy_in[NX], double *fade_table)
{ {
int i, j; int i, j;
for (i=0; i<NX; i++) for (i=0; i<NX; i++)
for (j=0; j<NY; j++) for (j=0; j<NY; j++)
fade_table[i*NY + j] = pow(tcc_table[i][j]/courant2, 0.25); {
if (xy_in[i][j] == 2) fade_table[i*NY + j] = 0.0;
else fade_table[i*NY + j] = pow(tcc_table[i][j]/courant2, 0.1);
// fade_table[i*NY + j] = tcc_table[i][j]/courant2; // fade_table[i*NY + j] = tcc_table[i][j]/courant2;
}
} }
void draw_wave_highres_palette(int size, double *phi[NX], double *psi[NX], double *total_energy[NX], double *average_energy[NX], double *fade_table, double *total_flux, double *color_scale[NX], short int *xy_in[NX], double scale, int time, int plot, int palette, int pnumber, int fade, double fade_value) void draw_wave_highres_palette(int size, double *phi[NX], double *psi[NX], double *total_energy[NX], double *average_energy[NX], double *fade_table, double *total_flux, double *color_scale[NX], short int *xy_in[NX], double scale, int time, int plot, int palette, int pnumber, int fade, double fade_value)

91
wave_comparison.c
View File

@ -44,12 +44,13 @@
#include <time.h> #include <time.h>
#define MOVIE 0 /* set to 1 to generate movie */ #define MOVIE 0 /* set to 1 to generate movie */
#define DOUBLE_MOVIE 0 /* set to 1 to produce movies for wave height and energy simultaneously */ #define 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 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 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 VARIABLE_IOR 1 /* set to 1 for a variable index of refraction */
#define IOR 10 /* choice of index of refraction, see list in global_pdes.c */ #define IOR 183 /* choice of index of refraction, see list in global_pdes.c */
#define IOR_B 183 /* 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 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 MANDEL_IOR_SCALE -0.05 /* parameter controlling dependence of IoR on Mandelbrot escape speed */
@ -61,10 +62,16 @@
#define NX 3840 /* number of grid points on x axis */ #define NX 3840 /* number of grid points on x axis */
#define NY 2300 /* number of grid points on y axis */ #define NY 2300 /* number of grid points on y axis */
#define YMID 1150 /* mid point of display */ #define YMID 1150 /* mid point of display */
#define XMIN -2.0
#define XMAX 2.0 /* x interval */ #define XMIN -2.2
#define YMIN -1.197916667 #define XMAX 3.8 /* x interval */
#define YMAX 1.197916667 /* y interval for 9/16 aspect ratio */ #define YMIN -1.796875
#define YMAX 1.796875 /* y interval for 9/16 aspect ratio */
// #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 HIGHRES 1 /* set to 1 if resolution of grid is double that of displayed image */ #define HIGHRES 1 /* set to 1 if resolution of grid is double that of displayed image */
@ -72,8 +79,8 @@
/* Choice of the billiard table */ /* Choice of the billiard table */
#define B_DOMAIN 94 /* choice of domain shape, see list in global_pdes.c */ #define B_DOMAIN 196 /* choice of domain shape, see list in global_pdes.c */
#define B_DOMAIN_B 94 /* choice of domain shape, see list in global_pdes.c */ #define B_DOMAIN_B 196 /* choice of domain shape, see list in global_pdes.c */
#define CIRCLE_PATTERN 13 /* pattern of circles, see list in global_pdes.c */ #define CIRCLE_PATTERN 13 /* pattern of circles, see list in global_pdes.c */
#define CIRCLE_PATTERN_B 13 /* pattern of circles, see list in global_pdes.c */ #define CIRCLE_PATTERN_B 13 /* pattern of circles, see list in global_pdes.c */
@ -92,14 +99,15 @@
#define HEX_NONUNIF_COMPRESSSION 0.15 /* compression factor for HEX_NONUNIF pattern */ #define HEX_NONUNIF_COMPRESSSION 0.15 /* compression factor for HEX_NONUNIF pattern */
#define HEX_NONUNIF_COMPRESSSION_B -0.15 /* compression factor for HEX_NONUNIF pattern */ #define HEX_NONUNIF_COMPRESSSION_B -0.15 /* compression factor for HEX_NONUNIF pattern */
#define LAMBDA 1.5 /* parameter controlling the dimensions of domain */ #define LAMBDA 1.0 /* parameter controlling the dimensions of domain */
#define MU 0.06 /* parameter controlling the dimensions of domain */ #define MU 0.06 /* parameter controlling the dimensions of domain */
#define MU_B 0.06 /* parameter controlling the dimensions of domain */ #define MU_B 0.06 /* parameter controlling the dimensions of domain */
#define MUB 0.06 /* parameter controlling the dimensions of domain */ #define MUB 0.06 /* parameter controlling the dimensions of domain */
#define NPOLY 3 /* number of sides of polygon */ #define NPOLY 3 /* number of sides of polygon */
#define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */ #define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */
#define APOLY_B 2.0 /* angle by which to turn polygon, in units of Pi/2 */ #define APOLY_B 2.0 /* angle by which to turn polygon, in units of Pi/2 */
#define MDEPTH 4 /* depth of computation of Menger gasket */ #define MDEPTH 5 /* depth of computation of Menger gasket */
#define MDEPTH_B 10 /* depth of computation of Menger gasket */
#define MRATIO 3 /* ratio defining Menger gasket */ #define MRATIO 3 /* ratio defining Menger gasket */
#define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ #define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */
#define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */ #define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */
@ -129,15 +137,16 @@
/* Physical parameters of wave equation */ /* 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_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 OSCILLATE_TOPBOT 1 /* set to 1 to enforce a planar wave on top and bottom boundary */
#define OSCIL_YMID -0.9 /* defines oscilling beam midpoint */
#define OMEGA 0.024 /* frequency of periodic excitation */ #define OMEGA 0.002 /* frequency of periodic excitation */
#define AMPLITUDE 1.0 /* amplitude of periodic excitation */ #define AMPLITUDE 1.0 /* amplitude of periodic excitation */
#define DAMPING 0.0 /* damping of periodic excitation */ #define DAMPING 0.0 /* damping of periodic excitation */
#define COURANT 0.1 /* Courant number */ #define COURANT 0.1 /* Courant number */
#define COURANTB 0.063 /* Courant number in medium B */ #define COURANTB 0.1 /* Courant number in medium B */
#define GAMMA 0.0 /* damping factor in wave equation */ #define GAMMA 0.0 /* damping factor in wave equation */
#define GAMMAB 0.0 /* damping factor in wave equation */ #define GAMMAB 0.0 /* damping factor in wave equation */
#define GAMMA_SIDES 1.0e-4 /* damping factor on boundary */ #define GAMMA_SIDES 1.0e-4 /* damping factor on boundary */
@ -150,7 +159,7 @@
/* Increasing COURANT speeds up the simulation, but decreases accuracy */ /* Increasing COURANT speeds up the simulation, but decreases accuracy */
/* For similar wave forms, COURANT^2*GAMMA should be kept constant */ /* 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 15.625 /* period of oscillating source */ #define OSCILLATING_SOURCE_PERIOD 15.625 /* period of oscillating source */
#define ALTERNATE_OSCILLATING_SOURCE 1 /* set to 1 to alternate sign of oscillating source */ #define ALTERNATE_OSCILLATING_SOURCE 1 /* set to 1 to alternate sign of oscillating source */
#define N_SOURCES 2 /* number of sources, for option draw_sources */ #define N_SOURCES 2 /* number of sources, for option draw_sources */
@ -160,11 +169,11 @@
/* Boundary conditions, see list in global_pdes.c */ /* Boundary conditions, see list in global_pdes.c */
#define B_COND 1 #define B_COND 4
/* Parameters for length and speed of simulation */ /* Parameters for length and speed of simulation */
#define NSTEPS 3900 /* number of frames of movie */ #define NSTEPS 2300 /* number of frames of movie */
#define NVID 7 /* number of iterations between images displayed on screen */ #define NVID 7 /* number of iterations between images displayed on screen */
#define NSEG 100 /* number of segments of boundary */ #define NSEG 100 /* number of segments of boundary */
#define INITIAL_TIME 0 /* time after which to start saving frames */ #define INITIAL_TIME 0 /* time after which to start saving frames */
@ -189,12 +198,12 @@
#define PLOT 0 #define PLOT 0
#define PLOT_B 5 #define PLOT_B 3
/* Color schemes */ /* Color schemes */
#define COLOR_PALETTE 14 /* Color palette, see list in global_pdes.c */ #define COLOR_PALETTE 11 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 14 /* 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 */ #define BLACK 1 /* background */
#define BLACK_TEXT 1 /* set to 1 to write text in black instead of white */ #define BLACK_TEXT 1 /* set to 1 to write text in black instead of white */
@ -207,9 +216,9 @@
#define PHASE_FACTOR 1.0 /* factor in computation of phase in color scheme P_3D_PHASE */ #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 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 ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */
#define VSHIFT_AMPLITUDE -0.2 /* additional shift for wave amplitude */ #define VSHIFT_AMPLITUDE -0.1 /* additional shift for wave amplitude */
#define VSCALE_AMPLITUDE 0.4 /* additional scaling factor for wave amplitude */ #define VSCALE_AMPLITUDE 1.0 /* additional scaling factor for wave amplitude */
#define E_SCALE 10.0 /* scaling factor for energy representation */ #define E_SCALE 300.0 /* scaling factor for energy representation */
#define LOG_SCALE 1.0 /* scaling factor for energy log representation */ #define LOG_SCALE 1.0 /* scaling factor for energy log representation */
#define LOG_SHIFT 2.0 /* shift of colors on log scale */ #define LOG_SHIFT 2.0 /* shift of colors on log scale */
#define FLUX_SCALE 5.0e3 /* scaling factor for enegy flux represtnation */ #define FLUX_SCALE 5.0e3 /* scaling factor for enegy flux represtnation */
@ -222,9 +231,9 @@
#define HUEMEAN 220.0 /* mean value of hue for color scheme C_HUE */ #define HUEMEAN 220.0 /* mean value of hue for color scheme C_HUE */
#define HUEAMP -220.0 /* amplitude of variation of hue for color scheme C_HUE */ #define HUEAMP -220.0 /* amplitude of variation of hue for color scheme C_HUE */
#define DRAW_COLOR_SCHEME 0 /* set to 1 to plot the color scheme */ #define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */
#define COLORBAR_RANGE 1.5 /* scale of color scheme bar */ #define COLORBAR_RANGE 1.2 /* scale of color scheme bar */
#define COLORBAR_RANGE_B 12.5 /* scale of color scheme bar for 2nd part */ #define COLORBAR_RANGE_B 5.0 /* scale of color scheme bar for 2nd part */
#define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ #define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */
@ -304,7 +313,7 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
/* phi is value of field at time t, psi at time t-1 */ /* phi is value of field at time t, psi at time t-1 */
{ {
int i, j, iplus, iminus, jplus, jminus, jmid = NY/2; int i, j, iplus, iminus, jplus, jminus, jmid = NY/2;
double delta, x, y, c, cc, gamma; double delta, x, y, c, cc, gamma, tb_shift;
static long time = 0; static long time = 0;
static double tc[NX][NY], tgamma[NX][NY]; static double tc[NX][NY], tgamma[NX][NY];
static short int first = 1; static short int first = 1;
@ -333,6 +342,8 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
} }
first = 0; first = 0;
} }
if (OSCILLATE_TOPBOT) tb_shift = (int)((XMAX - XMIN)*(double)NX/(XMAX - XMIN));
#pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta,x,y,c,cc,gamma) #pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta,x,y,c,cc,gamma)
/* evolution in the bulk */ /* evolution in the bulk */
@ -512,7 +523,15 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
x = phi_in[i][0]; x = phi_in[i][0];
y = psi_in[i][0]; y = psi_in[i][0];
switch (B_COND) { if ((OSCILLATE_TOPBOT)&&(i < tb_shift)&&(i<NX-1)&&(i>0))
{
iplus = i+1;
iminus = i-1; if (iminus < 0) iminus = 0;
delta = phi_in[iplus][0] + phi_in[iminus][0] + - 2.0*x;
phi_out[i][0] = -y + 2*x + tcc[i][0]*delta - KAPPA*x - tgamma[i][0]*(x-y);
}
else switch (B_COND) {
case (BC_DIRICHLET): case (BC_DIRICHLET):
{ {
iplus = i+1; if (iplus == NX) iplus = NX-1; iplus = i+1; if (iplus == NX) iplus = NX-1;
@ -568,8 +587,16 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
if ((TWOSPEEDS)||(xy_in[i][NY-1] != 0)){ if ((TWOSPEEDS)||(xy_in[i][NY-1] != 0)){
x = phi_in[i][NY-1]; x = phi_in[i][NY-1];
y = psi_in[i][NY-1]; y = psi_in[i][NY-1];
switch (B_COND) { if ((OSCILLATE_TOPBOT)&&(i < tb_shift)&&(i<NX-1)&&(i>0))
{
iplus = i+1;
iminus = i-1; if (iminus < 0) iminus = 0;
delta = phi_in[iplus][NY-1] + phi_in[iminus][NY-1] + - 2.0*x;
phi_out[i][NY-1] = -y + 2*x + tcc[i][NY-1]*delta - KAPPA*x - tgamma[i][NY-1]*(x-y);
}
else switch (B_COND) {
case (BC_DIRICHLET): case (BC_DIRICHLET):
{ {
iplus = i+1; if (iplus == NX) iplus = NX-1; iplus = i+1; if (iplus == NX) iplus = NX-1;

3
wave_energy.c
View File

@ -90,6 +90,7 @@
#define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */ #define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */
#define APOLY_B 0.335 /* angle by which to turn polygon, in units of Pi/2 */ #define APOLY_B 0.335 /* angle by which to turn polygon, in units of Pi/2 */
#define MDEPTH 4 /* depth of computation of Menger gasket */ #define MDEPTH 4 /* depth of computation of Menger gasket */
#define MDEPTH_B 10 /* depth of computation of Menger gasket */
#define MRATIO 3 /* ratio defining Menger gasket */ #define MRATIO 3 /* ratio defining Menger gasket */
#define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ #define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */
#define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */ #define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */
@ -238,11 +239,13 @@
#define MAZE_WIDTH 0.02 /* half width of maze walls */ #define MAZE_WIDTH 0.02 /* half width of maze walls */
#define VARIABLE_IOR 0 /* set to 1 for a variable index of refraction */ #define VARIABLE_IOR 0 /* set to 1 for a variable index of refraction */
#define IOR 7 /* choice of index of refraction, see list in global_pdes.c */ #define IOR 7 /* choice of index of refraction, see list in global_pdes.c */
#define IOR_B 183 /* 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 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 MANDEL_IOR_SCALE -0.05 /* parameter controlling dependence of IoR on Mandelbrot escape speed */
#define WAVE_PACKET_SOURCE_TYPE 1 /* type of wave packet sources */ #define WAVE_PACKET_SOURCE_TYPE 1 /* type of wave packet sources */
#define N_WAVE_PACKETS 15 /* number of wave packets */ #define N_WAVE_PACKETS 15 /* number of wave packets */
#define OSCIL_LEFT_YSHIFT 40.0 /* y-dependence of left oscillation (for non-horizontal waves) */ #define OSCIL_LEFT_YSHIFT 40.0 /* y-dependence of left oscillation (for non-horizontal waves) */
#define OSCIL_YMID -0.9 /* defines oscilling beam midpoint */
#define DRAW_WAVE_PROFILE 0 /* set to 1 to draw a profile of the wave */ #define DRAW_WAVE_PROFILE 0 /* set to 1 to draw a profile of the wave */
#define MU_B 1.0 /* parameter controlling the dimensions of domain */ #define MU_B 1.0 /* parameter controlling the dimensions of domain */
#define DRAW_WAVE_PROFILE 0 /* set to 1 to draw a profile of the wave */ #define DRAW_WAVE_PROFILE 0 /* set to 1 to draw a profile of the wave */

3
wave_sphere.c
View File

@ -39,7 +39,7 @@
#include <omp.h> #include <omp.h>
#include <time.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 0 /* set to 1 to produce movies for wave height and energy simultaneously */ #define DOUBLE_MOVIE 0 /* 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 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 NO_EXTRA_BUFFER_SWAP 1 /* some OS require one less buffer swap when recording images */
@ -144,6 +144,7 @@
#define KAPPA_SIDES 5.0e-4 /* "elasticity" term on absorbing boundary */ #define KAPPA_SIDES 5.0e-4 /* "elasticity" term on absorbing boundary */
#define KAPPA_TOPBOT 0.0 /* "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) */ #define OSCIL_LEFT_YSHIFT 0.0 /* y-dependence of left oscillation (for non-horizontal waves) */
#define OSCIL_YMID -0.9 /* defines oscilling beam midpoint */
/* The Courant number is given by c*DT/DX, where DT is the time step and DX the lattice spacing */ /* 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 */ /* The physical damping coefficient is given by GAMMA/(DT)^2 */
/* Increasing COURANT speeds up the simulation, but decreases accuracy */ /* Increasing COURANT speeds up the simulation, but decreases accuracy */