/*********************************************************************************/ /* */ /* Animation of interacting particles in a planar domain */ /* */ /* N. Berglund, november 2021 */ /* */ /* UPDATE 24/04: distinction between damping and "elasticity" parameters */ /* UPDATE 27/04: new billiard shapes, bug in color scheme fixed */ /* UPDATE 28/04: code made more efficient, with help of Marco Mancini */ /* */ /* 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 lennardjones lennardjones.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_ljones */ /* It may be possible to increase parameter PAUSE */ /* */ /* create movie using */ /* ffmpeg -i lj.%05d.tif -vcodec libx264 lj.mp4 */ /* */ /*********************************************************************************/ #include #include #include #include #include #include #include /* Sam Leffler's libtiff library. */ #include #include #define MOVIE 0 /* set to 1 to generate movie */ #define DOUBLE_MOVIE 0 /* set to 1 to produce movies for wave height and energy simultaneously */ #define SAVE_MEMORY 1 /* set to 1 to save memory while saving frames */ #define NO_EXTRA_BUFFER_SWAP 1 /* some OS require one less buffer swap when recording images */ #define TIME_LAPSE 0 /* set to 1 to add a time-lapse movie at the end */ /* so far incompatible with double movie */ #define TIME_LAPSE_FACTOR 3 /* factor of time-lapse movie */ #define TIME_LAPSE_FIRST 1 /* set to 1 to show time-lapse version first */ #define SAVE_TIME_SERIES 0 /* set to 1 to save time series of particle positions */ /* General geometrical parameters */ #define WINWIDTH 1760 /* window width */ #define WINHEIGHT 990 /* window height */ #define XMIN -2.0 #define XMAX 2.0 /* x interval */ #define YMIN -1.125 #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ #define INITXMIN -2.1 #define INITXMAX -2.0 /* x interval for initial condition */ #define INITYMIN -0.7 #define INITYMAX -0.6 /* y interval for initial condition */ #define THERMOXMIN -1.25 #define THERMOXMAX 1.25 /* x interval for initial condition */ #define THERMOYMIN 0.0 #define THERMOYMAX 0.75 /* y interval for initial condition */ #define ADDXMIN -2.1 #define ADDXMAX -2.0 /* x interval for adding particles */ #define ADDYMIN -0.6 #define ADDYMAX -0.5 /* y interval for adding particles */ #define ADDRMIN 4.75 #define ADDRMAX 6.0 /* r interval for adding particles */ #define BCXMIN -2.2 #define BCXMAX 2.0 /* x interval for boundary condition */ #define BCYMIN -1.125 #define BCYMAX 1.4 /* y interval for boundary condition */ #define OBSXMIN -2.0 #define OBSXMAX 2.0 /* x interval for motion of obstacle */ #define CIRCLE_PATTERN 1 /* pattern of circles, see list in global_ljones.c */ #define ADD_INITIAL_PARTICLES 0 /* set to 1 to add a second type of particles */ #define CIRCLE_PATTERN_B 0 /* pattern of circles for additional particles */ #define ADD_FIXED_OBSTACLES 0 /* set to 1 do add fixed circular obstacles */ #define OBSTACLE_PATTERN 82 /* pattern of obstacles, see list in global_ljones.c */ #define RATTLE_OBSTACLES 1 /* set to 1 to rattle obstacles (for pattern O_SIEVE_B) */ #define ADD_FIXED_SEGMENTS 1 /* set to 1 to add fixed segments as obstacles */ #define SEGMENT_PATTERN 361 /* 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 BELT_SPEED1 10.0 /* speed of first conveyor belt */ #define BELT_SPEED2 15.0 /* speed of second conveyor belt */ #define BELT_SPEED3 6.0 /* speed of second conveyor belt */ #define OBSTACLE_OMEGA 300.0 /* obstacle rotation speed */ #define TWO_TYPES 0 /* set to 1 to have two types of particles */ #define TYPE_PROPORTION 0.5 /* proportion of particles of first type */ #define TWOTYPE_CONFIG 0 /* choice of types, see TTC_ list in global_ljones.c */ #define SYMMETRIZE_FORCE 1 /* set to 1 to symmetrize two-particle interaction, only needed if particles are not all the same */ #define CENTER_PX 0 /* set to 1 to center horizontal momentum */ #define CENTER_PY 0 /* set to 1 to center vertical momentum */ #define CENTER_PANGLE 0 /* set to 1 to center angular momentum */ #define INTERACTION 1 /* particle interaction, see list in global_ljones.c */ #define INTERACTION_B 1 /* particle interaction for second type of particle, see list in global_ljones.c */ #define SPIN_INTER_FREQUENCY 4.0 /* angular frequency of spin-spin interaction */ #define SPIN_INTER_FREQUENCY_B 4.0 /* angular frequency of spin-spin interaction for second particle type */ #define MOL_ANGLE_FACTOR 4.0 /* rotation angle for P_MOL_ANGLE color scheme */ #define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */ #define NPOISSON 100 /* number of points for Poisson C_RAND_POISSON arrangement */ #define PDISC_DISTANCE 1.0 /* minimal distance in Poisson disc process, controls density of particles */ #define PDISC_CANDIDATES 100 /* number of candidates in construction of Poisson disc process */ #define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */ #define LAMBDA 0.2 /* parameter controlling the dimensions of domain */ #define MU 0.02 /* parameter controlling radius of particles */ #define MU_B 0.03 /* parameter controlling radius of particles of second type */ #define NPOLY 3 /* number of sides of polygon */ #define APOLY 0.075 /* 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 1 /* number of grid point for grid of disks */ #define NGRIDY 1 /* 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 */ #define DAM_WIDTH 0.05 /* width of dam for S_DAM segment configuration */ #define NOBSX 10 #define NOBSY 5 /* obstacles for O_HEX obstacle pattern */ #define NTREES 15 /* number of trees in S_TREES */ #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 */ /* Parameters for length and speed of simulation */ #define NSTEPS 4800 /* number of frames of movie */ #define NVID 120 /* number of iterations between images displayed on screen */ #define NSEG 25 /* number of segments of boundary of circles */ #define INITIAL_TIME 0 /* time after which to start saving frames */ #define OBSTACLE_INITIAL_TIME 0 /* time after which to start moving obstacle */ #define BOUNDARY_WIDTH 1 /* width of particle boundary */ #define LINK_WIDTH 2 /* width of links between particles */ #define CONTAINER_WIDTH 2 /* width of container boundary */ #define PAUSE 1000 /* number of frames after which to pause */ #define PSLEEP 1 /* 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 */ /* Boundary conditions, see list in global_ljones.c */ #define BOUNDARY_COND 1 /* Plot type, see list in global_ljones.c */ #define PLOT 23 #define PLOT_B 13 /* plot type for second movie */ /* Background color depending on particle properties */ #define COLOR_BACKGROUND 0 /* set to 1 to color background */ #define BG_COLOR 2 /* type of background coloring, see list in global_ljones.c */ #define BG_COLOR_B 0 /* type of background coloring, see list in global_ljones.c */ #define DRAW_BONDS 0 /* 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 DRAW_CLUSTER_LINKS 0 /* set to 1 to draw links between particles in cluster */ #define ALTITUDE_LINES 0 /* set to 1 to add horizontal lines to show altitude */ #define COLOR_SEG_GROUPS 0 /* set to 1 to collor segment groups differently */ #define N_PARTICLE_COLORS 300 /* 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-averaging in P_EMEAN color scheme */ #define DRATIO 0.999 /* ratio for time-averaging in P_DIRECT_EMEAN color scheme */ /* Color schemes */ #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 0 /* Color palette for angle representation */ #define COLOR_PALETTE_DIRECTION 17 /* Color palette for direction representation */ #define COLOR_PALETTE_INITIAL_POS 10 /* Color palette for initial position representation */ #define COLOR_PALETTE_DIFFNEIGH 10 /* Color palette for different neighbours representation */ #define COLOR_PALETTE_PRESSURE 11 /* Color palette for different neighbours representation */ #define COLOR_PALETTE_CHARGE 18 /* Color palette for charge representation */ #define COLOR_PALETTE_CLUSTER 14 /* Color palette for cluster representation */ #define COLOR_PALETTE_CLUSTER_SIZE 13 /* Color palette for cluster size representation */ #define COLOR_PALETTE_CLUSTER_SELECTED 11 /* Color palette for selected cluster representation */ #define COLOR_HUE_CLUSTER_SELECTED 90.0 /* Color hue for selected cluster */ #define COLOR_HUE_CLUSTER_NOT_SELECTED 220.0 /* Color hue for selected cluster */ #define BLACK 1 /* background */ #define COLOR_SCHEME 3 /* choice of color scheme, see list in global_ljones.c */ #define SCALE 0 /* set to 1 to adjust color scheme to variance of field */ #define SLOPE 0.5 /* sensitivity of color on wave amplitude */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ #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 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 1 /* set to 1 to print certain parameters */ #define PRINT_TEMPERATURE 0 /* set to 1 to print current temperature */ #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 0 /* set to 1 to print force on segments */ #define PRINT_NPARTICLES 0 /* print number of active particles */ #define PRINT_TYPE_PROP 0 /* print type proportion */ #define FORCE_FACTOR 0.1 /* factor controlling length of force vector */ /* particle properties */ #define ENERGY_HUE_MIN 330.0 /* color of original particle */ #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_EMIN 10.0 /* energy of particle with coolest color */ #define PARTICLE_EMAX 2000.0 /* energy of particle with hottest color */ #define HUE_TYPE0 320.0 /* hue of particles of type 0 */ #define HUE_TYPE1 60.0 /* hue of particles of type 1 */ #define HUE_TYPE2 320.0 /* hue of particles of type 2 */ #define HUE_TYPE3 260.0 /* hue of particles of type 3 */ #define HUE_TYPE4 200.0 /* hue of particles of type 4 */ #define HUE_TYPE5 60.0 /* hue of particles of type 5 */ #define HUE_TYPE6 130.0 /* hue of particles of type 6 */ #define HUE_TYPE7 150.0 /* hue of particles of type 7 */ #define BG_FORCE_SLOPE 7.5e-8 /* contant in BG_FORCE backgound color scheme*/ #define RANDOM_RADIUS 1 /* set to 1 for random particle radius */ #define RANDOM_RADIUS_MIN 0.25 /* min of random particle radius (default 0.75) */ #define RANDOM_RADIUS_RANGE 1.5 /* range of random particle radius (default 0.5) */ #define ADAPT_MASS_TO_RADIUS 0 /* set to positive value to for mass prop to power of radius */ #define ADAPT_DAMPING_TO_RADIUS 0.5 /* set to positive value to for friction prop to power of radius */ #define ADAPT_DAMPING_FACTOR 0.5 /* factor by which damping is adapted to radius */ #define DT_PARTICLE 3.0e-6 /* time step for particle displacement */ #define KREPEL 50.0 /* constant in repelling force between particles */ #define EQUILIBRIUM_DIST 2.5 /* Lennard-Jones equilibrium distance */ #define EQUILIBRIUM_DIST_B 2.5 /* Lennard-Jones equilibrium distance for second type of particle */ #define SEGMENT_FORCE_EQR 1.0 /* equilibrium distance factor for force from segments (default 1.5) */ #define REPEL_RADIUS 25.0 /* radius in which repelling force acts (in units of particle radius) */ #define DAMPING 150.0 /* damping coefficient of particles */ #define INITIAL_DAMPING 1000.0 /* damping coefficient of particles during initial phase */ #define DAMPING_ROT 5.0 /* damping coefficient for rotation of particles */ #define PARTICLE_MASS 2.0 /* mass of particle of radius MU */ #define PARTICLE_MASS_B 2.0 /* mass of particle of radius MU_B */ #define PARTICLE_INERTIA_MOMENT 0.5 /* moment of inertia of particle */ #define PARTICLE_INERTIA_MOMENT_B 0.5 /* moment of inertia of second type of particle */ #define V_INITIAL 800.0 /* initial velocity range */ #define OMEGA_INITIAL 100.0 /* initial angular velocity range */ #define VICSEK_VMIN 1.0 /* minimal speed of particles in Vicsek model */ #define VICSEK_VMAX 40.0 /* minimal speed of particles in Vicsek model */ #define COULOMB_LJ_FACTOR 1.0 /* relative intensity of LJ interaction in I_COULOMB_LJ interaction (default: 0.01) */ #define V_INITIAL_TYPE 0 /* 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.004 /* initial inverse temperature */ #define MU_XI 0.005 /* friction constant in thermostat */ #define KSPRING_BOUNDARY 2.0e11 /* confining harmonic potential outside simulation region */ #define KSPRING_OBSTACLE 1.0e9 /* harmonic potential of obstacles */ #define NBH_DIST_FACTOR 4.0 /* radius in which to count neighbours */ #define GRAVITY 5000.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 1 /* type of gravity schedule, see list in global_ljones.c */ #define GRAVITY_FACTOR 10.0 /* factor by which to increase gravity */ #define GRAVITY_INITIAL_TIME 250 /* time at start of simulation with constant gravity */ #define GRAVITY_RESTORE_TIME 500 /* time at end of simulation with gravity restored to initial value */ #define KSPRING_VICSEK 0.2 /* spring constant for I_VICSEK_SPEED interaction */ #define VICSEK_REPULSION 10.0 /* repulsion between particles in Vicsek model */ #define ADD_EFIELD 0 /* set to 1 to add an electric field */ #define EFIELD 30000.0 /* value of electric field */ #define ADD_BFIELD 0 /* set to 1 to add a magnetic field */ #define BFIELD 225.0 /* value of magnetic field */ #define CHARGE 1.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 5000000.0 /* factor by which to increase electric field */ #define INCREASE_B 0 /* set to 1 to increase magnetic field */ #define BFIELD_FACTOR 1.0 /* factor by which to increase magnetic field */ #define CHARGE_OBSTACLES 0 /* set to 1 for obstacles to be charged */ #define OBSTACLE_CHARGE 3.0 /* charge of obstacles */ #define KCOULOMB_OBSTACLE 1000.0 /* Coulomb force constant for charged obstacles */ #define BFIELD_REGION 1 /* space-dependent magnetic field (0 for constant) */ #define ADD_WIND 1 /* set to 1 to add a "wind" friction force */ #define WIND_FORCE 1.35e6 /* force of wind */ #define WIND_YMIN -0.6 /* min altitude of region with wind */ #define ROTATION 1 /* set to 1 to include rotation of particles */ #define COUPLE_ANGLE_TO_THERMOSTAT 0 /* set to 1 to couple angular degrees of freedom to thermostat */ #define DIMENSION_FACTOR 1.0 /* scaling factor taking into account number of degrees of freedom */ #define KTORQUE 1.0e5 /* force constant in angular dynamics */ #define KTORQUE_BOUNDARY 1.0e5 /* constant in torque from the boundary */ #define KTORQUE_B 10.0 /* force constant in angular dynamics */ #define KTORQUE_DIFF 500.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 0 /* set to 1 to draw cross on particles of second type */ #define DRAW_MINUS 0 /* set to 1 to draw cross on particles of negative charge */ #define SPIN_RANGE 10.0 /* range of spin-spin interaction */ #define SPIN_RANGE_B 10.0 /* range of spin-spin interaction for second type of particle */ #define QUADRUPOLE_RATIO 0.6 /* anisotropy in quadrupole potential */ #define INCREASE_BETA 0 /* set to 1 to increase BETA during simulation */ #define BETA_SCHEDULE 3 /* type of temperature schedule, see TS_* in global_ljones */ #define BETA_FACTOR 0.06 /* factor by which to change BETA during simulation */ #define TS_SLOPE 8.5 /* controls speed of change of BETA for TS_TANH schedule (default 1.0) */ #define N_TOSCILLATIONS 1.0 /* number of temperature oscillations in BETA schedule */ #define NO_OSCILLATION 0 /* set to 1 to have exponential BETA change only */ #define MIDDLE_CONSTANT_PHASE 0 /* final phase in which temperature is constant */ #define FINAL_DECREASE_PHASE 0 /* final phase in which temperature decreases */ #define FINAL_CONSTANT_PHASE -1 /* final phase in which temperature is constant */ #define DECREASE_CONTAINER_SIZE 0 /* set to 1 to decrease size of container */ #define SMOOTH_CONTAINER_DECREASE 1 /* set to 1 to decrease size smoothly at each simulation step */ #define SYMMETRIC_DECREASE 0 /* set tp 1 to decrease container symmetrically */ #define COMPRESSION_RATIO 0.25 /* final size of container */ #define RESTORE_CONTAINER_SIZE 1 /* set to 1 to restore container to initial size at end of simulation */ #define RESTORE_TIME 800 /* time before end of sim at which to restore size */ #define MOVE_OBSTACLE 0 /* set to 1 to have a moving obstacle */ #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.018 /* 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 */ #define RECORD_PRESSURES 0 /* set to 1 to record pressures on obstacle */ #define N_PRESSURES 100 /* number of intervals to record pressure */ #define N_P_AVERAGE 100 /* size of pressure averaging window */ #define N_T_AVERAGE 1 /* size of temperature averaging window */ #define MAX_PRESSURE 3.0e10 /* pressure shown in "hottest" color */ #define PARTIAL_THERMO_COUPLING 0 /* set to 1 to couple only some particles to thermostat */ #define PARTIAL_THERMO_REGION 9 /* region for partial thermostat coupling (see list in global_ljones.c) */ #define PARTIAL_THERMO_SHIFT 0.2 /* distance from obstacle at the right of which particles are coupled to thermostat */ #define PARTIAL_THERMO_WIDTH 1.0 /* vertical size of partial thermostat coupling */ #define PARTIAL_THERMO_HEIGHT 0.0 /* vertical size of partial thermostat coupling */ #define PARTIAL_THERMO_RIN 0.5 /* initial radius of region without coupling */ #define PARTIAL_THERMO_RFIN 1.3 /* final radius of region without coupling */ #define INCREASE_KREPEL 0 /* set to 1 to increase KREPEL during simulation */ #define KREPEL_FACTOR 100.0 /* factor by which to change KREPEL during simulation */ #define PART_AT_BOTTOM 0 /* set to 1 to include "seed" particles at bottom */ #define MASS_PART_BOTTOM 10000.0 /* mass of particles at bottom */ #define NPART_BOTTOM 100 /* number of particles at the bottom */ #define ADD_PARTICLES 1 /* set to 1 to add particles */ #define ADD_REGION 0 /* shape of add regions, cf ADD_* in global_ljones */ #define ADD_TIME 0 /* time at which to add first particle */ #define ADD_PERIOD 15 /* time interval between adding further particles */ #define N_ADD_PARTICLES 1 /* number of particles to add */ #define FINAL_NOADD_PERIOD 500 /* final period where no particles are added */ #define SAFETY_FACTOR 4.0 /* no particles are added at distance less than MU*SAFETY_FACTOR of other particles */ #define TRACER_PARTICLE 1 /* set to 1 to have a tracer particle */ #define N_TRACER_PARTICLES 500 /* number of tracer particles */ #define TRAJECTORY_LENGTH 4800 /* length of recorded trajectory */ #define TRACER_LUM_FACTOR 5.0 /* controls luminosity decrease of trajectories with time */ #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 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 150 /* 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 */ #define DEACTIVATE_SEGMENT 1 /* set to 1 to deactivate last segment after a certain time */ #define SEGMENT_DEACTIVATION_TIME 200 /* 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_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 */ #define SEGMENT_GROUP_I 1000.0 /* moment of inertia of segment group */ #define SEGMENT_GROUP_DAMPING 0.0 /* damping of segment groups */ #define GROUP_REPULSION 0 /* set to 1 for groups of segments to repel each other */ #define KSPRING_GROUPS 5.0e11 /* harmonic potential between segment groups */ #define KSPRING_BELT 1.0e4 /* spring constant from belt */ #define GROUP_WIDTH 0.05 /* interaction width of groups */ #define GROUP_G_REPEL 0 /* set to 1 to add repulsion between centers of mass of groups */ #define GROUP_G_REPEL_RADIUS 1.2 /* radius within which centers of mass of groups repel each other */ #define TRACK_SEGMENT_GROUPS 0 /* set to 1 for view to track group of segments */ #define TRACK_X_PADDING 2.0 /* distance from x boundary where tracking starts */ #define POSITION_DEPENDENT_TYPE 0 /* set to 1 to make particle type depend on initial position */ #define POSITION_Y_DEPENDENCE 0 /* set to 1 for the separation between particles to be horizontal */ #define POSITION_DEP_SIGN -1.0 /* sign in position dependence condition */ #define POSITION_DEP_X -0.625 /* threshold value for position-dependent type */ #define PRINT_ENTROPY 0 /* set to 1 to compute entropy */ #define REACTION_DIFFUSION 0 /* set to 1 to simulate a chemical reaction (particles may change type) */ #define RD_REACTION 262 /* type of reaction, see list in global_ljones.c */ #define RD_TYPES 1 /* number of types in reaction-diffusion equation */ #define RD_INITIAL_COND 0 /* initial condition of particles */ #define REACTION_DIST 1.2 /* maximal distance for reaction to occur */ #define REACTION_PROB 1.0 /* probability controlling reaction term */ #define DISSOCIATION_PROB 0.0 /* probability controlling dissociation reaction */ #define KILLING_PROB 0.0015 /* probability of enzymes being killed */ #define DELTAMAX 0.1 /* max orientation difference for pairing polygons */ #define CENTER_COLLIDED_PARTICLES 0 /* set to 1 to recenter particles upon reaction (may interfere with thermostat) */ #define EXOTHERMIC 0 /* set to 1 to make reaction exo/endothermic */ #define DELTA_EKIN 2000.0 /* change of kinetic energy in reaction */ #define COLLISION_TIME 25 /* time during which collisions are shown */ #define COLLISION_RADIUS 2.0 /* radius of discs showing collisions, in units of MU */ #define DELTAVMAX 200.0 /* maximal deltav allowed for pairing molecules */ #define AGREGMAX 6 /* maximal number of partners for CHEM_AGGREGATION reaction */ #define AGREG_DECOUPLE 12 /* minimal number of partners to decouple from thermostat */ #define CLUSTER_PARTICLES 0 /* set to 1 for particles to form rigid clusters */ #define CLUSTER_MAXSIZE 1000 /* max size of clusters */ #define SMALL_CLUSTER_MAXSIZE 10 /* size limitation on smaller cluster */ #define SMALL_NP_MAXSIZE 2 /* limitation on number of partners of particle in smaller cluster */ #define NOTSELECTED_CLUSTER_MAXSIZE 0 /* limit on size of clusters that can merge with non-selected cluster */ #define REPAIR_CLUSTERS 0 /* set to 1 to repair alignment in clusters */ #define REPAIR_MIN_DIST 0.75 /* relative distance below which overlapping polygons are inactivated */ #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 1.0 /* expected final number of particles over initial number */ #define PRINT_LEFT 0 /* set to 1 to print certain parameters at the top left instead of right */ #define PLOT_SPEEDS 0 /* set to 1 to add a plot of obstacle speeds (e.g. for rockets) */ #define PLOT_TRAJECTORIES 0 /* set to 1 to add a plot of obstacle trajectories (e.g. for rockets) */ #define VMAX_PLOT_SPEEDS 0.25 /* vertical scale of plot of obstacle speeds */ #define EHRENFEST_COPY 0 /* set to 1 to add equal number of larger particles (for Ehrenfest model) */ #define LID_MASS 1000.0 /* mass of lid for BC_RECTANGLE_LID b.c. */ #define LID_WIDTH 0.1 /* width of lid for BC_RECTANGLE_LID b.c. */ #define WALL_MASS 2000.0 /* mass of wall for BC_RECTANGLE_WALL b.c. */ #define WALL_FRICTION 0.0 /* friction on wall for BC_RECTANGLE_WALL b.c. */ #define WALL_WIDTH 0.015 /* width of wall for BC_RECTANGLE_WALL b.c. */ #define WALL_VMAX 100.0 /* max speed of wall */ #define WALL_TIME 0 /* time during which to keep wall */ #define CHANGE_TYPES 0 /* set to 1 to change type proportion in course of simulation */ #define PROP_MIN 0.1 /* min proportion of type 1 particles */ #define PROP_MAX 0.9 /* max proportion of type 1 particles */ #define PAIR_PARTICLES 0 /* set to 1 to form particles pairs */ #define RANDOMIZE_ANGLE 0 /* set to 1 for random orientation */ #define DEACIVATE_CLOSE_PAIRS 1 /* set to 1 to test for closeness to other particles */ #define PAIR_SAFETY_FACTOR 1.2 /* distance to deactivate divided by sum of radii */ #define THIRD_TYPE_PROPORTION 1.0 /* proportion of third type pairings, for certain pairing types */ #define KSPRING_PAIRS 5.0e10 /* spring constant for pair interaction */ #define KTORQUE_PAIRS 1.0e10 /* constant for angular coupling in pair interaction */ #define KTORQUE_PAIR_ANGLE 0.0 /* constant for coupling between orientation in pairs */ #define NPARTNERS 4 /* number of partners of particles - for DNA, set NPARTNERS_DNA */ #define NPARTNERS_DNA 8 /* number of partners of particles, case of DNA, should be at least 8 */ #define NARMS 5 /* number of "arms" for certain paring types */ #define PAIRING_TYPE 9 /* type of pairing, see POLY_ in global_ljones.c */ #define PARTNER_ANGLE 104.45 /* angle (in degrees) between ions for POLY_WATER case */ #define PAIR_DRATIO 1.1 /* ratio between equilibrium distance and radius (default: 1.0) */ #define MU_C 0.035 /* radius of partner particle */ #define PARTICLE_MASS_C 2.0 /* mass or partner particle */ #define CHARGE_C 1.0 /* charge of partner particle */ #define CLUSTER_COLOR_FACTOR 40 /* factor for initialization of cluster colors */ #define ALTERNATE_POLY_CHARGE 1 /* set to 1 for alternating charges in molecule */ #define SECONDARY_PAIRING 0 /* set to 1 to pair with secondary partners, experimental */ #define DNA_RIGIDITY 0.5 /* controls rigidity for POLY_DNA_DOUBLE pairs, default = 1 */ #define PAIR_TYPEB_PARTICLES 1 /* set to 1 to pair particle of type 1 */ #define NPARTNERS_B 5 /* number of partners of particles */ #define NARMS_B 1 /* number of "arms" for certain paring types */ #define PAIRING_TYPE_B 81 /* type of pairing, see POLY_ in global_ljones.c */ #define MU_D 0.035 /* radius of partner particle */ #define PARTICLE_MASS_D 2.0 /* mass or partner particle */ #define CHARGE_D 1.0 /* charge of partner particle */ #define NXMAZE 12 /* width of maze */ #define NYMAZE 12 /* height of maze */ #define MAZE_MAX_NGBH 4 /* max number of neighbours of maze cell */ #define RAND_SHIFT 4 /* seed of random number generator */ #define MAZE_XSHIFT 0.5 /* horizontal shift of maze */ #define MAZE_WIDTH 0.01 /* width of maze walls */ #define FLOOR_FORCE 1 /* set to 1 to limit force on particle to FMAX */ #define FMAX 1.0e9 /* maximal force */ #define FLOOR_OMEGA 0 /* set to 1 to limit particle momentum to PMAX */ #define PMAX 1000.0 /* maximal force */ #define HASHX 40 /* size of hashgrid in x direction */ #define HASHY 20 /* 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 */ #define DRAW_COLOR_SCHEME 0 /* set to 1 to plot the color scheme */ #define COLORBAR_RANGE 8.0 /* scale of color scheme bar */ #define COLORBAR_RANGE_B 12.0 /* scale of color scheme bar for 2nd part */ #define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ #define LIMIT_ENERGY 0 /* set to 1 to limit energy, when there is no thermostat */ #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_LOG_EMEAN)||(PLOT_B == P_LOG_EMEAN)||(PLOT == P_DIRECT_EMEAN)||(PLOT_B == P_DIRECT_EMEAN)||(PLOT == P_EMEAN_DENSITY)||(PLOT_B == P_EMEAN_DENSITY)) #define COMPUTE_DIRMEAN ((PLOT == P_DIRECT_EMEAN)||(PLOT_B == P_DIRECT_EMEAN)) #define COUNT_PARTNER_TYPE ((RD_REACTION == CHEM_H2O_H_OH)||(RD_REACTION == CHEM_2H2O_H3O_OH)) #define PAIR_FORCE ((PAIR_PARTICLES)||((REACTION_DIFFUSION)&&((RD_REACTION == CHEM_AGGREGATION)||(RD_REACTION == CHEM_AGGREGATION_CHARGE)||(RD_REACTION == CHEM_AGGREGATION_NNEIGH)||(RD_REACTION == CHEM_POLYGON_AGGREGATION)))) #define COMPUTE_PAIR_TORQUE (KTORQUE_PAIR_ANGLE != 0.0) #define ADD_CONVEYOR_FORCE ((ADD_FIXED_SEGMENTS)&&((SEGMENT_PATTERN == S_CONVEYOR_BELT)||(SEGMENT_PATTERN == S_TWO_CONVEYOR_BELTS)||(SEGMENT_PATTERN == S_PERIODIC_CONVEYORS)||(SEGMENT_PATTERN == S_TEST_CONVEYORS)||(SEGMENT_PATTERN == S_CONVEYOR_SHOVELS)||(SEGMENT_PATTERN == S_CONVEYOR_MIXED)||(SEGMENT_PATTERN == S_CONVEYOR_SIEVE)||(SEGMENT_PATTERN == S_CONVEYOR_SIEVE_B)||(SEGMENT_PATTERN == S_CONVEYOR_SIEVE_LONG))) #define MOVE_CONVEYOR_BELT ((ADD_FIXED_SEGMENTS)&&((SEGMENT_PATTERN == S_CONVEYOR_SHOVELS)||(SEGMENT_PATTERN == S_CONVEYOR_MIXED)||(SEGMENT_PATTERN == S_CONVEYOR_SIEVE_LONG))) #define ROTATE_OBSTACLES ((ADD_FIXED_OBSTACLES)&&((OBSTACLE_PATTERN == O_SIEVE)||(OBSTACLE_PATTERN == O_SIEVE_B)||(OBSTACLE_PATTERN == O_SIEVE_LONG))) #define POLYGON_INTERACTION ((INTERACTION == I_POLYGON)||(INTERACTION == I_POLYGON_ALIGN)) double xshift = 0.0; /* x shift of shown window */ double xspeed = 0.0; /* x speed of obstacle */ double ylid = 0.9; /* y coordinate of lid (for BC_RECTANGLE_LID b.c.) */ double vylid = 0.0; /* y speed coordinate of lid (for BC_RECTANGLE_LID b.c.) */ double xwall = 0.0; /* x coordinate of wall (for BC_RECTANGLE_WALL b.c.) */ double vxwall = 0.0; /* x speed of wall (for BC_RECTANGLE_WALL b.c.) */ double angular_speed = 0.0; /* angular speed of rotating segments */ double xtrack = 0.0; /* traking coordinate */ double ytrack = 0.0; /* traking coordinate */ double xsegments[2] = {SEGMENTS_X0, -SEGMENTS_X0}; /* x coordinate of segments (for option MOVE_BOUNDARY) */ double ysegments[2] = {SEGMENTS_Y0, SEGMENTS_Y0}; /* y coordinate of segments (for option MOVE_BOUNDARY) */ 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))) #include "global_ljones.c" #include "sub_maze.c" #include "sub_hashgrid.c" #include "sub_lj.c" FILE *lj_time_series, *lj_final_position; /*********************/ /* animation part */ /*********************/ double repel_schedule(int i) { static double kexponent; static int first = 1; double krepel; if (first) { kexponent = log(KREPEL_FACTOR)/(double)(INITIAL_TIME + NSTEPS); first = 0; } krepel = KREPEL*exp(kexponent*(double)i); printf("krepel = %.3lg\n", krepel); return(krepel); } double efield_schedule(int i) { static double efactor; static int first = 1; double efield; if (first) { efactor = EFIELD_FACTOR/(double)(NSTEPS); first = 0; } if (i < INITIAL_TIME) efield = EFIELD; else efield = EFIELD*(double)(i-INITIAL_TIME)*efactor; printf("E = %.3lg\n", efield); return(efield); } double bfield_schedule(int i) { static double bfactor; static int first = 1; double bfield; if (first) { bfactor = BFIELD_FACTOR/(double)(NSTEPS); first = 0; } if (i < INITIAL_TIME) bfield = BFIELD; else bfield = BFIELD*(double)(i-INITIAL_TIME)*bfactor; printf("B = %.3lg\n", bfield); return(bfield); } double temperature_schedule(int i) { static double bexponent, omega, bexp2, factor2, logf, ac, bc; static int first = 1, t1, t2, t3; double beta, t; if (first) { t1 = NSTEPS - MIDDLE_CONSTANT_PHASE - FINAL_DECREASE_PHASE - FINAL_CONSTANT_PHASE; t2 = NSTEPS - FINAL_DECREASE_PHASE - FINAL_CONSTANT_PHASE; t3 = NSTEPS - FINAL_CONSTANT_PHASE; bexponent = log(BETA_FACTOR)/(double)(t1); omega = N_TOSCILLATIONS*DPI/(double)(t1); logf = log(BETA_FACTOR); switch (BETA_SCHEDULE) { case (TS_EXPONENTIAL): { factor2 = BETA_FACTOR; break; } case (TS_CYCLING): { factor2 = BETA_FACTOR*2.0/(1.0 + cos(N_TOSCILLATIONS*DPI)); break; } case (TS_PERIODIC): { factor2 = exp(logf*sin(N_TOSCILLATIONS*DPI)); break; } case (TS_LINEAR): { factor2 = BETA_FACTOR; break; } case (TS_COSINE): { factor2 = BETA_FACTOR; ac = 2.0*BETA*BETA_FACTOR/(1.0 + BETA_FACTOR); bc = (BETA_FACTOR - 1.0)/(1.0 + BETA_FACTOR); break; } case (TS_EXPCOS): { factor2 = BETA_FACTOR; bc = -0.5*log(BETA_FACTOR); break; } case (TS_ASYM_EXPCOS): { factor2 = BETA_FACTOR; bc = -0.5*log(BETA_FACTOR); break; } case (TS_ATAN): { factor2 = (1.0/BETA_FACTOR - 1.0)*2.0/PI; break; } case (TS_TANH): { factor2 = 1.0/BETA_FACTOR - 1.0; break; } } bexp2 = -log(factor2)/(double)(FINAL_DECREASE_PHASE); first = 0; // printf("t1 = %i, factor2 = %.3lg\n", t1, factor2); } if (i < INITIAL_TIME) beta = BETA; else if (i < INITIAL_TIME + t1) { switch (BETA_SCHEDULE) { case (TS_EXPONENTIAL): { beta = BETA*exp(bexponent*(double)(i - INITIAL_TIME)); break; } case (TS_CYCLING): { beta = BETA*exp(bexponent*(double)(i - INITIAL_TIME)); beta = beta*2.0/(1.0 + cos(omega*(double)(i - INITIAL_TIME))); break; } case (TS_PERIODIC): { beta = BETA*exp(logf*sin(omega*(double)(i - INITIAL_TIME))); break; } case (TS_LINEAR): { beta = BETA/(1.0 + (1.0/BETA_FACTOR - 1.0)*(double)(i - INITIAL_TIME)/(double)(t1)); // printf("i = %i, beta = %.3lg\n", i, beta); break; } case (TS_COSINE): { beta = ac/(1.0 + bc*cos(omega*(double)(i - INITIAL_TIME))); printf("i = %i, beta = %.3lg\n", i, beta); break; } case (TS_EXPCOS): { beta = BETA*exp(bc*(-1.0 + cos(omega*(double)(i - INITIAL_TIME)))); // printf("i = %i, beta = %.3lg\n", i, beta); break; } case (TS_ASYM_EXPCOS): { t = (double)(i - INITIAL_TIME)/(double)(t1); beta = BETA*exp(bc*(-1.0 + cos(N_TOSCILLATIONS*DPI*(t - 0.5*t*(1.0-t))))); break; } case (TS_ATAN): { beta = BETA/(1.0 + factor2*atan(2.0*(double)(i - INITIAL_TIME)/(double)(t1))); break; } case (TS_TANH): { beta = BETA/(1.0 + factor2*tanh(TS_SLOPE*(double)(i - INITIAL_TIME)/(double)(t1))); break; } } } else if (i < INITIAL_TIME + t2) beta = BETA*factor2; else if (i < INITIAL_TIME + t3) { beta = BETA*exp(bexp2*(double)(i - INITIAL_TIME - t3)); } else beta = BETA; printf("beta = %.3lg\n", beta); return(beta); } double container_size_schedule(int i) { if ((i < INITIAL_TIME)||(i > INITIAL_TIME + NSTEPS - RESTORE_TIME)) return(INITXMIN); else return(INITXMIN + (1.0-COMPRESSION_RATIO)*(INITXMAX-INITXMIN)*(double)(i-INITIAL_TIME)/(double)(NSTEPS-RESTORE_TIME)); } double container_size_schedule_smooth(int i, int j) { double t; if ((i < INITIAL_TIME)||(i > INITIAL_TIME + NSTEPS - RESTORE_TIME)) return(INITXMIN); else { t = (double)(i-INITIAL_TIME) + (double)j/(double)NVID; return(INITXMIN + (1.0-COMPRESSION_RATIO)*(INITXMAX-INITXMIN)*t/(double)(NSTEPS-RESTORE_TIME)); } } double obstacle_schedule_old(int i) { double time; static double t1 = 0.5, t2 = 0.75, t3 = 0.875; if (i < INITIAL_TIME) return(OBSTACLE_XMIN); else { time = (double)(i-INITIAL_TIME)/(double)(NSTEPS); if (time < t1) return(OBSTACLE_XMIN + (OBSTACLE_XMAX - OBSTACLE_XMIN)*time/t1); else if (time < t2) return(OBSTACLE_XMIN + (OBSTACLE_XMAX - OBSTACLE_XMIN)*(time - t1)/(t2 - t1)); else if (time < t3) return(OBSTACLE_XMIN + (OBSTACLE_XMAX - OBSTACLE_XMIN)*(time - t2)/(t3 - t2)); else return(OBSTACLE_XMAX); } } double obstacle_schedule(int i) { double time, acceleration = 40.0; double x; // static double t1 = 0.5, t2 = 0.75, t3 = 0.875; if (i < INITIAL_TIME) return(OBSTACLE_XMIN); else { time = (double)(i-INITIAL_TIME)/(double)(NSTEPS); x = OBSTACLE_XMIN + 0.5*acceleration*time*time; xspeed = acceleration*time; // while (x > OBSXMAX) x += OBSXMIN - OBSXMAX; return(x); } } double obstacle_schedule_smooth(int i, int j) { double time, acceleration = 50.0; double x; // static double t1 = 0.5, t2 = 0.75, t3 = 0.875; if (i < INITIAL_TIME) return(OBSTACLE_XMIN); else { time = ((double)(i-INITIAL_TIME) + (double)j/(double)NVID)/(double)(NSTEPS); x = OBSTACLE_XMIN + 0.5*acceleration*time*time; xspeed = acceleration*time; // while (x > OBSXMAX) x += OBSXMIN - OBSXMAX; return(x); } } double gravity_schedule(int i, int j) { double time, gravity, x, y; switch (GRAVITY_SCHEDULE){ case (G_INCREASE_RELEASE): { if ((i < INITIAL_TIME + GRAVITY_INITIAL_TIME)||(i > NSTEPS + INITIAL_TIME - GRAVITY_RESTORE_TIME)) return(GRAVITY); else { time = ((double)(i - INITIAL_TIME - GRAVITY_INITIAL_TIME) + (double)j/(double)NVID)/(double)(NSTEPS - GRAVITY_RESTORE_TIME - GRAVITY_INITIAL_TIME); gravity = GRAVITY*(1.0 + time*(GRAVITY_FACTOR - 1.0)); return(gravity); } break; } case (G_INCREASE_DECREASE): { if ((i < INITIAL_TIME + GRAVITY_INITIAL_TIME)||(i > NSTEPS + INITIAL_TIME - GRAVITY_RESTORE_TIME)) return(GRAVITY); else { time = ((double)(i - INITIAL_TIME - GRAVITY_INITIAL_TIME) + (double)j/(double)NVID)/(double)(NSTEPS - GRAVITY_RESTORE_TIME - GRAVITY_INITIAL_TIME); x = 2.0 - cos(DPI*time); y = 0.5*((GRAVITY_FACTOR - 1.0)*x + 3.0 - GRAVITY_FACTOR); gravity = GRAVITY*y; return(gravity); } break; } } } double rotation_angle(double phase) { /* case of rotating hourglass */ // while (phase > DPI) phase -= DPI; // return(phase - 0.5*sin(2.0*phase)); /* case of centrifuge */ // while (phase > 1.0) phase -= 1.0; // phase *= DPI; // angular_speed = 0.5*OMEGAMAX*(1.0 - cos(phase)); // return(0.5*OMEGAMAX*(phase - sin(phase))); /* case of centrifuge remaining at constant speed for a while */ 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); } } } double rotation_schedule(int i) { double phase; static int imin = INITIAL_TIME + ROTATE_INITIAL_TIME, imax = INITIAL_TIME + NSTEPS - ROTATE_FINAL_TIME; if (i < imin) { angular_speed = 0.0; return(0.0); } else { if (i > imax) i = imax; phase = (DPI/(double)PERIOD_ROTATE_BOUNDARY)*(double)(i - imin); return(rotation_angle(phase)); } } double rotation_schedule_smooth(int i, int j) { double phase, angle, phase1, angle1; static int imin = INITIAL_TIME + ROTATE_INITIAL_TIME, imax = INITIAL_TIME + NSTEPS - ROTATE_FINAL_TIME; if (i < imin) { angular_speed = 0.0; return(0.0); } else { if (i > imax) { angle = rotation_angle(1.0); angular_speed = 0.0; } else { phase = (1.0/(double)(imax - imin))*((double)(i - imin) + (double)j/(double)NVID); angle = rotation_angle(phase); phase1 = (1.0/(double)(imax - imin))*((double)(i + 1 - imin) + (double)j/(double)NVID); angle1 = rotation_angle(phase1); angular_speed = 25.0*(angle1 - angle); } return(angle); } } int thermostat_schedule(int i) { if (i < INITIAL_TIME) return(1); 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 *ncoupled, int initial_phase) { double a, totalenergy = 0.0, damping, damping1, damping_rot1, direction, dmean, dratio; static double b = 0.25*SIGMA*SIGMA*DT_PARTICLE/MU_XI, xi = 0.0; int j, move, ncoup; if (initial_phase) damping = INITIAL_DAMPING; else damping = DAMPING; // printf("Evolving particles\n"); #pragma omp parallel for private(j,xi,totalenergy,a,move) for (j=0; j 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; qx[j] = particle[j].xc + 0.5*DT_PARTICLE*px[j]*particle[j].mass_inv; qy[j] = particle[j].yc + 0.5*DT_PARTICLE*py[j]*particle[j].mass_inv; qangle[j] = particle[j].angle + 0.5*DT_PARTICLE*pangle[j]*particle[j].inertia_moment_inv; if ((THERMOSTAT_ON)&&(particle[j].thermostat)) { px[j] *= exp(- 0.5*DT_PARTICLE*xi); py[j] *= exp(- 0.5*DT_PARTICLE*xi); } if ((COUPLE_ANGLE_TO_THERMOSTAT)&&(particle[j].thermostat)) pangle[j] *= exp(- 0.5*DT_PARTICLE*xi); } /* compute kinetic energy */ // *nactive = 0; ncoup = 1; for (j=0; j 0.0) { damping1 = damping*particle[j].damping; damping_rot1 = DAMPING_ROT*particle[j].damping; } if ((THERMOSTAT_ON)&&(particle[j].thermostat)) { px[j] *= exp(- 0.5*DT_PARTICLE*xi); py[j] *= exp(- 0.5*DT_PARTICLE*xi); if (!COUPLE_ANGLE_TO_THERMOSTAT) pangle[j] *= exp(- DT_PARTICLE*damping_rot1); } else { px[j] *= exp(- DT_PARTICLE*damping1); py[j] *= exp(- DT_PARTICLE*damping1); pangle[j] *= exp(- DT_PARTICLE*damping_rot1); // printf("Damping particle angular velocity\n"); } if ((THERMOSTAT_ON)&&(COUPLE_ANGLE_TO_THERMOSTAT)&&(particle[j].thermostat)) pangle[j] *= exp(- 0.5*DT_PARTICLE*xi); particle[j].xc = qx[j] + 0.5*DT_PARTICLE*px[j]*particle[j].mass_inv; particle[j].yc = qy[j] + 0.5*DT_PARTICLE*py[j]*particle[j].mass_inv; particle[j].angle = qangle[j] + 0.5*DT_PARTICLE*pangle[j]*particle[j].inertia_moment_inv; // particle[j].vx = px[j] + 0.5*DT_PARTICLE*particle[j].fx; // particle[j].vy = py[j] + 0.5*DT_PARTICLE*particle[j].fy; // particle[j].omega = pangle[j] + 0.5*DT_PARTICLE*particle[j].torque; /* TO DO: move this to function wrap_particle */ if ((BOUNDARY_COND == BC_PERIODIC_CIRCLE)||(BOUNDARY_COND == BC_PERIODIC_FUNNEL)||(BOUNDARY_COND == BC_PERIODIC_TRIANGLE)) { // if (particle[j].xc < BCXMIN) if (particle[j].xc < xshift + BCXMIN) { particle[j].xc += BCXMAX - BCXMIN; if (RESAMPLE_Y) { *nmove++; if (resample_particle(j, NTRIALS, particle) == 1) { px[j] = particle[j].vx; py[j] = particle[j].vy; update_hashgrid(particle, hashgrid, 0); *nsuccess++; } } } // else if (particle[j].xc > BCXMAX) else if (particle[j].xc > xshift + BCXMAX) { particle[j].xc += BCXMIN - BCXMAX; } if (particle[j].yc > BCYMAX) particle[j].yc += BCYMIN - BCYMAX; else if (particle[j].yc < BCYMIN) particle[j].yc += BCYMAX - BCYMIN; } else if (!NO_WRAP_BC) { move += wrap_particle(&particle[j], &px[j], &py[j]); } // if (move > 0) // { // compute_relative_positions(particle, hashgrid); // update_hashgrid(particle, hashgrid, 0); /* REDUNDANT ? */ // } } *ncoupled = ncoup; return(totalenergy); } double evolve_clusters(t_particle particle[NMAXCIRCLES], t_cluster cluster[NMAXCIRCLES], t_hashgrid hashgrid[HASHX*HASHY], double cqx[NMAXCIRCLES], double cqy[NMAXCIRCLES], double cqangle[NMAXCIRCLES], double cpx[NMAXCIRCLES], double cpy[NMAXCIRCLES], double cpangle[NMAXCIRCLES], double beta, int *nactive, int *nsuccess, int *nmove, int *ncoupled, int initial_phase, int verbose) { double a, totalenergy = 0.0, damping, direction, dmean, newx, newy, newangle, deltax, deltay, deltaangle; static double b = 0.25*SIGMA*SIGMA*DT_PARTICLE/MU_XI, xi = 0.0; int j, move, ncoup; if (initial_phase) damping = INITIAL_DAMPING; else damping = DAMPING; #pragma omp parallel for private(j,xi,totalenergy,a,move) for (j=0; j direction + PI) dmean -= DPI; cluster[j].dirmean = DRATIO*dmean + (1.0-DRATIO)*direction; if (cluster[j].dirmean < 0.0) cluster[j].dirmean += DPI; else if (cluster[j].dirmean > DPI) cluster[j].dirmean -= DPI; } if ((COUPLE_ANGLE_TO_THERMOSTAT)&&(cluster[j].thermostat)) cluster[j].energy += cpangle[j]*cpangle[j]*cluster[j].inertia_moment_inv; cqx[j] = cluster[j].xg + 0.5*DT_PARTICLE*cpx[j]*cluster[j].mass_inv; cqy[j] = cluster[j].yg + 0.5*DT_PARTICLE*cpy[j]*cluster[j].mass_inv; cqangle[j] = cluster[j].angle + 0.5*DT_PARTICLE*cpangle[j]*cluster[j].inertia_moment_inv; if ((THERMOSTAT_ON)&&(cluster[j].thermostat)) { cpx[j] *= exp(- 0.5*DT_PARTICLE*xi); cpy[j] *= exp(- 0.5*DT_PARTICLE*xi); } if ((COUPLE_ANGLE_TO_THERMOSTAT)&&(cluster[j].thermostat)) cpangle[j] *= exp(- 0.5*DT_PARTICLE*xi); } /* compute kinetic energy */ // *nactive = 0; ncoup = 1; for (j=0; j 0) // { // compute_relative_positions(particle, hashgrid); // update_hashgrid(particle, hashgrid, 0); /* REDUNDANT ? */ // } } // sleep(1); *ncoupled = ncoup; return(totalenergy); } void evolve_lid(double fboundary) { double force; force = fboundary - GRAVITY*LID_MASS; if (ylid > BCYMAX + LID_WIDTH) force -= KSPRING_BOUNDARY*(ylid - BCYMAX - LID_WIDTH); vylid += force*DT_PARTICLE/LID_MASS; ylid += vylid*DT_PARTICLE; } void evolve_wall(double fboundary) { double force; force = fboundary; if (xwall > BCYMAX - WALL_WIDTH) force -= KSPRING_BOUNDARY*(xwall - BCYMAX + WALL_WIDTH); else if (xwall < BCYMIN + WALL_WIDTH) force += KSPRING_BOUNDARY*(BCYMIN + WALL_WIDTH - xwall); force -= vxwall*WALL_FRICTION; vxwall += fboundary*DT_PARTICLE/WALL_MASS; if (vxwall > WALL_VMAX) vxwall = WALL_VMAX; else if (vxwall < -WALL_VMAX) vxwall = -WALL_VMAX; xwall += vxwall*DT_PARTICLE; // printf("fboundary = %.3lg, xwall = %.3lg, vxwall = %.3lg\n", fboundary, xwall, vxwall); } void evolve_segments(t_segment segment[NMAXSEGMENTS], int time) { int i, nactive = 0, group; double fx[2] = {0.0, 0.0}, fy[2] = {0.0, 0.0}, x, y, padding = 3.0*MU, mass2 = SEGMENTS_MASS; if (SEGMENT_PATTERN == S_TWO_CIRCLES_EXT) mass2 = SEGMENTS_MASS*TWO_CIRCLES_RADIUS_RATIO; for (group=0; group<2; group++) { fx[group] = 0.0; fy[group] = 0.0; } for (i=0; i XMAX - padding) fx[group] -= KSPRING_BOUNDARY*(x - XMAX + padding); if (y < YMIN + padding) fy[group] += KSPRING_BOUNDARY*(YMIN + padding - y); else if (y > YMAX - padding) fy[group] -= KSPRING_BOUNDARY*(y - YMAX + padding); } else if ((BOUNDARY_COND == BC_REFLECT_ABS)||(BOUNDARY_COND == BC_REFLECT_ABS_BOTTOM)) /* add force from simulation boundary */ { y = 0.5*(segment[i].y1 + segment[i].y2); if (y < YMIN) fy[group] += KSPRING_BOUNDARY*(YMIN - y); } if (group == 0) fy[group] -= GRAVITY*SEGMENTS_MASS; else fy[group] -= GRAVITY*mass2; } if (nactive > 0) for (group=0; group<2; group++) { fx[group] = fx[group]/(double)nactive; fy[group] = fy[group]/(double)nactive; } if (FLOOR_FORCE) { if (fx[0] > FMAX) fx[0] = FMAX; else if (fx[0] < -FMAX) fx[0] = -FMAX; if (fy[0] > FMAX) fy[0] = FMAX; else if (fy[0] < -FMAX) fy[0] = -FMAX; } vxsegments[0] += fx[0]*DT_PARTICLE/SEGMENTS_MASS; vysegments[0] += fy[0]*DT_PARTICLE/SEGMENTS_MASS; xsegments[0] += vxsegments[0]*DT_PARTICLE; ysegments[0] += vysegments[0]*DT_PARTICLE; if (TWO_OBSTACLES) { if (FLOOR_FORCE) { if (fx[1] > FMAX) fx[1] = FMAX; else if (fx[1] < -FMAX) fx[1] = -FMAX; if (fy[1] > FMAX) fy[1] = FMAX; else if (fy[1] < -FMAX) fy[1] = -FMAX; } vxsegments[1] += fx[1]*DT_PARTICLE/mass2; vysegments[1] += fy[1]*DT_PARTICLE/mass2; xsegments[1] += vxsegments[1]*DT_PARTICLE; ysegments[1] += vysegments[1]*DT_PARTICLE; } /* add some damping if y coordinate is small (for lunar landing) */ if (DAMP_SEGS_AT_NEGATIVE_Y) for (group=0; group<2; group++) if (ysegments[group] < 0.1) { vysegments[group] *= exp(-DAMPING*DT_PARTICLE); vxsegments[group] *= exp(-DAMPING*DT_PARTICLE); } /* to avoid numerical instabilities */ for (group=0; group<2; group++) { if (xsegments[group] + 1.0 > BCXMAX) { xsegments[group] = BCXMAX - 1.0; vxsegments[group] = 0.0; } if ((RELEASE_ROCKET_AT_DEACTIVATION)&&((BOUNDARY_COND == BC_REFLECT_ABS)||(BOUNDARY_COND == BC_ABSORBING))) { // ysegments[group] = SEGMENTS_Y0; if (time < SEGMENT_DEACTIVATION_TIME) vysegments[group] = 0.0; else if ((ysegments[group] < SEGMENTS_Y0)&&(vysegments[group] < 0.0)) vysegments[group] = -0.5*vysegments[group]; } } translate_segments(segment, xsegments, ysegments); } void evolve_segment_groups(t_segment segment[NMAXSEGMENTS], int time, t_group_segments segment_group[NMAXGROUPS]) /* new version of evolve_segments that takes the group structure into account */ { double fx[NMAXGROUPS], fy[NMAXGROUPS], torque[NMAXGROUPS], dx[NMAXGROUPS], dy[NMAXGROUPS], dalpha[NMAXGROUPS]; double x, y, dx0, dy0, padding, proj, distance, f, xx[2], yy[2], xmean = 0.0, ymean = 0.0, ymax = 0.0; int i, j, k, group = 0; static double maxdepth, saturation_depth, xmax; static int first = 1; if (first) { xmax = XMAX - TRACK_X_PADDING; if ((PLOT_SPEEDS)||(PLOT_TRAJECTORIES)) xmax -= 1.8; first = 0; } maxdepth = 0.5*GROUP_WIDTH; saturation_depth = 0.1*GROUP_WIDTH; padding = 0.1; for (group=0; group 0)) { group = segment[i].group; fx[group] += segment[i].fx; fy[group] += segment[i].fy; torque[group] += segment[i].torque; dx0 = segment[i].xc - segment_group[group].xc; dy0 = segment[i].yc - segment_group[group].yc; torque[group] += dx0*segment[i].fy - dy0*segment[i].fx; if (BOUNDARY_COND == BC_SCREEN) /* add force from simulation boundary */ { x = 0.5*(segment[i].x1 + segment[i].x2); y = 0.5*(segment[i].y1 + segment[i].y2); if (x < XMIN + padding) fx[group] += KSPRING_BOUNDARY*(XMIN + padding - x); else if (x > XMAX - padding) fx[group] -= KSPRING_BOUNDARY*(x - XMAX + padding); if (y < YMIN + padding) fy[group] += KSPRING_BOUNDARY*(YMIN + padding - y); else if (y > YMAX - padding) fy[group] -= KSPRING_BOUNDARY*(y - YMAX + padding); } else if ((BOUNDARY_COND == BC_REFLECT_ABS)||(BOUNDARY_COND == BC_REFLECT_ABS_BOTTOM)) /* add force from simulation boundary */ { y = 0.5*(segment[i].y1 + segment[i].y2); if (y < YMIN) fy[group] += KSPRING_BOUNDARY*(YMIN - y); } /* repulsion between different groups */ if (GROUP_REPULSION) for (j=0; j 0.0)&&(proj < 1.0)) { distance = segment[i].nx*x + segment[i].ny*y - segment[i].c; if ((distance > -maxdepth)&&(distance < 0.0)) { if (distance < -saturation_depth) distance = -saturation_depth; f = KSPRING_GROUPS*(-distance); segment[j].fx += f*segment[i].nx; segment[j].fy += f*segment[i].ny; segment[j].torque += (x - segment[i].xc)*f*segment[i].ny - (y - segment[i].yc)*f*segment[i].nx; fx[group] -= f*segment[i].nx; fy[group] -= f*segment[i].ny; torque[group] -= (x - segment[i].xc)*f*segment[i].ny - (y - segment[i].yc)*f*segment[i].nx; } } } } } if (GROUP_G_REPEL) for (i=0; i FMAX) fx[group] = FMAX; else if (fx[group] < -FMAX) fx[group] = -FMAX; if (fy[group] > FMAX) fy[group] = FMAX; else if (fy[group] < -FMAX) fy[group] = -FMAX; } for (group=1; group 0)) { group = segment[i].group; translate_one_segment(segment, i, dx[group], dy[group]); rotate_one_segment(segment, i, dalpha[group], segment_group[group].xc, segment_group[group].yc); } if (TRACK_SEGMENT_GROUPS) { /* compute mean position */ for (group=1; group ymax) ymax = segment_group[group].yc; } xmean = xmean/((double)(ngroups-1)); ymean = ymean/((double)(ngroups-1)); /* bias towards ymax */ ymean = 0.75*ymax + 0.25*ymean; if (ymean > ytrack) ytrack = ymean; if (xmean > xmax) xtrack = xmean - xmax; else if (xmean < XMIN + TRACK_X_PADDING) xtrack = xmean - XMIN - TRACK_X_PADDING; } } void animation() { double time, scale, diss, rgb[3], dissip, gradient[2], x, y, dx, dy, dt, xleft, xright, a, b, length, fx, fy, force[2], totalenergy = 0.0, pos[2], prop, vx, xi = 0.0, torque, torque_ij, pleft = 0.0, pright = 0.0, entropy[2], speed_ratio, xmin, xmax, ymin, ymax, delta_energy, speed, ratio = 1.0, ratioc, cum_etot = 0.0, emean = 0.0, radius_ratio, t, angle, theta, sum, alpha, bfield; double *qx, *qy, *px, *py, *qangle, *pangle, *pressure, *obstacle_speeds; double *cqx, *cqy, *cpx, *cpy, *cqangle, *cpangle; int i, j, k, n, m, s, ij[2], i0, iplus, iminus, j0, jplus, jminus, p, q, p1, q1, p2, q2, total_neighbours = 0, cl, 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, ncoupled = 1, np, belt, obs; int *particle_numbers; static int imin, imax; static short int first = 1; t_particle *particle; t_cluster *cluster; t_obstacle *obstacle; t_segment *segment; t_group_segments *segment_group; t_tracer *trajectory; t_group_data *group_speeds; t_collision *collisions; t_hashgrid *hashgrid; t_molecule *molecule; t_lj_parameters params; t_belt *conveyor_belt; 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 (CLUSTER_PARTICLES) cluster = (t_cluster *)malloc(NMAXCIRCLES*sizeof(t_cluster)); /* particle clusters */ if (ADD_FIXED_OBSTACLES) obstacle = (t_obstacle *)malloc(NMAXOBSTACLES*sizeof(t_obstacle)); /* circular obstacles */ if (ADD_FIXED_SEGMENTS) { segment = (t_segment *)malloc(NMAXSEGMENTS*sizeof(t_segment)); /* linear obstacles */ segment_group = (t_group_segments *)malloc(NMAXGROUPS*sizeof(t_group_segments)); conveyor_belt = (t_belt *)malloc(NMAXBELTS*sizeof(t_belt)); } if (TRACER_PARTICLE) trajectory = (t_tracer *)malloc(TRAJECTORY_LENGTH*N_TRACER_PARTICLES*sizeof(t_tracer)); if (PAIR_PARTICLES) molecule = (t_molecule *)malloc(NMAXCIRCLES*sizeof(t_molecule)); /* molecules */ hashgrid = (t_hashgrid *)malloc(HASHX*HASHY*sizeof(t_hashgrid)); /* hashgrid */ qx = (double *)malloc(NMAXCIRCLES*sizeof(double)); qy = (double *)malloc(NMAXCIRCLES*sizeof(double)); px = (double *)malloc(NMAXCIRCLES*sizeof(double)); py = (double *)malloc(NMAXCIRCLES*sizeof(double)); qangle = (double *)malloc(NMAXCIRCLES*sizeof(double)); pangle = (double *)malloc(NMAXCIRCLES*sizeof(double)); pressure = (double *)malloc(N_PRESSURES*sizeof(double)); if (CLUSTER_PARTICLES) { cqx = (double *)malloc(NMAXCIRCLES*sizeof(double)); cqy = (double *)malloc(NMAXCIRCLES*sizeof(double)); cpx = (double *)malloc(NMAXCIRCLES*sizeof(double)); cpy = (double *)malloc(NMAXCIRCLES*sizeof(double)); cqangle = (double *)malloc(NMAXCIRCLES*sizeof(double)); cpangle = (double *)malloc(NMAXCIRCLES*sizeof(double)); } if (REACTION_DIFFUSION) { collisions = (t_collision *)malloc(2*NMAXCOLLISIONS*sizeof(t_collision)); for (i=0; i<2*NMAXCOLLISIONS; i++) collisions[i].time = 0; } if (SAVE_TIME_SERIES) { lj_time_series = fopen("lj_time_series.dat", "w"); lj_final_position = fopen("lj_final_position.dat", "w"); } lj_log = fopen("lj_logfile.txt", "w"); if (ADD_FIXED_OBSTACLES) init_obstacle_config(obstacle); if (ADD_FIXED_SEGMENTS) init_segment_config(segment, conveyor_belt); if ((MOVE_SEGMENT_GROUPS)&&(ADD_FIXED_SEGMENTS)) { for (i=0; i= 2000) exit(0); frame_time++; printf("Computing frame %d\n",i); if (INCREASE_KREPEL) params.krepel = repel_schedule(i); if (INCREASE_BETA) params.beta = temperature_schedule(i); if (INCREASE_E) params.efield = efield_schedule(i); else params.efield = EFIELD; if (INCREASE_B) params.bfield = bfield_schedule(i); else params.bfield = BFIELD; if ((PARTIAL_THERMO_COUPLING)&&(PARTIAL_THERMO_REGION == TH_RING_EXPAND)) params.thermo_radius = PARTIAL_THERMO_RIN + (double)i/(double)NSTEPS*(PARTIAL_THERMO_RFIN - PARTIAL_THERMO_RIN); if (DECREASE_CONTAINER_SIZE) { 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) { thermostat_on = thermostat_schedule(i); printf("Termostat: %i\n", thermostat_on); } /* deactivate some segments */ if ((ADD_FIXED_SEGMENTS)&&(DEACTIVATE_SEGMENT)&&(i == INITIAL_TIME + SEGMENT_DEACTIVATION_TIME + 1)) for (j=0; j xmax) xmax = particle[j].xc; if (particle[j].yc < ymin) ymin = particle[j].yc; if (particle[j].yc > ymax) ymax = particle[j].yc; } for (j=0; j 0.0)&&(speed < VICSEK_VMIN)) speed = VICSEK_VMIN; if ((INTERACTION == I_VICSEK_SHARK)&&(particle[j].type == 2)) speed *= 1.75; if (speed > VICSEK_VMAX) speed = 0.5*(speed + VICSEK_VMAX); particle[j].vx = speed*cos(particle[j].angle); particle[j].vy = speed*sin(particle[j].angle); speed = module2(px[j],py[j]); if ((VICSEK_VMIN > 0.0)&&(speed < VICSEK_VMIN)) speed = VICSEK_VMIN; if ((INTERACTION == I_VICSEK_SHARK)&&(particle[j].type == 2)) speed *= 1.75; if (speed > VICSEK_VMAX) speed = 0.5*(speed + VICSEK_VMAX); px[j] = speed*cos(particle[j].angle); py[j] = speed*sin(particle[j].angle); } /* TEST - adjust angles */ if ((REACTION_DIFFUSION)&&(RD_REACTION == CHEM_AGGREGATION_NNEIGH)) { // if (particle[j].npartners >= 1) // { // angle = 0.0; // theta = 0.99; // for (p=0; p= 1) { x = 0.0; y = 0.0; for (p=0; p= 2) // { // angle = particle[j].angle; // while (angle > alpha) angle -= alpha; // sum = angle; // for (p=0; p alpha) angle -= alpha; // sum += angle; // } // sum *= 1.0/(double)(np+1); // particle[j].angle = sum; // for (p=0; p WIND_YMIN)) { particle[j].fx += WIND_FORCE*particle[j].radius*particle[j].mass_inv; } if (FLOOR_FORCE) { if (particle[j].fx > FMAX) particle[j].fx = FMAX; if (particle[j].fx < -FMAX) particle[j].fx = -FMAX; if (particle[j].fy > FMAX) particle[j].fy = FMAX; if (particle[j].fy < -FMAX) particle[j].fy = -FMAX; if (particle[j].torque > FMAX) particle[j].torque = FMAX; if (particle[j].torque < -FMAX) particle[j].torque = -FMAX; } } /* compute force and torque on clusters */ if (CLUSTER_PARTICLES) compute_cluster_force(cluster, particle); /* timestep of thermostat algorithm */ if (CLUSTER_PARTICLES) { totalenergy = evolve_clusters(particle, cluster, hashgrid, cqx, cqy, cqangle, cpx, cpy, cpangle, params.beta, ¶ms.nactive, &nsuccess, &nmove, &ncoupled, i < INITIAL_TIME, n == 0); /* FIXME: update particle data */ for (j=0; j= 2)) repair_cluster(cl, particle, cluster, 1000/NVID, 1); /* evolution of lid coordinate */ if (BOUNDARY_COND == BC_RECTANGLE_LID) evolve_lid(params.fboundary); if (BOUNDARY_COND == BC_RECTANGLE_WALL) { 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); if (ADD_CONVEYOR_FORCE) for (belt = 0; belt < nbelts; belt++) conveyor_belt[belt].position += conveyor_belt[belt].speed*DT_PARTICLE; if (MOVE_CONVEYOR_BELT) update_conveyor_belts(segment, conveyor_belt); if (ROTATE_OBSTACLES) for (obs = 0; obs < nobstacles; obs++) { /* TEST */ obstacle[obs].omega = obstacle[obs].omega0*(1.0 + 0.5*cos(i*DPI/200.0)); obstacle[obs].angle -= obstacle[obs].omega*DT_PARTICLE; if ((RATTLE_OBSTACLES)&&(obstacle[obs].oscillate)) { theta = obstacle[obs].phase + DPI*(double)i/(double)obstacle[obs].period; obstacle[obs].xc = obstacle[obs].xc0 + obstacle[obs].amplitude*cos(theta); obstacle[obs].yc = obstacle[obs].yc0 + obstacle[obs].amplitude*sin(theta); } } if ((MOVE_SEGMENT_GROUPS)&&(i > INITIAL_TIME + SEGMENT_DEACTIVATION_TIME)) evolve_segment_groups(segment, i, segment_group); // if ((MOVE_SEGMENT_GROUPS)&&(i > OBSTACLE_INITIAL_TIME)) evolve_segment_groups(segment, i, segment_group); } /* end of for (n=0; nINITIAL_TIME)&&(SAVE_TIME_SERIES)) { n_total_active = 0; for (j=0; jN_T_AVERAGE)) { nthermo = partial_thermostat_coupling(particle, xshift + PARTIAL_THERMO_SHIFT, segment, params); printf("%i particles coupled to thermostat out of %i active\n", nthermo, params.nactive); params.mean_energy = compute_mean_energy(particle); } else params.mean_energy = totalenergy/(double)ncircles; if (CENTER_PX) center_momentum(px); if (CENTER_PY) center_momentum(py); if (CENTER_PANGLE) center_momentum(pangle); if (FLOOR_OMEGA) floor = floor_momentum(pangle); // printf("pressure left %.5lg, pressure right %.5lg\n", pleft, pright); // for (j=0; j DPI) particle[j].angle -= DPI; while (particle[j].angle < 0.0) particle[j].angle += DPI; } /* update tracer particle trajectory */ if ((TRACER_PARTICLE)&&(i > INITIAL_TIME)) { for (j=0; j= TRAJECTORY_LENGTH) traj_position = 0; traj_length++; if (traj_length >= TRAJECTORY_LENGTH) traj_length = TRAJECTORY_LENGTH - 1; // for (j=0; j 0.0) { emean = cum_etot/(double)(i+1); if (totalenergy > 10.0*emean) { reset_energy(particle, px, py, totalenergy, emean); totalenergy = 0.0; for (j=0; j max_nb) max_nb = particle[j].neighb; if (particle[j].neighb < min_nb) min_nb = particle[j].neighb; } // printf("Mean number of neighbours: %.3f\n", (double)total_neighbours/(double)ncircles); // printf("Min number of neighbours: %i\n", min_nb); // printf("Max number of neighbours: %i\n", max_nb); blank(); /* case of reaction-diffusion equation */ if ((i > INITIAL_TIME)&&(REACTION_DIFFUSION)) { ncollisions = update_types(particle, molecule, cluster, collisions, ncollisions, particle_numbers, i - INITIAL_TIME - 1, &delta_energy); if (EXOTHERMIC) params.beta *= 1.0/(1.0 + delta_energy/totalenergy); params.nactive = 0; for (j=0; j ADD_TIME)&&((i - INITIAL_TIME - ADD_TIME)%ADD_PERIOD == 1)&&(i < NSTEPS - FINAL_NOADD_PERIOD)) { /* add enzymes */ if ((REACTION_DIFFUSION)&&((RD_REACTION == CHEM_DNA_ENZYME)||(RD_REACTION == CHEM_DNA_ENZYME_REPAIR))) { for (k=0; k INITIAL_TIME + WALL_TIME)&&(PRINT_ENTROPY)) { compute_entropy(particle, entropy); printf("Entropy 1 = %.5lg, Entropy 2 = %.5lg\n", entropy[0], entropy[1]); print_entropy(entropy); } /* these should be moved to draw_frame */ if (PRINT_SEGMENTS_SPEEDS) { if (MOVE_BOUNDARY) print_segments_speeds(vxsegments, vysegments); else print_segment_group_speeds(segment_group); } if ((i > INITIAL_TIME)&&(PLOT_PARTICLE_NUMBER)) count_particle_number(particle, particle_numbers, i - INITIAL_TIME); draw_frame(i, PLOT, BG_COLOR, ncollisions, traj_position, traj_length, wall, pressure, pleft, pright, particle_numbers, 1, params, particle, cluster, collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n); if (!((NO_EXTRA_BUFFER_SWAP)&&(MOVIE))) glutSwapBuffers(); if (MOVIE) { if (i >= INITIAL_TIME) { if ((TIME_LAPSE)&&(TIME_LAPSE_FIRST)) { if ((TIME_LAPSE)&&((i - INITIAL_TIME)%TIME_LAPSE_FACTOR == 0)&&(!DOUBLE_MOVIE)) { save_frame_lj(); } save_frame_lj_counter(NSTEPS/TIME_LAPSE_FACTOR + MID_FRAMES + i - INITIAL_TIME); } else { save_frame_lj(); if ((TIME_LAPSE)&&((i - INITIAL_TIME)%TIME_LAPSE_FACTOR == 0)&&(!DOUBLE_MOVIE)) { save_frame_lj_counter(NSTEPS + END_FRAMES + (i - INITIAL_TIME)/TIME_LAPSE_FACTOR); } } } else printf("Initial phase time %i of %i\n", i, INITIAL_TIME); if ((i >= INITIAL_TIME)&&(DOUBLE_MOVIE)) { draw_frame(i, PLOT_B, BG_COLOR_B, ncollisions, traj_position, traj_length, wall, pressure, pleft, pright, particle_numbers, 0, params, particle, cluster, collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n); glutSwapBuffers(); save_frame_lj_counter(NSTEPS + MID_FRAMES + 1 + counter); 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 lj*.tif tif_ljones/"); } } fclose(lj_log); lj_log = fopen("lj_logfile.txt", "a"); } if (SAVE_TIME_SERIES) { n_total_active = 0; for (j=0; j