Compare commits
4 Commits
178da1f6a9
...
283ebb8ebf
Author | SHA1 | Date |
---|---|---|
Nils Berglund | 283ebb8ebf | |
Nils Berglund | d33acfabbf | |
Nils Berglund | 717e35ddf4 | |
Nils Berglund | 10bdaaf598 |
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
9750
Parameters.md
9750
Parameters.md
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
32
README.md
32
README.md
|
@ -6,7 +6,7 @@ C code for videos on YouTube Channel https://www.youtube.com/c/NilsBerglund
|
|||
|
||||
Parameter values used in specific simulations will be gradually added to file `Parameters.md`, `Parameters_June21.md` and so on.
|
||||
|
||||
There are four groups of 6 files, 19 files, 5 files and 4 files.
|
||||
There are four groups of 8 files, 19 files, 5 files and 4 files.
|
||||
In addition the following files handling color schemes have been included:
|
||||
|
||||
1. `hsluv.c`and `hsluv.h` from https://github.com/adammaj1/hsluv-color-gradient
|
||||
|
@ -17,29 +17,31 @@ The following file (beta version) provides support for creating mazes:
|
|||
|
||||
4. `sub_maze.c`
|
||||
|
||||
The file
|
||||
The files
|
||||
|
||||
5. `Earth_Map_Blue_Marble_2002_large.ppm.gz`
|
||||
5. `*.ppm.gz`
|
||||
|
||||
is required by `wave_sphere.c` and should be unzipped before compiling.
|
||||
are required by `wave_sphere.c` and should be unzipped before compiling.
|
||||
|
||||
### Simulations of classical particles in billiards.
|
||||
|
||||
1. *particle_billiard.c*: simulation of a collection of non-interacting particles in a billiard
|
||||
2. *drop_billiard.c*: simulation of an expanding front of particles
|
||||
3. *particle_pinball.c*: variant of `particle_billiard` with some extra statistics plots
|
||||
4. *global_particles.c*: global variables and parameters
|
||||
5. *sub_part_billiard.c*: drawing/computation routines common to `particle_billiard` and `drop_billiard`
|
||||
6. *sub_part_pinball.c*: additional drawing/computation routines for `particle_pinball`
|
||||
3. *particle_pinball.c*: variant of `particle_billiard` with some extra statistics plots
|
||||
4. *billiard_phasespace.c*: variant of `particle_billiard` for phase portraits (only works for certain shapes)
|
||||
5. *global_particles.c*: global variables and parameters
|
||||
6. *sub_part_billiard.c*: drawing/computation routines common to `particle_billiard` and `drop_billiard`
|
||||
7. *sub_part_pinball.c*: additional drawing/computation routines for `particle_pinball`
|
||||
8. *sub_billiard_phasespace.c*: additional drawing/computation routines for `billiard_phasespace`
|
||||
|
||||
|
||||
- Create subfolders `tif_part`, `tif_drop`
|
||||
- Customize constants at beginning of .c file
|
||||
- Compile with `make particle_billiard`, `make_drop_billiard`, etc, or
|
||||
|
||||
`gcc -o particle_billiard particle_billiard.c-O3 -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lXmu -lglut`
|
||||
`gcc -o particle_billiard particle_billiard.c-O3 -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lglut`
|
||||
|
||||
`gcc -o drop_billiard drop_billiard.c-O3 -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lXmu -lglut`
|
||||
`gcc -o drop_billiard drop_billiard.c-O3 -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 lglut`
|
||||
|
||||
- Many laptops claim to have 4 cores, but two of those are virtual. OMP acceleration may be more effective after executing
|
||||
|
||||
|
@ -78,13 +80,13 @@ in the shell before running the program
|
|||
- Customize constants at beginning of .c file
|
||||
- Compile with `make wave_billiard`, etc, or
|
||||
|
||||
`gcc -o wave_billiard wave_billiard.c -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lXmu -lglut -O3 -fopenmp`
|
||||
`gcc -o wave_billiard wave_billiard.c -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lglut -O3 -fopenmp`
|
||||
|
||||
`gcc -o wave_comparison wave_comparison.c -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lXmu -lglut -O3 -fopenmp`
|
||||
`gcc -o wave_comparison wave_comparison.c -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lglut -O3 -fopenmp`
|
||||
|
||||
`gcc -o heat heat.c -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lXmu -lglut -O3 -fopenmp`
|
||||
|
||||
`gcc -o schrodinger schrodinger.c -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lXmu -lglut -O3 -fopenmp`
|
||||
`gcc -o schrodinger schrodinger.c -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lglut -O3 -fopenmp`
|
||||
|
||||
- Many laptops claim to have 4 cores, but two of those are virtual. OMP acceleration may be more effective after executing
|
||||
|
||||
|
@ -109,7 +111,7 @@ in the shell before running the program
|
|||
- Customize constants at beginning of .c file
|
||||
- Compile with `make lennardjones` or
|
||||
|
||||
`gcc -o lennardjones lennardjones.c -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lXmu -lglut -O3 -fopenmp`
|
||||
`gcc -o lennardjones lennardjones.c -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lglut -O3 -fopenmp`
|
||||
|
||||
- Generate movie with
|
||||
|
||||
|
@ -126,7 +128,7 @@ in the shell before running the program
|
|||
- Customize constants at beginning of .c file
|
||||
- Compile with `make percolation` or
|
||||
|
||||
`gcc -o percolation percolation.c -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lXmu -lglut -O3 -fopenmp`
|
||||
`gcc -o percolation percolation.c -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lglut -O3 -fopenmp`
|
||||
|
||||
- Generate movie with
|
||||
|
||||
|
|
Binary file not shown.
|
@ -30,7 +30,7 @@
|
|||
#include <tiffio.h> /* Sam Leffler's libtiff library. */
|
||||
#include <time.h>
|
||||
|
||||
#define MOVIE 1 /* set to 1 to generate movie */
|
||||
#define MOVIE 0 /* set to 1 to generate movie */
|
||||
#define SAVE_MEMORY 1 /* set to 1 to save memory when writing tiff images */
|
||||
|
||||
#define WINWIDTH 1280 /* window width */
|
||||
|
@ -52,7 +52,7 @@
|
|||
|
||||
/* Choice of the billiard table, see global_particles.c */
|
||||
|
||||
#define B_DOMAIN 11 /* choice of domain shape */
|
||||
#define B_DOMAIN 1 /* choice of domain shape */
|
||||
|
||||
#define CIRCLE_PATTERN 6 /* pattern of circles */
|
||||
#define POLYLINE_PATTERN 4 /* pattern of polyline */
|
||||
|
@ -67,8 +67,8 @@
|
|||
#define NGOLDENSPIRAL 2000 /* max number of points for C_GOLDEN_SPIRAL arrandement */
|
||||
#define SDEPTH 2 /* Sierpinski gastket depth */
|
||||
|
||||
#define LAMBDAMIN 6.0 /* parameter controlling shape of domain */
|
||||
#define LAMBDA 20.0 /* parameter controlling shape of domain */
|
||||
#define LAMBDAMIN 0.0 /* parameter controlling shape of domain */
|
||||
#define LAMBDA 1.5 /* parameter controlling shape of domain */
|
||||
#define MU 1.0 /* second parameter controlling shape of billiard */
|
||||
#define FOCI 1 /* set to 1 to draw focal points of ellipse */
|
||||
#define NPOLY 6 /* number of sides of polygon */
|
||||
|
@ -96,7 +96,7 @@
|
|||
#define PRINT_TRAJECTORY_LENGTH 0 /* set to 1 to print length of trajectory 0 */
|
||||
#define PRINT_TRAJECTORY_PERIOD 0 /* set to 1 to print period of trajectory 0 */
|
||||
#define DRAW_LENGTHS_PLOT 0 /* set to 1 to plot trajectory lengths */
|
||||
#define LENGTHS_LOG_SCALE 1 /* set to 1 to use log scale for plot of lengths */
|
||||
#define LENGTHS_LOG_SCALE 0 /* set to 1 to use log scale for plot of lengths */
|
||||
#define MAX_ANGLE 90.0 /* range of angles of trajectory */
|
||||
|
||||
#define NSTEPS 4000 /* number of frames of movie */
|
||||
|
@ -141,7 +141,7 @@
|
|||
#define PLOT_LYAPUNOV 1 /* set to 1 to add plot of Lyapunov exponents */
|
||||
#define LOGSCALEX_LYAP 0 /* set to 1 to use log scale on parameter axis of Lyapunov exponents */
|
||||
#define LYAP_MAX 1.0 /* maximal Lyapunov exponent */
|
||||
#define ADAPT_TO_SYMMETRY 1 /* set to 1 to show only one symmetric part of phase space */
|
||||
#define ADAPT_TO_SYMMETRY 0 /* set to 1 to show only one symmetric part of phase space */
|
||||
#define SYMMETRY_FACTOR 3 /* proportion of phase space to be shown */
|
||||
|
||||
#define PAUSE 1000 /* number of frames after which to pause */
|
||||
|
|
Binary file not shown.
14
global_3d.c
14
global_3d.c
|
@ -65,6 +65,14 @@
|
|||
|
||||
#define VP_HORIZONTAL 0 /* rotate in a horizontal plane */
|
||||
#define VP_ORBIT 1 /* rotate in a plane containing the origin */
|
||||
#define VP_ORBIT2 11 /* rotate in a plane specified by max latitude */
|
||||
#define VP_POLAR 2 /* polar orbit */
|
||||
|
||||
/* Type of digital elevation model */
|
||||
|
||||
#define DEM_EARTH 0 /* DEM of Earth */
|
||||
#define DEM_MARS 1 /* DEM of Mars */
|
||||
#define DEM_MOON 2 /* DEM of the Moon */
|
||||
|
||||
/* macros to avoid unnecessary computations in 3D plots */
|
||||
|
||||
|
@ -85,6 +93,9 @@
|
|||
|
||||
#define COMPUTE_TOTAL_ENERGY ((ZPLOT == P_3D_TOTAL_ENERGY)||(CPLOT == P_3D_TOTAL_ENERGY)||(ZPLOT == P_3D_LOG_TOTAL_ENERGY)||(CPLOT == P_3D_LOG_TOTAL_ENERGY)||(ZPLOT == P_3D_MEAN_ENERGY)||(CPLOT == P_3D_MEAN_ENERGY)||(ZPLOT == P_3D_LOG_MEAN_ENERGY)||(CPLOT == P_3D_LOG_MEAN_ENERGY)||(ZPLOT_B == P_3D_TOTAL_ENERGY)||(CPLOT_B == P_3D_TOTAL_ENERGY)||(ZPLOT_B == P_3D_LOG_TOTAL_ENERGY)||(CPLOT_B == P_3D_LOG_TOTAL_ENERGY)||(ZPLOT_B == P_3D_MEAN_ENERGY)||(CPLOT_B == P_3D_MEAN_ENERGY)||(ZPLOT_B == P_3D_LOG_MEAN_ENERGY)||(CPLOT_B == P_3D_LOG_MEAN_ENERGY))
|
||||
|
||||
#define PLANET ((B_DOMAIN == D_SPHERE_EARTH)||(B_DOMAIN == D_SPHERE_MARS)||(B_DOMAIN == D_SPHERE_MOON))
|
||||
#define OTHER_PLANET ((B_DOMAIN == D_SPHERE_MARS)||(B_DOMAIN == D_SPHERE_MOON))
|
||||
|
||||
#define NMAXCIRC_SPHERE 100 /* max number of circles on sphere */
|
||||
|
||||
int global_time = 0;
|
||||
|
@ -154,6 +165,9 @@ typedef struct
|
|||
double r, g, b; /* RGB values for image */
|
||||
short int indomain; /* has value 1 if lattice point is in domain */
|
||||
double x2d, y2d; /* x and y coordinates for 2D representation */
|
||||
double altitude; /* altitude in case of Earth with digital elevation model */
|
||||
double cos_angle; /* cosine of light angle */
|
||||
double cos_angle_sphere; /* cosing of light angle for perfect sphere */
|
||||
} t_wave_sphere;
|
||||
|
||||
|
||||
|
|
|
@ -75,8 +75,10 @@
|
|||
#define S_BIN_OPENING 20 /* bin containing particles opening at deactivation time */
|
||||
#define S_POLYGON_EXT 21 /* exterior of a regular polygon */
|
||||
#define S_WEDGE_EXT 22 /* exterior of a wedge */
|
||||
#define S_MIXER 23 /* exterior of a set of fan of rectangles */
|
||||
#define S_MIXER 23 /* exterior of a blender made of rectangles */
|
||||
#define S_AIRFOIL 24 /* exterior of an air foil */
|
||||
#define S_COANDA 25 /* wall for Coanda effect */
|
||||
#define S_COANDA_SHORT 26 /* shorter wall for Coanda effect */
|
||||
|
||||
/* particle interaction */
|
||||
|
||||
|
@ -92,6 +94,7 @@
|
|||
#define I_VICSEK_REPULSIVE 9 /* Vicsek-type interaction with harmonic repulsion */
|
||||
#define I_VICSEK_SPEED 10 /* Vicsek-type interaction with speed adjustment */
|
||||
#define I_VICSEK_SHARK 11 /* Vicsek-type interaction with speed adjustment, and one shark */
|
||||
#define I_COULOMB_LJ 12 /* Coulomb force regularised by Lennard-Jones repulsion */
|
||||
|
||||
/* Boundary conditions */
|
||||
|
||||
|
@ -178,6 +181,17 @@
|
|||
#define IC_TWOROCKETS 8 /* type 1 or 2 depending on rocket position */
|
||||
#define IC_TWOROCKETS_TWOFUELS 9 /* type 1 and 2 or 1 and 3 depending on rocket */
|
||||
|
||||
/* Initial conditions for option TWO_TYPES */
|
||||
|
||||
#define TTC_RANDOM 0 /* assign types randomly */
|
||||
#define TTC_CHESSBOARD 1 /* assign types according to chessboard, works with hex initial config */
|
||||
#define TTC_COANDA 2 /* type 1 in a band of width LAMBDA */
|
||||
|
||||
/* Initial speed distribution */
|
||||
|
||||
#define VI_RANDOM 0 /* random (Gaussian) initial speed distribution */
|
||||
#define VI_COANDA 1 /* nonzero speed in a band of width LAMBDA */
|
||||
|
||||
/* Plot types */
|
||||
|
||||
#define P_KINETIC 0 /* colors represent kinetic energy of particles */
|
||||
|
@ -197,6 +211,11 @@
|
|||
#define P_DIRECT_EMEAN 14 /* averaged version of P_DIRECT_ENERGY */
|
||||
#define P_NOPARTICLE 15 /* particles are not drawn (only the links between them) */
|
||||
|
||||
/* Rotation schedules */
|
||||
|
||||
#define ROT_SPEEDUP_SLOWDOWN 0 /* rotation speeds up and then slows down to zero */
|
||||
#define ROT_BACK_FORTH 1 /* rotation goes in one direction and then back */
|
||||
|
||||
/* Initial position dependence types */
|
||||
|
||||
#define IP_X 0 /* color depends on x coordinate of initial position */
|
||||
|
@ -241,6 +260,7 @@ typedef struct
|
|||
double vy; /* y velocity of particle */
|
||||
double omega; /* angular velocity of particle */
|
||||
double mass_inv; /* inverse of particle mass */
|
||||
double charge; /* electric charge */
|
||||
double inertia_moment_inv; /* inverse of moment of inertia */
|
||||
double fx; /* x component of force on particle */
|
||||
double fy; /* y component of force on particle */
|
||||
|
@ -313,6 +333,8 @@ typedef struct
|
|||
double angle01, angle02; /* initial values of angles in which concave corners repel */
|
||||
double fx, fy; /* x and y-components of force on segment */
|
||||
double torque; /* torque on segment with respect to its center */
|
||||
double pressure; /* pressure acting on segement */
|
||||
double avrg_pressure; /* time-averaged pressure */
|
||||
short int inactivate; /* set to 1 for segment to become inactive at time SEGMENT_DEACTIVATION_TIME */
|
||||
} t_segment;
|
||||
|
||||
|
@ -361,6 +383,7 @@ typedef struct
|
|||
double angle; /* orientation of obstacle */
|
||||
double omega; /* angular speed of obstacle */
|
||||
double bdry_fx, bdry_fy; /* components of boundary force */
|
||||
double efield, bfield; /* electric and magnetic field */
|
||||
} t_lj_parameters;
|
||||
|
||||
int ncircles, nobstacles, nsegments, ngroups = 1, counter = 0;
|
||||
|
|
|
@ -77,6 +77,7 @@
|
|||
#define D_ONE_FUNNEL 581 /* one funnel */
|
||||
#define D_LENSES_RING 59 /* several lenses forming a ring */
|
||||
#define D_MAZE_CIRCULAR 60 /* circular maze */
|
||||
#define D_LENS 61 /* symmetric lens made of circular faces */
|
||||
|
||||
#define D_WING 70 /* complement of wing-shaped domain */
|
||||
#define D_TESLA 71 /* Tesla valve */
|
||||
|
@ -89,6 +90,9 @@
|
|||
#define D_SPHERE_JULIA 82 /* Julia set on Riemann sphere */
|
||||
#define D_SPHERE_JULIA_INV 83 /* inverted Julia set on Riemann sphere */
|
||||
#define D_SPHERE_EARTH 84 /* map of the Earth */
|
||||
#define D_SPHERE_JULIA_CUBIC 85 /* Julia set for cubic polynomial on Riemann sphere */
|
||||
#define D_SPHERE_MARS 86 /* map of Mars */
|
||||
#define D_SPHERE_MOON 87 /* map of the Moon */
|
||||
|
||||
#define NMAXCIRCLES 10000 /* total number of circles/polygons (must be at least NCX*NCY for square grid) */
|
||||
#define NMAXPOLY 50000 /* maximal number of vertices of polygonal lines (for von Koch et al) */
|
||||
|
@ -149,6 +153,8 @@
|
|||
#define IOR_POISSON_WELLS 8 /* wells located on a random Poisson disc process */
|
||||
#define IOR_PPP_WELLS 9 /* wells located on a Poisson point process */
|
||||
|
||||
#define IOR_EARTH_DEM 10 /* digital elevation model (for waves on sphere) */
|
||||
|
||||
/* Boundary conditions */
|
||||
|
||||
#define BC_DIRICHLET 0 /* Dirichlet boundary conditions */
|
||||
|
@ -167,6 +173,18 @@
|
|||
#define OSC_WAVE_PACKET 2 /* Gaussian wave packet */
|
||||
#define OSC_CHIRP 3 /* chirp (linearly accelerating frequency) */
|
||||
|
||||
/* Wave packet types */
|
||||
|
||||
#define WP_RANDOM1 0 /* random, variant 1 */
|
||||
#define WP_RANDOM2 1 /* random, variant 2 */
|
||||
#define WP_PAIR 2 /* 2 sources */
|
||||
|
||||
/* Wave packet envelope types */
|
||||
|
||||
#define WE_SINE 0 /* sine envelope */
|
||||
#define WE_CUTOFF 1 /* sharp cut-off */
|
||||
#define WE_CONSTANT 2 /* constant envelope */
|
||||
|
||||
/* For debugging purposes only */
|
||||
// #define FLOOR 0 /* set to 1 to limit wave amplitude to VMAX */
|
||||
// #define VMAX 10.0 /* max value of wave amplitude */
|
||||
|
|
1
heat.c
1
heat.c
|
@ -212,6 +212,7 @@
|
|||
#define WAVE_PACKET_SOURCE_TYPE 1 /* type of wave packet sources */
|
||||
#define N_WAVE_PACKETS 15 /* number of wave packets */
|
||||
#define OSCIL_LEFT_YSHIFT 0.0 /* y-dependence of left oscillation (for non-horizontal waves) */
|
||||
#define OSCILLATING_SOURCE_PERIOD 20 /* period of oscillating source */
|
||||
/* end of constants only used by sub_wave and sub_maze */
|
||||
|
||||
#include "global_pdes.c"
|
||||
|
|
178
lennardjones.c
178
lennardjones.c
|
@ -22,7 +22,7 @@
|
|||
/* It may be possible to increase parameter PAUSE */
|
||||
/* */
|
||||
/* create movie using */
|
||||
/* ffmpeg -i lj.%05d.tif -vcodec libx264 lj.mp4 */
|
||||
/* ffmpeg -i lj.%05d.tif -vcodec libx264 lj.mp4 */
|
||||
/* */
|
||||
/*********************************************************************************/
|
||||
|
||||
|
@ -50,8 +50,8 @@
|
|||
|
||||
/* General geometrical parameters */
|
||||
|
||||
#define WINWIDTH 1280 /* window width */
|
||||
#define WINHEIGHT 720 /* window height */
|
||||
#define WINWIDTH 1600 /* window width */
|
||||
#define WINHEIGHT 900 /* window height */
|
||||
|
||||
#define XMIN -2.0
|
||||
#define XMAX 2.0 /* x interval */
|
||||
|
@ -59,7 +59,7 @@
|
|||
#define YMAX 1.125 /* y interval for 9/16 aspect ratio */
|
||||
|
||||
#define INITXMIN -2.0
|
||||
#define INITXMAX 2.0 /* x interval for initial condition */
|
||||
#define INITXMAX 2.03 /* x interval for initial condition */
|
||||
#define INITYMIN -1.1
|
||||
#define INITYMAX 1.1 /* y interval for initial condition */
|
||||
|
||||
|
@ -86,14 +86,15 @@
|
|||
#define OBSTACLE_PATTERN 4 /* pattern of obstacles, see list in global_ljones.c */
|
||||
|
||||
#define ADD_FIXED_SEGMENTS 1 /* set to 1 to add fixed segments as obstacles */
|
||||
#define SEGMENT_PATTERN 24 /* pattern of repelling segments, see list in global_ljones.c */
|
||||
#define SEGMENT_PATTERN 25 /* pattern of repelling segments, see list in global_ljones.c */
|
||||
#define ROCKET_SHAPE 3 /* shape of rocket combustion chamber, see list in global_ljones.c */
|
||||
#define ROCKET_SHAPE_B 3 /* shape of second rocket */
|
||||
#define NOZZLE_SHAPE 6 /* shape of nozzle, see list in global_ljones.c */
|
||||
#define NOZZLE_SHAPE_B 6 /* shape of nozzle for second rocket, see list in global_ljones.c */
|
||||
|
||||
#define TWO_TYPES 0 /* set to 1 to have two types of particles */
|
||||
#define TWO_TYPES 1 /* set to 1 to have two types of particles */
|
||||
#define TYPE_PROPORTION 0.5 /* proportion of particles of first type */
|
||||
#define TWOTYPE_CONFIG 2 /* choice of types, see TTC_ list in global_ljones.c */
|
||||
#define SYMMETRIZE_FORCE 1 /* set to 1 to symmetrize two-particle interaction, only needed if particles are not all the same */
|
||||
#define CENTER_PX 0 /* set to 1 to center horizontal momentum */
|
||||
#define CENTER_PY 0 /* set to 1 to center vertical momentum */
|
||||
|
@ -110,9 +111,9 @@
|
|||
#define PDISC_CANDIDATES 100 /* number of candidates in construction of Poisson disc process */
|
||||
#define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */
|
||||
|
||||
#define LAMBDA 0.75 /* parameter controlling the dimensions of domain */
|
||||
#define MU 0.012 /* parameter controlling radius of particles */
|
||||
#define MU_B 0.010 /* parameter controlling radius of particles of second type */
|
||||
#define LAMBDA 0.3 /* parameter controlling the dimensions of domain */
|
||||
#define MU 0.014 /* parameter controlling radius of particles */
|
||||
#define MU_B 0.01 /* parameter controlling radius of particles of second type */
|
||||
#define NPOLY 40 /* number of sides of polygon */
|
||||
#define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */
|
||||
#define AWEDGE 0.5 /* opening angle of wedge, in units of Pi/2 */
|
||||
|
@ -121,10 +122,8 @@
|
|||
#define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */
|
||||
#define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */
|
||||
#define FOCI 1 /* set to 1 to draw focal points of ellipse */
|
||||
// #define NGRIDX 140 /* number of grid point for grid of disks */
|
||||
// #define NGRIDY 70 /* number of grid point for grid of disks */
|
||||
#define NGRIDX 130 /* number of grid point for grid of disks */
|
||||
#define NGRIDY 65 /* number of grid point for grid of disks */
|
||||
#define NGRIDX 101 /* number of grid point for grid of disks */
|
||||
#define NGRIDY 47 /* number of grid point for grid of disks */
|
||||
#define EHRENFEST_RADIUS 0.9 /* radius of container for Ehrenfest urn configuration */
|
||||
#define EHRENFEST_WIDTH 0.035 /* width of tube for Ehrenfest urn configuration */
|
||||
#define TWO_CIRCLES_RADIUS_RATIO 0.8 /* ratio of radii for S_TWO_CIRCLES_EXT segment configuration */
|
||||
|
@ -137,9 +136,8 @@
|
|||
|
||||
/* Parameters for length and speed of simulation */
|
||||
|
||||
#define NSTEPS 4525 /* number of frames of movie */
|
||||
// #define NSTEPS 1000 /* number of frames of movie */
|
||||
#define NVID 30 /* number of iterations between images displayed on screen */
|
||||
#define NSTEPS 3200 /* number of frames of movie */
|
||||
#define NVID 15 /* number of iterations between images displayed on screen */
|
||||
#define NSEG 25 /* number of segments of boundary of circles */
|
||||
#define INITIAL_TIME 0 /* time after which to start saving frames */
|
||||
#define OBSTACLE_INITIAL_TIME 150 /* time after which to start moving obstacle */
|
||||
|
@ -152,7 +150,7 @@
|
|||
#define SLEEP1 1 /* initial sleeping time */
|
||||
#define SLEEP2 1 /* final sleeping time */
|
||||
#define MID_FRAMES 20 /* number of still frames between parts of two-part movie */
|
||||
#define END_FRAMES 100 /* number of still frames at end of movie */
|
||||
#define END_FRAMES 250 /* number of still frames at end of movie */
|
||||
|
||||
/* Boundary conditions, see list in global_ljones.c */
|
||||
|
||||
|
@ -160,18 +158,18 @@
|
|||
|
||||
/* Plot type, see list in global_ljones.c */
|
||||
|
||||
#define PLOT 13
|
||||
#define PLOT 5
|
||||
#define PLOT_B 11 /* plot type for second movie */
|
||||
|
||||
#define DRAW_BONDS 1 /* set to 1 to draw bonds between neighbours */
|
||||
#define COLOR_BONDS 1 /* set to 1 to color bonds according to length */
|
||||
#define FILL_TRIANGLES 0 /* set to 1 to fill triangles between neighbours */
|
||||
#define ALTITUDE_LINES 0 /* set to 1 to add horizontal lines to show altitude */
|
||||
#define COLOR_SEG_GROUPS 1 /* set to 1 to collor segment groups differently */
|
||||
#define COLOR_SEG_GROUPS 0 /* set to 1 to collor segment groups differently */
|
||||
#define N_PARTICLE_COLORS 200 /* number of colors for P_NUMBER color scheme */
|
||||
#define INITIAL_POS_TYPE 0 /* type of initial position dependence */
|
||||
#define ERATIO 0.995 /* ratio for time-averagin in P_EMEAN color scheme */
|
||||
#define DRATIO 0.995 /* ratio for time-averagin in P_DIRECT_EMEAN color scheme */
|
||||
#define ERATIO 0.995 /* ratio for time-averaging in P_EMEAN color scheme */
|
||||
#define DRATIO 0.995 /* ratio for time-averaging in P_DIRECT_EMEAN color scheme */
|
||||
|
||||
/* Color schemes */
|
||||
|
||||
|
@ -195,15 +193,15 @@
|
|||
#define LUMAMP 0.3 /* amplitude of luminosity variation for scheme C_LUM */
|
||||
#define HUEMEAN 220.0 /* mean value of hue for color scheme C_HUE */
|
||||
#define HUEAMP -50.0 /* amplitude of variation of hue for color scheme C_HUE */
|
||||
#define COLOR_HUESHIFT 0.5 /* shift in color hue (for some cyclic palettes) */
|
||||
#define COLOR_HUESHIFT 0.0 /* shift in color hue (for some cyclic palettes) */
|
||||
|
||||
#define PRINT_PARAMETERS 1 /* set to 1 to print certain parameters */
|
||||
#define PRINT_PARAMETERS 0 /* set to 1 to print certain parameters */
|
||||
#define PRINT_TEMPERATURE 0 /* set to 1 to print current temperature */
|
||||
#define PRINT_ANGLE 1 /* set to 1 to print obstacle orientation */
|
||||
#define PRINT_OMEGA 1 /* set to 1 to print angular speed */
|
||||
#define PRINT_ANGLE 0 /* set to 1 to print obstacle orientation */
|
||||
#define PRINT_OMEGA 0 /* set to 1 to print angular speed */
|
||||
#define PRINT_PARTICLE_SPEEDS 0 /* set to 1 to print average speeds/momenta of particles */
|
||||
#define PRINT_SEGMENTS_SPEEDS 0 /* set to 1 to print velocity of moving segments */
|
||||
#define PRINT_SEGMENTS_FORCE 1 /* set to 1 to print force on segments */
|
||||
#define PRINT_SEGMENTS_FORCE 0 /* set to 1 to print force on segments */
|
||||
#define FORCE_FACTOR 0.1 /* factor controlling length of force vector */
|
||||
|
||||
/* particle properties */
|
||||
|
@ -212,7 +210,7 @@
|
|||
#define ENERGY_HUE_MAX 50.0 /* color of saturated particle */
|
||||
#define PARTICLE_HUE_MIN 359.0 /* color of original particle */
|
||||
#define PARTICLE_HUE_MAX 0.0 /* color of saturated particle */
|
||||
#define PARTICLE_EMAX 5000.0 /* energy of particle with hottest color */
|
||||
#define PARTICLE_EMAX 1000.0 /* energy of particle with hottest color */
|
||||
#define HUE_TYPE0 280.0 /* hue of particles of type 0 */
|
||||
#define HUE_TYPE1 60.0 /* hue of particles of type 1 */
|
||||
#define HUE_TYPE2 140.0 /* hue of particles of type 2 */
|
||||
|
@ -221,35 +219,35 @@
|
|||
#define RANDOM_RADIUS 0 /* set to 1 for random circle radius */
|
||||
#define DT_PARTICLE 3.0e-6 /* time step for particle displacement */
|
||||
#define KREPEL 1.0 /* constant in repelling force between particles */
|
||||
#define EQUILIBRIUM_DIST 5.0 /* Lennard-Jones equilibrium distance */
|
||||
#define EQUILIBRIUM_DIST_B 5.0 /* Lennard-Jones equilibrium distance for second type of particle */
|
||||
#define EQUILIBRIUM_DIST 7.0 /* Lennard-Jones equilibrium distance */
|
||||
#define EQUILIBRIUM_DIST_B 6.5 /* Lennard-Jones equilibrium distance for second type of particle */
|
||||
#define REPEL_RADIUS 15.0 /* radius in which repelling force acts (in units of particle radius) */
|
||||
#define DAMPING 100.0 /* damping coefficient of particles */
|
||||
#define DAMPING 0.02 /* damping coefficient of particles */
|
||||
#define OMEGA_INITIAL 10.0 /* initial angular velocity range */
|
||||
#define INITIAL_DAMPING 5.0 /* damping coefficient of particles during initial phase */
|
||||
#define DAMPING_ROT 100.0 /* dampint coefficient for rotation of particles */
|
||||
#define PARTICLE_MASS 1.0 /* mass of particle of radius MU */
|
||||
#define PARTICLE_MASS_B 1.0 /* mass of particle of radius MU */
|
||||
#define PARTICLE_MASS_B 6.5 /* mass of particle of radius MU */
|
||||
#define PARTICLE_INERTIA_MOMENT 0.2 /* moment of inertia of particle */
|
||||
#define PARTICLE_INERTIA_MOMENT_B 0.02 /* moment of inertia of second type of particle */
|
||||
#define V_INITIAL 0.0 /* initial velocity range */
|
||||
#define V_INITIAL 1000.0 /* initial velocity range */
|
||||
#define OMEGA_INITIAL 10.0 /* initial angular velocity range */
|
||||
#define VICSEK_VMIN 1.0 /* minimal speed of particles in Vicsek model */
|
||||
#define VICSEK_VMAX 40.0 /* minimal speed of particles in Vicsek model */
|
||||
|
||||
#define V_INITIAL_TYPE 1 /* type of initial speed distribution (see VI_ in global_ljones.c) */
|
||||
|
||||
#define THERMOSTAT 0 /* set to 1 to switch on thermostat */
|
||||
#define VARY_THERMOSTAT 0 /* set to 1 for time-dependent thermostat schedule */
|
||||
#define SIGMA 5.0 /* noise intensity in thermostat */
|
||||
#define BETA 0.02 /* initial inverse temperature */
|
||||
#define BETA 0.0002 /* initial inverse temperature */
|
||||
#define MU_XI 0.01 /* friction constant in thermostat */
|
||||
#define KSPRING_BOUNDARY 1.0e7 /* confining harmonic potential outside simulation region */
|
||||
#define KSPRING_OBSTACLE 1.0e11 /* harmonic potential of obstacles */
|
||||
// #define NBH_DIST_FACTOR 2.3 /* radius in which to count neighbours */
|
||||
// #define NBH_DIST_FACTOR 2.7 /* radius in which to count neighbours */
|
||||
#define NBH_DIST_FACTOR 3.8 /* radius in which to count neighbours */
|
||||
// #define NBH_DIST_FACTOR 4.0 /* radius in which to count neighbours */
|
||||
#define GRAVITY 0.0 /* gravity acting on all particles */
|
||||
#define GRAVITY_X 5000.0 /* horizontal gravity acting on all particles */
|
||||
#define GRAVITY 0.0 /* gravity acting on all particles */
|
||||
#define GRAVITY_X 0.0 /* horizontal gravity acting on all particles */
|
||||
#define CIRCULAR_GRAVITY 0 /* set to 1 to have gravity directed to center */
|
||||
#define INCREASE_GRAVITY 0 /* set to 1 to increase gravity during the simulation */
|
||||
#define GRAVITY_SCHEDULE 2 /* type of gravity schedule, see list in global_ljones.c */
|
||||
#define GRAVITY_FACTOR 100.0 /* factor by which to increase gravity */
|
||||
|
@ -258,6 +256,15 @@
|
|||
#define KSPRING_VICSEK 0.2 /* spring constant for I_VICSEK_SPEED interaction */
|
||||
#define VICSEK_REPULSION 10.0 /* repulsion between particles in Vicsek model */
|
||||
|
||||
#define ADD_EFIELD 1 /* set to 1 to add an electric field */
|
||||
#define EFIELD 10000.0 /* value of electric field */
|
||||
#define ADD_BFIELD 0 /* set to 1 to add a magnetic field */
|
||||
#define BFIELD 1000.0 /* value of magnetic field */
|
||||
#define CHARGE 0.0 /* charge of particles of first type */
|
||||
#define CHARGE_B 1.0 /* charge of particles of second type */
|
||||
#define INCREASE_E 0 /* set to 1 to increase electric field */
|
||||
#define EFIELD_FACTOR 1000.0 /* factor by which to increase electric field */
|
||||
|
||||
#define ROTATION 0 /* set to 1 to include rotation of particles */
|
||||
#define COUPLE_ANGLE_TO_THERMOSTAT 1 /* set to 1 to couple angular degrees of freedom to thermostat */
|
||||
#define DIMENSION_FACTOR 1.0 /* scaling factor taking into account number of degrees of freedom */
|
||||
|
@ -267,7 +274,7 @@
|
|||
#define KTORQUE_DIFF 150.0 /* force constant in angular dynamics for different particles */
|
||||
#define DRAW_SPIN 0 /* set to 1 to draw spin vectors of particles */
|
||||
#define DRAW_SPIN_B 0 /* set to 1 to draw spin vectors of particles */
|
||||
#define DRAW_CROSS 1 /* set to 1 to draw cross on particles of second type */
|
||||
#define DRAW_CROSS 0 /* set to 1 to draw cross on particles of second type */
|
||||
#define SPIN_RANGE 10.0 /* range of spin-spin interaction */
|
||||
#define SPIN_RANGE_B 5.0 /* range of spin-spin interaction for second type of particle */
|
||||
#define QUADRUPOLE_RATIO 0.6 /* anisotropy in quadrupole potential */
|
||||
|
@ -290,7 +297,7 @@
|
|||
#define CENTER_VIEW_ON_OBSTACLE 0 /* set to 1 to center display on moving obstacle */
|
||||
#define RESAMPLE_Y 0 /* set to 1 to resample y coordinate of moved particles (for shock waves) */
|
||||
#define NTRIALS 2000 /* number of trials when resampling */
|
||||
#define OBSTACLE_RADIUS 0.25 /* radius of obstacle for circle boundary conditions */
|
||||
#define OBSTACLE_RADIUS 0.15 /* radius of obstacle for circle boundary conditions */
|
||||
#define FUNNEL_WIDTH 0.25 /* funnel width for funnel boundary conditions */
|
||||
#define OBSTACLE_XMIN 0.0 /* initial position of obstacle */
|
||||
#define OBSTACLE_XMAX 3.0 /* final position of obstacle */
|
||||
|
@ -325,24 +332,28 @@
|
|||
#define TRACER_PARTICLE_MASS 4.0 /* relative mass of tracer particle */
|
||||
#define TRAJECTORY_WIDTH 3 /* width of tracer particle trajectory */
|
||||
|
||||
#define ROTATE_BOUNDARY 1 /* set to 1 to rotate the repelling segments */
|
||||
#define ROTATE_BOUNDARY 0 /* set to 1 to rotate the repelling segments */
|
||||
#define SMOOTH_ROTATION 1 /* set to 1 to update segments at each time step (rather than at each movie frame) */
|
||||
#define ROTATION_SCHEDULE 0 /* time-dependence of rotation angle, see ROT_* in global_ljones.c */
|
||||
#define PERIOD_ROTATE_BOUNDARY 1000 /* period of rotating boundary */
|
||||
#define ROTATE_INITIAL_TIME 300 /* initial time without rotation */
|
||||
#define ROTATE_FINAL_TIME 300 /* final time without rotation */
|
||||
#define ROTATE_CHANGE_TIME 0.5 /* relative duration of acceleration/deceleration phases */
|
||||
#define OMEGAMAX -2.0*PI /* maximal rotation speed */
|
||||
#define OMEGAMAX -2.0*PI /* maximal rotation speed */
|
||||
|
||||
#define MOVE_BOUNDARY 0 /* set to 1 to move repelling segments, due to force from particles */
|
||||
#define SEGMENTS_MASS 40.0 /* mass of collection of segments */
|
||||
#define DEACTIVATE_SEGMENT 1 /* set to 1 to deactivate last segment after a certain time */
|
||||
#define SEGMENT_DEACTIVATION_TIME 20 /* time at which to deactivate last segment */
|
||||
#define RELEASE_ROCKET_AT_DEACTIVATION 0 /* set to 1 to limit segments velocity before segment release */
|
||||
#define SEGMENTS_X0 1.5 /* initial position of segments */
|
||||
#define SEGMENTS_Y0 0.0 /* initial position of segments */
|
||||
#define SEGMENTS_X0 0.0 /* initial position of segments */
|
||||
#define SEGMENTS_Y0 -0.6 /* initial position of segments */
|
||||
#define SEGMENTS_VX0 0.0 /* initial velocity of segments */
|
||||
#define SEGMENTS_VY0 0.0 /* initial velocity of segments */
|
||||
#define DAMP_SEGS_AT_NEGATIVE_Y 0 /* set to 1 to dampen segments when y coordinate is negative */
|
||||
#define SHOW_SEGMENTS_PRESSURE 0 /* set to 1 to show (averaged) pressure acting on segments */
|
||||
#define SEGMENT_PMAX 7.5e7 /* pressure of segment with hottest color */
|
||||
#define P_AVRG_FACTOR 0.02 /* factor in computation of mean pressure */
|
||||
|
||||
#define MOVE_SEGMENT_GROUPS 0 /* set to 1 to group segments into moving units */
|
||||
#define SEGMENT_GROUP_MASS 500.0 /* mass of segment group */
|
||||
|
@ -407,8 +418,8 @@
|
|||
#define FLOOR_OMEGA 0 /* set to 1 to limit particle momentum to PMAX */
|
||||
#define PMAX 1000.0 /* maximal force */
|
||||
|
||||
#define HASHX 120 /* size of hashgrid in x direction */
|
||||
#define HASHY 60 /* size of hashgrid in y direction */
|
||||
#define HASHX 60 /* size of hashgrid in x direction */
|
||||
#define HASHY 30 /* size of hashgrid in y direction */
|
||||
#define HASHMAX 100 /* maximal number of particles per hashgrid cell */
|
||||
#define HASHGRID_PADDING 0.1 /* padding of hashgrid outside simulation window */
|
||||
|
||||
|
@ -473,6 +484,23 @@ double repel_schedule(int i)
|
|||
}
|
||||
|
||||
|
||||
double efield_schedule(int i)
|
||||
{
|
||||
static double efactor;
|
||||
static int first = 1;
|
||||
double efield;
|
||||
|
||||
if (first)
|
||||
{
|
||||
efactor = EFIELD_FACTOR/(double)(INITIAL_TIME + NSTEPS);
|
||||
first = 0;
|
||||
}
|
||||
efield = EFIELD*(double)i*efactor;
|
||||
printf("E = %.3lg\n", efield);
|
||||
return(efield);
|
||||
}
|
||||
|
||||
|
||||
double temperature_schedule(int i)
|
||||
{
|
||||
static double bexponent, omega, bexp2;
|
||||
|
@ -614,20 +642,26 @@ double rotation_angle(double phase)
|
|||
// return(0.5*OMEGAMAX*(phase - sin(phase)));
|
||||
|
||||
/* case of centrifuge remaining at constant speed for a while */
|
||||
if (phase < ROTATE_CHANGE_TIME)
|
||||
{
|
||||
// angular_speed = 0.5*OMEGAMAX*(1.0 - cos(phase*PI/ROTATE_CHANGE_TIME));
|
||||
return(0.5*OMEGAMAX*(phase - (ROTATE_CHANGE_TIME/PI)*sin(phase*PI/ROTATE_CHANGE_TIME)));
|
||||
}
|
||||
else if (phase < 1.0 - ROTATE_CHANGE_TIME)
|
||||
{
|
||||
// angular_speed = OMEGAMAX;
|
||||
return(0.5*OMEGAMAX*(2.0*phase - ROTATE_CHANGE_TIME));
|
||||
}
|
||||
else
|
||||
{
|
||||
// angular_speed = 0.5*OMEGAMAX*(1.0 + cos((phase - 1.0 + ROTATE_CHANGE_TIME)*PI/ROTATE_CHANGE_TIME));
|
||||
return(0.5*OMEGAMAX*(2.0 - 2.0*ROTATE_CHANGE_TIME + phase - 1.0 + (ROTATE_CHANGE_TIME/PI)*sin((1.0-phase)*PI/ROTATE_CHANGE_TIME)));
|
||||
switch (ROTATION_SCHEDULE) {
|
||||
case (ROT_SPEEDUP_SLOWDOWN):
|
||||
{
|
||||
if (phase < ROTATE_CHANGE_TIME)
|
||||
{
|
||||
return(0.5*OMEGAMAX*(phase - (ROTATE_CHANGE_TIME/PI)*sin(phase*PI/ROTATE_CHANGE_TIME)));
|
||||
}
|
||||
else if (phase < 1.0 - ROTATE_CHANGE_TIME)
|
||||
{
|
||||
return(0.5*OMEGAMAX*(2.0*phase - ROTATE_CHANGE_TIME));
|
||||
}
|
||||
else
|
||||
{
|
||||
return(0.5*OMEGAMAX*(2.0 - 2.0*ROTATE_CHANGE_TIME + phase - 1.0 + (ROTATE_CHANGE_TIME/PI)*sin((1.0-phase)*PI/ROTATE_CHANGE_TIME)));
|
||||
}
|
||||
}
|
||||
case (ROT_BACK_FORTH):
|
||||
{
|
||||
return(OMEGAMAX*(1.0 - cos(DPI*phase))/DPI);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
@ -1252,7 +1286,8 @@ void animation()
|
|||
printf("Computing frame %d\n",i);
|
||||
|
||||
if (INCREASE_KREPEL) params.krepel = repel_schedule(i);
|
||||
if (INCREASE_BETA) params.beta = temperature_schedule(i);
|
||||
if (INCREASE_BETA) params.beta = temperature_schedule(i);
|
||||
if (INCREASE_E) params.efield = efield_schedule(i);
|
||||
if (DECREASE_CONTAINER_SIZE)
|
||||
{
|
||||
params.xmincontainer = container_size_schedule(i);
|
||||
|
@ -1358,12 +1393,33 @@ void animation()
|
|||
|
||||
/* add gravity */
|
||||
if (INCREASE_GRAVITY) particle[j].fy -= params.gravity/particle[j].mass_inv;
|
||||
else if (CIRCULAR_GRAVITY)
|
||||
{
|
||||
particle[j].fx -= GRAVITY*particle[j].xc/particle[j].mass_inv;
|
||||
particle[j].fy -= GRAVITY*particle[j].yc/particle[j].mass_inv;
|
||||
}
|
||||
else
|
||||
{
|
||||
particle[j].fy -= GRAVITY/particle[j].mass_inv;
|
||||
particle[j].fx += GRAVITY_X/particle[j].mass_inv;
|
||||
}
|
||||
|
||||
/* add electric force */
|
||||
if (ADD_EFIELD)
|
||||
{
|
||||
if (INCREASE_E) particle[j].fx += params.efield*particle[j].charge;
|
||||
else particle[j].fx += EFIELD*particle[j].charge;
|
||||
}
|
||||
|
||||
/* add magnetic force */
|
||||
if (ADD_BFIELD)
|
||||
{
|
||||
// particle[j].fx += BFIELD*particle[j].charge*particle[j].vy;
|
||||
// particle[j].fy -= BFIELD*particle[j].charge*particle[j].vx;
|
||||
particle[j].fx += BFIELD*particle[j].charge*particle[j].vy*particle[j].mass_inv;
|
||||
particle[j].fy -= BFIELD*particle[j].charge*particle[j].vx*particle[j].mass_inv;
|
||||
}
|
||||
|
||||
if (FLOOR_FORCE)
|
||||
{
|
||||
if (particle[j].fx > FMAX) particle[j].fx = FMAX;
|
||||
|
|
2
makefile
2
makefile
|
@ -1,4 +1,4 @@
|
|||
LIBS = -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lXmu -lglut -fopenmp
|
||||
LIBS = -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lglut -fopenmp
|
||||
|
||||
.c:
|
||||
gcc -o $@ $< -O3 $(LIBS)
|
||||
|
|
|
@ -251,16 +251,11 @@
|
|||
#define IOR 7 /* choice of index of refraction, see list in global_pdes.c */
|
||||
#define IOR_TOTAL_TURNS 1.5 /* total angle of rotation for IOR_PERIODIC_WELLS_ROTATING */
|
||||
#define MANDEL_IOR_SCALE -0.05 /* parameter controlling dependence of IoR on Mandelbrot escape speed */
|
||||
#define COURANT 0.04 /* Courant number */
|
||||
#define COURANTB 0.0 /* Courant number in medium B */
|
||||
#define INITIAL_AMP 0.5 /* amplitude of initial condition */
|
||||
#define INITIAL_VARIANCE 0.0003 /* variance of initial condition */
|
||||
#define INITIAL_WAVELENGTH 0.015 /* wavelength of initial condition */
|
||||
#define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */
|
||||
#define WAVE_PACKET_SOURCE_TYPE 1 /* type of wave packet sources */
|
||||
#define N_WAVE_PACKETS 15 /* number of wave packets */
|
||||
#define OSCIL_LEFT_YSHIFT 0.0 /* y-dependence of left oscillation (for non-horizontal waves) */
|
||||
#define DRAW_WAVE_PROFILE 0 /* set to 1 to draw a profile of the wave */
|
||||
#define OSCILLATING_SOURCE_PERIOD 20 /* period of oscillating source */
|
||||
/* end of constants only used by sub_wave and sub_maze */
|
||||
|
||||
#include "global_pdes.c"
|
||||
|
|
Binary file not shown.
|
@ -192,6 +192,7 @@
|
|||
#define WAVE_PACKET_SOURCE_TYPE 1 /* type of wave packet sources */
|
||||
#define N_WAVE_PACKETS 15 /* number of wave packets */
|
||||
#define OSCIL_LEFT_YSHIFT 0.0 /* y-dependence of left oscillation (for non-horizontal waves) */
|
||||
#define OSCILLATING_SOURCE_PERIOD 20 /* period of oscillating source */
|
||||
/* end of constants only used by sub_wave and sub_maze */
|
||||
|
||||
|
||||
|
|
263
sub_lj.c
263
sub_lj.c
|
@ -34,8 +34,8 @@ int writetiff(char *filename, char *description, int x, int y, int width, int he
|
|||
glPixelStorei(GL_PACK_ALIGNMENT, 1);
|
||||
|
||||
glReadPixels(x, y, width, height, GL_RGB, GL_UNSIGNED_BYTE, image);
|
||||
TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32) width);
|
||||
TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32) height);
|
||||
TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32_t) width);
|
||||
TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32_t) height);
|
||||
TIFFSetField(file, TIFFTAG_BITSPERSAMPLE, 8);
|
||||
TIFFSetField(file, TIFFTAG_COMPRESSION, compression);
|
||||
TIFFSetField(file, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB);
|
||||
|
@ -622,15 +622,18 @@ void init_particle_config(t_particle particles[NMAXCIRCLES])
|
|||
case (C_SQUARE):
|
||||
{
|
||||
ncircles = NGRIDX*NGRIDY;
|
||||
dy = (YMAX - YMIN)/((double)NGRIDY);
|
||||
dx = (INITXMAX - INITXMIN)/((double)NGRIDX);
|
||||
dy = (INITYMAX - INITYMIN)/((double)NGRIDY);
|
||||
for (i = 0; i < NGRIDX; i++)
|
||||
for (j = 0; j < NGRIDY; j++)
|
||||
{
|
||||
n = NGRIDY*i + j;
|
||||
particles[n].xc = ((double)(i-NGRIDX/2) + 0.5)*dy;
|
||||
particles[n].yc = YMIN + ((double)j + 0.5)*dy;
|
||||
particles[n].xc = INITXMIN + ((double)i - 0.5)*dx;
|
||||
particles[n].yc = INITYMIN + ((double)j - 0.5)*dy;
|
||||
particles[n].radius = MU;
|
||||
particles[n].active = 1;
|
||||
/* activate only circles that intersect the domain */
|
||||
if ((particles[n].yc < INITYMAX + MU)&&(particles[n].yc > INITYMIN - MU)&&(particles[n].xc < INITXMAX + MU)&&(particles[n].xc > INITXMIN - MU)) particles[n].active = 1;
|
||||
else particles[n].active = 0;
|
||||
}
|
||||
break;
|
||||
}
|
||||
|
@ -1616,7 +1619,7 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
|
|||
/* initialise linear obstacle configuration */
|
||||
{
|
||||
int i, j, cycle = 0, iminus, iplus, nsides, n, concave = 1;
|
||||
double angle, angle2, dangle, dx, width, height, a, b, length, xmid = 0.5*(BCXMIN + BCXMAX), lpocket, r, x, x1, y1, x2, y2, nozx, nozy, y, dy, ca, sa;
|
||||
double angle, angle2, dangle, dx, width, height, a, b, length, xmid = 0.5*(BCXMIN + BCXMAX), lpocket, r, x, x1, y1, x2, y2, nozx, nozy, y, yy, dy, ca, sa, padding;
|
||||
|
||||
switch (SEGMENT_PATTERN) {
|
||||
case (S_RECTANGLE):
|
||||
|
@ -2399,6 +2402,82 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
|
|||
|
||||
break;
|
||||
}
|
||||
case (S_COANDA):
|
||||
{
|
||||
y1 = SEGMENTS_Y0;
|
||||
width = 0.05;
|
||||
padding = 0.01;
|
||||
|
||||
height = 0.25;
|
||||
dx = (XMAX - XMIN + 2.0*padding)/(double)(NPOLY-1);
|
||||
|
||||
for (i=0; i<NPOLY; i++)
|
||||
{
|
||||
x = XMIN - padding + (double)i*dx;
|
||||
y = y1 + height*cos(PI*x/XMAX);
|
||||
|
||||
segment[i].x1 = x;
|
||||
segment[i].y1 = y + width;
|
||||
segment[i].x2 = x + dx;
|
||||
segment[i].x2 = y1 + height*cos(PI*(x+dx)/XMAX) + width;
|
||||
|
||||
segment[i].concave = 1;
|
||||
}
|
||||
|
||||
for (i=0; i<NPOLY; i++)
|
||||
{
|
||||
x = XMAX + padding - (double)i*dx;
|
||||
y = y1 + height*cos(PI*x/XMAX);
|
||||
|
||||
segment[NPOLY+i].x1 = x;
|
||||
segment[NPOLY+i].y1 = y - width;
|
||||
segment[NPOLY+i].x2 = x - dx;
|
||||
segment[NPOLY+i].x2 = y1 + height*cos(PI*(x-dx)/XMAX) - width;
|
||||
|
||||
segment[NPOLY+i].concave = 1;
|
||||
}
|
||||
|
||||
nsegments = 2*NPOLY;
|
||||
|
||||
cycle = 1;
|
||||
concave = 1;
|
||||
|
||||
for (i=0; i<nsegments; i++)
|
||||
{
|
||||
segment[i].inactivate = 0;
|
||||
}
|
||||
|
||||
break;
|
||||
}
|
||||
case (S_COANDA_SHORT):
|
||||
{
|
||||
y1 = SEGMENTS_Y0;
|
||||
width = 0.05;
|
||||
padding = 0.1;
|
||||
x1 = XMIN + padding;
|
||||
|
||||
height = 0.2;
|
||||
length = 0.85*(XMAX - XMIN - 2.0*padding);
|
||||
dx = length/(double)(NPOLY-1);
|
||||
|
||||
for (i=0; i<NPOLY; i++)
|
||||
{
|
||||
x = x1 + (double)i*dx;
|
||||
y = y1 - height*cos(DPI*(x-x1)/length);
|
||||
yy = y1 - height*cos(DPI*(x+dx-x1)/length);
|
||||
add_rotated_angle_to_segments(x, y, x+dx, yy, width, 1, segment, 0);
|
||||
}
|
||||
|
||||
cycle = 0;
|
||||
concave = 1;
|
||||
|
||||
for (i=0; i<nsegments; i++)
|
||||
{
|
||||
segment[i].inactivate = 0;
|
||||
}
|
||||
|
||||
break;
|
||||
}
|
||||
default:
|
||||
{
|
||||
printf("Function init_segment_config not defined for this pattern \n");
|
||||
|
@ -2491,6 +2570,10 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
|
|||
segment[i].angle02 = segment[i].angle2;
|
||||
}
|
||||
|
||||
/* initialize averaged pressure */
|
||||
if (SHOW_SEGMENTS_PRESSURE) for (i=0; i<nsegments; i++)
|
||||
segment[i].avrg_pressure = 0.0;
|
||||
|
||||
// for (i=0; i<nsegments; i++)
|
||||
// {
|
||||
// printf("Segment %i: (x1, y1) = (%.3lg,%.3lg), (x2, y2) = (%.3lg,%.3lg)\n (nx, ny) = (%.3lg,%.3lg), c = %.3lg, length = %.3lg\n", i, segment[i].x1, segment[i].y1, segment[i].x2, segment[i].y2, segment[i].nx, segment[i].ny, segment[i].c, segment[i].length);
|
||||
|
@ -2550,7 +2633,7 @@ int in_segment_region(double x, double y, t_segment segment[NMAXSEGMENTS])
|
|||
/* returns 1 if (x,y) is inside region delimited by obstacle segments */
|
||||
{
|
||||
int i;
|
||||
double angle, dx, height, width, theta, lx, ly, x1, y1, x2, y2, padding, ca, sa, r;
|
||||
double angle, dx, height, width, length, theta, lx, ly, x1, y1, x2, y2, padding, ca, sa, r;
|
||||
|
||||
if (x >= BCXMAX) return(0);
|
||||
if (x <= BCXMIN) return(0);
|
||||
|
@ -2750,6 +2833,17 @@ int in_segment_region(double x, double y, t_segment segment[NMAXSEGMENTS])
|
|||
y1 += 0.5*x1*x1;
|
||||
return(x1*x1 + 25.0*y1*y1 > LAMBDA*LAMBDA*(1.0 + padding));
|
||||
}
|
||||
case (S_COANDA):
|
||||
{
|
||||
return(vabs(y - SEGMENTS_Y0 - 0.25*cos(PI*x/XMAX)) > 0.1);
|
||||
}
|
||||
case (S_COANDA_SHORT):
|
||||
{
|
||||
x1 = XMIN + 0.1;
|
||||
length = 0.85*(XMAX - XMIN - 0.2);
|
||||
if (x > x1 + length + 0.1) return(1);
|
||||
return(vabs(y - SEGMENTS_Y0 + 0.2*cos(DPI*(x-x1)/length)) > 0.05);
|
||||
}
|
||||
case (S_MAZE):
|
||||
{
|
||||
for (i=0; i<nsegments; i++)
|
||||
|
@ -3384,7 +3478,7 @@ int compute_particle_interaction(int i, int k, double force[2], double *torque,
|
|||
/* compute repelling force and torque of particle #k on particle #i */
|
||||
/* returns 1 if distance between particles is smaller than NBH_DIST_FACTOR*MU */
|
||||
{
|
||||
double x1, y1, x2, y2, r, f, angle, aniso, fx, fy, ff[2], dist_scaled, spin_f, ck, sk, ck_rel, sk_rel, alpha, amp;
|
||||
double x1, y1, x2, y2, r, f, angle, aniso, fx, fy, ff[2], dist_scaled, spin_f, ck, sk, ck_rel, sk_rel, alpha, amp, charge;
|
||||
static double dxhalf = 0.5*(BCXMAX - BCXMIN), dyhalf = 0.5*(BCYMAX - BCYMIN);
|
||||
int wwrapx, wwrapy, twrapx, twrapy;
|
||||
|
||||
|
@ -3514,6 +3608,17 @@ int compute_particle_interaction(int i, int k, double force[2], double *torque,
|
|||
}
|
||||
break;
|
||||
}
|
||||
case (I_COULOMB_LJ):
|
||||
{
|
||||
charge = particle[i].charge*particle[k].charge;
|
||||
f = -100.0*krepel*charge/(1.0e-12 + distance*distance);
|
||||
if (charge < 0.0)
|
||||
f += 0.01*krepel*lennard_jones_force(distance, particle[k]);
|
||||
force[0] = f*ca;
|
||||
force[1] = f*sa;
|
||||
break;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
if (ROTATION)
|
||||
|
@ -3998,7 +4103,7 @@ void compute_particle_colors(t_particle particle, int plot, double rgb[3], doubl
|
|||
if (hue > DPI) hue -= DPI;
|
||||
hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*(hue)/DPI;
|
||||
ej = particle.emean;
|
||||
if (ej < 0.1*PARTICLE_EMAX) lum = 10.0*ej/PARTICLE_EMAX;
|
||||
if (ej < 0.5*PARTICLE_EMAX) lum = 2.0*ej/PARTICLE_EMAX;
|
||||
else lum = 1.0;
|
||||
*radius = particle.radius;
|
||||
*width = BOUNDARY_WIDTH;
|
||||
|
@ -4127,6 +4232,21 @@ void set_segment_group_color(int group, double lum, double rgb[3])
|
|||
glColor3f(lum*rgb[0], lum*rgb[1], lum*rgb[2]);
|
||||
}
|
||||
|
||||
void set_segment_pressure_color(double pressure, double lum, double rgb[3])
|
||||
{
|
||||
double hue;
|
||||
|
||||
// if (pressure < 0.0) pressure = 0.0;
|
||||
hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*pressure/SEGMENT_PMAX;
|
||||
if (hue > PARTICLE_HUE_MIN) hue = PARTICLE_HUE_MIN;
|
||||
else if (hue < PARTICLE_HUE_MAX) hue = PARTICLE_HUE_MAX;
|
||||
|
||||
hsl_to_rgb_turbo(hue, 0.9, lum, rgb);
|
||||
|
||||
glColor3f(rgb[0], rgb[1], rgb[2]);
|
||||
}
|
||||
|
||||
|
||||
void draw_altitude_lines()
|
||||
{
|
||||
int i, i1, i2;
|
||||
|
@ -4625,8 +4745,18 @@ int draw_special_particle(t_particle particle, double xc1, double yc1, double ra
|
|||
void draw_one_particle(t_particle particle, double xc, double yc, double radius, double angle, int nsides, double width, double rgb[3])
|
||||
/* draw one of the particles */
|
||||
{
|
||||
double ca, sa, x1, x2, y1, y2, xc1, yc1, wangle, newradius = radius;
|
||||
double x1, x2, y1, y2, xc1, yc1, wangle, newradius = radius;
|
||||
int wsign, cont = 1, draw = 1;
|
||||
static double ca, sa;
|
||||
static int first = 1;
|
||||
|
||||
if (first)
|
||||
{
|
||||
angle = APOLY*PID;
|
||||
ca = cos(angle);
|
||||
sa = sin(angle);
|
||||
first = 0;
|
||||
}
|
||||
|
||||
if (CENTER_VIEW_ON_OBSTACLE) xc1 = xc - xshift;
|
||||
else xc1 = xc;
|
||||
|
@ -4651,14 +4781,16 @@ void draw_one_particle(t_particle particle, double xc, double yc, double radius,
|
|||
|
||||
/* draw crosses on particles of second type */
|
||||
if ((TWO_TYPES)&&(DRAW_CROSS))
|
||||
if (particle.type == 1)
|
||||
if (particle.type == 2)
|
||||
{
|
||||
if (ROTATION) angle = angle + APOLY*PID;
|
||||
else angle = APOLY*PID;
|
||||
ca = cos(angle);
|
||||
sa = sin(angle);
|
||||
glLineWidth(3);
|
||||
glColor3f(0.0, 0.0, 0.0);
|
||||
if (ROTATION)
|
||||
{
|
||||
angle = angle + APOLY*PID;
|
||||
ca = cos(angle);
|
||||
sa = sin(angle);
|
||||
}
|
||||
glLineWidth(2);
|
||||
glColor3f(1.0, 1.0, 1.0);
|
||||
x1 = xc1 - MU_B*ca;
|
||||
y1 = yc1 - MU_B*sa;
|
||||
x2 = xc1 + MU_B*ca;
|
||||
|
@ -5002,6 +5134,7 @@ void draw_container(double xmin, double xmax, t_obstacle obstacle[NMAXOBSTACLES]
|
|||
for (i = 0; i < nsegments; i++) if (segment[i].active)
|
||||
{
|
||||
if (COLOR_SEG_GROUPS) set_segment_group_color(segment[i].group, 1.0, rgb);
|
||||
else if (SHOW_SEGMENTS_PRESSURE) set_segment_pressure_color(segment[i].avrg_pressure, 1.0, rgb);
|
||||
draw_line(segment[i].x1 - xtrack, segment[i].y1 - ytrack, segment[i].x2 - xtrack, segment[i].y2 - ytrack);
|
||||
}
|
||||
}
|
||||
|
@ -5248,6 +5381,12 @@ void print_omega(double angle, double angular_speed, double fx, double fy)
|
|||
glColor3f(1.0, 1.0, 1.0);
|
||||
sprintf(message, "Fy = %.4f", fy);
|
||||
write_text(xrighttext + 0.1, y1, message);
|
||||
|
||||
y1 -= 0.1;
|
||||
erase_area_hsl(xrightbox, y1 + 0.025, 0.42, 0.05, 0.0, 0.9, 0.0);
|
||||
glColor3f(1.0, 1.0, 1.0);
|
||||
sprintf(message, "Fy/Fx = %.4f", fy/fx);
|
||||
write_text(xrighttext + 0.1, y1, message);
|
||||
}
|
||||
}
|
||||
|
||||
|
@ -5389,10 +5528,17 @@ void print_parameters(t_lj_parameters params, short int left, double pressure[N_
|
|||
}
|
||||
else if (INCREASE_KREPEL) /* print force constant */
|
||||
{
|
||||
erase_area_hsl(xbox, y + 0.025*scale, 0.22*scale, 0.05*scale, 0.0, 0.9, 0.0);
|
||||
erase_area_hsl(xbox, y + 0.025*scale, 0.35*scale, 0.05*scale, 0.0, 0.9, 0.0);
|
||||
glColor3f(1.0, 1.0, 1.0);
|
||||
sprintf(message, "Force %.0f", params.krepel);
|
||||
write_text(xtext + 0.28, y, message);
|
||||
sprintf(message, "Force constant %.0f", params.krepel);
|
||||
write_text(xtext + 0.03, y, message);
|
||||
}
|
||||
else if (INCREASE_E) /* print electric field */
|
||||
{
|
||||
erase_area_hsl(xbox, y + 0.025*scale, 0.27*scale, 0.05*scale, 0.0, 0.9, 0.0);
|
||||
glColor3f(1.0, 1.0, 1.0);
|
||||
sprintf(message, "E field = %.2f", 25.0*NVID*DT_PARTICLE*params.efield);
|
||||
write_text(xtext + 0.08, y, message);
|
||||
}
|
||||
|
||||
if (RECORD_PRESSURES)
|
||||
|
@ -5869,9 +6015,7 @@ void print_particles_speeds(t_particle particle[NMAXCIRCLES])
|
|||
write_text(xrighttext + 0.1, y, message);
|
||||
}
|
||||
|
||||
double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NMAXOBSTACLES],
|
||||
t_segment segment[NMAXSEGMENTS],
|
||||
double xleft, double xright, double *pleft, double *pright, double pressure[N_PRESSURES], int wall)
|
||||
double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NMAXOBSTACLES], t_segment segment[NMAXSEGMENTS], double xleft, double xright, double *pleft, double *pright, double pressure[N_PRESSURES], int wall)
|
||||
{
|
||||
int i, k;
|
||||
double xmin, xmax, ymin, ymax, padding, r, rp, r2, cphi, sphi, f, fperp = 0.0, x, y, xtube, distance, dx, dy, width, ybin, angle, x1, x2, h, ytop, norm, dleft, dplus, dminus, tmp_pleft = 0.0, tmp_pright = 0.0, proj, pscal, pvect, pvmin;
|
||||
|
@ -5914,8 +6058,6 @@ double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacl
|
|||
particle[j].fx += f*segment[i].nx;
|
||||
particle[j].fy += f*segment[i].ny;
|
||||
|
||||
|
||||
|
||||
if ((MOVE_BOUNDARY)||(MOVE_SEGMENT_GROUPS)||(PRINT_SEGMENTS_FORCE))
|
||||
{
|
||||
segment[i].fx -= f*segment[i].nx;
|
||||
|
@ -5923,6 +6065,14 @@ double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacl
|
|||
segment[i].torque -= (x - segment[i].xc)*f*segment[i].ny - (y - segment[i].yc)*f*segment[i].nx;
|
||||
// printf("Segment %i: f = (%.3lg, %.3lg)\n", i, segment[i].fx, segment[i].fy);
|
||||
}
|
||||
|
||||
if (SHOW_SEGMENTS_PRESSURE)
|
||||
{
|
||||
segment[i].pressure = f/segment[i].length;
|
||||
segment[i].avrg_pressure *= (1.0 - P_AVRG_FACTOR);
|
||||
segment[i].avrg_pressure += P_AVRG_FACTOR*segment[i].pressure;
|
||||
// printf("Segment %i: pressure = %.3lg\n", i, segment[i].avrg_pressure);
|
||||
}
|
||||
}
|
||||
if ((VICSEK_INT)&&(vabs(distance) < 1.5*r))
|
||||
{
|
||||
|
@ -6480,21 +6630,51 @@ int reorder_particles(t_particle particle[NMAXCIRCLES], double py[NMAXCIRCLES],
|
|||
}
|
||||
|
||||
|
||||
int twotype_config(int i, t_particle particle[NMAXCIRCLES])
|
||||
/* assign different particle types */
|
||||
{
|
||||
switch (TWOTYPE_CONFIG) {
|
||||
case (TTC_RANDOM): return((double)rand()/RAND_MAX > TYPE_PROPORTION);
|
||||
case (TTC_CHESSBOARD):
|
||||
{
|
||||
switch (CIRCLE_PATTERN) {
|
||||
case (C_SQUARE): return((i%NGRIDY + i/NGRIDY)%2 == 0);
|
||||
case (C_HEX): return(i%2 == 0);
|
||||
default: return(i%2 == 0);
|
||||
}
|
||||
}
|
||||
case (TTC_COANDA):
|
||||
{
|
||||
if (vabs(particle[i].yc) < LAMBDA) return(1);
|
||||
if (vabs(particle[i].yc) < LAMBDA + MU) return(-1);
|
||||
return(0);
|
||||
}
|
||||
default: return(0);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
int initialize_configuration(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HASHX*HASHY], t_obstacle obstacle[NMAXOBSTACLES], double px[NMAXCIRCLES], double py[NMAXCIRCLES], double pangle[NMAXCIRCLES], int tracer_n[N_TRACER_PARTICLES],
|
||||
t_segment segment[NMAXSEGMENTS])
|
||||
/* initialize all particles, obstacles, and the hashgrid */
|
||||
{
|
||||
int i, j, k, n, nactive = 0;
|
||||
int i, j, k, n, type, nactive = 0;
|
||||
double x, y, h, xx, yy, rnd;
|
||||
|
||||
for (i=0; i < ncircles; i++)
|
||||
{
|
||||
/* set particle type */
|
||||
particle[i].type = 0;
|
||||
if ((TWO_TYPES)&&((double)rand()/RAND_MAX > TYPE_PROPORTION))
|
||||
// if ((TWO_TYPES)&&((double)rand()/RAND_MAX > TYPE_PROPORTION))
|
||||
if ((TWO_TYPES)&&(twotype_config(i, particle)))
|
||||
{
|
||||
particle[i].type = 2;
|
||||
particle[i].radius = MU_B;
|
||||
type = twotype_config(i, particle);
|
||||
if (type == 1)
|
||||
{
|
||||
particle[i].type = 2;
|
||||
particle[i].radius = MU_B;
|
||||
}
|
||||
else if (type == -1) particle[i].active = 0;
|
||||
}
|
||||
if ((INTERACTION == I_VICSEK_SHARK)&&(i==1))
|
||||
{
|
||||
|
@ -6524,6 +6704,7 @@ t_segment segment[NMAXSEGMENTS])
|
|||
particle[i].spin_freq = SPIN_INTER_FREQUENCY;
|
||||
particle[i].mass_inv = 1.0/PARTICLE_MASS;
|
||||
particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT;
|
||||
particle[i].charge = CHARGE;
|
||||
}
|
||||
else
|
||||
{
|
||||
|
@ -6533,10 +6714,25 @@ t_segment segment[NMAXSEGMENTS])
|
|||
particle[i].spin_freq = SPIN_INTER_FREQUENCY_B;
|
||||
particle[i].mass_inv = 1.0/PARTICLE_MASS_B;
|
||||
particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT_B;
|
||||
particle[i].charge = CHARGE_B;
|
||||
}
|
||||
|
||||
switch (V_INITIAL_TYPE)
|
||||
{
|
||||
case (VI_RANDOM):
|
||||
{
|
||||
particle[i].vx = V_INITIAL*gaussian();
|
||||
particle[i].vy = V_INITIAL*gaussian();
|
||||
break;
|
||||
}
|
||||
case (VI_COANDA):
|
||||
{
|
||||
if (vabs(particle[i].yc) < LAMBDA) particle[i].vx = V_INITIAL;
|
||||
else particle[i].vx = 0.0;
|
||||
particle[i].vy = 0.0;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
particle[i].vx = V_INITIAL*gaussian();
|
||||
particle[i].vy = V_INITIAL*gaussian();
|
||||
|
||||
if ((INTERACTION == I_VICSEK_SHARK)&&(i==1))
|
||||
{
|
||||
|
@ -6596,6 +6792,7 @@ t_segment segment[NMAXSEGMENTS])
|
|||
particle[i].dirmean = 0.0;
|
||||
particle[i].mass_inv = 1.0/PARTICLE_MASS;
|
||||
particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT;
|
||||
particle[i].charge = 0.0;
|
||||
particle[i].vx = 0.0;
|
||||
particle[i].vy = 0.0;
|
||||
px[i] = 0.0;
|
||||
|
|
|
@ -110,8 +110,8 @@ int writetiff(char *filename, char *description, int x, int y, int width, int he
|
|||
glPixelStorei(GL_PACK_ALIGNMENT, 1);
|
||||
|
||||
glReadPixels(x, y, width, height, GL_RGB, GL_UNSIGNED_BYTE, image);
|
||||
TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32) width);
|
||||
TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32) height);
|
||||
TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32_t) width);
|
||||
TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32_t) height);
|
||||
TIFFSetField(file, TIFFTAG_BITSPERSAMPLE, 8);
|
||||
TIFFSetField(file, TIFFTAG_COMPRESSION, compression);
|
||||
TIFFSetField(file, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB);
|
||||
|
|
|
@ -112,8 +112,8 @@ int writetiff(char *filename, char *description, int x, int y, int width, int he
|
|||
glPixelStorei(GL_PACK_ALIGNMENT, 1);
|
||||
|
||||
glReadPixels(x, y, width, height, GL_RGB, GL_UNSIGNED_BYTE, image);
|
||||
TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32) width);
|
||||
TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32) height);
|
||||
TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32_t) width);
|
||||
TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32_t) height);
|
||||
TIFFSetField(file, TIFFTAG_BITSPERSAMPLE, 8);
|
||||
TIFFSetField(file, TIFFTAG_COMPRESSION, compression);
|
||||
TIFFSetField(file, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB);
|
||||
|
|
|
@ -50,8 +50,8 @@ int writetiff(char *filename, char *description, int x, int y, int width, int he
|
|||
glPixelStorei(GL_PACK_ALIGNMENT, 1);
|
||||
|
||||
glReadPixels(x, y, width, height, GL_RGB, GL_UNSIGNED_BYTE, image);
|
||||
TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32) width);
|
||||
TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32) height);
|
||||
TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32_t) width);
|
||||
TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32_t) height);
|
||||
TIFFSetField(file, TIFFTAG_BITSPERSAMPLE, 8);
|
||||
TIFFSetField(file, TIFFTAG_COMPRESSION, compression);
|
||||
TIFFSetField(file, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB);
|
||||
|
|
835
sub_sphere.c
835
sub_sphere.c
File diff suppressed because it is too large
Load Diff
68
sub_wave.c
68
sub_wave.c
|
@ -30,8 +30,8 @@ int writetiff_new(char *filename, char *description, int x, int y, int width, in
|
|||
glPixelStorei(GL_PACK_ALIGNMENT, 1);
|
||||
|
||||
glReadPixels(x, y, width, height, GL_RGB, GL_UNSIGNED_BYTE, image);
|
||||
TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32) width);
|
||||
TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32) height);
|
||||
TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32_t) width);
|
||||
TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32_t) height);
|
||||
TIFFSetField(file, TIFFTAG_BITSPERSAMPLE, 8);
|
||||
TIFFSetField(file, TIFFTAG_COMPRESSION, compression);
|
||||
TIFFSetField(file, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB);
|
||||
|
@ -85,8 +85,8 @@ int writetiff(char *filename, char *description, int x, int y, int width, int he
|
|||
glPixelStorei(GL_PACK_ALIGNMENT, 1);
|
||||
|
||||
glReadPixels(x, y, width, height, GL_RGB, GL_UNSIGNED_BYTE, image);
|
||||
TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32) width);
|
||||
TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32) height);
|
||||
TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32_t) width);
|
||||
TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32_t) height);
|
||||
TIFFSetField(file, TIFFTAG_BITSPERSAMPLE, 8);
|
||||
TIFFSetField(file, TIFFTAG_COMPRESSION, compression);
|
||||
TIFFSetField(file, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB);
|
||||
|
@ -3130,6 +3130,13 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_
|
|||
}
|
||||
return(1);
|
||||
}
|
||||
case (D_LENS):
|
||||
{
|
||||
if (vabs(y) > MU) return(1);
|
||||
a = LAMBDA - 0.25;
|
||||
if (x > 0) return ((x+a)*(x+a) + y*y > LAMBDA*LAMBDA);
|
||||
return ((x-a)*(x-a) + y*y > LAMBDA*LAMBDA);
|
||||
}
|
||||
case (D_MENGER):
|
||||
{
|
||||
x1 = 0.5*(x+1.0);
|
||||
|
@ -4537,6 +4544,22 @@ void draw_billiard(int fade, double fade_value) /* draws the billiard bound
|
|||
}
|
||||
break;
|
||||
}
|
||||
case (D_LENS):
|
||||
{
|
||||
a = LAMBDA - 0.25;
|
||||
if (first)
|
||||
{
|
||||
arcangle = asin(MU/LAMBDA);
|
||||
ll = LAMBDA*cos(arcangle) - a;
|
||||
first = 0;
|
||||
}
|
||||
|
||||
draw_circle_arc(-a, 0.0, LAMBDA, -arcangle, 2.0*arcangle, NSEG);
|
||||
draw_line(ll, MU, -ll, MU);
|
||||
draw_circle_arc(a, 0.0, LAMBDA, PI-arcangle, 2.0*arcangle, NSEG);
|
||||
draw_line(-ll, -MU, ll, -MU);
|
||||
break;
|
||||
}
|
||||
case (D_MENGER):
|
||||
{
|
||||
glLineWidth(3);
|
||||
|
@ -6330,7 +6353,7 @@ void init_wave_packets(t_wave_packet *packet, int radius)
|
|||
|
||||
printf("Initialising wave packets\n");
|
||||
switch (WAVE_PACKET_SOURCE_TYPE) {
|
||||
case (0):
|
||||
case (WP_RANDOM1):
|
||||
{
|
||||
nx = (int)sqrt((double)N_WAVE_PACKETS);
|
||||
ny = N_WAVE_PACKETS/nx;
|
||||
|
@ -6356,15 +6379,12 @@ void init_wave_packets(t_wave_packet *packet, int radius)
|
|||
|
||||
break;
|
||||
}
|
||||
case (1):
|
||||
case (WP_RANDOM2):
|
||||
{
|
||||
for (i=0; i<N_WAVE_PACKETS; i++)
|
||||
{
|
||||
// j = i/nx;
|
||||
// k = i%nx;
|
||||
packet[i].xc = XMIN + 0.15*(XMAX - XMIN)*(double)rand()/RAND_MAX;
|
||||
packet[i].yc = 0.4*(YMAX - YMIN)*((double)rand()/RAND_MAX - 0.5);
|
||||
// packet[i].period = 30.0*(1.0 + 0.5*(double)rand()/RAND_MAX);
|
||||
packet[i].period = 50.0*(1.0 + 0.5*(double)rand()/RAND_MAX);
|
||||
packet[i].amp = INITIAL_AMP;
|
||||
packet[i].phase = DPI*(double)rand()/RAND_MAX;
|
||||
|
@ -6377,6 +6397,27 @@ void init_wave_packets(t_wave_packet *packet, int radius)
|
|||
packet[i].iy = ij[1];
|
||||
}
|
||||
|
||||
break;
|
||||
}
|
||||
case (WP_PAIR):
|
||||
{
|
||||
for (i=0; i<2; i++)
|
||||
{
|
||||
packet[i].xc = -1.25;
|
||||
packet[i].yc = 0.06;
|
||||
if (i==1) packet[i].yc *= -1.0;
|
||||
packet[i].period = OSCILLATING_SOURCE_PERIOD;
|
||||
packet[i].amp = INITIAL_AMP;
|
||||
packet[i].phase = 0.0;
|
||||
packet[i].var_envelope = 100000.0;
|
||||
packet[i].time_shift = 0;
|
||||
|
||||
xy_to_ij(packet[i].xc, packet[i].yc, ij);
|
||||
if(ij[0] <= radius) ij[0] = radius+1;
|
||||
packet[i].ix = ij[0];
|
||||
packet[i].iy = ij[1];
|
||||
}
|
||||
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
@ -6390,18 +6431,23 @@ double wave_packet_height(double t, t_wave_packet packet, int type_envelope)
|
|||
cwave = packet.amp*cos(DPI*t/packet.period + packet.phase);
|
||||
|
||||
switch (type_envelope) {
|
||||
case (0):
|
||||
case (WE_SINE):
|
||||
{
|
||||
envelope = 0.1 + sin(DPI*t/packet.var_envelope);
|
||||
envelope = envelope*envelope;
|
||||
break;
|
||||
}
|
||||
case (1):
|
||||
case (WE_CUTOFF):
|
||||
{
|
||||
if (t < (double)packet.var_envelope) envelope = 1.0;
|
||||
else envelope = 0.0;
|
||||
break;
|
||||
}
|
||||
case (WE_CONSTANT):
|
||||
{
|
||||
envelope = 1.0;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
return(cwave*envelope);
|
||||
|
|
|
@ -219,16 +219,11 @@
|
|||
#define IOR 7 /* choice of index of refraction, see list in global_pdes.c */
|
||||
#define IOR_TOTAL_TURNS 1.5 /* total angle of rotation for IOR_PERIODIC_WELLS_ROTATING */
|
||||
#define MANDEL_IOR_SCALE -0.05 /* parameter controlling dependence of IoR on Mandelbrot escape speed */
|
||||
#define COURANT 0.04 /* Courant number */
|
||||
#define COURANTB 0.0 /* Courant number in medium B */
|
||||
#define INITIAL_AMP 0.5 /* amplitude of initial condition */
|
||||
#define INITIAL_VARIANCE 0.0003 /* variance of initial condition */
|
||||
#define INITIAL_WAVELENGTH 0.015 /* wavelength of initial condition */
|
||||
#define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */
|
||||
#define WAVE_PACKET_SOURCE_TYPE 1 /* type of wave packet sources */
|
||||
#define N_WAVE_PACKETS 15 /* number of wave packets */
|
||||
#define OSCIL_LEFT_YSHIFT 0.0 /* y-dependence of left oscillation (for non-horizontal waves) */
|
||||
#define DRAW_WAVE_PROFILE 0 /* set to 1 to draw a profile of the wave */
|
||||
#define OSCILLATING_SOURCE_PERIOD 20 /* period of oscillating source */
|
||||
/* end of constants only used by sub_wave and sub_maze */
|
||||
|
||||
#include "global_pdes.c" /* constants and global variables */
|
||||
|
|
142
wave_sphere.c
142
wave_sphere.c
|
@ -53,6 +53,7 @@
|
|||
|
||||
#define DPOLE 20 /* safety distance to poles */
|
||||
#define SMOOTHPOLE 0.1 /* smoothing coefficient at poles */
|
||||
#define ZERO_MERIDIAN 0.0 /* choice of zero meridian (will be at left/right boundary of 2d plot) */
|
||||
|
||||
#define XMIN -2.0
|
||||
#define XMAX 2.0 /* x interval */
|
||||
|
@ -61,10 +62,10 @@
|
|||
|
||||
#define HIGHRES 0 /* set to 1 if resolution of grid is double that of displayed image */
|
||||
|
||||
#define JULIA_SCALE 0.25 /* scaling for Julia sets */
|
||||
#define JULIA_ROT 90.0 /* rotation of Julia set, in degrees */
|
||||
#define JULIA_RE -0.77145
|
||||
#define JULIA_IM -0.10295 /* parameters for Julia sets */
|
||||
#define JULIA_SCALE 0.6 /* scaling for Julia sets */
|
||||
#define JULIA_ROT -20.0 /* rotation of Julia set, in degrees */
|
||||
#define JULIA_RE 0.5
|
||||
#define JULIA_IM 0.462 /* parameters for Julia sets */
|
||||
|
||||
/* Choice of the billiard table */
|
||||
|
||||
|
@ -77,7 +78,7 @@
|
|||
#define CIRCLE_PATTERN_B 0 /* second pattern of circles or polygons */
|
||||
|
||||
#define VARIABLE_IOR 0 /* set to 1 for a variable index of refraction */
|
||||
#define IOR 9 /* choice of index of refraction, see list in global_pdes.c */
|
||||
#define IOR 10 /* choice of index of refraction, see list in global_pdes.c */
|
||||
#define IOR_TOTAL_TURNS 1.0 /* total angle of rotation for IOR_PERIODIC_WELLS_ROTATING */
|
||||
#define MANDEL_IOR_SCALE -0.05 /* parameter controlling dependence of IoR on Mandelbrot escape speed */
|
||||
|
||||
|
@ -114,7 +115,7 @@
|
|||
|
||||
/* Physical parameters of wave equation */
|
||||
|
||||
#define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */
|
||||
#define TWOSPEEDS 1 /* set to 1 to replace hardcore boundary by medium with different speed */
|
||||
#define OSCILLATE_LEFT 0 /* set to 1 to add oscilating boundary condition on the left */
|
||||
#define OSCILLATE_TOPBOT 0 /* set to 1 to enforce a planar wave on top and bottom boundary */
|
||||
#define OSCILLATION_SCHEDULE 3 /* oscillation schedule, see list in global_pdes.c */
|
||||
|
@ -124,7 +125,7 @@
|
|||
#define ACHIRP 0.2 /* acceleration coefficient in chirp */
|
||||
#define DAMPING 0.0 /* damping of periodic excitation */
|
||||
#define COURANT 0.05 /* Courant number */
|
||||
#define COURANTB 0.01 /* Courant number in medium B */
|
||||
#define COURANTB 0.002 /* Courant number in medium B */
|
||||
#define GAMMA 0.0 /* damping factor in wave equation */
|
||||
#define GAMMAB 1.0e-6 /* damping factor in wave equation */
|
||||
#define GAMMA_SIDES 1.0e-4 /* damping factor on boundary */
|
||||
|
@ -155,7 +156,7 @@
|
|||
|
||||
/* Parameters for length and speed of simulation */
|
||||
|
||||
#define NSTEPS 3600 /* number of frames of movie */
|
||||
#define NSTEPS 2500 /* number of frames of movie */
|
||||
#define NVID 4 /* number of iterations between images displayed on screen */
|
||||
#define NSEG 1000 /* number of segments of boundary */
|
||||
#define INITIAL_TIME 0 /* time after which to start saving frames */
|
||||
|
@ -167,15 +168,14 @@
|
|||
#define SLEEP1 1 /* initial sleeping time */
|
||||
#define SLEEP2 1 /* final sleeping time */
|
||||
#define MID_FRAMES 100 /* number of still frames between parts of two-part movie */
|
||||
#define END_FRAMES 100 /* number of still frames at end of movie */
|
||||
#define END_FRAMES 500 /* number of still frames at end of movie */
|
||||
#define FADE 1 /* set to 1 to fade at end of movie */
|
||||
#define ROTATE_VIEW_WHILE_FADE 1 /* set to 1 to keep rotating viewpoint during fade */
|
||||
|
||||
/* Parameters of initial condition */
|
||||
|
||||
#define INITIAL_AMP 0.75 /* amplitude of initial condition */
|
||||
#define INITIAL_VARIANCE 0.0005 /* variance of initial condition */
|
||||
#define INITIAL_WAVELENGTH 0.025 /* wavelength of initial condition */
|
||||
#define INITIAL_AMP 0.5 /* amplitude of initial condition */
|
||||
#define INITIAL_VARIANCE 0.001 /* variance of initial condition */
|
||||
#define INITIAL_WAVELENGTH 0.05 /* wavelength of initial condition */
|
||||
|
||||
/* Plot type, see list in global_pdes.c */
|
||||
|
||||
|
@ -188,15 +188,15 @@
|
|||
#define CHANGE_LUMINOSITY 1 /* set to 1 to let luminosity depend on energy flux intensity */
|
||||
#define FLUX_WINDOW 30 /* size of averaging window of flux intensity */
|
||||
#define AMPLITUDE_HIGH_RES 1 /* set to 1 to increase resolution of plot */
|
||||
#define SHADE_3D 1 /* set to 1 to change luminosity according to normal vector */
|
||||
#define SHADE_2D 0 /* set to 1 to change luminosity according to normal vector to plane */
|
||||
#define SHADE_3D 0 /* set to 1 to change luminosity according to normal vector */
|
||||
#define SHADE_2D 1 /* set to 1 to change luminosity according to normal vector to plane */
|
||||
#define SHADE_WAVE 1 /* set to 1 to have luminosity depend on wave height */
|
||||
#define NON_DIRICHLET_BC 0 /* set to 1 to draw only facets in domain, if field is not zero on boundary */
|
||||
#define FLOOR_ZCOORD 1 /* set to 1 to draw only facets with z not too negative */
|
||||
#define DRAW_BILLIARD 0 /* set to 1 to draw boundary */
|
||||
#define DRAW_BILLIARD_FRONT 0 /* set to 1 to draw front of boundary after drawing wave */
|
||||
#define DRAW_CONSTRUCTION_LINES 0 /* set to 1 to draw construction lines of certain domains */
|
||||
#define FADE_IN_OBSTACLE 1 /* set to 1 to fade color inside obstacles */
|
||||
#define FADE_IN_OBSTACLE 0 /* set to 1 to fade color inside obstacles */
|
||||
#define DRAW_OUTSIDE_GRAY 0 /* experimental, draw outside of billiard in gray */
|
||||
#define SHADE_SCALE_2D 10.0 /* controls "depth" of 2D shading */
|
||||
#define COS_LIGHT_MIN 0.0 /* controls angle-dependence of 2D shading */
|
||||
|
@ -208,16 +208,19 @@
|
|||
/* 3D representation */
|
||||
|
||||
#define REPRESENTATION_3D 1 /* choice of 3D representation */
|
||||
#define PLOT_2D 0 /* switch to 2D representation, equirectangular projection */
|
||||
#define PLOT_2D 1 /* switch to 2D representation, equirectangular projection */
|
||||
#define PHISHIFT 0.0 /* shift of phi in 2D plot (in degrees) */
|
||||
#define FLOODING 1 /* set to 1 to draw waves above altitude (for Earth representations) */
|
||||
|
||||
#define REP_AXO_3D 0 /* linear projection (axonometry) */
|
||||
#define REP_PROJ_3D 1 /* projection on plane orthogonal to observer line of sight */
|
||||
|
||||
#define ROTATE_VIEW 1 /* set to 1 to rotate position of observer */
|
||||
#define ROTATE_ANGLE -360.0 /* total angle of rotation during simulation */
|
||||
#define ROTATE_ANGLE 75.0 /* total angle of rotation during simulation */
|
||||
#define ROTATE_VIEW_WHILE_FADE 1 /* set to 1 to keep rotating viewpoint during fade */
|
||||
|
||||
#define VIEWPOINT_TRAJ 1 /* type of viewpoint trajectory */
|
||||
#define VIEWPOINT_TRAJ 2 /* type of viewpoint trajectory */
|
||||
#define MAX_LATITUDE 45.0 /* maximal latitude for viewpoint trajectory VP_ORBIT2 */
|
||||
|
||||
/* Color schemes */
|
||||
|
||||
|
@ -233,18 +236,18 @@
|
|||
|
||||
#define SCALE 0 /* set to 1 to adjust color scheme to variance of field */
|
||||
#define SLOPE 1.0 /* sensitivity of color on wave amplitude */
|
||||
#define VSCALE_AMPLITUDE 1.0 /* additional scaling factor for color scheme P_3D_AMPLITUDE */
|
||||
#define VSCALE_ENERGY 10.0 /* additional scaling factor for color scheme P_3D_ENERGY */
|
||||
#define VSCALE_AMPLITUDE 1.5 /* additional scaling factor for color scheme P_3D_AMPLITUDE */
|
||||
#define VSCALE_ENERGY 4.0 /* additional scaling factor for color scheme P_3D_ENERGY */
|
||||
#define PHASE_FACTOR 20.0 /* factor in computation of phase in color scheme P_3D_PHASE */
|
||||
#define PHASE_SHIFT 0.0 /* shift of phase in color scheme P_3D_PHASE */
|
||||
#define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */
|
||||
#define E_SCALE 100.0 /* scaling factor for energy representation */
|
||||
#define E_SCALE 150.0 /* scaling factor for energy representation */
|
||||
#define LOG_SCALE 0.75 /* scaling factor for energy log representation */
|
||||
#define LOG_SHIFT 0.5 /* shift of colors on log scale */
|
||||
#define LOG_ENERGY_FLOOR -10.0 /* floor value for log of (total) energy */
|
||||
#define LOG_MEAN_ENERGY_SHIFT 1.0 /* additional shift for log of mean energy */
|
||||
#define FLUX_SCALE 600.0 /* scaling factor for energy flux representation */
|
||||
#define FLUX_CSCALE 5.0 /* scaling factor for color in energy flux representation */
|
||||
#define FLUX_SCALE 1200.0 /* scaling factor for energy flux representation */
|
||||
#define FLUX_CSCALE 2.0 /* scaling factor for color in energy flux representation */
|
||||
#define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */
|
||||
|
||||
#define COLORHUE 260 /* initial hue of water color for scheme C_LUM */
|
||||
|
@ -261,9 +264,9 @@
|
|||
#define MAZE_XSHIFT 0.0 /* horizontal shift of maze */
|
||||
#define MAZE_WIDTH 0.02 /* half width of maze walls */
|
||||
|
||||
#define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */
|
||||
#define COLORBAR_RANGE 3.0 /* scale of color scheme bar */
|
||||
#define COLORBAR_RANGE_B 2.0 /* 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 3.0 /* scale of color scheme bar for 2nd part */
|
||||
#define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */
|
||||
#define CIRC_COLORBAR 0 /* set to 1 to draw circular color scheme */
|
||||
#define CIRC_COLORBAR_B 0 /* set to 1 to draw circular color scheme */
|
||||
|
@ -276,7 +279,6 @@
|
|||
#define POT_FACT 20.0
|
||||
/* end of constants only used by sub_wave and sub_maze */
|
||||
|
||||
|
||||
/* For debugging purposes only */
|
||||
#define FLOOR 1 /* set to 1 to limit wave amplitude to VMAX */
|
||||
#define VMAX 10.0 /* max value of wave amplitude */
|
||||
|
@ -286,18 +288,27 @@
|
|||
double u_3d[2] = {0.75, -0.45}; /* projections of basis vectors for REP_AXO_3D representation */
|
||||
double v_3d[2] = {-0.75, -0.45};
|
||||
double w_3d[2] = {0.0, 0.015};
|
||||
double light[3] = {0.816496581, -0.40824829, 0.40824829}; /* vector of "light" direction for P_3D_ANGLE color scheme */
|
||||
double observer[3] = {6.0, 8.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] = {-5.0, 8.0, -7.0}; /* location of observer for REP_PROJ_3D representation */
|
||||
int reset_view = 0; /* switch to reset 3D view parameters (for option ROTATE_VIEW) */
|
||||
|
||||
#define ADD_DEM 1 /* add DEM (digital elevation model) */
|
||||
#define ADD_NEGATIVE_DEM 0 /* add DEM with bathymetric data */
|
||||
#define RSCALE_DEM 0.1 /* scaling factor of radial component for DEM */
|
||||
#define SMOOTH_DEM 0 /* set to 1 to smoothen DEM (to make altitude less constant) */
|
||||
#define DEM_SMOOTH_STEPS 10 /* number of smoothening steps */
|
||||
#define DEM_SMOOTH_HEIGHT 0.5 /* relative height below which to smoothen */
|
||||
#define DEM_MAXHEIGHT 8000 /* max height of DEM (estimated from Everest) */
|
||||
#define PLANET_SEALEVEL 2500.0 /* sea level for flooded planet */
|
||||
|
||||
#define RSCALE 0.01 /* scaling factor of radial component */
|
||||
#define RMAX 10.0 /* max value of radial component */
|
||||
#define Z_SCALING_FACTOR 0.8 /* overall scaling factor of z axis for REP_PROJ_3D representation */
|
||||
#define Z_SCALING_FACTOR 0.85 /* overall scaling factor of z axis for REP_PROJ_3D representation */
|
||||
#define XY_SCALING_FACTOR 2.0 /* overall scaling factor for on-screen (x,y) coordinates after projection */
|
||||
#define ZMAX_FACTOR 1.0 /* max value of z coordinate for REP_PROJ_3D representation */
|
||||
#define XSHIFT_3D 0.0 /* overall x shift for REP_PROJ_3D representation */
|
||||
#define YSHIFT_3D 0.0 /* overall y shift for REP_PROJ_3D representation */
|
||||
#define COS_VISIBLE -0.75 /* limit on cosine of normal to shown facets */
|
||||
#define COS_VISIBLE -0.5 /* limit on cosine of normal to shown facets */
|
||||
|
||||
#include "global_pdes.c" /* constants and global variables */
|
||||
#include "global_3d.c" /* constants and global variables */
|
||||
|
@ -498,25 +509,36 @@ void evolve_wave(double phi[NX*NY], double psi[NX*NY], double tmp[NX*NY], short
|
|||
void draw_color_bar_palette(int plot, double range, int palette, int circular, int fade, double fade_value)
|
||||
{
|
||||
// double width = 0.2;
|
||||
double width = 0.14;
|
||||
double width = 0.12;
|
||||
// double width = 0.2;
|
||||
|
||||
width *= (double)NX/(double)WINWIDTH;
|
||||
|
||||
if (ROTATE_COLOR_SCHEME)
|
||||
draw_color_scheme_palette_3d(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range, palette, fade, fade_value);
|
||||
else if (circular)
|
||||
draw_circular_color_scheme_palette_3d(XMAX - 2.0*width, YMIN + 2.0*width, 1.5*width, 1.3*width, plot, -range, range, palette, fade, fade_value);
|
||||
else
|
||||
draw_color_scheme_palette_3d(XMAX - 1.5*width, YMIN + 0.1, XMAX - 0.5*width, YMAX - 0.1, plot, -range, range, palette, fade, fade_value);
|
||||
if (PLOT_2D)
|
||||
{
|
||||
if (ROTATE_COLOR_SCHEME)
|
||||
draw_color_scheme_palette_fade(XMIN + 0.6, YMIN + 0.1, XMAX - 0.6, YMIN + 0.1 + width, plot, -range, range, palette, fade, fade_value);
|
||||
else
|
||||
draw_color_scheme_palette_fade(XMAX - 1.5*width, YMIN + 0.1, XMAX - 0.5*width, YMAX - 0.1, plot, -range, range, palette, fade, fade_value);
|
||||
}
|
||||
else
|
||||
{
|
||||
if (ROTATE_COLOR_SCHEME)
|
||||
draw_color_scheme_palette_3d(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range, palette, fade, fade_value);
|
||||
else if (circular)
|
||||
draw_circular_color_scheme_palette_3d(XMAX - 2.0*width, YMIN + 2.0*width, 1.5*width, 1.3*width, plot, -range, range, palette, fade, fade_value);
|
||||
else
|
||||
draw_color_scheme_palette_3d(XMAX - 1.5*width, YMIN + 0.1, XMAX - 0.5*width, YMAX - 0.1, plot, -range, range, palette, fade, fade_value);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
void viewpoint_schedule(int i)
|
||||
/* change position of observer */
|
||||
{
|
||||
int j;
|
||||
double angle, ca, sa, r1;
|
||||
static double observer_initial[3], r, ratio;
|
||||
double angle, ca, sa, r1, interpolate, rho;
|
||||
static double observer_initial[3], r, ratio, rho0, zmax;
|
||||
static int first = 1;
|
||||
|
||||
if (first)
|
||||
|
@ -525,10 +547,14 @@ void viewpoint_schedule(int i)
|
|||
r1 = observer[0]*observer[0] + observer[1]*observer[1];
|
||||
r = sqrt(r1 + observer[2]*observer[2]);
|
||||
ratio = r/sqrt(r1);
|
||||
rho0 = module2(observer[0], observer[1]);
|
||||
if (vabs(rho0) < 0.001) rho0 = 0.001;
|
||||
zmax = r*sin(MAX_LATITUDE*PI/180.0);
|
||||
first = 0;
|
||||
}
|
||||
|
||||
angle = (ROTATE_ANGLE*DPI/360.0)*(double)i/(double)NSTEPS;
|
||||
interpolate = (double)i/(double)NSTEPS;
|
||||
angle = (ROTATE_ANGLE*DPI/360.0)*interpolate;
|
||||
ca = cos(angle);
|
||||
sa = sin(angle);
|
||||
switch (VIEWPOINT_TRAJ)
|
||||
|
@ -546,6 +572,21 @@ void viewpoint_schedule(int i)
|
|||
observer[2] = ca*observer_initial[2];
|
||||
break;
|
||||
}
|
||||
case (VP_ORBIT2):
|
||||
{
|
||||
observer[0] = ca*observer_initial[0] - sa*observer_initial[1]*ratio;
|
||||
observer[1] = ca*observer_initial[1] + sa*observer_initial[0]*ratio;
|
||||
observer[2] = sa*zmax;
|
||||
break;
|
||||
}
|
||||
case (VP_POLAR):
|
||||
{
|
||||
rho = -sa*observer_initial[2] + ca*rho0;
|
||||
observer[0] = observer_initial[0]*rho/rho0;
|
||||
observer[1] = observer_initial[1]*rho/rho0;
|
||||
observer[2] = ca*observer_initial[2] + sa*rho0;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
printf("Angle %.3lg, Observer position (%.3lg, %.3lg, %.3lg)\n", angle, observer[0], observer[1], observer[2]);
|
||||
|
@ -630,10 +671,7 @@ void animation()
|
|||
courant2 = COURANT*COURANT;
|
||||
courantb2 = COURANTB*COURANTB;
|
||||
c = COURANT*(XMAX - XMIN)/(double)NX;
|
||||
// a = 0.015;
|
||||
// b = 0.0003;
|
||||
// a = 0.04;
|
||||
// b = 0.0018;
|
||||
|
||||
a = 0.05;
|
||||
b = 0.0016;
|
||||
|
||||
|
@ -674,16 +712,18 @@ void animation()
|
|||
// for (j=1; j<NPOLY; j++)
|
||||
// add_circular_wave_mod(1.0, lambda1*cos(((double)j+0.5)*angle), lambda1*sin(((double)j+0.5)*angle), phi, psi, xy_in);
|
||||
|
||||
init_circular_wave_sphere(0.7, 0.5, phi, psi, xy_in, wsphere);
|
||||
// init_wave_flat_sphere(phi, psi, xy_in, wsphere);
|
||||
// init_circular_wave_sphere(0.25 + PID, 0.0, phi, psi, xy_in, wsphere);
|
||||
// add_circular_wave_sphere(-1.0, 0.25 + 3.0*PID, 0.0, phi, psi, xy_in, wsphere);
|
||||
|
||||
init_circular_wave_sphere(123.3*PI/180.0, -48.8*PI/180.0, phi, psi, xy_in, wsphere);
|
||||
|
||||
// add_circular_wave_sphere(1.0, 1.0 - PI/9.0 + DPI/3.0, 0.15, phi, psi, xy_in, wsphere);
|
||||
// add_circular_wave_sphere(1.0, 1.0 - PI/9.0 + 2.0*DPI/3.0, 0.15, phi, psi, xy_in, wsphere);
|
||||
|
||||
// printf("Wave initialized\n");
|
||||
|
||||
/* initialize table of wave speeds/dissipation */
|
||||
init_speed_dissipation(xy_in, tc, tcc, tgamma);
|
||||
init_speed_dissipation_sphere(xy_in, tc, tcc, tgamma, wsphere);
|
||||
|
||||
|
||||
|
||||
/* initialze potential to add to z coordinate */
|
||||
if (ADD_POTENTIAL)
|
||||
|
|
Loading…
Reference in New Issue