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

Add files via upload

Browse Source 
Main changes to lennardjones.c, wave_billiard.c, wave_comparison.c, wave_3d.c
This commit is contained in:
Nils Berglund authored and GitHub committed 2026-08-10 18:49:37 +02:00
1 parent 840dc80aab
commit 849a9833ca
18 files changed
+2768 -346
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

global_ljones.c
+19
View File
@@ -89,6 +89,10 @@
#define S_CYLINDRICAL 100 /* two lines at top and bottom */
#define S_BOTTOM 110 /* one line at bottom */
#define S_VERTICAL 111 /* vertical double line */
#define S_FOUR_VERTICALS 112 /* four vertical double lines */
#define S_TWO_VERTICALS 113 /* four vertical double lines */
#define S_TWO_DIAGONALS 114 /* two slanted rectangles */
#define S_TWO_DIAGONALS_ROT 115 /* two slanted rectangles, rotated by 90° */
#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 */
@@ -107,6 +111,8 @@
#define S_DAM_WITH_HOLE_AND_RAMP 14 /* segments forming a dam in which a hole can open */
#define S_MAZE 15 /* segments forming a maze */
#define S_MAZE_DIAG 151 /* segments forming a maze with diagonally opposed exits */
#define S_MAZE_SPHERE 152 /* segments forming a maze on the sphere with fewer cells near poles */
#define S_MAZE_SPHERE_REG 153 /* segments forming a regular maze on the sphere */
#define S_EXT_RECTANGLE 16 /* particles outside a rectangle */
#define S_DAM_BRICKS 17 /* dam made of several bricks */
#define S_HLINE_HOLE 18 /* horizontal line with a hole in the bottom */
@@ -116,6 +122,8 @@
#define S_BIN_OPENING 20 /* bin containing particles opening at deactivation time */
#define S_BIN_LARGE 201 /* larger bin */
#define S_POLYGON_EXT 21 /* exterior of a regular polygon */
#define S_TWO_POLYGONS_EXT 211 /* exterior of a regular polygon */
#define S_TWO_POLYGONS_EXT_SHIFTED 212 /* exterior of a regular polygon */
#define S_WEDGE_EXT 22 /* exterior of a wedge */
#define S_MIXER 23 /* exterior of a blender made of rectangles */
#define S_AIRFOIL 24 /* exterior of an air foil */
@@ -404,6 +412,7 @@
#define BG_CURRENTX 9 /* background color depends on x-component of current */
#define BG_DIRECTION 10 /* background color depends on particle orientation */
#define BG_POTENTIAL 12 /* potential added to sphere */
#define BG_DEM 13 /* color given by DEM */
/* Obstacle color schemes */
@@ -437,6 +446,15 @@
#define SPP_WELLS 0 /* potential wells located at points in table wells[] */
#define SPP_CUBE 1 /* potential having shape of a cube */
#define SPP_PLANET 2 /* potential depending on planet DEM */
/* Type of planet */
#define PLANET_EARTH 0 /* DEM uses planet Earth */
#define PLANET_MARS 1 /* DEM uses planet Mars */
#define PLANET_MOON 2 /* DEM uses the Moon */
#define PLANET_VENUS 3 /* DEM uses planet Venus */
#define PLANET_MERCURY 4 /* DEM uses planet Mercury */
/* Type of sphere gravity */
@@ -749,6 +767,7 @@ typedef struct
double radius; /* radius with wave height */
double initial_radius; /* radius before particles are added */
double r, g, b; /* RGB values for image */
double r0, g0, b0; /* default RGB values */
double cos_angle; /* cosine of light angle */
double cos_angle_sphere; /* cosine of light angle for perfect sphere */
double cos_angle_pot; /* cosine of light angle for 2D potential */
global_pdes.c
+25 -1
View File
@@ -49,6 +49,7 @@
#define D_CIRCLES 20 /* several circles */
#define D_CIRCLES_IN_RECT 201 /* several circles in a rectangle */
#define D_CIRCLES_NEUMANN 202 /* several circles, with Dirichlet b.c. in top halr and Neumann in bottom half */
#define D_FOUR_PARABOLAS 31 /* four parabolas with axes in NSEW directions */
#define D_POLY_PARABOLAS 32 /* polygon with parabolic sides */
@@ -122,6 +123,10 @@
#define D_POLYCIRCLES_ANGLED 792 /* variant of D_POLYCIRCLES with angled small cavities */
#define D_TWOCIRCLES_ELLIPSE 793 /* ellipse with two adjacent circular cavities */
#define D_CIRCLES_ELLIPSE 794 /* ellipse with NPOLY adjacent circles */
#define D_TWO_ELLIPSES 795 /* two disjoint ellipses */
#define D_TWO_RECTANGLES 796 /* two disjoint rectangles */
#define D_TWO_POLYGONS 797 /* two disjoint polygons */
#define D_TWO_POLY_PARABOLAS 798 /* two disjoint parabolic resonators */
#define D_CARDIOID 90 /* cardioid */
#define D_NEPHROID 91 /* nephroid */
@@ -175,7 +180,7 @@
#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 NMAXCIRCLES 1000 /* 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) */
#define NMAXLINES 1000 /* maximal number of lines (for bounadries) */
#define NMAXSOURCES 30 /* maximal number of sources */
@@ -186,7 +191,9 @@
#define C_RAND_PERCOL 3 /* random percolation arrangement */
#define C_RAND_POISSON 4 /* random Poisson point process */
#define C_CLOAK 5 /* invisibility cloak */
#define C_CLOAK_OBJECT 51 /* invisibility cloak with an object inside */
#define C_CLOAK_A 6 /* first optimized invisibility cloak */
#define C_CLOAK_AB 61 /* first optimized invisibility cloak with smaller discs */
#define C_LASER 7 /* laser fight in a room of mirrors */
#define C_POISSON_DISC 8 /* Poisson disc sampling */
@@ -194,6 +201,7 @@
#define C_GOLDEN_SPIRAL 11 /* spiral pattern based on golden mean */
#define C_SQUARE_HEX 12 /* alternating between square and hexagonal/triangular */
#define C_HEX_NONUNIF 13 /* triangular grid with non-constant column distance */
#define C_SQUARE_TWO_SPACINGS 14 /* square grid but with two spacings */
#define C_RINGS 20 /* obstacles arranged in concentric rings */
#define C_RINGS_T 201 /* obstacles arranged in concentric rings, triangular lattice */
@@ -264,6 +272,7 @@
#define BC_ABSORBING 2 /* absorbing boundary conditions (beta version) */
#define BC_VPER_HABS 3 /* vertically periodic and horizontally absorbing boundary conditions */
#define BC_ABS_REFLECT 4 /* absorbing boundary conditions, except reflecting at y=0, for comparisons */
// #define BC_NEUMANN 5 /* Neumann boundary conditions */
#define BC_LSHAPE 10 /* L-shaped boundary conditions (surface of genus 2) */
// #define BC_OSCILL_ABSORB 5 /* oscillating boundary condition on the left, absorbing on other walls */
@@ -284,6 +293,7 @@
#define OSC_BEAM_SINE_TWOPERIODS 51 /* sum of two periodic oscillations modulated by y cut-off */
#define OSC_TWO_WAVES 6 /* two linear waves at an angle, separate */
#define OSC_TWO_WAVES_ADDED 61 /* two linear waves at an angle, superimposed */
#define OSC_TWOFREQ_TOPBOT 7 /* two different frequencies in top and bottom half */
/* Wave packet types */
@@ -505,6 +515,19 @@ typedef struct
int sign;
} t_wave_source;
typedef struct
{
double nx, ny; /* coordinates of normal vector */
short int signx1, signy1; /* relative position of neighboring grid points */
short int signx2, signy2; /* relative position of neighboring grid points */
int nbi1, nbj1, nbk1; /* absolute position of neighboring grid points */
int nbi2, nbj2, nbk2; /* absolute position of neighboring grid points */
int i, j, k; /* grid coordinates of boundary point */
short int top; /* indicates point is in top half, for comparisons */
short int active; /* allows to make problematic points (e.g. corners) inactive */
} t_neumann_bc;
#define N_NEUMANN_POINTS 80000 /* number of Neumann boundary points */
int ncircles = NMAXCIRCLES; /* actual number of circles, can be decreased e.g. for random patterns */
int npolyline = NMAXPOLY; /* actual length of polyline */
@@ -529,6 +552,7 @@ int npolyline_b = NMAXPOLY; /* actual length of polyline */
int npolyrect_b = NMAXPOLY; /* actual number of polyrect */
int npolyrect_rot_b = NMAXPOLY; /* actual number of rotated polyrect */
int npolyarc_b = NMAXPOLY; /* actual number of arcs */
int n_neumann_points = N_NEUMANN_POINTS; /* actual number of Neumann bc points */
t_circle circles_b[NMAXCIRCLES]; /* circular scatterers */
t_polygon polygons_b[NMAXCIRCLES]; /* polygonal scatterers */
t_vertex polyline_b[NMAXPOLY]; /* vertices of polygonal line */
heat.c
+1
View File
@@ -190,6 +190,7 @@
#define ACHIRP 0.2 /* acceleration coefficient in chirp */
#define DAMPING 0.0 /* damping of periodic excitation */
#define OMEGA 0.001 /* frequency of periodic excitation */
#define OMEGA_B 0.001 /* frequency of periodic excitation */
#define AMPLITUDE 0.8 /* amplitude of periodic excitation */
/* end of constants only used by wave_billiard and wave_3d */
lennardjones.c
+116 -77
View File
@@ -58,8 +58,8 @@
#define YMIN 0.0
#define YMAX 3.141592654 /* y interval for 9/16 aspect ratio */
#define INITXMIN 0.0
#define INITXMAX 6.28 /* x interval for initial condition */
#define INITXMIN 0.1
#define INITXMAX 6.18 /* x interval for initial condition */
#define INITYMIN 0.0
#define INITYMAX 3.14 /* y interval for initial condition */
@@ -77,8 +77,8 @@
#define BCXMIN 0.0
#define BCXMAX 6.283185307 /* x interval for boundary condition */
#define BCYMIN 0.3
#define BCYMAX 2.841592654 /* y interval for boundary condition */
#define BCYMIN 0.1
#define BCYMAX 3.041592654 /* y interval for boundary condition */
#define OBSXMIN -2.0
#define OBSXMAX 2.0 /* x interval for motion of obstacle */
@@ -108,7 +108,7 @@
#define COUPLE_MINLENGTH 0.5 /* length at which bonds decouple */
#define ADD_FIXED_SEGMENTS 0 /* set to 1 to add fixed segments as obstacles */
#define SEGMENT_PATTERN 111 /* pattern of repelling segments, see list in global_ljones.c */
#define SEGMENT_PATTERN 100 /* pattern of repelling segments, see list in global_ljones.c */
#define ROCKET_SHAPE 3 /* shape of rocket combustion chamber, see list in global_ljones.c */
#define ROCKET_SHAPE_B 3 /* shape of second rocket */
#define NOZZLE_SHAPE 6 /* shape of nozzle, see list in global_ljones.c */
@@ -119,31 +119,31 @@
#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.5 /* proportion of particles of first type */
#define TYPE_PROPORTION 0.3333333 /* proportion of particles of first type */
#define TWOTYPE_CONFIG 0 /* choice of types, see TTC_ list in global_ljones.c */
#define SYMMETRIZE_FORCE 1 /* set to 1 to symmetrize two-particle interaction, only needed if particles are not all the same */
#define CENTER_PX 0 /* set to 1 to center horizontal momentum */
#define CENTER_PY 0 /* set to 1 to center vertical momentum */
#define CENTER_PANGLE 0 /* set to 1 to center angular momentum */
#define INTERACTION 1 /* particle interaction, see list in global_ljones.c */
#define INTERACTION_B 1 /* particle interaction for second type of particle, see list in global_ljones.c */
#define SPIN_INTER_FREQUENCY 4.0 /* angular frequency of spin-spin interaction */
#define SPIN_INTER_FREQUENCY_B 4.0 /* angular frequency of spin-spin interaction for second particle type */
#define INTERACTION 171 /* particle interaction, see list in global_ljones.c */
#define INTERACTION_B 171 /* particle interaction for second type of particle, see list in global_ljones.c */
#define SPIN_INTER_FREQUENCY 6.0 /* angular frequency of spin-spin interaction */
#define SPIN_INTER_FREQUENCY_B 6.0 /* angular frequency of spin-spin interaction for second particle type */
#define MOL_ANGLE_FACTOR 1.0 /* rotation angle for P_MOL_ANGLE color scheme */
#define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */
#define NPOISSON 100 /* number of points for Poisson C_RAND_POISSON arrangement */
#define PDISC_DISTANCE 1.8 /* minimal distance in Poisson disc process, controls density of particles */
#define PDISC_DISTANCE 5.0 /* minimal distance in Poisson disc process, controls density of particles */
#define PDISC_CANDIDATES 100 /* number of candidates in construction of Poisson disc process */
#define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */
#define LAMBDA 0.2 /* parameter controlling the dimensions of domain */
#define MU 0.018 /* parameter controlling radius of particles */
#define MU_B 0.018 /* parameter controlling radius of particles of second type */
#define LAMBDA 0.75 /* 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 NPOLY 3 /* number of sides of polygon */
#define NPOLY 6 /* number of sides of polygon */
#define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */
#define AWEDGE 0.5 /* opening angle of wedge, in units of Pi/2 */
#define MDEPTH 4 /* depth of computation of Menger gasket */
@@ -171,8 +171,8 @@
/* Parameters for length and speed of simulation */
#define NSTEPS 1750 /* number of frames of movie */
#define NVID 100 /* number of iterations between images displayed on screen */
#define NSTEPS 2200 /* number of frames of movie */
#define NVID 80 /* 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 */
@@ -193,7 +193,7 @@
/* Plot type, see list in global_ljones.c */
#define PLOT 14
#define PLOT 4
#define PLOT_B 13 /* plot type for second movie */
/* Background color depending on particle properties */
@@ -225,7 +225,7 @@
#define COLOR_PALETTE 10 /* Color palette, see list in global_ljones.c */
#define COLOR_PALETTE_EKIN 10 /* Color palette for kinetic energy */
#define COLOR_PALETTE_ANGLE 17 /* Color palette for angle representation */
#define COLOR_PALETTE_DIRECTION 0 /* Color palette for direction representation */
#define COLOR_PALETTE_DIRECTION 10 /* Color palette for direction representation */
#define COLOR_PALETTE_INITIAL_POS 10 /* Color palette for initial position representation */
#define COLOR_PALETTE_DIFFNEIGH 10 /* Color palette for different neighbours representation */
#define COLOR_PALETTE_PRESSURE 11 /* Color palette for different neighbours representation */
@@ -235,6 +235,7 @@
#define COLOR_PALETTE_CLUSTER_SELECTED 11 /* Color palette for selected cluster representation */
#define COLOR_PALETTE_ANGULAR_MOMENTUM 17 /* Color palette for angular momentum */
#define COLOR_PALETTE_CURRENT 17 /* Color palette for current */
#define COLOR_PALETTE_POTENTIAL 11 /* Color palette for direction representation */
#define COLOR_HUE_CLUSTER_SELECTED 90.0 /* Color hue for selected cluster */
#define COLOR_HUE_CLUSTER_NOT_SELECTED 220.0 /* Color hue for selected cluster */
@@ -274,15 +275,15 @@
#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 4000.0 /* energy of particle with hottest color */
#define PARTICLE_EMAX 100000.0 /* energy of particle with hottest color */
#define PARTICLE_DMIN 200.0 /* energy of particle with largest local density */
#define PARTICLE_DMAX 500.0 /* energy of particle with largest local density */
#define SEGMENT_HUE_MIN 275.0 /* color of original segment */
#define SEGMENT_HUE_MAX 30.0 /* color of saturated segment */
#define OBSTACLE_EMAX 1000000.0 /* energy of obstacle with hottest color */
#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 40.0 /* hue of particles of type 1 */
#define HUE_TYPE0 280.0 /* hue of particles of type 0 */
#define HUE_TYPE1 70.0 /* hue of particles of type 1 */
#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 180.0 /* hue of particles of type 4 */
@@ -304,28 +305,28 @@
#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 2.0e-6 /* time step for particle displacement */
#define KREPEL 40.0 /* constant in repelling force between particles */
#define EQUILIBRIUM_DIST 3.5 /* Lennard-Jones equilibrium distance */
#define EQUILIBRIUM_DIST_B 4.5 /* Lennard-Jones equilibrium distance for second type of particle */
#define KREPEL 600.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 25.0 /* damping coefficient of particles */
#define DAMPING 400.0 /* damping coefficient of particles */
#define INITIAL_DAMPING 0.0 /* damping coefficient of particles during initial phase */
#define DAMPING_ROT 0.0 /* damping coefficient for rotation of particles */
#define DAMPING_ROT 50.0 /* damping coefficient for rotation of particles */
#define DAMPING_PAIRS 0.0 /* damping between paired 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 PARTICLE_INERTIA_MOMENT 0.05 /* moment of inertia of particle */
#define PARTICLE_INERTIA_MOMENT_B 0.05 /* moment of inertia of second type of particle */
#define V_INITIAL 50.0 /* initial velocity range */
#define V_INITIAL_ADD 4500.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 500.0 /* relative intensity of Coulomb interaction in I_COULOMB_LJ (default: 100.0) */
#define KCOULOMB_FACTOR 100.0 /* relative intensity of Coulomb interaction in I_COULOMB_LJ (default: 100.0) */
#define COULOMB_ALWAYS_REPEL 1 /* set to 1 to always repel with I_COULOMB_IMAGINARY */
#define OBSTACLE_DAMPING 0.0 /* damping of oscillating obstacles */
@@ -334,28 +335,28 @@
#define THERMOSTAT 1 /* set to 1 to switch on thermostat */
#define VARY_THERMOSTAT 0 /* set to 1 for time-dependent thermostat schedule */
#define SIGMA 5.0 /* noise intensity in thermostat */
#define BETA 0.002 /* initial inverse temperature */
#define BETA 0.02 /* 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_BOUNDARY 5.0e11 /* confining harmonic potential outside simulation region */
#define KSPRING_OBSTACLE 5.0e11 /* harmonic potential of obstacles */
#define NBH_DIST_FACTOR 6.0 /* radius in which to count neighbours */
#define BOND_DIST_FACTOR 6.0 /* radius in which to draw bonds */
#define GRAVITY 0.0 /* gravity acting on all particles */
#define GRAVITY 20000.0 /* gravity acting on all particles */
#define GRAVITY_X 0.0 /* horizontal gravity acting on all particles */
#define CIRCULAR_GRAVITY 0 /* set to 1 to have gravity directed to center */
#define SPHERE_GRAVITY 0 /* set to 1 to have gravity at constant angle wrt sphere */
#define INCREASE_GRAVITY 0 /* set to 1 to increase gravity during the simulation */
#define SPHERE_GRAVITY 1 /* set to 1 to have gravity at constant angle wrt sphere */
#define INCREASE_GRAVITY 1 /* set to 1 to increase gravity during the simulation */
#define GRAVITY_SCHEDULE 1 /* type of gravity schedule, see list in global_ljones.c */
#define GRAVITY_FACTOR 2000.0 /* factor by which to increase gravity */
#define GRAVITY_INITIAL_TIME 100 /* time at start of simulation with constant gravity */
#define GRAVITY_RESTORE_TIME 750 /* time at end of simulation with gravity restored to initial value */
#define GRAVITY_INITIAL_ANGLE 0.0 /* initial angle for SPHERE_GRAVITY */
#define GRAVITY_DELTA_ANGLE 1440.0 /* increase of angle for SPHERE_GRAVITY */
#define GRAVITY_FACTOR 10.0 /* factor by which to increase gravity */
#define GRAVITY_INITIAL_TIME 0 /* time at start of simulation with constant gravity */
#define GRAVITY_RESTORE_TIME 600 /* time at end of simulation with gravity restored to initial value */
#define GRAVITY_INITIAL_ANGLE 90.0 /* initial angle for SPHERE_GRAVITY */
#define GRAVITY_DELTA_ANGLE 0.0 /* increase of angle for SPHERE_GRAVITY */
#define KSPRING_VICSEK 0.2 /* spring constant for I_VICSEK_SPEED interaction */
#define VICSEK_REPULSION 10.0 /* repulsion between particles in Vicsek model */
#define ADD_EFIELD 0 /* set to 1 to add an electric field */
#define EFIELD 20000.0 /* value of electric field */
#define EFIELD 2500000.0 /* value of electric field */
#define EFIELD_Y 0.0 /* value of electric field */
#define ADD_BFIELD 0 /* set to 1 to add a magnetic field */
#define BFIELD 20000.0 /* value of magnetic field */
@@ -366,7 +367,7 @@
#define INCREASE_E 0 /* set to 1 to increase electric field */
#define OSCILLATE_E 0 /* set to 1 for oscillating electric field */
#define E_PERIOD 1000 /* period of oscillating electric field */
#define EFIELD_FACTOR 1000.0 /* factor by which to increase electric field */
#define EFIELD_FACTOR 2500.0 /* factor by which to increase electric field */
#define INCREASE_B 0 /* set to 1 to increase magnetic field */
#define BFIELD_FACTOR 1000.0 /* factor by which to increase magnetic field */
#define CHARGE_OBSTACLES 1 /* set to 1 for obstacles to be charged */
@@ -389,12 +390,12 @@
#define WIND_FORCE 1.35e6 /* force of wind */
#define WIND_YMIN -0.6 /* min altitude of region with wind */
#define ROTATE_SPHERE 1 /* set to 1 to add Coriolis and centripetal force */
#define OMEGA_SPHERE 4.0 /* angular frequency of rotating sphere */
#define CHANGE_OMEGA_SPHERE 1 /* set to 1 to change sphere rotation frequency */
#define OMEGA_SPHERE_FACTOR 200.0 /* change factor of sphere rotation frequency */
#define ROTATE_SPHERE 0 /* set to 1 to add Coriolis and centripetal force */
#define OMEGA_SPHERE 10.0 /* angular frequency of rotating sphere */
#define CHANGE_OMEGA_SPHERE 0 /* set to 1 to change sphere rotation frequency */
#define OMEGA_SPHERE_FACTOR 30.0 /* change factor of sphere rotation frequency */
#define ROTATION 0 /* set to 1 to include rotation of particles */
#define ROTATION 1 /* set to 1 to include rotation of particles */
#define COUPLE_ANGLE_TO_THERMOSTAT 0 /* set to 1 to couple angular degrees of freedom to thermostat */
#define DIMENSION_FACTOR 0.25 /* scaling factor taking into account number of degrees of freedom */
#define KTORQUE 2.0e3 /* force constant in angular dynamics */
@@ -449,7 +450,7 @@
#define PARTIAL_THERMO_RFIN 1.3 /* final radius of region without coupling */
#define INCREASE_KREPEL 0 /* set to 1 to increase KREPEL during simulation */
#define KREPEL_FACTOR 100.0 /* factor by which to change KREPEL during simulation */
#define KREPEL_FACTOR 0.005 /* factor by which to change KREPEL during simulation */
#define PART_AT_BOTTOM 0 /* set to 1 to include "seed" particles at bottom */
#define MASS_PART_BOTTOM 10000.0 /* mass of particles at bottom */
@@ -467,12 +468,12 @@
#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 0 /* set to 1 to have a tracer particle */
#define N_TRACER_PARTICLES 5500 /* number of tracer particles */
#define TRACER_PARTICLE 1 /* set to 1 to have a tracer particle */
#define N_TRACER_PARTICLES 1200 /* 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_LENGTH 4000 /* length of recorded trajectory */
#define TRAJECTORY_DRAW_LENGTH 1000 /* length of drawn trajectory */
#define TRACER_LUM_FACTOR 50.0 /* controls luminosity decrease of trajectories with time */
#define TRACER_LUM_FACTOR 100.0 /* controls luminosity decrease of trajectories with time */
#define TRACER_PARTICLE_MASS 1.0 /* relative mass of tracer particle */
#define TRAJECTORY_WIDTH 2 /* width of tracer particle trajectory */
@@ -580,10 +581,10 @@
#define WALL_TIME 0 /* time during which to keep wall */
#define CHANGE_TYPES 0 /* set to 1 to change type proportion in course of simulation */
#define PROP_MIN 0.0 /* min proportion of type 1 particles */
#define PROP_MAX 1.0 /* max proportion of type 1 particles */
#define PROP_TINITIAL 250 /* initial time without change */
#define PROP_TFINAL 250 /* final time without change */
#define PROP_MIN 1.0 /* min proportion of type 1 particles */
#define PROP_MAX 0.0 /* max proportion of type 1 particles */
#define PROP_TINITIAL 50 /* initial time without change */
#define PROP_TFINAL 50 /* final time without change */
#define PAIR_PARTICLES 0 /* set to 1 to form particles pairs */
#define RANDOMIZE_ANGLE 0 /* set to 1 for random orientation */
@@ -625,17 +626,18 @@
#define ADD_POTENTIAL_SPHERE 0 /* add potential for gradient force field on sphere */
#define DRAW_POTENTIAL_SPHERE 1 /* draw sphere radius depending on potential */
#define SPHERE_POTENTIAL 1 /* type of sphere potential */
#define SPHERE_POTENTIAL 2 /* type of sphere potential */
#define SPHERE_POT_PATTERN 3 /* pattern of local minma of SPP_WELLS sphere potential */
#define PLANET_DEM 4 /* planet DEM used for SPP_PLANET */
#define POT_SPHERE_AMP 1.0 /* amplitude in definition of potential on sphere */
#define POT_SPHERE_RADIUS 0.1 /* radius in definition of potential on sphere */
#define POT_SPHERE_SMOOTH 0.5 /* smoothing of potential well */
#define POT_SPHERE_STRENGTH 1.0e3 /* coefficient of gradient force */
#define POT_SPHERE_STRENGTH 2.5e4 /* coefficient of gradient force */
#define NXMAZE 16 /* width of maze */
#define NYMAZE 16 /* height of maze */
#define NXMAZE 18 /* width of maze */
#define NYMAZE 10 /* height of maze */
#define MAZE_MAX_NGBH 4 /* max number of neighbours of maze cell */
#define RAND_SHIFT 5 /* seed of random number generator */
#define RAND_SHIFT 53 /* seed of random number generator */
#define MAZE_XSHIFT 0.0 /* horizontal shift of maze */
#define MAZE_WIDTH 0.015 /* width of maze walls */
@@ -643,9 +645,11 @@
#define FMAX 1.0e8 /* maximal force */
#define FLOOR_OMEGA 0 /* set to 1 to limit particle momentum to PMAX */
#define PMAX 1000.0 /* maximal force */
#define FLOOR_ENERGY 0 /* set to 1 to limit mean energy to EMAX, use sparingly! */
#define EMAX 250.0 /* maximal mean energy */
#define HASHX 120 /* size of hashgrid in x direction */
#define HASHY 60 /* size of hashgrid in y direction */
#define HASHX 100 /* size of hashgrid in x direction */
#define HASHY 50 /* size of hashgrid in y direction */
#define HASHMAX 100 /* maximal number of particles per hashgrid cell */
#define HASHGRID_PADDING 0.1 /* padding of hashgrid outside simulation window */
@@ -662,24 +666,25 @@
#define POLAR_PADDING 0.01 /* region around poles that belong to the same hashcell */
#define DRAW_SPHERE 1 /* set to 1 to draw 3D sphere */
#define DRAW_ELLIPSES_ON_SPHERE 1 /* set to 1 to draw ellipses instead of circles on sphere in 2D */
#define NX_SPHERE 3000
#define NY_SPHERE 1500 /* number of points on sphere */
#define NX_SPHERE 1800
#define NY_SPHERE 1350 /* 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 2.0 /* overall scaling factor for on-screen (x,y) coordinates after projection */
#define XY_SCALING_FACTOR 1.9 /* overall scaling factor for on-screen (x,y) coordinates after projection */
#define 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 RSCALE_POTENTIAL 1.0 /* radial scaling of potential */
#define RSCALE_POTENTIAL 0.15 /* radial scaling of potential */
#define ROTATE_VIEW 1 /* set to 1 to rotate position of observer */
#define ROTATE_ANGLE 360.0 /* total angle of rotation during simulation */
#define VIEWPOINT_TRAJ 1 /* type of viewpoint trajectory */
#define ROTATE_ANGLE 45.0 /* total angle of rotation during simulation */
#define VIEWPOINT_TRAJ 2 /* type of viewpoint trajectory */
#define MAX_LATITUDE 45.0 /* maximal latitude for viewpoint trajectory VP_ORBIT2 */
#define DRAW_POLAR_AXIS 1 /* set to 1 to draw polar axis */
double light[3] = {-0.40824829, 0.816496581, 0.40824829}; /* vector of "light" direction for P_3D_ANGLE color scheme */
double observer[3] = {-3.0, 2.0, 2.5}; /* location of observer for REP_PROJ_3D representation */
double light[3] = {-0.816496581, 0.40824829, 0.40824829}; /* vector of "light" direction for P_3D_ANGLE color scheme */
double observer[3] = {-2.0, 0.3, -1.0}; /* location of observer for REP_PROJ_3D representation */
#define NO_WRAP_BC ((BOUNDARY_COND != BC_PERIODIC)&&(BOUNDARY_COND != BC_PERIODIC_CIRCLE)&&(BOUNDARY_COND != BC_PERIODIC_TRIANGLE)&&(BOUNDARY_COND != BC_KLEIN)&&(BOUNDARY_COND != BC_PERIODIC_FUNNEL)&&(BOUNDARY_COND != BC_BOY)&&(BOUNDARY_COND != BC_GENUS_TWO))
#define PERIODIC_BC ((BOUNDARY_COND == BC_PERIODIC)||(BOUNDARY_COND == BC_PERIODIC_CIRCLE)||(BOUNDARY_COND == BC_PERIODIC_FUNNEL)||(BOUNDARY_COND == BC_PERIODIC_TRIANGLE))
@@ -1132,7 +1137,7 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NM
double beta, int *nactive, int *nsuccess, int *nmove, int *ncoupled, int initial_phase)
{
double a, totalenergy = 0.0, damping, damping1, damping_rot1, direction, dmean, dratio;
double ctheta, stheta, stheta_reg, ffx, ffy;
double ctheta, stheta, stheta_reg, ffx, ffy, vratio, vratio1, energy;
static double b = 0.25*SIGMA*SIGMA*DT_PARTICLE/MU_XI, xi = 0.0;
int j, move, ncoup, k, p;
@@ -1311,6 +1316,27 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NM
// wrap_particles(particle, molecule, cluster);
*ncoupled = ncoup;
/* TEST */
if ((FLOOR_ENERGY)&&(totalenergy > (double)*ncoupled*EMAX))
{
// printf("Flooring total energy\n");
vratio = sqrt((double)*ncoupled*EMAX/totalenergy);
for (p=0; p<ncircles; p++)
{
vratio1 = vratio;
energy = particle[p].energy;
if (energy > 2.0*EMAX) vratio1 = sqrt(2.0*EMAX/energy);
if (vabs(particle[p].energy) > 0.1*EMAX)
{
particle[p].vx *= vratio1;
particle[p].vy *= vratio1;
px[p] *= vratio1;
py[p] *= vratio1;
}
}
}
return(totalenergy);
}
@@ -1887,8 +1913,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, dist, dphi_inv, dtheta_inv, sincos;
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, dist, dphi_inv, dtheta_inv, sincos, mass, vratio;
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,
@@ -2375,10 +2400,11 @@ void animation()
{
/* centrifugal force */
sincos = cos(particle[j].yc)*sin(particle[j].yc);
particle[j].fy += params.omega_sphere*params.omega_sphere*sincos;
mass = 1.0/particle[j].mass_inv;
particle[j].fy += params.omega_sphere*params.omega_sphere*sincos*mass;
/* Coriolis force */
particle[j].fx += 2.0*params.omega_sphere*cos(particle[j].yc)*particle[j].vy;
particle[j].fy += 2.0*params.omega_sphere*sincos*particle[j].vx;
particle[j].fx += 2.0*params.omega_sphere*cos(particle[j].yc)*particle[j].vy*mass;
particle[j].fy += 2.0*params.omega_sphere*sincos*particle[j].vx*mass;
}
/* add potential force on sphere */
@@ -2596,6 +2622,19 @@ void animation()
printf("Mean kinetic energy: %.3f\n", totalenergy/(double)ncircles);
printf("Kinetic energy by coupled particle: %.3f\n", totalenergy/(double)ncoupled);
// if ((FLOOR_ENERGY)&&(totalenergy > (double)ncoupled*EMAX))
// {
// printf("Flooring total energy\n");
// vratio = sqrt((double)ncoupled*EMAX/totalenergy);
// for (p=0; p<ncircles; p++)
// {
// particle[p].vx *= vratio;
// particle[p].vy *= vratio;
// px[p] *= vratio;
// py[p] *= vratio;
// }
// }
if ((!THERMOSTAT)&&(LIMIT_ENERGY))
{
if (cum_etot > 0.0)
mangrove.c
+1
View File
@@ -109,6 +109,7 @@
#define X_SHIFT -0.9 /* x range on which to apply OSCILLATE_TOPBOT */
#define OMEGA 0.00133333333 /* frequency of periodic excitation */
#define OMEGA_B 0.001 /* frequency of periodic excitation */
#define K_BC 3.0 /* spatial period of periodic excitation in y direction */
#define KX_BC 10.0 /* spatial period of periodic excitation in x direction */
#define KY_BC 3.3333 /* spatial period of periodic excitation in y direction */
rde.c
+27 -22
View File
@@ -40,13 +40,14 @@
#include <time.h>
#define MOVIE 0 /* set to 1 to generate movie */
#define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */
#define DOUBLE_MOVIE 0 /* set to 1 to produce movies for wave height and energy simultaneously */
#define SAVE_MEMORY 1 /* set to 1 to save memory when writing tiff images */
#define NO_EXTRA_BUFFER_SWAP 1 /* some OS require one less buffer swap when recording images */
/* General geometrical parameters */
#define WINWIDTH 1920 /* window width */
// #define WINWIDTH 1920 /* window width */
#define WINWIDTH 1150 /* window width */
#define WINHEIGHT 1150 /* window height */
#define NX 1600 /* number of grid points on x axis */
#define NY 960 /* number of grid points on y axis */
@@ -54,8 +55,10 @@
#define XMIN -2.0
#define XMAX 2.0 /* x interval */
#define YMIN -1.041666667
#define YMAX 1.041666667 /* y interval for 9/16 aspect ratio */
#define YMIN -2.0
#define YMAX 2.0 /* y interval for 9/16 aspect ratio */
// #define YMIN -1.041666667
// #define YMAX 1.041666667 /* y interval for 9/16 aspect ratio */
/* Choice of simulated equation */
@@ -102,7 +105,7 @@
#define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */
#define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */
#define PDISC_FACTOR 3.25 /* controls density of Poisson disc process (default: 3.25) */
#define PDISC_FACTOR 4.0 /* controls density of Poisson disc process (default: 3.25) */
#define RANDOM_POLY_ANGLE 0 /* 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 */
@@ -172,8 +175,8 @@
#define KS_CHIMAX_FACTOR 30.0 /* limiter on the value of chi/D_KSU */
/* Gray-Scott model */
#define A_GS 0.029 /* coupling constant in Gray-Scott model */
#define B_GS 0.06 /* coupling constant in Gray-Scott model */
#define A_GS 0.014 /* coupling constant in Gray-Scott model */
#define B_GS 0.054 /* coupling constant in Gray-Scott model */
#define D_GSU 1.0 /* u diffusion in Gray-Scott model */
#define D_GSV 2.0 /* v diffusion in Gray-Scott model */
#define GS_RSCALE 0.1 /* scaling factor for reaction term of Gray-Scott model */
@@ -226,16 +229,16 @@
#define VISCOSITY_MAX 2.0 /* max value of viscosity beyond which NVID is increased and integration step is decrase, for numerical stability */
#define CHANGE_RPSLZB 0 /* set to 1 to change second parameter in Rock-Paper-Scissors-Lizard-Spock equation */
#define RPSLZB_CHANGE 0.75 /* factor by which to rpslzb parameter */
#define RPSLZB_CHANGE 0.75 /* amount by which to rpslzb parameter */
#define RPSLZB_INITIAL_TIME 0 /* initial time during which rpslzb remains constant */
#define RPSLZB_FINAL_TIME 500 /* final time during which rpslzb remains constant */
#define CHANGE_A_GS 0 /* set to 1 to change first parameter in Gray-Scott model */
#define CHANGE_A_GS 1 /* set to 1 to change first parameter in Gray-Scott model */
#define A_GS_CHANGE -0.02 /* amount by which to decrease a_gs parameter */
#define A_GS_INITIAL_TIME 0 /* initial time during which a_gs remains constant */
#define A_GS_FINAL_TIME 100 /* final time during which a_gs remains constant */
#define CHANGE_B_GS 0 /* set to 1 to change second parameter in Gray-Scott model */
#define CHANGE_B_GS 1 /* set to 1 to change second parameter in Gray-Scott model */
#define B_GS_CHANGE -0.01 /* amount by which to decrease b_gs parameter */
#define B_GS_INITIAL_TIME 0 /* initial time during which b_gs remains constant */
#define B_GS_FINAL_TIME 100 /* final time during which b_gs remains constant */
@@ -273,8 +276,8 @@
/* Parameters for length and speed of simulation */
#define NSTEPS 1400 /* number of frames of movie */
#define NVID 30 /* number of iterations between images displayed on screen */
#define NSTEPS 1800 /* number of frames of movie */
#define NVID 20 /* number of iterations between images displayed on screen */
#define ACCELERATION_FACTOR 1.0 /* factor by which to increase NVID in course of simulation */
#define DT_ACCELERATION_FACTOR 1.0 /* factor by which to increase time step in course of simulation */
#define MAX_DT 0.024 /* maximal value of integration step */
@@ -312,12 +315,12 @@
/* Plot type - color scheme */
#define CPLOT 0
#define CPLOT_B 80
#define CPLOT_B 1
/* Plot type - height of 3D plot */
#define ZPLOT 80 /* z coordinate in 3D plot */
#define ZPLOT_B 80 /* z coordinate in second 3D plot */
#define ZPLOT 0 /* z coordinate in 3D plot */
#define ZPLOT_B 1 /* z coordinate in second 3D plot */
#define AMPLITUDE_HIGH_RES 1 /* set to 1 to increase resolution of P_3D_AMPLITUDE plot */
#define NON_DIRICHLET_BC 0 /* set to 1 to draw only facets in domain, if field is not zero on boundary */
@@ -365,8 +368,8 @@
/* Color schemes, see list in global_pdes.c */
#define COLOR_PALETTE 17 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 14 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE 11 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 11 /* Color palette, see list in global_pdes.c */
#define BLACK 1 /* black background */
#define COLOR_OUT_R 1.0 /* color outside domain */
@@ -380,7 +383,7 @@
#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 VSHIFT_AMPLITUDE -0.4 /* additional shift for wave amplitude */
#define VSHIFT_AMPLITUDE -0.25 /* additional shift for wave amplitude */
#define VSCALE_AMPLITUDE 2.5 /* additional scaling factor for color scheme P_3D_AMPLITUDE */
#define VSCALE_TIMESLICE 1.6 /* additional scaling factor for color scheme Z_TIMESLICE */
#define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */
@@ -438,6 +441,7 @@
#define IOR_TOTAL_TURNS 1.5 /* total angle of rotation for IOR_PERIODIC_WELLS_ROTATING */
#define MANDEL_IOR_SCALE -0.05 /* parameter controlling dependence of IoR on Mandelbrot escape speed */
#define OMEGA 0.005 /* frequency of periodic excitation */
#define OMEGA_B 0.005 /* frequency of periodic excitation */
#define OSCIL_YMAX 0.2 /* defines oscillation range */
#define OSCIL_YMID -0.75 /* defines oscilling beam midpoint */
#define COURANT 0.08 /* Courant number */
@@ -491,7 +495,7 @@ double u_3d[2] = {0.75, -0.45}; /* projections of basis vectors for REP_AXO_
double v_3d[2] = {-0.75, -0.45};
double w_3d[2] = {0.0, 0.015};
double light[3] = {0.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 observer[3] = {-3.0, -3.0, 4.0}; /* location of observer for REP_PROJ_3D representation */
int reset_view = 0; /* switch to reset 3D view parameters (for option ROTATE_VIEW) */
/* constants for simulations on planets */
@@ -506,8 +510,8 @@ int reset_view = 0; /* switch to reset 3D view parameters (for option RO
#define PLANET_SEALEVEL 0.0 /* sea level for flooded planet */
#define VENUS_NODATA_FACTOR 0.5 /* altitude to assign to DEM points without data (fraction of mean altitude) */
#define Z_SCALING_FACTOR 0.75 /* overall scaling factor of z axis for REP_PROJ_3D representation */
#define XY_SCALING_FACTOR 1.7 /* overall scaling factor for on-screen (x,y) coordinates after projection */
#define Z_SCALING_FACTOR 1.0 /* overall scaling factor of z axis for REP_PROJ_3D representation */
#define XY_SCALING_FACTOR 3.0 /* overall scaling factor for on-screen (x,y) coordinates after projection */
#define ZMAX_FACTOR 1.0 /* max value of z coordinate for REP_PROJ_3D representation */
#define XSHIFT_3D 0.0 /* overall x shift for REP_PROJ_3D representation */
#define YSHIFT_3D 0.0 /* overall y shift for REP_PROJ_3D representation */
@@ -516,7 +520,7 @@ int reset_view = 0; /* switch to reset 3D view parameters (for option RO
#define ZMAX 3.0 /* maximal value of z coordinate */
#define ZMIN -3.0 /* maximal value of z coordinate */
#define RSCALE 0.25 /* scaling factor of radial component */
#define RSCALE 0.3 /* scaling factor of radial component */
#define RSCALE_B 1.00 /* experimental, additional radial scaling factor in ij_to_sphere */
#define RSHIFT 0.0 /* shift in radial component */
#define RMAX 4.0 /* max value of radial component */
@@ -528,6 +532,7 @@ int reset_view = 0; /* switch to reset 3D view parameters (for option RO
#define VMAX 10.0 /* max value of wave amplitude */
#define TEST_GRADIENT 0 /* print norm squared of gradient */
#define REFRESH_B (ZPLOT_B != ZPLOT)||(CPLOT_B != CPLOT) /* to save computing time, to be improved */
#define COMPUTE_WRAP_ANGLE ((WRAP_ANGLE)&&((cplot == Z_ANGLE_GRADIENT)||(cplot == Z_ANGLE_GRADIENTX)||(cplot == Z_ARGUMENT)||(cplot == Z_ANGLE_GRADIENTX)||(cplot == Z_EULER_DIRECTION_SPEED)||(cplot == Z_SWATER_DIRECTION_SPEED)||(cplot == Z_ANGLE)))
#define PRINT_PARAMETERS ((PRINT_TIME)||(PRINT_VISCOSITY)||(PRINT_RPSLZB)||(PRINT_A_GS)||(PRINT_B_GS)||(PRINT_PROBABILITIES)||(PRINT_NOISE)||(PRINT_FLOW_SPEED)||(PRINT_AVERAGE_SPEED))
schrodinger.c
+1
View File
@@ -170,6 +170,7 @@
#define ACHIRP 0.2 /* acceleration coefficient in chirp */
#define DAMPING 0.0 /* damping of periodic excitation */
#define OMEGA 0.001 /* frequency of periodic excitation */
#define OMEGA_B 0.001 /* frequency of periodic excitation */
#define AMPLITUDE 0.8 /* amplitude of periodic excitation */
/* end of constants only used by wave_billiard and wave_3d */
sub_lj.c
+601 -18
View File
@@ -409,6 +409,42 @@ void draw_colored_polygon(double x, double y, double r, int nsides, double angle
glEnd();
}
void draw_scaled_polygon(double x, double y, double rx, double ry, int nsides, double angle)
{
int i;
double pos[2], alpha, dalpha;
dalpha = DPI/(double)nsides;
glColor3f(1.0, 1.0, 1.0);
glBegin(GL_LINE_LOOP);
for (i=0; i<=nsides; i++)
{
alpha = angle + (double)i*dalpha;
glVertex2d(x + rx*cos(alpha), y + ry*sin(alpha));
}
glEnd();
}
void draw_scaled_colored_polygon(double x, double y, double rx, double ry, int nsides, double angle, double rgb[3])
{
int i;
double pos[2], alpha, dalpha;
dalpha = DPI/(double)nsides;
glColor3f(rgb[0], rgb[1], rgb[2]);
glBegin(GL_TRIANGLE_FAN);
glVertex2d(x, y);
for (i=0; i<=nsides; i++)
{
alpha = angle + (double)i*dalpha;
glVertex2d(x + rx*cos(alpha), y + ry*sin(alpha));
}
glEnd();
}
void draw_rhombus(double x, double y, double r, double angle)
{
int i;
@@ -2856,6 +2892,302 @@ int init_maze_segments(t_segment segment[NMAXSEGMENTS], int diag)
free(maze);
}
int init_maze_segments_sphere_reg(t_segment segment[NMAXSEGMENTS])
/* init segments forming a maze */
{
t_maze* maze;
int i, j, n;
double x1, y1, x2, y2, dx, dy, padding = 0.5, width = MAZE_WIDTH, phishift;
maze = (t_maze *)malloc(NXMAZE*NYMAZE*sizeof(t_maze));
init_cylinder_maze(maze);
/* build walls of maze */
dx = DPI/(double)(NXMAZE);
dy = (YMAX - YMIN - 2.0*padding)/(double)(NYMAZE);
phishift = 0.5*dx;
for (i=0; i<NXMAZE; i++)
for (j=0; j<NYMAZE; j++)
{
n = nmaze(i, j);
x1 = phishift + (double)i*dx;
y1 = YMIN + padding + (double)j*dy;
if (maze[n].west)
{
x2 = x1 - width;
y2 = y1 + dy;
add_rectangle_to_segments(x1, y1, x1 - width, y1 + dy, segment, 0);
if (x2 > DPI) add_rectangle_to_segments(x1 - DPI, y1, x1 - DPI - width, y1 + dy, segment, 0);
}
if (((j>0)||(i!=NXMAZE*3/4))&&(maze[n].south))
{
x2 = x1 + dx;
y2 = y1 - width;
add_rectangle_to_segments(x1, y1, x1 + dx, y1 - width, segment, 0);
if (x2 > DPI) add_rectangle_to_segments(x1 - DPI, y1, x1 - DPI + dx, y1 - width, segment, 0);
}
}
/* top side of maze */
for (i=0; i<NXMAZE; i++)
{
x1 = phishift + (double)i*dx;
y1 = YMAX - padding;
if (i!=NXMAZE/4)
{
add_rectangle_to_segments(x1, y1, x1 + dx, y1 - width, segment, 0);
if (x1 + dx > DPI) add_rectangle_to_segments(x1 - DPI, y1, x1 - DPI + dx, y1 - width, segment, 0);
}
}
free(maze);
}
int init_maze_segments_sphere(t_segment segment[NMAXSEGMENTS])
/* init segments forming a maze */
{
t_maze* maze;
int nblocks, block, i, j, k, n, p, q, np, na, inrect;
double x1, y1, x2, y2, dx, dy, padding = 0.02;
double rmin, rmax, angle, r, dr, phi, dphi, phi1, ww, width = 0.03, theta, theta1, phaseshift, phishift;
maze = (t_maze *)malloc(NXMAZE*NYMAZE*sizeof(t_maze));
init_circular_maze(maze);
/* build walls of maze */
nblocks = NYMAZE/NXMAZE;
rmin = 0.4;
rmax = 1.5*PID + 0.2;
angle = DPI/(double)nblocks;
dr = (rmax - rmin)/(double)(NXMAZE);
phishift = 0.1;
/* add north-south walls */
for (block = 0; block < nblocks; block++)
{
printf("Block %i\n", block);
dphi = angle;
/* first circle */
n = nmaze(0, block*NXMAZE);
r = rmin - 0.5*width;
phi = phishift + (double)block*angle;
if (maze[n].south)
{
x1 = phi - 0.5*width;
y1 = YMAX - r;
x2 = phi + 0.5*width;
y2 = YMAX - r - dr - width;
if (x1 > XMAX) x1 += XMIN - XMAX;
if (x2 > XMAX) x2 += XMIN - XMAX;
add_rectangle_to_segments(x1, y1, x2, y2, segment, 0);
np++;
}
/* second circle */
r = rmin + dr - 0.5*width;
dphi *= 0.5;
for (q=0; q<2; q++)
{
n = nmaze(1, block*NXMAZE + q);
phi = phishift + (double)(block)*angle + (double)q*dphi;
if (maze[n].south)
{
x1 = phi - 0.5*width;
y1 = YMAX - r;
x2 = phi + 0.5*width;
y2 = YMAX - r - dr - width;
if (x1 > XMAX) x1 += XMIN - XMAX;
if (x2 > XMAX) x2 += XMIN - XMAX;
add_rectangle_to_segments(x1, y1, x2, y2, segment, 0);
np++;
}
}
/* other circles */
ww = 2;
i = 2;
while (ww < NXMAZE)
{
dphi *= 0.5;
for (p = 0; p < ww; p++)
{
r = rmin + (double)i*dr - 0.5*width;
for (q = 0; q < 2*ww; q++)
{
j = block*NXMAZE + q;
n = nmaze(i,j);
phi = phishift + (double)(block)*angle + (double)q*dphi;
if (maze[n].south)
{
x1 = phi - 0.5*width;
y1 = YMAX - r;
x2 = phi + 0.5*width;
y2 = YMAX - r - dr - width;
if (x1 > XMAX) x1 += XMIN - XMAX;
if (x2 > XMAX) x2 += XMIN - XMAX;
add_rectangle_to_segments(x1, y1, x2, y2, segment, 0);
np++;
}
}
i++;
}
ww *= 2;
}
}
/* add east-west walls */
for (block = 0; block < nblocks; block++)
{
dphi = angle;
/* first circle */
n = nmaze(0, block*NXMAZE);
r = rmin;
phi = phishift + (double)block*angle;
if ((block > 0)&&(maze[n].west))
{
// x1 = phi;
// y1 = YMAX - (r - 0.5*width);
// x2 = phi + dphi + 0.05*width;
// y2 = YMAX - (r + 0.5*width);
//
// add_rectangle_to_segments(x1, y1, x2, y2, segment, 0);
//
// if (x2 > XMAX) add_rectangle_to_segments(x1 - DPI, y1, x2 - DPI, y2, segment, 0);
// if (x1 < XMIN) add_rectangle_to_segments(x1 + DPI, y1, x2 + DPI, y2, segment, 0);
// np++;
y1 = YMAX - (r - 0.5*width);
y2 = YMAX - (r + 0.5*width);
x1 = phi;
while (x1 < phi + dphi + 0.05*width)
{
x2 = x1 + 0.25*dphi;
add_rectangle_to_segments(x1, y1, x2, y2, segment, 0);
if (x2 > XMAX) add_rectangle_to_segments(x1 - DPI, y1, x2 - DPI, y2, segment, 0);
if (x1 < XMIN) add_rectangle_to_segments(x1 + DPI, y1, x2 + DPI, y2, segment, 0);
np++;
x1 += 0.25*dphi;
}
}
/* second circle */
r = rmin + dr;
dphi *= 0.5;
for (q=0; q<2; q++)
{
n = nmaze(1, block*NXMAZE + q);
phi = phishift + (double)(block)*angle + (double)q*dphi;
if (maze[n].west)
{
// x1 = phi;
// y1 = YMAX - (r - 0.5*width);
// x2 = phi + dphi + 0.05*width;
// y2 = YMAX - (r + 0.5*width);
//
// add_rectangle_to_segments(x1, y1, x2, y2, segment, 0);
//
// if (x2 > XMAX) add_rectangle_to_segments(x1 - DPI, y1, x2 - DPI, y2, segment, 0);
// if (x1 < XMIN) add_rectangle_to_segments(x1 + DPI, y1, x2 + DPI, y2, segment, 0);
//
// np++;
y1 = YMAX - (r - 0.5*width);
y2 = YMAX - (r + 0.5*width);
x1 = phi;
while (x1 < dphi + 0.05*width)
{
x2 = x1 + 0.25*dphi;
add_rectangle_to_segments(x1, y1, x2, y2, segment, 0);
if (x2 > XMAX) add_rectangle_to_segments(x1 - DPI, y1, x2 - DPI, y2, segment, 0);
if (x1 < XMIN) add_rectangle_to_segments(x1 + DPI, y1, x2 + DPI, y2, segment, 0);
np++;
x1 += 0.25*dphi;
}
}
}
/* other circles */
ww = 2;
i = 2;
while (ww < NXMAZE)
{
dphi *= 0.5;
for (p = 0; p < ww; p++)
{
r = rmin + (double)i*dr;
for (q = 0; q < 2*ww; q++)
{
j = block*NXMAZE + q;
n = nmaze(i,j);
phi = phishift + (double)(block)*angle + (double)q*dphi;
if (maze[n].west)
{
x1 = phi;
y1 = YMAX - (r - 0.5*width);
x2 = phi + dphi + 0.05*width;
y2 = YMAX - (r + 0.5*width);
add_rectangle_to_segments(x1, y1, x2, y2, segment, 0);
if (x2 > XMAX) add_rectangle_to_segments(x1 - DPI, y1, x2 - DPI, y2, segment, 0);
if (x1 < XMIN) add_rectangle_to_segments(x1 + DPI, y1, x2 + DPI, y2, segment, 0);
np++;
}
}
i++;
}
ww *= 2;
}
}
/* outer boundary of maze */
y1 = YMAX - (rmax - 0.5*width);
y2 = YMAX - (rmax + 0.5*width);
x1 = phishift + dphi;
while (x1 < DPI)
{
x2 = x1 + dphi;
add_rectangle_to_segments(x1, y1, x2, y2, segment, 0);
if (x2 > XMAX) add_rectangle_to_segments(x1 - DPI, y1, x2 - DPI, y2, segment, 0);
x1 += dphi;
}
np++;
free(maze);
}
void translate_segments(t_segment segment[NMAXSEGMENTS], double deltax[2], double deltay[2])
/* rotates the repelling segments by given angle */
{
@@ -3283,6 +3615,100 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS], t_belt conveyor_belt[N
break;
}
case (S_FOUR_VERTICALS):
{
width = 2.0*MU;
y = 0.5*(YMIN + YMAX);
add_rectangle_to_segments(width, y, -width, YMAX + 0.1, segment, 0);
add_rectangle_to_segments(XMAX+width, y, XMAX-width, YMAX + 0.1, segment, 0);
x = 0.5*(XMIN + XMAX);
add_rectangle_to_segments(x+width, y, x-width, YMAX + 0.1, segment, 0);
x = 0.25*(3.0*XMIN + XMAX);
add_rectangle_to_segments(x+width, YMIN - 0.1, x-width, y, segment, 0);
x = 0.25*(XMIN + 3.0*XMAX);
add_rectangle_to_segments(x+width, YMIN - 0.1, x-width, y, segment, 0);
cycle = 0;
concave = 0; /* add_rectangle_to_segments already deals with concave corners */
for (i=0; i<nsegments; i++)
{
// segment[i].concave = 1;
segment[i].group = 0;
segment[i].inactivate = 0;
}
break;
}
case (S_TWO_VERTICALS):
{
width = 2.0*MU;
y = 0.5*(YMIN + YMAX);
x = 0.25*(3.0*XMIN + XMAX);
add_rectangle_to_segments(x+width, YMIN - 0.1, x-width, y, segment, 0);
x = 0.25*(XMIN + 3.0*XMAX);
add_rectangle_to_segments(x+width, YMIN - 0.1, x-width, y, segment, 0);
cycle = 0;
concave = 0; /* add_rectangle_to_segments already deals with concave corners */
for (i=0; i<nsegments; i++)
{
// segment[i].concave = 1;
segment[i].group = 0;
segment[i].inactivate = 0;
}
break;
}
case (S_TWO_DIAGONALS):
{
width = 2.0*MU;
y = 0.5*(YMIN + YMAX);
x = 0.25*(3.0*XMIN + XMAX);
add_rotated_angle_to_segments(x, y, x - LAMBDA, y + LAMBDA, width, 1, segment, 0);
x = 0.25*(XMIN + 3.0*XMAX);
add_rotated_angle_to_segments(x - LAMBDA, y - LAMBDA, x, y, width, 1, segment, 0);
cycle = 0;
concave = 0; /* add_rectangle_to_segments already deals with concave corners */
for (i=0; i<nsegments; i++)
{
// segment[i].concave = 1;
segment[i].group = 0;
segment[i].inactivate = 0;
}
break;
}
case (S_TWO_DIAGONALS_ROT):
{
width = 2.0*MU;
y = 0.5*(YMIN + YMAX);
x = 0.25*(3.0*XMIN + XMAX);
add_rotated_angle_to_segments(x, y, x + LAMBDA, y + LAMBDA, width, 1, segment, 0);
x = 0.25*(XMIN + 3.0*XMAX);
add_rotated_angle_to_segments(x + LAMBDA, y - LAMBDA, x, y, width, 1, segment, 0);
cycle = 0;
concave = 0; /* add_rectangle_to_segments already deals with concave corners */
for (i=0; i<nsegments; i++)
{
// segment[i].concave = 1;
segment[i].group = 0;
segment[i].inactivate = 0;
}
break;
}
case (S_CUP):
{
angle = APOLY*PID;
@@ -3778,6 +4204,32 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS], t_belt conveyor_belt[N
break;
}
case (S_MAZE_SPHERE):
{
init_maze_segments_sphere(segment);
cycle = 0;
for (i=0; i<nsegments; i++)
{
segment[i].group = 0;
segment[i].inactivate = 0;
}
break;
}
case (S_MAZE_SPHERE_REG):
{
init_maze_segments_sphere_reg(segment);
cycle = 0;
for (i=0; i<nsegments; i++)
{
segment[i].group = 0;
segment[i].inactivate = 0;
}
break;
}
case (S_EXT_RECTANGLE):
{
width = 0.1*LAMBDA;
@@ -4026,7 +4478,7 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS], t_belt conveyor_belt[N
}
case (S_POLYGON_EXT):
{
add_circle_to_segments(0.0, 0.0, LAMBDA, NPOLY, APOLY, segment, 0);
add_circle_to_segments(0.5*(XMIN+XMAX), 0.5*(YMIN+YMAX), LAMBDA, NPOLY, APOLY, segment, 0);
cycle = 0;
concave = 0;
@@ -4041,6 +4493,42 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS], t_belt conveyor_belt[N
break;
}
case (S_TWO_POLYGONS_EXT):
{
add_circle_to_segments(XMIN + 0.25*(XMAX-XMIN), 0.5*(YMIN+YMAX), LAMBDA, NPOLY, APOLY, segment, 0);
add_circle_to_segments(XMIN + 0.75*(XMAX-XMIN), 0.5*(YMIN+YMAX), LAMBDA, NPOLY, APOLY, segment, 0);
cycle = 0;
concave = 0;
ngroups = 1;
for (i=0; i<nsegments; i++)
{
segment[i].concave = 1;
segment[i].inactivate = 0;
}
break;
}
case (S_TWO_POLYGONS_EXT_SHIFTED):
{
y1 = 0.4*YMIN + 0.6*YMAX;
add_circle_to_segments(XMIN + 0.25*(XMAX-XMIN), y1, LAMBDA, NPOLY, APOLY, segment, 0);
y1 = 0.6*YMIN + 0.4*YMAX;
add_circle_to_segments(XMIN + 0.75*(XMAX-XMIN), y1, LAMBDA, NPOLY, APOLY, segment, 0);
cycle = 0;
concave = 0;
ngroups = 1;
for (i=0; i<nsegments; i++)
{
segment[i].concave = 1;
segment[i].inactivate = 0;
}
break;
}
case (S_WEDGE_EXT):
{
angle = AWEDGE*PID;
@@ -4675,6 +5163,19 @@ int in_segment_region(double x, double y, t_segment segment[NMAXSEGMENTS])
if (y <= BCYMIN) return(0);
switch (SEGMENT_PATTERN) {
case (S_TWO_VERTICALS):
{
y1 = 0.5*(YMIN + YMAX);
width = 4.0*MU;
if (y < y1 + MU)
{
x1 = 0.25*(3.0*XMIN + XMAX);
if (vabs(x - x1) < width) return(0);
x1 = 0.25*(XMIN + 3.0*XMAX);
if (vabs(x - x1) < width) return(0);
}
return(1);
}
case (S_CUP):
{
angle = APOLY*PID;
@@ -4682,6 +5183,7 @@ int in_segment_region(double x, double y, t_segment segment[NMAXSEGMENTS])
if (y < BCYMAX - (BCYMAX - BCYMIN)*(x - BCXMIN)/dx) return(0);
if (y < BCYMAX - (BCYMAX - BCYMIN)*(BCXMAX - x)/dx) return(0);
break;
}
case (S_HOURGLASS):
{
@@ -4823,7 +5325,23 @@ int in_segment_region(double x, double y, t_segment segment[NMAXSEGMENTS])
case (S_POLYGON_EXT):
{
padding = 3.0*MU;
if (in_polygon(x, y, LAMBDA + padding, NPOLY, APOLY)) return(0);
if (in_polygon(x - 0.5*(XMIN+XMAX), y - 0.5*(YMIN+YMAX), LAMBDA + padding, NPOLY, APOLY)) return(0);
else return(1);
}
case (S_TWO_POLYGONS_EXT):
{
padding = 3.0*MU;
if (in_polygon(x - 0.25*XMIN - 0.75*XMAX, y - 0.5*(YMIN+YMAX), LAMBDA + padding, NPOLY, APOLY)) return(0);
else if (in_polygon(x - 0.75*XMIN - 0.25*XMAX, y - 0.5*(YMIN+YMAX), LAMBDA + padding, NPOLY, APOLY)) return(0);
else return(1);
}
case (S_TWO_POLYGONS_EXT_SHIFTED):
{
padding = 3.0*MU;
y1 = 0.6*YMIN + 0.4*YMAX;
y2 = 0.4*YMIN + 0.6*YMAX;
if (in_polygon(x - 0.25*XMIN - 0.75*XMAX, y - y1, LAMBDA + padding, NPOLY, APOLY)) return(0);
else if (in_polygon(x - 0.75*XMIN - 0.25*XMAX, y - y2, LAMBDA + padding, NPOLY, APOLY)) return(0);
else return(1);
}
case (S_WEDGE_EXT):
@@ -8186,8 +8704,13 @@ void draw_one_particle(t_particle particle, double xc, double yc, double radius,
draw_colored_ellipse_precomp(xc1, yc1, radius/(sin(yc1) + SIN_THETA_REG), radius, rgb);
else draw_colored_circle_precomp(xc1, yc1, radius, rgb);
}
else
{
if ((SPHERE)&&(DRAW_ELLIPSES_ON_SPHERE))
draw_scaled_colored_polygon(xc1, yc1, radius/(sin(yc1) + SIN_THETA_REG), radius, nsides, angle + APOLY*PID, rgb);
else draw_colored_polygon(xc1, yc1, radius, nsides, angle + APOLY*PID, rgb);
}
}
/* draw crosses/bars on charged particles */
// if (TWO_TYPES)
@@ -8254,8 +8777,13 @@ void draw_one_particle(t_particle particle, double xc, double yc, double radius,
draw_ellipse_precomp(xc1, yc1, radius/(sin(yc1) + SIN_THETA_REG), radius);
else draw_circle_precomp(xc1, yc1, radius);
}
else
{
if ((SPHERE)&&(DRAW_ELLIPSES_ON_SPHERE))
draw_scaled_polygon(xc1, yc1, radius/(sin(yc1) + SIN_THETA_REG), radius, nsides, angle + APOLY*PID);
else draw_polygon(xc1, yc1, radius, nsides, angle + APOLY*PID);
}
}
if (particle.interaction == I_LJ_WATER) for (wsign = -1; wsign <= 1; wsign+=2)
{
@@ -8888,7 +9416,53 @@ void draw_potential(t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE])
value = atan(BG_POTENTIAL_SLOPE*pot)/PI + 0.5;
// value = 0.5*(tanh(BG_POTENTIAL_SLOPE*pot)/PI + 1.0);
hue = ENERGY_HUE_MIN + (ENERGY_HUE_MAX - ENERGY_HUE_MIN)*value;
hsl_to_rgb_turbo(hue, 0.9, 0.5, rgb);
// hsl_to_rgb_turbo(hue, 0.9, 0.5, rgb);
hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_POTENTIAL);
if (SHADE_BG_COLOR_2D)
for (k=0; k<3; k++) rgb[k] *= wsphere[cell].cos_angle_pot;
phi = XMIN + (double)i*dphi;
theta = YMIN + (double)j*dtheta;
glColor3f(rgb[0], rgb[1], rgb[2]);
glVertex2d(phi, theta);
glVertex2d(phi + dphi_res, theta);
glVertex2d(phi + dphi_res, theta + dtheta_res);
glVertex2d(phi, theta + dtheta_res);
}
glEnd();
}
void draw_dem(t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE])
/* color background according to potential on sphere */
{
int i, j, k, cell, res = 3;
double pot, value, phi, theta, hue, rgb[3];
static double dphi, dtheta, dphi_res, dtheta_res;
static int first = 1;
if (first)
{
dphi = (XMAX-XMIN)/(double)NX_SPHERE;
dtheta = (YMAX-YMIN)/(double)NY_SPHERE;
dphi_res = (double)res*dphi;
dtheta_res = (double)res*dtheta;
first = 0;
}
glBegin(GL_QUADS);
for (i=0; i<NX_SPHERE; i+=res)
for (j=0; j<NY_SPHERE; j+=res)
{
cell = i*NY_SPHERE + j;
rgb[0] = wsphere[cell].r0;
rgb[1] = wsphere[cell].g0;
rgb[2] = wsphere[cell].b0;
// printf("rgb = (%.3f, %.3f, %.3f\n", rgb[0], rgb[1], rgb[2]);
if (SHADE_BG_COLOR_2D)
for (k=0; k<3; k++) rgb[k] *= wsphere[cell].cos_angle_pot;
@@ -9928,7 +10502,7 @@ void print_parameters(t_lj_parameters params, short int left, double pressure[N_
{
char message[100];
int i, j, k;
double density, hue, rgb[3], logratio, x, y, meanpress[N_PRESSURES], phi, sphi, dphi, pprint, mean_temp, lengthcontainer, boundary_force, fx, fy, r1, r2;
double density, hue, rgb[3], logratio, x, y, meanpress[N_PRESSURES], phi, sphi, dphi, pprint, mean_temp, lengthcontainer, boundary_force, fx, fy, r1, r2, rbox;
static double xbox, xtext, xmid, xmidtext, xxbox, xxtext, pressures[N_P_AVERAGE], meanpressure = 0.0, maxpressure = 0.0, mean_fx, mean_fy;
static double press[N_PRESSURES][N_P_AVERAGE], temp[N_T_AVERAGE], scale;
static int first = 1, i_pressure, i_temp;
@@ -9939,19 +10513,20 @@ void print_parameters(t_lj_parameters params, short int left, double pressure[N_
if (DRAW_SPHERE)
{
scale = 2.0/2.5;
rbox = 2.0*(double)WINWIDTH/1760.0;
if (left)
{
xbox = -2.0 + 0.45*scale;
xtext = -2.0 + 0.12*scale;
xxbox = 2.0 - 0.39*scale;
xxtext = 2.0 - 0.73*scale;
xbox = -rbox + 0.45*scale;
xtext = -rbox + 0.12*scale;
xxbox = rbox - 0.39*scale;
xxtext = rbox - 0.73*scale;
}
else
{
xbox = 2.0 - 0.45*scale;
xtext = 2.0 - 0.73*scale;
xxbox = -2.0 + 0.4*scale;
xxtext = -2.0 + 0.08*scale;
xbox = rbox - 0.45*scale;
xtext = rbox - 0.73*scale;
xxbox = -rbox + 0.4*scale;
xxtext = -rbox + 0.08*scale;
}
xmid = xbox;
xmidtext = xtext;
@@ -10215,7 +10790,7 @@ void print_parameters(t_lj_parameters params, short int left, double pressure[N_
y -= 0.1;
}
if (SPHERE_GRAVITY > 0)
if ((SPHERE_GRAVITY > 0)&&(GRAVITY_DELTA_ANGLE != 0.0))
{
erase_area_hsl(xmid, y + 0.025*scale, 0.22*scale, 0.05*scale, 0.0, 0.9, 0.0);
glColor3f(1.0, 1.0, 1.0);
@@ -11992,7 +12567,7 @@ t_segment segment[NMAXSEGMENTS], t_molecule molecule[NMAXCIRCLES])
{
while (!particle[i].active) i++;
tracer_n[n_tracers] = i;
printf("%i th tracer particle %i\n", n_tracers, i);
// printf("%i th tracer particle %i\n", n_tracers, i);
i++;
n_tracers++;
}
@@ -12133,7 +12708,6 @@ t_segment segment[NMAXSEGMENTS], t_molecule molecule[NMAXCIRCLES])
particle[i].active = 0;
}
if (ADD_FIXED_OBSTACLES)
{
for (i=0; i< ncircles; i++) for (j=0; j < nobstacles; j++)
@@ -18059,9 +18633,14 @@ void draw_frame_2d(int i, int plot, int bg_color, int ncollisions, int traj_posi
compute_all_particle_colors(particle, cluster, plot);
if ((COLOR_BACKGROUND)&&(bg_color > 0))
{
compute_background_color(particle, segment, obstacle, bg_color, hashgrid, wsphere);
// compute_background_color(particle, segment, obstacle, bg_color, hashgrid, wsphere);
if (bg_color == BG_POTENTIAL) draw_potential(wsphere);
else draw_background(particle, segment, obstacle, bg_color, hashgrid);
else if (bg_color == BG_DEM) draw_dem(wsphere);
else
{
compute_background_color(particle, segment, obstacle, bg_color, hashgrid, wsphere);
draw_background(particle, segment, obstacle, bg_color, hashgrid);
}
}
// else if (!TRACER_PARTICLE) blank();
@@ -18114,16 +18693,20 @@ void draw_frame_3d(int i, int plot, int bg_color, int ncollisions, int traj_posi
// }
/* initialize background colors */
if ((COLOR_BACKGROUND)&&(bg_color > 0))
if ((COLOR_BACKGROUND)&&(bg_color > 0)&&(bg_color != BG_DEM))
compute_background_color(particle, segment, obstacle, bg_color, hashgrid, wsphere);
// printf("1\n");
/* compute all particle colors */
compute_all_particle_colors(particle, cluster, plot);
// printf("2\n");
/* initialize sphere radius and color */
init_sphere_radius(particle, hashgrid, cluster, trajectory, segment, traj_position, traj_length, wsphere, tracer_n, plot, bg_color, absorber);
// printf("3\n");
/* initialize light angle (cf compute_light_angle_sphere_rde() in sub_rde.c) */
compute_light_angle_sphere(wsphere);
/* draw sphere with particles (cf draw_wave_sphere_3d_rde() in sub_rde.c) */
// printf("4\n");
draw_lj_sphere_3d(wsphere);
if (PRINT_PARAMETERS) print_parameters(params, PRINT_LEFT, pressure, refresh);
sub_lj_sphere.c
+349 -7
View File
@@ -1070,7 +1070,8 @@ void compute_background_color(t_particle particle[NMAXCIRCLES], t_segment segmen
value = atan(BG_POTENTIAL_SLOPE*wsphere[cell].potential)/PI + 0.5;
hue = ENERGY_HUE_MIN + (ENERGY_HUE_MAX - ENERGY_HUE_MIN)*value;
hsl_to_rgb_turbo(hue, 0.9, 0.5, rgb);
// hsl_to_rgb_turbo(hue, 0.9, 0.5, rgb);
hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_POTENTIAL);
if (SHADE_BG_COLOR_2D)
for (k=0; k<3; k++) rgb[k] *= wsphere[cell].cos_angle_pot;
@@ -1316,7 +1317,7 @@ double test_distance(double phi, double psi, int i, t_particle* particle)
/* determine particle shape */
{
double d, ca, sa, theta;
static double c1, c2, c1inv;
static double c1, c2, c1inv, clinvhalf;
static int first = 1;
if (first)
@@ -1335,6 +1336,22 @@ double test_distance(double phi, double psi, int i, t_particle* particle)
c2 = 1.0/cos(PI/5.0);
break;
}
case (I_POLYGON):
{
c1 = (double)NPOLY/DPI;
c1inv = DPI/(double)NPOLY;
c2 = 1.0/cos(PI/(double)NPOLY);
clinvhalf = PI/(double)NPOLY;
break;
}
case (I_POLYGON_ALIGN):
{
c1 = (double)NPOLY/DPI;
c1inv = DPI/(double)NPOLY;
c2 = 1.0/cos(PI/(double)NPOLY);
clinvhalf = PI/(double)NPOLY;
break;
}
}
first = 0;
}
@@ -1364,6 +1381,24 @@ double test_distance(double phi, double psi, int i, t_particle* particle)
theta = (argument(ca, sa) - particle[i].angle);
return(2.0*d*(1.0 - 0.5*cos(2.0*theta)));
}
case (I_POLYGON):
{
d = dist_point_to_particle_sphere(i, phi, psi, particle, &ca, &sa);
theta = (argument(ca, sa) - particle[i].angle + 2.0*DPI)*c1;
theta -= (double)((int)theta);
theta *= c1inv;
theta -= clinvhalf;
return(d*cos(theta)*c2);
}
case (I_POLYGON_ALIGN):
{
d = dist_point_to_particle_sphere(i, phi, psi, particle, &ca, &sa);
theta = (argument(ca, sa) - particle[i].angle + 2.0*DPI)*c1;
theta -= (double)((int)theta);
theta *= c1inv;
theta -= clinvhalf;
return(d*cos(theta)*c2);
}
default:
{
if (DRAW_SPIN)
@@ -1564,6 +1599,23 @@ void add_particle_to_sphere(int i, int j, int i0, int j0, int part, t_particle p
}
}
}
else if (POLYGON_INTERACTION)
{
dist = test_distance(x, y, part, particle);
if (dist < 1.2*r)
{
h = wsphere[i0*NY_SPHERE+j0].initial_radius;
wsphere[i*NY_SPHERE+j].radius = h;
if (dist < r) wsphere[i*NY_SPHERE+j].radius += r;
if (!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];
}
}
}
if (draw_normal_particle)
{
@@ -1585,7 +1637,7 @@ void add_particle_to_sphere(int i, int j, int i0, int j0, int part, t_particle p
}
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)
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, int bg_color)
/* draw tracer particle trajectory */
{
int i, j, i0, j0, time, p, q, width, imin, cell, i1, j1;
@@ -1625,6 +1677,15 @@ void draw_trajectory_sphere(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTI
lum1 = 1.0 - lum;
if (COLOR_BACKGROUND)
{
if (BG_COLOR == BG_DEM)
{
cell = i0*NY_SPHERE+j0;
rgb_bg[0] = wsphere[cell].r;
rgb_bg[1] = wsphere[cell].g;
rgb_bg[2] = wsphere[cell].b;
}
else
{
cell = hash_cell(x1, y1);
if (!wsphere[cell].locked)
@@ -1634,6 +1695,7 @@ void draw_trajectory_sphere(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTI
rgb_bg[2] = hashgrid[cell].b;
}
}
}
if ((x2 != x1)||(y2 != y1)) for (p=-1; p<2; p++)
for (q=-1; q<2; q++)
@@ -1777,7 +1839,7 @@ void draw_segments_sphere(t_segment segment[NMAXSEGMENTS], t_lj_sphere wsphere[N
/* 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, ca;
double x1, y1, x2, y2, x, y, dt, length, ca, xtmp;
static double dphi, dtheta;
static int first = 1;
@@ -1795,6 +1857,7 @@ void draw_segments_sphere(t_segment segment[NMAXSEGMENTS], t_lj_sphere wsphere[N
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;
@@ -1924,8 +1987,11 @@ void compute_gradient_potential_sphere(t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE])
void init_potential_sphere(t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE], t_absorber wells[NMAX_ABSORBERS])
/* initialise the potential and its gradient, for option ADD_POTENTIAL_SPHERE */
{
int i, j;
double phi, theta, phi1, theta1, dist, pot, x, y, z, max, nablax, nablay, norm, pscal, ca, shade_scale2;
int i, j, k, ii, jj, nmaxpixels, rnx, rny, nx, ny, maxrgb, scan, hmin, hmax, rgbval;
int *rgb_values;
double phi, theta, phi1, theta1, dist, pot, x, y, z, max, nablax, nablay, norm, pscal, ca, shade_scale2, hsea, cratio, rx, ry;
double *height_values;
FILE *image_file;
printf("Initialising potential on sphere\n");
for (i=0; i<NX_SPHERE*NY_SPHERE; i++) wsphere[i].potential = 0.0;
@@ -1963,7 +2029,249 @@ void init_potential_sphere(t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE], t_absorber
}
break;
}
case (SPP_PLANET):
{
switch (PLANET_DEM)
{
case (PLANET_EARTH):
{
nmaxpixels = 4915200;
image_file = fopen("digital_elevation_model_large.ppm", "r");
hsea = 12.0; /* sea level #0c0c0c */
break;
}
case (PLANET_MARS):
{
nmaxpixels = 8388608;
image_file = fopen("marscyl2.ppm", "r");
// hsea = 6.0 + 255.0*PLANET_SEALEVEL/DEM_MAXHEIGHT;
hsea = 6.0;
break;
}
case (PLANET_MOON):
{
nmaxpixels = 2097152;
image_file = fopen("Moon_LRO_LOLA_global_LDEM_1024.ppm", "r");
hsea = 0.0;
break;
}
case (PLANET_VENUS):
{
nmaxpixels = 4096*2048;
image_file = fopen("Venus_Magellan_Topography_Global_4641m_v02_scaled2.ppm", "r");
hsea = 0.0;
break;
}
case (PLANET_MERCURY):
{
nmaxpixels = 1280*380;
image_file = fopen("Mercury_Messenger_DEM_Global_665m_1024_1_cropped.ppm", "r");
hsea = 0.0;
break;
}
}
rnx = NX_SPHERE;
rny = NY_SPHERE;
printf("Reading digital elevation model\n");
rgb_values = (int *)malloc(3*nmaxpixels*sizeof(int));
height_values = (double *)malloc(rnx*rny*sizeof(double));
// height_values_tmp = (double *)malloc(rnx*rny*sizeof(double));
scan = fscanf(image_file,"%i %i\n", &nx, &ny);
scan = fscanf(image_file,"%i\n", &maxrgb);
hmin = maxrgb;
hmax = 0;
if (nx*ny > nmaxpixels)
{
printf("DEM too large, increase nmaxpixels in init_potential_sphere()\n");
exit(0);
}
/* read rgb values */
for (j=0; j<ny; j++)
for (i=0; i<nx; i++)
for (k=0; k<3; k++)
{
scan = fscanf(image_file,"%i\n", &rgbval);
rgb_values[3*(j*nx+i)+k] = rgbval;
if (rgbval < hmin) hmin = rgbval;
if (rgbval > hmax) hmax = rgbval;
}
cratio = 1.0/(double)(hmax-hmin);
rx = (double)nx/(double)(rnx);
ry = (double)ny/(double)(rny);
/* build height table */
for (i=0; i<rnx; i++)
for (j=0; j<rny; j++)
{
ii = nx/2 - (int)(rx*(double)(rnx-1 - i));
if (ii > nx-1) ii -= nx;
if (ii < 0) ii += nx;
jj = (int)(ry*(double)(rny - j));
if (jj > ny-1) jj = ny-1;
if (jj < 0) jj = 0;
height_values[i*rny+j] = ((double)rgb_values[3*(jj*nx+ii)]-hsea)*cratio;
wsphere[i*rny+j].potential = ((double)rgb_values[3*(jj*nx+ii)]-hsea)*cratio;
// wsphere[i*rny+j].r = ((double)rgb_values[3*(jj*nx+ii)])/256.0;
// wsphere[i*rny+j].g = ((double)rgb_values[3*(jj*nx+ii)+1])/256.0;
// wsphere[i*rny+j].b = ((double)rgb_values[3*(jj*nx+ii)+2])/256.0;
}
/* smooth values in case of high resolution */
if ((rnx > nx)||(rny > ny))
{
for (i=1; i<rnx-1; i++)
for (j=1; j<rny-1; j++)
{
height_values[i*rny+j] *= 0.2;
height_values[i*rny+j] += 0.2*height_values[(i+1)*rny+j];
height_values[i*rny+j] += 0.2*height_values[(i-1)*rny+j];
height_values[i*rny+j] += 0.2*height_values[i*rny+j-1];
height_values[i*rny+j] += 0.2*height_values[i*rny+j+1];
wsphere[i*rny+j].potential *= 0.2;
wsphere[i*rny+j].potential += 0.2*wsphere[(i+1)*rny+j].potential;
wsphere[i*rny+j].potential += 0.2*wsphere[(i-1)*rny+j].potential;
wsphere[i*rny+j].potential += 0.2*wsphere[i*rny+j-1].potential;
wsphere[i*rny+j].potential += 0.2*wsphere[i*rny+j+1].potential;
}
/* i = 0 */
for (j=1; j<rny-1; j++)
{
height_values[j] *= 0.2;
height_values[j] += 0.2*height_values[rny+j];
height_values[j] += 0.2*height_values[(rnx-1)*rny+j];
height_values[j] += 0.2*height_values[j-1];
height_values[j] += 0.2*height_values[j+1];
wsphere[j].potential *= 0.2;
wsphere[j].potential += 0.2*wsphere[rny+j].potential;
wsphere[j].potential += 0.2*wsphere[(rnx-1)*rny+j].potential;
wsphere[j].potential += 0.2*wsphere[j-1].potential;
wsphere[j].potential += 0.2*wsphere[j+1].potential;
}
/* i = rny-1 */
for (j=1; j<rny-1; j++)
{
height_values[(rny-1)*rny+j] *= 0.2;
height_values[(rny-1)*rny+j] += 0.2*height_values[j];
height_values[(rny-1)*rny+j] += 0.2*height_values[(rny-2)*rny+j];
height_values[(rny-1)*rny+j] += 0.2*height_values[(rny-1)*rny+j-1];
height_values[(rny-1)*rny+j] += 0.2*height_values[(rny-1)*rny+j+1];
wsphere[(rny-1)*rny+j].potential *= 0.2;
wsphere[(rny-1)*rny+j].potential += 0.2*wsphere[j].potential;
wsphere[(rny-1)*rny+j].potential += 0.2*wsphere[(rny-2)*rny+j].potential;
wsphere[(rny-1)*rny+j].potential += 0.2*wsphere[(rny-1)*rny+j-1].potential;
wsphere[(rny-1)*rny+j].potential += 0.2*wsphere[(rny-1)*rny+j+1].potential;
}
}
printf("Closing rgb_values\n");
fclose(image_file);
if (BG_COLOR == BG_DEM)
{
switch (PLANET_DEM)
{
case (PLANET_EARTH):
{
printf("Reading Blue Marble map\n");
image_file = fopen("Earth_Map_Blue_Marble_2002_large.ppm", "r");
break;
}
case (PLANET_MARS):
{
printf("Reading Mars map\n");
nmaxpixels = 8388608;
image_file = fopen("Mars_Viking_ClrMosaic_global_925m_scaled.ppm", "r");
break;
}
case (PLANET_MOON):
{
printf("Reading Moon map\n");
nmaxpixels = 2048*1024;
image_file = fopen("Moon_photo_map.ppm", "r");
break;
}
case (PLANET_VENUS):
{
printf("Reading Venus map\n");
nmaxpixels = 1440*720;
image_file = fopen("Venus_map_NASA_JPL_Magellan-Venera-Pioneer.ppm", "r");
break;
}
case (PLANET_MERCURY):
{
printf("Reading Mercury map\n");
nmaxpixels = 2304*1152;
image_file = fopen("Mercury_color_photo.ppm", "r");
break;
}
}
scan = fscanf(image_file,"%i %i\n", &nx, &ny);
scan = fscanf(image_file,"%i\n", &maxrgb);
free(rgb_values);
rgb_values = (int *)malloc(3*nmaxpixels*sizeof(int));
if (nx*ny > nmaxpixels)
{
printf("Image too large, increase nmaxpixels in init_potential_sphere()\n");
exit(0);
}
/* read rgb values */
for (j=0; j<ny; j++)
for (i=0; i<nx; i++)
for (k=0; k<3; k++)
{
scan = fscanf(image_file,"%i\n", &rgbval);
rgb_values[3*(j*nx+i)+k] = rgbval;
}
cratio = 1.0/(double)maxrgb;
rx = (double)nx/(double)(NX_SPHERE);
ry = (double)ny/(double)(NY_SPHERE);
/* build wave table */
for (i=0; i<NX_SPHERE; i++)
for (j=0; j<NY_SPHERE; j++)
{
ii = nx/2 - (int)(rx*(double)(rnx-1 - i));
if (ii > nx-1) ii -= nx;
if (ii < 0) ii += nx;
jj = (int)(ry*(double)(NY_SPHERE - j));
if (jj > ny-1) jj = ny-1;
if (jj < 0) jj = 0;
wsphere[i*NY_SPHERE+j].r0 = (double)rgb_values[3*(jj*nx+ii)]*cratio;
wsphere[i*NY_SPHERE+j].g0 = (double)rgb_values[3*(jj*nx+ii)+1]*cratio;
wsphere[i*NY_SPHERE+j].b0 = (double)rgb_values[3*(jj*nx+ii)+2]*cratio;
}
printf("Closing rgb_values\n");
fclose(image_file);
}
printf("1\n");
free(rgb_values);
free(height_values);
break;
}
}
printf("2\n");
for (i=0; i<NX_SPHERE*NY_SPHERE; i++)
{
@@ -2014,12 +2322,21 @@ void init_sphere_radius(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA
wsphere[i].radius = wsphere[i].initial_radius;
if ((COLOR_BACKGROUND)&&(bg_color > 0))
{
if (BG_COLOR == BG_DEM)
{
wsphere[i].r = wsphere[i].r0;
wsphere[i].g = wsphere[i].g0;
wsphere[i].b = wsphere[i].b0;
}
else
{
cell = hash_cell(wsphere[i].phi, wsphere[i].theta);
wsphere[i].r = hashgrid[cell].r;
wsphere[i].g = hashgrid[cell].g;
wsphere[i].b = hashgrid[cell].b;
}
}
else
{
wsphere[i].r = 1.0;
@@ -2030,7 +2347,7 @@ void init_sphere_radius(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA
/* add tracer trajectories */
if (TRACER_PARTICLE)
draw_trajectory_sphere(trajectory, hashgrid, 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, bg_color);
// if (ADD_FIXED_SEGMENTS)
// draw_segments_sphere(segment, wsphere);
@@ -2284,6 +2601,25 @@ void compute_light_angle_sphere(t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE])
}
}
void draw_axis(int north)
{
double xyz[3], z;
z = (double)north;
glColor3f(1.0, 1.0, 1.0);
glBegin(GL_LINE_STRIP);
xyz[0] = 0.0;
xyz[1] = 0.0;
xyz[2] = z;
draw_vertex_sphere(xyz);
xyz[0] = 0.0;
xyz[1] = 0.0;
xyz[2] = 2.0*z;
draw_vertex_sphere(xyz);
glEnd();
}
void draw_lj_sphere_ij(int i, int iplus, int j, int jplus, t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE], int fade, double fade_value)
/* draw particles at simulation grid point (i,j) */
{
@@ -2343,6 +2679,7 @@ void draw_lj_sphere_3d(t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE])
if (observer[2] > 0.0)
{
if (DRAW_POLAR_AXIS) draw_axis(-1.0);
if (imin < imax)
{
for (i=imax; i>imid; i--)
@@ -2436,9 +2773,12 @@ void draw_lj_sphere_3d(t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE])
draw_lj_sphere_ij(i, i+1, NY_SPHERE-3, NY_SPHERE-1, wsphere, fade, fade_value);
draw_lj_sphere_ij(NX_SPHERE-1, 0, NY_SPHERE-3, NY_SPHERE-1, wsphere, fade, fade_value);
if (DRAW_POLAR_AXIS) draw_axis(1.0);
}
else
{
if (DRAW_POLAR_AXIS) draw_axis(1.0);
if (imin < imax)
{
for (i=imax; i>imid; i--)
@@ -2536,6 +2876,8 @@ void draw_lj_sphere_3d(t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE])
for (j=2; j>0; j--)
draw_lj_sphere_ij(NX_SPHERE-1, 0, j-1, j, wsphere, fade, fade_value);
if (DRAW_POLAR_AXIS) draw_axis(-1.0);
}
}
sub_maze.c
+135 -1
View File
@@ -8,6 +8,7 @@
#define MAZE_TYPE_CIRCLE 1 /* circular maze */
#define MAZE_TYPE_HEX 2 /* honeycomb maze */
#define MAZE_TYPE_OCT 3 /* maze with octagonal and square cells */
#define MAZE_TYPE_CYLINDER 4 /* square maze with periodic boundary conditions */
typedef struct
{
@@ -139,6 +140,128 @@ void init_maze_graph(t_maze maze[NXMAZE*NYMAZE])
}
}
void init_cylinder_maze_graph(t_maze maze[NXMAZE*NYMAZE])
{
int i, j, k, n;
printf("Initializing maze\n");
/* initialize neighbours */
/* in the bulk */
for (i=1; i<NXMAZE-1; i++)
for (j=1; j<NYMAZE-1; j++)
{
n = nmaze(i, j);
maze[n].nneighb = 4;
maze[n].neighb[0] = nmaze(i, j+1);
maze[n].neighb[1] = nmaze(i+1, j);
maze[n].neighb[2] = nmaze(i, j-1);
maze[n].neighb[3] = nmaze(i-1, j);
for (k=0; k<4; k++) maze[n].directions[k] = k;
}
/* left side */
for (j=1; j<NYMAZE-1; j++)
{
n = nmaze(0, j);
maze[n].nneighb = 4;
maze[n].neighb[0] = nmaze(0, j+1);
maze[n].neighb[1] = nmaze(1, j);
maze[n].neighb[2] = nmaze(0, j-1);
maze[n].neighb[3] = nmaze(NXMAZE-1, j);
for (k=0; k<4; k++) maze[n].directions[k] = k;
}
/* right side */
for (j=1; j<NYMAZE-1; j++)
{
n = nmaze(NXMAZE-1, j);
maze[n].nneighb = 4;
maze[n].neighb[0] = nmaze(NXMAZE-1, j+1);
maze[n].neighb[1] = nmaze(NXMAZE-2, j);
maze[n].neighb[2] = nmaze(NXMAZE-1, j-1);
maze[n].neighb[3] = nmaze(0,j);
maze[n].directions[0] = 0;
maze[n].directions[1] = 3;
maze[n].directions[2] = 2;
maze[n].directions[3] = 1;
}
/* bottom side */
for (i=1; i<NXMAZE-1; i++)
{
n = nmaze(i, 0);
maze[n].nneighb = 3;
maze[n].neighb[0] = nmaze(i, 1);
maze[n].neighb[1] = nmaze(i+1, 0);
maze[n].neighb[2] = nmaze(i-1, 0);
maze[n].directions[0] = 0;
maze[n].directions[1] = 1;
maze[n].directions[2] = 3;
}
/* top side */
for (i=1; i<NXMAZE-1; i++)
{
n = nmaze(i, NYMAZE-1);
maze[n].nneighb = 3;
maze[n].neighb[0] = nmaze(i, NYMAZE-2);
maze[n].neighb[1] = nmaze(i+1, NYMAZE-1);
maze[n].neighb[2] = nmaze(i-1, NYMAZE-1);
maze[n].directions[0] = 2;
maze[n].directions[1] = 1;
maze[n].directions[2] = 3;
}
/* corners */
n = nmaze(0,0);
maze[n].nneighb = 3;
maze[n].neighb[0] = nmaze(1,0);
maze[n].neighb[1] = nmaze(0,1);
maze[n].neighb[2] = nmaze(NXMAZE-1,0);
maze[n].directions[0] = 1;
maze[n].directions[1] = 0;
maze[n].directions[2] = 3;
n = nmaze(NXMAZE-1,0);
maze[n].nneighb = 3;
maze[n].neighb[0] = nmaze(NXMAZE-2,0);
maze[n].neighb[1] = nmaze(NXMAZE-1,1);
maze[n].neighb[2] = nmaze(0,0);
maze[n].directions[0] = 3;
maze[n].directions[1] = 0;
maze[n].directions[2] = 1;
n = nmaze(0,NYMAZE-1);
maze[n].nneighb = 3;
maze[n].neighb[0] = nmaze(1,NYMAZE-1);
maze[n].neighb[1] = nmaze(0,NYMAZE-2);
maze[n].neighb[2] = nmaze(NXMAZE-1,NYMAZE-1);
maze[n].directions[0] = 1;
maze[n].directions[1] = 2;
maze[n].directions[2] = 3;
n = nmaze(NXMAZE-1,NYMAZE-1);
maze[n].nneighb = 3;
maze[n].neighb[0] = nmaze(NXMAZE-2,NYMAZE-1);
maze[n].neighb[1] = nmaze(NXMAZE-1,NYMAZE-2);
maze[n].neighb[2] = nmaze(0,NYMAZE-1);
maze[n].directions[0] = 3;
maze[n].directions[1] = 2;
maze[n].directions[2] = 1;
/* initialize other parameters */
for (i=0; i<NXMAZE; i++)
for (j=0; j<NYMAZE; j++)
{
n = nmaze(i, j);
maze[n].active = 0;
maze[n].tested = 0;
maze[n].connected = 0;
maze[n].closed = 0;
maze[n].north = 1;
maze[n].east = 1;
maze[n].south = 1;
maze[n].west = 1;
}
}
void init_circular_maze_graph(t_maze maze[NXMAZE*NYMAZE])
/* initialise graph of circular maze */
@@ -738,7 +861,7 @@ int find_maze_path(t_maze maze[NXMAZE*NYMAZE], int n0, int *path, int *pathlengt
inext = next_table[rand()%nnext];
nextcell = maze[n].neighb[inext];
/* square and circular maze */
if (mazetype < MAZE_TYPE_HEX) switch(maze[n].directions[inext]){
if ((mazetype < MAZE_TYPE_HEX)||(mazetype == MAZE_TYPE_CYLINDER)) switch(maze[n].directions[inext]){
case(0):
{
printf("Moving north\n");
@@ -957,6 +1080,11 @@ void init_maze_oftype(t_maze maze[NXMAZE*NYMAZE], int type)
init_oct_maze_graph(maze);
break;
}
case (MAZE_TYPE_CYLINDER):
{
init_cylinder_maze_graph(maze);
break;
}
}
for (i=0; i<RAND_SHIFT; i++) rand();
@@ -1036,6 +1164,12 @@ void init_oct_maze(t_maze maze[NXMAZE*NYMAZE])
init_maze_oftype(maze, MAZE_TYPE_OCT);
}
void init_cylinder_maze(t_maze maze[NXMAZE*NYMAZE])
/* init a maze with hexagonal cells */
{
init_maze_oftype(maze, MAZE_TYPE_CYLINDER);
}
void init_maze_exit(int nx, int ny, t_maze maze[NXMAZE*NYMAZE])
/* init a maze with exit at (nx, ny) */
{
sub_wave.c
+851 -3
View File
@@ -303,6 +303,12 @@ void write_text( double x, double y, char *st)
return(alph);
}
short int signp(double x) /* sign with sign(0) = 1 */
{
if (x >= 0.0) return(1);
else return(-1);
}
double iabs(int i) /* absolute value */
{
int res;
@@ -890,6 +896,27 @@ int init_circle_config_pattern(t_circle circles[NMAXCIRCLES], int circle_pattern
}
break;
}
case (C_CLOAK_OBJECT):
{
for (i = 0; i < 40; i++)
for (j = 0; j < 5; j++)
{
n = 5*i + j;
phi = (double)i*DPI/40.0;
r = LAMBDA*0.5*(1.0 + (double)j/5.0);
circles[n].xc = r*cos(phi);
circles[n].yc = r*sin(phi);
circles[n].radius = MU;
circles[n].active = 1;
}
n = ncircles;
circles[n].xc = 0.0;
circles[n].yc = 0.0;
circles[n].radius = 0.25*LAMBDA;
circles[n].active = 1;
ncircles++;
break;
}
case (C_CLOAK_A): /* optimized model A1 by C. Jo et al */
{
ncircles = 200;
@@ -911,6 +938,27 @@ int init_circle_config_pattern(t_circle circles[NMAXCIRCLES], int circle_pattern
}
break;
}
case (C_CLOAK_AB): /* optimized model A1 by C. Jo et al with smaller discs */
{
ncircles = 200;
ra[0] = 0.0731; sa[0] = 1.115;
ra[1] = 0.0768; sa[1] = 1.292;
ra[2] = 0.0652; sa[2] = 1.464;
ra[3] = 0.056; sa[3] = 1.633;
ra[4] = 0.0375; sa[4] = 1.794;
for (i = 0; i < 40; i++)
for (j = 0; j < 5; j++)
{
n = 5*i + j;
phi = (double)i*DPI/40.0;
r = LAMBDA*sa[j];
circles[n].xc = r*cos(phi);
circles[n].yc = r*sin(phi);
circles[n].radius = 0.6*LAMBDA*ra[j];
circles[n].active = 1;
}
break;
}
case (C_LASER):
{
ncircles = 17;
@@ -978,6 +1026,11 @@ int init_circle_config_pattern(t_circle circles[NMAXCIRCLES], int circle_pattern
{
/* new circle is far away from circle k */
far = far*((x - circles[k].xc)*(x - circles[k].xc) + (y - circles[k].yc)*(y - circles[k].yc) >= dpoisson*dpoisson);
if ((B_COND == BC_PERIODIC)||(B_COND == BC_VPER_HABS))
{
far *= ((x - circles[k].xc)*(x - circles[k].xc) + (y + YMAX - YMIN - circles[k].yc)*(y + YMAX - YMIN - circles[k].yc) >= dpoisson*dpoisson);
far *= ((x - circles[k].xc)*(x - circles[k].xc) + (y - YMAX + YMIN - circles[k].yc)*(y - YMAX + YMIN - circles[k].yc) >= dpoisson*dpoisson);
}
/* new circle is in domain */
far = far*(vabs(x) < LAMBDA)*(y < YMAX)*(y > YMIN);
}
@@ -1004,6 +1057,33 @@ int init_circle_config_pattern(t_circle circles[NMAXCIRCLES], int circle_pattern
printf("%i active circles\n", n_p_active);
}
/* deal with periodic b.c. */
if ((B_COND == BC_PERIODIC)||(B_COND == BC_VPER_HABS))
{
ncirc0 = ncircles;
for (i=0; i<ncirc0; i++)
{
x = circles[i].xc;
y = circles[i].yc;
if (y > YMAX - MU)
{
circles[ncircles].xc = x;
circles[ncircles].yc = y - YMAX + YMIN;
circles[ncircles].radius = MU;
circles[ncircles].active = 1;
ncircles++;
}
else if (y < YMIN + MU)
{
circles[ncircles].xc = x;
circles[ncircles].yc = y + YMAX - YMIN;
circles[ncircles].radius = MU;
circles[ncircles].active = 1;
ncircles++;
}
}
}
printf("Generated %i circles\n", ncircles);
break;
}
@@ -6145,6 +6225,18 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_
}
return(1);
}
case (D_CIRCLES_NEUMANN):
{
for (i = 0; i < ncirc; i++)
if (circles[i].active)
{
x1 = circles[i].xc;
y1 = circles[i].yc;
r2 = circles[i].radius*circles[i].radius;
if ((x-x1)*(x-x1) + (y-y1)*(y-y1) < r2) return(0);
}
return(1);
}
case (D_CIRCLES_IN_RECT): /* returns 2 inside circles, 0 outside rectangle */
{
for (i = 0; i < ncirc; i++)
@@ -6153,9 +6245,10 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_
x1 = circles[i].xc;
y1 = circles[i].yc;
r2 = circles[i].radius*circles[i].radius;
if ((x-x1)*(x-x1) + (y-y1)*(y-y1) < r2) return(2);
if ((x-x1)*(x-x1) + (y-y1)*(y-y1) < r2) return(0);
// if ((x-x1)*(x-x1) + (y-y1)*(y-y1) < r2) return(2);
}
if ((vabs(x) >= LAMBDA)||(vabs(y) >= 1.0)) return(0);
if ((vabs(x) >= LAMBDA)||(vabs(y) >= YMAX - XMAX + LAMBDA)) return(0);
else return(1);
}
case (D_POLYGONS):
@@ -6956,6 +7049,37 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_
}
return(0);
}
case (D_TWO_ELLIPSES):
{
x1 = vabs(x);
return((x1 - LAMBDA)*(x1 - LAMBDA)/(MU*MU) + y*y < 1.0);
}
case (D_TWO_RECTANGLES):
{
x1 = vabs(x);
return((vabs(x1 - LAMBDA) < MU)&&(vabs(y) < 1.0));
}
case (D_TWO_POLYGONS):
{
// if (x < 0.0) y = -y;
x1 = vabs(x);
return(in_polygon(x1 - LAMBDA, y, MU, NPOLY, APOLY));
}
case (D_TWO_POLY_PARABOLAS):
{
condition = 1;
omega = DPI/((double)NPOLY);
if (x > 0.0) x -= LAMBDA;
else x += LAMBDA;
for (k=0; k<NPOLY; k++)
{
angle = APOLY*PID + ((double)k+0.5)*omega;
x1 = x*cos(angle) + y*sin(angle);
y1 = -x*sin(angle) + y*cos(angle);
condition = condition*(x1 < MU - 0.25*y1*y1/MU);
}
return(condition);
}
case (D_CARDIOID):
{
y1 = y - MU;
@@ -9350,12 +9474,19 @@ void draw_billiard(int fade, double fade_value) /* draws the billiard bound
if (circles[i].active) draw_circle(circles[i].xc, circles[i].yc, circles[i].radius, NSEG);
break;
}
case (D_CIRCLES_NEUMANN):
{
glLineWidth(BOUNDARY_WIDTH);
for (i = 0; i < ncircles; i++)
if (circles[i].active) draw_circle(circles[i].xc, circles[i].yc, circles[i].radius, NSEG);
break;
}
case (D_CIRCLES_IN_RECT):
{
glLineWidth(BOUNDARY_WIDTH);
for (i = 0; i < ncircles; i++)
if (circles[i].active) draw_circle(circles[i].xc, circles[i].yc, circles[i].radius, NSEG);
draw_rectangle(-LAMBDA, -1.0, LAMBDA, 1.0);
draw_rectangle(-LAMBDA, -(YMAX - XMAX + LAMBDA), LAMBDA, YMAX - XMAX + LAMBDA);
if ((FOCI)&&(CIRCLE_PATTERN == C_LASER))
{
glColor3f(0.3, 0.3, 0.3);
@@ -9976,6 +10107,109 @@ void draw_billiard(int fade, double fade_value) /* draws the billiard bound
break;
}
case (D_TWO_ELLIPSES):
{
glBegin(GL_LINE_LOOP);
for (i=0; i<NSEG; i++)
{
phi = DPI*(double)i/(double)NSEG;
xy_to_pos(LAMBDA + MU*cos(phi), sin(phi), pos);
glVertex2d(pos[0], pos[1]);
}
glEnd();
glBegin(GL_LINE_LOOP);
for (i=0; i<NSEG; i++)
{
phi = DPI*(double)i/(double)NSEG;
xy_to_pos(-LAMBDA + MU*cos(phi), sin(phi), pos);
glVertex2d(pos[0], pos[1]);
}
glEnd();
/* draw foci */
if (FOCI)
{
if (fade) glColor3f(0.3*fade_value, 0.3*fade_value, 0.3*fade_value);
else glColor3f(0.3, 0.3, 0.3);
x0 = sqrt(1.0 - MU*MU);
glLineWidth(2);
glEnable(GL_LINE_SMOOTH);
draw_circle(LAMBDA, x0, 0.01, NSEG);
draw_circle(LAMBDA, -x0, 0.01, NSEG);
draw_circle(-LAMBDA, x0, 0.01, NSEG);
draw_circle(-LAMBDA, -x0, 0.01, NSEG);
}
break;
}
case (D_TWO_RECTANGLES):
{
draw_rectangle(LAMBDA - MU, -1.0, LAMBDA + MU, 1.0);
draw_rectangle(-LAMBDA - MU, -1.0, -LAMBDA + MU, 1.0);
break;
}
case (D_TWO_POLYGONS):
{
omega = DPI/((double)NPOLY);
glBegin(GL_LINE_LOOP);
for (i=0; i<=NPOLY; i++)
{
x = LAMBDA + MU*cos(i*omega + APOLY*PID);
y = MU*sin(i*omega + APOLY*PID);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
glBegin(GL_LINE_LOOP);
for (i=0; i<=NPOLY; i++)
{
x = -LAMBDA + MU*cos(i*omega + APOLY*PID);
y = MU*sin(i*omega + APOLY*PID);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
break;
}
case (D_TWO_POLY_PARABOLAS):
{
omega = PI/((double)NPOLY);
a = 0.25/MU;
b = 1.0/tan(omega);
c = MU;
ymax = (-b + sqrt(b*b + 4.0*a*c))/(2.0*a);
dy = 2.0*ymax/(double)NSEG;
for (k1 = -1; k1 < 2; k1+=2)
{
glBegin(GL_LINE_LOOP);
for (k=0; k<NPOLY; k++)
{
alpha = APOLY*PID + (2.0*(double)k+1.0)*omega;
for (i = 0; i < NSEG+1; i++)
{
y1 = -ymax + dy*(double)i;
x1 = MU - 0.25*y1*y1/MU;
x = (double)k1*LAMBDA + x1*cos(alpha) - y1*sin(alpha);
y = x1*sin(alpha) + y1*cos(alpha);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
}
glEnd ();
}
if (FOCI)
{
glColor3f(0.3, 0.3, 0.3);
draw_circle(LAMBDA, 0.0, r, NSEG);
draw_circle(-LAMBDA, 0.0, r, NSEG);
}
break;
}
case (D_CARDIOID):
{
glBegin(GL_LINE_LOOP);
@@ -11915,6 +12149,15 @@ double oscillating_bc(int time, int j)
if (t1 > 0.0) amp += cos(t1);
return(AMPLITUDE*amp);
}
case (OSC_TWOFREQ_TOPBOT):
{
if (j > NY/2)
t = (double)time*OMEGA + (double)(NY-j)*OSCIL_LEFT_YSHIFT/(double)NY;
else
t = (double)time*OMEGA_B + (double)(NY-j)*OSCIL_LEFT_YSHIFT/(double)NY;
return(AMPLITUDE*cos(t)*exp(-(double)t*DAMPING));
}
}
}
@@ -13138,5 +13381,610 @@ int init_input_signal()
printf("Tail frames: %i\n", tailcounter);
}
int compute_normal_vector(int i, int j, int n, t_neumann_bc neumann_bc[N_NEUMANN_POINTS], t_circle circles[NMAXCIRCLES])
/* compute normal vector to domain at boundary point */
{
int on_boundary = 1, side, k, circle;
double xy[2], v1, v2, w1, w2, sum, angle, theta, x1, y1, c, w, distance, dist0, dx, dy;
dx = (XMAX-XMIN)/(double)NX;
dy = (YMAX-YMIN)/(double)NY;
/* on_boundary = 0 is used to optionally declare some boundary points as Dirichlet */
ij_to_xy(i, j, xy);
switch (B_DOMAIN) {
case (D_ELLIPSE) :
{
v1 = 2.0*xy[0]/(LAMBDA*LAMBDA);
v2 = 2.0*xy[1];
break;
}
case (D_POLYGON):
{
theta = DPI/(double)NPOLY;
angle = argument(xy[0], xy[1]) - PID*APOLY + DPI;
side = (int)(angle/theta);
angle = PID*APOLY + 0.5*theta + (double)side*theta;
v1 = cos(angle);
v2 = sin(angle);
break;
}
case (D_POLY_PARABOLAS):
{
theta = DPI/(double)NPOLY;
angle = argument(xy[0], xy[1]) - PID*APOLY;
if (angle < 0.0) angle += DPI;
side = (int)((angle+DPI)/theta);
side = side%NPOLY;
angle = PID*APOLY + 0.5*theta + (double)side*theta;
y1 = -(xy[0]-LAMBDA)*sin(angle) + xy[1]*cos(angle);
w1 = 1.0;
w2 = y1/(2.0*MU);
v1 = w1*cos(angle) - w2*sin(angle);
v2 = w1*sin(angle) + w2*cos(angle);
break;
}
case (D_PENROSE):
{
x1 = xy[0];
y1 = xy[1];
c = sqrt(LAMBDA*LAMBDA - (1.0-MU)*(1.0-MU));
w = 0.1*MU;
/* upper ellipse */
if (y1 > MU)
{
v1 = 2.0*x1/(LAMBDA*LAMBDA);
v2 = 2.0*(y1-MU)/((1.0-MU)*(1.0-MU));
}
/* lower ellipse */
else if (y1 < -MU)
{
v1 = 2.0*x1/(LAMBDA*LAMBDA);
v2 = 2.0*(y1+MU)/((1.0-MU)*(1.0-MU));
}
/* right ellipse */
else if ((x1 > 0.0)&&(x1 < LAMBDA))
{
v1 = -8.0*(x1-c)/(MU*MU);
v2 = -2.0*y1/(MU*MU);
}
/* left ellipse */
else if ((x1 < 0.0)&&(x1 > -LAMBDA))
{
v1 = -8.0*(x1+c)/(MU*MU);
v2 = -2.0*y1/(MU*MU);
}
else if (y1 > 0.0)
{
/* top right straight parts */
if (x1 > 0.0)
{
if (y1 <= w)
{
v1 = 0.0;
v2 = -1.0;
}
else if (x1 >= 0.5*(LAMBDA+c))
{
v1 = 1.0;
v2 = 0.0;
}
else
{
v1 = -1.0;
v2 = 0.0;
}
}
/* top left straight parts */
else
{
if (y1 <= w)
{
v1 = 0.0;
v2 = -1.0;
}
else if (x1 <= -0.5*(LAMBDA+c))
{
v1 = -1.0;
v2 = 0.0;
}
else
{
v1 = 1.0;
v2 = 0.0;
}
}
}
else
{
/* bottom right straight parts */
if (x1 > 0.0)
{
if (y1 >= -w)
{
v1 = 0.0;
v2 = 1.0;
}
else if (x1 >= 0.5*(LAMBDA+c))
{
v1 = 1.0;
v2 = 0.0;
}
else
{
v1 = -1.0;
v2 = 0.0;
}
}
/* bottom left straight parts */
else
{
if (y1 >= -w)
{
v1 = 0.0;
v2 = 1.0;
}
else if (x1 <= -0.5*(LAMBDA+c))
{
v1 = -1.0;
v2 = 0.0;
}
else
{
v1 = 1.0;
v2 = 0.0;
}
}
}
break;
}
case (D_CIRCLES):
{
/* find closest circle */
distance = 10.0;
for (k=0; k<ncircles; k++) if (circles[k].active)
{
dist0 = module2(xy[0] - circles[k].xc, xy[1] - circles[k].yc);
if (dist0 < distance)
{
distance = dist0;
circle = k;
}
}
v1 = circles[circle].xc - xy[0];
v2 = circles[circle].yc - xy[1];
break;
}
case (D_CIRCLES_NEUMANN):
{
if (xy[1] < 0.0)
{
/* find closest circle */
distance = 10.0;
for (k=0; k<ncircles; k++) if (circles[k].active)
{
dist0 = module2(xy[0] - circles[k].xc, xy[1] - circles[k].yc);
if (dist0 < distance)
{
distance = dist0;
circle = k;
}
}
v1 = circles[circle].xc - xy[0];
v2 = circles[circle].yc - xy[1];
}
else on_boundary = 0;
break;
}
case (D_CIRCLES_IN_RECT):
{
/* points on outer boundary */
if (xy[0] <= -LAMBDA + dx)
{
v1 = -1.0;
v2 = 0.0;
}
else if (xy[0] >= LAMBDA - dx)
{
v1 = 1.0;
v2 = 0.0;
}
else if (xy[1] <= -(YMAX - XMAX + LAMBDA) + dy)
{
v1 = 0.0;
v2 = -1.0;
}
else if (xy[1] >= YMAX - XMAX + LAMBDA - dy)
{
v1 = 0.0;
v2 = 1.0;
}
else
{
/* find closest circle */
distance = 10.0;
for (k=0; k<ncircles; k++)
{
dist0 = module2(xy[0] - circles[k].xc, xy[1] - circles[k].yc);
if (dist0 < distance)
{
distance = dist0;
circle = k;
}
}
v1 = circles[circle].xc - xy[0];
v2 = circles[circle].yc - xy[1];
}
break;
}
case (D_TWO_ELLIPSES):
{
if (xy[0] > 0.0)
{
v1 = 2.0*(xy[0] - LAMBDA)/(MU*MU);
v2 = 2.0*xy[1];
}
else on_boundary = 0;
break;
}
case (D_TWO_RECTANGLES):
{
if (xy[0] > 0.0)
{
if (vabs(xy[1]) < 1.0)
{
if (xy[0] > LAMBDA) v1 = 1.0;
else v1 = -1.0;
v2 = 0.0;
}
else
{
v1 = 0.0;
if (xy[1] > 0.0) v2 = 1.0;
else v2 = -1.0;
}
}
else on_boundary = 0;
break;
}
case (D_TWO_POLYGONS):
{
if (xy[0] > 0.0)
{
theta = DPI/(double)NPOLY;
angle = argument(xy[0] - LAMBDA, xy[1]) - PID*APOLY + DPI;
side = (int)(angle/theta);
angle = PID*APOLY + 0.5*theta + (double)side*theta;
v1 = cos(angle);
v2 = sin(angle);
}
else on_boundary = 0;
break;
}
case (D_TWO_POLY_PARABOLAS):
{
if (xy[0] > 0.0)
{
theta = DPI/(double)NPOLY;
angle = argument(xy[0] - LAMBDA, xy[1]) - PID*APOLY;
if (angle < 0.0) angle += DPI;
printf("(x,y) = (%.3lg,%.3lg), angle = %.3lg ", xy[0] - LAMBDA, xy[1], angle*180.0/PI);
side = (int)((angle+DPI)/theta);
side = side%NPOLY;
printf("side %i ", side);
// angle += (double)side*theta;
angle = PID*APOLY + 0.5*theta + (double)side*theta;
// angle = PID*APOLY + (double)side*theta;
printf("angle = %.5lg\n", angle*180.0/PI);
// theta *= (double)side;
// x1 = (xy[0]-LAMBDA)*cos(angle) + xy[1]*sin(angle);
y1 = -(xy[0]-LAMBDA)*sin(angle) + xy[1]*cos(angle);
w1 = 1.0;
w2 = y1/(2.0*MU);
v1 = w1*cos(angle) - w2*sin(angle);
v2 = w1*sin(angle) + w2*cos(angle);
printf("(w1, w2) = (%.3lg, %.3lg) (v1, v2) = (%.3lg, %.3lg)\n", w1, w2, v1, v2);
// angle = argument(xy[0] - LAMBDA, xy[1]);
// v1 = cos(angle);
// v2 = sin(angle);
}
else on_boundary = 0;
break;
}
default:
{
printf("Function compute_normal vector not defined for this domain");
}
}
if (on_boundary)
{
neumann_bc[n].signx1 = -signp(v1);
neumann_bc[n].signy1 = 0;
neumann_bc[n].signx2 = 0;
neumann_bc[n].signy2 = -signp(v2);
v1 = vabs(v1);
v2 = vabs(v2);
sum = v1 + v2;
neumann_bc[n].nx = v1/sum;
neumann_bc[n].ny = v2/sum;
// printf("Point (%.3lg, %.3lg): n = (%.3lg, %.3lg), sign = (%i,%i)\n", xy[0], xy[1], v1/sum, v2/sum, neumann_bc[n].signx, neumann_bc[n].signy);
}
neumann_bc[n].top = (xy[1] > 0.0);
return (on_boundary);
}
void build_neumann_data(t_neumann_bc neumann_bc[N_NEUMANN_POINTS], int n)
/* compute correct values of Neumann boundary data */
{
int p, i, j;
short int signx1, signy2;
double n1, n2;
for (p = 0; p < n; p++)
{
i = neumann_bc[p].i;
j = neumann_bc[p].j;
signx1 = neumann_bc[p].signx1;
signy2 = neumann_bc[p].signy2;
n1 = neumann_bc[p].nx;
n2 = neumann_bc[p].ny;
if (n1 > n2) /* boundary is closer to vertical */
{
neumann_bc[p].nx = (n1-n2)/n1;
neumann_bc[p].ny = n2/n1;
neumann_bc[p].signx2 = signx1;
}
else if (n2 > n1) /* boundary is closer to horizontal */
{
neumann_bc[p].nx = n1/n2;
neumann_bc[p].ny = (n2-n1)/n2;
neumann_bc[p].signy1 = signy2;
}
/* compute absolute coordinates of neighbours */
neumann_bc[p].nbi1 = i + neumann_bc[p].signx1;
neumann_bc[p].nbj1 = j + neumann_bc[p].signy1;
neumann_bc[p].nbi2 = i + neumann_bc[p].signx2;
neumann_bc[p].nbj2 = j + neumann_bc[p].signy2;
/* make points active */
neumann_bc[p].active = 1;
// printf("Neighbouring points in domain: (%i, %i), n = (%.3lg, %.3lg)\n", xy_in[i+signx1][j+signy1], xy_in[i+signx2][j+signy2], neumann_bc[p].nx, neumann_bc[p].ny);
}
}
void adapt_neumann_data_to_bc(t_neumann_bc neumann_bc[N_NEUMANN_POINTS], int n, short int comparison)
/* adapt Neumann boundary data to boundary conditions */
{
int p, i, j, jmid = NY/2;
/* x direction */
if (B_COND == BC_PERIODIC) for (p=0; p<n; p++)
{
if (neumann_bc[p].nbi1 == NX) neumann_bc[p].nbi1 == 0;
else if (neumann_bc[p].nbi1 == -1) neumann_bc[p].nbi1 == NX-1;
if (neumann_bc[p].nbi2 == NX) neumann_bc[p].nbi2 == 0;
else if (neumann_bc[p].nbi2 == -1) neumann_bc[p].nbi2 == NX-1;
}
else for (p=0; p<n; p++)
{
if (neumann_bc[p].nbi1 == NX) neumann_bc[p].nbi1 == NX-1;
else if (neumann_bc[p].nbi1 == -1) neumann_bc[p].nbi1 == 0;
if (neumann_bc[p].nbi2 == NX) neumann_bc[p].nbi2 == NX-1;
else if (neumann_bc[p].nbi2 == -1) neumann_bc[p].nbi2 == 0;
}
/* y direction */
if ((B_COND == BC_PERIODIC)||(B_COND == BC_VPER_HABS)||(B_COND == BC_ABS_REFLECT))
{
if (comparison) for (p=0; p<n; p++)
{
if (neumann_bc[p].nbj1 == NY) neumann_bc[p].nbj1 == jmid;
else if (neumann_bc[p].nbj1 == -1) neumann_bc[p].nbj1 == jmid-1;
// else if ((neumann_bc[p].nbj1 == jmid-1)&&(neumann_bc[p].signy1 == -1))
// neumann_bc[p].nbj1 == NY-1;
// else if ((neumann_bc[p].nbj1 == jmid)&&(neumann_bc[p].signy1 == 1))
// neumann_bc[p].nbj1 == 0;
if (neumann_bc[p].nbj2 == NY) neumann_bc[p].nbj2 == jmid;
else if (neumann_bc[p].nbj2 == -1) neumann_bc[p].nbj2 == jmid-1;
// else if ((neumann_bc[p].nbj2 == jmid-1)&&(neumann_bc[p].signy1 == -1))
// neumann_bc[p].nbj2 == NY-1;
// else if ((neumann_bc[p].nbj2 == jmid)&&(neumann_bc[p].signy1 == 1))
// neumann_bc[p].nbj2 == 0;
/* TEST */
if ((neumann_bc[p].nbj1 == jmid)||(neumann_bc[p].nbj1 == jmid-1))
neumann_bc[p].active = 0;
if ((neumann_bc[p].nbj2 == jmid)||(neumann_bc[p].nbj2 == jmid-1))
neumann_bc[p].active = 0;
}
else for (p=0; p<n; p++)
{
if (neumann_bc[p].nbj1 == NY) neumann_bc[p].nbj1 == 0;
else if (neumann_bc[p].nbj1 == -1) neumann_bc[p].nbj1 == NY-1;
if (neumann_bc[p].nbj2 == NY) neumann_bc[p].nbj2 == 0;
else if (neumann_bc[p].nbj2 == -1) neumann_bc[p].nbj2 == NY-1;
}
}
else
{
if (neumann_bc[p].nbj1 == NY) neumann_bc[p].nbj1 == NY-1;
else if (neumann_bc[p].nbj1 == -1) neumann_bc[p].nbj1 == 0;
if (neumann_bc[p].nbj2 == NY) neumann_bc[p].nbj2 == NY-1;
else if (neumann_bc[p].nbj2 == -1) neumann_bc[p].nbj2 == 0;
if (comparison)
{
if ((neumann_bc[p].nbj1 == jmid-1)&&(neumann_bc[p].signy1 == -1))
neumann_bc[p].nbj1 == jmid;
else if ((neumann_bc[p].nbj1 == jmid)&&(neumann_bc[p].signy1 == 1))
neumann_bc[p].nbj1 == jmid-1;
}
}
}
void extend_neumann_data_flat(t_neumann_bc neumann_bc[N_NEUMANN_POINTS], int n)
/* extend Neumann boundary data to flat table */
{
int p;
for (p=0; p<n; p++)
{
neumann_bc[p].k = (neumann_bc[p].i)*NY + neumann_bc[p].j;
neumann_bc[p].nbk1 = (neumann_bc[p].nbi1)*NY + neumann_bc[p].nbj1;
neumann_bc[p].nbk2 = (neumann_bc[p].nbi2)*NY + neumann_bc[p].nbj2;
}
}
void test_neumann_data(short int * xy_in[NX], t_neumann_bc neumann_bc[N_NEUMANN_POINTS], int n)
/* test values of Neumann boundary data, for debugging */
{
int p, i, j, errors = 0;
short int signx1, signy1, signx2, signy2;
double xy[2];
for (p = 0; p < n; p++)
{
i = neumann_bc[p].i;
j = neumann_bc[p].j;
signx1 = neumann_bc[p].signx1;
signy1 = neumann_bc[p].signy1;
signx2 = neumann_bc[p].signx2;
signy2 = neumann_bc[p].signy2;
ij_to_xy(i, j, xy);
// printf("Signs: (%i, %i), (%i, %i)\n", signx1, signy1, signx2, signy2);
// printf("pnt %i of %i, coords = (%.3lg, %.3lg), nbh pts in domain: (%i, %i), n = (%.3lg, %.3lg), sign1 = (%i, %i), sign2 = (%i, %i)\n", p, n, xy[0], xy[1], xy_in[i+signx1][j+signy1], xy_in[i+signx2][j+signy2], neumann_bc[p].nx, neumann_bc[p].ny, signx1, signy1, signx2, signy2);
if (xy_in[i+signx1][j+signy1] == -1)
{
neumann_bc[p].active = 0;
errors++;
}
if (xy_in[i+signx2][j+signy2] == -1)
{
neumann_bc[p].active = 0;
errors++;
}
}
printf("%i errors among %i points\n", errors, n);
}
void test_neumann_data_flat(short int xy_in[NX*NY], t_neumann_bc neumann_bc[N_NEUMANN_POINTS], int n)
/* test values of Neumann boundary data, for debugging */
{
int p, i, j, errors = 0;
short int signx1, signy1, signx2, signy2;
double xy[2];
for (p = 0; p < n; p++)
{
i = neumann_bc[p].i;
j = neumann_bc[p].j;
signx1 = neumann_bc[p].signx1;
signy1 = neumann_bc[p].signy1;
signx2 = neumann_bc[p].signx2;
signy2 = neumann_bc[p].signy2;
ij_to_xy(i, j, xy);
// printf("Signs: (%i, %i), (%i, %i)\n", signx1, signy1, signx2, signy2);
printf("pnt %i of %i, coords = (%.3lg, %.3lg), nbh pts in domain: (%i, %i), n = (%.3lg, %.3lg), sign1 = (%i, %i), sign2 = (%i, %i)\n", p, n, xy[0], xy[1], xy_in[(i+signx1)*NY+j+signy1], xy_in[(i+signx2)*NY+j+signy2], neumann_bc[p].nx, neumann_bc[p].ny, signx1, signy1, signx2, signy2);
if (xy_in[(i+signx1)*NY+j+signy1] == -1) errors++;
if (xy_in[(i+signx2)*NY+j+signy2] == -1) errors++;
}
printf("%i errors among %i points\n", errors, n);
}
int init_neumann_bc(short int * xy_in[NX], t_neumann_bc neumann_bc[N_NEUMANN_POINTS], t_circle circles[NMAXCIRCLES], short int comparison)
/* initialize data for Neumann boundary conditions */
/* points are numbered in order of appearance, and treated as such by evolve_wave() */
{
int i, j, n = 0;
double vn[2];
printf("Building Neumann boundary data\n");
/* search points in the bulk that are not in the domain */
for (i=1; i<NX-1; i++)
for (j=1; j<NY-1; j++) if (!xy_in[i][j])
{
if ((xy_in[i-1][j]==1)||(xy_in[i+1][j]==1)||(xy_in[i][j-1]==1)||(xy_in[i][j+1]==1))
/* the point has a neighbour in the domain */
{
if (compute_normal_vector(i, j, n, neumann_bc, circles))
{
xy_in[i][j] = -1;
neumann_bc[n].i = i;
neumann_bc[n].j = j;
n++;
}
if (n >= N_NEUMANN_POINTS)
{
printf("Error: N_NEUMANN_POINTS should be larger than %i\n", n);
exit(0);
}
}
}
printf("Found %i boundary points\n", n);
/* compute the correct coefficients and points */
build_neumann_data(neumann_bc, n);
adapt_neumann_data_to_bc(neumann_bc, n, comparison);
test_neumann_data(xy_in, neumann_bc, n);
return(n);
}
int init_neumann_bc_flat(short int xy_in[NX*NY], t_neumann_bc neumann_bc[N_NEUMANN_POINTS], t_circle circles[NMAXCIRCLES], int comparison)
/* initialize data for Neumann boundary conditions, version for flat tables */
/* points are numbered in order of appearance, and treated as such by evolve_wave() */
{
int i, j, n = 0;
double vn[2];
printf("Building Neumann boundary data\n");
/* search points in the bulk that are not in the domain */
for (i=1; i<NX-1; i++)
for (j=1; j<NY-1; j++) if (!xy_in[i*NY+j])
{
if ((xy_in[(i-1)*NY+j]==1)||(xy_in[(i+1)*NY+j]==1)||(xy_in[i*NY+j-1]==1)||(xy_in[i*NY+j+1]==1))
/* the point has a neighbour in the domain */
{
if (compute_normal_vector(i, j, n, neumann_bc, circles))
{
xy_in[i*NY+j] = -1;
neumann_bc[n].i = i;
neumann_bc[n].j = j;
// neumann_bc[n].k = i*NY+j;
n++;
}
if (n >= N_NEUMANN_POINTS)
{
printf("Error: N_NEUMANN_POINTS should be larger than %i\n", n);
exit(0);
}
}
}
printf("Found %i boundary points\n", n);
/* compute the correct coefficients and points */
build_neumann_data(neumann_bc, n);
adapt_neumann_data_to_bc(neumann_bc, n, comparison);
extend_neumann_data_flat(neumann_bc, n);
test_neumann_data_flat(xy_in, neumann_bc, n);
return(n);
}
sub_wave_3d.c
+1 -1
View File
@@ -1133,7 +1133,7 @@ void compute_field_color(double value, double value2, int cplot, int palette, do
switch (cplot) {
case (P_3D_AMP_ANGLE):
{
color_scheme_palette(COLOR_SCHEME, palette, VSCALE_AMPLITUDE*value, 1.0, 0, rgb);
color_scheme_palette(COLOR_SCHEME, palette, VSHIFT_AMPLITUDE + VSCALE_AMPLITUDE*value, 1.0, 0, rgb);
break;
}
case (P_3D_ENERGY):
sub_wave_comp.c
+256 -11
View File
@@ -8,11 +8,21 @@ void init_circle_config_half(int pattern, int top, t_circle circles[NMAXCIRCLES]
/* initialise the arrays circlex, circley, circlerad and circleactive */
/* for billiard shape D_CIRCLES */
{
int i, j, k, n, ncirc0, n_p_active, ncandidates=5000, naccepted, nnew;
double dx, dy, p, phi, r, r0, ra[5], sa[5], height, x, y = 0.0, gamma, ymean, ytop, ybottom, dpoisson = 3.05*MU, xmax;
int i, j, k, n, ncirc0, n_p_active, ncandidates=5000, naccepted, nnew, ncirc1;
double dx, dy, p, phi, r, r0, ra[5], sa[5], height, x, y = 0.0, gamma, ymean, ytop, ybottom, dpoisson = PDISC_FACTOR*MU, xmax;
short int active_poisson[NMAXCIRCLES], far;
ymean = 0.5*(YMIN + YMAX);
if (top)
{
ytop = YMAX;
ybottom = ymean;
}
else
{
ytop = ymean;
ybottom = YMIN;
}
switch (pattern) {
case (C_SQUARE):
{
@@ -34,6 +44,41 @@ void init_circle_config_half(int pattern, int top, t_circle circles[NMAXCIRCLES]
ncircles += NGRIDX*NGRIDY/2;
break;
}
case (C_SQUARE_TWO_SPACINGS):
{
dy = (YMAX - YMIN)/((double)NGRIDY);
for (i = 0; i < NGRIDX/2; i++)
for (j = 0; j < NGRIDY/2; j++)
{
printf("i = %i, j = %i, n = %i\n", i, j, n);
n = ncircles;
circles[n].xc = ((double)(i-NGRIDX/2) + 0.5)*dy;
y = ((double)j + 0.5)*dy;
if (top) circles[n].yc = ymean + y;
else circles[n].yc = ymean - y;
if (top) circles[n].radius = MU;
else circles[n].radius = MUB;
circles[n].active = 1;
circletop[n] = top;
ncircles++;
}
for (i = 0; i < NGRIDX/2; i++)
for (j = 0; j < NGRIDY*3/2; j++)
{
printf("i = %i, j = %i, n = %i\n", i, j, n);
n = ncircles;
circles[n].xc = ((double)i + 0.5)*dy;
y = ((double)j + 0.5)*dy/1.5;
if (top) circles[n].yc = ymean + y;
else circles[n].yc = ymean - y;
if (top) circles[n].radius = MU;
else circles[n].radius = MUB;
circles[n].active = 1;
circletop[n] = top;
ncircles++;
}
break;
}
case (C_HEX):
{
dy = (YMAX - YMIN)/((double)NGRIDY);
@@ -43,7 +88,7 @@ void init_circle_config_half(int pattern, int top, t_circle circles[NMAXCIRCLES]
{
n = ncircles + (NGRIDY+2)*i/2 + j;
circles[n].xc = ((double)(i-NGRIDX/2) + 0.5)*dy;
y = ((double)j - 0.5)*dy;
y = ((double)j - 0.5)*dy + 0.25*dy;
if ((i+NGRIDX)%2 == 1) y += 0.5*dy;
if (top) circles[n].yc = ymean + 0.5*dy + y;
else circles[n].yc = ymean - 0.5*dy - y;
@@ -97,7 +142,37 @@ void init_circle_config_half(int pattern, int top, t_circle circles[NMAXCIRCLES]
circletop[n] = top;
printf("n = %i, x = %.3lg\n", n, circles[n].xc);
}
ncircles += NGRIDX*NGRIDY/2;
ncirc0 = ncircles + NGRIDX*NGRIDY/2;
ncirc1 = 0;
/* deal with periodic b.c. */
if ((B_COND == BC_PERIODIC)||(B_COND == BC_VPER_HABS))
{
for (i=0; i<NGRIDX*NGRIDY/2; i++)
{
x = circles[ncircles + i].xc;
y = circles[ncircles + i].yc;
r = circles[ncircles + i].radius;
if (y > ytop - r)
{
circles[ncirc0 + ncirc1].xc = x;
circles[ncirc0 + ncirc1].yc = y - ytop + ybottom;
circles[ncirc0 + ncirc1].radius = r;
circles[ncirc0 + ncirc1].active = 1;
ncirc1++;
}
else if (y < ybottom + r)
{
circles[ncirc0 + ncirc1].xc = x;
circles[ncirc0 + ncirc1].yc = y + ytop - ybottom;
circles[ncirc0 + ncirc1].radius = r;
circles[ncirc0 + ncirc1].active = 1;
ncirc1++;
}
}
}
ncircles = ncirc0 + ncirc1;
break;
}
case (C_RAND_PERCOL):
@@ -145,19 +220,51 @@ void init_circle_config_half(int pattern, int top, t_circle circles[NMAXCIRCLES]
for (i = 0; i < 40; i++)
for (j = 0; j < 5; j++)
{
n = ncircles + 5*i + j;
phi = (double)i*DPI/40.0;
n = ncircles;
phi = (0.5 + (double)i)*DPI/40.0;
r = LAMBDA*0.5*(1.0 + (double)j/5.0);
circles[n].xc = r*cos(phi);
y = r*sin(phi);
if ( ((top)&&(y > 0.0))||((!top)&&(y <= 0.0)) )
{
circles[n].yc = y;
if (top) circles[n].radius = MU;
else circles[n].radius = MUB;
circles[n].active = 1;
circletop[n] = top;
ncircles++;
}
ncircles += 200;
}
break;
}
case (C_CLOAK_OBJECT):
{
for (i = 0; i < 40; i++)
for (j = 0; j < 5; j++)
{
n = ncircles;
phi = (0.5 + (double)i)*DPI/40.0;
r = LAMBDA*0.5*(1.0 + (double)j/5.0);
circles[n].xc = r*cos(phi);
y = r*sin(phi);
if ( ((top)&&(y > 0.0))||((!top)&&(y <= 0.0)) )
{
circles[n].yc = y;
if (top) circles[n].radius = MU;
else circles[n].radius = MUB;
circles[n].active = 1;
circletop[n] = top;
ncircles++;
}
}
n = ncircles;
circles[n].xc = 0.0;
if (top) circles[n].yc = MU;
else circles[n].yc = -MU;
circles[n].radius = 0.25*LAMBDA;
circles[n].active = 1;
circletop[n] = top;
ncircles++;
break;
}
case (C_CLOAK_A): /* optimized model A1 by C. Jo et al */
@@ -197,14 +304,60 @@ void init_circle_config_half(int pattern, int top, t_circle circles[NMAXCIRCLES]
ncircles += 1;
break;
}
case (C_CLOAK_AB): /* optimized model A1 by C. Jo et al */
{
ra[0] = 0.0731; sa[0] = 1.115;
ra[1] = 0.0768; sa[1] = 1.292;
ra[2] = 0.0652; sa[2] = 1.464;
ra[3] = 0.056; sa[3] = 1.633;
ra[4] = 0.0375; sa[4] = 1.794;
for (i = 0; i < 21; i++)
for (j = 0; j < 5; j++)
{
n = ncircles + 5*i + j;
phi = ((double)i + 0.5)*DPI/40.0;
r = 0.5*LAMBDA*sa[j];
circles[n].xc = r*cos(phi);
circles[n].yc = r*sin(phi);
if (top) y = r*sin(phi);
else y = -r*sin(phi);
circles[n].yc = y;
circles[n].radius = 0.3*LAMBDA*ra[j];
circles[n].active = 1;
circletop[n] = top;
}
ncircles += 105;
/* add circle in the center */
// circles[ncircles].xc = 0.0;
// circles[ncircles].yc = 0.0;
// if (top) circles[ncircles].radius = MU;
// else circles[ncircles].radius = MUB;
// circles[ncircles].active = 2;
// circletop[ncircles] = top;
// ncircles += 1;
break;
}
case (C_POISSON_DISC):
{
printf("Generating Poisson disc sample\n");
/* generate first circle */
n = ncircles;
ncirc1 = ncircles;
circles[n].xc = LAMBDA*(2.0*(double)rand()/RAND_MAX - 1.0);
if (NO_CIRCLES_ON_BOUNDARY)
{
if (top) y = ymean + MU + (YMAX-ymean-2.0*MU)*(double)rand()/RAND_MAX;
else y = ymean - MU_B + (YMIN-ymean-2.0*MU_B)*(double)rand()/RAND_MAX;
}
else
{
if (top) y = ymean + (YMAX-ymean)*(double)rand()/RAND_MAX;
else y = ymean + (YMIN-ymean)*(double)rand()/RAND_MAX;
}
circles[n].yc = y;
if (top) circles[n].radius = MU;
else circles[n].radius = MUB;
@@ -235,10 +388,24 @@ void init_circle_config_half(int pattern, int top, t_circle circles[NMAXCIRCLES]
{
/* new circle is far away from circle k */
far = far*((x - circles[k].xc)*(x - circles[k].xc) + (y - circles[k].yc)*(y - circles[k].yc) >= dpoisson*dpoisson);
if ((B_COND == BC_PERIODIC)||(B_COND == BC_VPER_HABS))
{
far *= ((x - circles[k].xc)*(x - circles[k].xc) + (y + ytop - ybottom - circles[k].yc)*(y + ytop - ybottom - circles[k].yc) >= dpoisson*dpoisson);
far *= ((x - circles[k].xc)*(x - circles[k].xc) + (y - ytop + ybottom - circles[k].yc)*(y - ytop + ybottom - circles[k].yc) >= dpoisson*dpoisson);
}
/* new circle is in domain */
if (NO_CIRCLES_ON_BOUNDARY)
{
if (top) far = far*(vabs(x) < LAMBDA)*(y < YMAX - MU)*(y > ymean + MU);
else far = far*(vabs(x) < LAMBDA)*(y > YMIN + MU_B)*(y < ymean - MU_B);
}
else
{
if (top) far = far*(vabs(x) < LAMBDA)*(y < YMAX)*(y > 0.0);
else far = far*(vabs(x) < LAMBDA)*(y > YMIN)*(y < 0.0);
}
}
if (far) /* accept new circle */
{
printf("New circle at (%.3f,%.3f) accepted\n", x, y);
@@ -267,6 +434,32 @@ void init_circle_config_half(int pattern, int top, t_circle circles[NMAXCIRCLES]
// sleep(1);
}
/* deal with periodic b.c. */
if ((B_COND == BC_PERIODIC)||(B_COND == BC_VPER_HABS))
{
for (i=0; i<ncirc0; i++)
{
x = circles[ncirc1 + i].xc;
y = circles[ncirc1 + i].yc;
if (y > ytop - MU)
{
circles[ncircles].xc = x;
circles[ncircles].yc = y - ytop + ybottom;
circles[ncircles].radius = MU;
circles[ncircles].active = 1;
ncircles++;
}
else if (y < ybottom + MU)
{
circles[ncircles].xc = x;
circles[ncircles].yc = y + ytop - ybottom;
circles[ncircles].radius = MU;
circles[ncircles].active = 1;
ncircles++;
}
}
}
printf("Already existing: %i circles\n", ncircles);
ncircles += ncirc0;
printf("Generated %i circles\n", ncirc0);
@@ -415,13 +608,39 @@ void init_circle_config_half(int pattern, int top, t_circle circles[NMAXCIRCLES]
}
void mirror_circle_config(t_circle circles[NMAXCIRCLES])
/* mirror circles in top half to bottom half */
{
int i, ncircles0 = ncircles;
printf("Mirroring circles, ncircles = %i\n", ncircles);
for (i=0; i<ncircles0; i++)
{
// printf("Mirroring circle %i to %i\n", i, ncircles);
circles[ncircles].xc = circles[i].xc;
circles[ncircles].yc = -circles[i].yc;
circles[ncircles].radius = circles[i].radius;
circles[ncircles].active = circles[i].active;
circles[ncircles].top = 0;
// printf("Circle[%i].top = %i\n", ncircles, circles[ncircles].top);
ncircles++;
}
}
void init_circle_config_comp(t_circle circles[NMAXCIRCLES])
/* initialise the arrays circlex, circley, circlerad and circleactive */
/* for billiard shape D_CIRCLES */
{
ncircles = 0;
init_circle_config_half(CIRCLE_PATTERN, 1, circles);
init_circle_config_half(CIRCLE_PATTERN_B, 0, circles);
ncircles++;
if (SYMMETRIC_CIRCLE_PATTERNS)
mirror_circle_config(circles);
else init_circle_config_half(CIRCLE_PATTERN_B, 0, circles);
circles[126].top = 0;
}
void init_circle_config_energy(t_circle circles[NMAXCIRCLES])
@@ -531,6 +750,23 @@ int xy_in_billiard_half(double x, double y, int domain, int pattern, int top)
}
return(1);
}
case (D_CIRCLES_NEUMANN):
{
for (i = 0; i < ncircles; i++)
if (circles[i].active != 0)
{
/* choose specific type according to value of circles[i].active */
if (circles[i].active == 1) type = 0;
else type = circles[i].active;
x1 = circles[i].xc;
y1 = circles[i].yc;
r2 = circles[i].radius*circles[i].radius;
if ((top)&&(circletop[i])&&(y > 0.0)&&((x-x1)*(x-x1) + (y-y1)*(y-y1) < r2)) return(type);
else if ((!top)&&(!circletop[i])&&(y < 0.0)&&((x-x1)*(x-x1) + (y-y1)*(y-y1) < r2)) return(type);
}
return(1);
}
case (D_POLYGONS):
{
for (i = 0; i < ncircles; i++)
@@ -769,10 +1005,10 @@ void draw_billiard_half(int domain, int pattern, int top, int fade, double fade_
glEnable(GL_SCISSOR_TEST);
// if (top) glScissor(0.0, YMID, NX, YMID);
// if (top) glScissor(0.0, YMID, NX, YMAX);
// if (top) glScissor(0.0, YMID, NX, YMID);
// else glScissor(0.0, 0.0, NX, YMID);
if (top) glScissor(0.0, YMID/(HIGHRES+1), NX, YMID/(HIGHRES+1));
if (top) glScissor(0.0, YMID/(HIGHRES+1.0), NX, YMID/(HIGHRES+1.0));
// if (top) glScissor(0.0, YMID, NX, YMAX);
else glScissor(0.0, 0.0, NX, YMID/(HIGHRES+1));
@@ -925,6 +1161,13 @@ void draw_billiard_half(int domain, int pattern, int top, int fade, double fade_
}
break;
}
case (D_CIRCLES_NEUMANN):
{
glLineWidth(BOUNDARY_WIDTH);
for (i = 0; i < ncircles; i++)
if (circles[i].active) draw_circle(circles[i].xc, circles[i].yc, circles[i].radius, NSEG);
break;
}
case (D_POLYGONS):
{
glLineWidth(BOUNDARY_WIDTH);
@@ -1413,7 +1656,7 @@ void draw_wave_comp(double *phi[NX], double *psi[NX], short int *xy_in[NX], doub
}
void draw_wave_comp_highres_palette(int size, double *phi[NX], double *psi[NX], double *total_energy[NX], double *average_energy[NX], double *total_flux, double *color_scale[NX], short int *xy_in[NX], double scale, int time, int plot, int palette, int fade, double fade_value)
void draw_wave_comp_highres_palette(int size, double *phi[NX], double *psi[NX], double *total_energy[NX], double *average_energy[NX], double *total_flux, double *color_scale[NX], short int *xy_in[NX], double scale, int time, int plot, int palette, int pnumber, int fade, double fade_value)
/* draw the field */
{
int i, j, iplus, iminus, jplus, jminus, k;
@@ -1500,6 +1743,8 @@ void draw_wave_comp_highres_palette(int size, double *phi[NX], double *psi[NX],
glVertex2d(pos[0], pos[1]);
glEnd();
if (DRAW_WAVE_PROFILE) draw_wave_profile(values, size, pnumber, fade, fade_value);
free(values);
free(rgbvals);
}
wave_3d.c
+187 -81
View File
@@ -44,7 +44,7 @@
#include <time.h>
#define MOVIE 0 /* set to 1 to generate movie */
#define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */
#define DOUBLE_MOVIE 0 /* set to 1 to produce movies for wave height and energy simultaneously */
#define SAVE_MEMORY 1 /* set to 1 to save memory when writing tiff images */
#define NO_EXTRA_BUFFER_SWAP 1 /* some OS require one less buffer swap when recording images */
@@ -52,15 +52,15 @@
#define WINWIDTH 1920 /* window width */
#define WINHEIGHT 1150 /* window height */
// #define NX 1920 /* number of grid points on x axis */
// #define NY 1150 /* number of grid points on y axis */
// #define NX 3840 /* number of grid points on x axis */
// #define NY 2300 /* number of grid points on y axis */
#define NX 2880 /* number of grid points on x axis */
#define NY 1725 /* number of grid points on y axis */
// #define NX 2400 /* number of grid points on x axis */
// #define NY 1400 /* number of grid points on y axis */
#define XMIN -2.0
#define XMAX 2.0 /* x interval */
// #define XMIN -2.0
// #define XMAX 2.0 /* x interval */
#define XMIN -0.75
#define XMAX 3.25 /* x interval */
#define YMIN -1.197916667
#define YMAX 1.197916667 /* y interval for 9/16 aspect ratio */
@@ -70,17 +70,17 @@
/* Choice of the billiard table */
#define B_DOMAIN 192 /* choice of domain shape, see list in global_pdes.c */
#define B_DOMAIN 202 /* choice of domain shape, see list in global_pdes.c */
#define CIRCLE_PATTERN 2 /* pattern of circles or polygons, see list in global_pdes.c */
#define CIRCLE_PATTERN 5 /* pattern of circles or polygons, see list in global_pdes.c */
#define COMPARISON 0 /* set to 1 to compare two different patterns */
#define B_DOMAIN_B 20 /* second domain shape, for comparisons */
#define CIRCLE_PATTERN_B 0 /* second pattern of circles or polygons */
#define COMPARISON 1 /* set to 1 to compare two different patterns */
#define B_DOMAIN_B 202 /* second domain shape, for comparisons */
#define CIRCLE_PATTERN_B 5 /* second pattern of circles or polygons */
#define IMAGE_FILE 5 /* for option D_IMAGE */
#define VARIABLE_IOR 0 /* set to 1 for a variable index of refraction */
#define IOR 182 /* choice of index of refraction, see list in global_pdes.c */
#define IOR 17 /* choice of index of refraction, see list in global_pdes.c */
#define IOR_TOTAL_TURNS 1.0 /* total angle of rotation for IOR_PERIODIC_WELLS_ROTATING */
#define MANDEL_IOR_SCALE -0.05 /* parameter controlling dependence of IoR on Mandelbrot escape speed */
@@ -90,20 +90,20 @@
#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 2.3 /* parameter controlling the dimensions of domain */
#define MU 1.15 /* parameter controlling the dimensions of domain */
#define MU_B 1.0 /* parameter controlling the dimensions of domain */
#define NPOLY 3 /* number of sides of polygon */
#define APOLY -1.0 /* angle by which to turn polygon, in units of Pi/2 */
#define MDEPTH 6 /* depth of computation of Menger gasket */
#define LAMBDA 0.4 /* parameter controlling the dimensions of domain */
#define MU 0.008 /* parameter controlling the dimensions of domain */
#define MU_B 0.08 /* parameter controlling the dimensions of domain */
#define NPOLY 5 /* number of sides of polygon */
#define APOLY 0.2 /* angle by which to turn polygon, in units of Pi/2 */
#define MDEPTH 5 /* depth of computation of Menger gasket */
#define MRATIO 3 /* ratio defining Menger gasket */
#define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */
#define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */
#define FOCI 1 /* set to 1 to draw focal points of ellipse */
#define NGRIDX 12 /* number of grid point for grid of disks */
#define NGRIDX 4 /* number of grid point for grid of disks */
#define NGRIDY 8 /* 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) */
@@ -128,19 +128,20 @@
#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.1 /* defines oscilling beam range */
#define OSCIL_YMID -0.75 /* defines oscilling beam midpoint */
#define OSCILLATE_TOPBOT 1 /* set to 1 to enforce a planar wave on top and bottom boundary */
#define OSCILLATION_SCHEDULE 0 /* oscillation schedule, see list in global_pdes.c */
#define 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.02 /* frequency of periodic excitation */
#define OMEGA 0.007 /* frequency of periodic excitation */
#define OMEGA_B 0.01 /* frequency of periodic excitation */
#define AMPLITUDE 1.0 /* amplitude of periodic excitation */
#define ACHIRP 0.2 /* acceleration coefficient in chirp */
#define DAMPING 0.0 /* damping of periodic excitation */
#define COURANT 0.15 /* Courant number in medium B */
#define COURANTB 0.1 /* Courant number */
#define COURANT 0.09 /* Courant number in medium B */
#define COURANTB 0.2 /* Courant number */
#define GAMMA 0.0 /* damping factor in wave equation */
#define GAMMAB 1.0e-6 /* damping factor in wave equation */
#define GAMMA_SIDES 1.0e-4 /* damping factor on boundary */
@@ -155,9 +156,9 @@
/* 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 18 /* period of oscillating source */
#define OSCILLATING_SOURCE_PERIOD 12 /* period of oscillating source */
#define ALTERNATE_OSCILLATING_SOURCE 1 /* set to 1 to alternate sign of oscillating source */
#define MAX_PULSING_TIME 1500 /* max time for adding pulses */
#define MAX_PULSING_TIME 23 /* 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 1 /* type of wave packet sources */
@@ -166,16 +167,17 @@
/* Boundary conditions, see list in global_pdes.c */
#define B_COND 2
#define B_COND 4
#define BC_NEUMANN 1 /* set to 1 to use Neumann boundary conditions on domain */
#define PRECOMPUTE_BC 0 /* set to 1 to compute neighbours for Laplacian in advance */
/* Parameters for length and speed of simulation */
#define NSTEPS 2000 /* number of frames of movie */
#define NSTEPS 2600 /* number of frames of movie */
#define NVID 10 /* number of iterations between images displayed on screen */
#define NSEG 200 /* number of segments of boundary */
#define INITIAL_TIME 100 /* time after which to start saving frames */
#define INITIAL_TIME 0 /* time after which to start saving frames */
#define BOUNDARY_WIDTH 2 /* width of billiard boundary */
#define PRINT_SPEED 0 /* set to 1 to print speed of moving source */
@@ -190,9 +192,9 @@
/* Parameters of initial condition */
#define INITIAL_AMP 0.75 /* amplitude of initial condition */
#define INITIAL_VARIANCE 0.0005 /* variance of initial condition */
#define INITIAL_WAVELENGTH 0.025 /* wavelength of initial condition */
#define INITIAL_AMP 1.0 /* amplitude of initial condition */
#define INITIAL_VARIANCE 0.001 /* variance of initial condition */
#define INITIAL_WAVELENGTH 0.007 /* wavelength of initial condition */
/* Plot type, see list in global_pdes.c */
@@ -206,12 +208,12 @@
#define FLUX_WINDOW 30 /* size of averaging window of flux intensity */
#define AMPLITUDE_HIGH_RES 1 /* set to 1 to increase resolution of plot */
#define SHADE_3D 1 /* set to 1 to change luminosity according to normal vector */
#define NON_DIRICHLET_BC 0 /* set to 1 to draw only facets in domain, if field is not zero on boundary */
#define NON_DIRICHLET_BC 1 /* set to 1 to draw only facets in domain, if field is not zero on boundary */
#define FLOOR_ZCOORD 1 /* set to 1 to draw only facets with z not too negative */
#define DRAW_BILLIARD 0 /* set to 1 to draw boundary */
#define DRAW_BILLIARD_FRONT 0 /* set to 1 to draw front of boundary after drawing wave */
#define DRAW_CONSTRUCTION_LINES 0 /* set to 1 to draw construction lines of certain domains */
#define FADE_IN_OBSTACLE 0 /* set to 1 to fade color inside obstacles */
#define FADE_IN_OBSTACLE 1 /* set to 1 to fade color inside obstacles */
#define DRAW_OUTSIDE_GRAY 0 /* experimental, draw outside of billiard in gray */
#define PLOT_SCALE_ENERGY 0.4 /* vertical scaling in energy plot */
@@ -225,11 +227,11 @@
#define REP_PROJ_3D 1 /* projection on plane orthogonal to observer line of sight */
#define ROTATE_VIEW 0 /* set to 1 to rotate position of observer */
#define ROTATE_ANGLE 360.0 /* total angle of rotation during simulation */
#define ROTATE_ANGLE -360.0 /* total angle of rotation during simulation */
/* Color schemes */
#define COLOR_PALETTE 11 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE 17 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 13 /* Color palette, see list in global_pdes.c */
#define BLACK 1 /* background */
@@ -240,12 +242,12 @@
#define SLOPE 1.0 /* sensitivity of color on wave amplitude */
#define COLOR_RANGE 1.0 /* max range of color (default: 1.0) */
#define VSCALE_AMPLITUDE 1.0 /* additional scaling factor for color scheme P_3D_AMPLITUDE */
#define VSHIFT_AMPLITUDE 0.0 /* additional shift for wave amplitude */
#define VSCALE_ENERGY 2.5 /* additional scaling factor for color scheme P_3D_ENERGY */
#define VSHIFT_AMPLITUDE -0.0 /* additional shift for wave amplitude */
#define VSCALE_ENERGY 0.75 /* additional scaling factor for color scheme P_3D_ENERGY */
#define PHASE_FACTOR 20.0 /* factor in computation of phase in color scheme P_3D_PHASE */
#define PHASE_SHIFT 0.0 /* shift of phase in color scheme P_3D_PHASE */
#define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */
#define E_SCALE 100.0 /* scaling factor for energy representation */
#define E_SCALE 150.0 /* scaling factor for energy representation */
#define LOG_SCALE 1.5 /* scaling factor for energy log representation */
#define LOG_SHIFT 0.25 /* shift of colors on log scale */
#define LOG_ENERGY_FLOOR -10.0 /* floor value for log of (total) energy */
@@ -276,9 +278,9 @@
#define MAZE_XSHIFT 0.0 /* horizontal shift of maze */
#define MAZE_WIDTH 0.02 /* half width of maze walls */
#define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */
#define COLORBAR_RANGE 3.0 /* scale of color scheme bar */
#define COLORBAR_RANGE_B 1.5 /* scale of color scheme bar for 2nd part */
#define DRAW_COLOR_SCHEME 0 /* set to 1 to plot the color scheme */
#define COLORBAR_RANGE 2.0 /* scale of color scheme bar */
#define COLORBAR_RANGE_B 8.0 /* scale of color scheme bar for 2nd part */
#define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */
#define SAVE_TIME_SERIES 0 /* set to 1 to save wave time series at a point */
@@ -287,7 +289,6 @@
#define POTENTIAL 10
#define POT_FACT 200.0
#define DRAW_WAVE_PROFILE 0 /* set to 1 to draw a profile of the wave */
#define MU_B 1.0 /* parameter controlling the dimensions of domain */
#define HORIZONTAL_WAVE_PROFILE 0 /* set to 1 to draw wave profile vertically */
#define VERTICAL_WAVE_PROFILE 0 /* set to 1 to draw wave profile vertically */
#define DRAW_WAVE_TIMESERIES 0 /* set to 1 to draw a time series of the wave */
@@ -300,7 +301,8 @@
#define HRES 1 /* dummy, only used by rde.c */
#define SHADE_2D 0 /* set to 1 to add pseudo-3d shading effect */
#define SHADE_SCALE_2D 0.05 /* lower value increases sensitivity of shading */
#define N_SOURCES 1 /* number of sources, for option draw_sources */
#define N_SOURCES 4 /* number of sources, for option draw_sources */
#define ALTERNATE_SOURCE_PHASES 0 /* set to 1 to alternate initial phases of sources */
#define XYIN_INITIALISED (B_DOMAIN == D_IMAGE)
/* end of constants only used by sub_wave and sub_maze */
@@ -318,11 +320,12 @@ double light[3] = {0.816496581, -0.40824829, 0.40824829}; /* vector of "lig
double observer[3] = {8.0, 6.0, 6.0}; /* location of observer for REP_PROJ_3D representation */
int reset_view = 0; /* switch to reset 3D view parameters (for option ROTATE_VIEW) */
#define Z_SCALING_FACTOR 0.1 /* overall scaling factor of z axis for REP_PROJ_3D representation */
#define XY_SCALING_FACTOR 2.3 /* overall scaling factor for on-screen (x,y) coordinates after projection */
#define Z_SCALING_FACTOR 0.07 /* overall scaling factor of z axis for REP_PROJ_3D representation */
#define XY_SCALING_FACTOR 1.8 /* overall scaling factor for on-screen (x,y) coordinates after projection */
#define 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.25 /* overall y shift for REP_PROJ_3D representation */
#define XSHIFT_3D 0.2 /* overall x shift for REP_PROJ_3D representation */
#define YSHIFT_3D 0.6 /* overall y shift for REP_PROJ_3D representation */
#define BEAM_BC ((OSCILLATION_SCHEDULE == OSC_BEAM_SINE)||(OSCILLATION_SCHEDULE == OSC_BEAM_SINE_DECREASING)||(OSCILLATION_SCHEDULE == OSC_BEAM_SINE_TWOPERIODS)||(OSCILLATION_SCHEDULE == OSC_BEAM_SINE_CHIRP))
@@ -433,6 +436,11 @@ void evolve_wave_half_new(double phi_in[NX*NY], double psi_in[NX*NY], double phi
phi_out[j] = x - tc[j]*(x - phi_in[NY+j]) - KAPPA_SIDES*x - GAMMA_SIDES*(x-y);
break;
}
case (BC_ABS_REFLECT):
{
phi_out[j] = x - tc[j]*(x - phi_in[NY+j]) - KAPPA_SIDES*x - GAMMA_SIDES*(x-y);
break;
}
}
}
}
@@ -464,6 +472,11 @@ void evolve_wave_half_new(double phi_in[NX*NY], double psi_in[NX*NY], double phi
phi_out[(NX-1)*NY+j] = x - tc[(NX-1)*NY+j]*(x - phi_in[(NX-2)*NY+j]) - KAPPA_SIDES*x - GAMMA_SIDES*(x-y);
break;
}
case (BC_ABS_REFLECT):
{
phi_out[(NX-1)*NY+j] = x - tc[(NX-1)*NY+j]*(x - phi_in[(NX-2)*NY+j]) - KAPPA_SIDES*x - GAMMA_SIDES*(x-y);
break;
}
}
}
}
@@ -504,6 +517,12 @@ void evolve_wave_half_new(double phi_in[NX*NY], double psi_in[NX*NY], double phi
else phi_out[i*NY+NY-1] = -y + 2*x + tcc[i*NY+NY-1]*delta[i*NY+NY-1] - KAPPA*x - tgamma[i*NY+NY-1]*(x-y);
break;
}
case (BC_ABS_REFLECT):
{
phi_out[i*NY+NY-1] = x - tc[i*NY+NY-1]*(x - phi_in[i*NY+NY-2]) - KAPPA_TOPBOT*x - GAMMA_TOPBOT*(x-y);
break;
}
}
}
}
@@ -542,6 +561,12 @@ void evolve_wave_half_new(double phi_in[NX*NY], double psi_in[NX*NY], double phi
else phi_out[i*NY] = -y + 2*x + tcc[i*NY]*delta[i*NY] - KAPPA*x - tgamma[i*NY]*(x-y);
break;
}
case (BC_ABS_REFLECT):
{
phi_out[i*NY] = x - tc[i*NY]*(x - phi_in[i*NY+1]) - KAPPA_TOPBOT*x - GAMMA_TOPBOT*(x-y);
break;
}
}
}
}
@@ -572,13 +597,12 @@ void evolve_wave_half_new(double phi_in[NX*NY], double psi_in[NX*NY], double phi
}
void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out[NX*NY],
short int xy_in[NX*NY], double tc[NX*NY], double tcc[NX*NY], double tgamma[NX*NY])
void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out[NX*NY], short int xy_in[NX*NY], double tc[NX*NY], double tcc[NX*NY], double tgamma[NX*NY], t_neumann_bc neumann_bc[N_NEUMANN_POINTS])
/* time step of field evolution */
/* phi is value of field at time t, psi at time t-1 */
/* this version of the function has been rewritten in order to minimize the number of if-branches */
{
int i, j, iplus, iminus, jplus, jminus, jtop, jbot;
int i, j, k, iplus, iminus, jplus, jminus, jtop, jbot, nc = 0;
double delta, x, y, c, cc, gamma;
static long time = 0;
static short int first = 1;
@@ -597,7 +621,7 @@ void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out
/* evolution in the bulk */
for (i=1; i<NX-1; i++){
for (j=1; j<NY-1; j++){
if ((TWOSPEEDS)||(xy_in[i*NY+j] != 0)){
if ((TWOSPEEDS)||(xy_in[i*NY+j] > 0)){
x = phi_in[i*NY+j];
y = psi_in[i*NY+j];
@@ -610,6 +634,18 @@ void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out
}
}
if (BC_NEUMANN)
{
#pragma omp parallel for private(nc,k)
for (nc= 0; nc < n_neumann_points; nc++) if (neumann_bc[nc].active)
{
k = neumann_bc[nc].k;
phi_out[k] = neumann_bc[nc].nx*phi_in[neumann_bc[nc].nbk1];
phi_out[k] += neumann_bc[nc].ny*phi_in[neumann_bc[nc].nbk2];
}
}
/* left boundary */
// if (OSCILLATE_LEFT) for (j=1; j<NY-1; j++) phi_out[j] = AMPLITUDE*cos((double)time*OMEGA);
if (OSCILLATE_LEFT)
@@ -647,6 +683,12 @@ void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out
phi_out[j] = x - tc[j]*(x - phi_in[NY+j]) - KAPPA_SIDES*x - GAMMA_SIDES*(x-y);
break;
}
case (BC_ABS_REFLECT):
{
delta = phi_in[NY+j] + phi_in[j+1] + phi_in[j-1] - 3.0*x;
phi_out[j] = x - tc[j]*(x - phi_in[NY+j]) - KAPPA_SIDES*x - GAMMA_SIDES*(x-y);
break;
}
}
}
}
@@ -682,6 +724,12 @@ void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out
phi_out[(NX-1)*NY+j] = x - tc[(NX-1)*NY+j]*(x - phi_in[(NX-2)*NY+j]) - KAPPA_SIDES*x - GAMMA_SIDES*(x-y);
break;
}
case (BC_ABS_REFLECT):
{
delta = phi_in[(NX-2)*NY+j] + phi_in[(NX-1)*NY+j+1] + phi_in[(NX-1)*NY+j-1] - 3.0*x;
phi_out[(NX-1)*NY+j] = x - tc[(NX-1)*NY+j]*(x - phi_in[(NX-2)*NY+j]) - KAPPA_SIDES*x - GAMMA_SIDES*(x-y);
break;
}
}
}
}
@@ -741,6 +789,15 @@ void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out
else phi_out[i*NY+NY-1] = -y + 2*x + tcc[i*NY+NY-1]*delta - KAPPA*x - tgamma[i*NY+NY-1]*(x-y);
break;
}
case (BC_ABS_REFLECT):
{
iplus = (i+1); if (iplus == NX) iplus = NX-1;
iminus = (i-1); if (iminus == -1) iminus = 0;
delta = phi_in[iplus*NY+NY-1] + phi_in[iminus*NY+NY-1] + phi_in[i*NY+NY-2] - 3.0*x;
phi_out[i*NY+NY-1] = x - tc[i*NY+NY-1]*(x - phi_in[i*NY+NY-2]) - KAPPA_TOPBOT*x - GAMMA_TOPBOT*(x-y);
break;
}
}
}
}
@@ -800,6 +857,15 @@ void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out
else phi_out[i*NY] = -y + 2*x + tcc[i*NY]*delta - KAPPA*x - tgamma[i*NY]*(x-y);
break;
}
case (BC_ABS_REFLECT):
{
iplus = (i+1); if (iplus == NX) iplus = NX-1;
iminus = (i-1); if (iminus == -1) iminus = 0;
delta = phi_in[iplus*NY] + phi_in[iminus*NY] + phi_in[i*NY+1] - 3.0*x;
phi_out[i*NY] = x - tc[i*NY]*(x - phi_in[i*NY+1]) - KAPPA_TOPBOT*x - GAMMA_TOPBOT*(x-y);
break;
}
}
}
}
@@ -849,6 +915,15 @@ void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out
else phi_out[i*NY+jtop] = -y + 2*x + tcc[i*NY+jtop]*delta - KAPPA*x - tgamma[i*NY+jtop]*(x-y);
break;
}
case (BC_ABS_REFLECT):
{
iplus = i+1; if (iplus == NX) iplus = NX-1;
iminus = i-1; if (iminus == -1) iminus = 0;
delta = phi_in[iplus*NY+jtop] + phi_in[iminus*NY+jtop] + phi_in[i*NY+jtop-1] - 3.0*x;
phi_out[i*NY+jtop] = -y + 2*x + tcc[i*NY+jtop]*delta - KAPPA*x - tgamma[i*NY+jtop]*(x-y);
break;
}
}
}
}
@@ -898,6 +973,14 @@ void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out
else phi_out[i*NY+jbot] = -y + 2*x + tcc[i*NY+jbot]*delta - KAPPA*x - tgamma[i*NY+jbot]*(x-y);
break;
}
case (BC_ABS_REFLECT):
{
iplus = (i+1); if (iplus == NX) iplus = NX-1;
iminus = (i-1); if (iminus == -1) iminus = 0;
delta = phi_in[iplus*NY+jbot] + phi_in[iminus*NY+jbot] + phi_in[i*NY+jbot+1] - 3.0*x;
phi_out[i*NY+jbot] = -y + 2*x + tcc[i*NY+jbot]*delta - KAPPA*x - tgamma[i*NY+jbot]*(x-y);
}
}
}
}
@@ -925,8 +1008,7 @@ void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out
void evolve_wave(double phi[NX*NY], double psi[NX*NY], double tmp[NX*NY], short int xy_in[NX*NY],
double tc[NX*NY], double tcc[NX*NY], double tgamma[NX*NY], t_laplacian laplace[NX*NY], t_laplacian laplace1[NX*NY],
t_laplacian laplace2[NX*NY])
double tc[NX*NY], double tcc[NX*NY], double tgamma[NX*NY], t_laplacian laplace[NX*NY], t_laplacian laplace1[NX*NY], t_laplacian laplace2[NX*NY], t_neumann_bc neumann_bc[N_NEUMANN_POINTS])
/* time step of field evolution */
/* phi is value of field at time t, psi at time t-1 */
{
@@ -942,9 +1024,9 @@ void evolve_wave(double phi[NX*NY], double psi[NX*NY], double tmp[NX*NY], short
{
// evolve_wave_half(phi, psi, phi_tmp, psi_tmp, xy_in, tc, tcc, tgamma);
// evolve_wave_half(phi_tmp, psi_tmp, phi, psi, xy_in, tc, tcc, tgamma);
evolve_wave_half(phi, psi, tmp, xy_in, tc, tcc, tgamma);
evolve_wave_half(tmp, phi, psi, xy_in, tc, tcc, tgamma);
evolve_wave_half(psi, tmp, phi, xy_in, tc, tcc, tgamma);
evolve_wave_half(phi, psi, tmp, xy_in, tc, tcc, tgamma, neumann_bc);
evolve_wave_half(tmp, phi, psi, xy_in, tc, tcc, tgamma, neumann_bc);
evolve_wave_half(psi, tmp, phi, xy_in, tc, tcc, tgamma, neumann_bc);
}
}
@@ -988,15 +1070,16 @@ void viewpoint_schedule(int i)
void animation()
{
double time, scale, ratio, startleft[2], startright[2], sign = 1.0, r2, xy[2], fade_value, yshift, speed = 0.0, a, b, c, angle = 0.0, lambda1, y, x1, sign1, omega, phase_shift, source_amp[N_SOURCES];
double time, scale, ratio, startleft[2], startright[2], sign[N_SOURCES], r2, xy[2], fade_value, yshift, speed = 0.0, a, b, c, angle = 0.0, lambda1, y, x1, sign1, omega, phase_shift, source_amp[N_SOURCES], source_periods[N_SOURCES], dperiod;
double *phi, *psi, *tmp, *color_scale, *tc, *tcc, *tgamma;
// double *total_energy;
short int *xy_in;
int i, j, s, sample_left[2], sample_right[2], period = 0, fade, source_counter = 0, k, p, q, source, source_periods[N_SOURCES];
int i, j, s, sample_left[2], sample_right[2], period = 0, fade, source_counter = 0, k, p, q, source, phase, add_counter[N_SOURCES];
static int counter = 0, first_source = 1;
long int wave_value;
t_wave *wave;
t_laplacian *laplace, *laplace1, *laplace2;
t_neumann_bc *neumann_bc;
t_wave_source wave_source[25];
if (SAVE_TIME_SERIES)
@@ -1025,24 +1108,33 @@ void animation()
/* initialise positions and radii of circles */
if (COMPARISON)
{
if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN == D_CIRCLES_IN_RECT))
if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN == D_CIRCLES_IN_RECT)||(B_DOMAIN == D_CIRCLES_NEUMANN))
ncircles = init_circle_config_pattern(circles, CIRCLE_PATTERN);
else if (B_DOMAIN == D_POLYGONS) ncircles = init_polygon_config_pattern(polygons, CIRCLE_PATTERN);
if ((B_DOMAIN_B == D_CIRCLES)||(B_DOMAIN_B == D_CIRCLES_IN_RECT))
if ((B_DOMAIN_B == D_CIRCLES)||(B_DOMAIN_B == D_CIRCLES_IN_RECT)||(B_DOMAIN_B == D_CIRCLES_NEUMANN))
ncircles_b = init_circle_config_pattern(circles_b, CIRCLE_PATTERN_B);
else if (B_DOMAIN_B == D_POLYGONS) ncircles_b = init_polygon_config_pattern(polygons_b, CIRCLE_PATTERN_B);
/* TO DO: adapt to different polygon patterns */
}
else
{
if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN == D_CIRCLES_IN_RECT)) ncircles = init_circle_config(circles);
if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN == D_CIRCLES_IN_RECT)||(B_DOMAIN == D_CIRCLES_NEUMANN))
ncircles = init_circle_config(circles);
else if (B_DOMAIN == D_POLYGONS) ncircles = init_polygon_config(polygons);
}
printf("Polygons initialized\n");
if (XYIN_INITIALISED) init_xyin_from_image_3d(xy_in);
if (BC_NEUMANN)
neumann_bc = (t_neumann_bc *)malloc(N_NEUMANN_POINTS*sizeof(t_neumann_bc));
for (i=0; i<N_SOURCES; i++) add_counter[i] = 0;
if (ALTERNATE_SOURCE_PHASES) for (i=0; i<N_SOURCES; i++) sign[i] = pow(-1.0,(double)i);
else for (i=0; i<N_SOURCES; i++) sign[i] = 1.0;
/* initialise polyline for von Koch and similar domains */
// npolyline = init_polyline(MDEPTH, polyline);
//
@@ -1118,8 +1210,7 @@ void animation()
// add_circular_wave_mod(1.0, lambda1*cos(((double)j+0.5)*angle), lambda1*sin(((double)j+0.5)*angle), phi, psi, xy_in);
init_wave_flat_mod(phi, psi, xy_in);
// init_circular_wave_mod(-0.5, 0.0, phi, psi, xy_in);
// add_circular_wave_mod(1.0, 1.0, 0.0, phi, psi, xy_in);
// printf("Wave initialized\n");
@@ -1127,6 +1218,10 @@ void animation()
/* initialize table of wave speeds/dissipation */
init_speed_dissipation(xy_in, tc, tcc, tgamma);
/* initialize Neumann boundary data */
if (BC_NEUMANN)
n_neumann_points = init_neumann_bc_flat(xy_in, neumann_bc, circles, COMPARISON);
/* initialze potential to add to z coordinate */
if (ADD_POTENTIAL)
{
@@ -1181,7 +1276,7 @@ void animation()
draw_wave_3d(0, phi, psi, xy_in, wave, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1);
for (j=0; j<NVID; j++)
{
evolve_wave(phi, psi, tmp, xy_in, tc, tcc, tgamma, laplace, laplace1, laplace2);
evolve_wave(phi, psi, tmp, xy_in, tc, tcc, tgamma, laplace, laplace1, laplace2, neumann_bc);
if (SAVE_TIME_SERIES)
{
wave_value = (long int)(phi[sample_left[0]*NY+sample_left[1]]*1.0e16);
@@ -1199,16 +1294,24 @@ void animation()
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, fade, fade_value);
/* add oscillating waves */
wave_source_x[0] = 0.0;
wave_source_y[0] = 0.0;
source_periods[0] = OSCILLATING_SOURCE_PERIOD;
source_amp[0] = INITIAL_AMP;
for (source = 0; source < N_SOURCES; source++)
{
if ((ADD_OSCILLATING_SOURCE)&&(i%(OSCILLATING_SOURCE_PERIOD) == 1)&&(i<MAX_PULSING_TIME))
wave_source_x[source] = (2.0*(double)source - 1.0)*LAMBDA;
wave_source_y[source] = 0.0;
source_amp[source] = INITIAL_AMP;
}
source_periods[0] = OSCILLATING_SOURCE_PERIOD;
source_periods[1] = OSCILLATING_SOURCE_PERIOD;
for (source = 0; source < N_SOURCES; source++)
{
if (ALTERNATE_OSCILLATING_SOURCE) sign = -sign;
add_circular_wave_mod(-sign*INITIAL_AMP, wave_source_x[source], wave_source_y[source], phi, psi, xy_in);
dperiod = source_periods[source];
phase = i - (int)(dperiod*(double)((int)((double)i/dperiod)));
if ((ADD_OSCILLATING_SOURCE)&&(phase == 1)&&(i<MAX_PULSING_TIME))
{
printf("Source %i: Adding pulse %i at (%.3lg, %.3lg)\n", source, add_counter[source], wave_source_x[source], wave_source_y[source]);
add_counter[source]++;
if (ALTERNATE_OSCILLATING_SOURCE) sign[source] = -sign[source];
add_circular_wave_mod(sign[source]*INITIAL_AMP, wave_source_x[source], wave_source_y[source], phi, psi, xy_in);
}
}
@@ -1270,7 +1373,7 @@ void animation()
{
if ((ROTATE_VIEW)&&(ROTATE_VIEW_WHILE_FADE))
{
viewpoint_schedule(NSTEPS - INITIAL_TIME + i);
viewpoint_schedule(NSTEPS + i);
reset_view = 1;
}
fade_value = 1.0 - (double)i/(double)MID_FRAMES;
@@ -1283,7 +1386,7 @@ void animation()
if ((ROTATE_VIEW)&&(ROTATE_VIEW_WHILE_FADE))
{
viewpoint_schedule(NSTEPS - INITIAL_TIME);
viewpoint_schedule(NSTEPS);
reset_view = 1;
}
draw_wave_3d(1, phi, psi, xy_in, wave, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 0, 1.0, 1);
@@ -1296,7 +1399,8 @@ void animation()
{
if ((ROTATE_VIEW)&&(ROTATE_VIEW_WHILE_FADE))
{
viewpoint_schedule(NSTEPS - INITIAL_TIME + i);
// viewpoint_schedule(NSTEPS + MID_FRAMES + i);
viewpoint_schedule(NSTEPS + i);
reset_view = 1;
}
fade_value = 1.0 - (double)i/(double)END_FRAMES;
@@ -1314,7 +1418,7 @@ void animation()
{
if ((ROTATE_VIEW)&&(ROTATE_VIEW_WHILE_FADE))
{
viewpoint_schedule(NSTEPS - INITIAL_TIME + i);
viewpoint_schedule(NSTEPS + i);
reset_view = 1;
}
fade_value = 1.0 - (double)i/(double)END_FRAMES;
@@ -1344,6 +1448,8 @@ void animation()
free(laplace1);
free(laplace2);
if (BC_NEUMANN) free(neumann_bc);
if (SAVE_TIME_SERIES)
{
wave_billiard.c
+93 -59
View File
@@ -60,10 +60,12 @@
#define NX 3840 /* number of grid points on x axis */
#define NY 2300 /* number of grid points on y axis */
#define XMIN -2.0
#define XMAX 2.0 /* x interval */
#define YMIN -1.297916667
#define YMAX 1.097916667 /* y interval for 9/16 aspect ratio */
// #define XMIN -2.0
// #define XMAX 2.0 /* x interval */
#define XMIN -0.75
#define XMAX 3.25 /* x interval */
#define YMIN -1.197916667
#define YMAX 1.197916667 /* y interval for 9/16 aspect ratio */
#define HIGHRES 1 /* set to 1 if resolution of grid is double that of displayed image */
#define HRES 1 /* dummy, only used by rde.c */
@@ -72,9 +74,9 @@
/* Choice of the billiard table */
#define B_DOMAIN 991 /* choice of domain shape, see list in global_pdes.c */
#define B_DOMAIN 20 /* choice of domain shape, see list in global_pdes.c */
#define CIRCLE_PATTERN 202 /* pattern of circles or polygons, see list in global_pdes.c */
#define CIRCLE_PATTERN 61 /* pattern of circles or polygons, see list in global_pdes.c */
#define IMAGE_FILE 5 /* for option D_IMAGE */
#define COMPARISON 0 /* set to 1 to compare two different patterns (beta) */
@@ -83,22 +85,22 @@
#define P_PERCOL 0.15 /* probability of having a circle in C_RAND_PERCOL arrangement */
#define NPOISSON 1000 /* number of points for Poisson C_RAND_POISSON arrangement */
#define PDISC_FACTOR 2.3 /* controls density of Poisson disc process (default: 3.25) */
#define PDISC_FACTOR 5.0 /* controls density of Poisson disc process (default: 3.25) */
#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 0.6 /* parameter controlling the dimensions of domain */
#define MU 0.15 /* parameter controlling the dimensions of domain */
#define MU_B 0.45 /* parameter controlling the dimensions of domain */
#define NPOLY 5 /* number of sides of polygon */
#define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */
#define LAMBDA 0.2 /* parameter controlling the dimensions of domain */
#define MU 0.008 /* parameter controlling the dimensions of domain */
#define MU_B 0.08 /* parameter controlling the dimensions of domain */
#define NPOLY 9 /* number of sides of polygon */
#define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */
#define MDEPTH 6 /* depth of computation of Menger gasket */
#define MRATIO 3 /* ratio defining Menger gasket */
#define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */
#define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */
#define FOCI 1 /* set to 1 to draw focal points of ellipse */
#define NGRIDX 4 /* number of grid point for grid of disks */
#define NGRIDY 14 /* number of grid point for grid of disks */
#define NGRIDX 10 /* number of grid point for grid of disks */
#define NGRIDY 10 /* number of grid point for grid of disks */
#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 */
@@ -122,19 +124,20 @@
/* Physical parameters of wave equation */
#define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */
#define OSCILLATE_LEFT 0 /* set to 1 to add oscilating boundary condition on the left */
#define OSCILLATE_LEFT 1 /* set to 1 to add oscilating boundary condition on the left */
#define OSCILLATE_TOPBOT 0 /* set to 1 to enforce a planar wave on top and bottom boundary */
#define OSCILLATION_SCHEDULE 42 /* oscillation schedule, see list in global_pdes.c */
#define OSCIL_YMAX 0.0375 /* defines oscilling beam range */
#define OSCILLATION_SCHEDULE 0 /* oscillation schedule, see list in global_pdes.c */
#define OSCIL_YMAX 1.2 /* 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.006 /* frequency of periodic excitation */
#define AMPLITUDE 2.0 /* amplitude of periodic excitation */
#define OMEGA 0.005 /* frequency of periodic excitation */
#define OMEGA_B 0.01 /* frequency of periodic excitation */
#define AMPLITUDE 1.0 /* amplitude of periodic excitation */
#define ACHIRP 0.25 /* acceleration coefficient in chirp */
#define DAMPING 0.0 /* damping of periodic excitation */
#define COURANT 0.12 /* Courant number in medium B */
#define COURANT 0.09 /* 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 */
@@ -149,12 +152,12 @@
/* Increasing COURANT speeds up the simulation, but decreases accuracy */
/* For similar wave forms, COURANT^2*GAMMA should be kept constant */
#define ADD_OSCILLATING_SOURCE 1 /* set to 1 to add an oscillating wave source */
#define OSCILLATING_SOURCE_PERIOD 12.0 /* period of oscillating source */
#define ADD_OSCILLATING_SOURCE 0 /* set to 1 to add an oscillating wave source */
#define OSCILLATING_SOURCE_PERIOD 7.0 /* period of oscillating source */
#define ALTERNATE_OSCILLATING_SOURCE 1 /* set to 1 to alternate sign of oscillating source */
#define N_SOURCES 4 /* number of sources, for option draw_sources */
#define N_SOURCES 1 /* number of sources, for option draw_sources */
#define ALTERNATE_SOURCE_PHASES 0 /* set to 1 to alternate initial phases of sources */
#define MAX_PULSING_TIME 5000 /* max time for adding pulses */
#define MAX_PULSING_TIME 13 /* 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 */
@@ -165,14 +168,15 @@
/* Boundary conditions, see list in global_pdes.c */
#define B_COND 2
#define B_COND 3
#define BC_NEUMANN 1 /* set to 1 to use Neumann boundary conditions on domain */
/* Parameters for length and speed of simulation */
#define NSTEPS 3200 /* number of frames of movie */
#define NVID 12 /* number of iterations between images displayed on screen */
#define NSTEPS 3000 /* number of frames of movie */
#define NVID 10 /* number of iterations between images displayed on screen */
#define NSEG 1000 /* number of segments of boundary */
#define INITIAL_TIME 0 /* time after which to start saving frames */
#define INITIAL_TIME 50 /* 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) */
@@ -187,9 +191,9 @@
/* Parameters of initial condition */
#define INITIAL_AMP 0.5 /* amplitude of initial condition */
#define INITIAL_VARIANCE 0.00025 /* variance of initial condition */
#define INITIAL_WAVELENGTH 0.015 /* wavelength of initial condition */
#define INITIAL_AMP 1.0 /* amplitude of initial condition */
#define INITIAL_VARIANCE 0.0003 /* variance of initial condition */
#define INITIAL_WAVELENGTH 0.01 /* wavelength of initial condition */
/* Plot type, see list in global_pdes.c */
@@ -199,8 +203,8 @@
/* Color schemes */
#define COLOR_PALETTE 15 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 18 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE 18 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 17 /* Color palette, see list in global_pdes.c */
#define BLACK 1 /* background */
@@ -208,12 +212,12 @@
#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 0.75 /* max range of color (default: 1.0) */
#define COLOR_RANGE 0.8 /* 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.0 /* additional shift for wave amplitude */
#define VSCALE_AMPLITUDE 0.5 /* additional scaling factor for wave amplitude */
#define VSHIFT_AMPLITUDE -0.2 /* additional shift for wave amplitude */
#define VSCALE_AMPLITUDE 1.0 /* additional scaling factor for wave amplitude */
#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 */
@@ -222,7 +226,7 @@
#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 */
#define SHADE_SCALE_2D 0.01 /* lower value increases sensitivity of shading */
#define SHADE_SCALE_2D 0.5 /* lower value increases sensitivity of shading */
#define COLORHUE 260 /* initial hue of water color for scheme C_LUM */
#define COLORDRIFT 0.0 /* how much the color hue drifts during the whole simulation */
@@ -238,10 +242,10 @@
#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 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 DRAW_WAVE_PROFILE 1 /* 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 WAVE_PROFILE_X 1.75 /* value of x 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*/
@@ -272,7 +276,7 @@
/* For debugging purposes only */
#define FLOOR 0 /* set to 1 to limit wave amplitude to VMAX */
#define VMAX 10.0 /* max value of wave amplitude */
#define VMAX 2.5 /* max value of wave amplitude */
#define MEAN_FLUX (PLOT == P_TOTAL_ENERGY_FLUX)||(PLOT_B == P_TOTAL_ENERGY_FLUX)
@@ -296,12 +300,12 @@ FILE *time_series_left, *time_series_right;
// void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX], double *psi_out[NX],
// short int *xy_in[NX])
void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX],
short int *xy_in[NX], double *tcc[NX], double *tc[NX], double *tgamma[NX])
short int *xy_in[NX], double *tcc[NX], double *tc[NX], double *tgamma[NX], t_neumann_bc neumann_bc[N_NEUMANN_POINTS])
/* time step of field evolution */
/* phi is value of field at time t, psi at time t-1 */
/* this version of the function has been rewritten in order to minimize the number of if-branches */
{
int i, j, iplus, iminus, jplus, jminus, ij[2];
int i, j, iplus, iminus, jplus, jminus, ij[2], nc = 0;
double delta, x, y, c, cc, gamma, tb_shift;
static long time = 0;
// static double tc[NX][NY];
@@ -358,7 +362,7 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
/* evolution in the bulk */
for (i=1; i<NX-1; i++){
for (j=1; j<NY-1; j++){
if ((TWOSPEEDS)||(xy_in[i][j] != 0)){
if ((TWOSPEEDS)||(xy_in[i][j] > 0)){
x = phi_in[i][j];
y = psi_in[i][j];
@@ -372,6 +376,20 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
}
}
if (BC_NEUMANN)
{
#pragma omp parallel for private(i,j,nc)
for (nc= 0; nc < n_neumann_points; nc++) if (neumann_bc[nc].active)
{
i = neumann_bc[nc].i;
j = neumann_bc[nc].j;
phi_out[i][j] = neumann_bc[nc].nx*phi_in[neumann_bc[nc].nbi1][neumann_bc[nc].nbj1];
phi_out[i][j] += neumann_bc[nc].ny*phi_in[neumann_bc[nc].nbi2][neumann_bc[nc].nbj2];
}
}
/* left boundary */
// if (OSCILLATE_LEFT) for (j=bc_jmin+1; j<bc_jmax; j++) phi_out[0][j] = oscillating_bc(time, j);
// else
@@ -572,21 +590,27 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
}
/* for debugging purposes/if there is a risk of blow-up */
if (FLOOR) for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
if (FLOOR)
{
#pragma omp parallel for private(i,j)
for (i=0; i<NX; i++)
{
for (j=0; j<NY; j++)
{
if (xy_in[i][j] != 0)
{
if (phi_out[i][j] > VMAX) phi_out[i][j] = VMAX;
if (phi_out[i][j] < -VMAX) phi_out[i][j] = -VMAX;
else if (phi_out[i][j] < -VMAX) phi_out[i][j] = -VMAX;
// if (psi_out[i][j] > VMAX) psi_out[i][j] = VMAX;
// if (psi_out[i][j] < -VMAX) psi_out[i][j] = -VMAX;
}
}
}
}
}
void evolve_wave(double *phi[NX], double *psi[NX], double *tmp[NX], short int *xy_in[NX], double *tcc_table[NX], double *tc_table[NX], double *tgamma_table[NX])
void evolve_wave(double *phi[NX], double *psi[NX], double *tmp[NX], short int *xy_in[NX], double *tcc_table[NX], double *tc_table[NX], double *tgamma_table[NX], t_neumann_bc neumann_bc[N_NEUMANN_POINTS])
/* time step of field evolution */
/* phi is value of field at time t, psi at time t-1 */
{
@@ -596,11 +620,11 @@ void evolve_wave(double *phi[NX], double *psi[NX], double *tmp[NX], short int *x
// At the beginning w[t] is saved in phi, w[t-1] in psi and tmp is space
// for the next wave state w[t+1]. Take w[t] and w[t-1] to calculate the
// next wave state. Write this new state in temp
evolve_wave_half(phi, psi, tmp, xy_in, tcc_table, tc_table, tgamma_table);
evolve_wave_half(phi, psi, tmp, xy_in, tcc_table, tc_table, tgamma_table, neumann_bc);
// now w[t] is saved in tmp, w[t-1] in phi and the result is written to psi
evolve_wave_half(tmp, phi, psi, xy_in, tcc_table, tc_table, tgamma_table);
evolve_wave_half(tmp, phi, psi, xy_in, tcc_table, tc_table, tgamma_table, neumann_bc);
// now w[t] is saved in psi, w[t-1] in tmp and the result is written to phi
evolve_wave_half(psi, tmp, phi, xy_in, tcc_table, tc_table, tgamma_table);
evolve_wave_half(psi, tmp, phi, xy_in, tcc_table, tc_table, tgamma_table, neumann_bc);
// now w[t] is saved in phi, w[t-1] in psi and tmp is free again to take
// the new wave state w[t+1] in the next call to this function, thus
// matching the given parameter names again
@@ -637,6 +661,7 @@ void animation()
int image_counter = 0, add_counter[N_SOURCES];
long int wave_value;
t_wave_packet *packet;
t_neumann_bc *neumann_bc;
t_wave_source wave_source[25];
if (SAVE_TIME_SERIES)
@@ -670,11 +695,14 @@ void animation()
init_wave_packets(packet, WAVE_PACKET_RADIUS);
}
if (BC_NEUMANN)
neumann_bc = (t_neumann_bc *)malloc(N_NEUMANN_POINTS*sizeof(t_neumann_bc));
if (ALTERNATE_SOURCE_PHASES) for (i=0; i<N_SOURCES; i++) sign[i] = pow(-1.0,(double)i);
else for (i=0; i<N_SOURCES; i++) sign[i] = 1.0;
/* initialise positions and radii of circles */
if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN == D_CIRCLES_IN_RECT)) ncircles = init_circle_config(circles);
if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN == D_CIRCLES_IN_RECT)||(B_DOMAIN == D_CIRCLES_NEUMANN)) ncircles = init_circle_config(circles);
else if (B_DOMAIN == D_POLYGONS) ncircles = init_polygon_config(polygons);
printf("Polygons initialized\n");
@@ -767,6 +795,10 @@ void animation()
// printf("xleft = (%.3f, %.3f) xright = (%.3f, %.3f)\n", xin_left, yin_left, xin_right, yin_right);
init_wave_flat(phi, psi, xy_in);
// init_circular_wave(LAMBDA, 0.0, phi, psi, xy_in);
// add_circular_wave(1.0, -LAMBDA, 0.0, phi, psi, xy_in);
init_ior_2d(xy_in, tcc_table, tc_table, tgamma_table, ior_angle);
if (FADE_IN_OBSTACLE) init_fade_table(tcc_table, xy_in, fade_table);
@@ -819,6 +851,9 @@ void animation()
// add_drop_to_wave(1.0, -0.7, 0.0, phi, psi);
// add_drop_to_wave(1.0, 0.0, -0.7, phi, psi);
/* initialize Neumann boundary data */
if (BC_NEUMANN) n_neumann_points = init_neumann_bc(xy_in, neumann_bc, circles, 0);
blank();
glColor3f(0.0, 0.0, 0.0);
// draw_wave(phi, psi, xy_in, 1.0, 0, PLOT);
@@ -861,7 +896,7 @@ void animation()
else draw_wave_epalette(phi, psi, total_energy, average_energy, total_flux, color_scale, xy_in, scale, i, PLOT, COLOR_PALETTE, 0, 1.0);
for (j=0; j<NVID; j++)
{
evolve_wave(phi, psi, tmp, xy_in, tcc_table, tc_table, tgamma_table);
evolve_wave(phi, psi, tmp, xy_in, tcc_table, tc_table, tgamma_table, neumann_bc);
if (SAVE_TIME_SERIES)
{
wave_value = (long int)(phi[sample_left[0]][sample_left[1]]*1.0e16);
@@ -880,14 +915,11 @@ void animation()
/* add oscillating waves */
for (source = 0; source < N_SOURCES; source++)
{
wave_source_x[source] = polyarc[source].xc;
wave_source_y[source] = polyarc[source].yc;
wave_source_x[source] = -2.0;
wave_source_y[source] = 0.0;
source_amp[source] = INITIAL_AMP;
}
source_periods[0] = OSCILLATING_SOURCE_PERIOD;
source_periods[1] = OSCILLATING_SOURCE_PERIOD*3/4;
source_periods[2] = OSCILLATING_SOURCE_PERIOD*2/3;
source_periods[3] = OSCILLATING_SOURCE_PERIOD/2;
for (source = 0; source < N_SOURCES; source++)
{
dperiod = source_periods[source];
@@ -1032,6 +1064,8 @@ void animation()
if (ADD_WAVE_PACKET_SOURCES) free(packet);
if (BC_NEUMANN) free(neumann_bc);
if (SAVE_TIME_SERIES)
{
fclose(time_series_left);
wave_comparison.c
+101 -65
View File
@@ -44,11 +44,11 @@
#include <time.h>
#define MOVIE 0 /* set to 1 to generate movie */
#define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */
#define DOUBLE_MOVIE 0 /* set to 1 to produce movies for wave height and energy simultaneously */
#define SAVE_MEMORY 1 /* set to 1 to save memory when writing tiff images */
#define 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 183 /* choice of index of refraction, see list in global_pdes.c */
#define IOR_B 183 /* choice of index of refraction, see list in global_pdes.c */
#define IOR_TOTAL_TURNS 1.5 /* total angle of rotation for IOR_PERIODIC_WELLS_ROTATING */
@@ -63,15 +63,12 @@
#define NY 2300 /* number of grid points on y axis */
#define YMID 1150 /* mid point of display */
#define XMIN -2.2
#define XMAX 3.8 /* x interval */
#define YMIN -1.796875
#define YMAX 1.796875 /* y interval for 9/16 aspect ratio */
// #define XMIN -2.0
// #define XMAX 2.0 /* x interval */
// #define YMIN -1.197916667
// #define YMAX 1.197916667 /* y interval for 9/16 aspect ratio */
#define XMIN -0.75
#define XMAX 3.25 /* x interval */
#define YMIN -1.197916667
#define YMAX 1.197916667 /* y interval for 9/16 aspect ratio */
#define HIGHRES 1 /* set to 1 if resolution of grid is double that of displayed image */
@@ -79,16 +76,18 @@
/* Choice of the billiard table */
#define B_DOMAIN 196 /* choice of domain shape, see list in global_pdes.c */
#define B_DOMAIN_B 196 /* choice of domain shape, see list in global_pdes.c */
#define B_DOMAIN 20 /* choice of domain shape, see list in global_pdes.c */
#define B_DOMAIN_B 20 /* choice of domain shape, see list in global_pdes.c */
#define CIRCLE_PATTERN 13 /* pattern of circles, see list in global_pdes.c */
#define CIRCLE_PATTERN_B 13 /* pattern of circles, see list in global_pdes.c */
#define CIRCLE_PATTERN 5 /* pattern of circles, see list in global_pdes.c */
#define CIRCLE_PATTERN_B 61 /* pattern of circles, see list in global_pdes.c */
#define SYMMETRIC_CIRCLE_PATTERNS 0 /* set to 1 to have symmetric patterns in top and bottom */
#define NO_CIRCLES_ON_BOUNDARY 0 /* set to 1 to avoid circles overlapping boundary */
#define IMAGE_FILE 5 /* for option D_IMAGE */
#define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */
#define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */
#define PDISC_FACTOR 3.25 /* controls density of Poisson disc process (default: 3.25) */
#define PDISC_FACTOR 3.5 /* controls density of Poisson disc process (default: 3.25) */
#define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */
#define RANDOM_POLY_ANGLE_B 0 /* 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 */
@@ -99,10 +98,10 @@
#define HEX_NONUNIF_COMPRESSSION 0.15 /* compression factor for HEX_NONUNIF pattern */
#define HEX_NONUNIF_COMPRESSSION_B -0.15 /* compression factor for HEX_NONUNIF pattern */
#define LAMBDA 1.0 /* parameter controlling the dimensions of domain */
#define MU 0.06 /* parameter controlling the dimensions of domain */
#define MU_B 0.06 /* parameter controlling the dimensions of domain */
#define MUB 0.06 /* parameter controlling the dimensions of domain */
#define LAMBDA 0.4 /* parameter controlling the dimensions of domain */
#define MU 0.008 /* parameter controlling the dimensions of domain */
#define MU_B 0.008 /* parameter controlling the dimensions of domain */
#define MUB 0.008 /* parameter controlling the dimensions of domain */
#define NPOLY 3 /* number of sides of polygon */
#define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */
#define APOLY_B 2.0 /* angle by which to turn polygon, in units of Pi/2 */
@@ -112,8 +111,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 17 /* number of grid point for grid of disks */
#define NGRIDY 12 /* number of grid point for grid of disks */
#define NGRIDX 14 /* number of grid point for grid of disks */
#define NGRIDY 16 /* number of grid point for grid of disks */
#define WALL_WIDTH 0.012 /* width of channels/wall separating lenses */
#define WALL_WIDTH_B 0.012 /* width of channels/wall separating lenses */
#define WALL_WIDTH_RND 0.0 /* proportion of width of width for random arrangements */
@@ -137,12 +136,14 @@
/* 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 1 /* set to 1 to enforce a planar wave on top and bottom boundary */
#define OSCIL_YMID -0.9 /* defines oscilling beam midpoint */
#define OSCIL_YMID 0.0 /* defines oscilling beam midpoint */
#define OSCILLATION_SCHEDULE 0 /* oscillation schedule, see list in global_pdes.c */
#define OMEGA 0.002 /* frequency of periodic excitation */
#define OMEGA 0.01 /* frequency of periodic excitation */
#define OMEGA_B 0.015 /* frequency of periodic excitation */
#define AMPLITUDE 1.0 /* amplitude of periodic excitation */
#define DAMPING 0.0 /* damping of periodic excitation */
#define COURANT 0.1 /* Courant number */
@@ -170,13 +171,14 @@
/* Boundary conditions, see list in global_pdes.c */
#define B_COND 4
#define BC_NEUMANN 1 /* set to 1 to use Neumann boundary conditions on domain */
/* Parameters for length and speed of simulation */
#define NSTEPS 2300 /* number of frames of movie */
#define NVID 7 /* number of iterations between images displayed on screen */
#define NSTEPS 3000 /* number of frames of movie */
#define NVID 10 /* number of iterations between images displayed on screen */
#define NSEG 100 /* number of segments of boundary */
#define INITIAL_TIME 0 /* time after which to start saving frames */
#define INITIAL_TIME 100 /* time after which to start saving frames */
#define COMPUTE_ENERGIES 0 /* set to 1 to compute and print energies */
#define BOUNDARY_WIDTH 2 /* width of billiard boundary */
@@ -202,8 +204,8 @@
/* Color schemes */
#define COLOR_PALETTE 11 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 12 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE 18 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 11 /* Color palette, see list in global_pdes.c */
#define BLACK 1 /* background */
#define BLACK_TEXT 1 /* set to 1 to write text in black instead of white */
@@ -212,11 +214,11 @@
#define SCALE 0 /* set to 1 to adjust color scheme to variance of field */
#define SLOPE 1.0 /* sensitivity of color on wave amplitude */
#define COLOR_RANGE 1.0 /* max range of color (default: 1.0) */
#define COLOR_RANGE 0.8 /* 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.0 /* additional shift for wave amplitude */
#define VSCALE_AMPLITUDE 1.0 /* additional scaling factor for wave amplitude */
#define E_SCALE 300.0 /* scaling factor for energy representation */
#define LOG_SCALE 1.0 /* scaling factor for energy log representation */
@@ -231,21 +233,29 @@
#define HUEMEAN 220.0 /* mean value of hue for color scheme C_HUE */
#define HUEAMP -220.0 /* amplitude of variation of hue for color scheme C_HUE */
#define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */
#define DRAW_COLOR_SCHEME 0 /* set to 1 to plot the color scheme */
#define COLORBAR_RANGE 1.2 /* scale of color scheme bar */
#define COLORBAR_RANGE_B 5.0 /* scale of color scheme bar for 2nd part */
#define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */
/* For debugging purposes only */
#define FLOOR 0 /* set to 1 to limit wave amplitude to VMAX */
#define VMAX 5.0 /* max value of wave amplitude */
#define FLOOR 1 /* set to 1 to limit wave amplitude to VMAX */
#define VMAX 2.0 /* max value of wave amplitude */
/* the following constants are only used by wave_billiard and wave_3d so far */
#define COMPARISON 0 /* set to 1 to compare two different patterns */
#define OSCILLATION_SCHEDULE 3 /* oscillation schedule, see list in global_pdes.c */
#define ACHIRP 0.2 /* acceleration coefficient in chirp */
#define DAMPING 0.0 /* damping of periodic excitation */
#define DRAW_WAVE_PROFILE 1 /* 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 WAVE_PROFILE_X 1.75 /* value of x 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*/
/* end of constants only used by wave_billiard and wave_3d */
/* for compatibility with sub_wave and sub_maze */
@@ -261,8 +271,6 @@
#define WAVE_PACKET_SOURCE_TYPE 1 /* type of wave packet sources */
#define N_WAVE_PACKETS 15 /* number of wave packets */
#define OSCIL_LEFT_YSHIFT 0.0 /* y-dependence of left oscillation (for non-horizontal waves) */
#define DRAW_WAVE_PROFILE 0 /* set to 1 to draw a profile of the wave */
#define VERTICAL_WAVE_PROFILE 0 /* set to 1 to draw wave profile vertically */
#define DRAW_WAVE_TIMESERIES 0 /* set to 1 to draw a time series of the wave */
#define RADIUS_FACTOR 0.3 /* controls inner radius for C_RING arrangements */
#define OSCIL_YMAX 0.35 /* defines oscillation range */
@@ -273,11 +281,6 @@
#define MESSAGE_LSPACE 48 /* length of space for Morse code message */
#define MESSAGE_INITIAL_TIME 100 /* initial time before starting message for Morse code message */
#define AVRG_E_FACTOR 0.95 /* controls time window size in P_AVERAGE_ENERGY scheme */
#define HORIZONTAL_WAVE_PROFILE 0 /* set to 1 to draw wave profile vertically */
#define WAVE_PROFILE_X 2.1 /* value of x to sample wave profile */
#define WAVE_PROFILE_Y -1.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 1 /* set to 1 to draw time-average of wave profile squared*/
#define TIMESERIES_NVALUES 400 /* number of values plotted in time series */
#define DRAW_WAVE_SOURCE 0 /* set to 1 to draw source of wave at (wave_source_x, wave_source_y) */
#define HRES 1 /* dummy, only used by rde.c */
@@ -285,7 +288,7 @@
#define WAVE_PACKET_SHIFT 200.0 /* time shift between wave packets (in oscillation periods) */
#define FADE_IN_OBSTACLE 0 /* set to 1 to fade color inside obstacles */
#define SHADE_2D 1 /* set to 1 to add pseudo-3d shading effect */
#define SHADE_SCALE_2D 0.05 /* lower value increases sensitivity of shading */
#define SHADE_SCALE_2D 0.025 /* lower value increases sensitivity of shading */
#define MEAN_FLUX (PLOT == P_TOTAL_ENERGY_FLUX)||(PLOT_B == P_TOTAL_ENERGY_FLUX)
#define XYIN_INITIALISED (B_DOMAIN == D_IMAGE)
@@ -308,11 +311,11 @@ FILE *monitor_sources;
/*********************/
void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX],
short int *xy_in[NX], double *tcc[NX])
short int *xy_in[NX], double *tcc[NX], t_neumann_bc neumann_bc[N_NEUMANN_POINTS])
/* time step of field evolution */
/* phi is value of field at time t, psi at time t-1 */
{
int i, j, iplus, iminus, jplus, jminus, jmid = NY/2;
int i, j, iplus, iminus, jplus, jminus, jmid = NY/2, nc, j1, j2;
double delta, x, y, c, cc, gamma, tb_shift;
static long time = 0;
static double tc[NX][NY], tgamma[NX][NY];
@@ -374,10 +377,24 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
}
}
if (BC_NEUMANN)
{
#pragma omp parallel for private(i,j,nc)
for (nc= 0; nc < n_neumann_points; nc++) if (neumann_bc[nc].active)
{
i = neumann_bc[nc].i;
j = neumann_bc[nc].j;
phi_out[i][j] = neumann_bc[nc].nx*phi_in[neumann_bc[nc].nbi1][neumann_bc[nc].nbj1];
phi_out[i][j] += neumann_bc[nc].ny*phi_in[neumann_bc[nc].nbi2][neumann_bc[nc].nbj2];
}
}
/* left boundary */
if (OSCILLATE_LEFT) {
for (j=1; j<jmid-1; j++) phi_out[0][j] = AMPLITUDE*cos((double)time*OMEGA);
for (j=jmid+1; j<NY-1; j++) phi_out[0][j] = AMPLITUDE*cos((double)time*OMEGA);
for (j=1; j<jmid-1; j++) phi_out[0][j] = oscillating_bc(time, j);
for (j=jmid+1; j<NY-1; j++) phi_out[0][j] = oscillating_bc(time, j);
// for (j=1; j<jmid-1; j++) phi_out[0][j] = AMPLITUDE*cos((double)time*OMEGA);
// for (j=jmid+1; j<NY-1; j++) phi_out[0][j] = AMPLITUDE*cos((double)time*OMEGA);
}
else for (j=1; j<NY-1; j++) if ((j!=jmid-1)&&(j!=jmid)) {
if ((TWOSPEEDS)||(xy_in[0][j] != 0)){
@@ -714,12 +731,18 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
}
/* for debugging purposes/if there is a risk of blow-up */
if (FLOOR) for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
if (FLOOR)
{
#pragma omp parallel for private(i,j)
for (i=0; i<NX; i++)
{
for (j=0; j<NY; j++)
{
if (xy_in[i][j] != 0)
{
if (phi_out[i][j] > VMAX) phi_out[i][j] = VMAX;
if (phi_out[i][j] < -VMAX) phi_out[i][j] = -VMAX;
else if (phi_out[i][j] < -VMAX) phi_out[i][j] = -VMAX;
}
}
}
}
@@ -727,13 +750,13 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
}
void evolve_wave(double *phi[NX], double *psi[NX], double *tmp[NX], short int *xy_in[NX], double *tcc_table[NX])
void evolve_wave(double *phi[NX], double *psi[NX], double *tmp[NX], short int *xy_in[NX], double *tcc_table[NX], t_neumann_bc neumann_bc[N_NEUMANN_POINTS])
/* time step of field evolution */
/* phi is value of field at time t, psi at time t-1 */
{
evolve_wave_half(phi, psi, tmp, xy_in, tcc_table);
evolve_wave_half(tmp, phi, psi, xy_in, tcc_table);
evolve_wave_half(psi, tmp, phi, xy_in, tcc_table);
evolve_wave_half(phi, psi, tmp, xy_in, tcc_table, neumann_bc);
evolve_wave_half(tmp, phi, psi, xy_in, tcc_table, neumann_bc);
evolve_wave_half(psi, tmp, phi, xy_in, tcc_table, neumann_bc);
}
@@ -761,6 +784,7 @@ void animation()
short int *xy_in[NX];
int i, j, s, counter = 0, k, first_source = 1, fade, source, resol = HIGHRES + 1, phase, add_counter[N_SOURCES];
t_wave_source wave_source[NSOURCES];
t_neumann_bc *neumann_bc;
monitor_sources = fopen("monitor_sources.dat", "w");
@@ -780,13 +804,19 @@ void animation()
if (MEAN_FLUX) total_flux = (double *)malloc(4*NX*NY*sizeof(double));
if (BC_NEUMANN)
neumann_bc = (t_neumann_bc *)malloc(N_NEUMANN_POINTS*sizeof(t_neumann_bc));
if (ALTERNATE_SOURCE_PHASES) for (i=0; i<N_SOURCES; i++) sign[i] = pow(-1.0,(double)i);
else for (i=0; i<N_SOURCES; i++) sign[i] = 1.0;
/* initialise positions and radii of circles */
printf("initializing circle configuration\n");
if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN_B == D_CIRCLES)) init_circle_config_comp(circles);
if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN_B == D_CIRCLES)||(B_DOMAIN == D_CIRCLES_NEUMANN)) init_circle_config_comp(circles);
if ((B_DOMAIN == D_POLYGONS)|(B_DOMAIN_B == D_POLYGONS)) init_polygon_config_comp(polygons);
circles[126].top = 0;
// for (i=0; i<ncircles; i++) printf("polygon %i at (%.3f, %.3f) radius %.3f\n", i, polygons[i].xc, polygons[i].yc, polygons[i].radius);
/* initialise polyline for von Koch and similar domains */
@@ -800,8 +830,11 @@ void animation()
// for (i=0; i<npolyarc; i++) printf("polyarc %i: center (%.3f, %.3f) angles (%.3f, %.3f)\n", i, polyarc[i].xc, polyarc[i].yc, polyarc[i].angle1*180.0/PI, (polyarc[i].angle1 + polyarc[i].dangle)*180.0/PI);
// for (i=0; i<npolyarc_b; i++) printf("polyarc_b %i: center (%.3f, %.3f) angles (%.3f, %.3f)\n", i, polyarc_b[i].xc, polyarc_b[i].yc, polyarc_b[i].angle1*180.0/PI, (polyarc_b[i].angle1 + polyarc_b[i].dangle)*180.0/PI);
for (i=0; i<ncircles; i++) printf("circle %i: (%.3f, %.3f) radius %.3f\n", i, circles[i].xc, circles[i].yc, circles[i].radius);
for (i=0; i<ncircles_b; i++) printf("circle %i: (%.3f, %.3f) radius %.3f\n", i, circles_b[i].xc, circles_b[i].yc, circles_b[i].radius);
for (i=0; i<ncircles; i++) if (circles[i].active)
printf("circle %i: (%.3f, %.3f) radius %.3f top = %i\n", i, circles[i].xc, circles[i].yc, circles[i].radius, circletop[i]);
for (i=0; i<ncircles_b; i++) if (circles_b[i].active)
printf("circle %i: (%.3f, %.3f) radius %.3f\n", i, circles_b[i].xc, circles_b[i].yc, circles_b[i].radius);
/* initialize total energy table */
if ((PLOT == P_MEAN_ENERGY)||(PLOT_B == P_MEAN_ENERGY)||(PLOT == P_LOG_MEAN_ENERGY)||(PLOT_B == P_LOG_MEAN_ENERGY))
@@ -839,11 +872,14 @@ void animation()
printf("computed energies\n");
}
/* initialize Neumann boundary data */
if (BC_NEUMANN) n_neumann_points = init_neumann_bc(xy_in, neumann_bc, circles, 1);
blank();
glColor3f(0.0, 0.0, 0.0);
printf("drawing wave\n");
// draw_wave_comp(phi, psi, xy_in, 1.0, 0, PLOT);
draw_wave_comp_highres_palette(2, phi, psi, total_energy, average_energy, total_flux, color_scale, xy_in, 1.0, 0, PLOT, COLOR_PALETTE, 0, 1.0);
draw_wave_comp_highres_palette(2, phi, psi, total_energy, average_energy, total_flux, color_scale, xy_in, 1.0, 0, PLOT, COLOR_PALETTE, 1, 0, 1.0);
printf("drawing billiard\n");
draw_billiard_comp(0, 1.0);
@@ -869,11 +905,11 @@ void animation()
else scale = 1.0;
// draw_wave_comp(phi, psi, xy_in, scale, i, PLOT);
draw_wave_comp_highres_palette(2, phi, psi, total_energy, average_energy, total_flux, color_scale, xy_in, 1.0, i, PLOT, COLOR_PALETTE, 0, 1.0);
draw_wave_comp_highres_palette(2, phi, psi, total_energy, average_energy, total_flux, color_scale, xy_in, 1.0, i, PLOT, COLOR_PALETTE, 0, 0, 1.0);
for (j=0; j<NVID; j++)
{
evolve_wave(phi, psi, tmp, xy_in, tcc_table);
evolve_wave(phi, psi, tmp, xy_in, tcc_table, neumann_bc);
// if (i % 10 == 9) oscillate_linear_wave(0.2*scale, 0.15*(double)(i*NVID + j), -1.5, YMIN, -1.5, YMAX, phi, psi);
}
@@ -982,7 +1018,7 @@ void animation()
else /*if (DOUBLE_MOVIE)*/
{
// draw_wave_comp(phi, psi, xy_in, scale, i, PLOT_B);
draw_wave_comp_highres_palette(2, phi, psi, total_energy, average_energy, total_flux, color_scale, xy_in, 1.0, i, PLOT_B, COLOR_PALETTE_B, 0, 1.0);
draw_wave_comp_highres_palette(2, phi, psi, total_energy, average_energy, total_flux, color_scale, xy_in, 1.0, i, PLOT_B, COLOR_PALETTE_B, 0, 0, 1.0);
draw_billiard_comp(0, 1.0);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0);
glutSwapBuffers();
@@ -1013,7 +1049,7 @@ void animation()
if (DOUBLE_MOVIE)
{
// draw_wave_comp(phi, psi, xy_in, scale, i, PLOT);
draw_wave_comp_highres_palette(2, phi, psi, total_energy, average_energy, total_flux, color_scale, xy_in, 1.0, NSTEPS, PLOT, COLOR_PALETTE, 0, 1.0);
draw_wave_comp_highres_palette(2, phi, psi, total_energy, average_energy, total_flux, color_scale, xy_in, 1.0, NSTEPS, PLOT, COLOR_PALETTE, 0, 0, 1.0);
draw_billiard_comp(0, 1.0);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, 0, 1.0);
glutSwapBuffers();
@@ -1022,7 +1058,7 @@ void animation()
else for (i=0; i<MID_FRAMES; i++)
{
fade_value = 1.0 - (double)i/(double)MID_FRAMES;
draw_wave_comp_highres_palette(2, phi, psi, total_energy, average_energy, total_flux, color_scale, xy_in, 1.0, NSTEPS, PLOT, COLOR_PALETTE, 1, fade_value);
draw_wave_comp_highres_palette(2, phi, psi, total_energy, average_energy, total_flux, color_scale, xy_in, 1.0, NSTEPS, PLOT, COLOR_PALETTE, 0, 1, fade_value);
draw_billiard_comp(1, fade_value);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, 1, fade_value);
if (!NO_EXTRA_BUFFER_SWAP) glutSwapBuffers();
@@ -1031,7 +1067,7 @@ void animation()
if (DOUBLE_MOVIE)
{
// draw_wave_comp(phi, psi, xy_in, scale, i, PLOT_B);
draw_wave_comp_highres_palette(2, phi, psi, total_energy, average_energy, total_flux, color_scale, xy_in, 1.0, NSTEPS, PLOT_B, COLOR_PALETTE_B, 0, 1.0);
draw_wave_comp_highres_palette(2, phi, psi, total_energy, average_energy, total_flux, color_scale, xy_in, 1.0, NSTEPS, PLOT_B, COLOR_PALETTE_B, 1, 0, 1.0);
draw_billiard_comp(0, 1.0);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0);
glutSwapBuffers();
@@ -1040,7 +1076,7 @@ void animation()
else for (i=0; i<END_FRAMES; i++)
{
fade_value = 1.0 - (double)i/(double)END_FRAMES;
draw_wave_comp_highres_palette(2, phi, psi, total_energy, average_energy, total_flux, color_scale, xy_in, 1.0, NSTEPS, PLOT_B, COLOR_PALETTE_B, 1, fade_value);
draw_wave_comp_highres_palette(2, phi, psi, total_energy, average_energy, total_flux, color_scale, xy_in, 1.0, NSTEPS, PLOT_B, COLOR_PALETTE_B, 1, 1, fade_value);
draw_billiard_comp(1, fade_value);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 1, fade_value);
glutSwapBuffers();
@@ -1068,7 +1104,7 @@ void animation()
if (MEAN_FLUX) free(total_flux);
if (BC_NEUMANN) free(neumann_bc);
}
wave_energy.c
+3
View File
@@ -126,6 +126,7 @@
#define OSCILLATE_TOPBOT 1 /* set to 1 to enforce a planar wave on top and bottom boundary */
#define OMEGA 0.0125 /* frequency of periodic excitation */
#define OMEGA_B 0.01 /* frequency of periodic excitation */
#define AMPLITUDE 0.5 /* amplitude of periodic excitation */
#define COURANT 0.1 /* Courant number */
#define COURANTB 0.035 /* Courant number in medium B */
@@ -275,6 +276,8 @@
#define SHADE_SCALE_2D 0.01 /* lower value increases sensitivity of shading */
#define N_SOURCES 1 /* number of sources, for option draw_sources */
#define XYIN_INITIALISED (B_DOMAIN == D_IMAGE)
#define SYMMETRIC_CIRCLE_PATTERNS 0 /* set to 1 to have symmetric patterns in top and bottom */
#define NO_CIRCLES_ON_BOUNDARY 0 /* set to 1 to avoid circles overlapping boundary */
double light[2] = {0.40824829, 0.816496581}; /* location of light source for SHADE_2D option*/
/* end of constants only used by sub_wave and sub_maze */
wave_sphere.c
+1
View File
@@ -131,6 +131,7 @@
#define OSCILLATION_SCHEDULE 3 /* oscillation schedule, see list in global_pdes.c */
#define OMEGA 0.001 /* frequency of periodic excitation */
#define OMEGA_B 0.001 /* frequency of periodic excitation */
#define AMPLITUDE 0.8 /* amplitude of periodic excitation */
#define ACHIRP 0.2 /* acceleration coefficient in chirp */
#define DAMPING 0.0 /* damping of periodic excitation */