diff --git a/Earth_Map_Blue_Marble_2002_large.ppm.gz b/Earth_Map_Blue_Marble_2002_large.ppm.gz index ea9834e..e2c3764 100644 Binary files a/Earth_Map_Blue_Marble_2002_large.ppm.gz and b/Earth_Map_Blue_Marble_2002_large.ppm.gz differ diff --git a/Mars_Viking_ClrMosaic_global_925m_scaled.ppm.gz b/Mars_Viking_ClrMosaic_global_925m_scaled.ppm.gz new file mode 100644 index 0000000..9d07f36 Binary files /dev/null and b/Mars_Viking_ClrMosaic_global_925m_scaled.ppm.gz differ diff --git a/Moon_LRO_LOLA_global_LDEM_1024.ppm.gz b/Moon_LRO_LOLA_global_LDEM_1024.ppm.gz new file mode 100644 index 0000000..34286aa Binary files /dev/null and b/Moon_LRO_LOLA_global_LDEM_1024.ppm.gz differ diff --git a/Moon_photo_map.ppm.gz b/Moon_photo_map.ppm.gz new file mode 100644 index 0000000..2e9634b Binary files /dev/null and b/Moon_photo_map.ppm.gz differ diff --git a/bathymetry_gebco_2560_1280_mod2_color.ppm.gz b/bathymetry_gebco_2560_1280_mod2_color.ppm.gz new file mode 100644 index 0000000..6eb5ee6 Binary files /dev/null and b/bathymetry_gebco_2560_1280_mod2_color.ppm.gz differ diff --git a/billiard_phase_space.c b/billiard_phase_space.c index dfcf908..ce0d50d 100644 --- a/billiard_phase_space.c +++ b/billiard_phase_space.c @@ -30,7 +30,7 @@ #include /* Sam Leffler's libtiff library. */ #include -#define MOVIE 1 /* set to 1 to generate movie */ +#define MOVIE 0 /* set to 1 to generate movie */ #define SAVE_MEMORY 1 /* set to 1 to save memory when writing tiff images */ #define WINWIDTH 1280 /* window width */ @@ -52,7 +52,7 @@ /* Choice of the billiard table, see global_particles.c */ -#define B_DOMAIN 11 /* choice of domain shape */ +#define B_DOMAIN 1 /* choice of domain shape */ #define CIRCLE_PATTERN 6 /* pattern of circles */ #define POLYLINE_PATTERN 4 /* pattern of polyline */ @@ -67,8 +67,8 @@ #define NGOLDENSPIRAL 2000 /* max number of points for C_GOLDEN_SPIRAL arrandement */ #define SDEPTH 2 /* Sierpinski gastket depth */ -#define LAMBDAMIN 6.0 /* parameter controlling shape of domain */ -#define LAMBDA 20.0 /* parameter controlling shape of domain */ +#define LAMBDAMIN 0.0 /* parameter controlling shape of domain */ +#define LAMBDA 1.5 /* parameter controlling shape of domain */ #define MU 1.0 /* second parameter controlling shape of billiard */ #define FOCI 1 /* set to 1 to draw focal points of ellipse */ #define NPOLY 6 /* number of sides of polygon */ @@ -96,7 +96,7 @@ #define PRINT_TRAJECTORY_LENGTH 0 /* set to 1 to print length of trajectory 0 */ #define PRINT_TRAJECTORY_PERIOD 0 /* set to 1 to print period of trajectory 0 */ #define DRAW_LENGTHS_PLOT 0 /* set to 1 to plot trajectory lengths */ -#define LENGTHS_LOG_SCALE 1 /* set to 1 to use log scale for plot of lengths */ +#define LENGTHS_LOG_SCALE 0 /* set to 1 to use log scale for plot of lengths */ #define MAX_ANGLE 90.0 /* range of angles of trajectory */ #define NSTEPS 4000 /* number of frames of movie */ @@ -141,7 +141,7 @@ #define PLOT_LYAPUNOV 1 /* set to 1 to add plot of Lyapunov exponents */ #define LOGSCALEX_LYAP 0 /* set to 1 to use log scale on parameter axis of Lyapunov exponents */ #define LYAP_MAX 1.0 /* maximal Lyapunov exponent */ -#define ADAPT_TO_SYMMETRY 1 /* set to 1 to show only one symmetric part of phase space */ +#define ADAPT_TO_SYMMETRY 0 /* set to 1 to show only one symmetric part of phase space */ #define SYMMETRY_FACTOR 3 /* proportion of phase space to be shown */ #define PAUSE 1000 /* number of frames after which to pause */ diff --git a/digital_elevation_model_large.ppm.gz b/digital_elevation_model_large.ppm.gz new file mode 100644 index 0000000..01b6432 Binary files /dev/null and b/digital_elevation_model_large.ppm.gz differ diff --git a/global_3d.c b/global_3d.c index 16ca4f8..7e15f8f 100644 --- a/global_3d.c +++ b/global_3d.c @@ -65,6 +65,14 @@ #define VP_HORIZONTAL 0 /* rotate in a horizontal plane */ #define VP_ORBIT 1 /* rotate in a plane containing the origin */ +#define VP_ORBIT2 11 /* rotate in a plane specified by max latitude */ +#define VP_POLAR 2 /* polar orbit */ + +/* Type of digital elevation model */ + +#define DEM_EARTH 0 /* DEM of Earth */ +#define DEM_MARS 1 /* DEM of Mars */ +#define DEM_MOON 2 /* DEM of the Moon */ /* macros to avoid unnecessary computations in 3D plots */ @@ -85,6 +93,9 @@ #define COMPUTE_TOTAL_ENERGY ((ZPLOT == P_3D_TOTAL_ENERGY)||(CPLOT == P_3D_TOTAL_ENERGY)||(ZPLOT == P_3D_LOG_TOTAL_ENERGY)||(CPLOT == P_3D_LOG_TOTAL_ENERGY)||(ZPLOT == P_3D_MEAN_ENERGY)||(CPLOT == P_3D_MEAN_ENERGY)||(ZPLOT == P_3D_LOG_MEAN_ENERGY)||(CPLOT == P_3D_LOG_MEAN_ENERGY)||(ZPLOT_B == P_3D_TOTAL_ENERGY)||(CPLOT_B == P_3D_TOTAL_ENERGY)||(ZPLOT_B == P_3D_LOG_TOTAL_ENERGY)||(CPLOT_B == P_3D_LOG_TOTAL_ENERGY)||(ZPLOT_B == P_3D_MEAN_ENERGY)||(CPLOT_B == P_3D_MEAN_ENERGY)||(ZPLOT_B == P_3D_LOG_MEAN_ENERGY)||(CPLOT_B == P_3D_LOG_MEAN_ENERGY)) +#define PLANET ((B_DOMAIN == D_SPHERE_EARTH)||(B_DOMAIN == D_SPHERE_MARS)||(B_DOMAIN == D_SPHERE_MOON)) +#define OTHER_PLANET ((B_DOMAIN == D_SPHERE_MARS)||(B_DOMAIN == D_SPHERE_MOON)) + #define NMAXCIRC_SPHERE 100 /* max number of circles on sphere */ int global_time = 0; @@ -154,6 +165,9 @@ typedef struct double r, g, b; /* RGB values for image */ short int indomain; /* has value 1 if lattice point is in domain */ double x2d, y2d; /* x and y coordinates for 2D representation */ + double altitude; /* altitude in case of Earth with digital elevation model */ + double cos_angle; /* cosine of light angle */ + double cos_angle_sphere; /* cosing of light angle for perfect sphere */ } t_wave_sphere; diff --git a/global_ljones.c b/global_ljones.c index 4d76d63..fe6d9b5 100644 --- a/global_ljones.c +++ b/global_ljones.c @@ -75,8 +75,10 @@ #define S_BIN_OPENING 20 /* bin containing particles opening at deactivation time */ #define S_POLYGON_EXT 21 /* exterior of a regular polygon */ #define S_WEDGE_EXT 22 /* exterior of a wedge */ -#define S_MIXER 23 /* exterior of a set of fan of rectangles */ +#define S_MIXER 23 /* exterior of a blender made of rectangles */ #define S_AIRFOIL 24 /* exterior of an air foil */ +#define S_COANDA 25 /* wall for Coanda effect */ +#define S_COANDA_SHORT 26 /* shorter wall for Coanda effect */ /* particle interaction */ @@ -92,6 +94,7 @@ #define I_VICSEK_REPULSIVE 9 /* Vicsek-type interaction with harmonic repulsion */ #define I_VICSEK_SPEED 10 /* Vicsek-type interaction with speed adjustment */ #define I_VICSEK_SHARK 11 /* Vicsek-type interaction with speed adjustment, and one shark */ +#define I_COULOMB_LJ 12 /* Coulomb force regularised by Lennard-Jones repulsion */ /* Boundary conditions */ @@ -178,6 +181,17 @@ #define IC_TWOROCKETS 8 /* type 1 or 2 depending on rocket position */ #define IC_TWOROCKETS_TWOFUELS 9 /* type 1 and 2 or 1 and 3 depending on rocket */ +/* Initial conditions for option TWO_TYPES */ + +#define TTC_RANDOM 0 /* assign types randomly */ +#define TTC_CHESSBOARD 1 /* assign types according to chessboard, works with hex initial config */ +#define TTC_COANDA 2 /* type 1 in a band of width LAMBDA */ + +/* Initial speed distribution */ + +#define VI_RANDOM 0 /* random (Gaussian) initial speed distribution */ +#define VI_COANDA 1 /* nonzero speed in a band of width LAMBDA */ + /* Plot types */ #define P_KINETIC 0 /* colors represent kinetic energy of particles */ @@ -197,6 +211,11 @@ #define P_DIRECT_EMEAN 14 /* averaged version of P_DIRECT_ENERGY */ #define P_NOPARTICLE 15 /* particles are not drawn (only the links between them) */ +/* Rotation schedules */ + +#define ROT_SPEEDUP_SLOWDOWN 0 /* rotation speeds up and then slows down to zero */ +#define ROT_BACK_FORTH 1 /* rotation goes in one direction and then back */ + /* Initial position dependence types */ #define IP_X 0 /* color depends on x coordinate of initial position */ @@ -241,6 +260,7 @@ typedef struct double vy; /* y velocity of particle */ double omega; /* angular velocity of particle */ double mass_inv; /* inverse of particle mass */ + double charge; /* electric charge */ double inertia_moment_inv; /* inverse of moment of inertia */ double fx; /* x component of force on particle */ double fy; /* y component of force on particle */ @@ -313,6 +333,8 @@ typedef struct double angle01, angle02; /* initial values of angles in which concave corners repel */ double fx, fy; /* x and y-components of force on segment */ double torque; /* torque on segment with respect to its center */ + double pressure; /* pressure acting on segement */ + double avrg_pressure; /* time-averaged pressure */ short int inactivate; /* set to 1 for segment to become inactive at time SEGMENT_DEACTIVATION_TIME */ } t_segment; @@ -361,6 +383,7 @@ typedef struct double angle; /* orientation of obstacle */ double omega; /* angular speed of obstacle */ double bdry_fx, bdry_fy; /* components of boundary force */ + double efield, bfield; /* electric and magnetic field */ } t_lj_parameters; int ncircles, nobstacles, nsegments, ngroups = 1, counter = 0; diff --git a/global_pdes.c b/global_pdes.c index 4df65a6..3045d69 100644 --- a/global_pdes.c +++ b/global_pdes.c @@ -77,6 +77,7 @@ #define D_ONE_FUNNEL 581 /* one funnel */ #define D_LENSES_RING 59 /* several lenses forming a ring */ #define D_MAZE_CIRCULAR 60 /* circular maze */ +#define D_LENS 61 /* symmetric lens made of circular faces */ #define D_WING 70 /* complement of wing-shaped domain */ #define D_TESLA 71 /* Tesla valve */ @@ -89,6 +90,9 @@ #define D_SPHERE_JULIA 82 /* Julia set on Riemann sphere */ #define D_SPHERE_JULIA_INV 83 /* inverted Julia set on Riemann sphere */ #define D_SPHERE_EARTH 84 /* map of the Earth */ +#define D_SPHERE_JULIA_CUBIC 85 /* Julia set for cubic polynomial on Riemann sphere */ +#define D_SPHERE_MARS 86 /* map of Mars */ +#define D_SPHERE_MOON 87 /* map of the Moon */ #define NMAXCIRCLES 10000 /* total number of circles/polygons (must be at least NCX*NCY for square grid) */ #define NMAXPOLY 50000 /* maximal number of vertices of polygonal lines (for von Koch et al) */ @@ -149,6 +153,8 @@ #define IOR_POISSON_WELLS 8 /* wells located on a random Poisson disc process */ #define IOR_PPP_WELLS 9 /* wells located on a Poisson point process */ +#define IOR_EARTH_DEM 10 /* digital elevation model (for waves on sphere) */ + /* Boundary conditions */ #define BC_DIRICHLET 0 /* Dirichlet boundary conditions */ @@ -167,6 +173,18 @@ #define OSC_WAVE_PACKET 2 /* Gaussian wave packet */ #define OSC_CHIRP 3 /* chirp (linearly accelerating frequency) */ +/* Wave packet types */ + +#define WP_RANDOM1 0 /* random, variant 1 */ +#define WP_RANDOM2 1 /* random, variant 2 */ +#define WP_PAIR 2 /* 2 sources */ + +/* Wave packet envelope types */ + +#define WE_SINE 0 /* sine envelope */ +#define WE_CUTOFF 1 /* sharp cut-off */ +#define WE_CONSTANT 2 /* constant envelope */ + /* For debugging purposes only */ // #define FLOOR 0 /* set to 1 to limit wave amplitude to VMAX */ // #define VMAX 10.0 /* max value of wave amplitude */ diff --git a/heat.c b/heat.c index 55eb8f6..99b9b0c 100644 --- a/heat.c +++ b/heat.c @@ -212,6 +212,7 @@ #define WAVE_PACKET_SOURCE_TYPE 1 /* type of wave packet sources */ #define N_WAVE_PACKETS 15 /* number of wave packets */ #define OSCIL_LEFT_YSHIFT 0.0 /* y-dependence of left oscillation (for non-horizontal waves) */ +#define OSCILLATING_SOURCE_PERIOD 20 /* period of oscillating source */ /* end of constants only used by sub_wave and sub_maze */ #include "global_pdes.c" diff --git a/lennardjones.c b/lennardjones.c index 75014c3..721c688 100644 --- a/lennardjones.c +++ b/lennardjones.c @@ -22,7 +22,7 @@ /* It may be possible to increase parameter PAUSE */ /* */ /* create movie using */ -/* ffmpeg -i lj.%05d.tif -vcodec libx264 lj.mp4 */ +/* ffmpeg -i lj.%05d.tif -vcodec libx264 lj.mp4 */ /* */ /*********************************************************************************/ @@ -50,8 +50,8 @@ /* General geometrical parameters */ -#define WINWIDTH 1280 /* window width */ -#define WINHEIGHT 720 /* window height */ +#define WINWIDTH 1600 /* window width */ +#define WINHEIGHT 900 /* window height */ #define XMIN -2.0 #define XMAX 2.0 /* x interval */ @@ -59,7 +59,7 @@ #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ #define INITXMIN -2.0 -#define INITXMAX 2.0 /* x interval for initial condition */ +#define INITXMAX 2.03 /* x interval for initial condition */ #define INITYMIN -1.1 #define INITYMAX 1.1 /* y interval for initial condition */ @@ -86,14 +86,15 @@ #define OBSTACLE_PATTERN 4 /* pattern of obstacles, see list in global_ljones.c */ #define ADD_FIXED_SEGMENTS 1 /* set to 1 to add fixed segments as obstacles */ -#define SEGMENT_PATTERN 24 /* pattern of repelling segments, see list in global_ljones.c */ +#define SEGMENT_PATTERN 25 /* pattern of repelling segments, see list in global_ljones.c */ #define ROCKET_SHAPE 3 /* shape of rocket combustion chamber, see list in global_ljones.c */ #define ROCKET_SHAPE_B 3 /* shape of second rocket */ #define NOZZLE_SHAPE 6 /* shape of nozzle, see list in global_ljones.c */ #define NOZZLE_SHAPE_B 6 /* shape of nozzle for second rocket, see list in global_ljones.c */ -#define TWO_TYPES 0 /* set to 1 to have two types of particles */ +#define TWO_TYPES 1 /* set to 1 to have two types of particles */ #define TYPE_PROPORTION 0.5 /* proportion of particles of first type */ +#define TWOTYPE_CONFIG 2 /* choice of types, see TTC_ list in global_ljones.c */ #define SYMMETRIZE_FORCE 1 /* set to 1 to symmetrize two-particle interaction, only needed if particles are not all the same */ #define CENTER_PX 0 /* set to 1 to center horizontal momentum */ #define CENTER_PY 0 /* set to 1 to center vertical momentum */ @@ -110,9 +111,9 @@ #define PDISC_CANDIDATES 100 /* number of candidates in construction of Poisson disc process */ #define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */ -#define LAMBDA 0.75 /* parameter controlling the dimensions of domain */ -#define MU 0.012 /* parameter controlling radius of particles */ -#define MU_B 0.010 /* parameter controlling radius of particles of second type */ +#define LAMBDA 0.3 /* parameter controlling the dimensions of domain */ +#define MU 0.014 /* parameter controlling radius of particles */ +#define MU_B 0.01 /* parameter controlling radius of particles of second type */ #define NPOLY 40 /* number of sides of polygon */ #define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */ #define AWEDGE 0.5 /* opening angle of wedge, in units of Pi/2 */ @@ -121,10 +122,8 @@ #define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ #define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */ #define FOCI 1 /* set to 1 to draw focal points of ellipse */ -// #define NGRIDX 140 /* number of grid point for grid of disks */ -// #define NGRIDY 70 /* number of grid point for grid of disks */ -#define NGRIDX 130 /* number of grid point for grid of disks */ -#define NGRIDY 65 /* number of grid point for grid of disks */ +#define NGRIDX 101 /* number of grid point for grid of disks */ +#define NGRIDY 47 /* number of grid point for grid of disks */ #define EHRENFEST_RADIUS 0.9 /* radius of container for Ehrenfest urn configuration */ #define EHRENFEST_WIDTH 0.035 /* width of tube for Ehrenfest urn configuration */ #define TWO_CIRCLES_RADIUS_RATIO 0.8 /* ratio of radii for S_TWO_CIRCLES_EXT segment configuration */ @@ -137,9 +136,8 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 4525 /* number of frames of movie */ -// #define NSTEPS 1000 /* number of frames of movie */ -#define NVID 30 /* number of iterations between images displayed on screen */ +#define NSTEPS 3200 /* number of frames of movie */ +#define NVID 15 /* number of iterations between images displayed on screen */ #define NSEG 25 /* number of segments of boundary of circles */ #define INITIAL_TIME 0 /* time after which to start saving frames */ #define OBSTACLE_INITIAL_TIME 150 /* time after which to start moving obstacle */ @@ -152,7 +150,7 @@ #define SLEEP1 1 /* initial sleeping time */ #define SLEEP2 1 /* final sleeping time */ #define MID_FRAMES 20 /* number of still frames between parts of two-part movie */ -#define END_FRAMES 100 /* number of still frames at end of movie */ +#define END_FRAMES 250 /* number of still frames at end of movie */ /* Boundary conditions, see list in global_ljones.c */ @@ -160,18 +158,18 @@ /* Plot type, see list in global_ljones.c */ -#define PLOT 13 +#define PLOT 5 #define PLOT_B 11 /* plot type for second movie */ #define DRAW_BONDS 1 /* set to 1 to draw bonds between neighbours */ #define COLOR_BONDS 1 /* set to 1 to color bonds according to length */ #define FILL_TRIANGLES 0 /* set to 1 to fill triangles between neighbours */ #define ALTITUDE_LINES 0 /* set to 1 to add horizontal lines to show altitude */ -#define COLOR_SEG_GROUPS 1 /* set to 1 to collor segment groups differently */ +#define COLOR_SEG_GROUPS 0 /* set to 1 to collor segment groups differently */ #define N_PARTICLE_COLORS 200 /* number of colors for P_NUMBER color scheme */ #define INITIAL_POS_TYPE 0 /* type of initial position dependence */ -#define ERATIO 0.995 /* ratio for time-averagin in P_EMEAN color scheme */ -#define DRATIO 0.995 /* ratio for time-averagin in P_DIRECT_EMEAN color scheme */ +#define ERATIO 0.995 /* ratio for time-averaging in P_EMEAN color scheme */ +#define DRATIO 0.995 /* ratio for time-averaging in P_DIRECT_EMEAN color scheme */ /* Color schemes */ @@ -195,15 +193,15 @@ #define LUMAMP 0.3 /* amplitude of luminosity variation for scheme C_LUM */ #define HUEMEAN 220.0 /* mean value of hue for color scheme C_HUE */ #define HUEAMP -50.0 /* amplitude of variation of hue for color scheme C_HUE */ -#define COLOR_HUESHIFT 0.5 /* shift in color hue (for some cyclic palettes) */ +#define COLOR_HUESHIFT 0.0 /* shift in color hue (for some cyclic palettes) */ -#define PRINT_PARAMETERS 1 /* set to 1 to print certain parameters */ +#define PRINT_PARAMETERS 0 /* set to 1 to print certain parameters */ #define PRINT_TEMPERATURE 0 /* set to 1 to print current temperature */ -#define PRINT_ANGLE 1 /* set to 1 to print obstacle orientation */ -#define PRINT_OMEGA 1 /* set to 1 to print angular speed */ +#define PRINT_ANGLE 0 /* set to 1 to print obstacle orientation */ +#define PRINT_OMEGA 0 /* set to 1 to print angular speed */ #define PRINT_PARTICLE_SPEEDS 0 /* set to 1 to print average speeds/momenta of particles */ #define PRINT_SEGMENTS_SPEEDS 0 /* set to 1 to print velocity of moving segments */ -#define PRINT_SEGMENTS_FORCE 1 /* set to 1 to print force on segments */ +#define PRINT_SEGMENTS_FORCE 0 /* set to 1 to print force on segments */ #define FORCE_FACTOR 0.1 /* factor controlling length of force vector */ /* particle properties */ @@ -212,7 +210,7 @@ #define ENERGY_HUE_MAX 50.0 /* color of saturated particle */ #define PARTICLE_HUE_MIN 359.0 /* color of original particle */ #define PARTICLE_HUE_MAX 0.0 /* color of saturated particle */ -#define PARTICLE_EMAX 5000.0 /* energy of particle with hottest color */ +#define PARTICLE_EMAX 1000.0 /* energy of particle with hottest color */ #define HUE_TYPE0 280.0 /* hue of particles of type 0 */ #define HUE_TYPE1 60.0 /* hue of particles of type 1 */ #define HUE_TYPE2 140.0 /* hue of particles of type 2 */ @@ -221,35 +219,35 @@ #define RANDOM_RADIUS 0 /* set to 1 for random circle radius */ #define DT_PARTICLE 3.0e-6 /* time step for particle displacement */ #define KREPEL 1.0 /* constant in repelling force between particles */ -#define EQUILIBRIUM_DIST 5.0 /* Lennard-Jones equilibrium distance */ -#define EQUILIBRIUM_DIST_B 5.0 /* Lennard-Jones equilibrium distance for second type of particle */ +#define EQUILIBRIUM_DIST 7.0 /* Lennard-Jones equilibrium distance */ +#define EQUILIBRIUM_DIST_B 6.5 /* Lennard-Jones equilibrium distance for second type of particle */ #define REPEL_RADIUS 15.0 /* radius in which repelling force acts (in units of particle radius) */ -#define DAMPING 100.0 /* damping coefficient of particles */ +#define DAMPING 0.02 /* damping coefficient of particles */ #define OMEGA_INITIAL 10.0 /* initial angular velocity range */ #define INITIAL_DAMPING 5.0 /* damping coefficient of particles during initial phase */ #define DAMPING_ROT 100.0 /* dampint coefficient for rotation of particles */ #define PARTICLE_MASS 1.0 /* mass of particle of radius MU */ -#define PARTICLE_MASS_B 1.0 /* mass of particle of radius MU */ +#define PARTICLE_MASS_B 6.5 /* mass of particle of radius MU */ #define PARTICLE_INERTIA_MOMENT 0.2 /* moment of inertia of particle */ #define PARTICLE_INERTIA_MOMENT_B 0.02 /* moment of inertia of second type of particle */ -#define V_INITIAL 0.0 /* initial velocity range */ +#define V_INITIAL 1000.0 /* initial velocity range */ #define OMEGA_INITIAL 10.0 /* initial angular velocity range */ #define VICSEK_VMIN 1.0 /* minimal speed of particles in Vicsek model */ #define VICSEK_VMAX 40.0 /* minimal speed of particles in Vicsek model */ +#define V_INITIAL_TYPE 1 /* type of initial speed distribution (see VI_ in global_ljones.c) */ + #define THERMOSTAT 0 /* set to 1 to switch on thermostat */ #define VARY_THERMOSTAT 0 /* set to 1 for time-dependent thermostat schedule */ #define SIGMA 5.0 /* noise intensity in thermostat */ -#define BETA 0.02 /* initial inverse temperature */ +#define BETA 0.0002 /* initial inverse temperature */ #define MU_XI 0.01 /* friction constant in thermostat */ #define KSPRING_BOUNDARY 1.0e7 /* confining harmonic potential outside simulation region */ #define KSPRING_OBSTACLE 1.0e11 /* harmonic potential of obstacles */ -// #define NBH_DIST_FACTOR 2.3 /* radius in which to count neighbours */ -// #define NBH_DIST_FACTOR 2.7 /* radius in which to count neighbours */ #define NBH_DIST_FACTOR 3.8 /* radius in which to count neighbours */ -// #define NBH_DIST_FACTOR 4.0 /* radius in which to count neighbours */ -#define GRAVITY 0.0 /* gravity acting on all particles */ -#define GRAVITY_X 5000.0 /* horizontal gravity acting on all particles */ +#define GRAVITY 0.0 /* gravity acting on all particles */ +#define GRAVITY_X 0.0 /* horizontal gravity acting on all particles */ +#define CIRCULAR_GRAVITY 0 /* set to 1 to have gravity directed to center */ #define INCREASE_GRAVITY 0 /* set to 1 to increase gravity during the simulation */ #define GRAVITY_SCHEDULE 2 /* type of gravity schedule, see list in global_ljones.c */ #define GRAVITY_FACTOR 100.0 /* factor by which to increase gravity */ @@ -258,6 +256,15 @@ #define KSPRING_VICSEK 0.2 /* spring constant for I_VICSEK_SPEED interaction */ #define VICSEK_REPULSION 10.0 /* repulsion between particles in Vicsek model */ +#define ADD_EFIELD 1 /* set to 1 to add an electric field */ +#define EFIELD 10000.0 /* value of electric field */ +#define ADD_BFIELD 0 /* set to 1 to add a magnetic field */ +#define BFIELD 1000.0 /* value of magnetic field */ +#define CHARGE 0.0 /* charge of particles of first type */ +#define CHARGE_B 1.0 /* charge of particles of second type */ +#define INCREASE_E 0 /* set to 1 to increase electric field */ +#define EFIELD_FACTOR 1000.0 /* factor by which to increase electric field */ + #define ROTATION 0 /* set to 1 to include rotation of particles */ #define COUPLE_ANGLE_TO_THERMOSTAT 1 /* set to 1 to couple angular degrees of freedom to thermostat */ #define DIMENSION_FACTOR 1.0 /* scaling factor taking into account number of degrees of freedom */ @@ -267,7 +274,7 @@ #define KTORQUE_DIFF 150.0 /* force constant in angular dynamics for different particles */ #define DRAW_SPIN 0 /* set to 1 to draw spin vectors of particles */ #define DRAW_SPIN_B 0 /* set to 1 to draw spin vectors of particles */ -#define DRAW_CROSS 1 /* set to 1 to draw cross on particles of second type */ +#define DRAW_CROSS 0 /* set to 1 to draw cross on particles of second type */ #define SPIN_RANGE 10.0 /* range of spin-spin interaction */ #define SPIN_RANGE_B 5.0 /* range of spin-spin interaction for second type of particle */ #define QUADRUPOLE_RATIO 0.6 /* anisotropy in quadrupole potential */ @@ -290,7 +297,7 @@ #define CENTER_VIEW_ON_OBSTACLE 0 /* set to 1 to center display on moving obstacle */ #define RESAMPLE_Y 0 /* set to 1 to resample y coordinate of moved particles (for shock waves) */ #define NTRIALS 2000 /* number of trials when resampling */ -#define OBSTACLE_RADIUS 0.25 /* radius of obstacle for circle boundary conditions */ +#define OBSTACLE_RADIUS 0.15 /* radius of obstacle for circle boundary conditions */ #define FUNNEL_WIDTH 0.25 /* funnel width for funnel boundary conditions */ #define OBSTACLE_XMIN 0.0 /* initial position of obstacle */ #define OBSTACLE_XMAX 3.0 /* final position of obstacle */ @@ -325,24 +332,28 @@ #define TRACER_PARTICLE_MASS 4.0 /* relative mass of tracer particle */ #define TRAJECTORY_WIDTH 3 /* width of tracer particle trajectory */ -#define ROTATE_BOUNDARY 1 /* set to 1 to rotate the repelling segments */ +#define ROTATE_BOUNDARY 0 /* set to 1 to rotate the repelling segments */ #define SMOOTH_ROTATION 1 /* set to 1 to update segments at each time step (rather than at each movie frame) */ +#define ROTATION_SCHEDULE 0 /* time-dependence of rotation angle, see ROT_* in global_ljones.c */ #define PERIOD_ROTATE_BOUNDARY 1000 /* period of rotating boundary */ #define ROTATE_INITIAL_TIME 300 /* initial time without rotation */ #define ROTATE_FINAL_TIME 300 /* final time without rotation */ #define ROTATE_CHANGE_TIME 0.5 /* relative duration of acceleration/deceleration phases */ -#define OMEGAMAX -2.0*PI /* maximal rotation speed */ +#define OMEGAMAX -2.0*PI /* maximal rotation speed */ #define MOVE_BOUNDARY 0 /* set to 1 to move repelling segments, due to force from particles */ #define SEGMENTS_MASS 40.0 /* mass of collection of segments */ #define DEACTIVATE_SEGMENT 1 /* set to 1 to deactivate last segment after a certain time */ #define SEGMENT_DEACTIVATION_TIME 20 /* time at which to deactivate last segment */ #define RELEASE_ROCKET_AT_DEACTIVATION 0 /* set to 1 to limit segments velocity before segment release */ -#define SEGMENTS_X0 1.5 /* initial position of segments */ -#define SEGMENTS_Y0 0.0 /* initial position of segments */ +#define SEGMENTS_X0 0.0 /* initial position of segments */ +#define SEGMENTS_Y0 -0.6 /* initial position of segments */ #define SEGMENTS_VX0 0.0 /* initial velocity of segments */ #define SEGMENTS_VY0 0.0 /* initial velocity of segments */ #define DAMP_SEGS_AT_NEGATIVE_Y 0 /* set to 1 to dampen segments when y coordinate is negative */ +#define SHOW_SEGMENTS_PRESSURE 0 /* set to 1 to show (averaged) pressure acting on segments */ +#define SEGMENT_PMAX 7.5e7 /* pressure of segment with hottest color */ +#define P_AVRG_FACTOR 0.02 /* factor in computation of mean pressure */ #define MOVE_SEGMENT_GROUPS 0 /* set to 1 to group segments into moving units */ #define SEGMENT_GROUP_MASS 500.0 /* mass of segment group */ @@ -407,8 +418,8 @@ #define FLOOR_OMEGA 0 /* set to 1 to limit particle momentum to PMAX */ #define PMAX 1000.0 /* maximal force */ -#define HASHX 120 /* size of hashgrid in x direction */ -#define HASHY 60 /* size of hashgrid in y direction */ +#define HASHX 60 /* size of hashgrid in x direction */ +#define HASHY 30 /* size of hashgrid in y direction */ #define HASHMAX 100 /* maximal number of particles per hashgrid cell */ #define HASHGRID_PADDING 0.1 /* padding of hashgrid outside simulation window */ @@ -473,6 +484,23 @@ double repel_schedule(int i) } +double efield_schedule(int i) +{ + static double efactor; + static int first = 1; + double efield; + + if (first) + { + efactor = EFIELD_FACTOR/(double)(INITIAL_TIME + NSTEPS); + first = 0; + } + efield = EFIELD*(double)i*efactor; + printf("E = %.3lg\n", efield); + return(efield); +} + + double temperature_schedule(int i) { static double bexponent, omega, bexp2; @@ -614,20 +642,26 @@ double rotation_angle(double phase) // return(0.5*OMEGAMAX*(phase - sin(phase))); /* case of centrifuge remaining at constant speed for a while */ - if (phase < ROTATE_CHANGE_TIME) - { -// angular_speed = 0.5*OMEGAMAX*(1.0 - cos(phase*PI/ROTATE_CHANGE_TIME)); - return(0.5*OMEGAMAX*(phase - (ROTATE_CHANGE_TIME/PI)*sin(phase*PI/ROTATE_CHANGE_TIME))); - } - else if (phase < 1.0 - ROTATE_CHANGE_TIME) - { -// angular_speed = OMEGAMAX; - return(0.5*OMEGAMAX*(2.0*phase - ROTATE_CHANGE_TIME)); - } - else - { -// angular_speed = 0.5*OMEGAMAX*(1.0 + cos((phase - 1.0 + ROTATE_CHANGE_TIME)*PI/ROTATE_CHANGE_TIME)); - return(0.5*OMEGAMAX*(2.0 - 2.0*ROTATE_CHANGE_TIME + phase - 1.0 + (ROTATE_CHANGE_TIME/PI)*sin((1.0-phase)*PI/ROTATE_CHANGE_TIME))); + switch (ROTATION_SCHEDULE) { + case (ROT_SPEEDUP_SLOWDOWN): + { + if (phase < ROTATE_CHANGE_TIME) + { + return(0.5*OMEGAMAX*(phase - (ROTATE_CHANGE_TIME/PI)*sin(phase*PI/ROTATE_CHANGE_TIME))); + } + else if (phase < 1.0 - ROTATE_CHANGE_TIME) + { + return(0.5*OMEGAMAX*(2.0*phase - ROTATE_CHANGE_TIME)); + } + else + { + return(0.5*OMEGAMAX*(2.0 - 2.0*ROTATE_CHANGE_TIME + phase - 1.0 + (ROTATE_CHANGE_TIME/PI)*sin((1.0-phase)*PI/ROTATE_CHANGE_TIME))); + } + } + case (ROT_BACK_FORTH): + { + return(OMEGAMAX*(1.0 - cos(DPI*phase))/DPI); + } } } @@ -1252,7 +1286,8 @@ void animation() printf("Computing frame %d\n",i); if (INCREASE_KREPEL) params.krepel = repel_schedule(i); - if (INCREASE_BETA) params.beta = temperature_schedule(i); + if (INCREASE_BETA) params.beta = temperature_schedule(i); + if (INCREASE_E) params.efield = efield_schedule(i); if (DECREASE_CONTAINER_SIZE) { params.xmincontainer = container_size_schedule(i); @@ -1358,12 +1393,33 @@ void animation() /* add gravity */ if (INCREASE_GRAVITY) particle[j].fy -= params.gravity/particle[j].mass_inv; + else if (CIRCULAR_GRAVITY) + { + particle[j].fx -= GRAVITY*particle[j].xc/particle[j].mass_inv; + particle[j].fy -= GRAVITY*particle[j].yc/particle[j].mass_inv; + } else { particle[j].fy -= GRAVITY/particle[j].mass_inv; particle[j].fx += GRAVITY_X/particle[j].mass_inv; } + /* add electric force */ + if (ADD_EFIELD) + { + if (INCREASE_E) particle[j].fx += params.efield*particle[j].charge; + else particle[j].fx += EFIELD*particle[j].charge; + } + + /* add magnetic force */ + if (ADD_BFIELD) + { +// particle[j].fx += BFIELD*particle[j].charge*particle[j].vy; +// particle[j].fy -= BFIELD*particle[j].charge*particle[j].vx; + particle[j].fx += BFIELD*particle[j].charge*particle[j].vy*particle[j].mass_inv; + particle[j].fy -= BFIELD*particle[j].charge*particle[j].vx*particle[j].mass_inv; + } + if (FLOOR_FORCE) { if (particle[j].fx > FMAX) particle[j].fx = FMAX; diff --git a/makefile b/makefile index 54fa8ba..f4f9a32 100644 --- a/makefile +++ b/makefile @@ -1,4 +1,4 @@ -LIBS = -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lXmu -lglut -fopenmp +LIBS = -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lglut -fopenmp .c: gcc -o $@ $< -O3 $(LIBS) diff --git a/mangrove.c b/mangrove.c index d901517..1253439 100644 --- a/mangrove.c +++ b/mangrove.c @@ -251,16 +251,11 @@ #define IOR 7 /* choice of index of refraction, see list in global_pdes.c */ #define IOR_TOTAL_TURNS 1.5 /* total angle of rotation for IOR_PERIODIC_WELLS_ROTATING */ #define MANDEL_IOR_SCALE -0.05 /* parameter controlling dependence of IoR on Mandelbrot escape speed */ -#define COURANT 0.04 /* Courant number */ -#define COURANTB 0.0 /* Courant number in medium B */ -#define INITIAL_AMP 0.5 /* amplitude of initial condition */ -#define INITIAL_VARIANCE 0.0003 /* variance of initial condition */ -#define INITIAL_WAVELENGTH 0.015 /* wavelength of initial condition */ -#define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */ #define WAVE_PACKET_SOURCE_TYPE 1 /* type of wave packet sources */ #define N_WAVE_PACKETS 15 /* number of wave packets */ #define OSCIL_LEFT_YSHIFT 0.0 /* y-dependence of left oscillation (for non-horizontal waves) */ #define DRAW_WAVE_PROFILE 0 /* set to 1 to draw a profile of the wave */ +#define OSCILLATING_SOURCE_PERIOD 20 /* period of oscillating source */ /* end of constants only used by sub_wave and sub_maze */ #include "global_pdes.c" diff --git a/marscyl2.ppm.gz b/marscyl2.ppm.gz new file mode 100644 index 0000000..faa9fab Binary files /dev/null and b/marscyl2.ppm.gz differ diff --git a/schrodinger.c b/schrodinger.c index fef5209..9b8bc45 100644 --- a/schrodinger.c +++ b/schrodinger.c @@ -192,6 +192,7 @@ #define WAVE_PACKET_SOURCE_TYPE 1 /* type of wave packet sources */ #define N_WAVE_PACKETS 15 /* number of wave packets */ #define OSCIL_LEFT_YSHIFT 0.0 /* y-dependence of left oscillation (for non-horizontal waves) */ +#define OSCILLATING_SOURCE_PERIOD 20 /* period of oscillating source */ /* end of constants only used by sub_wave and sub_maze */ diff --git a/sub_lj.c b/sub_lj.c index 28c988e..42d559d 100644 --- a/sub_lj.c +++ b/sub_lj.c @@ -34,8 +34,8 @@ int writetiff(char *filename, char *description, int x, int y, int width, int he glPixelStorei(GL_PACK_ALIGNMENT, 1); glReadPixels(x, y, width, height, GL_RGB, GL_UNSIGNED_BYTE, image); - TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32) width); - TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32) height); + TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32_t) width); + TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32_t) height); TIFFSetField(file, TIFFTAG_BITSPERSAMPLE, 8); TIFFSetField(file, TIFFTAG_COMPRESSION, compression); TIFFSetField(file, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB); @@ -622,15 +622,18 @@ void init_particle_config(t_particle particles[NMAXCIRCLES]) case (C_SQUARE): { ncircles = NGRIDX*NGRIDY; - dy = (YMAX - YMIN)/((double)NGRIDY); + dx = (INITXMAX - INITXMIN)/((double)NGRIDX); + dy = (INITYMAX - INITYMIN)/((double)NGRIDY); for (i = 0; i < NGRIDX; i++) for (j = 0; j < NGRIDY; j++) { n = NGRIDY*i + j; - particles[n].xc = ((double)(i-NGRIDX/2) + 0.5)*dy; - particles[n].yc = YMIN + ((double)j + 0.5)*dy; + particles[n].xc = INITXMIN + ((double)i - 0.5)*dx; + particles[n].yc = INITYMIN + ((double)j - 0.5)*dy; particles[n].radius = MU; - particles[n].active = 1; + /* activate only circles that intersect the domain */ + if ((particles[n].yc < INITYMAX + MU)&&(particles[n].yc > INITYMIN - MU)&&(particles[n].xc < INITXMAX + MU)&&(particles[n].xc > INITXMIN - MU)) particles[n].active = 1; + else particles[n].active = 0; } break; } @@ -1616,7 +1619,7 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS]) /* initialise linear obstacle configuration */ { int i, j, cycle = 0, iminus, iplus, nsides, n, concave = 1; - double angle, angle2, dangle, dx, width, height, a, b, length, xmid = 0.5*(BCXMIN + BCXMAX), lpocket, r, x, x1, y1, x2, y2, nozx, nozy, y, dy, ca, sa; + double angle, angle2, dangle, dx, width, height, a, b, length, xmid = 0.5*(BCXMIN + BCXMAX), lpocket, r, x, x1, y1, x2, y2, nozx, nozy, y, yy, dy, ca, sa, padding; switch (SEGMENT_PATTERN) { case (S_RECTANGLE): @@ -2399,6 +2402,82 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS]) break; } + case (S_COANDA): + { + y1 = SEGMENTS_Y0; + width = 0.05; + padding = 0.01; + + height = 0.25; + dx = (XMAX - XMIN + 2.0*padding)/(double)(NPOLY-1); + + for (i=0; i= BCXMAX) return(0); if (x <= BCXMIN) return(0); @@ -2750,6 +2833,17 @@ int in_segment_region(double x, double y, t_segment segment[NMAXSEGMENTS]) y1 += 0.5*x1*x1; return(x1*x1 + 25.0*y1*y1 > LAMBDA*LAMBDA*(1.0 + padding)); } + case (S_COANDA): + { + return(vabs(y - SEGMENTS_Y0 - 0.25*cos(PI*x/XMAX)) > 0.1); + } + case (S_COANDA_SHORT): + { + x1 = XMIN + 0.1; + length = 0.85*(XMAX - XMIN - 0.2); + if (x > x1 + length + 0.1) return(1); + return(vabs(y - SEGMENTS_Y0 + 0.2*cos(DPI*(x-x1)/length)) > 0.05); + } case (S_MAZE): { for (i=0; i DPI) hue -= DPI; hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*(hue)/DPI; ej = particle.emean; - if (ej < 0.1*PARTICLE_EMAX) lum = 10.0*ej/PARTICLE_EMAX; + if (ej < 0.5*PARTICLE_EMAX) lum = 2.0*ej/PARTICLE_EMAX; else lum = 1.0; *radius = particle.radius; *width = BOUNDARY_WIDTH; @@ -4127,6 +4232,21 @@ void set_segment_group_color(int group, double lum, double rgb[3]) glColor3f(lum*rgb[0], lum*rgb[1], lum*rgb[2]); } +void set_segment_pressure_color(double pressure, double lum, double rgb[3]) +{ + double hue; + +// if (pressure < 0.0) pressure = 0.0; + hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*pressure/SEGMENT_PMAX; + if (hue > PARTICLE_HUE_MIN) hue = PARTICLE_HUE_MIN; + else if (hue < PARTICLE_HUE_MAX) hue = PARTICLE_HUE_MAX; + + hsl_to_rgb_turbo(hue, 0.9, lum, rgb); + + glColor3f(rgb[0], rgb[1], rgb[2]); +} + + void draw_altitude_lines() { int i, i1, i2; @@ -4625,8 +4745,18 @@ int draw_special_particle(t_particle particle, double xc1, double yc1, double ra void draw_one_particle(t_particle particle, double xc, double yc, double radius, double angle, int nsides, double width, double rgb[3]) /* draw one of the particles */ { - double ca, sa, x1, x2, y1, y2, xc1, yc1, wangle, newradius = radius; + double x1, x2, y1, y2, xc1, yc1, wangle, newradius = radius; int wsign, cont = 1, draw = 1; + static double ca, sa; + static int first = 1; + + if (first) + { + angle = APOLY*PID; + ca = cos(angle); + sa = sin(angle); + first = 0; + } if (CENTER_VIEW_ON_OBSTACLE) xc1 = xc - xshift; else xc1 = xc; @@ -4651,14 +4781,16 @@ void draw_one_particle(t_particle particle, double xc, double yc, double radius, /* draw crosses on particles of second type */ if ((TWO_TYPES)&&(DRAW_CROSS)) - if (particle.type == 1) + if (particle.type == 2) { - if (ROTATION) angle = angle + APOLY*PID; - else angle = APOLY*PID; - ca = cos(angle); - sa = sin(angle); - glLineWidth(3); - glColor3f(0.0, 0.0, 0.0); + if (ROTATION) + { + angle = angle + APOLY*PID; + ca = cos(angle); + sa = sin(angle); + } + glLineWidth(2); + glColor3f(1.0, 1.0, 1.0); x1 = xc1 - MU_B*ca; y1 = yc1 - MU_B*sa; x2 = xc1 + MU_B*ca; @@ -5002,6 +5134,7 @@ void draw_container(double xmin, double xmax, t_obstacle obstacle[NMAXOBSTACLES] for (i = 0; i < nsegments; i++) if (segment[i].active) { if (COLOR_SEG_GROUPS) set_segment_group_color(segment[i].group, 1.0, rgb); + else if (SHOW_SEGMENTS_PRESSURE) set_segment_pressure_color(segment[i].avrg_pressure, 1.0, rgb); draw_line(segment[i].x1 - xtrack, segment[i].y1 - ytrack, segment[i].x2 - xtrack, segment[i].y2 - ytrack); } } @@ -5248,6 +5381,12 @@ void print_omega(double angle, double angular_speed, double fx, double fy) glColor3f(1.0, 1.0, 1.0); sprintf(message, "Fy = %.4f", fy); write_text(xrighttext + 0.1, y1, message); + + y1 -= 0.1; + erase_area_hsl(xrightbox, y1 + 0.025, 0.42, 0.05, 0.0, 0.9, 0.0); + glColor3f(1.0, 1.0, 1.0); + sprintf(message, "Fy/Fx = %.4f", fy/fx); + write_text(xrighttext + 0.1, y1, message); } } @@ -5389,10 +5528,17 @@ void print_parameters(t_lj_parameters params, short int left, double pressure[N_ } else if (INCREASE_KREPEL) /* print force constant */ { - erase_area_hsl(xbox, y + 0.025*scale, 0.22*scale, 0.05*scale, 0.0, 0.9, 0.0); + erase_area_hsl(xbox, y + 0.025*scale, 0.35*scale, 0.05*scale, 0.0, 0.9, 0.0); glColor3f(1.0, 1.0, 1.0); - sprintf(message, "Force %.0f", params.krepel); - write_text(xtext + 0.28, y, message); + sprintf(message, "Force constant %.0f", params.krepel); + write_text(xtext + 0.03, y, message); + } + else if (INCREASE_E) /* print electric field */ + { + erase_area_hsl(xbox, y + 0.025*scale, 0.27*scale, 0.05*scale, 0.0, 0.9, 0.0); + glColor3f(1.0, 1.0, 1.0); + sprintf(message, "E field = %.2f", 25.0*NVID*DT_PARTICLE*params.efield); + write_text(xtext + 0.08, y, message); } if (RECORD_PRESSURES) @@ -5869,9 +6015,7 @@ void print_particles_speeds(t_particle particle[NMAXCIRCLES]) write_text(xrighttext + 0.1, y, message); } -double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NMAXOBSTACLES], - t_segment segment[NMAXSEGMENTS], - double xleft, double xright, double *pleft, double *pright, double pressure[N_PRESSURES], int wall) +double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NMAXOBSTACLES], t_segment segment[NMAXSEGMENTS], double xleft, double xright, double *pleft, double *pright, double pressure[N_PRESSURES], int wall) { int i, k; double xmin, xmax, ymin, ymax, padding, r, rp, r2, cphi, sphi, f, fperp = 0.0, x, y, xtube, distance, dx, dy, width, ybin, angle, x1, x2, h, ytop, norm, dleft, dplus, dminus, tmp_pleft = 0.0, tmp_pright = 0.0, proj, pscal, pvect, pvmin; @@ -5914,8 +6058,6 @@ double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacl particle[j].fx += f*segment[i].nx; particle[j].fy += f*segment[i].ny; - - if ((MOVE_BOUNDARY)||(MOVE_SEGMENT_GROUPS)||(PRINT_SEGMENTS_FORCE)) { segment[i].fx -= f*segment[i].nx; @@ -5923,6 +6065,14 @@ double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacl segment[i].torque -= (x - segment[i].xc)*f*segment[i].ny - (y - segment[i].yc)*f*segment[i].nx; // printf("Segment %i: f = (%.3lg, %.3lg)\n", i, segment[i].fx, segment[i].fy); } + + if (SHOW_SEGMENTS_PRESSURE) + { + segment[i].pressure = f/segment[i].length; + segment[i].avrg_pressure *= (1.0 - P_AVRG_FACTOR); + segment[i].avrg_pressure += P_AVRG_FACTOR*segment[i].pressure; +// printf("Segment %i: pressure = %.3lg\n", i, segment[i].avrg_pressure); + } } if ((VICSEK_INT)&&(vabs(distance) < 1.5*r)) { @@ -6480,21 +6630,51 @@ int reorder_particles(t_particle particle[NMAXCIRCLES], double py[NMAXCIRCLES], } +int twotype_config(int i, t_particle particle[NMAXCIRCLES]) +/* assign different particle types */ +{ + switch (TWOTYPE_CONFIG) { + case (TTC_RANDOM): return((double)rand()/RAND_MAX > TYPE_PROPORTION); + case (TTC_CHESSBOARD): + { + switch (CIRCLE_PATTERN) { + case (C_SQUARE): return((i%NGRIDY + i/NGRIDY)%2 == 0); + case (C_HEX): return(i%2 == 0); + default: return(i%2 == 0); + } + } + case (TTC_COANDA): + { + if (vabs(particle[i].yc) < LAMBDA) return(1); + if (vabs(particle[i].yc) < LAMBDA + MU) return(-1); + return(0); + } + default: return(0); + } +} + + int initialize_configuration(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HASHX*HASHY], t_obstacle obstacle[NMAXOBSTACLES], double px[NMAXCIRCLES], double py[NMAXCIRCLES], double pangle[NMAXCIRCLES], int tracer_n[N_TRACER_PARTICLES], t_segment segment[NMAXSEGMENTS]) /* initialize all particles, obstacles, and the hashgrid */ { - int i, j, k, n, nactive = 0; + int i, j, k, n, type, nactive = 0; double x, y, h, xx, yy, rnd; for (i=0; i < ncircles; i++) { /* set particle type */ particle[i].type = 0; - if ((TWO_TYPES)&&((double)rand()/RAND_MAX > TYPE_PROPORTION)) +// if ((TWO_TYPES)&&((double)rand()/RAND_MAX > TYPE_PROPORTION)) + if ((TWO_TYPES)&&(twotype_config(i, particle))) { - particle[i].type = 2; - particle[i].radius = MU_B; + type = twotype_config(i, particle); + if (type == 1) + { + particle[i].type = 2; + particle[i].radius = MU_B; + } + else if (type == -1) particle[i].active = 0; } if ((INTERACTION == I_VICSEK_SHARK)&&(i==1)) { @@ -6524,6 +6704,7 @@ t_segment segment[NMAXSEGMENTS]) particle[i].spin_freq = SPIN_INTER_FREQUENCY; particle[i].mass_inv = 1.0/PARTICLE_MASS; particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT; + particle[i].charge = CHARGE; } else { @@ -6533,10 +6714,25 @@ t_segment segment[NMAXSEGMENTS]) particle[i].spin_freq = SPIN_INTER_FREQUENCY_B; particle[i].mass_inv = 1.0/PARTICLE_MASS_B; particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT_B; + particle[i].charge = CHARGE_B; + } + + switch (V_INITIAL_TYPE) + { + case (VI_RANDOM): + { + particle[i].vx = V_INITIAL*gaussian(); + particle[i].vy = V_INITIAL*gaussian(); + break; + } + case (VI_COANDA): + { + if (vabs(particle[i].yc) < LAMBDA) particle[i].vx = V_INITIAL; + else particle[i].vx = 0.0; + particle[i].vy = 0.0; + break; + } } - - particle[i].vx = V_INITIAL*gaussian(); - particle[i].vy = V_INITIAL*gaussian(); if ((INTERACTION == I_VICSEK_SHARK)&&(i==1)) { @@ -6596,6 +6792,7 @@ t_segment segment[NMAXSEGMENTS]) particle[i].dirmean = 0.0; particle[i].mass_inv = 1.0/PARTICLE_MASS; particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT; + particle[i].charge = 0.0; particle[i].vx = 0.0; particle[i].vy = 0.0; px[i] = 0.0; diff --git a/sub_part_billiard.c b/sub_part_billiard.c index ef8fa4f..02a594f 100644 --- a/sub_part_billiard.c +++ b/sub_part_billiard.c @@ -110,8 +110,8 @@ int writetiff(char *filename, char *description, int x, int y, int width, int he glPixelStorei(GL_PACK_ALIGNMENT, 1); glReadPixels(x, y, width, height, GL_RGB, GL_UNSIGNED_BYTE, image); - TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32) width); - TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32) height); + TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32_t) width); + TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32_t) height); TIFFSetField(file, TIFFTAG_BITSPERSAMPLE, 8); TIFFSetField(file, TIFFTAG_COMPRESSION, compression); TIFFSetField(file, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB); diff --git a/sub_part_billiard_phasespace.c b/sub_part_billiard_phasespace.c index bda0376..0c22f7c 100644 --- a/sub_part_billiard_phasespace.c +++ b/sub_part_billiard_phasespace.c @@ -112,8 +112,8 @@ int writetiff(char *filename, char *description, int x, int y, int width, int he glPixelStorei(GL_PACK_ALIGNMENT, 1); glReadPixels(x, y, width, height, GL_RGB, GL_UNSIGNED_BYTE, image); - TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32) width); - TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32) height); + TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32_t) width); + TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32_t) height); TIFFSetField(file, TIFFTAG_BITSPERSAMPLE, 8); TIFFSetField(file, TIFFTAG_COMPRESSION, compression); TIFFSetField(file, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB); diff --git a/sub_perco.c b/sub_perco.c index a320e60..e034837 100644 --- a/sub_perco.c +++ b/sub_perco.c @@ -50,8 +50,8 @@ int writetiff(char *filename, char *description, int x, int y, int width, int he glPixelStorei(GL_PACK_ALIGNMENT, 1); glReadPixels(x, y, width, height, GL_RGB, GL_UNSIGNED_BYTE, image); - TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32) width); - TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32) height); + TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32_t) width); + TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32_t) height); TIFFSetField(file, TIFFTAG_BITSPERSAMPLE, 8); TIFFSetField(file, TIFFTAG_COMPRESSION, compression); TIFFSetField(file, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB); diff --git a/sub_sphere.c b/sub_sphere.c index b3e15e7..79c9bc6 100644 --- a/sub_sphere.c +++ b/sub_sphere.c @@ -53,6 +53,8 @@ void init_wave_sphere(t_wave_sphere wsphere[NX*NY]) 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; + wsphere[i*NY+j].radius = 1.0; + 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)); @@ -64,6 +66,8 @@ void init_wave_sphere(t_wave_sphere wsphere[NX*NY]) wsphere[i*NY+j].x2d = xy[0]; wsphere[i*NY+j].y2d = xy[1]; + + wsphere[i*NY+j].cos_angle_sphere = wsphere[i*NY+j].x*light[0] + wsphere[i*NY+j].y*light[1] + wsphere[i*NY+j].z*light[2]; } /* cotangent, taking care of not dividing by zero */ @@ -74,10 +78,524 @@ void init_wave_sphere(t_wave_sphere wsphere[NX*NY]) } +void read_negative_dem_values(double *height_values, t_wave_sphere wsphere[NX*NY]) +/* init bathymetric data */ +{ + int i, j, k, ii, jj, nx, ny, maxrgb, nmaxpixels = 6480000, hmin, hmax, ishift, nshift, sshift, rgbval, scan, rgbtot; + int *rgb_values, *int_height_values; + int hcont = 50, rgbdiff; + double cratio, rx, ry, height; + double *height_values_tmp, *height_values_tmp2; + FILE *image_file; + + printf("Reading bathymetric data\n"); + + rgb_values = (int *)malloc(3*nmaxpixels*sizeof(int)); + int_height_values = (int *)malloc(3*nmaxpixels*sizeof(int)); + height_values_tmp = (double *)malloc(NX*NY*sizeof(double)); + height_values_tmp2 = (double *)malloc(NX*NY*sizeof(double)); + +// image_file = fopen("bathymetry_gebco_3600x1800_color.ppm", "r"); + image_file = fopen("bathymetry_gebco_2560_1280_mod2_color.ppm", "r"); + scan = fscanf(image_file,"%i %i\n", &nx, &ny); + scan = fscanf(image_file,"%i\n", &maxrgb); + + for (i=0; i nmaxpixels) + { + printf("bathymetric data file too large, increase nmaxpixels in read_negative_dem_values()\n"); + exit(0); + } + + /* shift due to min/max latitudes of image */ + sshift = 0 + DPOLE; + nshift = 0 + DPOLE; + + /* choice of zero meridian */ + ishift = (int)(nx*ZERO_MERIDIAN/360.0); + + /* read rgb values */ + for (j=0; j hcont)) hmin = rgbtot; + if (rgbtot > hmax) hmax = rgbtot; + int_height_values[3*(j*nx+i)] = rgbtot; + } + + printf("hmin = %i, hmax = %i\n", hmin, hmax); + + /* remove remaining black continents */ + for (i=0; i<3*nx*ny; i++) if (int_height_values[i] < hcont) int_height_values[i] = hmax; + + cratio = 1.0/(double)(hmax-hmin); + rx = (double)nx/(double)NX; + ry = (double)ny/(double)(NY - sshift - nshift); + + /* option: set continents color to white */ +// for (i=0; i nx-1) ii -= nx; + jj = (int)(ry*(double)(NY-nshift - j)); + if (jj > ny-1) jj = ny-1; + if (jj < 0) jj = 0; + + if (wsphere[i*NY+j].indomain) + { + /* set height to zero if color is black (due to black patches in bathymetric map) */ + if (int_height_values[3*(jj*nx+ii)] < hcont) height = 0.0; + else height = -1.0 + (double)(int_height_values[3*(jj*nx+ii)])*cratio; + if (height > 0.0) height = 0.0; + height_values_tmp[i*NY+j] = height; + wsphere[i*NY+j].altitude = height; + + if (int_height_values[3*(jj*nx+ii)] > hcont) + { + wsphere[i*NY+j].r = 0.9*(double)rgb_values[3*(jj*nx+ii)]*cratio; + wsphere[i*NY+j].g = 0.9*(double)rgb_values[3*(jj*nx+ii)+1]*cratio; + wsphere[i*NY+j].b = 0.9*(double)rgb_values[3*(jj*nx+ii)+2]*cratio; + } + else + { + wsphere[i*NY+j].r = 0.29; + wsphere[i*NY+j].g = 0.29; + wsphere[i*NY+j].b = 0.29; + } + } + else + { + height_values_tmp[i*NY+j] = 0.0; + height_values_tmp2[i*NY+j] = 0.0; + } + } + + /* smoothen values at low depth */ + if (SMOOTH_DEM) for (k=1; k= -0.25)) + { + height_values_tmp2[i*NY+j] = height_values_tmp[i*NY+j] + 0.1*(height_values_tmp[(i+1)*NY+j] + height_values_tmp[(i-1)*NY+j] + height_values_tmp[i*NY+j+1] + height_values_tmp[i*NY+j-1] - 4.0*height_values_tmp[i*NY+j]); + + height_values_tmp[i*NY+j] = height_values_tmp2[i*NY+j] + 0.1*(height_values_tmp2[(i+1)*NY+j] + height_values_tmp2[(i-1)*NY+j] + height_values_tmp2[i*NY+j+1] + height_values_tmp2[i*NY+j-1] - 4.0*height_values_tmp2[i*NY+j]); + } + } + + if (SMOOTH_DEM) for (i=1; i= -0.25)) + { + wsphere[i*NY+j].altitude = height_values_tmp[i*NY+j]; + } + + for (i=0; i nmaxpixels) + { + printf("DEM too large, increase nmaxpixels in init_dem()\n"); + exit(0); + } + + /* shift due to min/max latitudes of image */ + sshift = 0 + DPOLE; + nshift = 0 + DPOLE; + + /* choice of zero meridian */ + ishift = (int)(nx*ZERO_MERIDIAN/360.0); + + printf("Reading RGB values\n"); + + /* read rgb values */ + for (j=0; j hmax) hmax = rgbval; + } + + printf("hmin = %i, hmax = %i, hsea = %i\n", hmin, hmax, (int)hsea); + + cratio = 1.0/(double)(hmax-hmin); + rx = (double)nx/(double)NX; + ry = (double)ny/(double)(NY - sshift - nshift); + + /* build height table */ + for (i=0; i nx-1) ii -= nx; + if (ii < 0) ii = 0; + jj = (int)(ry*(double)(NY-nshift - j)); + if (jj > ny-1) jj = ny-1; + if (jj < 0) jj = 0; + height_values[i*NY+j] = ((double)rgb_values[3*(jj*nx+ii)]-hsea)*cratio; + wsphere[i*NY+j].altitude = ((double)rgb_values[3*(jj*nx+ii)]-hsea)*cratio; + + if (OTHER_PLANET) + wsphere[i*NY+j].indomain = (wsphere[i*NY+j].altitude < 0.0); + +// if (wsphere[i*NY+j].indomain) printf("rgb = %i, altitude = %.3lg\n", rgb_values[3*(jj*nx+ii)], height_values[i*NY+j]); + } + + /* smooth values in case of high resolution */ + if ((NX > nx)||(NY > ny)) + { + for (i=1; i 0.0)&&(wsphere[i*NY+j].altitude <= 0.15)) +// printf("Smoothed altitude at (%i, %i): %.3lg\n", i, j, wsphere[i*NY+j].altitude); + + if ((ADD_NEGATIVE_DEM)&&(dem_number == DEM_EARTH)) + read_negative_dem_values(height_values, wsphere); + + /* set radius */ + for (i=0; i nmaxpixels) { - printf("Image too large, increase nmaxpixels in choose_colors()\n"); + printf("Image too large, increase nmaxpixels in init_earth_map()\n"); exit(0); } @@ -100,6 +618,9 @@ void init_earth_map(t_wave_sphere wsphere[NX*NY]) sshift = 0 + DPOLE; nshift = 0 + DPOLE; + /* choice of zero meridian */ + ishift = (int)(nx*ZERO_MERIDIAN/360.0); + /* read rgb values */ for (j=0; j nx-1) ii -= nx; // jj = (int)(-ry*(double)j + cy); // jj = (int)(ry*(double)(NY-nshift - j)) + sshift; @@ -133,13 +654,157 @@ void init_earth_map(t_wave_sphere wsphere[NX*NY]) 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); + wsphere[i*NY+j].indomain = (diff < 15); } + /* smooth colors in case of high resolution */ + if ((NX > nx)||(NY > ny)) + for (i=1; i nmaxpixels) + { + printf("Image too large, increase nmaxpixels in init_planet_map()\n"); + exit(0); + } + + /* shift due to min/max latitudes of image */ + sshift = 0 + DPOLE; + nshift = 0 + DPOLE; + + /* choice of zero meridian */ + ishift = (int)(nx*ZERO_MERIDIAN/360.0); + + /* read rgb values */ + for (j=0; j nx-1) ii -= nx; +// jj = (int)(-ry*(double)j + cy); +// jj = (int)(ry*(double)(NY-nshift - j)) + sshift; + jj = (int)(ry*(double)(NY-nshift - j)); + if (jj > ny-1) jj = ny-1; + if (jj < 0) jj = 0; + wsphere[i*NY+j].r = (double)rgb_values[3*(jj*nx+ii)]*cratio; + wsphere[i*NY+j].g = (double)rgb_values[3*(jj*nx+ii)+1]*cratio; + wsphere[i*NY+j].b = (double)rgb_values[3*(jj*nx+ii)+2]*cratio; + +// printf("RGB = (%.2f, %.2f, %.2f)\n", wsphere[i*NY+j].r, wsphere[i*NY+j].g, wsphere[i*NY+j].b); + + /* decide which points are in the Sea */ + wsphere[i*NY+j].indomain = 1; +// wsphere[i*NY+j].indomain = 0; +// wsphere[i*NY+j].indomain = (diff < 15); + } + + /* smooth colors in case of high resolution */ + if ((NX > nx)||(NY > ny)) + for (i=1; i RMAX) r = RMAX; @@ -568,6 +1333,7 @@ void compute_light_angle_sphere(short int xy_in[NX*NY], t_wave wave[NX*NY], t_wa wave[i*NY+j].cos_angle = pscal/norm; } +// else wave[i*NY+j].cos_angle = wsphere[i*NY+j].cos_angle; else { pscal = wsphere[i*NY+j].x*light[0] + wsphere[i*NY+j].y*light[1] + wsphere[i*NY+j].z*light[2]; @@ -596,10 +1362,11 @@ void compute_light_angle_sphere(short int xy_in[NX*NY], t_wave wave[NX*NY], t_wa } } -void compute_light_angle_sphere_2d(short int xy_in[NX*NY], t_wave wave[NX*NY], int movie) +void compute_light_angle_sphere_2d(short int xy_in[NX*NY], double phi[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; + short int draw; double gradx, grady, norm, pscal; static double dx, dy, vscale2; static int first = 1; @@ -616,7 +1383,11 @@ void compute_light_angle_sphere_2d(short int xy_in[NX*NY], t_wave wave[NX*NY], i for (i=1; i= wsphere[i*NY+j].altitude + 0.001)) + if (FLOODING) draw = ((wsphere[i*NY+j].indomain)||(phi[i*NY+j] >= wsphere[i*NY+j].altitude + 0.001)); + else draw = ((TWOSPEEDS)||(xy_in[i*NY+j])); + if (draw) { 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; @@ -626,12 +1397,14 @@ void compute_light_angle_sphere_2d(short int xy_in[NX*NY], t_wave wave[NX*NY], i wave[i*NY+j].cos_angle = pscal/norm; } + else wave[i*NY+j].cos_angle = wsphere[i*NY+j].cos_angle; } /* i=0 */ for (j=1; j= wsphere[j].altitude + 0.001)) { 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; @@ -641,12 +1414,14 @@ void compute_light_angle_sphere_2d(short int xy_in[NX*NY], t_wave wave[NX*NY], i wave[j].cos_angle = pscal/norm; } + else wave[j].cos_angle = wsphere[j].cos_angle; } /* i=NX-1 */ for (j=1; j= wsphere[(NX-1)*NY+j].altitude + 0.001)) { 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; @@ -656,6 +1431,7 @@ void compute_light_angle_sphere_2d(short int xy_in[NX*NY], t_wave wave[NX*NY], i wave[(NX-1)*NY+j].cos_angle = pscal/norm; } + else wave[(NX-1)*NY+j].cos_angle = wsphere[(NX-1)*NY+j].cos_angle; } } @@ -694,6 +1470,7 @@ void compute_cfield_sphere(short int xy_in[NX*NY], int cplot, int palette, t_wav 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; + short int draw; double x, y, ca; static int ishift; static double dx, dy; @@ -702,7 +1479,7 @@ void draw_wave_sphere_2D(int movie, double phi[NX*NY], double psi[NX*NY], short { 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); + else if (SHADE_2D) compute_light_angle_sphere_2d(xy_in, phi, wave, wsphere, 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); @@ -714,14 +1491,16 @@ void draw_wave_sphere_2D(int movie, double phi[NX*NY], double psi[NX*NY], short for (i=0; i= wsphere[i*NY+j].altitude + 0.001); + else draw = (TWOSPEEDS)||(xy_in[i*NY+j]); + if (draw) 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) + if (PLANET) 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); } @@ -744,16 +1523,22 @@ void draw_wave_sphere_ij(int i, int iplus, int j, int jplus, int jcolor, int mov /* draw wave at simulation grid point (i,j) */ { int k, l; + short int draw; 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]); +// if ((TWOSPEEDS)||(xy_in[i*NY+j])) +// if (wsphere[i*NY+j].indomain) +// if ((wsphere[i*NY+j].indomain)||(phi[i*NY+j] >= wsphere[i*NY+j].altitude + 0.001)) + if (FLOODING) draw = (wsphere[i*NY+j].indomain)||(phi[i*NY+j] >= wsphere[i*NY+j].altitude + 0.001); + else draw = ((TWOSPEEDS)||(xy_in[i*NY+j])); + if (draw) 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 = wave[i*NY+j].cos_angle; + ca = wsphere[i*NY+j].cos_angle; ca = (ca + 1.0)*0.4 + 0.2; if (fade) ca *= fade_value; - if (B_DOMAIN == D_SPHERE_EARTH) + if (PLANET) 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); } @@ -781,7 +1566,7 @@ void draw_wave_sphere_3D(int movie, double phi[NX*NY], double psi[NX*NY], short { 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); + else if (SHADE_2D) compute_light_angle_sphere_2d(xy_in, phi, wave, wsphere, movie); compute_cfield_sphere(xy_in, cplot, palette, wave, fade, fade_value, movie); } diff --git a/sub_wave.c b/sub_wave.c index c775d96..dd4bdd2 100644 --- a/sub_wave.c +++ b/sub_wave.c @@ -30,8 +30,8 @@ int writetiff_new(char *filename, char *description, int x, int y, int width, in glPixelStorei(GL_PACK_ALIGNMENT, 1); glReadPixels(x, y, width, height, GL_RGB, GL_UNSIGNED_BYTE, image); - TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32) width); - TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32) height); + TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32_t) width); + TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32_t) height); TIFFSetField(file, TIFFTAG_BITSPERSAMPLE, 8); TIFFSetField(file, TIFFTAG_COMPRESSION, compression); TIFFSetField(file, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB); @@ -85,8 +85,8 @@ int writetiff(char *filename, char *description, int x, int y, int width, int he glPixelStorei(GL_PACK_ALIGNMENT, 1); glReadPixels(x, y, width, height, GL_RGB, GL_UNSIGNED_BYTE, image); - TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32) width); - TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32) height); + TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32_t) width); + TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32_t) height); TIFFSetField(file, TIFFTAG_BITSPERSAMPLE, 8); TIFFSetField(file, TIFFTAG_COMPRESSION, compression); TIFFSetField(file, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB); @@ -3130,6 +3130,13 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_ } return(1); } + case (D_LENS): + { + if (vabs(y) > MU) return(1); + a = LAMBDA - 0.25; + if (x > 0) return ((x+a)*(x+a) + y*y > LAMBDA*LAMBDA); + return ((x-a)*(x-a) + y*y > LAMBDA*LAMBDA); + } case (D_MENGER): { x1 = 0.5*(x+1.0); @@ -4537,6 +4544,22 @@ void draw_billiard(int fade, double fade_value) /* draws the billiard bound } break; } + case (D_LENS): + { + a = LAMBDA - 0.25; + if (first) + { + arcangle = asin(MU/LAMBDA); + ll = LAMBDA*cos(arcangle) - a; + first = 0; + } + + draw_circle_arc(-a, 0.0, LAMBDA, -arcangle, 2.0*arcangle, NSEG); + draw_line(ll, MU, -ll, MU); + draw_circle_arc(a, 0.0, LAMBDA, PI-arcangle, 2.0*arcangle, NSEG); + draw_line(-ll, -MU, ll, -MU); + break; + } case (D_MENGER): { glLineWidth(3); @@ -6330,7 +6353,7 @@ void init_wave_packets(t_wave_packet *packet, int radius) printf("Initialising wave packets\n"); switch (WAVE_PACKET_SOURCE_TYPE) { - case (0): + case (WP_RANDOM1): { nx = (int)sqrt((double)N_WAVE_PACKETS); ny = N_WAVE_PACKETS/nx; @@ -6356,15 +6379,12 @@ void init_wave_packets(t_wave_packet *packet, int radius) break; } - case (1): + case (WP_RANDOM2): { for (i=0; i