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

Add files via upload

Browse Source 
Bugs corrected in lennardjones.c
New domains for wave_billiard.c
This commit is contained in:
Nils Berglund 2025-12-30 18:32:47 +01:00 committed by GitHub
parent aa023d3799
commit b0dbe991d1
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
8 changed files with 2068 additions and 301 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

37
global_ljones.c
View File

@ -26,6 +26,7 @@
#define NMAXSHOVELS 50 /* max number of shovels */
#define NMAX_TRIANGLES_PER_OBSTACLE 10 /* max number of triangles per obstacle */
#define NMAX_TRIANGLES_PER_FACET 4 /* max number of triangles per facet between obstacles */
#define NMAX_ABSORBERS 10 /* max number of absorbing discs */
#define C_SQUARE 0 /* square grid of circles */
#define C_HEX 1 /* hexagonal/triangular grid of circles */
@ -36,6 +37,7 @@
#define C_CLOAK_A 6 /* first optimized invisibility cloak */
#define C_LASER 7 /* laser fight in a room of mirrors */
#define C_POISSON_DISC 8 /* Poisson disc sampling */
#define C_POISSON_DISC_SPHERE 81 /* Poisson disc sampling in spherical geometry */
#define C_GOLDEN_MEAN 10 /* pattern based on vertical shifts by golden mean */
#define C_GOLDEN_SPIRAL 11 /* spiral pattern based on golden mean */
@ -43,6 +45,8 @@
#define C_POOL_TABLE 20 /* pool table initial position */
#define C_SQUARE_SPHERE 30 /* "square" grid of circles on a sphere */
#define C_ONE 97 /* one single circle, as for Sinai */
#define C_TWO 98 /* two concentric circles of different type */
#define C_NOTHING 99 /* no circle at all, for comparisons */
@ -83,6 +87,7 @@
/* pattern of additional repelling segments */
#define S_RECTANGLE 0 /* segments forming a rectangle */
#define S_CYLINDRICAL 100 /* two lines at top and bottom */
#define S_BOTTOM 110 /* one line at bottom */
#define S_CUP 1 /* segments forming a cup (for increasing gravity) */
#define S_HOURGLASS 2 /* segments forming an hour glass */
#define S_PENTA 3 /* segments forming a pentagon with 3 angles of 120° and 2 right angles */
@ -252,6 +257,7 @@
#define CHEM_H2O_H_OH 20 /* H2O <-> H+ + OH- */
#define CHEM_2H2O_H3O_OH 21 /* 2 H2O <-> H3O+ + OH- */
#define CHEM_AGGREGATION 22 /* agregation of molecules coming close */
#define CHEM_AGGREGATION_MAX 221 /* agregation of molecules coming close, with max total partner number */
#define CHEM_AGGREGATION_CHARGE 23 /* agregation of charged molecules coming close */
#define CHEM_AGGREGATION_CHARGE_NOTRIANGLE 231 /* agregation of charged molecules coming close, no multiple pairings */
#define CHEM_AGGREGATION_NNEIGH 24 /* agregation of molecules with limitation on neighbours */
@ -409,6 +415,14 @@
#define VP_ORBIT2 11 /* rotate in a plane specified by max latitude */
#define VP_POLAR 2 /* polar orbit */
/* Absorber pattern */
#define AP_ONE 0 /* one single absorber */
#define AP_TWO 1 /* two absorbers, symmetric wrp x-axis */
#define AP_CUBE 2 /* absorbers on the vertices of a cube */
#define AP_DODECA 3 /* absorbers on the vertices of a dodecahedron */
#define AP_ICOSA 4 /* absorbers on the vertices of an icosahedron */
/* Color schemes */
#define C_LUM 0 /* color scheme modifies luminosity (with slow drift of hue) */
@ -571,6 +585,12 @@ typedef struct
short int chessboard; /* has value 1 on a chessboard, for some arrangements */
} t_obstacle;
typedef struct
{
double xc, yc; /* center of absorber */
double radius; /* radius of absorber */
} t_absorber;
typedef struct
{
int i[3]; /* indices of obstacles forming triangle */
@ -678,6 +698,7 @@ typedef struct
typedef struct
{
int nactive; /* number of active particles */
int nabsorbed; /* number of absorbed particles */
double beta; /* inverse temperature */
double mean_energy; /* mean energy */
double krepel; /* force constant */
@ -704,24 +725,12 @@ typedef struct
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 */
short int locked; /* has value 1 if color of pixel is fixed */
} 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, nabsorbers = 0;
FILE *lj_log;

12
global_pdes.c
View File

@ -41,6 +41,11 @@
#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_TWO_PARABOLAS_SEMITRANS 193 /* 192 with a semitransparent part */
#define D_TWO_PARABOLAS_LENS 194 /* 193 with a lens */
#define D_TWO_PARABOLAS_CLOSED_SEMITRANS 195 /* cavity with two parabolic sides and semitransparent part */
#define D_TWO_PARABOLAS_CLOSED_LENS 1951 /* 105 with a lens as window */
#define D_TWO_PARABOLAS_CLOSED_LENS_R 1952 /* 105 with a lens of adaptable radius as window */
#define D_CIRCLES 20 /* several circles */
#define D_CIRCLES_IN_RECT 201 /* several circles in a rectangle */
@ -57,6 +62,7 @@
#define D_POLYGONS 40 /* several polygons */
#define D_VONKOCH 41 /* von Koch snowflake fractal */
#define D_STAR 42 /* star shape */
#define D_STAR_CHANNEL 421 /* star shape with an attached channel */
#define D_FRESNEL 43 /* Fresnel lens */
#define D_NOISEPANEL 44 /* zigzag noise insulating panel */
#define D_NOISEPANEL_RECT 441 /* comparison between zigzag noise insulating panel and flat walls */
@ -102,6 +108,7 @@
#define D_TESLA_FOUR 72 /* four Tesla valves */
#define D_TREE 73 /* Christmas tree, to use with IOR_TREE */
#define D_TREE_BEAM 731 /* Christmas tree with an attached wave guide */
#define D_MICHELSON 74 /* Michelson interferometer, to use with IOR_MICHELSON */
#define D_MICHELSON_MOVING 741 /* moving Michelson interferometer, to use with IOR_MICHELSON */
#define D_RITCHEY_CHRETIEN_SPHERICAL 75 /* Ritchey-Chrétien telescope with spherical mirrors */
@ -136,7 +143,9 @@
#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 */
#define D_DISC_WAVEGUIDE_ELLIPSE 982 /* an ellipse with a shifted attached waveguide */
#define D_WAVEGUIDE_ELLIPSE_OBLIQUE 983 /* an ellipse with an angled attached waveguide */
#define D_CIRCLE_PAIRS 99 /* pairs of circles connected by channels */
/* for wave_sphere.c */
@ -163,6 +172,7 @@
#define IM_HAPPY_NEW_YEAR 3 /* "Happy New Year 2025" */
#define IM_DICKSON 4 /* Dickson fjord and others nearbys in Greenland */
#define IM_DICKSON_ZOOM 5 /* zoom on Dickson fjord in Greenland */
#define IM_HAPPY_NEW_YEAR_TWOSIX 6 /* "Happy New Year 2026" */
#define NMAXCIRCLES 10000 /* total number of circles/polygons (must be at least NCX*NCY for square grid) */
#define NMAXPOLY 50000 /* maximal number of vertices of polygonal lines (for von Koch et al) */

303
lennardjones.c
View File

@ -58,34 +58,34 @@
#define YMIN 0.0
#define YMAX 3.141592654 /* y interval for 9/16 aspect ratio */
#define INITXMIN 0.025
#define INITXMAX 6.255 /* x interval for initial condition */
#define INITYMIN 0.4
#define INITYMAX 2.74 /* y interval for initial condition */
#define INITXMIN 2.95
#define INITXMAX 3.4 /* x interval for initial condition */
#define INITYMIN 1.37
#define INITYMAX 1.77 /* y interval for initial condition */
#define THERMOXMIN -1.25
#define THERMOXMAX 1.25 /* x interval for initial condition */
#define THERMOYMIN 0.0
#define THERMOYMAX 0.75 /* y interval for initial condition */
#define ADDXMIN -1.9
#define ADDXMAX 1.9 /* x interval for adding particles */
#define ADDYMIN 1.2
#define ADDYMAX 1.3 /* y interval for adding particles */
#define ADDXMIN 0.0
#define ADDXMAX 0.1 /* x interval for adding particles */
#define ADDYMIN 1.57
#define ADDYMAX 1.57 /* y interval for adding particles */
#define ADDRMIN 2.0
#define ADDRMAX 2.1 /* r interval for adding particles */
#define BCXMIN 0.0
#define BCXMAX 6.283185307 /* x interval for boundary condition */
#define BCYMIN 0.05
#define BCYMAX 3.091592654 /* y interval for boundary condition */
#define BCYMIN 0.3
#define BCYMAX 2.841592654 /* y interval for boundary condition */
#define OBSXMIN -2.0
#define OBSXMAX 2.0 /* x interval for motion of obstacle */
#define OBSYMIN -1.125
#define OBSYMAX 1.125 /* x interval for motion of obstacle */
#define CIRCLE_PATTERN 0 /* pattern of circles, see list in global_ljones.c */
#define CIRCLE_PATTERN 20 /* pattern of circles, see list in global_ljones.c */
#define ADD_INITIAL_PARTICLES 0 /* set to 1 to add a second type of particles */
#define CIRCLE_PATTERN_B 0 /* pattern of circles for additional particles */
@ -119,9 +119,9 @@
#define OBSTACLE_OMEGA 300.0 /* obstacle rotation speed */
#define TWO_TYPES 0 /* set to 1 to have two types of particles */
#define TYPE_PROPORTION 0.75 /* proportion of particles of first type */
#define TYPE_PROPORTION 0.5 /* proportion of particles of first type */
#define TWOTYPE_CONFIG 0 /* choice of types, see TTC_ list in global_ljones.c */
#define SYMMETRIZE_FORCE 0 /* set to 1 to symmetrize two-particle interaction, only needed if particles are not all the same */
#define SYMMETRIZE_FORCE 1 /* set to 1 to symmetrize two-particle interaction, only needed if particles are not all the same */
#define CENTER_PX 0 /* set to 1 to center horizontal momentum */
#define CENTER_PY 0 /* set to 1 to center vertical momentum */
#define CENTER_PANGLE 0 /* set to 1 to center angular momentum */
@ -139,10 +139,10 @@
#define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */
#define LAMBDA 0.2 /* parameter controlling the dimensions of domain */
#define MU 0.022 /* parameter controlling radius of particles */
#define MU_B 0.022 /* parameter controlling radius of particles of second type */
#define MU_ADD 0.022 /* parameter controlling radius of added particles */
#define MU_ADD_B 0.022 /* parameter controlling radius of added particles */
#define MU 0.02 /* parameter controlling radius of particles */
#define MU_B 0.02 /* parameter controlling radius of particles of second type */
#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 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 */
@ -151,8 +151,8 @@
#define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */
#define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */
#define FOCI 1 /* set to 1 to draw focal points of ellipse */
#define NGRIDX 40 /* number of grid point for grid of disks */
#define NGRIDY 20 /* number of grid point for grid of disks */
#define NGRIDX 15 /* number of grid point for grid of disks */
#define NGRIDY 15 /* number of grid point for grid of disks */
#define EHRENFEST_RADIUS 0.9 /* radius of container for Ehrenfest urn configuration */
#define EHRENFEST_WIDTH 0.035 /* width of tube for Ehrenfest urn configuration */
#define TWO_CIRCLES_RADIUS_RATIO 0.8 /* ratio of radii for S_TWO_CIRCLES_EXT segment configuration */
@ -160,6 +160,9 @@
#define NOBSX 24
#define NOBSY 14 /* obstacles for O_HEX obstacle pattern */
#define NTREES 15 /* number of trees in S_TREES */
#define OFSSET_TREES 0.5 /* vertical offset in S_TREES_B */
#define SLOPE_TREES 0.015 /* slope in S_TREES_B (default: 0.3) */
#define SLOPE_TREES_B 0.015 /* slope in S_TREES_B (default: 0.25) */
#define X_SHOOTER -0.2
#define Y_SHOOTER -0.6
@ -168,8 +171,8 @@
/* Parameters for length and speed of simulation */
#define NSTEPS 2300 /* number of frames of movie */
#define NVID 100 /* number of iterations between images displayed on screen */
#define NSTEPS 3800 /* number of frames of movie */
#define NVID 100 /* number of iterations between images displayed on screen */
#define NSEG 25 /* number of segments of boundary of circles */
#define INITIAL_TIME 0 /* time after which to start saving frames */
#define OBSTACLE_INITIAL_TIME 0 /* time after which to start moving obstacle */
@ -190,13 +193,13 @@
/* Plot type, see list in global_ljones.c */
#define PLOT 5
#define PLOT 13
#define PLOT_B 13 /* plot type for second movie */
/* Background color depending on particle properties */
#define COLOR_BACKGROUND 0 /* set to 1 to color background */
#define BG_COLOR 10 /* type of background coloring, see list in global_ljones.c */
#define BG_COLOR 1 /* 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 */
@ -250,8 +253,8 @@
#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 PRINT_PARAMETERS 0 /* set to 1 to print certain parameters */
#define PRINT_TEMPERATURE 1 /* set to 1 to print current temperature */
#define PRINT_PARAMETERS 1 /* set to 1 to print certain parameters */
#define PRINT_TEMPERATURE 0 /* set to 1 to print current temperature */
#define PRINT_ANGLE 0 /* set to 1 to print obstacle orientation */
#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 */
@ -259,6 +262,7 @@
#define PRINT_SEGMENTS_FORCE 0 /* set to 1 to print force on segments */
#define PRINT_NPARTICLES 0 /* print number of active particles */
#define PRINT_TYPE_PROP 0 /* print type proportion */
#define PRINT_NABSORBED 1 /* print number of absorbed particles */
#define FORCE_FACTOR 0.1 /* factor controlling length of force vector */
/* particle properties */
@ -268,7 +272,7 @@
#define PARTICLE_HUE_MIN 359.0 /* color of original particle */
#define PARTICLE_HUE_MAX 0.0 /* color of saturated particle */
#define PARTICLE_EMIN 0.0 /* energy of particle with coolest color */
#define PARTICLE_EMAX 2000.0 /* energy of particle with hottest color */
#define PARTICLE_EMAX 250000.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 */
@ -277,15 +281,17 @@
#define OBSTACLE_VMAX 4.0 /* speed of obstacle with largest luminosity */
#define HUE_TYPE0 320.0 /* hue of particles of type 0 */
#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 100.0 /* hue of particles of type 2 */
#define HUE_TYPE3 140.0 /* hue of particles of type 3 */
#define HUE_TYPE4 200.0 /* hue of particles of type 4 */
#define HUE_TYPE5 60.0 /* hue of particles of type 5 */
#define HUE_TYPE6 130.0 /* hue of particles of type 6 */
#define HUE_TYPE7 200.0 /* hue of particles of type 7 */
#define HUE_TYPE4 180.0 /* hue of particles of type 4 */
#define HUE_TYPE5 220.0 /* hue of particles of type 5 */
#define HUE_TYPE6 260.0 /* hue of particles of type 6 */
#define HUE_TYPE7 300.0 /* hue of particles of type 7 */
#define HUE_TYPE8 330.0 /* hue of particles of type 7 */
#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 BG_CHARGE_SLOPE 1.0 /* constant in BG_CHARGE backgound color scheme (default: 0.5) */
#define CHARGE_HUE_RANGE 0.5 /* range of charge colors */
#define PARTICLE_LMAX 1.5e4 /* angular momentum particle with brightest color */
#define RANDOM_RADIUS 0 /* set to 1 for random particle radius */
@ -295,36 +301,36 @@
#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 DT_PARTICLE 1.0e-6 /* time step for particle displacement */
#define KREPEL 50.0 /* constant in repelling force between particles */
#define EQUILIBRIUM_DIST 3.5 /* Lennard-Jones equilibrium distance */
#define EQUILIBRIUM_DIST_B 3.5 /* Lennard-Jones equilibrium distance for second type of particle */
#define KREPEL 40.0 /* constant in repelling force between particles */
#define EQUILIBRIUM_DIST 2.5 /* Lennard-Jones equilibrium distance */
#define EQUILIBRIUM_DIST_B 2.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 REPEL_RADIUS 25.0 /* radius in which repelling force acts (in units of particle radius) */
#define DAMPING 0.0 /* damping coefficient of particles */
#define INITIAL_DAMPING 1000.0 /* damping coefficient of particles during initial phase */
#define DAMPING_ROT 0.0 /* damping coefficient for rotation of particles */
#define DAMPING_PAIRS 0.0 /* damping between paired particles */
#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_ADD_MASS 1.0 /* mass of added particles */
#define PARTICLE_MASS 2.0 /* mass of particle of radius MU */
#define PARTICLE_MASS_B 2.0 /* mass of particle of radius MU_B */
#define PARTICLE_ADD_MASS 2.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_B 0.1 /* moment of inertia of second type of particle */
#define V_INITIAL 25.0 /* initial velocity range */
#define V_INITIAL_ADD 0.0 /* initial velocity range for added particles */
#define V_INITIAL 0.0 /* initial velocity range */
#define V_INITIAL_ADD 5500.0 /* initial velocity range for added particles */
#define OMEGA_INITIAL 100.0 /* initial angular velocity range */
#define VICSEK_VMIN 1.0 /* minimal speed of particles in Vicsek model */
#define VICSEK_VMAX 40.0 /* minimal speed of particles in Vicsek model */
#define 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 KCOULOMB_FACTOR 500.0 /* relative intensity of Coulomb interaction in I_COULOMB_LJ (default: 100.0) */
#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 THERMOSTAT 1 /* set to 1 to switch on thermostat */
#define THERMOSTAT 0 /* set to 1 to switch on thermostat */
#define VARY_THERMOSTAT 0 /* set to 1 for time-dependent thermostat schedule */
#define SIGMA 5.0 /* noise intensity in thermostat */
#define BETA 0.0002 /* initial inverse temperature */
#define BETA 0.00005 /* initial inverse temperature */
#define MU_XI 0.005 /* friction constant in thermostat */
#define KSPRING_BOUNDARY 2.0e11 /* confining harmonic potential outside simulation region */
#define KSPRING_OBSTACLE 5.0e11 /* harmonic potential of obstacles */
@ -376,7 +382,7 @@
#define WIND_FORCE 1.35e6 /* force of wind */
#define WIND_YMIN -0.6 /* min altitude of region with wind */
#define ROTATION 1 /* set to 1 to include rotation of particles */
#define ROTATION 0 /* set to 1 to include rotation of particles */
#define COUPLE_ANGLE_TO_THERMOSTAT 0 /* set to 1 to couple angular degrees of freedom to thermostat */
#define DIMENSION_FACTOR 0.25 /* scaling factor taking into account number of degrees of freedom */
#define KTORQUE 2.0e3 /* force constant in angular dynamics */
@ -392,12 +398,12 @@
#define QUADRUPOLE_RATIO 0.6 /* anisotropy in quadrupole potential */
#define INCREASE_BETA 0 /* set to 1 to increase BETA during simulation */
#define BETA_SCHEDULE 5 /* type of temperature schedule, see TS_* in global_ljones */
#define BETA_FACTOR 50.0 /* factor by which to change BETA during simulation */
#define BETA_SCHEDULE 0 /* type of temperature schedule, see TS_* in global_ljones */
#define BETA_FACTOR 50000.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 N_TOSCILLATIONS 1.5 /* number of temperature oscillations in BETA schedule */
#define N_TOSCILLATIONS 1.0 /* number of temperature oscillations in BETA schedule */
#define NO_OSCILLATION 0 /* set to 1 to have exponential BETA change only */
#define INITIAL_CONSTANT_PHASE 500 /* initial phase in which temperature is constant */
#define INITIAL_CONSTANT_PHASE 200 /* initial phase in which temperature is constant */
#define MIDDLE_CONSTANT_PHASE 0 /* middle 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 */
@ -437,24 +443,24 @@
#define MASS_PART_BOTTOM 10000.0 /* mass of particles at bottom */
#define NPART_BOTTOM 100 /* number of particles at the bottom */
#define ADD_PARTICLES 0 /* set to 1 to add particles */
#define ADD_REGION 1 /* shape of add regions, cf ADD_* in global_ljones */
#define ADD_TIME 0 /* time at which to add first particle */
#define ADD_PERIOD 6 /* time interval between adding further particles */
#define ADD_PARTICLES 1 /* set to 1 to add particles */
#define ADD_REGION 0 /* shape of add regions, cf ADD_* in global_ljones */
#define ADD_TIME 20 /* time at which to add first particle */
#define ADD_PERIOD 10000 /* 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 FINAL_NOADD_PERIOD 1800 /* final period where no particles are added */
#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 0 /* 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 TRACER_PARTICLE 1 /* set to 1 to have a tracer particle */
#define N_TRACER_PARTICLES 3000 /* number of tracer particles */
#define N_TRACER_PARTICLES 1000 /* number of tracer particles */
#define TRACER_STEPS 5 /* number of tracer steps recorded between images */
#define TRAJECTORY_LENGTH 7000 /* length of recorded trajectory */
#define TRAJECTORY_DRAW_LENGTH 250 /* length of drawn trajectory */
#define TRACER_LUM_FACTOR 100.0 /* controls luminosity decrease of trajectories with time */
#define TRAJECTORY_DRAW_LENGTH 1000 /* length of drawn trajectory */
#define TRACER_LUM_FACTOR 50.0 /* controls luminosity decrease of trajectories with time */
#define TRACER_PARTICLE_MASS 1.0 /* relative mass of tracer particle */
#define TRAJECTORY_WIDTH 2 /* width of tracer particle trajectory */
@ -506,28 +512,29 @@
#define POSITION_DEP_X -0.625 /* threshold value for position-dependent type */
#define PRINT_ENTROPY 0 /* set to 1 to compute entropy */
#define SPECIAL_IC 1 /* set to 1 for choosing specaial initial condition RD_INITIAL_COND */
#define REACTION_DIFFUSION 1 /* set to 1 to simulate a chemical reaction (particles may change type) */
#define SPECIAL_IC 0 /* set to 1 for choosing special initial condition RD_INITIAL_COND */
#define REACTION_DIFFUSION 0 /* set to 1 to simulate a chemical reaction (particles may change type) */
#define REACTION_MAX_TIME 100000 /* time after which no reactions take place */
#define RD_REACTION 11 /* type of reaction, see list in global_ljones.c */
#define RD_TYPES 5 /* number of types in reaction-diffusion equation */
#define RD_PLOT_TYPES 4 /* number of types shown in graph */
#define RD_INITIAL_COND 5 /* initial condition of particles */
#define REACTION_DIST 2.0 /* maximal distance for reaction to occur */
#define REACTION_PROB 0.5 /* probability controlling reaction term */
#define RD_REACTION 22 /* type of reaction, see list in global_ljones.c */
#define RD_TYPES 8 /* number of types in reaction-diffusion equation */
#define RD_PLOT_TYPES 8 /* number of types shown in graph */
#define RD_INITIAL_COND 2 /* initial condition of particles */
#define REACTION_DIST 2.8 /* maximal distance for reaction to occur */
#define REACTION_PROB 1.0 /* probability controlling reaction term */
#define DISSOCIATION_PROB 0.0 /* probability controlling dissociation reaction */
#define KILLING_PROB 0.0015 /* probability of enzymes being killed */
#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 EXOTHERMIC 1 /* set to 1 to make reaction exo/endothermic */
#define EXOTHERMIC 0 /* set to 1 to make reaction exo/endothermic */
#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 COLLISION_TIME 25 /* time during which collisions are shown */
#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 AGREGMAX 4 /* maximal number of partners for CHEM_AGGREGATION reaction */
#define NUDGE_FACTOR 0.0005 /* factor by which to correct particle distance */
#define COLLISION_TIME 35 /* time during which collisions are shown */
#define COLLISION_RADIUS 3.0 /* radius of discs showing collisions, in units of MU */
#define DELTAVMAX 500.0 /* maximal deltav allowed for pairing molecules */
#define AGREGMAX 3 /* maximal number of partners for CHEM_AGGREGATION reaction */
#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 NEUTRALIZE_REACTING_PARTICLES 1 /* set to 1 for reacting particles to become neutral */
#define CLUSTER_PARTICLES 0 /* set to 1 for particles to form rigid clusters */
#define CLUSTER_MAXSIZE 2 /* max size of clusters */
#define SMALL_CLUSTER_MAXSIZE 2 /* size limitation on smaller cluster */
@ -540,7 +547,7 @@
#define MU_RATIO 0.666666667 /* ratio by which to increase radius */
#define PRINT_PARTICLE_NUMBER 0 /* set to 1 to print total number of particles */
#define PLOT_PARTICLE_NUMBER 1 /* set to 1 to make of plot of particle number over time */
#define PLOT_PARTICLE_NUMBER 0 /* 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 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) */
@ -573,27 +580,33 @@
#define KSPRING_PAIRS 5.0e9 /* spring constant for pair interaction */
#define KTORQUE_PAIRS 1.0e10 /* constant for angular coupling in pair interaction */
#define KTORQUE_PAIR_ANGLE 0.0 /* constant for coupling between orientation in pairs */
#define NPARTNERS 1 /* number of partners of particles - for DNA, set NPARTNERS_DNA */
#define NPARTNERS 2 /* number of partners of particles - for DNA, set NPARTNERS_DNA */
#define NPARTNERS_DNA 8 /* number of partners of particles, case of DNA, should be at least 8 */
#define NARMS 5 /* number of "arms" for certain paring types */
#define PAIRING_TYPE 99 /* type of pairing, see POLY_ in global_ljones.c */
#define NARMS 4 /* number of "arms" for certain paring types */
#define PAIRING_TYPE 99 /* type of pairing, see POLY_ in global_ljones.c */
#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 MU_C 0.035 /* radius of partner particle */
#define MU_C 0.022 /* radius of 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 ALTERNATE_POLY_CHARGE 1 /* set to 1 for alternating charges in molecule */
#define ALTERNATE_POLY_CHARGE 0 /* set to 1 for alternating charges in molecule */
#define SECONDARY_PAIRING 0 /* set to 1 to pair with secondary partners, experimental */
#define DNA_RIGIDITY 0.5 /* controls rigidity for POLY_DNA_DOUBLE pairs, default = 1 */
#define PAIR_TYPEB_PARTICLES 0 /* set to 1 to pair particle of type 1 */
#define NPARTNERS_B 5 /* number of partners of particles */
#define NPARTNERS_B 18 /* number of partners of particles */
#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 MU_D 0.035 /* radius of partner particle */
#define PAIRING_TYPE_B 5 /* type of pairing, see POLY_ in global_ljones.c */
#define MU_D 0.022 /* radius of 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 ADD_ABSORBERS 1 /* set to 1 to add absorbing discs */
#define ABSORBER_PATTERN 4 /* pattern of absorbers, see AP_* in global_ljones */
#define ABSORBER_X 0.0
#define ABSORBER_Y 0.0 /* coordinates of first absorber */
#define ABSORBER_R 0.16 /* radius of absorber */
#define NXMAZE 16 /* width of maze */
#define NYMAZE 16 /* height of maze */
@ -621,25 +634,26 @@
/* 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 SIN_THETA_REG 0.01 /* regularization of sin(theta) for motion on sphere */
#define POLAR_PADDING 0.05 /* 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 FLIPX -1.0 /* set to -1 to flip left/right */
#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 ROTATE_ANGLE 750.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 */
double light[3] = {-0.40824829, 0.816496581, 0.40824829}; /* vector of "light" direction for P_3D_ANGLE color scheme */
double observer[3] = {-3.0, 0.0, 1.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))
@ -863,7 +877,8 @@ double temperature_schedule(int i)
// else if ((i < INITIAL_TIME + t3)&&(t3 > t2)) beta = BETA*factor2;
// else if ((i < INITIAL_TIME + t4)&&(t4 > t3))
// beta = BETA*exp(bexp2*(double)(i - INITIAL_TIME - t3) + 1.0e-10);
else beta = BETA*1.0e10;
else beta = BETA*factor2;
// else beta = BETA*1.0e10;
// else beta = BETA;
printf("beta = %.3lg\n", beta);
return(beta);
@ -1082,7 +1097,7 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NM
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;
int j, move, ncoup;
int j, move, ncoup, k, p;
if (initial_phase) damping = INITIAL_DAMPING;
else damping = DAMPING;
@ -1098,16 +1113,20 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NM
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];
// ffx = -2.0*ctheta*px[j]*py[j];
// ffy = stheta*ctheta*px[j]*px[j];
ffx = -2.0*ctheta*px[j]*py[j]*particle[j].mass_inv;
ffy = stheta*ctheta*px[j]*px[j]*particle[j].mass_inv;
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;
// ffx = -2.0*ctheta*particle[j].vx*particle[j].vy;
// ffy = stheta*ctheta*particle[j].vx*particle[j].vx;
ffx = -2.0*ctheta*particle[j].vx*particle[j].vy*particle[j].mass_inv;
ffy = stheta*ctheta*particle[j].vx*particle[j].vx*particle[j].mass_inv;
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;
@ -1252,6 +1271,8 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NM
// }
}
// wrap_particles(particle, molecule, cluster);
*ncoupled = ncoup;
return(totalenergy);
}
@ -1266,6 +1287,7 @@ double evolve_clusters(t_particle particle[NMAXCIRCLES], t_cluster cluster[NMAXC
double beta, int *nactive, int *nsuccess, int *nmove, int *ncoupled, int initial_phase, int verbose)
{
double a, totalenergy = 0.0, damping, direction, dmean, newx, newy, newangle, deltax, deltay, deltaangle;
double ctheta, stheta, stheta_reg, ffx, ffy;
static double b = 0.25*SIGMA*SIGMA*DT_PARTICLE/MU_XI, xi = 0.0;
int j, k, move, ncoup, mol, cl;
short int moved[NMAXCIRCLES];
@ -1277,16 +1299,41 @@ double evolve_clusters(t_particle particle[NMAXCIRCLES], t_cluster cluster[NMAXC
#pragma omp parallel for private(j,xi,totalenergy,a,move)
for (j=0; j<ncircles; j++) if (cluster[j].active)
{
cluster[j].vx = cpx[j] + 0.5*DT_PARTICLE*cluster[j].fx;
cluster[j].vy = cpy[j] + 0.5*DT_PARTICLE*cluster[j].fy;
cluster[j].omega = cpangle[j] + 0.5*DT_PARTICLE*cluster[j].torque;
{
if (SPHERE) /* add inertial terms for equation on a sphere */
{
ctheta = cos(cluster[j].xg);
stheta = sin(cluster[j].yg);
stheta_reg = stheta + SIN_THETA_REG;
ffx = -2.0*ctheta*cpx[j]*cpy[j];
ffy = stheta*ctheta*cpx[j]*cpx[j];
cluster[j].vx = cpx[j] + 0.5*DT_PARTICLE*(cluster[j].fx + ffx)/stheta_reg;
cluster[j].vy = cpy[j] + 0.5*DT_PARTICLE*(cluster[j].fy + ffy);
cluster[j].omega = cpangle[j] + 0.5*DT_PARTICLE*cluster[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*cluster[j].vx*cluster[j].vy;
ffy = stheta*ctheta*cluster[j].vx*cluster[j].vx;
cpx[j] = cluster[j].vx + 0.5*DT_PARTICLE*(cluster[j].fx + ffx)/stheta_reg;
cpy[j] = cluster[j].vy + 0.5*DT_PARTICLE*(cluster[j].fy + ffy);
cpangle[j] = cluster[j].omega + 0.5*DT_PARTICLE*cluster[j].torque;
// pangle[j] = particle[j].omega + 0.5*DT_PARTICLE*particle[j].torque/stheta_reg;
cluster[j].energy = (stheta*stheta*cpx[j]*cpx[j] + cpy[j]*cpy[j])*cluster[j].mass_inv;
}
else
{
cluster[j].vx = cpx[j] + 0.5*DT_PARTICLE*cluster[j].fx;
cluster[j].vy = cpy[j] + 0.5*DT_PARTICLE*cluster[j].fy;
cluster[j].omega = cpangle[j] + 0.5*DT_PARTICLE*cluster[j].torque;
cpx[j] = cluster[j].vx + 0.5*DT_PARTICLE*cluster[j].fx;
cpy[j] = cluster[j].vy + 0.5*DT_PARTICLE*cluster[j].fy;
cpangle[j] = cluster[j].omega + 0.5*DT_PARTICLE*cluster[j].torque;
cpx[j] = cluster[j].vx + 0.5*DT_PARTICLE*cluster[j].fx;
cpy[j] = cluster[j].vy + 0.5*DT_PARTICLE*cluster[j].fy;
cpangle[j] = cluster[j].omega + 0.5*DT_PARTICLE*cluster[j].torque;
cluster[j].energy = (cpx[j]*cpx[j] + cpy[j]*cpy[j])*cluster[j].mass_inv;
cluster[j].energy = (cpx[j]*cpx[j] + cpy[j]*cpy[j])*cluster[j].mass_inv;
}
// if (j == 110) printf("2-Cluster 110 is at (%.3lg, %.3lg)\n", cluster[110].xg, cluster[110].yg);
@ -1804,7 +1851,7 @@ void evolve_obstacles(t_obstacle obstacle[NMAXOBSTACLES])
void animation()
{
double time, scale, diss, rgb[3], dissip, gradient[2], x, y, dx, dy, dt, xleft, xright,
a, b, length, fx, fy, force[2], totalenergy = 0.0, pos[2], prop, vx, xi = 0.0, torque, torque_ij, pleft = 0.0, pright = 0.0, entropy[2], speed_ratio, xmin, xmax, ymin, ymax, delta_energy, speed, ratio = 1.0, ratioc, cum_etot = 0.0, emean = 0.0, radius_ratio, t, angle, theta, sum, alpha, bfield, track_x0, track_y0, efield, efieldy;
a, b, length, fx, fy, force[2], totalenergy = 0.0, pos[2], prop, vx, xi = 0.0, torque, torque_ij, pleft = 0.0, pright = 0.0, entropy[2], speed_ratio, xmin, xmax, ymin, ymax, delta_energy, speed, ratio = 1.0, ratioc, cum_etot = 0.0, emean = 0.0, radius_ratio, t, angle, theta, sum, alpha, bfield, track_x0, track_y0, efield, efieldy, dist;
double *qx, *qy, *px, *py, *qangle, *pangle, *pressure, *obstacle_speeds, *currents;
double *cqx, *cqy, *cpx, *cpy, *cqangle, *cpangle;
int i, j, k, n, m, s, ij[2], i0, iplus, iminus, j0, jplus, jminus, p, q, p1, q1, p2, q2, total_neighbours = 0, cl,
@ -1829,6 +1876,7 @@ void animation()
t_otriangle *otriangle;
t_ofacet *ofacet;
t_lj_sphere *wsphere;
t_absorber *absorber;
char message[100];
ratioc = 1.0 - ratio;
@ -1841,6 +1889,7 @@ void animation()
params.fboundary = 0.0;
params.gravity = GRAVITY;
params.radius = MU;
params.nabsorbed = 0;
particle = (t_particle *)malloc(NMAXCIRCLES*sizeof(t_particle)); /* particles */
@ -1900,6 +1949,13 @@ void animation()
for (i=0; i<2*NMAXCOLLISIONS; i++) collisions[i].time = 0;
}
if (ADD_ABSORBERS)
{
absorber = (t_absorber *)malloc(2*NMAX_ABSORBERS*sizeof(t_absorber));
nabsorbers = init_absorbers(absorber);
if (DRAW_SPHERE) draw_absorbers_sphere(absorber, wsphere);
}
if (SAVE_TIME_SERIES)
{
lj_time_series = fopen("lj_time_series.dat", "w");
@ -1917,6 +1973,10 @@ void animation()
group_speeds = (t_group_data *)malloc(ngroups*(INITIAL_TIME + NSTEPS)*sizeof(t_group_data));
}
/* test temperature schedule */
// for (i=0; i<INITIAL_TIME + NSTEPS; i++)
// fprintf(lj_log, "T(%i) = %.3lg\n", i, 1.0/temperature_schedule(i));
/* initialise array of trig functions to speed up drawing particles */
init_angles();
@ -1947,6 +2007,11 @@ void animation()
params.nactive = initialize_configuration(particle, hashgrid, obstacle, px, py, pangle, tracer_n, segment, molecule);
/* TEST - correcting bug */
// if (PAIRING_TYPE == POLY_PLUSMINUS)
// for (i=0; i<NMAXCIRCLES; i++)
// if (particle[i].charge != 0.0) particle[i].radius = MU_C;
printf("%i active particles\n", params.nactive);
printf("%i tracers\n", n_tracers);
@ -2048,7 +2113,6 @@ void animation()
pright = 0.0;
if (RECORD_PRESSURES) for (j=0; j<N_PRESSURES; j++) pressure[j] = 0.0;
printf("1\n");
// printf("evolving particles\n");
for(n=0; n<NVID; n++)
@ -2260,9 +2324,7 @@ void animation()
}
else
totalenergy = evolve_particles(particle, obstacle, hashgrid, qx, qy, qangle, px, py, pangle, params.beta, &params.nactive, &nsuccess, &nmove, &ncoupled, i < INITIAL_TIME);
/* TEST */
/* repair clusters */
if ((CLUSTER_PARTICLES)&&(REPAIR_CLUSTERS)) for (cl=0; cl<ncircles; cl++)
if ((cluster[cl].active)&&(cluster[cl].nparticles >= 2))
@ -2317,6 +2379,22 @@ void animation()
}
} /* end of for (n=0; n<NVID; n++) */
/* absorb particles */
if (ADD_ABSORBERS)
{
for (j=0; j<ncircles; j++) if (particle[j].active)
for (k=0; k<nabsorbers; k++)
{
if (SPHERE) dist = dist_sphere(particle[j].xc, particle[j].yc, absorber[k].xc, absorber[k].yc);
else dist = module2(particle[j].xc - absorber[k].xc, particle[j].yc - absorber[k].yc);
if (dist < absorber[k].radius)
{
particle[j].active = 0;
params.nabsorbed++;
}
}
}
if (TRACK_PARTICLE)
{
if (i < TRACK_INITIAL_TIME)
@ -2444,6 +2522,7 @@ void animation()
cum_etot += totalenergy;
}
printf("1\n");
printf("Boundary force: %.3f\n", params.fboundary/(double)(ncircles*NVID));
if (RESAMPLE_Y) printf("%i succesful moves out of %i trials\n", nsuccess, nmove);
if (INCREASE_GRAVITY) printf("Gravity: %.3f\n", params.gravity);
@ -2479,7 +2558,7 @@ void animation()
/* case of reaction-diffusion equation */
if ((i > INITIAL_TIME)&&(REACTION_DIFFUSION)&&(i < REACTION_MAX_TIME))
{
ncollisions = update_types(particle, molecule, cluster, collisions, ncollisions, particle_numbers, i - INITIAL_TIME - 1, &delta_energy);
ncollisions = update_types(particle, molecule, cluster, collisions, ncollisions, particle_numbers, i - INITIAL_TIME - 1, &delta_energy, tracer_n);
if (EXOTHERMIC) params.beta *= 1.0/(1.0 + delta_energy/totalenergy);
params.nactive = 0;
for (j=0; j<ncircles; j++) if (particle[j].active)
@ -2568,7 +2647,7 @@ void animation()
draw_frame(i, PLOT, BG_COLOR, ncollisions, traj_position, traj_length,
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, wsphere);
collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n, otriangle, ofacet, wsphere, absorber);
if (!((NO_EXTRA_BUFFER_SWAP)&&(MOVIE))) glutSwapBuffers();
@ -2600,7 +2679,7 @@ void animation()
{
draw_frame(i, PLOT_B, BG_COLOR_B, ncollisions, traj_position, traj_length,
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, wsphere);
collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n, otriangle, ofacet, wsphere, absorber);
glutSwapBuffers();
save_frame_lj_counter(NSTEPS + MID_FRAMES + 1 + counter);
counter++;
@ -2641,7 +2720,7 @@ void animation()
blank();
draw_frame(NSTEPS, PLOT, BG_COLOR, ncollisions, traj_position, traj_length,
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, wsphere);
collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n, otriangle, ofacet, wsphere, absorber);
}
if (DOUBLE_MOVIE) for (i=0; i<MID_FRAMES; i++)
{
@ -2653,7 +2732,7 @@ void animation()
{
draw_frame(NSTEPS, PLOT_B, BG_COLOR_B, ncollisions, traj_position, traj_length,
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, wsphere);
collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n, otriangle, ofacet, wsphere, absorber);
if (!((NO_EXTRA_BUFFER_SWAP)&&(MOVIE))) glutSwapBuffers();
}
if ((TIME_LAPSE)&&(!DOUBLE_MOVIE))
@ -2733,6 +2812,8 @@ void animation()
if (PAIR_PARTICLES) free(molecule);
if (ADD_ABSORBERS) free(absorber);
if (SAVE_TIME_SERIES)
{
fclose(lj_time_series);

777
sub_lj.c
View File

File diff suppressed because it is too large Load Diff

340
sub_lj_sphere.c
View File

@ -92,6 +92,64 @@ double dist_sphere(double phi1, double theta1, double phi2, double theta2)
}
void compute_midpoint_sphere(double phi1, double theta1, double phi2, double theta2, double *phi, double *theta)
/* compute midpoint between two points on the sphere */
{
double x1, y1, z1, x2, y2, z2, x3, y3, z3, r;
x1 = cos(phi1)*sin(theta1);
y1 = sin(phi1)*sin(theta1);
z1 = -cos(theta1);
x2 = cos(phi2)*sin(theta2);
y2 = sin(phi2)*sin(theta2);
z2 = -cos(theta2);
x3 = x1 + x2;
y3 = y1 + y2;
z3 = z1 + z2;
r = sqrt(x3*x3 + y3*y3 + z3*z3);
*theta = acos(-z3/r);
*phi = argument(x3, y3);
if (*phi < 0.0) *phi += DPI;
}
void compute_convex_combination_sphere(double phi1, double theta1, double phi2, double theta2, double *phia, double *thetaa, double *phib, double *thetab, double lambda)
/* compute convex combination of two points on the sphere */
{
double x1, y1, z1, x2, y2, z2, x3, y3, z3, r, lam1;
lam1 = 1.0 - lambda;
x1 = cos(phi1)*sin(theta1);
y1 = sin(phi1)*sin(theta1);
z1 = -cos(theta1);
x2 = cos(phi2)*sin(theta2);
y2 = sin(phi2)*sin(theta2);
z2 = -cos(theta2);
x3 = lambda*x1 + lam1*x2;
y3 = lambda*y1 + lam1*y2;
z3 = lambda*z1 + lam1*z2;
r = sqrt(x3*x3 + y3*y3 + z3*z3);
*thetaa = acos(-z3/r);
*phia = argument(x3, y3);
if (*phia < 0.0) *phia += DPI;
x3 = lam1*x1 + lambda*x2;
y3 = lam1*y1 + lambda*y2;
z3 = lam1*z1 + lambda*z2;
r = sqrt(x3*x3 + y3*y3 + z3*z3);
*thetab = acos(-z3/r);
*phib = argument(x3, y3);
if (*phib < 0.0) *phib += DPI;
}
int in_polygon(double x, double y, double r, int npoly, double apoly)
/* test whether (x,y) is in regular polygon of npoly sides inscribed in circle of radius r, turned by apoly Pi/2 */
{
@ -176,6 +234,7 @@ double type_hue(int type)
if (RD_REACTION == CHEM_BZ) return(HUE_TYPE2);
else return(HUE_TYPE7);
}
case (8): return(HUE_TYPE8);
default:
{
if (RD_REACTION == CHEM_BZ)
@ -350,7 +409,7 @@ void compute_particle_colors(t_particle particle, t_cluster cluster[NMAXCIRCLES]
}
case (P_CHARGE):
{
hue = (-0.6*tanh(particle.charge)+1.0)*180.0;
hue = (-CHARGE_HUE_RANGE*tanh(BG_CHARGE_SLOPE*particle.charge)+1.0)*180.0;
break;
}
case (P_MOL_ANGLE):
@ -631,7 +690,7 @@ void compute_all_particle_colors(t_particle particle[NMAXCIRCLES], t_cluster clu
}
void compute_background_color(t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NMAXOBSTACLES], int bg_color, t_hashgrid hashgrid[HASHX*HASHY])
void compute_background_color(t_particle particle[NMAXCIRCLES], t_segment segment[NMAXSEGMENTS], t_obstacle obstacle[NMAXOBSTACLES], int bg_color, t_hashgrid hashgrid[HASHX*HASHY])
/* color background according to particle properties */
{
int i, j, k, n, p, q, m, nnb, number, avrg_fact, obs;
@ -1139,12 +1198,19 @@ double dist_point_to_particle_sphere(int i, double phi, double psi, t_particle*
void init_3d() /* initialisation of window */
{
double width, height;
width = 2.0*(WINWIDTH/1760.0);
height = WINHEIGHT/990.0;
glLineWidth(3);
glClearColor(0.0, 0.0, 0.0, 1.0);
glClear(GL_COLOR_BUFFER_BIT);
glOrtho(-2.0, 2.0, -1.0, 1.0 , -1.0, 1.0);
// glOrtho(-2.0, 2.0, -1.0, 1.0 , -1.0, 1.0);
glOrtho(-width, width, -height, height, -1.0, 1.0);
}
void xyz_to_xy(double x, double y, double z, double xy_out[2])
@ -1177,7 +1243,7 @@ void xyz_to_xy(double x, double y, double z, double xy_out[2])
xinter[1] = t*observer[1] + (1.0-t)*y;
xinter[2] = t*observer[2] + (1.0-t)*z;
xy_out[0] = XSHIFT_3D + XY_SCALING_FACTOR*(xinter[0]*h[0] + xinter[1]*h[1]);
xy_out[0] = XSHIFT_3D + FLIPX*XY_SCALING_FACTOR*(xinter[0]*h[0] + xinter[1]*h[1]);
xy_out[1] = YSHIFT_3D + XY_SCALING_FACTOR*(xinter[0]*v[0] + xinter[1]*v[1] + xinter[2]*v[2]);
}
@ -1222,6 +1288,8 @@ void init_lj_sphere(t_lj_sphere *wsphere)
wsphere[i*NY_SPHERE+j].radius = 1.0;
wsphere[i*NY_SPHERE+j].cos_angle_sphere = wsphere[i*NY_SPHERE+j].x*light[0] + wsphere[i*NY_SPHERE+j].y*light[1] + wsphere[i*NY_SPHERE+j].z*light[2];
wsphere[i*NY_SPHERE+j].locked = 0;
}
}
}
@ -1296,8 +1364,8 @@ double test_distance(double phi, double psi, int i, t_particle* particle)
void add_particle_to_sphere(int i, int j, int part, t_particle particle[NMAXCIRCLES], t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE])
/* compute the effect at point (i,j) of adding a particle */
{
int draw_normal_particle = 1;
double x, y, r, dist, ca, sa, x1, y1, x2, y2, d1, d2, r2;
int draw_normal_particle = 1, k;
double x, y, r, dist, ca, sa, x1, y1, x2, y2, d1, d2, r2, omega, h;
static double dphi, dtheta;
static int first = 1;
@ -1335,7 +1403,7 @@ void add_particle_to_sphere(int i, int j, int part, t_particle particle[NMAXCIRC
r2 = 0.49*r*r;
if (d1 < r2)
if ((d1 < r2)&&(!wsphere[i*NY_SPHERE+j].locked))
{
wsphere[i*NY_SPHERE+j].r = particle[part].rgb[0];
wsphere[i*NY_SPHERE+j].g = particle[part].rgb[1];
@ -1365,7 +1433,7 @@ void add_particle_to_sphere(int i, int j, int part, t_particle particle[NMAXCIRC
r2 = 0.49*r*r;
if (d1 < r2)
if ((d1 < r2)&&(!wsphere[i*NY_SPHERE+j].locked))
{
wsphere[i*NY_SPHERE+j].r = particle[part].rgb[0];
wsphere[i*NY_SPHERE+j].g = particle[part].rgb[1];
@ -1395,7 +1463,7 @@ void add_particle_to_sphere(int i, int j, int part, t_particle particle[NMAXCIRC
r2 = 0.9*r*r;
if (d1 < r2)
if ((d1 < r2)&&(!wsphere[i*NY_SPHERE+j].locked))
{
wsphere[i*NY_SPHERE+j].r = particle[part].rgb[0];
wsphere[i*NY_SPHERE+j].g = particle[part].rgb[1];
@ -1421,7 +1489,7 @@ void add_particle_to_sphere(int i, int j, int part, t_particle particle[NMAXCIRC
r2 = 0.7*r*r;
if (d1 < r2)
if ((d1 < r2)&&(!wsphere[i*NY_SPHERE+j].locked))
{
wsphere[i*NY_SPHERE+j].r = particle[part].rgb[0];
wsphere[i*NY_SPHERE+j].g = particle[part].rgb[1];
@ -1433,6 +1501,50 @@ void add_particle_to_sphere(int i, int j, int part, t_particle particle[NMAXCIRC
}
break;
}
case (CHEM_POLYMER):
{
dist = dist_point_to_particle_sphere(part, x, y, particle, &ca, &sa);
x2 = dist*ca;
y2 = dist*sa;
r = 1.2*MU;
if ((dist < r)&&(!wsphere[i*NY_SPHERE+j].locked))
{
wsphere[i*NY_SPHERE+j].r = particle[part].rgb[0];
wsphere[i*NY_SPHERE+j].g = particle[part].rgb[1];
wsphere[i*NY_SPHERE+j].b = particle[part].rgb[2];
wsphere[i*NY_SPHERE+j].radius += 1.5*sqrt(r*r - dist*dist);
}
if (particle[part].type > 2)
{
omega = DPI/(double)(particle[part].type - 2);
for (k=0; k<particle[part].type-2; k++)
{
x1 = 1.4*r*cos(particle[part].angle + (double)k*omega);
y1 = 1.4*r*sin(particle[part].angle + (double)k*omega);
d1 = (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2);
r2 = 0.9*r*r;
if (d1 < r2)
{
h = 1.0 + 1.5*sqrt(r2 - d1);
if ((h > wsphere[i*NY_SPHERE+j].radius)&&(!wsphere[i*NY_SPHERE+j].locked))
{
wsphere[i*NY_SPHERE+j].r = particle[part].rgb[0];
wsphere[i*NY_SPHERE+j].g = particle[part].rgb[1];
wsphere[i*NY_SPHERE+j].b = particle[part].rgb[2];
wsphere[i*NY_SPHERE+j].radius = h;
}
}
}
}
draw_normal_particle = 0;
break;
}
}
}
@ -1442,19 +1554,23 @@ void add_particle_to_sphere(int i, int j, int part, t_particle particle[NMAXCIRC
if (dist < r)
{
wsphere[i*NY_SPHERE+j].r = particle[part].rgb[0];
wsphere[i*NY_SPHERE+j].g = particle[part].rgb[1];
wsphere[i*NY_SPHERE+j].b = particle[part].rgb[2];
wsphere[i*NY_SPHERE+j].radius += 1.5*sqrt(r*r - dist*dist);
h = 1.0 + 1.5*sqrt(r*r - dist*dist);
if ((h > wsphere[i*NY_SPHERE+j].radius)&&(!wsphere[i*NY_SPHERE+j].locked))
{
wsphere[i*NY_SPHERE+j].r = particle[part].rgb[0];
wsphere[i*NY_SPHERE+j].g = particle[part].rgb[1];
wsphere[i*NY_SPHERE+j].b = particle[part].rgb[2];
wsphere[i*NY_SPHERE+j].radius = h;
}
}
}
}
void draw_trajectory_sphere(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTICLES], int traj_position, int traj_length, t_particle *particle, t_cluster *cluster, t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE], int *tracer_n, int plot)
void draw_trajectory_sphere(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTICLES], t_hashgrid hashgrid[HASHX*HASHY], int traj_position, int traj_length, t_particle *particle, t_cluster *cluster, t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE], int *tracer_n, int plot)
/* draw tracer particle trajectory */
{
int i, j, i0, j0, time, p, q, width, imin, cell, i1, j1;
double x1, x2, y1, y2, rgb[3], rgbx[3], rgby[3], radius, lum, lum1;
double x1, x2, y1, y2, rgb[3], rgbx[3], rgby[3], radius, lum, lum1, rgb_bg[3];
static double dphi, dtheta;
static int first = 1;
@ -1489,6 +1605,17 @@ void draw_trajectory_sphere(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTI
if (lum < 0.0) lum = 0.0;
lum1 = 1.0 - lum;
if (COLOR_BACKGROUND)
{
cell = hash_cell(x1, y1);
if (!wsphere[cell].locked)
{
rgb_bg[0] = hashgrid[cell].r;
rgb_bg[1] = hashgrid[cell].g;
rgb_bg[2] = hashgrid[cell].b;
}
}
if ((x2 != x1)||(y2 != y1)) for (p=-1; p<2; p++)
for (q=-1; q<2; q++)
{
@ -1499,9 +1626,21 @@ void draw_trajectory_sphere(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTI
if (j1 < 0) j1 = 0;
if (j1 >= NY_SPHERE) j1 = NY_SPHERE-1;
cell = i1*NY_SPHERE+j1;
wsphere[cell].r = particle[tracer_n[j]].rgb[0]*lum + lum1;
wsphere[cell].g = particle[tracer_n[j]].rgb[1]*lum + lum1;
wsphere[cell].b = particle[tracer_n[j]].rgb[2]*lum + lum1;
if (!wsphere[cell].locked)
{
if (COLOR_BACKGROUND)
{
wsphere[cell].r = particle[tracer_n[j]].rgb[0]*lum + lum1*rgb_bg[0];
wsphere[cell].g = particle[tracer_n[j]].rgb[1]*lum + lum1*rgb_bg[1];
wsphere[cell].b = particle[tracer_n[j]].rgb[2]*lum + lum1*rgb_bg[2];
}
else
{
wsphere[cell].r = particle[tracer_n[j]].rgb[0]*lum + lum1;
wsphere[cell].g = particle[tracer_n[j]].rgb[1]*lum + lum1;
wsphere[cell].b = particle[tracer_n[j]].rgb[2]*lum + lum1;
}
}
}
}
else
@ -1522,6 +1661,17 @@ void draw_trajectory_sphere(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTI
if (lum < 0.0) lum = 0.0;
lum1 = 1.0 - lum;
if (COLOR_BACKGROUND)
{
cell = hash_cell(x1, y1);
if (!wsphere[cell].locked)
{
rgb_bg[0] = hashgrid[cell].r;
rgb_bg[1] = hashgrid[cell].g;
rgb_bg[2] = hashgrid[cell].b;
}
}
if ((x2 != x1)||(y2 != y1)) for (p=-1; p<2; p++)
for (q=-1; q<2; q++)
{
@ -1532,9 +1682,21 @@ void draw_trajectory_sphere(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTI
if (j1 < 0) j1 = 0;
if (j1 >= NY_SPHERE) j1 = NY_SPHERE-1;
cell = i1*NY_SPHERE+j1;
wsphere[cell].r = particle[tracer_n[j]].rgb[0]*lum + lum1;
wsphere[cell].g = particle[tracer_n[j]].rgb[1]*lum + lum1;
wsphere[cell].b = particle[tracer_n[j]].rgb[2]*lum + lum1;
if (!wsphere[cell].locked)
{
if (COLOR_BACKGROUND)
{
wsphere[cell].r = particle[tracer_n[j]].rgb[0]*lum + lum1*rgb_bg[0];
wsphere[cell].g = particle[tracer_n[j]].rgb[1]*lum + lum1*rgb_bg[1];
wsphere[cell].b = particle[tracer_n[j]].rgb[2]*lum + lum1*rgb_bg[2];
}
else
{
wsphere[cell].r = particle[tracer_n[j]].rgb[0]*lum + lum1;
wsphere[cell].g = particle[tracer_n[j]].rgb[1]*lum + lum1;
wsphere[cell].b = particle[tracer_n[j]].rgb[2]*lum + lum1;
}
}
}
}
for (i=imin; i < traj_position-1; i++)
@ -1550,6 +1712,17 @@ void draw_trajectory_sphere(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTI
if (lum < 0.0) lum = 0.0;
lum1 = 1.0 - lum;
if (COLOR_BACKGROUND)
{
cell = hash_cell(x1, y1);
if (!wsphere[cell].locked)
{
rgb_bg[0] = hashgrid[cell].r;
rgb_bg[1] = hashgrid[cell].g;
rgb_bg[2] = hashgrid[cell].b;
}
}
if ((x2 != x1)||(y2 != y1)) for (p=-1; p<2; p++)
for (q=-1; q<2; q++)
{
@ -1560,16 +1733,119 @@ void draw_trajectory_sphere(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTI
if (j1 < 0) j1 = 0;
if (j1 >= NY_SPHERE) j1 = NY_SPHERE-1;
cell = i1*NY_SPHERE+j1;
wsphere[cell].r = particle[tracer_n[j]].rgb[0]*lum + lum1;
wsphere[cell].g = particle[tracer_n[j]].rgb[1]*lum + lum1;
wsphere[cell].b = particle[tracer_n[j]].rgb[2]*lum + lum1;
if (!wsphere[cell].locked)
{
if (COLOR_BACKGROUND)
{
wsphere[cell].r = particle[tracer_n[j]].rgb[0]*lum + lum1*rgb_bg[0];
wsphere[cell].g = particle[tracer_n[j]].rgb[1]*lum + lum1*rgb_bg[1];
wsphere[cell].b = particle[tracer_n[j]].rgb[2]*lum + lum1*rgb_bg[2];
}
else
{
wsphere[cell].r = particle[tracer_n[j]].rgb[0]*lum + lum1;
wsphere[cell].g = particle[tracer_n[j]].rgb[1]*lum + lum1;
wsphere[cell].b = particle[tracer_n[j]].rgb[2]*lum + lum1;
}
}
}
}
}
}
void init_sphere_radius(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HASHX*HASHY], t_cluster cluster[NMAXCIRCLES], t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTICLES], int traj_position, int traj_length, t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE], int *tracer_n, int plot)
void draw_segments_sphere(t_segment segment[NMAXSEGMENTS], t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE])
/* draw the repelling segments on the sphere */
{
int s, i, cell, npoints, i0, j0, i1, j1, p, q;
double x1, y1, x2, y2, x, y, dt, length;
static double dphi, dtheta;
static int first = 1;
if (first)
{
dphi = DPI/(double)NX_SPHERE;
dtheta = PI/(double)NY_SPHERE;
first = 0;
}
for (s=0; s<nsegments; s++)
{
x1 = segment[s].x1;
y1 = segment[s].y1;
x2 = segment[s].x2;
y2 = segment[s].y2;
length = dist_sphere(x1, y1, x2, y2);
npoints = (int)(length*300.0);
dt = 1.0/(double)npoints;
for (i=0; i<npoints; i++)
{
x = x1 + (x2 - x1)*(double)i*dt;
y = y1 + (y2 - y1)*(double)i*dt;
i0 = (int)(x/dphi);
j0 = (int)(y/dtheta);
for (p=-1; p<2; p++)
for (q=-1; q<2; q++)
{
i1 = i0 + p;
if (i1 < 0) i1 = NX_SPHERE-1;
if (i1 >= NX_SPHERE) i1 = 0;
j1 = j0 + q;
if (j1 < 0) j1 = 0;
if (j1 >= NY_SPHERE) j1 = NY_SPHERE-1;
cell = i1*NY_SPHERE+j1;
wsphere[cell].r = 0.0;
wsphere[cell].g = 0.0;
wsphere[cell].b = 0.0;
wsphere[cell].radius += 0.005;
}
}
}
}
void draw_absorbers_sphere(t_absorber absorber[NMAX_ABSORBERS], t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE])
/* draw the absorbing discs on the sphere */
{
int i, j, n, cell;
double x, y, dist;
static double dphi, dtheta;
static int first = 1;
if (first)
{
dphi = DPI/(double)NX_SPHERE;
dtheta = PI/(double)NY_SPHERE;
first = 0;
}
for (i=0; i<NX_SPHERE; i++)
for (j=0; j<NY_SPHERE; j++)
for (n=0; n<nabsorbers; n++)
{
x = XMIN + (double)i*dphi;
y = YMIN + (double)j*dtheta;
dist = dist_sphere(x, y, absorber[n].xc, absorber[n].yc);
if (dist < absorber[n].radius)
{
cell = i*NY_SPHERE+j;
wsphere[cell].r = 0.0;
wsphere[cell].g = 0.0;
wsphere[cell].b = 0.0;
}
if (dist < absorber[n].radius*1.02)
{
cell = i*NY_SPHERE+j;
wsphere[cell].locked = 1;
wsphere[cell].radius = 0.95;
}
}
}
void init_sphere_radius(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HASHX*HASHY], t_cluster cluster[NMAXCIRCLES], t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTICLES], t_segment segment[NMAXSEGMENTS], int traj_position, int traj_length, t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE], int *tracer_n, int plot, int bg_color, t_absorber absorber[NMAX_ABSORBERS])
/* initialize sphere radius and colors from particles */
{
int i, j, n, m, i0, j0, i1, j1, p, q, part, width, deltai, imin, imax, deltaj, jmin, jmax, cell;
@ -1586,11 +1862,11 @@ void init_sphere_radius(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA
/* default radius and color */
#pragma omp parallel for private(i)
for (i=0; i<NX_SPHERE*NY_SPHERE; i++)
for (i=0; i<NX_SPHERE*NY_SPHERE; i++) if (!wsphere[i].locked)
{
wsphere[i].radius = 1.0;
if (COLOR_BACKGROUND)
if ((COLOR_BACKGROUND)&&(bg_color > 0))
{
cell = hash_cell(wsphere[i].phi, wsphere[i].theta);
wsphere[i].r = hashgrid[cell].r;
@ -1607,7 +1883,13 @@ void init_sphere_radius(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA
/* add tracer trajectories */
if (TRACER_PARTICLE)
draw_trajectory_sphere(trajectory, traj_position, traj_length, particle, cluster, wsphere, tracer_n, plot);
draw_trajectory_sphere(trajectory, hashgrid, traj_position, traj_length, particle, cluster, wsphere, tracer_n, plot);
if (ADD_FIXED_SEGMENTS)
draw_segments_sphere(segment, wsphere);
// if (ADD_ABSORBERS)
// draw_absorbers_sphere(absorber, wsphere);
/* draw zero meridian, for debugging */
// for (j=0; j<NY_SPHERE; j++)

790
sub_wave.c
View File

@ -1713,6 +1713,51 @@ int compute_star_coordinates(t_vertex polyline[NMAXPOLY])
return(NPOLY);
}
int compute_star_channel_coordinates(t_vertex polyline[NMAXPOLY])
/* compute positions of vertices of star-shaped domain with a channel */
{
int i;
double alpha, r, x, y, t, pos[2];
alpha = DPI/(double)NPOLY;
polyline[0].x = XMIN;
polyline[0].y = 0.5*WALL_WIDTH;
polyline[NPOLY+2].x = XMIN;
polyline[NPOLY+2].y = -0.5*WALL_WIDTH;
y = (LAMBDA-MU)*sin(alpha);
t = 0.5*WALL_WIDTH/y;
polyline[1].x = -LAMBDA + t*(LAMBDA - (LAMBDA - MU)*cos(alpha));
polyline[1].y = 0.5*WALL_WIDTH;
polyline[NPOLY+1].x = polyline[1].x;
polyline[NPOLY+1].y = -polyline[1].y;
for (i=1; i<NPOLY; i++)
{
if (i%2 == 1) r = LAMBDA - MU;
else r = LAMBDA;
x = -r*cos(alpha*(double)i);
y = r*sin(alpha*(double)i);
polyline[i+1].x = x;
polyline[i+1].y = y;
}
/* add origin to compute xy_in_billiard */
polyline[NPOLY+3].x = 0.0;
polyline[NPOLY+3].y = 0.0;
for (i=0; i<=NPOLY+3; i++)
{
xy_to_pos(polyline[i].x, polyline[i].y, pos);
polyline[i].posi = pos[0];
polyline[i].posj = pos[1];
}
return(NPOLY+3);
}
int compute_fresnel_coordinates(t_vertex polyline[NMAXPOLY])
/* compute positions of vertices approximating Fresnel lens */
{
@ -3511,6 +3556,71 @@ int compute_disc_hex_lattice_strip_coordinates(t_vertex polyline[NMAXPOLY], t_ar
return(nlines);
}
int compute_circle_pairs_coordinates(t_vertex polyline[NMAXPOLY], t_arc polyarc[NMAXPOLY], t_circle circle[NMAXCIRCLES], int *npolyarc, int *ncircles, short int top)
/* initialise lines for D_CIRCLE_PAIRS domain */
{
double dx, alpha, alphab, h1, h2, h1b, h2b, x1, y1, pos[2], mu, beta;
int i, j, k, nlines, narcs, ncirc;
nlines = 0;
narcs = 0;
ncirc = 0;
dx = (XMAX - XMIN)/(double)NGRIDX;
h1 = sqrt(MU*MU - 0.25*WALL_WIDTH*WALL_WIDTH);
h2 = LAMBDA - h1;
alpha = atan(0.5*WALL_WIDTH/h1);
h1b = sqrt(MU_B*MU_B - 0.25*WALL_WIDTH*WALL_WIDTH);
h2b = LAMBDA - h1b;
alphab = atan(0.5*WALL_WIDTH/h1b);
for (i=0; i<NGRIDX; i++)
{
x1 = XMIN + ((double)i+0.5)*dx;
y1 = LAMBDA;
polyarc[narcs].xc = x1;
polyarc[narcs].yc = y1;
polyarc[narcs].r = MU;
polyarc[narcs].angle1 = -PID + alpha;
polyarc[narcs].dangle = DPI - 2.0*alpha;
narcs++;
y1 = -LAMBDA;
polyarc[narcs].xc = x1;
polyarc[narcs].yc = y1;
polyarc[narcs].r = MU_B;
polyarc[narcs].angle1 = PID + alphab;
polyarc[narcs].dangle = DPI - 2.0*alphab;
narcs++;
polyline[nlines].x = x1 - 0.5*WALL_WIDTH;
polyline[nlines].y = h2;
nlines++;
polyline[nlines].x = x1 - 0.5*WALL_WIDTH;
polyline[nlines].y = -h2b;
nlines++;
polyline[nlines].x = x1 + 0.5*WALL_WIDTH;
polyline[nlines].y = h2;
nlines++;
polyline[nlines].x = x1 + 0.5*WALL_WIDTH;
polyline[nlines].y = -h2b;
nlines++;
}
for (i=0; i<nlines; i++)
{
xy_to_pos(polyline[i].x, polyline[i].y, pos);
polyline[i].posi = pos[0];
polyline[i].posj = pos[1];
}
*npolyarc = narcs;
*ncircles = ncirc;
return(nlines);
}
int compute_disc_honeycomb_lattice_coordinates(t_vertex polyline[NMAXPOLY], t_rect_rotated polyrect[NMAXPOLY], t_arc polyarc[NMAXPOLY], int *npolyrect_rot, int *npolyarc)
/* initialise lines for D_CIRCLE_LATTICE_HEX domain */
{
@ -4695,6 +4805,10 @@ int init_polyline(int depth, t_vertex polyline[NMAXPOLY])
{
return(compute_star_coordinates(polyline));
}
case (D_STAR_CHANNEL):
{
return(compute_star_channel_coordinates(polyline));
}
case (D_FRESNEL):
{
return(compute_fresnel_coordinates(polyline));
@ -4833,6 +4947,10 @@ int init_poly(int depth, t_vertex polyline[NMAXPOLY], t_rectangle polyrect[NMAXP
{
return(compute_star_coordinates(polyline));
}
case (D_STAR_CHANNEL):
{
return(compute_star_channel_coordinates(polyline));
}
case (D_FRESNEL):
{
return(compute_fresnel_coordinates(polyline));
@ -4938,6 +5056,10 @@ 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));
}
case (D_CIRCLE_PAIRS):
{
return(compute_circle_pairs_coordinates(polyline, polyarc, circles, npolyarc, ncircles, top));
}
case (D_TWO_PARABOLAS_ASYM_GUIDE):
{
polyrect[0].x1 = XMIN - WALL_WIDTH;
@ -4953,6 +5075,36 @@ int init_poly(int depth, t_vertex polyline[NMAXPOLY], t_rectangle polyrect[NMAXP
*npolyrect = 2;
return(0);
}
case (D_TWO_PARABOLAS_SEMITRANS):
{
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);
}
case (D_TWO_PARABOLAS_LENS):
{
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:
{
if ((ADD_POTENTIAL)&&(POTENTIAL == POT_MAZE))
@ -5223,6 +5375,12 @@ int init_xyin_from_image(short int * xy_in[NX])
nmaxpixels = 1070*601;
break;
}
case (IM_HAPPY_NEW_YEAR_TWOSIX):
{
image_file = fopen("Happy_New_Year_2026.ppm", "r");
nmaxpixels = 2073600;
break;
}
}
scan = fscanf(image_file,"%i %i\n", &nx, &ny);
@ -5284,6 +5442,14 @@ int init_xyin_from_image(short int * xy_in[NX])
sign = -1;
break;
}
case (IM_HAPPY_NEW_YEAR_TWOSIX):
{
rgbmax = 3*maxrgb/2;
scale = scalex;
sign = -1;
break;
}
}
@ -5324,7 +5490,26 @@ int init_xyin_from_image(short int * xy_in[NX])
/* to avoid reflection on left boundary */
if (IMAGE_FILE == IM_HAPPY_NEW_YEAR)
for (j=0; j<NY; j++) xy_in[0][j] = 1;
/* to avoid leaking on left boundary */
else if (IMAGE_FILE == IM_HAPPY_NEW_YEAR_TWOSIX)
{
for (j=0; j<NY; j++)
{
for (i=0; i<100; i++)
if (xy_in[i][j] == 1) xy_in[i][j] = 2;
for (i=NX-50; i<NX; i++)
if (xy_in[i][j] == 1) xy_in[i][j] = 2;
}
for (i=0; i<NX; i++)
{
for (j=0; j<200; j++)
if (xy_in[i][j] == 1) xy_in[i][j] = 2;
for (j=NY-200; j<NX; j++)
if (xy_in[i][j] == 1) xy_in[i][j] = 2;
}
}
fclose(image_file);
free(rgb_values);
return(1);
@ -5364,7 +5549,7 @@ void init_epicyloid()
int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_circle *circles)
/* returns 1 if (x,y) represents a point in the billiard */
{
double l2, r1, r2, r2mu, omega, b, c, d, angle, z, x1, y1, x2, y2, y3, u, v, u1, v1, dx, dy, width, alpha, s, a, r, height, ca, sa, l, ht, xshift, zz[9][2], x5, x6, f, fp, h1, cb, sb, c1, c2, nx, ny, phi;
double l2, r1, r2, r2mu, omega, b, c, d, angle, z, x1, y1, x2, y2, x3, y3, u, v, u1, v1, dx, dy, width, alpha, s, a, r, height, ca, sa, l, ht, xshift, zz[9][2], x5, x6, f, fp, h1, cb, sb, c1, c2, nx, ny, phi;
int i, j, k, k1, k2, condition = 0, m;
static int first = 1, nsides;
static double h, hh, ra, rb, ll, salpha;
@ -5570,6 +5755,111 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_
}
}
case (D_TWO_PARABOLAS_SEMITRANS):
{
if (x > 0.0)
{
if (vabs(y) > 0.5*LAMBDA) return(1);
x1 = -y*y/LAMBDA + 0.25*LAMBDA;
if ((x < x1)||(x > x1 + WALL_WIDTH)) return(1);
if (vabs(y) < 0.1*LAMBDA) return(0);
return(2);
}
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(2);
if (vabs(y) > 0.5*MU) return(1);
x1 = y*y/MU - 0.25*MU;
if ((x > x1)||(x < x1 - WALL_WIDTH)) return(1);
return(2);
}
}
case (D_TWO_PARABOLAS_LENS):
{
if (x > 0.0)
{
if (vabs(y) > 0.5*LAMBDA) return(1);
if (vabs(y) > 0.1*LAMBDA)
{
x1 = -y*y/LAMBDA + 0.25*LAMBDA;
if ((x < x1)||(x > x1 + WALL_WIDTH)) return(1);
return(2);
}
x1 = 0.24*LAMBDA;
x2 = x1 + WALL_WIDTH;
r2 = 1.5*0.25*LAMBDA;
r2 = r2*r2;
x3 = sqrt(r2 - 0.01*LAMBDA*LAMBDA) - x2;
// r2 = x2*x2 + 0.01*LAMBDA*LAMBDA;
// if ((x*x + y*y < r2)&&((x-x1-x2)*(x-x1-x2) + y*y < r2)) return(0);
if (((x+x3)*(x+x3) + y*y < r2)&&(x > x1)) return(0);
// if (vabs(y) < 0.1*LAMBDA) return(0);
return(1);
}
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(2);
if (vabs(y) > 0.5*MU) return(1);
x1 = y*y/MU - 0.25*MU;
if ((x > x1)||(x < x1 - WALL_WIDTH)) return(1);
return(2);
}
}
case (D_TWO_PARABOLAS_CLOSED_SEMITRANS):
{
if (vabs(y) > 0.75 + WALL_WIDTH) return(1);
x1 = 0.5*LAMBDA + 0.25*MU - y*y/MU;
if ((vabs(y) > 0.75)&&(vabs(x) < x1 + WALL_WIDTH)) return(2);
if ((vabs(x) < x1)||(vabs(x) > x1 + WALL_WIDTH)) return(1);
if ((x > 0.0)&&(vabs(y) < 0.2)) return(0);
return(2);
}
case (D_TWO_PARABOLAS_CLOSED_LENS):
{
if (vabs(y) > 0.75 + WALL_WIDTH) return(1);
if ((x > 0.0)&&(vabs(y) < 0.2))
{
/* lens */
y1 = 0.2;
x1 = 0.5*LAMBDA + 0.25*MU - y1*y1/MU;
x2 = x1 + WALL_WIDTH;
r2 = 1.5*0.25*LAMBDA;
r2 = r2*r2;
x3 = sqrt(r2 - y1*y1) - x2;
if (((x+x3)*(x+x3) + y*y < r2)&&(x > x1)) return(0);
return(1);
}
x1 = 0.5*LAMBDA + 0.25*MU - y*y/MU;
if ((vabs(y) > 0.75)&&(vabs(x) < x1 + WALL_WIDTH)) return(2);
if ((vabs(x) < x1)||(vabs(x) > x1 + WALL_WIDTH)) return(1);
return(2);
}
case (D_TWO_PARABOLAS_CLOSED_LENS_R):
{
if (vabs(y) > 0.75 + WALL_WIDTH) return(1);
if ((x > 0.0)&&(vabs(y) < WALL_WIDTH_B))
{
/* lens */
y1 = WALL_WIDTH_B;
x1 = 0.5*LAMBDA + 0.25*MU - y1*y1/MU;
x2 = x1 + WALL_WIDTH;
r2 = MU_B;
r2 = r2*r2;
x3 = sqrt(r2 - y1*y1) - x2;
if (((x+x3)*(x+x3) + y*y < r2)&&(x > x1)) return(0);
return(1);
}
x1 = 0.5*LAMBDA + 0.25*MU - y*y/MU;
if ((vabs(y) > 0.75)&&(vabs(x) < x1 + WALL_WIDTH)) return(2);
if ((vabs(x) < x1)||(vabs(x) > x1 + WALL_WIDTH)) return(1);
return(2);
}
case (D_FOUR_PARABOLAS):
{
x1 = MU + LAMBDA - 0.25*y*y/MU;
@ -5812,6 +6102,13 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_
condition += xy_in_triangle_tvertex(x, y, polyline[NPOLY], polyline[i], polyline[i+1]);
return(condition >= 1);
}
case (D_STAR_CHANNEL):
{
condition = ((x < 0.0)&&(vabs(y) < 0.5*WALL_WIDTH));
for (i = 0; i < NPOLY+3; i++)
condition += xy_in_triangle_tvertex(x, y, polyline[NPOLY+3], polyline[i+1], polyline[i]);
return(condition >= 1);
}
case (D_FRESNEL):
{
if (vabs(y) > 0.9*vabs(LAMBDA)) return(1);
@ -6352,6 +6649,45 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_
return(1);
}
case (D_TREE_BEAM):
{
/* upper triangle */
h = 1.5*LAMBDA;
r = 1.5*LAMBDA;
x1 = vabs(x);
y1 = y - h;
zz[0][0] = r*cos(PI/6.0); zz[0][1] = -0.5*r;
zz[1][0] = 0.0; zz[1][1] = r;
zz[2][0] = 0.0; zz[2][1] = -0.5*r;
/* middle triangle */
h = 0.25*LAMBDA;
r = 2.0*LAMBDA;
y2 = y - h;
zz[3][0] = r*cos(PI/6.0); zz[3][1] = -0.5*r;
zz[4][0] = 0.0; zz[4][1] = r;
zz[5][0] = 0.0; zz[5][1] = -0.5*r;
/* bottom triangle */
h = -1.5*LAMBDA;
r = 3.0*LAMBDA;
y3 = y - h;
zz[6][0] = r*cos(PI/6.0); zz[6][1] = -0.5*r;
zz[7][0] = 0.0; zz[7][1] = r;
zz[8][0] = 0.0; zz[8][1] = -0.5*r;
if (xy_in_triangle(x1, y1, zz[0], zz[1], zz[2])) return(1);
if (xy_in_triangle(x1, y2, zz[3], zz[4], zz[5])) return(1);
if (xy_in_triangle(x1, y3, zz[6], zz[7], zz[8])) return(1);
if ((x < 0.0)&&(y > OSCIL_YMID - 0.5*WALL_WIDTH)&&(y < OSCIL_YMID + 0.5*WALL_WIDTH))
return(1);
if ((x < -zz[6][0] + 0.5*WALL_WIDTH)&&(y > OSCIL_YMID - 0.75*WALL_WIDTH)&&(y < OSCIL_YMID + 0.75*WALL_WIDTH))
return(2);
return(0);
}
case (D_MICHELSON):
{
if ((vabs(x) > LAMBDA)&&(vabs(y) > LAMBDA)) return(2);
@ -6850,6 +7186,28 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_
if (vabs(y1) < 0.75*WALL_WIDTH) return(2);
return(0);
}
case (D_WAVEGUIDE_ELLIPSE_OBLIQUE):
{
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 - (x - XMIN)*tan(APOLY*PID);
if (vabs(y1) < 0.5*WALL_WIDTH) return(1);
if (vabs(y1) < 0.75*WALL_WIDTH) return(2);
return(0);
}
case (D_CIRCLE_PAIRS):
{
dx = (XMAX - XMIN)/(double)NGRIDX;
for (k=0; k<NGRIDX; k++)
{
x1 = XMIN + dx*((double)k + 0.5);
y1 = LAMBDA;
if (module2(x-x1, y-y1) < MU) return(1);
if (module2(x-x1, y+y1) < MU_B) return(1);
if ((vabs(y) < LAMBDA)&&(vabs(x-x1) < 0.5*WALL_WIDTH)) return(1);
}
return(0);
}
case (D_MENGER):
{
x1 = 0.5*(x+1.0);
@ -7139,7 +7497,7 @@ void draw_rc_hyp()
void draw_billiard(int fade, double fade_value) /* draws the billiard boundary */
{
double x0, y0, x, y, x1, y1, x2, y2, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, alpha2, dphi, omega, z, l, width, a, b, c, d, r1, r2, ymax, height, xmin, xmax, ca, sa, xshift, x5, x6, f, fp, xratio, w, nx, ny, x3, y3, psi, psi0, psi2, ca2, sa2, ca3, sa3;
double x0, y0, x, y, x1, y1, x2, y2, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, alpha2, dphi, omega, z, l, width, a, b, c, d, r1, r2, ymax, height, xmin, xmax, ca, sa, xshift, x5, x6, f, fp, xratio, w, nx, ny, x3, y3, psi, psi0, psi2, ca2, sa2, ca3, sa3, a1, b1, c1;
int i, j, k, k1, k2, mr2, ntiles;
static int first = 1, nsides;
static double h, hh, sqr3, ll, salpha, arcangle;
@ -7655,6 +8013,339 @@ void draw_billiard(int fade, double fade_value) /* draws the billiard bound
}
break;
}
case (D_TWO_PARABOLAS_SEMITRANS):
{
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 ();
y = 0.1*LAMBDA;
x = 0.25*LAMBDA - y*y/LAMBDA;
draw_line(x, y, x + WALL_WIDTH, y);
draw_line(x, -y, x + WALL_WIDTH, -y);
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_TWO_PARABOLAS_LENS):
{
dy = 0.4*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.1*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 ();
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.1*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 ();
y = 0.1*LAMBDA;
x = 0.25*LAMBDA - y*y/LAMBDA + WALL_WIDTH;
r1 = 0.25*1.5*LAMBDA;
alpha = asin(0.2666667); /* sin(alpha) = 0.1*LAMBDA/r1 */
x2 = r1*cos(alpha) - x,
// r1 = module2(x,y);
// alpha = atan(y/x);
glColor3f(1.0, 1.0, 1.0);
draw_circle_arc(0.0 - x2, 0.0, r1, -alpha, 2.0*alpha, NSEG);
// draw_circle_arc(0.0, 0.0, r1, -alpha, 2.0*alpha, NSEG);
// draw_circle_arc(2.0*x - WALL_WIDTH, 0.0, r1, PI-alpha, 2.0*alpha, NSEG);
x = 0.25*LAMBDA - y*y/LAMBDA;
draw_line(x,y,x,-y);
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_TWO_PARABOLAS_CLOSED_SEMITRANS):
{
dy = 1.5/(double)NSEG;
glBegin(GL_LINE_LOOP);
for (i = 0; i < NSEG+1; i++)
{
y = -0.75 + dy*(double)i;
x = 0.5*LAMBDA + 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.75 - dy*(double)i;
x = 0.5*LAMBDA + 0.25*MU - y*y/MU;
xy_to_pos(-x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
dy = (1.5 + 2.0*WALL_WIDTH)/(double)NSEG;
glBegin(GL_LINE_LOOP);
for (i = 0; i < NSEG+1; i++)
{
y = -0.75 - WALL_WIDTH + dy*(double)i;
x = WALL_WIDTH + 0.5*LAMBDA + 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.75 + WALL_WIDTH - dy*(double)i;
x = WALL_WIDTH + 0.5*LAMBDA + 0.25*MU - y*y/MU;
xy_to_pos(-x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
y = 0.2;
x = 0.5*LAMBDA + 0.25*MU - y*y/MU;
draw_line(x, y, x + WALL_WIDTH, y);
draw_line(x, -y, x + WALL_WIDTH, -y);
if (FOCI)
{
glColor3f(0.3, 0.3, 0.3);
draw_circle(0.5*LAMBDA, 0.0, r, NSEG);
draw_circle(-0.5*LAMBDA, 0.0, r, NSEG);
}
break;
}
case (D_TWO_PARABOLAS_CLOSED_LENS):
{
/* inner boundary */
dy = 1.5/(double)NSEG;
glBegin(GL_LINE_LOOP);
for (i = 0; i < NSEG+1; i++)
{
y = -0.75 + dy*(double)i;
x = 0.5*LAMBDA + 0.25*MU - y*y/MU;
xy_to_pos(-x, y, pos);
glVertex2d(pos[0], pos[1]);
}
dy = 0.55/(double)NSEG;
for (i = 0; i < NSEG+1; i++)
{
y = 0.75 - dy*(double)i;
x = 0.5*LAMBDA + 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.2 - dy*(double)i;
x = 0.5*LAMBDA + 0.25*MU - y*y/MU;
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
/* outer boundary */
glBegin(GL_LINE_STRIP);
dy = (0.55 + WALL_WIDTH)/(double)NSEG;
for (i = 0; i < NSEG+1; i++)
{
y = -0.2 - dy*(double)i;
x = WALL_WIDTH + 0.5*LAMBDA + 0.25*MU - y*y/MU;
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
dy = (1.5 + 2.0*WALL_WIDTH)/(double)NSEG;
for (i = 0; i < NSEG+1; i++)
{
y = -0.75 - WALL_WIDTH + dy*(double)i;
x = WALL_WIDTH + 0.5*LAMBDA + 0.25*MU - y*y/MU;
xy_to_pos(-x, y, pos);
glVertex2d(pos[0], pos[1]);
}
dy = (0.55 + WALL_WIDTH)/(double)NSEG;
for (i = 0; i < NSEG+1; i++)
{
y = 0.75 + WALL_WIDTH - dy*(double)i;
x = WALL_WIDTH + 0.5*LAMBDA + 0.25*MU - y*y/MU;
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
/* lens */
y = 0.2;
x = 0.5*LAMBDA + 0.25*MU - y*y/MU;
draw_line(x, y, x + WALL_WIDTH, y);
draw_line(x, -y, x + WALL_WIDTH, -y);
r1 = 0.25*1.5*LAMBDA;
alpha = asin(0.2/r1);
x1 = sqrt(r1*r1 - y*y) - x;
draw_circle_arc(-x1 + WALL_WIDTH, 0.0, r1, -alpha, 2.0*alpha, NSEG);
if (FOCI)
{
glColor3f(0.3, 0.3, 0.3);
draw_circle(0.5*LAMBDA, 0.0, r, NSEG);
draw_circle(-0.5*LAMBDA, 0.0, r, NSEG);
}
break;
}
case (D_TWO_PARABOLAS_CLOSED_LENS_R):
{
/* inner boundary */
dy = 1.5/(double)NSEG;
glBegin(GL_LINE_LOOP);
for (i = 0; i < NSEG+1; i++)
{
y = -0.75 + dy*(double)i;
x = 0.5*LAMBDA + 0.25*MU - y*y/MU;
xy_to_pos(-x, y, pos);
glVertex2d(pos[0], pos[1]);
}
dy = (0.75 - WALL_WIDTH_B)/(double)NSEG;
for (i = 0; i < NSEG+1; i++)
{
y = 0.75 - dy*(double)i;
x = 0.5*LAMBDA + 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 = -WALL_WIDTH_B - dy*(double)i;
x = 0.5*LAMBDA + 0.25*MU - y*y/MU;
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
/* outer boundary */
glBegin(GL_LINE_STRIP);
dy = (0.75 - WALL_WIDTH_B + WALL_WIDTH)/(double)NSEG;
for (i = 0; i < NSEG+1; i++)
{
y = -WALL_WIDTH_B - dy*(double)i;
x = WALL_WIDTH + 0.5*LAMBDA + 0.25*MU - y*y/MU;
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
dy = (1.5 + 2.0*WALL_WIDTH)/(double)NSEG;
for (i = 0; i < NSEG+1; i++)
{
y = -0.75 - WALL_WIDTH + dy*(double)i;
x = WALL_WIDTH + 0.5*LAMBDA + 0.25*MU - y*y/MU;
xy_to_pos(-x, y, pos);
glVertex2d(pos[0], pos[1]);
}
dy = (0.75 - WALL_WIDTH_B + WALL_WIDTH)/(double)NSEG;
for (i = 0; i < NSEG+1; i++)
{
y = 0.75 + WALL_WIDTH - dy*(double)i;
x = WALL_WIDTH + 0.5*LAMBDA + 0.25*MU - y*y/MU;
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
/* lens */
y = WALL_WIDTH_B;
x = 0.5*LAMBDA + 0.25*MU - y*y/MU;
draw_line(x, y, x + WALL_WIDTH, y);
draw_line(x, -y, x + WALL_WIDTH, -y);
r1 = MU_B;
alpha = asin(y/r1);
x1 = sqrt(r1*r1 - y*y) - x;
draw_circle_arc(-x1 + WALL_WIDTH, 0.0, r1, -alpha, 2.0*alpha, NSEG);
if (FOCI)
{
glColor3f(0.3, 0.3, 0.3);
draw_circle(0.5*LAMBDA, 0.0, r, NSEG);
draw_circle(-0.5*LAMBDA, 0.0, r, NSEG);
}
break;
}
case (D_FOUR_PARABOLAS):
{
x1 = 2.0*(sqrt(MU*(2.0*MU + LAMBDA)) - MU);
@ -8081,6 +8772,14 @@ void draw_billiard(int fade, double fade_value) /* draws the billiard bound
glEnd();
break;
}
case (D_STAR_CHANNEL):
{
glLineWidth(BOUNDARY_WIDTH);
glBegin(GL_LINE_LOOP);
for (i=0; i<npolyline; i++) tvertex_lineto(polyline[i]);
glEnd();
break;
}
case (D_FRESNEL):
{
glLineWidth(BOUNDARY_WIDTH);
@ -8960,6 +9659,11 @@ void draw_billiard(int fade, double fade_value) /* draws the billiard bound
/* do nothing */
break;
}
case (D_TREE_BEAM):
{
/* TODO - use polyline */
break;
}
case (D_MICHELSON):
{
w = 0.25*sqrt(2.0)*WALL_WIDTH;
@ -9507,6 +10211,64 @@ void draw_billiard(int fade, double fade_value) /* draws the billiard bound
}
break;
}
case (D_WAVEGUIDE_ELLIPSE_OBLIQUE):
{
/* lower waveguide boundary */
a = tan(PID*APOLY);
b = -1.0 + MU - XMIN*a;
a1 = a*a*LAMBDA*LAMBDA + 1.0;
b1 = a*b*LAMBDA*LAMBDA;
c1 = LAMBDA*LAMBDA*(b*b-1);
x1 = -(b1 + sqrt(b1*b1 - a1*c1))/a1;
y1 = a*x1 + b;
alpha = asin(y1);
draw_line(XMIN, -1.0 + MU, x1, y1);
/* upper waveguide boundary */
b += WALL_WIDTH;
b1 = a*b*LAMBDA*LAMBDA;
c1 = LAMBDA*LAMBDA*(b*b-1);
x1 = -(b1 + sqrt(b1*b1 - a1*c1))/a1;
y1 = a*x1 + b;
alpha2 = asin(y1);
draw_line(XMIN, -1.0 + MU + WALL_WIDTH, x1, y1);
draw_ellipse_arc(0.0, 0.0, LAMBDA, 1.0, PI - alpha, DPI - alpha2 + alpha, NSEG);
/* 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_PAIRS):
{
glBegin(GL_LINES);
for (i=0; i<npolyline; i++)
{
glVertex2d(polyline[i].posi, polyline[i].posj);
}
glEnd();
for (i=0; i<npolyarc; i++)
{
draw_circle_arc(polyarc[i].xc, polyarc[i].yc, polyarc[i].r, polyarc[i].angle1, polyarc[i].dangle, NSEG);
}
break;
}
case (D_MENGER):
{
glLineWidth(3);
@ -10875,6 +11637,8 @@ double oscillating_bc(int time, int j)
{
a = 0.0025;
t = (double)time*OMEGA;
if (OSCIL_LEFT_YSHIFT > 0.0)
t -= (double)j*OSCIL_LEFT_YSHIFT/(double)NY;
phase = t - a*t*t;
// if (time%1000 == 0) printf("time = %i, phase = %.4lg\n", time, phase);
return(AMPLITUDE*cos(phase)*exp(-phase*DAMPING));
@ -10882,6 +11646,8 @@ double oscillating_bc(int time, int j)
case (OSC_WAVE_PACKET):
{
t = (double)time*OMEGA;
if (OSCIL_LEFT_YSHIFT > 0.0)
t -= (double)j*OSCIL_LEFT_YSHIFT/(double)NY;
a = sqrt(INITIAL_VARIANCE)/OMEGA;
phase = AMPLITUDE*cos(t);
envelope = exp(-(t-INITIAL_SHIFT)*(t-INITIAL_SHIFT)/(a*a))*sqrt(DPI/a);
@ -10890,6 +11656,8 @@ double oscillating_bc(int time, int j)
case (OSC_WAVE_PACKETS):
{
t = (double)time*OMEGA;
if (OSCIL_LEFT_YSHIFT > 0.0)
t -= (double)j*OSCIL_LEFT_YSHIFT/(double)NY;
t1 = t - (double)((int)(t/WAVE_PACKET_SHIFT + 0.5))*WAVE_PACKET_SHIFT;
a = sqrt(INITIAL_VARIANCE)/OMEGA;
phase = AMPLITUDE*cos(t);
@ -10900,6 +11668,8 @@ double oscillating_bc(int time, int j)
{
// a = 0.25;
t = (double)time*OMEGA;
if (OSCIL_LEFT_YSHIFT > 0.0)
t -= (double)j*OSCIL_LEFT_YSHIFT/(double)NY;
phase = t + ACHIRP*t*t;
return(AMPLITUDE*sin(phase)*exp(-phase*DAMPING));
}
@ -10913,6 +11683,8 @@ double oscillating_bc(int time, int j)
else dist2 = 0.0;
dist2 = 500.0*dist2*dist2;
t = (double)time*OMEGA;
if (OSCIL_LEFT_YSHIFT > 0.0)
t -= (double)j*OSCIL_LEFT_YSHIFT/(double)NY;
amp = AMPLITUDE*cos(t)*exp(-(double)t*DAMPING);
amp *= exp(-dist2/INITIAL_VARIANCE);
return(amp);
@ -10922,6 +11694,8 @@ double oscillating_bc(int time, int j)
dist2 = (double)(j - NY/2)*(YMAX-YMIN)/(double)NY;
dist2 = dist2*dist2;
t = (double)time*OMEGA;
if (OSCIL_LEFT_YSHIFT > 0.0)
t -= (double)j*OSCIL_LEFT_YSHIFT/(double)NY;
amp = AMPLITUDE*cos(t)*exp(-(double)t*DAMPING);
amp *= exp(-dist2/INITIAL_VARIANCE);
return(amp);
@ -10937,6 +11711,8 @@ double oscillating_bc(int time, int j)
if (j > jmax) return(0.0);
if (j < jmin) return(0.0);
t = (double)time*OMEGA;
if (OSCIL_LEFT_YSHIFT > 0.0)
t -= (double)j*OSCIL_LEFT_YSHIFT/(double)NY;
amp = AMPLITUDE*cos(t)*exp(-(double)t*DAMPING);
amp *= cos(PID*dist2);
return(amp);
@ -10952,6 +11728,8 @@ double oscillating_bc(int time, int j)
if (j > jmax) return(0.0);
if (j < jmin) return(0.0);
t = (double)time*OMEGA;
if (OSCIL_LEFT_YSHIFT > 0.0)
t -= (double)j*OSCIL_LEFT_YSHIFT/(double)NY;
a = INITIAL_VARIANCE/(OMEGA*OMEGA);
amp = AMPLITUDE*cos(t)*exp(-(double)(t-INITIAL_SHIFT)*(t-INITIAL_SHIFT)/a);
amp *= cos(PID*dist2);
@ -10967,6 +11745,8 @@ double oscillating_bc(int time, int j)
else dist2 = 0.0;
dist2 = 500.0*dist2*dist2;
t = (double)time*OMEGA;
if (OSCIL_LEFT_YSHIFT > 0.0)
t -= (double)j*OSCIL_LEFT_YSHIFT/(double)NY;
amp = AMPLITUDE*(0.5*cos(t) + 0.5*cos(sqrt(7.0)*t))*exp(-(double)t*DAMPING);
amp *= exp(-dist2/INITIAL_VARIANCE);
return(amp);
@ -10982,6 +11762,8 @@ double oscillating_bc(int time, int j)
if (j > jmax) return(0.0);
if (j < jmin) return(0.0);
t = (double)time*OMEGA;
if (OSCIL_LEFT_YSHIFT > 0.0)
t -= (double)j*OSCIL_LEFT_YSHIFT/(double)NY;
amp = AMPLITUDE*(0.5*cos(t) + 0.5*cos(sqrt(7.0)*t))*exp(-(double)t*DAMPING);
amp *= cos(PID*dist2);
return(amp);
@ -10997,6 +11779,8 @@ double oscillating_bc(int time, int j)
if (j > jmax) return(0.0);
if (j < jmin) return(0.0);
t = (double)time*OMEGA;
if (OSCIL_LEFT_YSHIFT > 0.0)
t -= (double)j*OSCIL_LEFT_YSHIFT/(double)NY;
phase = t + ACHIRP*t*t;
// if (phase > DPI) phase = DPI;
amp = AMPLITUDE*sin(phase)*exp(-phase*DAMPING);

29
sub_wave_3d.c
View File

@ -1665,8 +1665,26 @@ void init_speed_dissipation(short int xy_in[NX*NY], double tc[NX*NY], double tcc
}
case (IOR_MICHELSON):
{
printf("Case IOR_MICHELSON in init_speed_dissipation of sub_wave_3d needs to be updated\n");
exit(1);
for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
tc[i*NY+j] = COURANT;
if (xy_in[i*NY+j] == 0)
{
tcc[i*NY+j] = courant2;
tgamma[i*NY+j] = GAMMA;
}
else if (xy_in[i*NY+j] == 2)
{
tcc[i*NY+j] = 0.0;
tgamma[i*NY+j] = 1.0;
}
else
{
tcc[i*NY+j] = courantb2;
tgamma[i*NY+j] = GAMMAB;
}
}
}
break;
}
case (IOR_GRADIENT_INDEX_LENS):
@ -1975,7 +1993,12 @@ void draw_wave_3d_ij(int i, int j, int movie, double phi[NX*NY], double psi[NX*N
if (NON_DIRICHLET_BC)
draw = (xy_in[i*NY+j])&&(xy_in[(i+1)*NY+j])&&(xy_in[i*NY+j+1])&&(xy_in[(i+1)*NY+j+1]);
else draw = (TWOSPEEDS)||(xy_in[i*NY+j]);
// else draw = (TWOSPEEDS)||(xy_in[i*NY+j]);
else
{
if (TWOSPEEDS) draw = (xy_in[i*NY+j] <= 1);
else draw = (xy_in[i*NY+j] == 1);
}
if (FLOOR_ZCOORD)
draw = (draw)&&(*wave[i*NY+j].p_zfield[movie] > zfloor)&&(*wave[(i+1)*NY+j].p_zfield[movie] > zfloor)&&(*wave[i*NY+j+1].p_zfield[movie] > zfloor)&&(*wave[(i+1)*NY+j+1].p_zfield[movie] > zfloor);

81
wave_billiard.c
View File

@ -48,7 +48,7 @@
#define SAVE_MEMORY 1 /* set to 1 to save memory when writing tiff images */
#define NO_EXTRA_BUFFER_SWAP 1 /* some OS require one less buffer swap when recording images */
#define VARIABLE_IOR 1 /* set to 1 for a variable index of refraction */
#define VARIABLE_IOR 0 /* set to 1 for a variable index of refraction */
#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 MANDEL_IOR_SCALE -0.05 /* parameter controlling dependence of IoR on Mandelbrot escape speed */
@ -62,8 +62,8 @@
#define XMIN -2.0
#define XMAX 2.0 /* x interval */
#define YMIN -1.197916667
#define YMAX 1.197916667 /* y interval for 9/16 aspect ratio */
#define YMIN -1.297916667
#define YMAX 1.097916667 /* y interval for 9/16 aspect ratio */
#define HIGHRES 1 /* set to 1 if resolution of grid is double that of displayed image */
#define HRES 1 /* dummy, only used by rde.c */
@ -72,7 +72,7 @@
/* Choice of the billiard table */
#define B_DOMAIN 982 /* choice of domain shape, see list in global_pdes.c */
#define B_DOMAIN 421 /* 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 IMAGE_FILE 5 /* for option D_IMAGE */
@ -87,11 +87,11 @@
#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 LAMBDA 1.4 /* 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 NPOLY 3 /* number of sides of polygon */
#define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */
#define LAMBDA 0.9 /* parameter controlling the dimensions of domain */
#define MU -0.2124612 /* parameter controlling the dimensions of domain */
#define MU_B 0.3 /* parameter controlling the dimensions of domain */
#define NPOLY 10 /* number of sides of polygon */
#define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */
#define MDEPTH 6 /* depth of computation of Menger gasket */
#define MRATIO 3 /* ratio defining Menger gasket */
#define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */
@ -99,8 +99,8 @@
#define FOCI 1 /* set to 1 to draw focal points of ellipse */
#define NGRIDX 6 /* number of grid point for grid of disks */
#define NGRIDY 14 /* number of grid point for grid of disks */
#define WALL_WIDTH 0.1 /* width of wall separating lenses */
#define WALL_WIDTH_B 0.05 /* width of wall separating lenses */
#define WALL_WIDTH 0.075 /* width of wall separating lenses */
#define WALL_WIDTH_B 0.1 /* width of wall separating lenses */
#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 WALL_WIDTH_ASYM 0.75 /* asymmetry of wall width (D_CIRCLE_LATTICE_NONISO) */
@ -121,22 +121,22 @@
/* xy_in_billiard and draw_billiard below */
/* Physical parameters of wave equation */
#define TWOSPEEDS 1 /* 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 1 /* set to 1 to add oscilating boundary condition on the left */
#define OSCILLATE_TOPBOT 0 /* set to 1 to enforce a planar wave on top and bottom boundary */
#define OSCILLATION_SCHEDULE 42 /* oscillation schedule, see list in global_pdes.c */
#define OSCIL_YMAX 0.05 /* defines oscilling beam range */
#define OSCIL_YMID -0.9 /* defines oscilling beam midpoint */
#define OSCIL_YMAX 0.0375 /* defines oscilling beam range */
#define OSCIL_YMID 0.0 /* defines oscilling beam midpoint */
#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 OMEGA 0.01 /* frequency of periodic excitation */
#define OMEGA 0.006 /* frequency of periodic excitation */
#define AMPLITUDE 2.0 /* amplitude of periodic excitation */
#define ACHIRP 0.25 /* acceleration coefficient in chirp */
#define DAMPING 0.0 /* damping of periodic excitation */
#define COURANT 0.25 /* Courant number in medium B */
#define COURANTB 0.1 /* Courant number */
#define GAMMA 5.0e-6 /* damping factor in wave equation */
#define COURANT 0.08 /* Courant number in medium B */
#define COURANTB 0.25 /* Courant number */
#define GAMMA 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_TOPBOT 1.0e-7 /* damping factor on boundary */
@ -150,11 +150,11 @@
/* For similar wave forms, COURANT^2*GAMMA should be kept constant */
#define ADD_OSCILLATING_SOURCE 0 /* set to 1 to add an oscillating wave source */
#define OSCILLATING_SOURCE_PERIOD 25.0 /* period of oscillating source */
#define OSCILLATING_SOURCE_PERIOD 6.0 /* period 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 ALTERNATE_SOURCE_PHASES 0 /* set to 1 to alternate initial phases of sources */
#define MAX_PULSING_TIME 750 /* max time for adding pulses */
#define MAX_PULSING_TIME 5000 /* max time for adding pulses */
#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 */
@ -169,10 +169,10 @@
/* Parameters for length and speed of simulation */
#define NSTEPS 2500 /* number of frames of movie */
#define NSTEPS 2700 /* number of frames of movie */
#define NVID 12 /* number of iterations between images displayed on screen */
#define NSEG 1000 /* number of segments of boundary */
#define INITIAL_TIME 150 /* 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 PRINT_SPEED 0 /* print speed of moving source */
#define PRINT_FREQUENCY 0 /* print frequency (for phased array) */
@ -199,8 +199,8 @@
/* Color schemes */
#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 18 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 18 /* Color palette, see list in global_pdes.c */
#define BLACK 1 /* background */
@ -208,17 +208,17 @@
#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 COLOR_RANGE 1.0 /* max range of color (default: 1.0) */
#define COLOR_RANGE 0.75 /* max range of color (default: 1.0) */
#define PHASE_FACTOR 1.0 /* factor in computation of phase in color scheme P_3D_PHASE */
#define PHASE_SHIFT 0.0 /* shift of phase in color scheme P_3D_PHASE */
#define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */
#define VSHIFT_AMPLITUDE -0.1 /* additional shift for wave amplitude */
#define VSHIFT_AMPLITUDE -0.5 /* additional shift for wave amplitude */
#define VSCALE_AMPLITUDE 1.0 /* additional scaling factor for wave amplitude */
#define E_SCALE 50.0 /* scaling factor for energy representation */
#define E_SCALE 300.0 /* scaling factor for energy representation */
#define LOG_SCALE 0.75 /* scaling factor for energy log representation */
#define LOG_SHIFT 0.75 /* shift of colors on log scale */
#define FLUX_SCALE 250.0 /* scaling factor for energy flux represtnation */
#define AVRG_E_FACTOR 0.75 /* controls time window size in P_AVERAGE_ENERGY scheme */
#define AVRG_E_FACTOR 0.95 /* controls time window size in P_AVERAGE_ENERGY scheme */
#define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */
#define FADE_IN_OBSTACLE 1 /* set to 1 to fade color inside obstacles */
#define SHADE_2D 1 /* set to 1 to add pseudo-3d shading effect */
@ -233,16 +233,16 @@
#define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */
#define COLORBAR_RANGE 1.7 /* scale of color scheme bar */
#define COLORBAR_RANGE_B 0.8 /* scale of color scheme bar for 2nd part */
#define COLORBAR_RANGE_B 4.5 /* scale of color scheme bar for 2nd part */
#define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */
#define CIRC_COLORBAR 0 /* set to 1 to draw circular color scheme */
#define CIRC_COLORBAR_B 0 /* set to 1 to draw circular color scheme */
#define DRAW_WAVE_PROFILE 0 /* set to 1 to draw a profile of the wave */
#define 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 DRAW_WAVE_PROFILE 1 /* set to 1 to draw a profile of the wave */
#define HORIZONTAL_WAVE_PROFILE 1 /* set to 1 to draw wave profile vertically */
#define VERTICAL_WAVE_PROFILE 0 /* set to 1 to draw wave profile vertically */
#define WAVE_PROFILE_X 1.9 /* value of x to sample wave profile */
#define WAVE_PROFILE_Y 0.12 /* value of y to sample wave profile */
#define WAVE_PROFILE_Y 0.0 /* value of y to sample wave profile */
#define PROFILE_AT_BOTTOM 1 /* draw wave profile at bottom instead of top */
#define AVERAGE_WAVE_PROFILE 0 /* set to 1 to draw time-average of wave profile squared*/
#define DRAW_WAVE_TIMESERIES 0 /* set to 1 to draw a time series of the wave, 2 to also draw it at the top */
@ -691,9 +691,9 @@ void animation()
/* initialise polyline and similar for drawing some domains */
npolyline = init_poly(MDEPTH, polyline, polyrect, polyrectrot, polyarc, circles, &npolyrect, &npolyrect_rot, &npolyarc, &ncircles, 1);
for (i=0; i<ncircles; i++) printf("circle %i: (%.3f, %.3f)\n", i, circles[i].xc, circles[i].yc);
// for (i=0; i<ncircles; i++) printf("circle %i: (%.3f, %.3f)\n", i, circles[i].xc, circles[i].yc);
// for (i=0; i<npolyline; i++) printf("vertex %i: (%.3f, %.3f)\n", i, polyline[i].x, polyline[i].y);
for (i=0; i<npolyline; i++) printf("vertex %i: (%.3f, %.3f)\n", i, polyline[i].x, polyline[i].y);
// for (i=0; i<npolyrect; i++) printf("polyrect vertex %i: (%.3f, %.3f) - (%.3f, %.3f)\n", i, polyrect[i].x1, polyrect[i].y1, polyrect[i].x2, polyrect[i].y2);
@ -878,10 +878,13 @@ void animation()
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, fade, fade_value);
/* add oscillating waves */
wave_source_x[0] = -1.0;
wave_source_y[0] = -0.5;
source_periods[0] = OSCILLATING_SOURCE_PERIOD;
source_amp[0] = INITIAL_AMP;
for (source = 0; source < N_SOURCES; source++)
{
wave_source_x[source] = polyarc[2*source].xc;
wave_source_y[source] = polyarc[2*source].yc;
source_periods[source] = OSCILLATING_SOURCE_PERIOD/(1.0 + 0.5*(double)source);
source_amp[source] = INITIAL_AMP;
}
for (source = 0; source < N_SOURCES; source++)
{
dperiod = source_periods[source];