diff --git a/Earth_Map_Blue_Marble_2002_large.ppm.gz b/Earth_Map_Blue_Marble_2002_large.ppm.gz new file mode 100644 index 0000000..ea9834e Binary files /dev/null and b/Earth_Map_Blue_Marble_2002_large.ppm.gz differ diff --git a/global_3d.c b/global_3d.c index 14a3de8..16ca4f8 100644 --- a/global_3d.c +++ b/global_3d.c @@ -61,6 +61,11 @@ #define SH_COAST 3 /* depth varying with x-coordinate */ #define SH_COAST_MONOTONE 4 /* depth decreasing with x-coordinate */ +/* Type of rotating viewpoint */ + +#define VP_HORIZONTAL 0 /* rotate in a horizontal plane */ +#define VP_ORBIT 1 /* rotate in a plane containing the origin */ + /* macros to avoid unnecessary computations in 3D plots */ #define COMPUTE_THETA ((cplot == Z_POLAR)||(cplot == Z_NORM_GRADIENT)||(cplot == Z_ANGLE_GRADIENT)||(cplot == Z_NORM_GRADIENT_INTENSITY)||(cplot == Z_VORTICITY)||(cplot == Z_VORTICITY_ABS)) @@ -80,6 +85,7 @@ #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 NMAXCIRC_SPHERE 100 /* max number of circles on sphere */ int global_time = 0; double max_depth = 1.0; @@ -137,3 +143,26 @@ typedef struct double gradx, grady; /* gradient of water depth */ } t_swater_depth; + +typedef struct +{ + double phi, theta; /* phi, theta angles */ + double cphi, sphi; /* cos and sin of phi */ + double ctheta, stheta, cottheta; /* cos, sin and cotangent of theta */ + double x, y, z; /* x, y, z coordinates of point on sphere */ + double radius; /* radius with wave height */ + 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 */ +} t_wave_sphere; + + +typedef struct +{ + double phi, theta; /* longitude, latitude */ + double radius; /* radius */ + double x, y, z; /* x, y, z coordinates of point on sphere */ +} t_circles_sphere; + + +t_circles_sphere circ_sphere[NMAXCIRC_SPHERE]; /* circular scatterers on sphere */ diff --git a/global_ljones.c b/global_ljones.c index 7bb9adc..4d76d63 100644 --- a/global_ljones.c +++ b/global_ljones.c @@ -45,6 +45,7 @@ #define O_POOL_TABLE 3 /* obstacles around pockets of pool table */ #define O_HLINE_HOLE_SPOKES 181 /* tips of spokes for S_HLINE_HOLE_SPOKES segment pattern */ #define O_CIRCLE 4 /* one circle at the origin */ +#define O_FOUR_CIRCLES 5 /* four circles */ /* pattern of additional repelling segments */ #define S_RECTANGLE 0 /* segments forming a rectangle */ @@ -72,6 +73,10 @@ #define S_HLINE_HOLE_SPOKES 181 /* horizontal line with a hole in the bottom and extra spokes */ #define S_EXT_CIRCLE_RECT 19 /* particles outside a circle and a rectangle */ #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_AIRFOIL 24 /* exterior of an air foil */ /* particle interaction */ @@ -183,11 +188,19 @@ #define P_TYPE 5 /* colors represent type of particle */ #define P_DIRECTION 6 /* colors represent direction of velocity */ #define P_ANGULAR_SPEED 7 /* colors represent angular speed */ -#define P_DIRECT_ENERGY 8 /* hues represent direction, luminosity represents energy */ +#define P_DIRECT_ENERGY 8 /* hue represents direction, luminosity represents energy */ #define P_DIFF_NEIGHB 9 /* colors represent number of neighbours of different type */ #define P_THERMOSTAT 10 /* colors show which particles are coupled to the thermostat */ #define P_INITIAL_POS 11 /* colors depend on initial position of particle */ #define P_NUMBER 12 /* colors depend on particle number */ +#define P_EMEAN 13 /* averaged kinetic energy (with exponential damping) */ +#define P_DIRECT_EMEAN 14 /* averaged version of P_DIRECT_ENERGY */ +#define P_NOPARTICLE 15 /* particles are not drawn (only the links between them) */ + +/* Initial position dependence types */ + +#define IP_X 0 /* color depends on x coordinate of initial position */ +#define IP_Y 1 /* color depends on y coordinate of initial position */ /* Color schemes */ @@ -223,6 +236,7 @@ typedef struct double angle; /* angle of particle's "spin" */ short int active; /* circle is active */ double energy; /* dissipated energy */ + double emean; /* mean energy */ double vx; /* x velocity of particle */ double vy; /* y velocity of particle */ double omega; /* angular velocity of particle */ @@ -231,6 +245,7 @@ typedef struct double fx; /* x component of force on particle */ double fy; /* y component of force on particle */ double torque; /* torque on particle */ + double dirmean; /* time averaged direction */ int close_to_boundary; /* has value 1 if particle is close to a boundary */ short int thermostat; /* whether particle is coupled to thermostat */ int hashcell; /* hash cell in which particle is located */ @@ -331,5 +346,22 @@ typedef struct int color; /* color hue in case of different collisions */ } t_collision; + +typedef struct +{ + int nactive; /* number of active particles */ + double beta; /* inverse temperature */ + double mean_energy; /* mean energy */ + double krepel; /* force constant */ + double xmincontainer, xmaxcontainer; /* container size */ + double fboundary; /* boundary force */ + double pressure; /* pressure */ + double gravity; /* gravity */ + double radius; /* particle radius */ + double angle; /* orientation of obstacle */ + double omega; /* angular speed of obstacle */ + double bdry_fx, bdry_fy; /* components of boundary force */ +} t_lj_parameters; + int ncircles, nobstacles, nsegments, ngroups = 1, counter = 0; diff --git a/global_pdes.c b/global_pdes.c index 6674c62..4df65a6 100644 --- a/global_pdes.c +++ b/global_pdes.c @@ -82,6 +82,14 @@ #define D_TESLA 71 /* Tesla valve */ #define D_TESLA_FOUR 72 /* four Tesla valves */ +/* for wave_sphere.c */ + +#define D_LATITUDE 80 /* strip between two latitudes */ +#define D_SPHERE_CIRCLES 81 /* circles on the sphere */ +#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 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) */ @@ -104,6 +112,15 @@ #define C_RINGS_T 201 /* obstacles arranged in concentric rings, triangular lattice */ #define C_RINGS_SPIRAL 202 /* obstacles arranged on a "subflower" spiral, similar to C_GOLDEN_SPIRAL */ +#define C_HEX_BOTTOM 101 /* hex/triangular lattice in lower half */ +#define C_HEX_BOTTOM2 102 /* smaller hex/triangular lattice in lower half */ +#define C_SQUARE_BOTTOM 103 /* square lattice in lower half */ + +#define C_SPH_DODECA 30 /* dodecahedron (on sphere) */ +#define C_SPH_ICOSA 31 /* icosahedron (on sphere) */ +#define C_SPH_OCTA 32 /* octahedron (on sphere) */ +#define C_SPH_CUBE 33 /* cube (on sphere) */ + #define C_ONE 97 /* one single circle, as for Sinai */ #define C_TWO 98 /* two concentric circles of different type */ #define C_NOTHING 99 /* no circle at all, for comparisons */ diff --git a/heat.c b/heat.c index 57e2093..55eb8f6 100644 --- a/heat.c +++ b/heat.c @@ -211,7 +211,7 @@ #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) */ /* end of constants only used by sub_wave and sub_maze */ #include "global_pdes.c" diff --git a/lennardjones.c b/lennardjones.c index 5848c39..75014c3 100644 --- a/lennardjones.c +++ b/lennardjones.c @@ -62,8 +62,6 @@ #define INITXMAX 2.0 /* x interval for initial condition */ #define INITYMIN -1.1 #define INITYMAX 1.1 /* y interval for initial condition */ -// #define INITYMAX 7.0 /* y interval for initial condition */ -// #define INITYMAX 9.0 /* y interval for initial condition */ #define ADDXMIN -1.97 #define ADDXMAX -0.8 /* x interval for adding particles */ @@ -84,11 +82,11 @@ #define ADD_INITIAL_PARTICLES 0 /* set to 1 to add a second type of particles */ #define CIRCLE_PATTERN_B 1 /* pattern of circles for additional particles */ -#define ADD_FIXED_OBSTACLES 1 /* set to 1 do add fixed circular obstacles */ +#define ADD_FIXED_OBSTACLES 0 /* set to 1 do add fixed circular obstacles */ #define OBSTACLE_PATTERN 4 /* pattern of obstacles, see list in global_ljones.c */ -#define ADD_FIXED_SEGMENTS 0 /* set to 1 to add fixed segments as obstacles */ -#define SEGMENT_PATTERN 151 /* pattern of repelling segments, 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 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 */ @@ -101,8 +99,8 @@ #define CENTER_PY 0 /* set to 1 to center vertical momentum */ #define CENTER_PANGLE 0 /* set to 1 to center angular momentum */ -#define INTERACTION 2 /* particle interaction, see list in global_ljones.c */ -#define INTERACTION_B 2 /* particle interaction for second type of particle, see list in global_ljones.c */ +#define INTERACTION 1 /* particle interaction, see list in global_ljones.c */ +#define INTERACTION_B 1 /* particle interaction for second type of particle, see list in global_ljones.c */ #define 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 */ @@ -112,18 +110,21 @@ #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.8 /* parameter controlling the dimensions of domain */ -#define MU 0.015 /* parameter controlling radius of particles */ -#define MU_B 0.015 /* parameter controlling radius of particles of second type */ -#define NPOLY 25 /* number of sides of polygon */ +#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 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 */ #define MDEPTH 4 /* 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 140 /* number of grid point for grid of disks */ -#define NGRIDY 70 /* number of grid point for grid of disks */ +// #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 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 */ @@ -136,10 +137,10 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 2800 /* number of frames of movie */ -// #define NSTEPS 200 /* number of frames of movie */ -#define NVID 60 /* number of iterations between images displayed on screen */ -#define NSEG 250 /* number of segments of boundary */ +#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 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 */ #define BOUNDARY_WIDTH 1 /* width of particle boundary */ @@ -159,7 +160,7 @@ /* Plot type, see list in global_ljones.c */ -#define PLOT 4 +#define PLOT 13 #define PLOT_B 11 /* plot type for second movie */ #define DRAW_BONDS 1 /* set to 1 to draw bonds between neighbours */ @@ -168,12 +169,17 @@ #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 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 */ /* Color schemes */ -#define COLOR_PALETTE 18 /* Color palette, see list in global_ljones.c */ +#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 10 /* Color palette for angle representation */ -#define COLOR_PALETTE_INITIAL_POS 10 /* Color palette for initial position representation */ +#define COLOR_PALETTE_DIRECTION 17 /* Color palette for direction representation */ +#define COLOR_PALETTE_INITIAL_POS 0 /* Color palette for initial position representation */ #define BLACK 1 /* background */ @@ -189,9 +195,16 @@ #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 PRINT_PARAMETERS 0 /* set to 1 to print certain parameters */ +#define PRINT_PARAMETERS 1 /* set to 1 to print certain parameters */ #define PRINT_TEMPERATURE 0 /* set to 1 to print current temperature */ +#define PRINT_ANGLE 1 /* set to 1 to print obstacle orientation */ +#define PRINT_OMEGA 1 /* 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 FORCE_FACTOR 0.1 /* factor controlling length of force vector */ /* particle properties */ @@ -199,9 +212,9 @@ #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 6000.0 /* energy of particle with hottest color */ -#define HUE_TYPE0 60.0 /* hue of particles of type 0 */ -#define HUE_TYPE1 280.0 /* hue of particles of type 1 */ +#define PARTICLE_EMAX 5000.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 */ #define HUE_TYPE3 200.0 /* hue of particles of type 3 */ @@ -215,7 +228,7 @@ #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 2.5 /* mass of particle of radius MU */ +#define PARTICLE_MASS 1.0 /* mass of particle of radius MU */ #define PARTICLE_MASS_B 1.0 /* mass of particle of radius MU */ #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 */ @@ -231,7 +244,9 @@ #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.5 /* radius in which to count neighbours */ +// #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 */ @@ -243,7 +258,7 @@ #define KSPRING_VICSEK 0.2 /* spring constant for I_VICSEK_SPEED interaction */ #define VICSEK_REPULSION 10.0 /* repulsion between particles in Vicsek model */ -#define ROTATION 1 /* set to 1 to include rotation of particles */ +#define ROTATION 0 /* set to 1 to include rotation of particles */ #define COUPLE_ANGLE_TO_THERMOSTAT 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 */ #define KTORQUE 50.0 /* force constant in angular dynamics */ @@ -275,7 +290,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.35 /* radius of obstacle for circle boundary conditions */ +#define OBSTACLE_RADIUS 0.25 /* 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 */ @@ -310,16 +325,13 @@ #define TRACER_PARTICLE_MASS 4.0 /* relative mass of tracer particle */ #define TRAJECTORY_WIDTH 3 /* width of tracer particle trajectory */ -#define ROTATE_BOUNDARY 0 /* set to 1 to rotate the repelling segments */ +#define ROTATE_BOUNDARY 1 /* 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 PERIOD_ROTATE_BOUNDARY 1000 /* period of rotating boundary */ -#define ROTATE_INITIAL_TIME 0 /* initial time without rotation */ -#define ROTATE_FINAL_TIME 100 /* final time without rotation */ -#define ROTATE_CHANGE_TIME 0.33 /* relative duration of acceleration/deceleration phases */ -#define OMEGAMAX 100.0 /* maximal rotation speed */ -#define PRINT_OMEGA 0 /* set to 1 to print angular speed */ -#define PRINT_PARTICLE_SPEEDS 0 /* set to 1 to print average speeds/momenta of particles */ -#define PRINT_SEGMENTS_SPEEDS 1 /* set to 1 to print velocity of moving segments */ +#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 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 */ @@ -362,6 +374,9 @@ #define DELTA_EKIN 2000.0 /* change of kinetic energy in reaction */ #define COLLISION_TIME 15 /* time during which collisions are shown */ +#define CHANGE_RADIUS 0 /* set to 1 to change particle radius during simulation */ +#define MU_RATIO 0.666666667 /* ratio by which to increase radius */ + #define PRINT_PARTICLE_NUMBER 0 /* set to 1 to print total number of particles */ #define PLOT_PARTICLE_NUMBER 0 /* set to 1 to make of plot of particle number over time */ #define PARTICLE_NB_PLOT_FACTOR 0.5 /* expected final number of particles over initial number */ @@ -392,8 +407,8 @@ #define FLOOR_OMEGA 0 /* set to 1 to limit particle momentum to PMAX */ #define PMAX 1000.0 /* maximal force */ -#define HASHX 100 /* size of hashgrid in x direction */ -#define HASHY 50 /* size of hashgrid in y direction */ +#define HASHX 120 /* size of hashgrid in x direction */ +#define HASHY 60 /* 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 */ @@ -407,6 +422,8 @@ #define NO_WRAP_BC ((BOUNDARY_COND != BC_PERIODIC)&&(BOUNDARY_COND != BC_PERIODIC_CIRCLE)&&(BOUNDARY_COND != BC_PERIODIC_TRIANGLE)&&(BOUNDARY_COND != BC_KLEIN)&&(BOUNDARY_COND != BC_PERIODIC_FUNNEL)&&(BOUNDARY_COND != BC_BOY)&&(BOUNDARY_COND != BC_GENUS_TWO)) #define PERIODIC_BC ((BOUNDARY_COND == BC_PERIODIC)||(BOUNDARY_COND == BC_PERIODIC_CIRCLE)||(BOUNDARY_COND == BC_PERIODIC_FUNNEL)||(BOUNDARY_COND == BC_PERIODIC_TRIANGLE)) #define TWO_OBSTACLES ((SEGMENT_PATTERN == S_TWO_CIRCLES_EXT)||(SEGMENT_PATTERN == S_TWO_ROCKETS)) +#define COMPUTE_EMEAN ((PLOT == P_EMEAN)||(PLOT_B == P_EMEAN)||(PLOT == P_DIRECT_EMEAN)||(PLOT_B == P_DIRECT_EMEAN)) +#define COMPUTE_DIRMEAN ((PLOT == P_DIRECT_EMEAN)||(PLOT_B == P_DIRECT_EMEAN)) double xshift = 0.0; /* x shift of shown window */ double xspeed = 0.0; /* x speed of obstacle */ @@ -422,6 +439,8 @@ double ysegments[2] = {SEGMENTS_Y0, SEGMENTS_Y0}; /* y coordinate of segments ( double vxsegments[2] = {SEGMENTS_VX0, SEGMENTS_VX0}; /* vx coordinate of segments (for option MOVE_BOUNDARY) */ double vysegments[2] = {SEGMENTS_VY0, SEGMENTS_VY0}; /* vy coordinate of segments (for option MOVE_BOUNDARY) */ int thermostat_on = 1; /* thermostat switch used when VARY_THERMOSTAT is on */ +double cosangle[NSEG+1]; +double sinangle[NSEG+1]; /* precomputed trig functions of angles to draw circles faster */ #define THERMOSTAT_ON ((THERMOSTAT)&&((!VARY_THERMOSTAT)||(thermostat_on))) @@ -667,12 +686,17 @@ int thermostat_schedule(int i) else return(0); } +double radius_schedule(int i) +{ + return(1.0 + (MU_RATIO - 1.0)*(double)i/(double)NSTEPS); +} + double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HASHX*HASHY], double qx[NMAXCIRCLES], double qy[NMAXCIRCLES], double qangle[NMAXCIRCLES], double px[NMAXCIRCLES], double py[NMAXCIRCLES], double pangle[NMAXCIRCLES], double beta, int *nactive, int *nsuccess, int *nmove, int initial_phase) { - double a, totalenergy = 0.0, damping; + double a, totalenergy = 0.0, damping, direction, dmean; static double b = 0.25*SIGMA*SIGMA*DT_PARTICLE/MU_XI, xi = 0.0; int j, move; @@ -691,6 +715,21 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA pangle[j] = particle[j].omega + 0.5*DT_PARTICLE*particle[j].torque; particle[j].energy = (px[j]*px[j] + py[j]*py[j])*particle[j].mass_inv; + + if (COMPUTE_EMEAN) + particle[j].emean = ERATIO*particle[j].emean + (1.0-ERATIO)*particle[j].energy; + + if (COMPUTE_DIRMEAN) + { + direction = argument(particle[j].vx, particle[j].vy); + dmean = particle[j].dirmean; + if (dmean < direction - PI) dmean += DPI; + else if (dmean > direction + PI) dmean -= DPI; + particle[j].dirmean = DRATIO*dmean + (1.0-DRATIO)*direction; + if (particle[j].dirmean < 0.0) particle[j].dirmean += DPI; + else if (particle[j].dirmean > DPI) particle[j].dirmean -= DPI; + } + if ((COUPLE_ANGLE_TO_THERMOSTAT)&&(particle[j].thermostat)) particle[j].energy += pangle[j]*pangle[j]*particle[j].inertia_moment_inv; @@ -1091,10 +1130,10 @@ void evolve_segment_groups(t_segment segment[NMAXSEGMENTS], int time, t_group_se 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, krepel = KREPEL, pos[2], prop, vx, beta = BETA, xi = 0.0, xmincontainer = BCXMIN, xmaxcontainer = BCXMAX, torque, torque_ij, fboundary = 0.0, pleft = 0.0, pright = 0.0, entropy[2], mean_energy, gravity = GRAVITY, speed_ratio, ymin, ymax, delta_energy, speed, ratio = 1.0, ratioc, cum_etot = 0.0, emean = 0.0; + 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, ymin, ymax, delta_energy, speed, ratio = 1.0, ratioc, cum_etot = 0.0, emean = 0.0, radius_ratio; double *qx, *qy, *px, *py, *qangle, *pangle, *pressure, *obstacle_speeds; int i, j, k, n, m, s, ij[2], i0, iplus, iminus, j0, jplus, jminus, p, q, p1, q1, p2, q2, total_neighbours = 0, - min_nb, max_nb, close, wrapx = 0, wrapy = 0, nactive = 0, nadd_particle = 0, nmove = 0, nsuccess = 0, + min_nb, max_nb, close, wrapx = 0, wrapy = 0, nadd_particle = 0, nmove = 0, nsuccess = 0, tracer_n[N_TRACER_PARTICLES], traj_position = 0, traj_length = 0, move = 0, old, m0, floor, nthermo, wall = 0, group, gshift, n_total_active = 0, ncollisions = 0; int *particle_numbers; @@ -1108,10 +1147,20 @@ void animation() t_group_data *group_speeds; t_collision *collisions; t_hashgrid *hashgrid; + t_lj_parameters params; char message[100]; ratioc = 1.0 - ratio; + /* parameter values, grouped in a structure to simplify parameter printing */ + params.beta = BETA; + params.krepel = KREPEL; + params.xmincontainer = BCXMIN; + params.xmaxcontainer = BCXMAX; + params.fboundary = 0.0; + params.gravity = GRAVITY; + params.radius = MU; + particle = (t_particle *)malloc(NMAXCIRCLES*sizeof(t_particle)); /* particles */ if (ADD_FIXED_OBSTACLES) obstacle = (t_obstacle *)malloc(NMAXOBSTACLES*sizeof(t_obstacle)); /* circular obstacles */ if (ADD_FIXED_SEGMENTS) @@ -1154,6 +1203,9 @@ void animation() group_speeds = (t_group_data *)malloc(ngroups*(INITIAL_TIME + NSTEPS)*sizeof(t_group_data)); } + /* initialise array of trig functions to speed up drawing particles */ + init_angles(); + /* initialise positions and radii of circles */ init_particle_config(particle); @@ -1174,7 +1226,7 @@ void animation() printf("Initializing configuration\n"); - nactive = initialize_configuration(particle, hashgrid, obstacle, px, py, pangle, tracer_n, segment); + params.nactive = initialize_configuration(particle, hashgrid, obstacle, px, py, pangle, tracer_n, segment); // xi = 0.0; @@ -1199,12 +1251,12 @@ void animation() { printf("Computing frame %d\n",i); - if (INCREASE_KREPEL) krepel = repel_schedule(i); - if (INCREASE_BETA) beta = temperature_schedule(i); + if (INCREASE_KREPEL) params.krepel = repel_schedule(i); + if (INCREASE_BETA) params.beta = temperature_schedule(i); if (DECREASE_CONTAINER_SIZE) { - xmincontainer = container_size_schedule(i); - if (SYMMETRIC_DECREASE) xmaxcontainer = -container_size_schedule(i); + params.xmincontainer = container_size_schedule(i); + if (SYMMETRIC_DECREASE) params.xmaxcontainer = -container_size_schedule(i); } if ((ROTATE_BOUNDARY)&&(!SMOOTH_ROTATION)) rotate_segments(segment, rotation_schedule(i)); if (VARY_THERMOSTAT) @@ -1235,7 +1287,7 @@ void animation() blank(); - fboundary = 0.0; + params.fboundary = 0.0; pleft = 0.0; pright = 0.0; if (RECORD_PRESSURES) for (j=0; j OBSTACLE_INITIAL_TIME)) evolve_segments(segment, i); @@ -1378,10 +1436,10 @@ void animation() if ((PARTIAL_THERMO_COUPLING)&&(i>N_T_AVERAGE)) { nthermo = partial_thermostat_coupling(particle, xshift + PARTIAL_THERMO_SHIFT, segment); - printf("%i particles coupled to thermostat out of %i active\n", nthermo, nactive); - mean_energy = compute_mean_energy(particle); + printf("%i particles coupled to thermostat out of %i active\n", nthermo, params.nactive); + params.mean_energy = compute_mean_energy(particle); } - else mean_energy = totalenergy/(double)ncircles; + else params.mean_energy = totalenergy/(double)ncircles; if (CENTER_PX) center_momentum(px); if (CENTER_PY) center_momentum(py); @@ -1434,9 +1492,9 @@ void animation() cum_etot += totalenergy; } - printf("Boundary force: %.3f\n", fboundary/(double)(ncircles*NVID)); + printf("Boundary force: %.3f\n", params.fboundary/(double)(ncircles*NVID)); if (RESAMPLE_Y) printf("%i succesful moves out of %i trials\n", nsuccess, nmove); - if (INCREASE_GRAVITY) printf("Gravity: %.3f\n", gravity); + if (INCREASE_GRAVITY) printf("Gravity: %.3f\n", params.gravity); total_neighbours = 0; min_nb = 100; @@ -1456,11 +1514,11 @@ void animation() if ((i > INITIAL_TIME)&&(REACTION_DIFFUSION)) { ncollisions = update_types(particle, collisions, ncollisions, particle_numbers, i - INITIAL_TIME - 1, &delta_energy); - if (EXOTHERMIC) beta *= 1.0/(1.0 + delta_energy/totalenergy); - nactive = 0; + if (EXOTHERMIC) params.beta *= 1.0/(1.0 + delta_energy/totalenergy); + params.nactive = 0; for (j=0; j ADD_TIME)&&((i - INITIAL_TIME - ADD_TIME)%ADD_PERIOD == 1)&&(i < NSTEPS - FINAL_NOADD_PERIOD)) @@ -1480,11 +1538,22 @@ void animation() nadd_particle = add_particles(particle, px, py, nadd_particle); } + /* change particle radius */ + if (CHANGE_RADIUS) + { + radius_ratio = radius_schedule(i+1)/radius_schedule(i); + printf("Particle radius factor %.5lg\t", radius_schedule(i+1)); + for (j=0; j= INITIAL_TIME)&&(DOUBLE_MOVIE)) { if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length); - draw_particles(particle, PLOT_B, beta, collisions, ncollisions); - draw_container(xmincontainer, xmaxcontainer, obstacle, segment, wall); - if (PRINT_PARAMETERS) - print_parameters(beta, mean_energy, krepel, xmaxcontainer - xmincontainer, - fboundary/(double)(ncircles*NVID), PRINT_LEFT, pressure, gravity); + draw_particles(particle, PLOT_B, params.beta, collisions, ncollisions); + draw_container(params.xmincontainer, params.xmaxcontainer, obstacle, segment, wall); + if (PRINT_PARAMETERS) print_parameters(params, PRINT_LEFT, pressure, 0); if (PLOT_SPEEDS) draw_speed_plot(group_speeds, i); if (PLOT_TRAJECTORIES) draw_trajectory_plot(group_speeds, i); if (BOUNDARY_COND == BC_EHRENFEST) print_ehrenfest_parameters(particle, pleft, pright); @@ -1560,7 +1626,6 @@ void animation() print_particle_types_number(particle, RD_TYPES); else print_particle_number(ncircles); } - if (PRINT_OMEGA) print_omega(angular_speed); else if (PRINT_PARTICLE_SPEEDS) print_particles_speeds(particle); else if (PRINT_SEGMENTS_SPEEDS) print_segment_group_speeds(segment_group); // print_segments_speeds(vxsegments, vysegments); @@ -1601,11 +1666,9 @@ void animation() { blank(); if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length); - draw_particles(particle, PLOT, beta, collisions, ncollisions); - draw_container(xmincontainer, xmaxcontainer, obstacle, segment, wall); - if (PRINT_PARAMETERS) - print_parameters(beta, mean_energy, krepel, xmaxcontainer - xmincontainer, - fboundary/(double)(ncircles*NVID), PRINT_LEFT, pressure, gravity); + draw_particles(particle, PLOT, params.beta, collisions, ncollisions); + draw_container(params.xmincontainer, params.xmaxcontainer, obstacle, segment, wall); + if (PRINT_PARAMETERS) print_parameters(params, PRINT_LEFT, pressure, 1); if (PLOT_SPEEDS) draw_speed_plot(group_speeds, i); if (PLOT_TRAJECTORIES) draw_trajectory_plot(group_speeds, i); if (BOUNDARY_COND == BC_EHRENFEST) print_ehrenfest_parameters(particle, pleft, pright); @@ -1615,7 +1678,6 @@ void animation() print_particle_types_number(particle, RD_TYPES); else print_particle_number(ncircles); } - if (PRINT_OMEGA) print_omega(angular_speed); else if (PRINT_PARTICLE_SPEEDS) print_particles_speeds(particle); else if (PRINT_SEGMENTS_SPEEDS) print_segment_group_speeds(segment_group); // print_segments_speeds(vxsegments, vysegments); @@ -1630,11 +1692,9 @@ void animation() if (DOUBLE_MOVIE) { if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length); - draw_particles(particle, PLOT_B, beta, collisions, ncollisions); - draw_container(xmincontainer, xmaxcontainer, obstacle, segment, wall); - if (PRINT_PARAMETERS) - print_parameters(beta, mean_energy, krepel, xmaxcontainer - xmincontainer, - fboundary/(double)(ncircles*NVID), PRINT_LEFT, pressure, gravity); + draw_particles(particle, PLOT_B, params.beta, collisions, ncollisions); + draw_container(params.xmincontainer, params.xmaxcontainer, obstacle, segment, wall); + if (PRINT_PARAMETERS) print_parameters(params, PRINT_LEFT, pressure, 0); if (PLOT_SPEEDS) draw_speed_plot(group_speeds, i); if (PLOT_TRAJECTORIES) draw_trajectory_plot(group_speeds, i); if (BOUNDARY_COND == BC_EHRENFEST) print_ehrenfest_parameters(particle, pleft, pright); @@ -1644,7 +1704,6 @@ void animation() print_particle_types_number(particle, RD_TYPES); else print_particle_number(ncircles); } - if (PRINT_OMEGA) print_omega(angular_speed); else if (PRINT_PARTICLE_SPEEDS) print_particles_speeds(particle); else if (PRINT_SEGMENTS_SPEEDS) print_segment_group_speeds(segment_group); // print_segments_speeds(vxsegments, vysegments); @@ -1660,45 +1719,45 @@ void animation() s = system("mv lj*.tif tif_ljones/"); } - nactive = 0; - for (j=0; j YMAX) particles[n].yc += YMIN - YMAX; +// else if (particles[n].yc < YMIN) particles[n].yc += YMAX - YMIN; particles[n].radius = MU; /* 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; @@ -1149,6 +1228,36 @@ void init_obstacle_config(t_obstacle obstacle[NMAXOBSTACLES]) obstacle[n].active = 1; nobstacles = 1; break; + }case (O_FOUR_CIRCLES): + { + n = 0; + + obstacle[n].xc = -1.5; + obstacle[n].yc = -0.5; + obstacle[n].radius = OBSTACLE_RADIUS; + obstacle[n].active = 1; + n++; + + obstacle[n].xc = -0.5; + obstacle[n].yc = 0.5; + obstacle[n].radius = OBSTACLE_RADIUS; + obstacle[n].active = 1; + n++; + + obstacle[n].xc = 0.5; + obstacle[n].yc = -0.5; + obstacle[n].radius = OBSTACLE_RADIUS; + obstacle[n].active = 1; + n++; + + obstacle[n].xc = 1.5; + obstacle[n].yc = 0.5; + obstacle[n].radius = OBSTACLE_RADIUS; + obstacle[n].active = 1; + n++; + + nobstacles = 4; + break; } default: @@ -1158,7 +1267,7 @@ void init_obstacle_config(t_obstacle obstacle[NMAXOBSTACLES]) } } -void add_rotated_angle_to_segments(double x1, double y1, double x2, double y2, double width, t_segment segment[NMAXSEGMENTS], int group) +void add_rotated_angle_to_segments(double x1, double y1, double x2, double y2, double width, int center, t_segment segment[NMAXSEGMENTS], int group) /* add four segments forming a rectangle, specified by two adjacent corners and width */ { double tx, ty, ux, uy, norm, x3, y3, x4, y4; @@ -1169,6 +1278,13 @@ void add_rotated_angle_to_segments(double x1, double y1, double x2, double y2, d norm = module2(tx, ty); tx = tx/norm; ty = ty/norm; + if (center) + { + x2 -= 0.5*width*ty; + y2 += 0.5*width*tx; + x1 -= 0.5*width*ty; + y1 += 0.5*width*tx; + } x3 = x2 + width*ty; y3 = y2 - width*tx; x4 = x1 + width*ty; @@ -1258,7 +1374,7 @@ void add_rectangle_to_segments(double x1, double y1, double x2, double y2, t_seg void add_circle_to_segments(double x, double y, double r, int nsegs, double angle0, t_segment segment[NMAXSEGMENTS], int group) -/* add segments forming a circle to linear obstacle configuration */ +/* add segments forming a circle/polygon to linear obstacle configuration */ { int i, n = nsegments, nplus, nminus; double angle; @@ -1269,10 +1385,10 @@ void add_circle_to_segments(double x, double y, double r, int nsegs, double angl { for (i=0; i= BCXMAX) return(0); if (x <= BCXMIN) return(0); @@ -2487,6 +2703,53 @@ int in_segment_region(double x, double y, t_segment segment[NMAXSEGMENTS]) if ((y > 1.0 - LAMBDA + padding)&&(vabs(x) < LAMBDA - padding)) return(1); return(0); } + case (S_POLYGON_EXT): + { + padding = 3.0*MU; + if (in_polygon(x, y, LAMBDA + padding, NPOLY, APOLY)) return(0); + else return(1); + } + case (S_WEDGE_EXT): + { + padding = 3.0*MU; + angle = AWEDGE*PID; + ca = cos(APOLY*PID); + sa = sin(APOLY*PID); + x1 = x*ca - y*sa; + y1 = x*sa + y*ca; + if (vabs(y1) - padding > (LAMBDA-x1)*sin(angle)/(1.0+cos(angle))) return(1); + if (vabs(y1) + padding < -x1*tan(angle)) return(1); + return(0); + } + case (S_MIXER): + { + padding = 1.5*MU; + r = module2(x,y); + if (r > LAMBDA + padding) return(1); + if (r < 2.0*padding) return(0); + for (i=0; i 0.0)&&(vabs(y1) < 0.025 + padding)) return(0); + } + return(1); + } + case (S_AIRFOIL): + { + if (vabs(x) > LAMBDA + 4.0*MU) return(1); + padding = 0.35; + angle = APOLY*PID; + ca = cos(angle); + sa = sin(angle); + x1 = x*ca - y*sa; + y1 = x*sa + y*ca; + y1 += 0.5*x1*x1; + return(x1*x1 + 25.0*y1*y1 > LAMBDA*LAMBDA*(1.0 + padding)); + } case (S_MAZE): { for (i=0; i 0.0) + { + hue = ENERGY_HUE_MIN + (ENERGY_HUE_MAX - ENERGY_HUE_MIN)*ej/PARTICLE_EMAX; + if (hue > ENERGY_HUE_MIN) hue = ENERGY_HUE_MIN; + if (hue < ENERGY_HUE_MAX) hue = ENERGY_HUE_MAX; + } + *radius = particle.radius; + *width = BOUNDARY_WIDTH; + break; + } + case (P_DIRECT_EMEAN): + { + hue = particle.dirmean + COLOR_HUESHIFT*PI; + 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; + else lum = 1.0; + *radius = particle.radius; + *width = BOUNDARY_WIDTH; + break; + } + case (P_NOPARTICLE): + { + hue = 0.0; + lum = 1.0; + *radius = particle.radius; + *width = BOUNDARY_WIDTH; + break; + } } switch (plot) { case (P_KINETIC): { - hsl_to_rgb_turbo(hue, 0.9, 0.5, rgb); - hsl_to_rgb_turbo(hue, 0.9, 0.5, rgbx); - hsl_to_rgb_turbo(hue, 0.9, 0.5, rgby); +// hsl_to_rgb_turbo(hue, 0.9, 0.5, rgb); +// hsl_to_rgb_turbo(hue, 0.9, 0.5, rgbx); +// hsl_to_rgb_turbo(hue, 0.9, 0.5, rgby); + hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_EKIN); + hsl_to_rgb_palette(hue, 0.9, 0.5, rgbx, COLOR_PALETTE_EKIN); + hsl_to_rgb_palette(hue, 0.9, 0.5, rgby, COLOR_PALETTE_EKIN); break; } case (P_BONDS): @@ -3733,9 +4034,9 @@ void compute_particle_colors(t_particle particle, int plot, double rgb[3], doubl } case (P_DIRECTION): { - hsl_to_rgb_twilight(hue, 0.9, 0.5, rgb); - hsl_to_rgb_twilight(hue, 0.9, 0.5, rgbx); - hsl_to_rgb_twilight(hue, 0.9, 0.5, rgby); + hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_DIRECTION); + hsl_to_rgb_palette(hue, 0.9, 0.5, rgbx, COLOR_PALETTE_DIRECTION); + hsl_to_rgb_palette(hue, 0.9, 0.5, rgby, COLOR_PALETTE_DIRECTION); break; } case (P_ANGLE): @@ -3772,6 +4073,26 @@ void compute_particle_colors(t_particle particle, int plot, double rgb[3], doubl hsl_to_rgb_palette(hue, 0.9, 0.5, rgby, COLOR_PALETTE_INITIAL_POS); break; } + case (P_EMEAN): + { + hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_EKIN); + hsl_to_rgb_palette(hue, 0.9, 0.5, rgbx, COLOR_PALETTE_EKIN); + hsl_to_rgb_palette(hue, 0.9, 0.5, rgby, COLOR_PALETTE_EKIN); + break; + } + case (P_DIRECT_EMEAN): + { + hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_DIRECTION); + hsl_to_rgb_palette(hue, 0.9, 0.5, rgbx, COLOR_PALETTE_DIRECTION); + hsl_to_rgb_palette(hue, 0.9, 0.5, rgby, COLOR_PALETTE_DIRECTION); + for (i=0; i<3; i++) + { + rgb[i] *= lum; + rgbx[i] *= lum; + rgby[i] *= lum; + } + break; + } default: { hsl_to_rgb(hue, 0.9, 0.5, rgb); @@ -4322,8 +4643,12 @@ void draw_one_particle(t_particle particle, double xc, double yc, double radius, if ((particle.interaction == I_LJ_QUADRUPOLE)||(particle.interaction == I_LJ_DIPOLE)||(particle.interaction == I_VICSEK)||(particle.interaction == I_VICSEK_REPULSIVE)||(particle.interaction == I_VICSEK_SPEED)||(particle.interaction == I_VICSEK_SHARK)) draw_colored_rhombus(xc1, yc1, radius, angle + APOLY*PID, rgb); - else if (cont) draw_colored_polygon(xc1, yc1, radius, nsides, angle + APOLY*PID, rgb); - + else if (cont) + { + if (nsides == NSEG) draw_colored_circle_precomp(xc1, yc1, radius, rgb); + else draw_colored_polygon(xc1, yc1, radius, nsides, angle + APOLY*PID, rgb); + } + /* draw crosses on particles of second type */ if ((TWO_TYPES)&&(DRAW_CROSS)) if (particle.type == 1) @@ -4352,7 +4677,11 @@ void draw_one_particle(t_particle particle, double xc, double yc, double radius, if ((particle.interaction == I_LJ_QUADRUPOLE)||(particle.interaction == I_LJ_DIPOLE)||(particle.interaction == I_VICSEK)||(particle.interaction == I_VICSEK_REPULSIVE)||(particle.interaction == I_VICSEK_SPEED)||(particle.interaction == I_VICSEK_SHARK)) draw_rhombus(xc1, yc1, radius, angle + APOLY*PID); - else if (cont) draw_polygon(xc1, yc1, radius, nsides, angle + APOLY*PID); + else if (cont) + { + if (nsides == NSEG) draw_circle_precomp(xc1, yc1, radius); + else draw_polygon(xc1, yc1, radius, nsides, angle + APOLY*PID); + } if (particle.interaction == I_LJ_WATER) for (wsign = -1; wsign <= 1; wsign+=2) { @@ -4470,6 +4799,7 @@ void draw_particles(t_particle particle[NMAXCIRCLES], int plot, double beta, t_c char message[100]; // if (!TRACER_PARTICLE) blank(); + if (plot == P_NOPARTICLE) blank(); glColor3f(1.0, 1.0, 1.0); /* show region of partial thermostat */ @@ -4536,19 +4866,17 @@ void draw_particles(t_particle particle[NMAXCIRCLES], int plot, double beta, t_c angle = particle[j].angle + APOLY*DPI; - draw_one_particle(particle[j], particle[j].xc, particle[j].yc, radius, angle, nsides, width, rgb); - /* in case of periodic b.c., draw translates of particles */ - if (PERIODIC_BC) + if ((PERIODIC_BC)&&(plot != P_NOPARTICLE)) { x1 = particle[j].xc; y1 = particle[j].yc; - for (i=-2; i<3; i++) + for (i=-1; i<2; i++) for (k=-1; k<2; k++) - draw_one_particle(particle[j], x1 + (double)i*(BCXMAX - BCXMIN), y1 + (double)k*(BCYMAX - BCYMIN), radius, angle, nsides, width, rgb); + draw_one_particle(particle[j], x1 + (double)i*(BCXMAX - BCXMIN), y1 + (double)k*(BCYMAX - BCYMIN), radius, angle, nsides, width, rgb); } - else if (BOUNDARY_COND == BC_KLEIN) + else if ((BOUNDARY_COND == BC_KLEIN)&&(plot != P_NOPARTICLE)) { x1 = particle[j].xc; y1 = particle[j].yc; @@ -4563,7 +4891,7 @@ void draw_particles(t_particle particle[NMAXCIRCLES], int plot, double beta, t_c radius, angle1, nsides, width, rgb); } } - else if (BOUNDARY_COND == BC_BOY) + else if ((BOUNDARY_COND == BC_BOY)&&(plot != P_NOPARTICLE)) { x1 = particle[j].xc; y1 = particle[j].yc; @@ -4584,7 +4912,7 @@ void draw_particles(t_particle particle[NMAXCIRCLES], int plot, double beta, t_c sign*(y1 + (double)k*(BCYMAX - BCYMIN)), radius, angle1, nsides, width, rgb); } } - else if (BOUNDARY_COND == BC_GENUS_TWO) + else if ((BOUNDARY_COND == BC_GENUS_TWO)&&(plot != P_NOPARTICLE)) { x1 = particle[j].xc; y1 = particle[j].yc; @@ -4622,6 +4950,10 @@ void draw_particles(t_particle particle[NMAXCIRCLES], int plot, double beta, t_c draw_one_particle(particle[j], x, y, radius, angle, nsides, width, rgb); } } + else if (plot != P_NOPARTICLE) + draw_one_particle(particle[j], particle[j].xc, particle[j].yc, radius, angle, nsides, width, rgb); + + } // /* draw spin vectors */ @@ -4658,10 +4990,10 @@ void draw_container(double xmin, double xmax, t_obstacle obstacle[NMAXOBSTACLES] glLineWidth(CONTAINER_WIDTH); hsl_to_rgb(300.0, 0.1, 0.5, rgb); for (i = 0; i < nobstacles; i++) - draw_colored_circle(obstacle[i].xc - xtrack, obstacle[i].yc - ytrack, obstacle[i].radius, NSEG, rgb); + draw_colored_circle_precomp(obstacle[i].xc - xtrack, obstacle[i].yc - ytrack, obstacle[i].radius, rgb); glColor3f(1.0, 1.0, 1.0); for (i = 0; i < nobstacles; i++) - draw_circle(obstacle[i].xc - xtrack, obstacle[i].yc - ytrack, obstacle[i].radius, NSEG); + draw_circle_precomp(obstacle[i].xc - xtrack, obstacle[i].yc - ytrack, obstacle[i].radius); } if (ADD_FIXED_SEGMENTS) { @@ -4710,9 +5042,9 @@ void draw_container(double xmin, double xmax, t_obstacle obstacle[NMAXOBSTACLES] if (CENTER_VIEW_ON_OBSTACLE) x = 0.0; else x = xmin + (double)i*(OBSXMAX - OBSXMIN); - draw_colored_circle(x, 0.0, OBSTACLE_RADIUS, NSEG, rgb); + draw_colored_circle_precomp(x, 0.0, OBSTACLE_RADIUS, rgb); glColor3f(1.0, 1.0, 1.0); - draw_circle(x, 0.0, OBSTACLE_RADIUS, NSEG); + draw_circle_precomp(x, 0.0, OBSTACLE_RADIUS); glColor3f(0.0, 0.0, 0.0); sprintf(message, "Mach %.3f", xspeed/20.0); @@ -4730,9 +5062,9 @@ void draw_container(double xmin, double xmax, t_obstacle obstacle[NMAXOBSTACLES] if (CENTER_VIEW_ON_OBSTACLE) x = 0.0; else x = xmin + (double)i*(OBSXMAX - OBSXMIN); - draw_colored_circle(x, 0.0, OBSTACLE_RADIUS, NSEG, rgb); + draw_colored_circle_precomp(x, 0.0, OBSTACLE_RADIUS, rgb); glColor3f(1.0, 1.0, 1.0); - draw_circle(x, 0.0, OBSTACLE_RADIUS, NSEG); + draw_circle_precomp(x, 0.0, OBSTACLE_RADIUS); glColor3f(0.0, 0.0, 0.0); sprintf(message, "Mach %.2f", xspeed/20.0); @@ -4774,9 +5106,9 @@ void draw_container(double xmin, double xmax, t_obstacle obstacle[NMAXOBSTACLES] for (j=-1; j<2; j+=2) { - draw_colored_circle(x, (double)j*(FUNNEL_WIDTH + OBSTACLE_RADIUS), OBSTACLE_RADIUS, NSEG, rgb); + draw_colored_circle_precomp(x, (double)j*(FUNNEL_WIDTH + OBSTACLE_RADIUS), OBSTACLE_RADIUS, rgb); glColor3f(1.0, 1.0, 1.0); - draw_circle(x, (double)j*(FUNNEL_WIDTH + OBSTACLE_RADIUS), OBSTACLE_RADIUS, NSEG); + draw_circle_precomp(x, (double)j*(FUNNEL_WIDTH + OBSTACLE_RADIUS), OBSTACLE_RADIUS); } glColor3f(0.0, 0.0, 0.0); @@ -4869,12 +5201,77 @@ void draw_container(double xmin, double xmax, t_obstacle obstacle[NMAXOBSTACLES] } -void print_parameters(double beta, double temperature, double krepel, double lengthcontainer, double boundary_force, short int left, double pressure[N_PRESSURES], double gravity) +void print_omega(double angle, double angular_speed, double fx, double fy) +{ + char message[100]; + double rgb[3], y1, frac, absa; + static double xleftbox, xlefttext, xrightbox, xrighttext, y = YMAX - 0.1, ymin = YMIN + 0.05; + static int first = 1; + + if (first) + { + xrightbox = XMAX - 0.54; + xrighttext = xrightbox - 0.48; + first = 0; + } + + y1 = y; + + if (PRINT_ANGLE) + { + erase_area_hsl(xrightbox, y + 0.025, 0.42, 0.05, 0.0, 0.9, 0.0); + glColor3f(1.0, 1.0, 1.0); + angle = angle*360.0/DPI; + absa = vabs(angle); + frac = absa - (double)((int)absa); + sprintf(message, "Angle = %3d.%d degrees", (int)absa, (int)(10.0*frac)); + write_text(xrighttext + 0.1, y, message); + y1 -= 0.1; + } + if (PRINT_OMEGA) + { + 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, "Angular speed = %.4f", angular_speed); + write_text(xrighttext + 0.1, y1, message); + y1 -= 0.1; + } + if (PRINT_SEGMENTS_FORCE) + { + 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, "Fx = %.4f", fx); + 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 = %.4f", fy); + write_text(xrighttext + 0.1, y1, message); + } +} + +void compute_segments_force(t_lj_parameters *params, t_segment segment[NMAXSEGMENTS]) +{ + int i; + double fx = 0.0, fy = 0.0; + + for (i=0; ibdry_fx = fx; + params->bdry_fy = fy; + +} + +void print_parameters(t_lj_parameters params, short int left, double pressure[N_PRESSURES], short int refresh) { char message[100]; int i, j, k; - double density, hue, rgb[3], logratio, x, y, meanpress[N_PRESSURES], phi, sphi, dphi, pprint, mean_temp; - static double xbox, xtext, xmid, xmidtext, xxbox, xxtext, pressures[N_P_AVERAGE], meanpressure = 0.0, maxpressure = 0.0; + double density, hue, rgb[3], logratio, x, y, meanpress[N_PRESSURES], phi, sphi, dphi, pprint, mean_temp, lengthcontainer, boundary_force, fx, fy, r1, r2; + 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; @@ -4907,10 +5304,17 @@ void print_parameters(double beta, double temperature, double krepel, double len i_pressure = 0; i_temp = 0; for (i=0; i 1.0) logratio = 1.0; else if (logratio < 0.0) logratio = 0.0; if (BETA_FACTOR > 1.0) hue = PARTICLE_HUE_MAX - (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*logratio; @@ -4954,7 +5358,7 @@ void print_parameters(double beta, double temperature, double krepel, double len else erase_area_hsl_turbo(xmid + 0.1, y + 0.025*scale, 0.45*scale, 0.05*scale, hue, 0.9, 0.5); if ((hue < 90)||(hue > 270)) glColor3f(1.0, 1.0, 1.0); else glColor3f(0.0, 0.0, 0.0); - sprintf(message, "Temperature %.2f", 1.0/beta); + sprintf(message, "Temperature %.2f", 1.0/params.beta); if (PRINT_LEFT) write_text(xtext, y, message); else write_text(xmidtext, y, message); // y -= 0.1; @@ -4974,7 +5378,7 @@ void print_parameters(double beta, double temperature, double krepel, double len erase_area_hsl(xmid, y + 0.025*scale, 0.37*scale, 0.05*scale, 0.0, 0.9, 0.0); glColor3f(1.0, 1.0, 1.0); - sprintf(message, "Temperature %.2f", temperature); + sprintf(message, "Temperature %.2f", params.mean_energy); write_text(xmidtext, y, message); erase_area_hsl(xxbox, y + 0.025*scale, 0.37*scale, 0.05*scale, 0.0, 0.9, 0.0); @@ -4987,7 +5391,7 @@ void print_parameters(double beta, double temperature, double krepel, double len { erase_area_hsl(xbox, y + 0.025*scale, 0.22*scale, 0.05*scale, 0.0, 0.9, 0.0); glColor3f(1.0, 1.0, 1.0); - sprintf(message, "Force %.0f", krepel); + sprintf(message, "Force %.0f", params.krepel); write_text(xtext + 0.28, y, message); } @@ -5023,8 +5427,8 @@ void print_parameters(double beta, double temperature, double krepel, double len if ((PARTIAL_THERMO_COUPLING)&&(!INCREASE_BETA)&&(!EXOTHERMIC)) { - printf("Temperature %i in average: %.3lg\n", i_temp, temperature); - temp[i_temp] = temperature; + printf("Temperature %i in average: %.3lg\n", i_temp, params.mean_energy); + temp[i_temp] = params.mean_energy; i_temp++; if (i_temp >= N_T_AVERAGE) i_temp = 0; @@ -5045,9 +5449,44 @@ void print_parameters(double beta, double temperature, double krepel, double len { 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); - sprintf(message, "Gravity %.2f", gravity/GRAVITY); + sprintf(message, "Gravity %.2f", params.gravity/GRAVITY); write_text(xmidtext + 0.1, y, message); - } + } + + if (CHANGE_RADIUS) + { + erase_area_hsl(xmid, y + 0.025*scale, 0.3*scale, 0.05*scale, 0.0, 0.9, 0.0); + glColor3f(1.0, 1.0, 1.0); + sprintf(message, "Radius %.4f", params.radius); + write_text(xmidtext + 0.05, y, message); + } + + if (PRINT_SEGMENTS_FORCE) + { + glColor3f(1.0, 1.0, 1.0); + if (refresh) + { + fx = 0.01*params.bdry_fx/(double)params.nactive; + fy = 0.01*params.bdry_fy/(double)params.nactive; + /* average boundary force */ + mean_fx = r2*mean_fx + r1*fx; + mean_fy = r2*mean_fy + r1*fy; + } + + if ((PRINT_ANGLE)||(PRINT_OMEGA)) + draw_arrow(0.0, 0.0, FORCE_FACTOR*mean_fx, FORCE_FACTOR*mean_fy, 15.0, 0.1); + else + { + erase_area_hsl(xmid, 0.0, 0.2*scale, 0.05*scale, 0.0, 0.9, 0.0); + glColor3f(1.0, 1.0, 1.0); + sprintf(message, "Fy = %.2f", mean_fy); + write_text(xmidtext + 0.15*scale, -0.02*scale, message); + if (mean_fx*mean_fx + mean_fy*mean_fy > 5.0*FORCE_FACTOR) + draw_arrow(0.0, 0.0, FORCE_FACTOR*mean_fx, FORCE_FACTOR*mean_fy, 15.0, 0.1); + } + } + if ((PRINT_ANGLE)||(PRINT_OMEGA)) print_omega(params.angle, params.omega, mean_fx, mean_fy); + } @@ -5266,29 +5705,6 @@ void print_entropy(double entropy[2]) } -void print_omega(double angular_speed) -{ - char message[100]; - double rgb[3]; - static double xleftbox, xlefttext, xrightbox, xrighttext, y = YMAX - 0.1, ymin = YMIN + 0.05; - static int first = 1; - - if (first) - { - xrightbox = XMAX - 0.39; - xrighttext = xrightbox - 0.55; - first = 0; - } - - - erase_area_hsl(xrightbox, y + 0.025, 0.35, 0.05, 0.0, 0.9, 0.0); - glColor3f(1.0, 1.0, 1.0); - sprintf(message, "Angular speed = %.4f", angular_speed); -// sprintf(message, "Angular speed = %.4f", DPI*angular_speed*25.0/(double)(PERIOD_ROTATE_BOUNDARY)); - write_text(xrighttext + 0.1, y, message); - -} - void print_segments_speeds(double vx[2], double vy[2]) { char message[100]; @@ -5500,11 +5916,12 @@ double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacl - if ((MOVE_BOUNDARY)||(MOVE_SEGMENT_GROUPS)) + if ((MOVE_BOUNDARY)||(MOVE_SEGMENT_GROUPS)||(PRINT_SEGMENTS_FORCE)) { segment[i].fx -= f*segment[i].nx; segment[i].fy -= f*segment[i].ny; 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 ((VICSEK_INT)&&(vabs(distance) < 1.5*r)) @@ -6089,6 +6506,8 @@ t_segment segment[NMAXSEGMENTS]) particle[i].diff_neighb = 0; particle[i].thermostat = 1; particle[i].close_to_boundary = 0; + particle[i].emean = 0.0; + particle[i].dirmean = 0.0; // particle[i].energy = 0.0; // y = particle[i].yc; @@ -6126,6 +6545,8 @@ t_segment segment[NMAXSEGMENTS]) } particle[i].energy = (particle[i].vx*particle[i].vx + particle[i].vy*particle[i].vy)*particle[i].mass_inv; + particle[i].emean = particle[i].energy; + particle[i].dirmean = 0.0; px[i] = particle[i].vx; py[i] = particle[i].vy; @@ -6145,7 +6566,21 @@ t_segment segment[NMAXSEGMENTS]) pangle[i] = particle[i].omega; if ((PLOT == P_INITIAL_POS)||(PLOT_B == P_INITIAL_POS)) - particle[i].color_hue = 360.0*(particle[i].yc - INITYMIN)/(INITYMAX - INITYMIN); + { + switch (INITIAL_POS_TYPE) { + case (IP_X): + { + particle[i].color_hue = 360.0*(particle[i].xc - INITXMIN)/(INITXMAX - INITXMIN); + break; + } + case (IP_Y): + { + particle[i].color_hue = 360.0*(particle[i].yc - INITYMIN)/(INITYMAX - INITYMIN); + break; + } + } + + } else if ((PLOT == P_NUMBER)||(PLOT_B == P_NUMBER)) particle[i].color_hue = 360.0*(double)(i%N_PARTICLE_COLORS)/(double)N_PARTICLE_COLORS; } @@ -6157,6 +6592,8 @@ t_segment segment[NMAXSEGMENTS]) particle[i].neighb = 0; particle[i].thermostat = 0; particle[i].energy = 0.0; + particle[i].emean = 0.0; + 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].vx = 0.0; @@ -7825,7 +8262,7 @@ void draw_trajectory_plot(t_group_data *group_speeds, int i) x1 = x2; y1 = y2; } - if (i>0) draw_colored_circle(x1, y1, 0.015, NSEG, rgb); + if (i>0) draw_colored_circle_precomp(x1, y1, 0.015, rgb); } glColor3f(1.0, 1.0, 1.0); diff --git a/sub_rde.c b/sub_rde.c index 7b78421..c533574 100644 --- a/sub_rde.c +++ b/sub_rde.c @@ -794,7 +794,93 @@ double init_linear_blob(double x, double y, double vx, double vy, double amp, do } +double init_circular_vibration(double x, double y, double v, double radius, double amp, double var, double omega, double min, double *phi[NFIELDS], short int xy_in[NX*NY]) +/* initialise gaussian wave height for shallow water equation */ +{ + int i, j; + double xy[2], r, angle, a, height; + + for (i=0; i XMAX) xy[0] += XMIN - XMAX; + + xy[1] *= (double)NY/(double)(NY-2*DPOLE); + + wsphere[i*NY+j].x2d = xy[0]; + wsphere[i*NY+j].y2d = xy[1]; + } + + /* cotangent, taking care of not dividing by zero */ + for (j=1; j nmaxpixels) + { + printf("Image too large, increase nmaxpixels in choose_colors()\n"); + exit(0); + } + + /* shift due to min/max latitudes of image */ + sshift = 0 + DPOLE; + nshift = 0 + DPOLE; + + /* read rgb values */ + for (j=0; j nx-1) ii -= nx; +// jj = (int)(-ry*(double)j + cy); +// jj = (int)(ry*(double)(NY-nshift - j)) + sshift; + jj = (int)(ry*(double)(NY-nshift - j)); + if (jj > ny-1) jj = ny-1; + if (jj < 0) jj = 0; + wsphere[i*NY+j].r = (double)rgb_values[3*(jj*nx+ii)]*cratio; + wsphere[i*NY+j].g = (double)rgb_values[3*(jj*nx+ii)+1]*cratio; + wsphere[i*NY+j].b = (double)rgb_values[3*(jj*nx+ii)+2]*cratio; + + /* decide which points are in the Sea */ + diff = iabs(rgb_values[3*(jj*nx+ii)] - 10); + diff += iabs(rgb_values[3*(jj*nx+ii)+1] - 10); + diff += iabs(rgb_values[3*(jj*nx+ii)+2] - 51); + wsphere[i*NY+j].indomain = (diff < 10); + } + + free(rgb_values); + fclose(image_file); +} + + +int ij_to_sphere(int i, int j, double r, t_wave_sphere wsphere[NX*NY], double xyz[3]) +/* convert spherical to rectangular coordinates */ +{ + double pscal, newr; + static double norm_observer; + static int first = 1; + + if (first) + { + norm_observer = sqrt(observer[0]*observer[0] + observer[1]*observer[1] + observer[2]*observer[2]); + first = 0; + } + + newr = wsphere[i*NY+j].radius; + + xyz[0] = wsphere[i*NY+j].x; + xyz[1] = wsphere[i*NY+j].y; + xyz[2] = wsphere[i*NY+j].z; + + pscal = xyz[0]*observer[0] + xyz[1]*observer[1] + xyz[2]*observer[2]; + + xyz[0] *= newr; + xyz[1] *= newr; + xyz[2] *= newr; + + return(pscal/norm_observer > COS_VISIBLE); +} + +int init_circle_sphere(t_circles_sphere *circles, int circle_pattern) +/* initialise circles on sphere */ +/* for billiard shape D_SPHERE_CIRCLES */ +{ + int ncircles, k; + double alpha, beta, gamma; + + switch (circle_pattern) { + case (C_SPH_DODECA): + { + alpha = acos(sqrt(5.0)/3.0); + beta = acos(1.0/3.0); + gamma = asin(sqrt(3.0/8.0)); + + circles[0].theta = 0.0; + circles[0].phi = 0.0; + + for (k=0; k<3; k++) + { + circles[1+k].theta = alpha; + circles[1+k].phi = (double)k*DPI/3.0; + } + + for (k=0; k<3; k++) + { + circles[4+k].theta = beta; + circles[4+k].phi = (double)k*DPI/3.0 + gamma; + circles[7+k].theta = beta; + circles[7+k].phi = (double)k*DPI/3.0 - gamma; + } + + for (k=0; k<3; k++) + { + circles[10+k].theta = PI - beta; + circles[10+k].phi = (double)k*DPI/3.0 + PI/3.0 + gamma; + circles[13+k].theta = PI - beta; + circles[13+k].phi = (double)k*DPI/3.0 + PI/3.0 - gamma; + } + + for (k=0; k<3; k++) + { + circles[16+k].theta = PI - alpha; + circles[16+k].phi = (double)k*DPI/3.0 + PI/3.0; + } + + circles[19].theta = PI; + circles[19].phi = 0.0; + + for (k=0; k<20; k++) + { + circles[k].radius = MU; + circles[k].x = sin(circles[k].theta)*cos(circles[k].phi); + circles[k].y = sin(circles[k].theta)*sin(circles[k].phi); + circles[k].z = cos(circles[k].theta); + } + + ncircles = 20; + break; + } + case (C_SPH_ICOSA): + { + alpha = acos(1.0/sqrt(5.0)); + + circles[0].theta = 0.0; + circles[0].phi = 0.0; + + for (k=0; k<5; k++) + { + circles[1+k].theta = alpha; + circles[1+k].phi = (double)k*DPI/5.0; + } + + for (k=0; k<5; k++) + { + circles[6+k].theta = PI - alpha; + circles[6+k].phi = (double)k*DPI/5.0 + PI/5.0; + } + + circles[11].theta = PI; + circles[11].phi = 0.0; + + for (k=0; k<12; k++) + { + circles[k].radius = MU; + circles[k].x = sin(circles[k].theta)*cos(circles[k].phi); + circles[k].y = sin(circles[k].theta)*sin(circles[k].phi); + circles[k].z = cos(circles[k].theta); + } + + ncircles = 12; + break; + } + case (C_SPH_OCTA): + { + circles[0].theta = 0.0; + circles[0].phi = 0.0; + + for (k=0; k<4; k++) + { + circles[1+k].theta = PID; + circles[1+k].phi = (double)k*PID; + } + + circles[5].theta = PI; + circles[5].phi = 0.0; + + for (k=0; k<6; k++) + { + circles[k].radius = MU; + circles[k].x = sin(circles[k].theta)*cos(circles[k].phi); + circles[k].y = sin(circles[k].theta)*sin(circles[k].phi); + circles[k].z = cos(circles[k].theta); + } + + ncircles = 6; + break; + } + case (C_SPH_CUBE): + { + alpha = acos(1.0/3.0); + + circles[0].theta = 0.0; + circles[0].phi = 0.0; + + for (k=0; k<3; k++) + { + circles[1+k].theta = alpha; + circles[1+k].phi = (double)k*DPI/3.0; + circles[4+k].theta = PI - alpha; + circles[4+k].phi = (double)k*DPI/3.0 + PI/3.0; + } + + circles[7].theta = PI; + circles[7].phi = 0.0; + + for (k=0; k<8; k++) + { + circles[k].radius = MU; + circles[k].x = sin(circles[k].theta)*cos(circles[k].phi); + circles[k].y = sin(circles[k].theta)*sin(circles[k].phi); + circles[k].z = cos(circles[k].theta); + } + + ncircles = 8; + break; + } + default: + { + printf("Function init_circle_sphere not defined for this pattern \n"); + } + } + return(ncircles); +} + +int xy_in_billiard_sphere(int i, int j, t_wave_sphere wsphere[NX*NY]) +/* returns 1 if (x,y) represents a point in the billiard */ +{ + int k; + double pscal, dist, r, u, v, u1, v1; + static double cos_rot, sin_rot; + static int first = 1; + + if (first) + { + if (B_DOMAIN == D_SPHERE_EARTH) init_earth_map(wsphere); + else if ((B_DOMAIN == D_SPHERE_JULIA)||(B_DOMAIN == D_SPHERE_JULIA_INV)) + { + cos_rot = cos(JULIA_ROT*DPI/360.0); + sin_rot = sin(JULIA_ROT*DPI/360.0); + } + first = 0; + } + + switch (B_DOMAIN) { + case (D_NOTHING): + { + return(1); + } + case (D_LATITUDE): + { + return(vabs(wsphere[i*NY+j].theta - PID) < LAMBDA*PID); + } + case (D_SPHERE_CIRCLES): + { + for (k=0; k RMAX) r = RMAX; + wsphere[i*NY+j].radius = r; + } + + + if (SHADE_WAVE) + { + #pragma omp parallel for private(i,j,norm,pscal,deltar,deltai,deltaj,n) + for (i=0; i NX) ii -= NX; + + glVertex2i(ii, j); + glVertex2i(ii+1, j); + glVertex2i(ii+1, j+1); + glVertex2i(ii, j+1); + } + glEnd (); +} + +void draw_wave_sphere_ij(int i, int iplus, int j, int jplus, int jcolor, int movie, double phi[NX*NY], double psi[NX*NY], short int xy_in[NX*NY], t_wave wave[NX*NY], t_wave_sphere wsphere[NX*NY], int zplot, int cplot, int palette, int fade, double fade_value) +/* draw wave at simulation grid point (i,j) */ +{ + int k, l; + double xyz[3], ca; + + if ((TWOSPEEDS)||(xy_in[i*NY+j])) + glColor3f(wave[i*NY+jcolor].rgb[0], wave[i*NY+jcolor].rgb[1], wave[i*NY+jcolor].rgb[2]); + else + { + ca = wave[i*NY+j].cos_angle; + ca = (ca + 1.0)*0.4 + 0.2; + if (fade) ca *= fade_value; + if (B_DOMAIN == D_SPHERE_EARTH) + glColor3f(wsphere[i*NY+j].r*ca, wsphere[i*NY+j].g*ca, wsphere[i*NY+j].b*ca); + else glColor3f(COLOR_OUT_R*ca, COLOR_OUT_G*ca, COLOR_OUT_B*ca); + } + glBegin(GL_TRIANGLE_FAN); + if (ij_to_sphere(i, j, *wave[i*NY+j].p_zfield[movie], wsphere, xyz)) + draw_vertex_sphere(xyz); + if (ij_to_sphere(iplus, j, *wave[iplus*NY+j].p_zfield[movie], wsphere, xyz)) + draw_vertex_sphere(xyz); + if (ij_to_sphere(iplus, jplus, *wave[iplus*NY+j+1].p_zfield[movie], wsphere, xyz)) + draw_vertex_sphere(xyz); + if (ij_to_sphere(i, jplus, *wave[i*NY+j+1].p_zfield[movie], wsphere, xyz)) + draw_vertex_sphere(xyz); + glEnd (); +} + + +void draw_wave_sphere_3D(int movie, double phi[NX*NY], double psi[NX*NY], short int xy_in[NX*NY], t_wave wave[NX*NY], t_wave_sphere wsphere[NX*NY], int zplot, int cplot, int palette, int fade, double fade_value, int refresh) +{ + int i, j, imax, imin; + double observer_angle, angle2; + + blank(); + + if (refresh) + { + compute_wave_fields(phi, psi, xy_in, zplot, cplot, wave); + if (SHADE_3D) compute_light_angle_sphere(xy_in, wave, wsphere, movie); + else if (SHADE_2D) compute_light_angle_sphere_2d(xy_in, wave, movie); + compute_cfield_sphere(xy_in, cplot, palette, wave, fade, fade_value, movie); + } + + observer_angle = argument(observer[0], observer[1]); + if (observer_angle < 0.0) observer_angle += DPI; + + angle2 = observer_angle + PI; + if (angle2 > DPI) angle2 -= DPI; + + imin = (int)(observer_angle*(double)NX/DPI); + imax = (int)(angle2*(double)NX/DPI); + if (imin >= NX-1) imin = NX-2; + if (imax >= NX-1) imax = NX-2; + +// printf("Angle = %.5lg, angle2 = %.5lg, imin = %i, imax = %i\n", observer_angle, angle2, imin, imax); + + if (observer[2] > 0.0) + { + if (imin < imax) + { + for (i=imax; i>imin; i--) + for (j=0; j=0; i--) + for (j=0; j=imin; i--) + for (j=0; jimin; i--) + for (j=NY-3; j>=0; j--) + draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value); + + for (i=imax+1; i=0; j--) + draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value); + + for (j=NY-3; j>=0; j--) + draw_wave_sphere_ij(NX-1, 0, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value); + + for (i=0; i<=imin; i++) + for (j=NY-3; j>=0; j--) + draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value); + } + else + { + for (i=imax; i=0; j--) + draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value); + + for (i=imax-1; i>=0; i--) + for (j=NY-3; j>=0; j--) + draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value); + + for (j=NY-3; j>=0; j--) + draw_wave_sphere_ij(NX-1, 0, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value); + + for (i=NX-2; i>=imin; i--) + for (j=NY-3; j>=0; j--) + draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value); + } + + /* South pole */ + for (i=0; i0; j--) + draw_wave_sphere_ij(i, i+1, j-1, j, j, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value); + + for (j=2; j>0; j--) + draw_wave_sphere_ij(NX-1, 0, j-1, j, DPOLE, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value); + } +} + +void draw_wave_sphere(int movie, double phi[NX*NY], double psi[NX*NY], short int xy_in[NX*NY], t_wave wave[NX*NY], t_wave_sphere wsphere[NX*NY], int zplot, int cplot, int palette, int fade, double fade_value, int refresh) +{ + if (PLOT_2D) draw_wave_sphere_2D(movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade,fade_value, refresh); + else draw_wave_sphere_3D(movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade,fade_value, refresh); +} diff --git a/sub_wave.c b/sub_wave.c index 89252ed..c775d96 100644 --- a/sub_wave.c +++ b/sub_wave.c @@ -293,6 +293,15 @@ void write_text( double x, double y, char *st) return(alph); } + double iabs(int i) /* absolute value */ + { + int res; + + if (i<0) res = -i; + else res = i; + return(i); + } + double gaussian() /* returns standard normal random variable, using Box-Mueller algorithm */ @@ -446,6 +455,16 @@ void erase_area(double x, double y, double dx, double dy) glEnd(); } +void erase_area_ij(int imin, int jmin, int imax, int jmax) +{ + glColor3f(0.0, 0.0, 0.0); + glBegin(GL_QUADS); + glVertex2i(imin, jmin); + glVertex2i(imax, jmin); + glVertex2i(imax, jmax); + glVertex2i(imin, jmax); + glEnd(); +} void erase_area_rgb(double x, double y, double dx, double dy, double rgb[3]) { @@ -1004,6 +1023,63 @@ int init_circle_config_pattern(t_circle circles[NMAXCIRCLES], int circle_pattern } break; } + case (C_HEX_BOTTOM): + { + ncircles = NGRIDX*(NGRIDY+1); + dy = (- YMIN)/((double)NGRIDY); +// dx = dy*0.5*sqrt(3.0); + dx = (XMAX - XMIN)/((double)NGRIDX); + for (i = 0; i < NGRIDX; i++) + for (j = 0; j < NGRIDY+1; j++) + { + n = (NGRIDY+1)*i + j; + circles[n].xc = ((double)(i-NGRIDX/2) + 0.5)*dx; /* is +0.5 needed? */ + circles[n].yc = YMIN + ((double)j - 0.5)*dy; + if ((i+NGRIDX)%2 == 1) circles[n].yc += 0.5*dy; + circles[n].radius = MU; + /* activate only circles that intersect the domain */ + if ((circles[n].yc < YMAX + MU)&&(circles[n].yc > YMIN - MU)) circles[n].active = 1; + else circles[n].active = 0; + } + break; + } + case (C_HEX_BOTTOM2): + { + ncircles = NGRIDX*(NGRIDY+1); + dy = (- YMIN)/((double)NGRIDY); + dx = 2.0*LAMBDA/((double)NGRIDX); + for (i = 0; i < NGRIDX; i++) + for (j = 0; j < NGRIDY+1; j++) + { + n = (NGRIDY+1)*i + j; + circles[n].xc = ((double)(i-NGRIDX/2) + 0.5)*dx; /* is +0.5 needed? */ + circles[n].yc = YMIN + ((double)j - 0.5)*dy; + if ((i+NGRIDX)%2 == 1) circles[n].yc += 0.5*dy; + circles[n].radius = MU; + /* activate only circles that intersect the domain */ + if ((circles[n].yc < YMAX + MU)&&(circles[n].yc > YMIN - MU)) circles[n].active = 1; + else circles[n].active = 0; + } + break; + } + case (C_SQUARE_BOTTOM): + { + ncircles = NGRIDX*(NGRIDY+1); + dy = (- YMIN)/((double)NGRIDY); + dx = 2.0*LAMBDA/((double)NGRIDX); + for (i = 0; i < NGRIDX; i++) + for (j = 0; j < NGRIDY+1; j++) + { + n = (NGRIDY+1)*i + j; + circles[n].xc = ((double)(i-NGRIDX/2) + 0.5)*dx; + circles[n].yc = YMIN + ((double)j - 0.5)*dy; + circles[n].radius = MU; + /* activate only circles that intersect the domain */ + if ((circles[n].yc < YMAX + MU)&&(circles[n].yc > YMIN - MU)) circles[n].active = 1; + else circles[n].active = 0; + } + break; + } case (C_ONE): { circles[ncircles].xc = 0.0; @@ -5752,7 +5828,7 @@ void compute_laplacian(double phi[NX*NY], t_laplacian laplace[NX*NY], double del } } -double oscillating_bc(int time) +double oscillating_bc(int time, int j) { double t, phase, a, envelope, omega; @@ -5760,7 +5836,12 @@ double oscillating_bc(int time) { case (OSC_PERIODIC): { - return(AMPLITUDE*cos((double)time*OMEGA)*exp(-(double)time*DAMPING)); + if (OSCIL_LEFT_YSHIFT > 0.0) + t = (double)time*OMEGA - (double)j*OSCIL_LEFT_YSHIFT/(double)NY; + else + t = (double)time*OMEGA + (double)(NY-j)*OSCIL_LEFT_YSHIFT/(double)NY; + if (t < 0.0) return(0.0); + else return(AMPLITUDE*cos(t)*exp(-(double)t*DAMPING)); } case (OSC_SLOWING): { diff --git a/sub_wave_3d.c b/sub_wave_3d.c index 1da1380..b12a4b6 100644 --- a/sub_wave_3d.c +++ b/sub_wave_3d.c @@ -1850,6 +1850,7 @@ void draw_wave_3d_ij(int i, int j, int movie, double phi[NX*NY], double psi[NX*N { z_mid = compute_interpolated_colors_wave(i, j, xy_in, wave, palette, cplot, rgb_e, rgb_w, rgb_n, rgb_s, fade, fade_value, movie); + /* TODO: incorporate these in wave structure */ ij_to_xy(i, j, xy_sw); ij_to_xy(i+1, j, xy_se); ij_to_xy(i, j+1, xy_nw); @@ -1863,6 +1864,7 @@ void draw_wave_3d_ij(int i, int j, int movie, double phi[NX*NY], double psi[NX*N // z_mid += *wave[(i+1)*NY+j+1].potential*POT_FACT*0.25; // } + /* TODO: incorporate these in wave structure */ for (k=0; k<2; k++) xy_mid[k] = 0.25*(xy_sw[k] + xy_se[k] + xy_nw[k] + xy_ne[k]); if (AMPLITUDE_HIGH_RES == 1) @@ -2198,6 +2200,162 @@ void draw_color_scheme_palette_3d(double x1, double y1, double x2, double y2, in draw_rectangle_noscale(x1, y1, x2, y2); } +void draw_circular_color_scheme_palette_3d(double x1, double y1, double radiusx, double radiusy, int plot, double min, double max, int palette, int fade, double fade_value) +{ + int j, k, ij[2], jmin=0; + double x, y, dy, dy_e, dy_phase, rgb[3], value, lum, amp, dphi, pos[2], phi, xy[2]; + + printf("Drawing circular color scheme\n"); + + glBegin(GL_TRIANGLE_FAN); +// xy_to_pos(x1, y1, xy); +// glVertex2d(xy[0], xy[1]); + xy_to_ij(x1, y1, ij); + draw_vertex_ij(ij[0], ij[1]); + dy = (max - min)/360.0; + dy_e = max/360.0; + dy_phase = 1.0/360.0; + dphi = DPI/360.0; + + for (j = 0; j <= 360; j++) + { + switch (plot) { + case (P_3D_AMPLITUDE): + { + value = min + 1.0*dy*(double)(j - jmin); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + case (P_3D_ANGLE): + { + value = 1.0*dy*(double)(j - jmin); + color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 0, rgb); + break; + } + case (P_3D_AMP_ANGLE): + { + value = min + 1.0*dy*(double)(j - jmin); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + case (P_3D_ENERGY): + { + value = dy_e*(double)(j - jmin)*100.0/E_SCALE; + if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + break; + } + case (P_3D_LOG_ENERGY): + { + value = LOG_SCALE*dy_e*(double)(j - jmin)*100.0/E_SCALE; + color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + break; + } + case (P_3D_TOTAL_ENERGY): + { + value = dy_e*(double)(j - jmin)*100.0/E_SCALE; + if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + break; + } + case (P_3D_LOG_TOTAL_ENERGY): + { + value = LOG_SCALE*dy_e*(double)(j - jmin)*100.0/E_SCALE; + color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + break; + } + case (P_3D_MEAN_ENERGY): + { + value = dy_e*(double)(j - jmin)*100.0/E_SCALE; + if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + break; + } + case (P_3D_LOG_MEAN_ENERGY): + { + value = LOG_SCALE*dy_e*(double)(j - jmin)*100.0/E_SCALE; + color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + break; + } + case (P_3D_PHASE): + { + value = dy_phase*(double)(j - jmin); + color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb); + break; + } + case (P_3D_FLUX_INTENSITY): + { + value = dy_e*(double)(j - jmin)*100.0/E_SCALE; + if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + break; + } + case (P_3D_FLUX_DIRECTION): + { + value = dy_phase*(double)(j - jmin); + color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb); + break; + } + case (Z_EULER_VORTICITY): + { + value = min + 1.0*dy*(double)(j - jmin); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + case (Z_EULER_LOG_VORTICITY): + { + value = min + 1.0*dy*(double)(j - jmin); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + case (Z_EULER_VORTICITY_ASYM): + { + value = min + 1.0*dy*(double)(j - jmin); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + case (Z_EULER_LPRESSURE): + { + value = min + 1.0*dy*(double)(j - jmin); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + case (Z_EULERC_VORTICITY): + { + value = min + 1.0*dy*(double)(j - jmin); + printf("Palette value %.3lg\n", value); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + + } + if (fade) for (k=0; k<3; k++) rgb[k] *= fade_value; + glColor3f(rgb[0], rgb[1], rgb[2]); + + xy_to_ij(x1 + radiusx*cos(dphi*(double)j), y1 + radiusy*sin(dphi*(double)j), ij); + draw_vertex_ij(ij[0], ij[1]); + } + + xy_to_ij(x1 + radiusx*cos(dphi), y1 + radiusy*sin(dphi), ij); + draw_vertex_ij(ij[0], ij[1]); + glEnd (); + + if (fade) glColor3f(fade_value, fade_value, fade_value); + else glColor3f(1.0, 1.0, 1.0); + glLineWidth(BOUNDARY_WIDTH*3/2); + glEnable(GL_LINE_SMOOTH); + + dphi = DPI/(double)NSEG; + glBegin(GL_LINE_LOOP); + for (j = 0; j < NSEG; j++) + { + xy_to_ij(x1 + radiusx*cos(dphi*(double)j), y1 + radiusy*sin(dphi*(double)j), ij); + draw_vertex_ij(ij[0], ij[1]); + } + glEnd (); +} + + void print_speed_3d(double speed, int fade, double fade_value) { char message[100]; diff --git a/wave_3d.c b/wave_3d.c index a79ba10..8df21a5 100644 --- a/wave_3d.c +++ b/wave_3d.c @@ -156,6 +156,7 @@ #define KAPPA 0.0 /* "elasticity" term enforcing oscillations */ #define KAPPA_SIDES 5.0e-4 /* "elasticity" term on absorbing boundary */ #define KAPPA_TOPBOT 0.0 /* "elasticity" term on absorbing boundary */ +#define OSCIL_LEFT_YSHIFT 0.0 /* y-dependence of left oscillation (for non-horizontal waves) */ /* The Courant number is given by c*DT/DX, where DT is the time step and DX the lattice spacing */ /* The physical damping coefficient is given by GAMMA/(DT)^2 */ /* Increasing COURANT speeds up the simulation, but decreases accuracy */ @@ -289,6 +290,7 @@ // #define POT_MAZE 7 #define POTENTIAL 10 #define POT_FACT 20.0 +#define DRAW_WAVE_PROFILE 0 /* set to 1 to draw a profile of the wave */ /* end of constants only used by sub_wave and sub_maze */ @@ -364,8 +366,8 @@ void evolve_wave_half_new(double phi_in[NX*NY], double psi_in[NX*NY], double phi // if (OSCILLATE_LEFT) for (j=1; j= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, scale, time, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, value, scale, time, rgb); + break; + } + case (P_MIXED): + { + if (j > NY/2) value = phi[i][j]; + else value = compute_energy(phi, psi, xy_in, i, j); + color_scheme_palette(COLOR_SCHEME, palette, value, scale, time, rgb); + break; + } + case (P_MEAN_ENERGY): + { + energy = compute_energy(phi, psi, xy_in, i, j); + total_energy[i][j] += energy; + value = total_energy[i][j]/(double)(time+1); + if (COLOR_PALETTE >= COL_TURBO) + color_scheme_asym_palette(COLOR_SCHEME, palette, value, scale, time, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, value, scale, time, rgb); + break; + } + case (P_LOG_ENERGY): + { + energy = compute_energy(phi, psi, xy_in, i, j); +// energy = LOG_SHIFT + LOG_SCALE*log(energy); +// if (energy < 0.0) energy = 0.0; + value = LOG_SHIFT + LOG_SCALE*log(energy); + color_scheme_palette(COLOR_SCHEME, palette, value, scale, time, rgb); + break; + } + case (P_LOG_MEAN_ENERGY): + { + energy = compute_energy(phi, psi, xy_in, i, j); + if (energy == 0.0) energy = 1.0e-20; + total_energy[i][j] += energy; + value = LOG_SHIFT + LOG_SCALE*log(total_energy[i][j]/(double)(time+1)); + color_scheme_palette(COLOR_SCHEME, palette, value, scale, time, rgb); + break; + } + case (P_ENERGY_FLUX): + { + compute_energy_flux(phi, psi, xy_in, i, j, &gx, &gy, &arg, &mod); +// color_scheme_palette(C_ONEDIM_LINEAR, palette, arg/DPI, 1.0, 1, rgb); +// flux_factor = tanh(mod*E_SCALE); +// for (k=0; k<3; k++) rgb[k] *= flux_factor; + value = mod*FLUX_SCALE; + color_scheme_asym_palette(COLOR_SCHEME, palette, value, scale, time, rgb); + break; + } + case (P_TOTAL_ENERGY_FLUX): + { + compute_energy_flux(phi, psi, xy_in, i, j, &gx, &gy, &arg, &mod); + total_flux[2*j*NX + 2*i] *= 0.99; + total_flux[2*j*NX + 2*i + 1] *= 0.99; + total_flux[2*j*NX + 2*i] += gx; + total_flux[2*j*NX + 2*i + 1] += gy; +// mgx = total_flux[2*j*NX + 2*i]/(double)(time+1); +// mgy = total_flux[2*j*NX + 2*i + 1]/(double)(time+1); + mgx = total_flux[2*j*NX + 2*i]; + mgy = total_flux[2*j*NX + 2*i + 1]; +// mgx = total_flux[2*j*NX + 2*i]/log((double)(time+2)); +// mgy = total_flux[2*j*NX + 2*i + 1]/log((double)(time+2)); + mod = module2(mgx, mgy); + arg = argument(mgx, mgy); + if (arg < 0.0) arg += DPI; + value = arg/DPI; + color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb); + flux_factor = tanh(mod*E_SCALE); + for (k=0; k<3; k++) rgb[k] *= flux_factor; + break; + } + } + + return(value); +} + + +void draw_wave_profile(double *values, int size, int fade, double fade_value) +/* draw a profile of the wave, if option DRAW_WAVE_PROFILE is active */ +{ + int i, k; + double vmin, vmax, deltav, a, b, y; + static int imin, imax, jmin, jmax, jmid, d, first = 1; + static double deltaj, ymin; + + if (first) + { + imin = 100; + imax = NX - 250; + jmin = NY - 150; + jmax = NY - 50; + jmid = (jmin + jmax)/2; + jmid -= (jmid%size); + d = (jmax - jmin)/10 + 1; + deltaj = (double)(jmax - jmin - 2*d); + ymin = (double)(jmin + d); + first = 0; + } + + vmin = 1.0e10; + vmax = -vmin; + for (i=imin; i vmax) vmax = values[i*NY+jmin]; + if (values[i*NY+jmin] < vmin) vmin = values[i*NY+jmin]; + } + if ((vmin < 0.0)&&(vmax > 0.0)) + { + if (vmax > -vmin) vmin = -vmax; + else if (vmax < -vmin) vmax = -vmin; + } +// printf("vmin = %.3lg, vmax = %.3lg\n", vmin, vmax); + deltav = vmax-vmin; + if (deltav == 0.0) deltav = 0.01; + a = deltaj/deltav; + b = ymin - a*vmin; + + erase_area_ij(imin, jmin, imax, jmax); + + if (fade) glColor3f(fade_value, fade_value, fade_value); + else glColor3f(1.0, 1.0, 1.0); + glBegin(GL_LINE_STRIP); + for (i=imin; i= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, energy, scale, time, rgb); - else color_scheme_palette(COLOR_SCHEME, palette, energy, scale, time, rgb); - break; - } - case (P_MIXED): - { - if (j > NY/2) color_scheme_palette(COLOR_SCHEME, palette, phi[i][j], scale, time, rgb); - else color_scheme_palette(COLOR_SCHEME, palette, compute_energy(phi, psi, xy_in, i, j), scale, time, rgb); - break; - } - case (P_MEAN_ENERGY): - { - energy = compute_energy(phi, psi, xy_in, i, j); - total_energy[i][j] += energy; - if (COLOR_PALETTE >= COL_TURBO) - color_scheme_asym_palette(COLOR_SCHEME, palette, total_energy[i][j]/(double)(time+1), scale, time, rgb); - else color_scheme_palette(COLOR_SCHEME, palette, total_energy[i][j]/(double)(time+1), scale, time, rgb); - break; - } - case (P_LOG_ENERGY): - { - energy = compute_energy(phi, psi, xy_in, i, j); -// energy = LOG_SHIFT + LOG_SCALE*log(energy); -// if (energy < 0.0) energy = 0.0; - color_scheme_palette(COLOR_SCHEME, palette, LOG_SHIFT + LOG_SCALE*log(energy), scale, time, rgb); - break; - } - case (P_LOG_MEAN_ENERGY): - { - energy = compute_energy(phi, psi, xy_in, i, j); - if (energy == 0.0) energy = 1.0e-20; - total_energy[i][j] += energy; - color_scheme_palette(COLOR_SCHEME, palette, LOG_SHIFT + LOG_SCALE*log(total_energy[i][j]/(double)(time+1)), scale, time, rgb); - break; - } - case (P_ENERGY_FLUX): - { - compute_energy_flux(phi, psi, xy_in, i, j, &gx, &gy, &arg, &mod); -// color_scheme_palette(C_ONEDIM_LINEAR, palette, arg/DPI, 1.0, 1, rgb); -// flux_factor = tanh(mod*E_SCALE); -// for (k=0; k<3; k++) rgb[k] *= flux_factor; - color_scheme_asym_palette(COLOR_SCHEME, palette, mod*FLUX_SCALE, scale, time, rgb); - break; - } - case (P_TOTAL_ENERGY_FLUX): - { - compute_energy_flux(phi, psi, xy_in, i, j, &gx, &gy, &arg, &mod); - total_flux[2*j*NX + 2*i] *= 0.99; - total_flux[2*j*NX + 2*i + 1] *= 0.99; - total_flux[2*j*NX + 2*i] += gx; - total_flux[2*j*NX + 2*i + 1] += gy; -// mgx = total_flux[2*j*NX + 2*i]/(double)(time+1); -// mgy = total_flux[2*j*NX + 2*i + 1]/(double)(time+1); - mgx = total_flux[2*j*NX + 2*i]; - mgy = total_flux[2*j*NX + 2*i + 1]; -// mgx = total_flux[2*j*NX + 2*i]/log((double)(time+2)); -// mgy = total_flux[2*j*NX + 2*i + 1]/log((double)(time+2)); - mod = module2(mgx, mgy); - arg = argument(mgx, mgy); - if (arg < 0.0) arg += DPI; - color_scheme_palette(C_ONEDIM_LINEAR, palette, arg/DPI, 1.0, 1, rgb); - flux_factor = tanh(mod*E_SCALE); - for (k=0; k<3; k++) rgb[k] *= flux_factor; - break; - } - - } + value = wave_value(i, j, phi, psi, total_energy, total_flux, xy_in, scale, time, plot, palette, rgb); + + if (DRAW_WAVE_PROFILE) values[i*NY+j] = value; + if (fade) for (k=0; k<3; k++) rgb[k] *= fade_value; glColor3f(rgb[0], rgb[1], rgb[2]); @@ -960,6 +1043,12 @@ void draw_wave_highres_palette(int size, double *phi[NX], double *psi[NX], doubl } glEnd (); + + if (DRAW_WAVE_PROFILE) + { + draw_wave_profile(values, size, fade, fade_value); + free(values); + } } /* modified function for "flattened" wave tables */ diff --git a/wave_comparison.c b/wave_comparison.c index 5118942..5c62861 100644 --- a/wave_comparison.c +++ b/wave_comparison.c @@ -43,7 +43,7 @@ #include #include -#define MOVIE 1 /* set to 1 to generate movie */ +#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 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 */ @@ -263,6 +263,8 @@ // #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 */ /* end of constants only used by sub_wave and sub_maze */ diff --git a/wave_energy.c b/wave_energy.c index 23e4ce6..e5bf07a 100644 --- a/wave_energy.c +++ b/wave_energy.c @@ -219,8 +219,16 @@ #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 */ /* end of constants only used by sub_wave and sub_maze */ #include "global_pdes.c" /* constants and global variables */ diff --git a/wave_sphere.c b/wave_sphere.c new file mode 100644 index 0000000..3ab6305 --- /dev/null +++ b/wave_sphere.c @@ -0,0 +1,944 @@ +/*********************************************************************************/ +/* */ +/* Animation of wave equation on a sphere */ +/* */ +/* N. Berglund, july 2023 */ +/* */ +/* Feel free to reuse, but if doing so it would be nice to drop a */ +/* line to nils.berglund@univ-orleans.fr - Thanks! */ +/* */ +/* compile with */ +/* gcc -o wave_sphere wave_sphere.c */ +/* -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lXmu -lglut -O3 -fopenmp */ +/* */ +/* OMP acceleration may be more effective after executing */ +/* export OMP_NUM_THREADS=2 in the shell before running the program */ +/* */ +/* To make a video, set MOVIE to 1 and create subfolder tif_wave */ +/* It may be possible to increase parameter PAUSE */ +/* */ +/* create movie using */ +/* ffmpeg -i wave.%05d.tif -vcodec libx264 wave.mp4 */ +/* */ +/*********************************************************************************/ + +/*********************************************************************************/ +/* */ +/* NB: The algorithm used to simulate the wave equation is highly paralellizable */ +/* One could make it much faster by using a GPU */ +/* */ +/*********************************************************************************/ + +#include +#include +#include +#include +#include +#include +#include /* Sam Leffler's libtiff library. */ +#include +#include + +#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 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 WINHEIGHT 1150 /* window height */ +#define NX 2560 /* number of grid points on x axis */ +#define NY 1280 /* number of grid points on y axis */ + +#define DPOLE 20 /* safety distance to poles */ +#define SMOOTHPOLE 0.1 /* smoothing coefficient at poles */ + +#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 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 */ + +/* Choice of the billiard table */ + +#define B_DOMAIN 84 /* choice of domain shape, see list in global_pdes.c */ + +#define CIRCLE_PATTERN 33 /* 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 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_TOTAL_TURNS 1.0 /* total angle of rotation for IOR_PERIODIC_WELLS_ROTATING */ +#define MANDEL_IOR_SCALE -0.05 /* parameter controlling dependence of IoR on Mandelbrot escape speed */ + +#define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */ +#define NPOISSON 1000 /* number of points for Poisson C_RAND_POISSON arrangement */ +#define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */ + +#define LAMBDA 0.75 /* parameter controlling the dimensions of domain */ +#define MU 0.1 /* parameter controlling the dimensions of domain */ +#define NPOLY 6 /* number of sides of polygon */ +#define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */ +#define MDEPTH 7 /* depth of computation of Menger gasket */ +#define MRATIO 3 /* ratio defining Menger gasket */ +#define MANDELLEVEL 2000 /* iteration level for Mandelbrot set */ +#define MANDELLIMIT 20.0 /* limit value for approximation of Mandelbrot set */ +#define FOCI 1 /* set to 1 to draw focal points of ellipse */ +#define NGRIDX 30 /* number of grid point for grid of disks */ +#define NGRIDY 18 /* number of grid point for grid of disks */ + +#define X_SHOOTER -0.2 +#define Y_SHOOTER -0.6 +#define X_TARGET 0.4 +#define Y_TARGET 0.7 /* shooter and target positions in laser fight */ + +#define ISO_XSHIFT_LEFT -2.9 +#define ISO_XSHIFT_RIGHT 1.4 +#define ISO_YSHIFT_LEFT -0.15 +#define ISO_YSHIFT_RIGHT -0.15 +#define ISO_SCALE 0.5 /* coordinates for isospectral billiards */ + + +/* You can add more billiard tables by adapting the functions */ +/* xy_in_billiard and draw_billiard below */ + +/* 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_TOPBOT 0 /* set to 1 to enforce a planar wave on top and bottom boundary */ +#define OSCILLATION_SCHEDULE 3 /* oscillation schedule, see list in global_pdes.c */ + +#define OMEGA 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 */ +#define COURANT 0.05 /* Courant number */ +#define COURANTB 0.01 /* 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 */ +#define GAMMA_TOPBOT 1.0e-7 /* damping factor on boundary */ +#define KAPPA 0.0 /* "elasticity" term enforcing oscillations */ +#define KAPPA_SIDES 5.0e-4 /* "elasticity" term on absorbing boundary */ +#define KAPPA_TOPBOT 0.0 /* "elasticity" term on absorbing boundary */ +#define OSCIL_LEFT_YSHIFT 0.0 /* y-dependence of left oscillation (for non-horizontal waves) */ +/* The Courant number is given by c*DT/DX, where DT is the time step and DX the lattice spacing */ +/* The physical damping coefficient is given by GAMMA/(DT)^2 */ +/* Increasing COURANT speeds up the simulation, but decreases accuracy */ +/* 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 30 /* period of oscillating source */ +#define ALTERNATE_OSCILLATING_SOURCE 1 /* set to 1 to alternate sign of oscillating source */ + +#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 */ +#define N_WAVE_PACKETS 15 /* number of wave packets */ +#define WAVE_PACKET_RADIUS 20 /* radius of wave packets */ + +/* Boundary conditions, see list in global_pdes.c */ + +#define B_COND 2 + +#define PRECOMPUTE_BC 0 /* set to 1 to compute neighbours for Laplacian in advance */ + +/* Parameters for length and speed of simulation */ + +#define NSTEPS 3600 /* 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 */ +#define BOUNDARY_WIDTH 2 /* width of billiard boundary */ +#define PRINT_SPEED 0 /* set to 1 to print speed of moving source */ + +#define PAUSE 100 /* number of frames after which to pause */ +#define PSLEEP 3 /* sleep time during pause */ +#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 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 */ + +/* Plot type, see list in global_pdes.c */ + +#define ZPLOT 103 /* wave height */ +#define CPLOT 103 /* color scheme */ + +#define ZPLOT_B 108 +#define CPLOT_B 108 /* plot type for second movie */ + +#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_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 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 */ +#define COS_LIGHT_MAX 0.8 /* controls angle-dependence of 2D shading */ + +#define PLOT_SCALE_ENERGY 0.4 /* vertical scaling in energy plot */ +#define PLOT_SCALE_LOG_ENERGY 0.5 /* vertical scaling in log energy plot */ + +/* 3D representation */ + +#define REPRESENTATION_3D 1 /* choice of 3D representation */ +#define PLOT_2D 0 /* switch to 2D representation, equirectangular projection */ +#define PHISHIFT 0.0 /* shift of phi in 2D plot (in degrees) */ + +#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 VIEWPOINT_TRAJ 1 /* type of viewpoint trajectory */ + +/* Color schemes */ + +#define COLOR_PALETTE 11 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE_B 16 /* Color palette, see list in global_pdes.c */ + +#define BLACK 1 /* background */ +#define COLOR_OUT_R 1.0 /* color outside domain */ +#define COLOR_OUT_G 1.0 +#define COLOR_OUT_B 1.0 + +#define COLOR_SCHEME 3 /* choice of color scheme, see list in global_pdes.c */ + +#define SCALE 0 /* set to 1 to adjust color scheme to variance of field */ +#define 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 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 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 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 */ +#define COLORDRIFT 0.0 /* how much the color hue drifts during the whole simulation */ +#define LUMMEAN 0.5 /* amplitude of luminosity variation for scheme C_LUM */ +#define LUMAMP 0.3 /* amplitude of luminosity variation for scheme C_LUM */ +#define HUEMEAN 240.0 /* mean value of hue for color scheme C_HUE */ +#define HUEAMP -200.0 /* amplitude of variation of hue for color scheme C_HUE */ + +#define NXMAZE 8 /* width of maze */ +#define NYMAZE 32 /* height of maze */ +#define MAZE_MAX_NGBH 5 /* max number of neighbours of maze cell */ +#define RAND_SHIFT 5 /* seed of random number generator */ +#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 ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ +#define CIRC_COLORBAR 0 /* set to 1 to draw circular color scheme */ +#define CIRC_COLORBAR_B 0 /* set to 1 to draw circular color scheme */ + +#define DRAW_WAVE_PROFILE 0 /* set to 1 to draw a profile of the wave */ +#define SAVE_TIME_SERIES 0 /* set to 1 to save wave time series at a point */ + +#define ADD_POTENTIAL 0 /* set to 1 to add potential to z coordinate */ +#define POTENTIAL 10 +#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 */ + +/* Parameters controlling 3D projection */ + +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 */ +int reset_view = 0; /* switch to reset 3D view parameters (for option ROTATE_VIEW) */ + +#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 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 */ + +#include "global_pdes.c" /* constants and global variables */ +#include "global_3d.c" /* constants and global variables */ +#include "sub_maze.c" /* support for generating mazes */ +#include "sub_wave.c" /* common functions for wave_billiard, heat and schrodinger */ +#include "wave_common.c" /* common functions for wave_billiard, wave_comparison, etc */ + +#include "sub_wave_3d.c" /* graphical functions specific to wave_3d */ +#include "sub_sphere.c" /* graphical functions specific to wave_sphere */ + +FILE *time_series_left, *time_series_right, *image_file; + +double courant2, courantb2; /* Courant parameters squared */ + + +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]) +/* 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; + double delta, x, y, c, cc, gamma, sintheta, cottheta, invstheta, sum, avrg; + static long time = 0; + static short int first = 1; + static double dphi, dtheta, cphiphi, ctheta; + + if (first) + { + dphi = DPI/(double)NX; + dtheta = PI/(double)NY; + cphiphi = dphi*dphi/(dtheta*dtheta); + ctheta = dphi*dphi/(2.0*dtheta); + + printf("dphi = %.5lg, dtheta = %.5lg, cphiphi = %.5lg, ctheta = %.5lg\n", dphi, dtheta, cphiphi, ctheta); + + first = 0; + } + + time++; + + #pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta,x,y) + /* evolution in the bulk */ + for (j=DPOLE; j VMAX) phi_out[i*NY+j] = VMAX; + if (phi_out[i*NY+j] < -VMAX) phi_out[i*NY+j] = -VMAX; + } + } + } +} + + +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]) +/* 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, 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); +} + + +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.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); +} + +void viewpoint_schedule(int i) +/* change position of observer */ +{ + int j; + double angle, ca, sa, r1; + static double observer_initial[3], r, ratio; + static int first = 1; + + if (first) + { + for (j=0; j<3; j++) observer_initial[j] = observer[j]; + r1 = observer[0]*observer[0] + observer[1]*observer[1]; + r = sqrt(r1 + observer[2]*observer[2]); + ratio = r/sqrt(r1); + first = 0; + } + + angle = (ROTATE_ANGLE*DPI/360.0)*(double)i/(double)NSTEPS; + ca = cos(angle); + sa = sin(angle); + switch (VIEWPOINT_TRAJ) + { + case (VP_HORIZONTAL): + { + observer[0] = ca*observer_initial[0] - sa*observer_initial[1]; + observer[1] = sa*observer_initial[0] + ca*observer_initial[1]; + break; + } + case (VP_ORBIT): + { + observer[0] = ca*observer_initial[0] - sa*observer_initial[1]*ratio; + observer[1] = ca*observer_initial[1] + sa*observer_initial[0]*ratio; + observer[2] = ca*observer_initial[2]; + break; + } + } + + printf("Angle %.3lg, Observer position (%.3lg, %.3lg, %.3lg)\n", angle, observer[0], observer[1], observer[2]); +} + + +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, theta, amp; + 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; + static int counter = 0, first_source = 1; + long int wave_value; + t_wave *wave; + t_laplacian *laplace, *laplace1, *laplace2; + t_wave_source wave_source[25]; + t_wave_sphere *wsphere; + + if (SAVE_TIME_SERIES) + { + time_series_left = fopen("wave_left.dat", "w"); + time_series_right = fopen("wave_right.dat", "w"); + } + + /* Since NX and NY are big, it seemed wiser to use some memory allocation here */ + xy_in = (short int *)malloc(NX*NY*sizeof(short int)); + phi = (double *)malloc(NX*NY*sizeof(double)); + psi = (double *)malloc(NX*NY*sizeof(double)); + tmp = (double *)malloc(NX*NY*sizeof(double)); +// total_energy = (double *)malloc(NX*NY*sizeof(double)); + color_scale = (double *)malloc(NX*NY*sizeof(double)); + tc = (double *)malloc(NX*NY*sizeof(double)); + tcc = (double *)malloc(NX*NY*sizeof(double)); + tgamma = (double *)malloc(NX*NY*sizeof(double)); + + wave = (t_wave *)malloc(NX*NY*sizeof(t_wave)); + wsphere = (t_wave_sphere *)malloc(NX*NY*sizeof(t_wave_sphere)); + + laplace = (t_laplacian *)malloc(NX*NY*sizeof(t_laplacian)); + laplace1 = (t_laplacian *)malloc(NX*NY*sizeof(t_laplacian)); + laplace2 = (t_laplacian *)malloc(NX*NY*sizeof(t_laplacian)); + + /* initialise positions and radii of circles */ + if (COMPARISON) + { + if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN == D_CIRCLES_IN_RECT)) + 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)) + 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); + else if (B_DOMAIN == D_POLYGONS) ncircles = init_polygon_config(polygons); + } + + printf("Polygons initialized\n"); + + /* initialise polyline for von Koch and similar domains */ + npolyline = init_polyline(MDEPTH, polyline); + for (i=0; i= INITIAL_TIME) save_frame(); +// if (i >= INITIAL_TIME) save_frame_counter(i); +// if (i >= INITIAL_TIME) save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter); + else printf("Initial phase time %i of %i\n", i, INITIAL_TIME); + + if ((i >= INITIAL_TIME)&&(DOUBLE_MOVIE)) + { + draw_wave_sphere(1, phi, psi, xy_in, wave, wsphere, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 0, 1.0, 1); + if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, CIRC_COLORBAR_B, 0, 1.0); + if (PRINT_SPEED) print_speed_3d(speed, 0, 1.0); + glutSwapBuffers(); + save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter); +// save_frame_counter(i); + counter++; + } + else if (NO_EXTRA_BUFFER_SWAP) glutSwapBuffers(); + + /* it seems that saving too many files too fast can cause trouble with the file system */ + /* so this is to make a pause from time to time - parameter PAUSE may need adjusting */ + if (i % PAUSE == PAUSE - 1) + { + printf("Making a short pause\n"); + sleep(PSLEEP); + s = system("mv wave*.tif tif_wave/"); + } + } + else printf("Computing frame %i\n", i); + + } + + if (MOVIE) + { + if (DOUBLE_MOVIE) + { + draw_wave_sphere(0, phi, psi, xy_in, wave, wsphere, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1); + if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, 0, 1.0); + if (PRINT_SPEED) print_speed_3d(speed, 0, 1.0); + glutSwapBuffers(); + + if (!FADE) for (i=0; i