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

Add files via upload

Browse Source
This commit is contained in:
Nils Berglund 2023-09-02 11:36:15 +02:00 committed by GitHub
parent ca83ff16e0
commit c7de52bfa4
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
20 changed files with 3227 additions and 365 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
Earth_Map_Blue_Marble_2002_large.ppm.gz Normal file
View File

Binary file not shown.

29
global_3d.c
View File

@ -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 */

34
global_ljones.c
View File

@ -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;

17
global_pdes.c
View File

@ -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 */

2
heat.c
View File

@ -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"

277
lennardjones.c
View File

@ -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<N_PRESSURES; j++) pressure[j] = 0.0;
@ -1245,16 +1297,22 @@ void animation()
{
if (MOVE_OBSTACLE)
{
xmincontainer = obstacle_schedule_smooth(i, n);
xshift = xmincontainer;
params.xmincontainer = obstacle_schedule_smooth(i, n);
xshift = params.xmincontainer;
}
if ((ROTATE_BOUNDARY)&&(SMOOTH_ROTATION))
rotate_segments(segment, rotation_schedule_smooth(i,n));
if (ROTATE_BOUNDARY)
{
params.omega = angular_speed;
params.angle = rotation_schedule_smooth(i,n);
}
if ((ROTATE_BOUNDARY)&&(SMOOTH_ROTATION)) rotate_segments(segment, rotation_schedule_smooth(i,n));
if (INCREASE_GRAVITY) gravity = gravity_schedule(i,n);
if (INCREASE_GRAVITY) params.gravity = gravity_schedule(i,n);
if ((BOUNDARY_COND == BC_RECTANGLE_WALL)&&(i < INITIAL_TIME + WALL_TIME)) wall = 1;
else wall = 0;
if ((MOVE_BOUNDARY)||(MOVE_SEGMENT_GROUPS)) for (j=0; j<nsegments; j++)
if ((MOVE_BOUNDARY)||(MOVE_SEGMENT_GROUPS)||(PRINT_SEGMENTS_FORCE)) for (j=0; j<nsegments; j++)
{
segment[j].fx = 0.0;
segment[j].fy = 0.0;
@ -1272,10 +1330,10 @@ void animation()
particle[j].torque = 0.0;
/* compute force from other particles */
compute_particle_force(j, krepel, particle, hashgrid);
compute_particle_force(j, params.krepel, particle, hashgrid);
/* take care of boundary conditions */
fboundary += compute_boundary_force(j, particle, obstacle, segment, xmincontainer, xmaxcontainer, &pleft, &pright, pressure, wall);
params.fboundary += compute_boundary_force(j, particle, obstacle, segment, params.xmincontainer, params.xmaxcontainer, &pleft, &pright, pressure, wall);
/* align velocities in case of Vicsek models */
// if (VICSEK_INT)
@ -1299,7 +1357,7 @@ void animation()
}
/* add gravity */
if (INCREASE_GRAVITY) particle[j].fy -= gravity/particle[j].mass_inv;
if (INCREASE_GRAVITY) particle[j].fy -= params.gravity/particle[j].mass_inv;
else
{
particle[j].fy -= GRAVITY/particle[j].mass_inv;
@ -1318,14 +1376,14 @@ void animation()
}
/* timestep of thermostat algorithm */
totalenergy = evolve_particles(particle, hashgrid, qx, qy, qangle, px, py, pangle, beta, &nactive, &nsuccess, &nmove, i < INITIAL_TIME);
totalenergy = evolve_particles(particle, hashgrid, qx, qy, qangle, px, py, pangle, params.beta, &params.nactive, &nsuccess, &nmove, i < INITIAL_TIME);
/* evolution of lid coordinate */
if (BOUNDARY_COND == BC_RECTANGLE_LID) evolve_lid(fboundary);
if (BOUNDARY_COND == BC_RECTANGLE_LID) evolve_lid(params.fboundary);
if (BOUNDARY_COND == BC_RECTANGLE_WALL)
{
if (i < INITIAL_TIME + WALL_TIME) evolve_wall(fboundary);
if (i < INITIAL_TIME + WALL_TIME) evolve_wall(params.fboundary);
else xwall = 0.0;
}
if ((MOVE_BOUNDARY)&&(i > 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<ncircles; j++) if (particle[j].active)
{
nactive++;
params.nactive++;
qx[j] = particle[j].xc;
qy[j] = particle[j].yc;
px[j] = particle[j].vx;
@ -1470,8 +1528,8 @@ void animation()
}
if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length);
draw_particles(particle, PLOT, beta, collisions, ncollisions);
draw_container(xmincontainer, xmaxcontainer, obstacle, segment, wall);
draw_particles(particle, PLOT, params.beta, collisions, ncollisions);
draw_container(params.xmincontainer, params.xmaxcontainer, obstacle, segment, wall);
/* add a particle */
if ((ADD_PARTICLES)&&(i > 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<ncircles; j++) particle[j].radius *= radius_ratio;
printf("Particle 0 radius %.5lg\n", particle[0].radius);
params.radius *= radius_ratio;
}
/* compute force on segments */
if (PRINT_SEGMENTS_FORCE) compute_segments_force(&params, segment);
update_hashgrid(particle, hashgrid, 1);
if (PRINT_PARAMETERS)
print_parameters(beta, mean_energy, krepel, xmaxcontainer - xmincontainer,
fboundary/(double)(ncircles*NVID), PRINT_LEFT, pressure, gravity);
if (PRINT_PARAMETERS) print_parameters(params, PRINT_LEFT, pressure, 1);
if ((BOUNDARY_COND == BC_EHRENFEST)||(BOUNDARY_COND == BC_RECTANGLE_WALL))
print_ehrenfest_parameters(particle, pleft, pright);
else if (PRINT_PARTICLE_NUMBER)
@ -1504,7 +1573,6 @@ void animation()
if (PLOT_SPEEDS) draw_speed_plot(group_speeds, i);
if (PLOT_TRAJECTORIES) draw_trajectory_plot(group_speeds, i);
if (PRINT_OMEGA) print_omega(angular_speed);
else if (PRINT_PARTICLE_SPEEDS) print_particles_speeds(particle);
else if (PRINT_SEGMENTS_SPEEDS)
{
@ -1546,11 +1614,9 @@ void animation()
if ((i >= 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<ncircles; j++) if (particle[j].active) nactive++;
printf("%i active particles\n", nactive);
params.nactive = 0;
for (j=0; j<ncircles; j++) if (particle[j].active) params.nactive++;
printf("%i active particles\n", params.nactive);
printf("1\n");
// printf("1\n");
free(particle);
printf("2\n");
// printf("2\n");
if (ADD_FIXED_OBSTACLES) free(obstacle);
printf("3\n");
// printf("3\n");
if (ADD_FIXED_SEGMENTS)
{
free(segment);
free(segment_group);
}
printf("4\n");
// printf("4\n");
if (MOVE_SEGMENT_GROUPS) free(group_speeds);
printf("5\n");
// printf("5\n");
if (TRACER_PARTICLE) free(trajectory);
printf("6\n");
// printf("6\n");
if (PLOT_SPEEDS) free(obstacle_speeds);
printf("7\n");
// printf("7\n");
if (PLOT_PARTICLE_NUMBER) free(particle_numbers);
printf("8\n");
// printf("8\n");
free(hashgrid);
printf("9\n");
// printf("9\n");
free(qx);
printf("10\n");
// printf("10\n");
free(qy);
printf("11\n");
// printf("11\n");
free(px);
printf("12\n");
// printf("12\n");
free(py);
printf("13\n");
// printf("13\n");
free(qangle);
printf("14\n");
// printf("14\n");
free(pangle);
printf("15\n");
// printf("15\n");
free(pressure);
printf("16\n");
// printf("16\n");
if (REACTION_DIFFUSION) free(collisions);
if (SAVE_TIME_SERIES)

2
mangrove.c
View File

@ -259,6 +259,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 */
#include "global_pdes.c"

73
rde.c
View File

@ -49,13 +49,13 @@
#define WINWIDTH 1920 /* window width */
#define WINHEIGHT 1150 /* window height */
// #define NY 575 /* number of grid points on y axis */
// #define NX 960 /* number of grid points on x axis */
// #define NY 575 /* number of grid points on y axis */
#define NX 960 /* number of grid points on x axis */
#define NY 575 /* number of grid points on y axis */
// #define WINWIDTH 1280 /* window width */
// #define WINHEIGHT 720 /* window height */
#define NX 640 /* number of grid points on x axis */
#define NY 360 /* number of grid points on y axis */
// #define NX 640 /* number of grid points on x axis */
// #define NY 360 /* number of grid points on y axis */
// #define NX 320 /* number of grid points on x axis */
// #define NY 180 /* number of grid points on y axis */
@ -76,16 +76,16 @@
#define ADD_POTENTIAL 0 /* set to 1 to add a potential (for Schrodinger equation) */
#define ADD_MAGNETIC_FIELD 0 /* set to 1 to add a magnetic field (for Schrodinger equation) - then set POTENTIAL 1 */
#define ADD_FORCE_FIELD 0 /* set to 1 to add a foce field (for compressible Euler equation) */
#define ADD_FORCE_FIELD 1 /* set to 1 to add a foce field (for compressible Euler equation) */
#define POTENTIAL 7 /* type of potential or vector potential, see list in global_3d.c */
#define FORCE_FIELD 5 /* type of force field, see list in global_3d.c */
#define ADD_CORIOLIS_FORCE 0 /* set to 1 to add Coriolis force (quasigeostrophic Euler equations) */
#define VARIABLE_DEPTH 1 /* set to 1 for variable depth in shallow water equation */
#define VARIABLE_DEPTH 0 /* set to 1 for variable depth in shallow water equation */
#define SWATER_DEPTH 4 /* variable depth in shallow water equation */
#define ANTISYMMETRIZE_WAVE_FCT 0 /* set tot 1 to make wave function antisymmetric */
#define ADAPT_STATE_TO_BC 0 /* to smoothly adapt initial state to obstacles */
#define OBSTACLE_GEOMETRY 7 /* geometry of obstacles, as in B_DOMAIN */
#define ADAPT_STATE_TO_BC 1 /* to smoothly adapt initial state to obstacles */
#define OBSTACLE_GEOMETRY 1 /* geometry of obstacles, as in B_DOMAIN */
// #define BC_STIFFNESS 100.0 /* controls region of boundary condition control */
#define BC_STIFFNESS 50.0 /* controls region of boundary condition control */
#define CHECK_INTEGRAL 1 /* set to 1 to check integral of first field */
@ -102,7 +102,7 @@
#define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */
#define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */
#define LAMBDA 0.5 /* parameter controlling the dimensions of domain */
#define LAMBDA 1.2 /* parameter controlling the dimensions of domain */
#define MU 0.06 /* parameter controlling the dimensions of domain */
#define NPOLY 5 /* number of sides of polygon */
#define APOLY 2.0 /* angle by which to turn polygon, in units of Pi/2 */
@ -129,7 +129,7 @@
/* You can add more billiard tables by adapting the functions */
/* xy_in_billiard and draw_billiard in sub_wave.c */
/* Physical patameters of wave equation */
/* Physical parameters of wave equation */
#define DT 0.00000025
@ -164,6 +164,10 @@
#define SMOOTHEN_PERIOD 10 /* period between smoothenings */
#define SMOOTH_FACTOR 0.15 /* factor by which to smoothen */
#define ADD_OSCILLATING_SOURCE 1 /* set to 1 to add an oscillating wave source */
#define OSCILLATING_SOURCE_PERIOD 1 /* period of oscillating source */
#define OSCILLATING_SOURCE_OMEGA 0.2 /* frequency of oscillating source */
#define ADD_TRACERS 1 /* set to 1 to add tracer particles (for Euler equations) */
#define N_TRACERS 1000 /* number of tracer particles */
#define TRACERS_STEP 0.005 /* step size in tracer evolution */
@ -201,7 +205,7 @@
#define LAMINAR_FLOW_YPERIOD 1.0 /* period of laminar flow in y direction */
#define PRESSURE_GRADIENT 0.2 /* amplitude of pressure gradient for Euler equation */
#define SWATER_MIN_HEIGHT 2.0 /* min height of initial condition for shallow water equation */
#define SWATER_MIN_HEIGHT 0.5 /* min height of initial condition for shallow water equation */
// #define DEPTH_FACTOR 0.0125 /* proportion of min height in variable depth */
// #define DEPTH_FACTOR 0.001 /* proportion of min height in variable depth */
// #define DEPTH_FACTOR 0.0012 /* proportion of min height in variable depth */
@ -216,20 +220,20 @@
/* Boundary conditions, see list in global_pdes.c */
#define B_COND 1
#define B_COND 0
#define B_COND_LEFT 0
#define B_COND_RIGHT 0
#define B_COND_TOP 1
#define B_COND_BOTTOM 1
#define B_COND_TOP 0
#define B_COND_BOTTOM 0
/* Parameters for length and speed of simulation */
#define NSTEPS 2100 /* number of frames of movie */
#define NSTEPS 1300 /* number of frames of movie */
// #define NSTEPS 500 /* number of frames of movie */
// #define NVID 100 /* number of iterations between images displayed on screen */
#define NVID 30 /* number of iterations between images displayed on screen */
#define NVID 45 /* number of iterations between images displayed on screen */
#define ACCELERATION_FACTOR 1.0 /* factor by which to increase NVID in course of simulation */
#define DT_ACCELERATION_FACTOR 1.0 /* factor by which to increase time step in course of simulation */
#define MAX_DT 0.024 /* maximal value of integration step */
@ -249,7 +253,7 @@
#define PLOT_3D 1 /* controls whether plot is 2D or 3D */
#define ROTATE_VIEW 0 /* set to 1 to rotate position of observer */
#define ROTATE_VIEW 1 /* set to 1 to rotate position of observer */
#define ROTATE_ANGLE 360.0 /* total angle of rotation during simulation */
#define DRAW_PERIODICISED 0 /* set to 1 to repeat wave periodically in x and y directions */
@ -273,8 +277,8 @@
#define DRAW_OUTSIDE_GRAY 0 /* experimental - draw outside of billiard in gray */
#define ADD_POTENTIAL_TO_Z 0 /* set to 1 to add the external potential to z-coordinate of plot */
#define ADD_POT_CONSTANT 0.35 /* constant added to wave height */
#define ADD_HEIGHT_CONSTANT -2.0 /* constant added to wave height */
#define DRAW_DEPTH 1 /* set to 1 to draw water depth */
#define ADD_HEIGHT_CONSTANT -0.5 /* constant added to wave height */
#define DRAW_DEPTH 0 /* set to 1 to draw water depth */
#define DEPTH_SCALE 0.75 /* vertical scaling of depth plot */
#define DEPTH_SHIFT -0.015 /* vertical shift of depth plot */
@ -306,9 +310,9 @@
/* Color schemes, see list in global_pdes.c */
#define COLOR_PALETTE 11 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 10 /* Color palette, see list in global_pdes.c */
// #define COLOR_PALETTE_B 0 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE 14 /* Color palette, see list in global_pdes.c */
// #define COLOR_PALETTE_B 17 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 0 /* Color palette, see list in global_pdes.c */
#define BLACK 1 /* black background */
@ -318,11 +322,11 @@
#define SCALE 0 /* set to 1 to adjust color scheme to variance of field */
#define SLOPE 1.0 /* sensitivity of color on wave amplitude */
#define VSCALE_AMPLITUDE 10.0 /* additional scaling factor for color scheme P_3D_AMPLITUDE */
#define VSCALE_AMPLITUDE 15.0 /* additional scaling factor for color scheme P_3D_AMPLITUDE */
#define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */
#define CURL_SCALE 0.000015 /* scaling factor for curl representation */
#define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */
#define SLOPE_SCHROD_LUM 200.0 /* sensitivity of luminosity on module, for color scheme Z_ARGUMENT */
#define SLOPE_SCHROD_LUM 400.0 /* sensitivity of luminosity on module, for color scheme Z_ARGUMENT */
#define MIN_SCHROD_LUM 0.2 /* minimal luminosity in color scheme Z_ARGUMENT*/
#define VSCALE_PRESSURE 2.0 /* additional scaling factor for color scheme Z_EULER_PRESSURE */
#define PRESSURE_SHIFT 10.0 /* shift for color scheme Z_EULER_PRESSURE */
@ -347,7 +351,7 @@
#define VSCALE_VORTICITY 20.0 /* additional scaling factor for color scheme Z_EULERC_VORTICITY */
#define VORTICITY_SHIFT 0.0 /* vertical shift of vorticity */
#define ZSCALE_SPEED 0.3 /* additional scaling factor for z-coord Z_EULER_SPEED and Z_SWATER_SPEED */
#define VSCALE_SWATER 300.0 /* additional scaling factor for color scheme Z_EULER_DENSITY */
#define VSCALE_SWATER 250.0 /* additional scaling factor for color scheme Z_EULER_DENSITY */
#define NXMAZE 7 /* width of maze */
#define NYMAZE 7 /* height of maze */
@ -390,6 +394,8 @@
#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 */
#define OSCIL_LEFT_YSHIFT 25.0 /* y-dependence of left oscillation (for non-horizontal waves) */
#define DRAW_WAVE_PROFILE 1 /* set to 1 to draw a profile of the wave */
/* end of constants added only for compatibility with wave_common.c */
@ -400,8 +406,8 @@ double light[3] = {0.816496581, 0.40824829, 0.40824829}; /* vector of "ligh
double observer[3] = {8.0, 8.0, 7.0}; /* location of observer for REP_PROJ_3D representation */
int reset_view = 0; /* switch to reset 3D view parameters (for option ROTATE_VIEW) */
#define Z_SCALING_FACTOR 35.0 /* overall scaling factor of z axis for REP_PROJ_3D representation */
#define XY_SCALING_FACTOR 1.7 /* overall scaling factor for on-screen (x,y) coordinates after projection */
#define Z_SCALING_FACTOR 150.0 /* overall scaling factor of z axis for REP_PROJ_3D representation */
#define XY_SCALING_FACTOR 2.5 /* 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.1 /* overall y shift for REP_PROJ_3D representation */
@ -1565,7 +1571,7 @@ void viewpoint_schedule(int i)
void animation()
{
double time = 0.0, scale, dx, var, jangle, cosj, sinj, sqrintstep,
intstep0, viscosity_printed, fade_value, noise = NOISE_INTENSITY, x, y, sign;
intstep0, viscosity_printed, fade_value, noise = NOISE_INTENSITY, x, y, sign, phase;
double *phi[NFIELDS], *phi_tmp[NFIELDS], *potential_field, *vector_potential_field, *tracers, *gfield, *bc_field, *bc_field2;
short int *xy_in;
int i, j, k, s, nvid, field;
@ -1672,10 +1678,12 @@ void animation()
// init_gaussian_wave(-1.0, 0.0, 0.005, 0.25, SWATER_MIN_HEIGHT, phi, xy_in);
init_linear_wave(-1.0, 0.01, 5.0e-8, 0.25, SWATER_MIN_HEIGHT, phi, xy_in);
// init_linear_wave(-1.0, 0.01, 5.0e-8, 0.25, SWATER_MIN_HEIGHT, phi, xy_in);
// init_linear_blob(0.0, -0.75, 0.025, 0.0, 0.001, 0.25, 0.5, SWATER_MIN_HEIGHT, phi, xy_in);
init_elliptical_vibration(0.0, 0.0, 0.0075, LAMBDA, 1.0, 0.0003, 0.1, 1.0, SWATER_MIN_HEIGHT, phi, xy_in);
// add_gaussian_wave(-1.6, -0.5, 0.015, 0.25, SWATER_MIN_HEIGHT, phi, xy_in);
if (ADAPT_STATE_TO_BC) adapt_state_to_bc(phi, bc_field, xy_in);
@ -1759,6 +1767,13 @@ void animation()
if (ADAPT_STATE_TO_BC) adapt_state_to_bc(phi, bc_field, xy_in);
if ((ADD_OSCILLATING_SOURCE)&&(i%OSCILLATING_SOURCE_PERIOD == 0))
{
phase = (double)i*OSCILLATING_SOURCE_OMEGA;
printf("Adding vibration with phase %.3lg\n", phase);
add_elliptical_vibration(0.0, 0.0, 0.00007*cos(phase), LAMBDA, 1.0, 0.00004*cos(phase), 0.1, 1.0, SWATER_MIN_HEIGHT, phi, xy_in);
}
if (ADD_TRACERS)
{
printf("Evolving tracer particles\n");

1
schrodinger.c
View File

@ -191,6 +191,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 */

641
sub_lj.c
View File

File diff suppressed because it is too large Load Diff

86
sub_rde.c
View File

@ -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<NX; i++)
for (j=0; j<NY; j++)
{
ij_to_xy(i, j, xy);
xy_in[i*NY+j] = xy_in_billiard(xy[0],xy[1]);
if (xy_in[i*NY+j])
{
r = module2(xy[0]-x, xy[1]-y);
angle = argument(xy[0]-x, xy[1]-y);
height = exp(-(r - radius)*(r - radius)/var)*cos(omega*angle);
phi[0][i*NY+j] = min + amp*height;
phi[1][i*NY+j] = v*cos(angle)*height;
phi[2][i*NY+j] = v*sin(angle)*height;
}
else
{
phi[0][i*NY+j] = 0.0;
phi[1][i*NY+j] = 0.0;
phi[2][i*NY+j] = 0.0;
}
}
}
double init_elliptical_vibration(double x, double y, double v, double radiusx, double radiusy, 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<NX; i++)
for (j=0; j<NY; j++)
{
ij_to_xy(i, j, xy);
xy_in[i*NY+j] = xy_in_billiard(xy[0],xy[1]);
if (xy_in[i*NY+j])
{
r = module2((xy[0]-x)/radiusx, (xy[1]-y)/radiusy);
angle = argument((xy[0]-x)/radiusx, (xy[1]-y)/radiusy);
height = exp(-(r - 1.0)*(r - 1.0)/var)*cos(omega*angle);
phi[0][i*NY+j] = min + amp*height;
phi[1][i*NY+j] = v*cos(angle)*height*radiusx;
phi[2][i*NY+j] = v*sin(angle)*height*radiusy;
}
else
{
phi[0][i*NY+j] = 0.0;
phi[1][i*NY+j] = 0.0;
phi[2][i*NY+j] = 0.0;
}
}
}
double add_elliptical_vibration(double x, double y, double v, double radiusx, double radiusy, 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<NX; i++)
for (j=0; j<NY; j++)
{
ij_to_xy(i, j, xy);
xy_in[i*NY+j] = xy_in_billiard(xy[0],xy[1]);
if (xy_in[i*NY+j])
{
r = module2((xy[0]-x)/radiusx, (xy[1]-y)/radiusy);
if (r < 1.0 + var)
{
angle = argument((xy[0]-x)/radiusx, (xy[1]-y)/radiusy);
height = exp(-(r - 1.0)*(r - 1.0)/var)*cos(omega*angle);
phi[0][i*NY+j] += amp*height;
phi[1][i*NY+j] += v*cos(angle)*height*radiusx;
phi[2][i*NY+j] += v*sin(angle)*height*radiusy;
}
}
}
}
double add_gaussian_wave(double x, double y, double amp, double radius, double min, double *phi[NFIELDS], short int xy_in[NX*NY])
/* initialise gaussian wave height for shallow water equation */
{

895
sub_sphere.c Normal file
View File

@ -0,0 +1,895 @@
void init_sphere_2D() /* initialisation of window */
{
glLineWidth(3);
glClearColor(0.0, 0.0, 0.0, 1.0);
glClear(GL_COLOR_BUFFER_BIT);
glOrtho(0.0, NX, DPOLE, NY-DPOLE, -1.0, 1.0);
}
void draw_vertex_sphere(double xyz[3])
{
double xy_screen[2];
xyz_to_xy(xyz[0], xyz[1], xyz[2], xy_screen);
glVertex2d(xy_screen[0], xy_screen[1]);
}
void init_wave_sphere(t_wave_sphere wsphere[NX*NY])
/* initialize sphere data */
{
int i, j;
double dphi, dtheta, theta0, xy[2], phishift;
printf("Initializing wsphere\n");
dphi = DPI/(double)NX;
// dtheta = PI/(double)NY;
dtheta = PI/(double)(NY-2*(DPOLE));
theta0 = (double)(DPOLE)*dtheta;
phishift = PHISHIFT*(XMAX-XMIN)/360.0;
#pragma omp parallel for private(i,j)
for (i=0; i<NX; i++)
{
for (j=DPOLE; j<NY-DPOLE; j++)
wsphere[i*NY+j].theta = (double)j*dtheta - theta0;
for (j=0; j<DPOLE; j++) wsphere[i*NY+j].theta = 0.0;
for (j=NY-DPOLE; j<NY; j++) wsphere[i*NY+j].theta = PI;
for (j=0; j<NY; j++)
{
wsphere[i*NY+j].phi = (double)i*dphi;
// wsphere[i*NY+j].theta = (double)j*dtheta;
wsphere[i*NY+j].cphi = cos(wsphere[i*NY+j].phi);
wsphere[i*NY+j].sphi = sin(wsphere[i*NY+j].phi);
wsphere[i*NY+j].ctheta = cos(wsphere[i*NY+j].theta);
wsphere[i*NY+j].stheta = sin(wsphere[i*NY+j].theta);
wsphere[i*NY+j].x = wsphere[i*NY+j].cphi*wsphere[i*NY+j].stheta;
wsphere[i*NY+j].y = wsphere[i*NY+j].sphi*wsphere[i*NY+j].stheta;
wsphere[i*NY+j].z = -wsphere[i*NY+j].ctheta;
ij_to_xy(NX-1-i,j,xy);
// xy[0] = XMIN + ((double)(NX-i-1))*(XMAX-XMIN)/((double)NX);
// xy[1] = YMIN + ((double)(j-DPOLE))*(YMAX-YMIN)/((double)(NY-2*DPOLE));
xy[0] += phishift;
if (xy[0] > 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<NY-1; j++) wsphere[i*NY+j].cottheta = wsphere[i*NY+j].ctheta/wsphere[i*NY+j].stheta;
wsphere[i*NY].cottheta = wsphere[i*NY+1].cottheta;
wsphere[i*NY + NY-1].cottheta = wsphere[i*NY+1].cottheta;
}
}
void init_earth_map(t_wave_sphere wsphere[NX*NY])
/* init file from earth map */
{
int i, j, ii, jj, k, nx, ny, maxrgb, nmaxpixels = 4915200, scan, rgbval, diff, sshift, nshift;
int *rgb_values;
double cratio, rx, ry, cy;
FILE *image_file;
printf("Reading Earth map\n");
rgb_values = (int *)malloc(3*nmaxpixels*sizeof(int));
image_file = fopen("Earth_Map_Blue_Marble_2002_large.ppm", "r");
scan = fscanf(image_file,"%i %i\n", &nx, &ny);
scan = fscanf(image_file,"%i\n", &maxrgb);
if (nx*ny > 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<ny; j++)
for (i=0; i<nx; i++)
for (k=0; k<3; k++)
{
scan = fscanf(image_file,"%i\n", &rgbval);
rgb_values[3*(j*nx+i)+k] = rgbval;
}
cratio = 1.0/(double)maxrgb;
rx = (double)nx/(double)NX;
ry = (double)ny/(double)(NY - sshift - nshift);
// cy = rx*(double)(NY - nshift);
/* build wave table */
for (i=0; i<NX; i++)
for (j=0; j<NY; j++)
{
ii = (int)(rx*(double)(NX-1 - i)) + nx/2;
if (ii > 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<ncircles; k++)
{
pscal = cos(wsphere[i*NY+j].phi - circ_sphere[k].phi)*wsphere[i*NY+j].stheta*sin(circ_sphere[k].theta);
pscal += wsphere[i*NY+j].ctheta*cos(circ_sphere[k].theta);
dist = acos(pscal);
if (dist < circ_sphere[k].radius) return(0);
}
return(1);
}
case (D_SPHERE_JULIA):
{
if (wsphere[i*NY+j].z == 1.0) return(1);
if (wsphere[i*NY+j].z == -1.0) return(0);
r = (1.0 + wsphere[i*NY+j].ctheta)/wsphere[i*NY+j].stheta;
u1 = r*wsphere[i*NY+j].cphi*JULIA_SCALE;
v1 = r*wsphere[i*NY+j].sphi*JULIA_SCALE;
u = u1*cos_rot + v1*sin_rot;
v = -u1*sin_rot + v1*cos_rot;
i = 0;
while ((i<MANDELLEVEL)&&(u*u+v*v < 1000.0*MANDELLIMIT))
{
u1 = u*u - v*v + JULIA_RE;
v = 2.0*u*v + JULIA_IM;
u = u1;
i++;
}
if (u*u + v*v < MANDELLIMIT) return(1);
return(0);
}
case (D_SPHERE_JULIA_INV):
{
if (wsphere[i*NY+j].z == 1.0) return(1);
if (wsphere[i*NY+j].z == -1.0) return(0);
r = (1.0 - wsphere[i*NY+j].ctheta)/wsphere[i*NY+j].stheta;
u1 = r*wsphere[i*NY+j].cphi*JULIA_SCALE;
v1 = r*wsphere[i*NY+j].sphi*JULIA_SCALE;
u = u1*cos_rot + v1*sin_rot;
v = -u1*sin_rot + v1*cos_rot;
i = 0;
while ((i<MANDELLEVEL)&&(u*u+v*v < 1000.0*MANDELLIMIT))
{
u1 = u*u - v*v + JULIA_RE;
v = 2.0*u*v + JULIA_IM;
u = u1;
i++;
}
if (u*u + v*v < MANDELLIMIT) return(0);
return(1);
}
case (D_SPHERE_EARTH):
{
return(wsphere[i*NY+j].indomain);
}
default:
{
printf("Function xy_in_billiard_sphere not defined for this billiard \n");
return(0);
}
}
}
void init_wave_flat_sphere(double phi[NX*NY], double psi[NX*NY], short int xy_in[NX*NY], t_wave_sphere wsphere[NX*NY])
/* initialise field with drop - phi is wave height, psi is phi at time t-1 */
/* phi0 is longidude, theta0 is latitude */
{
int i, j;
printf("Initializing wave\n");
// #pragma omp parallel for private(i,j,xy,dist2)
for (i=0; i<NX; i++)
{
if (i%100 == 0) printf("Initializing column %i of %i\n", i, NX);
for (j=0; j<NY; j++)
{
xy_in[i*NY+j] = xy_in_billiard_sphere(i, j, wsphere);
phi[i*NY+j] = 0.0;
psi[i*NY+j] = 0.0;
}
}
}
void init_circular_wave_sphere(double phi0, double theta0, double phi[NX*NY], double psi[NX*NY], short int xy_in[NX*NY], t_wave_sphere wsphere[NX*NY])
/* initialise field with drop - phi is wave height, psi is phi at time t-1 */
/* phi0 is longidude, theta0 is latitude */
{
int i, j;
double xy[2], pscal, dist, dist2, stheta, ctheta;
ctheta = cos(theta0 + PID);
stheta = sin(theta0 + PID);
printf("Initializing wave\n");
// #pragma omp parallel for private(i,j,xy,dist2)
for (i=0; i<NX; i++)
{
if (i%100 == 0) printf("Initializing column %i of %i\n", i, NX);
for (j=0; j<NY; j++)
{
pscal = cos(wsphere[i*NY+j].phi - phi0)*wsphere[i*NY+j].stheta*stheta;
pscal += wsphere[i*NY+j].ctheta*ctheta;
dist = acos(pscal);
dist2 = dist*dist;
xy_in[i*NY+j] = xy_in_billiard_sphere(i, j, wsphere);
if ((xy_in[i*NY+j])||(TWOSPEEDS))
phi[i*NY+j] = INITIAL_AMP*exp(-dist2/INITIAL_VARIANCE)*cos(-sqrt(dist2)/INITIAL_WAVELENGTH);
else phi[i*NY+j] = 0.0;
psi[i*NY+j] = 0.0;
}
}
}
void add_circular_wave_sphere(double amp, double phi0, double theta0, double phi[NX*NY], double psi[NX*NY], short int xy_in[NX*NY], t_wave_sphere wsphere[NX*NY])
/* initialise field with drop - phi is wave height, psi is phi at time t-1 */
/* phi0 is longidude, theta0 is latitude */
{
int i, j;
double xy[2], pscal, dist, dist2, stheta, ctheta;
ctheta = cos(theta0 + PID);
stheta = sin(theta0 + PID);
printf("Initializing wave\n");
// #pragma omp parallel for private(i,j,xy,dist2)
for (i=0; i<NX; i++)
{
if (i%100 == 0) printf("Initializing column %i of %i\n", i, NX);
for (j=0; j<NY; j++)
{
pscal = cos(wsphere[i*NY+j].phi - phi0)*wsphere[i*NY+j].stheta*stheta;
pscal += wsphere[i*NY+j].ctheta*ctheta;
dist = acos(pscal);
dist2 = dist*dist;
if ((xy_in[i*NY+j])||(TWOSPEEDS))
phi[i*NY+j] += amp*INITIAL_AMP*exp(-dist2/INITIAL_VARIANCE)*cos(-sqrt(dist2)/INITIAL_WAVELENGTH);
else phi[i*NY+j] = 0.0;
psi[i*NY+j] = 0.0;
}
}
}
void compute_light_angle_sphere(short int xy_in[NX*NY], t_wave wave[NX*NY], t_wave_sphere wsphere[NX*NY], int movie)
/* computes cosine of angle between normal vector and vector light */
{
int i, j;
double x, y, z, norm, pscal, deltai[3], deltaj[3], deltar, n[3], r;
static double dphi, dtheta;
static int first = 1;
if (first)
{
dphi = DPI/(double)NX;
dtheta = PI/(double)NY;
first = 0;
}
#pragma omp parallel for private(i,j)
for (i=0; i<NX; i++)
for (j=0; j<NY; j++)
{
r = 1.0 + RSCALE*(*wave[i*NY+j].p_zfield[movie]);
if (r > 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-1; i++)
for (j=0; j<NY-1; j++)
{
if ((TWOSPEEDS)||(xy_in[i*NY+j]))
{
/* computation of tangent vectors */
deltar = (wsphere[(i+1)*NY+j].radius - wsphere[i*NY+j].radius)/dphi;
deltai[0] = -wsphere[i*NY+j].radius*wsphere[i*NY+j].sphi;
deltai[0] += deltar*wsphere[i*NY+j].cphi;
deltai[1] = wsphere[i*NY+j].radius*wsphere[i*NY+j].cphi;
deltai[1] += deltar*wsphere[i*NY+j].sphi;
deltai[2] = -deltar*wsphere[i*NY+j].cottheta;
deltar = (wsphere[i*NY+j+1].radius - wsphere[i*NY+j].radius)/dtheta;
deltaj[0] = wsphere[i*NY+j].radius*wsphere[i*NY+j].cphi*wsphere[i*NY+j].ctheta;
deltaj[0] += deltar*wsphere[i*NY+j].cphi*wsphere[i*NY+j].stheta;
deltaj[1] = wsphere[i*NY+j].radius*wsphere[i*NY+j].sphi*wsphere[i*NY+j].ctheta;
deltaj[1] += deltar*wsphere[i*NY+j].sphi*wsphere[i*NY+j].stheta;
deltaj[2] = wsphere[i*NY+j].radius*wsphere[i*NY+j].stheta;
deltaj[2] += -deltar*wsphere[i*NY+j].ctheta;
/* computation of normal vector */
n[0] = deltai[1]*deltaj[2] - deltai[2]*deltaj[1];
n[1] = deltai[2]*deltaj[0] - deltai[0]*deltaj[2];
n[2] = deltai[0]*deltaj[1] - deltai[1]*deltaj[0];
norm = sqrt(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]);
pscal = n[0]*light[0] + n[1]*light[1] + n[2]*light[2];
wave[i*NY+j].cos_angle = pscal/norm;
}
else
{
pscal = wsphere[i*NY+j].x*light[0] + wsphere[i*NY+j].y*light[1] + wsphere[i*NY+j].z*light[2];
wave[i*NY+j].cos_angle = pscal;
}
}
for (i=0; i<NX-1; i++) wave[i*NY+NY-1].cos_angle = wave[i*NY+NY-2].cos_angle;
for (j=0; j<NY-1; j++) wave[(NX-1)*NY+j].cos_angle = wave[(NX-2)*NY+j].cos_angle;
wave[(NX-1)*NY+NY-1].cos_angle = wave[(NX-1)*NY+NY-2].cos_angle;
}
else
{
#pragma omp parallel for private(i,j,pscal)
for (i=0; i<NX; i++)
for (j=0; j<NY; j++)
{
if ((TWOSPEEDS)||(xy_in[i*NY+j]))
{
pscal = wsphere[i*NY+j].x*light[0] + wsphere[i*NY+j].y*light[1] + wsphere[i*NY+j].z*light[2];
wave[i*NY+j].cos_angle = pscal;
}
}
}
}
void compute_light_angle_sphere_2d(short int xy_in[NX*NY], t_wave wave[NX*NY], int movie)
/* computes cosine of angle between normal vector and vector light */
{
int i, j;
double gradx, grady, norm, pscal;
static double dx, dy, vscale2;
static int first = 1;
if (first)
{
dx = 2.0*(XMAX - XMIN)/(double)NX;
dy = 2.0*(YMAX - YMIN)/(double)NY;
vscale2 = SHADE_SCALE_2D*SHADE_SCALE_2D;
first = 0;
}
#pragma omp parallel for private(i,j,gradx, grady, norm, pscal)
for (i=1; i<NX-1; i++)
for (j=1; j<NY-1; j++)
{
if ((TWOSPEEDS)||(xy_in[i*NY+j]))
{
gradx = (*wave[(i+1)*NY+j].p_zfield[movie] - *wave[(i-1)*NY+j].p_zfield[movie])/dx;
grady = (*wave[i*NY+j+1].p_zfield[movie] - *wave[i*NY+j-1].p_zfield[movie])/dy;
norm = sqrt(vscale2 + gradx*gradx + grady*grady);
pscal = -gradx*light[0] - grady*light[1] + SHADE_SCALE_2D;
wave[i*NY+j].cos_angle = pscal/norm;
}
}
/* i=0 */
for (j=1; j<NY-1; j++)
{
if ((TWOSPEEDS)||(xy_in[j]))
{
gradx = (*wave[NY+j].p_zfield[movie] - *wave[(NY-1)*NY+j].p_zfield[movie])/dx;
grady = (*wave[j+1].p_zfield[movie] - *wave[j-1].p_zfield[movie])/dy;
norm = sqrt(vscale2 + gradx*gradx + grady*grady);
pscal = -gradx*light[0] - grady*light[1] + SHADE_SCALE_2D;
wave[j].cos_angle = pscal/norm;
}
}
/* i=NX-1 */
for (j=1; j<NY-1; j++)
{
if ((TWOSPEEDS)||(xy_in[(NX-1)*NY+j]))
{
gradx = (*wave[j].p_zfield[movie] - *wave[(NX-2)*NY+j].p_zfield[movie])/dx;
grady = (*wave[(NX-1)*NY+j+1].p_zfield[movie] - *wave[(NX-1)*NY+j-1].p_zfield[movie])/dy;
norm = sqrt(vscale2 + gradx*gradx + grady*grady);
pscal = -gradx*light[0] - grady*light[1] + SHADE_SCALE_2D;
wave[(NX-1)*NY+j].cos_angle = pscal/norm;
}
}
}
void compute_cfield_sphere(short int xy_in[NX*NY], int cplot, int palette, t_wave wave[NX*NY], int fade, double fade_value, int movie)
/* compute the colors */
{
int i, j, k;
double ca;
static double fact_max;
static int first = 1;
if (first)
{
fact_max = 0.5*COS_LIGHT_MAX;
first = 0;
}
#pragma omp parallel for private(i,j,ca)
for (i=0; i<NX; i++) for (j=0; j<NY; j++)
{
compute_field_color(*wave[i*NY+j].p_cfield[movie], *wave[i*NY+j].p_cfield[movie+2], cplot, palette, wave[i*NY+j].rgb);
if ((SHADE_3D)||(SHADE_2D))
{
ca = wave[i*NY+j].cos_angle;
// ca = (ca + 1.0)*0.4 + 0.2;
ca = (ca + 1.0)*fact_max + COS_LIGHT_MIN;
if ((FADE_IN_OBSTACLE)&&(!xy_in[i*NY+j])) ca = (ca + 0.1)*1.6;
for (k=0; k<3; k++) wave[i*NY+j].rgb[k] *= ca;
}
if (fade) for (k=0; k<3; k++) wave[i*NY+j].rgb[k] *= fade_value;
}
}
void draw_wave_sphere_2D(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, ii;
double x, y, ca;
static int ishift;
static double dx, dy;
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);
dx = (XMAX - XMIN)/(double)NX;
dy = (YMAX - YMIN)/(double)(NY-2*DPOLE);
ishift = (int)((double)NX*PHISHIFT/360.0);
}
glBegin(GL_QUADS);
for (i=0; i<NX; i++)
for (j=0; j<NY; j++)
{
if ((TWOSPEEDS)||(xy_in[i*NY+j]))
glColor3f(wave[i*NY+j].rgb[0], wave[i*NY+j].rgb[1], wave[i*NY+j].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);
}
x = wsphere[i*NY+j].x2d;
y = wsphere[i*NY+j].y2d;
ii = NX-i-1+ishift;
if (ii > 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<NY-2; 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<NX-1; i++)
for (j=0; j<NY-2; 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=0; j<NY-2; 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=0; j<NY-2; 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<imin; i++)
for (j=0; j<NY-2; 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=0; j<NY-2; 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=0; j<NY-2; 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=0; j<NY-2; 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);
}
/* North pole */
for (i=0; i<NX-1; i++) for (j=NY-2; j<NY-1; 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=NY-2; j<NY-1; j++)
draw_wave_sphere_ij(NX-1, 0, j-1, j, NY-DPOLE, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
}
else
{
if (imin < imax)
{
for (i=imax; 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);
for (i=imax+1; i<NX-1; 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=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<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);
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; i<NX-1; i++) for (j=2; j>0; 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);
}

85
sub_wave.c
View File

@ -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):
{

158
sub_wave_3d.c
View File

@ -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];

12
wave_3d.c
View File

@ -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<NY-1; j++) phi_out[j] = AMPLITUDE*cos((double)time*OMEGA);
if (OSCILLATE_LEFT)
{
for (j=1; j<NY-1; j++) phi_out[j] = oscillating_bc(time);
printf("Boundary condition %.3lg\n", oscillating_bc(time));
for (j=1; j<NY-1; j++) phi_out[j] = oscillating_bc(time, j);
printf("Boundary condition %.3lg\n", oscillating_bc(time, 0));
}
else for (j=1; j<NY-1; j++){
if ((TWOSPEEDS)||(xy_in[j] != 0)){
@ -574,7 +576,7 @@ void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out
// if (OSCILLATE_LEFT) for (j=1; j<NY-1; j++) phi_out[j] = AMPLITUDE*cos((double)time*OMEGA);
if (OSCILLATE_LEFT)
{
for (j=1; j<NY-1; j++) phi_out[j] = oscillating_bc(time);
for (j=1; j<NY-1; j++) phi_out[j] = oscillating_bc(time, j);
// printf("Boundary condition %.3lg\n", oscillating_bc(time));
}
else for (j=1; j<NY-1; j++){
@ -865,8 +867,8 @@ void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out
/* add oscillating boundary condition on the left corners */
if (OSCILLATE_LEFT)
{
phi_out[0] = oscillating_bc(time);
phi_out[NY-1] = oscillating_bc(time);
phi_out[0] = oscillating_bc(time, 0);
phi_out[NY-1] = oscillating_bc(time, NY-1);
}

65
wave_billiard.c
View File

@ -24,7 +24,7 @@
/* create movie using */
/* ffmpeg -i wave.%05d.tif -vcodec libx264 wave.mp4 */
/* */
/******************************F***************************************************/
/**********************************************************************************/
/*********************************************************************************/
/* */
@ -61,10 +61,12 @@
#define NX 3840 /* number of grid points on x axis */
#define NY 2300 /* number of grid points on y axis */
#define XMIN -1.0
#define XMAX 3.0 /* x interval */
#define YMIN -1.197916667
#define YMAX 1.197916667 /* y interval for 9/16 aspect ratio */
#define XMIN -2.0
#define XMAX 2.0 /* x interval */
#define YMIN -0.397916667
#define YMAX 1.997916667 /* y interval for 9/16 aspect ratio */
// #define YMIN -1.197916667
// #define YMAX 1.197916667 /* y interval for 9/16 aspect ratio */
#define HIGHRES 1 /* set to 1 if resolution of grid is double that of displayed image */
@ -87,9 +89,9 @@
/* Choice of the billiard table */
#define B_DOMAIN 10 /* choice of domain shape, see list in global_pdes.c */
#define B_DOMAIN 20 /* choice of domain shape, see list in global_pdes.c */
#define CIRCLE_PATTERN 1 /* pattern of circles or polygons, see list in global_pdes.c */
#define CIRCLE_PATTERN 103 /* pattern of circles or polygons, see list in global_pdes.c */
#define COMPARISON 0 /* set to 1 to compare two different patterns (beta) */
#define B_DOMAIN_B 20 /* second domain shape, for comparisons */
@ -99,8 +101,8 @@
#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.1197916667 /* parameter controlling the dimensions of domain */
#define MU 0.035 /* parameter controlling the dimensions of domain */
#define LAMBDA 1.0 /* parameter controlling the dimensions of domain */
#define MU 0.005 /* 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 6 /* depth of computation of Menger gasket */
@ -108,8 +110,8 @@
#define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */
#define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */
#define FOCI 1 /* set to 1 to draw focal points of ellipse */
#define NGRIDX 30 /* number of grid point for grid of disks */
#define NGRIDY 30 /* number of grid point for grid of disks */
#define NGRIDX 60 /* number of grid point for grid of disks */
#define NGRIDY 10 /* number of grid point for grid of disks */
// #define NGRIDY 18 /* number of grid point for grid of disks */
#define X_SHOOTER -0.2
@ -129,11 +131,11 @@
/* Physical parameters of wave equation */
#define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */
#define OSCILLATE_LEFT 0 /* set to 1 to add oscilating boundary condition on the left */
#define OSCILLATE_LEFT 1 /* set to 1 to add oscilating boundary condition on the left */
#define OSCILLATE_TOPBOT 0 /* set to 1 to enforce a planar wave on top and bottom boundary */
#define OSCILLATION_SCHEDULE 0 /* oscillation schedule, see list in global_pdes.c */
#define OMEGA 0.004 /* frequency of periodic excitation */
#define OMEGA 0.024 /* frequency of periodic excitation */
#define AMPLITUDE 1.0 /* amplitude of periodic excitation */
#define ACHIRP 0.25 /* acceleration coefficient in chirp */
#define DAMPING 0.0 /* damping of periodic excitation */
@ -146,12 +148,13 @@
#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 -400.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 1 /* set to 1 to add an oscillating wave source */
#define ADD_OSCILLATING_SOURCE 0 /* set to 1 to add an oscillating wave source */
#define OSCILLATING_SOURCE_PERIOD 30 /* period of oscillating source */
#define ALTERNATE_OSCILLATING_SOURCE 1 /* set to 1 to alternate sign of oscillating source */
@ -166,11 +169,11 @@
/* Parameters for length and speed of simulation */
#define NSTEPS 2400 /* number of frames of movie */
// #define NSTEPS 1300 /* number of frames of movie */
#define NVID 6 /* number of iterations between images displayed on screen */
#define NSTEPS 3600 /* number of frames of movie */
// #define NSTEPS 500 /* number of frames of movie */
#define NVID 7 /* number of iterations between images displayed on screen */
#define NSEG 1000 /* number of segments of boundary */
#define INITIAL_TIME 0 /* time after which to start saving frames */
#define INITIAL_TIME 700 /* time after which to start saving frames */
#define BOUNDARY_WIDTH 2 /* width of billiard boundary */
#define PRINT_SPEED 0 /* print speed of moving source */
#define PRINT_FREQUENCY 0 /* print frequency (for phased array) */
@ -199,9 +202,9 @@
/* Color schemes */
// #define COLOR_PALETTE 15 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE 17 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 13 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE 18 /* Color palette, see list in global_pdes.c */
// #define COLOR_PALETTE 17 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 12 /* Color palette, see list in global_pdes.c */
#define BLACK 1 /* background */
@ -212,9 +215,9 @@
#define PHASE_FACTOR 1.0 /* factor in computation of phase in color scheme P_3D_PHASE */
#define PHASE_SHIFT 0.0 /* shift of phase in color scheme P_3D_PHASE */
#define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */
#define E_SCALE 100.0 /* scaling factor for energy representation */
#define E_SCALE 50.0 /* scaling factor for energy representation */
#define LOG_SCALE 1.0 /* scaling factor for energy log representation */
#define LOG_SHIFT 1.0 /* shift of colors on log scale */
#define LOG_SHIFT 1.5 /* shift of colors on log scale */
#define FLUX_SCALE 5.0e3 /* scaling factor for enegy flux represtnation */
#define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */
@ -226,12 +229,14 @@
#define HUEAMP -180.0 /* amplitude of variation of hue for color scheme C_HUE */
#define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */
#define COLORBAR_RANGE 2.0 /* scale of color scheme bar */
#define COLORBAR_RANGE_B 2.0 /* scale of color scheme bar for 2nd part */
#define COLORBAR_RANGE 1.0 /* scale of color scheme bar */
#define COLORBAR_RANGE_B 0.4 /* 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 1 /* 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 NXMAZE 8 /* width of maze */
@ -328,7 +333,7 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
}
/* left boundary */
if (OSCILLATE_LEFT) for (j=1; j<NY-1; j++) phi_out[0][j] = oscillating_bc(time);
if (OSCILLATE_LEFT) for (j=1; j<NY-1; j++) phi_out[0][j] = oscillating_bc(time, j);
else for (j=1; j<NY-1; j++){
if ((TWOSPEEDS)||(xy_in[0][j] != 0)){
x = phi_in[0][j];
@ -519,8 +524,8 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
/* add oscillating boundary condition on the left corners */
if (OSCILLATE_LEFT)
{
phi_out[0][0] = oscillating_bc(time);
phi_out[0][NY-1] = oscillating_bc(time);
phi_out[0][0] = oscillating_bc(time, 0);
phi_out[0][NY-1] = oscillating_bc(time, NY-1);
}
/* for debugging purposes/if there is a risk of blow-up */
@ -676,9 +681,9 @@ void animation()
// xy_to_ij(startright[0], startright[1], sample_right);
// printf("xleft = (%.3f, %.3f) xright = (%.3f, %.3f)\n", xin_left, yin_left, xin_right, yin_right);
// init_wave_flat(phi, psi, xy_in);
init_wave_flat(phi, psi, xy_in);
init_circular_wave(-0.5, 0.0, phi, psi, xy_in);
// init_circular_wave(-0.5, 0.0, phi, psi, xy_in);
// x = XMIN + (XMAX - XMIN)*rand()/RAND_MAX;
// y = YMIN + (YMAX - YMIN)*rand()/RAND_MAX;
// init_circular_wave(0.0, -0.8, phi, psi, xy_in);

149
wave_common.c
View File

@ -853,52 +853,42 @@ void draw_wave_epalette(double *phi[NX], double *psi[NX], double *total_energy[N
}
void draw_wave_highres_palette(int size, double *phi[NX], double *psi[NX], double *total_energy[NX], double *total_flux, short int *xy_in[NX], double scale, int time, int plot, int palette, int fade, double fade_value)
/* same as draw_wave_highres, but with color scheme option */
double wave_value(int i, int j, double *phi[NX], double *psi[NX], double *total_energy[NX], double *total_flux, short int *xy_in[NX], double scale, int time, int plot, int palette, double rgb[3])
/* compute value of wave and color */
{
int i, j, k, iplus, iminus, jplus, jminus;
double rgb[3], xy[2], x1, y1, x2, y2, velocity, energy, gradientx2, gradienty2, arg, mod, flux_factor, gx, gy, mgx, mgy;
static double dtinverse = ((double)NX)/(COURANT*(XMAX-XMIN)), dx = (XMAX-XMIN)/((double)NX);
int k;
double value, velocity, energy, gradientx2, gradienty2, arg, mod, flux_factor, gx, gy, mgx, mgy;
glBegin(GL_QUADS);
// printf("dtinverse = %.5lg\n", dtinverse);
for (i=0; i<NX; i+=size)
for (j=0; j<NY; j+=size)
{
if ((TWOSPEEDS)||(xy_in[i][j]))
{
switch (plot) {
case (P_AMPLITUDE):
{
// /* make wave luminosity larger inside obstacles */
// if (!(xy_in[i][j])) color_scheme_lum(COLOR_SCHEME, phi[i][j], scale, time, 0.7, rgb);
// else
color_scheme_palette(COLOR_SCHEME, palette, phi[i][j], scale, time, rgb);
value = phi[i][j];
color_scheme_palette(COLOR_SCHEME, palette, value, scale, time, rgb);
break;
}
case (P_ENERGY):
{
energy = compute_energy(phi, psi, xy_in, i, j);
value = compute_energy(phi, psi, xy_in, i, j);
/* adjust energy to color palette */
if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, energy, scale, time, rgb);
else color_scheme_palette(COLOR_SCHEME, palette, energy, scale, time, rgb);
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_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);
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, 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);
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):
@ -906,7 +896,8 @@ void draw_wave_highres_palette(int size, double *phi[NX], double *psi[NX], doubl
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);
value = LOG_SHIFT + LOG_SCALE*log(energy);
color_scheme_palette(COLOR_SCHEME, palette, value, scale, time, rgb);
break;
}
case (P_LOG_MEAN_ENERGY):
@ -914,7 +905,8 @@ void draw_wave_highres_palette(int size, double *phi[NX], double *psi[NX], doubl
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);
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):
@ -923,7 +915,8 @@ void draw_wave_highres_palette(int size, double *phi[NX], double *psi[NX], doubl
// 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);
value = mod*FLUX_SCALE;
color_scheme_asym_palette(COLOR_SCHEME, palette, value, scale, time, rgb);
break;
}
case (P_TOTAL_ENERGY_FLUX):
@ -942,13 +935,103 @@ void draw_wave_highres_palette(int size, double *phi[NX], double *psi[NX], doubl
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);
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<imax; i+=size)
{
if (values[i*NY+jmin] > 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<imax; i+=size)
{
y = a*values[i*NY+jmin] + b;
glVertex2d((double)i, y);
}
glEnd();
glBegin(GL_LINE_LOOP);
glVertex2i(imin, jmin);
glVertex2i(imax, jmin);
glVertex2i(imax, jmax);
glVertex2i(imin, jmax);
glEnd();
}
void draw_wave_highres_palette(int size, double *phi[NX], double *psi[NX], double *total_energy[NX], double *total_flux, short int *xy_in[NX], double scale, int time, int plot, int palette, int fade, double fade_value)
/* same as draw_wave_highres, but with color scheme option */
{
int i, j, k, iplus, iminus, jplus, jminus;
double value, rgb[3], xy[2], x1, y1, x2, y2, velocity, energy, gradientx2, gradienty2, arg, mod, flux_factor, gx, gy, mgx, mgy;
// double vmin, vmax, deltav;
static double dtinverse = ((double)NX)/(COURANT*(XMAX-XMIN)), dx = (XMAX-XMIN)/((double)NX);
double *values;
if (DRAW_WAVE_PROFILE) values = (double *)malloc(NX*NY*sizeof(double));
glBegin(GL_QUADS);
// printf("dtinverse = %.5lg\n", dtinverse);
for (i=0; i<NX; i+=size)
for (j=0; j<NY; j+=size)
{
if ((TWOSPEEDS)||(xy_in[i][j]))
{
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 */

4
wave_comparison.c
View File

@ -43,7 +43,7 @@
#include <omp.h>
#include <time.h>
#define MOVIE 1 /* set to 1 to generate movie */
#define MOVIE 0 /* set to 1 to generate movie */
#define DOUBLE_MOVIE 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 */

8
wave_energy.c
View File

@ -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 */

944
wave_sphere.c Normal file
View File

@ -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 <math.h>
#include <string.h>
#include <GL/glut.h>
#include <GL/glu.h>
#include <unistd.h>
#include <sys/types.h>
#include <tiffio.h> /* Sam Leffler's libtiff library. */
#include <omp.h>
#include <time.h>
#define MOVIE 0 /* set to 1 to generate movie */
#define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */
#define 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<NY-DPOLE; j++){
sintheta = sin(j*dtheta);
// invstheta = 1.0/(sintheta*sintheta);
invstheta = 1.0/(sintheta*sintheta + SMOOTHPOLE*SMOOTHPOLE);
// cottheta = ctheta*cos(j*dtheta)/sintheta;
cottheta = ctheta*cos(j*dtheta)/(sintheta + SMOOTHPOLE);
for (i=1; i<NX-1; i++){
if ((TWOSPEEDS)||(xy_in[i*NY+j] != 0)){
x = phi_in[i*NY+j];
y = psi_in[i*NY+j];
/* discretized Laplacian */
/* 2nd phi derivative */
delta = invstheta*(phi_in[(i+1)*NY+j] + phi_in[(i-1)*NY+j] - 2.0*x);
/* 2nd theta derivative */
delta += cphiphi*(phi_in[i*NY+j+1] + phi_in[i*NY+j-1] - 2.0*x);
/* first theta derivative */
delta += cottheta*(phi_in[i*NY+j+1] - phi_in[i*NY+j-1]);
/* evolve phi */
phi_out[i*NY+j] = -y + 2*x + tcc[i*NY+j]*delta - KAPPA*x - tgamma[i*NY+j]*(x-y);
}
}
}
/* evolution at longitude zero */
for (j=DPOLE; j<NY-DPOLE; j++){
sintheta = sin(j*dtheta);
invstheta = 1.0/(sintheta*sintheta);
cottheta = ctheta*cos(j*dtheta)/sintheta;
/* i = 0 */
if ((TWOSPEEDS)||(xy_in[j] != 0)){
x = phi_in[j];
y = psi_in[j];
/* discretized Laplacian */
/* 2nd phi derivative */
delta = invstheta*(phi_in[NY+j] + phi_in[(NX-1)*NY+j] - 2.0*x);
/* 2nd theta derivative */
delta += cphiphi*(phi_in[j+1] + phi_in[j-1] - 2.0*x);
/* first theta derivative */
delta += cottheta*(phi_in[j+1] - phi_in[j-1]);
/* evolve phi */
phi_out[j] = -y + 2*x + tcc[j]*delta - KAPPA*x - tgamma[j]*(x-y);
}
/* i = NX-1 */
if ((TWOSPEEDS)||(xy_in[(NX-1)*NY+j] != 0)){
x = phi_in[(NX-1)*NY+j];
y = psi_in[(NX-1)*NY+j];
/* discretized Laplacian */
/* 2nd phi derivative */
delta = invstheta*(phi_in[j] + phi_in[(NX-2)*NY+j] - 2.0*x);
/* 2nd theta derivative */
delta += cphiphi*(phi_in[(NX-1)*NY+j+1] + phi_in[(NX-1)*NY+j-1] - 2.0*x);
/* first theta derivative */
delta += cottheta*(phi_in[(NX-1)*NY+j+1] - phi_in[(NX-1)*NY+j-1]);
/* evolve phi */
phi_out[(NX-1)*NY+j] = -y + 2*x + tcc[(NX-1)*NY+j]*delta - KAPPA*x - tgamma[(NX-1)*NY+j]*(x-y);
}
}
/* compute average at north pole */
sum = 0.0;
for (i=0; i<NX; i++) sum += phi_out[i*NY + DPOLE];
avrg = sum/(double)NX;
for (i=0; i<NX; i++) for (j=0; j<DPOLE; j++)
phi_out[i*NY + j] = avrg;
// {
// x = phi_in[(NX-1)*NY+j];
// y = psi_in[(NX-1)*NY+j];
// phi_out[(NX-1)*NY+j] = -y + 2*x + tcc[(NX-1)*NY+j]*avrg - KAPPA*x - tgamma[(NX-1)*NY+j]*(x-y);
// //
// }
/* compute average at south pole */
sum = 0.0;
for (i=0; i<NX; i++) sum += phi_out[i*NY + NY-1-DPOLE];
avrg = sum/(double)NX;
for (i=0; i<NX; i++) for (j=NY-DPOLE; j<NY; j++)
phi_out[i*NY + j] = avrg;
// {
// x = phi_in[(NX-1)*NY+j];
// y = psi_in[(NX-1)*NY+j];
// phi_out[(NX-1)*NY+j] = -y + 2*x + tcc[(NX-1)*NY+j]*avrg - KAPPA*x - tgamma[(NX-1)*NY+j]*(x-y);
// //
// }
// /* north pole, j = 0 */
// if ((TWOSPEEDS)||(xy_in[0] != 0)){
// x = phi_in[0];
// y = psi_in[0];
//
// /* discretized Laplacian */
// delta = cphiphi*(phi_in[1] + phi_in[(NX/4)*NY+1] + phi_in[(NX/2)*NY+1] + phi_in[(3*NX/4)*NY+1] - 4.0*x);
//
// /* evolve phi */
// phi_out[0] = -y + 2*x + tcc[0]*delta - KAPPA*x - tgamma[0]*(x-y);
//
// /* set same values for all i */
// for (i=1; i<NX; i++) phi_out[i*NY] = phi_out[0];
// }
//
// /* south pole, j = NY-1 */
// if ((TWOSPEEDS)||(xy_in[NY-1] != 0)){
// x = phi_in[NY-1];
// y = psi_in[NY-1];
//
// /* discretized Laplacian */
// delta = cphiphi*(phi_in[NY-2] + phi_in[(NX/4)*NY+NY-2] + phi_in[(NX/2)*NY+NY-2] + phi_in[(3*NX/4)*NY+NY-2] - 4.0*x);
//
// /* evolve phi */
// phi_out[NY-1] = -y + 2*x + tcc[0]*delta - KAPPA*x - tgamma[0]*(x-y);
//
// /* set same values for all i */
// for (i=1; i<NX; i++) phi_out[i*NY+NY-1] = phi_out[NY-1];
// }
/* for debugging purposes/if there is a risk of blow-up */
if (FLOOR) for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
if (xy_in[i*NY+j] != 0)
{
if (phi_out[i*NY+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<npolyline; i++) printf("vertex %i: (%.3f, %.3f)\n", i, polyline[i].x, polyline[i].y);
if (COMPARISON) npolyline_b = init_polyline(MDEPTH, polyline);
npolyrect = init_polyrect(polyrect);
for (i=0; i<npolyrect; i++) printf("polyrect vertex %i: (%.3f, %.3f) - (%.3f, %.3f)\n", i, polyrect[i].x1, polyrect[i].y1, polyrect[i].x2, polyrect[i].y2);
init_polyrect_arc(polyrectrot, polyarc, &npolyrect_rot, &npolyarc);
printf("Rotated rectangles and arcs initialized\n");
printf("%i rotated rectangles, %i arcs\n", npolyrect_rot, npolyarc);
if ((B_DOMAIN == D_SPHERE_CIRCLES)||(B_DOMAIN_B == D_SPHERE_CIRCLES))
{
ncircles = init_circle_sphere(circ_sphere, CIRCLE_PATTERN);
}
courant2 = COURANT*COURANT;
courantb2 = COURANTB*COURANTB;
c = COURANT*(XMAX - XMIN)/(double)NX;
// a = 0.015;
// b = 0.0003;
// a = 0.04;
// b = 0.0018;
a = 0.05;
b = 0.0016;
/* initialize color scale, for option RESCALE_COLOR_IN_CENTER */
if (RESCALE_COLOR_IN_CENTER)
{
for (i=0; i<NX; i++)
for (j=0; j<NY; j++)
{
ij_to_xy(i, j, xy);
r2 = xy[0]*xy[0] + xy[1]*xy[1];
color_scale[i*NY+j] = 1.0 - exp(-4.0*r2/LAMBDA*LAMBDA);
}
}
/* initialize wave with a drop at one point, zero elsewhere */
// init_circular_wave(0.0, -LAMBDA, phi, psi, xy_in);
init_wave_fields(wave);
init_wave_sphere(wsphere);
/* initialize total energy table - no longer needed */
// if ((ZPLOT == P_MEAN_ENERGY)||(ZPLOT_B == P_MEAN_ENERGY)||(ZPLOT == P_LOG_MEAN_ENERGY)||(ZPLOT_B == P_LOG_MEAN_ENERGY))
// for (i=0; i<NX; i++)
// for (j=0; j<NY; j++)
// total_energy[i*NY+j] = 0.0;
ratio = (XMAX - XMIN)/8.4; /* for Tokarsky billiard */
// init_circular_wave_mod(polyline[85].x, polyline[85].y, phi, psi, xy_in);
// init_circular_wave_mod(LAMBDA*cos(APOLY*PID), LAMBDA*sin(APOLY*PID), phi, psi, xy_in);
lambda1 = LAMBDA;
angle = DPI/(double)NPOLY;
// init_circular_wave_mod(lambda1*cos(0.5*angle), lambda1*sin(0.5*angle), phi, psi, xy_in);
// for (j=1; j<NPOLY; j++)
// add_circular_wave_mod(1.0, lambda1*cos(((double)j+0.5)*angle), lambda1*sin(((double)j+0.5)*angle), phi, psi, xy_in);
init_circular_wave_sphere(0.7, 0.5, phi, psi, xy_in, wsphere);
// init_wave_flat_sphere(phi, psi, xy_in, wsphere);
// init_circular_wave_sphere(0.25 + PID, 0.0, phi, psi, xy_in, wsphere);
// add_circular_wave_sphere(-1.0, 0.25 + 3.0*PID, 0.0, phi, psi, xy_in, wsphere);
// printf("Wave initialized\n");
/* initialize table of wave speeds/dissipation */
init_speed_dissipation(xy_in, tc, tcc, tgamma);
/* initialze potential to add to z coordinate */
if (ADD_POTENTIAL)
{
if (POTENTIAL == POT_IOR)
for (i=0; i<NX*NY; i++)
wave[i].potential = &tcc[i];
}
init_zfield(phi, psi, xy_in, ZPLOT, wave, 0);
init_cfield(phi, psi, xy_in, CPLOT, wave, 0);
if (DOUBLE_MOVIE)
{
init_zfield(phi, psi, xy_in, ZPLOT_B, wave, 1);
init_cfield(phi, psi, xy_in, CPLOT_B, wave, 1);
}
blank();
glColor3f(0.0, 0.0, 0.0);
draw_wave_sphere(0, phi, psi, xy_in, wave, wsphere, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1);
// draw_billiard();
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, 0, 1.0);
glutSwapBuffers();
sleep(SLEEP1);
for (i=0; i<=INITIAL_TIME + NSTEPS; i++)
{
global_time++;
//printf("%d\n",i);
/* compute the variance of the field to adjust color scheme */
/* the color depends on the field divided by sqrt(1 + variance) */
if (SCALE)
{
scale = sqrt(1.0 + compute_variance_mod(phi,psi, xy_in));
// printf("Scaling factor: %5lg\n", scale);
}
else scale = 1.0;
if (ROTATE_VIEW)
{
viewpoint_schedule(i - INITIAL_TIME);
reset_view = 1;
}
draw_wave_sphere(0, phi, psi, xy_in, wave, wsphere, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1);
for (j=0; j<NVID; j++)
{
evolve_wave(phi, psi, tmp, xy_in, tc, tcc, tgamma, laplace, laplace1, laplace2);
if (SAVE_TIME_SERIES)
{
wave_value = (long int)(phi[sample_left[0]*NY+sample_left[1]]*1.0e16);
fprintf(time_series_left, "%019ld\n", wave_value);
wave_value = (long int)(phi[sample_right[0]*NY+sample_right[1]]*1.0e16);
fprintf(time_series_right, "%019ld\n", wave_value);
if ((j == 0)&&(i%10 == 0)) printf("Frame %i of %i\n", i, NSTEPS);
// fprintf(time_series_right, "%.15f\n", phi[sample_right[0]][sample_right[1]]);
}
// if (i % 10 == 9) oscillate_linear_wave(0.2*scale, 0.15*(double)(i*NVID + j), -1.5, YMIN, -1.5, YMAX, phi, psi);
}
// draw_billiard();
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, fade, fade_value);
/* add oscillating waves */
// if ((ADD_OSCILLATING_SOURCE)&&(i%OSCILLATING_SOURCE_PERIOD == OSCILLATING_SOURCE_PERIOD - 1))
if ((ADD_OSCILLATING_SOURCE)&&(i%OSCILLATING_SOURCE_PERIOD == 1))
{
if (ALTERNATE_OSCILLATING_SOURCE) sign = -sign;
add_circular_wave_mod(sign, -0.5, 0.0, phi, psi, xy_in);
}
if (PRINT_SPEED) print_speed_3d(speed, 0, 1.0);
if (!((NO_EXTRA_BUFFER_SWAP)&&(MOVIE))) glutSwapBuffers();
if (MOVIE)
{
if (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<MID_FRAMES; i++) save_frame();
else for (i=0; i<MID_FRAMES; i++)
{
if ((ROTATE_VIEW)&&(ROTATE_VIEW_WHILE_FADE))
{
viewpoint_schedule(NSTEPS - INITIAL_TIME + i);
reset_view = 1;
}
fade_value = 1.0 - (double)i/(double)MID_FRAMES;
draw_wave_sphere(0, phi, psi, xy_in, wave, wsphere, ZPLOT, CPLOT, COLOR_PALETTE, 1, fade_value, 1);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, 1, fade_value);
if (PRINT_SPEED) print_speed_3d(speed, 1, fade_value);
if (!NO_EXTRA_BUFFER_SWAP) glutSwapBuffers();
save_frame_counter(NSTEPS + i + 1);
}
if ((ROTATE_VIEW)&&(ROTATE_VIEW_WHILE_FADE))
{
viewpoint_schedule(NSTEPS - INITIAL_TIME);
reset_view = 1;
}
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();
if (!FADE) for (i=0; i<END_FRAMES; i++) save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter + i);
else for (i=0; i<END_FRAMES; i++)
{
if ((ROTATE_VIEW)&&(ROTATE_VIEW_WHILE_FADE))
{
viewpoint_schedule(NSTEPS - INITIAL_TIME + i);
reset_view = 1;
}
fade_value = 1.0 - (double)i/(double)END_FRAMES;
draw_wave_sphere(1, phi, psi, xy_in, wave, wsphere, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 1, fade_value, 1);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, CIRC_COLORBAR_B, 1, fade_value);
if (PRINT_SPEED) print_speed_3d(speed, 1, fade_value);
glutSwapBuffers();
save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter + i);
}
}
else
{
if (!FADE) for (i=0; i<END_FRAMES; i++) save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter + i);
else for (i=0; i<END_FRAMES; i++)
{
if ((ROTATE_VIEW)&&(ROTATE_VIEW_WHILE_FADE))
{
viewpoint_schedule(NSTEPS - INITIAL_TIME + i);
reset_view = 1;
}
fade_value = 1.0 - (double)i/(double)END_FRAMES;
draw_wave_sphere(0, phi, psi, xy_in, wave, wsphere, ZPLOT, CPLOT, COLOR_PALETTE, 1, fade_value, 1);
if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, 1, fade_value);
if (PRINT_SPEED) print_speed_3d(speed, 1, fade_value);
glutSwapBuffers();
save_frame_counter(NSTEPS + 1 + counter + i);
}
}
s = system("mv wave*.tif tif_wave/");
}
free(xy_in);
free(phi);
free(psi);
free(tmp);
// free(total_energy);
free(color_scale);
free(tc);
free(tcc);
free(tgamma);
free(wave);
free(laplace);
free(laplace1);
free(laplace2);
if (SAVE_TIME_SERIES)
{
fclose(time_series_left);
fclose(time_series_right);
}
}
void display(void)
{
time_t rawtime;
struct tm * timeinfo;
time(&rawtime);
timeinfo = localtime(&rawtime);
glPushMatrix();
blank();
glutSwapBuffers();
blank();
glutSwapBuffers();
animation();
sleep(SLEEP2);
glPopMatrix();
glutDestroyWindow(glutGetWindow());
printf("Start local time and date: %s", asctime(timeinfo));
time(&rawtime);
timeinfo = localtime(&rawtime);
printf("Current local time and date: %s", asctime(timeinfo));
}
int main(int argc, char** argv)
{
glutInit(&argc, argv);
glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH);
glutInitWindowSize(WINWIDTH,WINHEIGHT);
glutCreateWindow("Wave equation in a planar domain");
if (PLOT_2D) init_sphere_2D();
else init_3d();
glutDisplayFunc(display);
glutMainLoop();
return 0;
}