diff --git a/global_3d.c b/global_3d.c index c4e0bbd..51109a6 100644 --- a/global_3d.c +++ b/global_3d.c @@ -37,6 +37,7 @@ #define GF_AIRFOIL 3 /* curved repelling ellipse */ #define GF_WING 4 /* wing shape */ #define GF_COMPUTE_FROM_BC 5 /* compute force field as gradient of bc_field2 */ +#define GF_EARTH 6 /* field depends on altitude on continents */ /* Choice of water depth for shallow water equation */ @@ -82,7 +83,10 @@ #define PLANET ((B_DOMAIN == D_SPHERE_EARTH)||(B_DOMAIN == D_SPHERE_MARS)||(B_DOMAIN == D_SPHERE_MOON)||(B_DOMAIN == D_SPHERE_VENUS)||(B_DOMAIN == D_SPHERE_MERCURY)) #define OTHER_PLANET ((B_DOMAIN == D_SPHERE_MARS)||(B_DOMAIN == D_SPHERE_MOON)||(B_DOMAIN == D_SPHERE_VENUS)||(B_DOMAIN == D_SPHERE_MERCURY)) +#define RDE_PLANET ((ADAPT_STATE_TO_BC)&&(OBSTACLE_GEOMETRY == D_SPHERE_EARTH)) + #define NMAXCIRC_SPHERE 100 /* max number of circles on sphere */ +#define NMAX_TRACER_PTS 20 /* max number of tracer points recorded per cell */ int global_time = 0; double max_depth = 1.0; @@ -132,6 +136,10 @@ typedef struct double depth; /* water depth */ double cos_depth_angle; /* cos of angle of depth profile */ double gradx, grady; /* gradient of water depth */ + short int tracer; /* has value 1 if cell contains a tracer */ + short int n_tracer_pts; /* number of recorded tracer points per cell */ + double tracerx[NMAX_TRACER_PTS], tracery[NMAX_TRACER_PTS]; /* coordinates of tracer */ + int prev_cell; /* cell where tracer was previously */ } t_rde; @@ -147,6 +155,7 @@ typedef struct double phi, theta; /* phi, theta angles */ double cphi, sphi; /* cos and sin of phi */ double ctheta, stheta, cottheta; /* cos, sin and cotangent of theta */ + double reg_cottheta; /* regularized cotangent of theta */ double x, y, z; /* x, y, z coordinates of point on sphere */ double radius; /* radius with wave height */ double radius_dem; /* radius with digital elevation model */ diff --git a/global_ljones.c b/global_ljones.c index 9fa3d4d..e376dfb 100644 --- a/global_ljones.c +++ b/global_ljones.c @@ -18,6 +18,7 @@ #define NMAXGROUPS 50 /* max number of groups of segments */ #define NMAXCOLLISIONS 200000 /* max number of collisions */ #define NMAXPARTNERS 30 /* max number of partners in molecule */ +#define NMAXPARTNERMOLECULES 10 /* max number of partners of a molecule */ #define C_SQUARE 0 /* square grid of circles */ #define C_HEX 1 /* hexagonal/triangular grid of circles */ @@ -48,6 +49,8 @@ #define O_CIRCLE 4 /* one circle at the origin */ #define O_FOUR_CIRCLES 5 /* four circles */ #define O_HEX 6 /* hexagonal lattice */ +#define O_SIDES 7 /* grid along the sides of the simulation rectangle */ +#define O_SIDES_B 71 /* finer grid along the sides of the simulation rectangle */ /* pattern of additional repelling segments */ #define S_RECTANGLE 0 /* segments forming a rectangle */ @@ -101,6 +104,10 @@ #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 */ +#define I_COULOMB_PENTA 13 /* Lennard-Jones force with or without pentagonal symmetry depending on charge */ +#define I_COULOMB_IMAGINARY 14 /* Coulomb interaction with "imaginary charge" */ +#define I_DNA_CHARGED 15 /* Coulomb-type interaction between end points of DNA nucleotides */ +#define I_DNA_CHARGED_B 151 /* stronger Coulomb-type interaction between end points of DNA nucleotides */ /* Boundary conditions */ @@ -144,6 +151,8 @@ #define TS_COSINE 4 /* periodic time dependence, cosine */ #define TS_EXPCOS 5 /* periodic time dependence, exponential of cosine */ #define TS_ASYM_EXPCOS 6 /* periodic time dependence, asymmetric exponential of cosine */ +#define TS_ATAN 7 /* atan approaching asymptotic value */ +#define TS_TANH 8 /* tanh approaching asymptotic value */ /* Gravity schedules */ @@ -191,6 +200,14 @@ #define CHEM_2H2O_H3O_OH 21 /* 2 H2O <-> H3O+ + OH- */ #define CHEM_AGGREGATION 22 /* agregation of molecules coming close */ #define CHEM_AGGREGATION_CHARGE 23 /* agregation of charged molecules coming close */ +#define CHEM_AGGREGATION_NNEIGH 24 /* agregation of molecules with limitation on neighbours */ +#define CHEM_DNA 25 /* aggregation of DNA molecules */ +#define CHEM_DNA_ALT 251 /* aggregation of DNA molecules with constraints on connections */ +#define CHEM_DNA_DOUBLE 252 /* aggregation of DNA molecules with different ends */ +#define CHEM_DNA_DSPLIT 253 /* aggregation/splitting of DNA molecules with different ends */ +#define CHEM_DNA_BASE_SPLIT 254 /* aggregation/splitting of DNA molecules when base pairs don't match */ +#define CHEM_DNA_ENZYME 255 /* aggregation/splitting of DNA molecules in presence of enzymes */ +#define CHEM_DNA_ENZYME_REPAIR 256 /* aggregation/splitting of DNA molecules in presence of enzymes and additional repairing of bad connections */ /* Initial conditions for chemical reactions */ @@ -206,6 +223,8 @@ #define IC_SIGNX 7 /* type 1 or 2 depending on sign of x */ #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 */ +#define IC_DNA_POLYMERASE 10 /* initial condition for DNA polymerase */ +#define IC_DNA_POLYMERASE_REC 11 /* initial condition for DNA polymerase with recombination */ /* Initial conditions for option TWO_TYPES */ @@ -259,10 +278,15 @@ #define POLY_WATER 2 /* star-shaped with a 120° separation between anions */ #define POLY_SOAP 3 /* polymers with all-to-all coupling and polar end */ #define POLY_SOAP_B 4 /* polymers with pairwise coupling and polar end */ +#define POLY_SOAP_N 41 /* polymers with pairwise coupling and neutral polar end */ +#define POLY_SOAP_NMIX 42 /* polymers mixing neutral polar and neutral end */ #define POLY_PLUSMINUS 5 /* polymers with ends of opposite charge */ #define POLY_HYDRA 6 /* star-shaped with longer arms */ #define POLY_HYDRA_RIGID 61 /* star-shaped with longer arms and rigid first ring */ -// #define POLY_GLUE 99 /* dummy value for option CHEM_AGGREGATION */ +#define POLY_DNA 7 /* simplified model for DNA */ +#define POLY_DNA_ALT 71 /* simplified model for DNA with different short ends */ +#define POLY_DNA_DOUBLE 72 /* simplified model for DNA with double ends for rigidity */ +#define POLY_DNA_FLEX 73 /* simplified model for DNA with less backbone rigidity (beta) */ /* Background color schemes */ @@ -342,8 +366,15 @@ typedef struct short int npartners; /* number of partner particles */ double partner_eqd[NMAXPARTNERS]; /* equilibrium distances between partners */ int p0, p1; /* numbers of two first partners (for P_MOL_ANGLE color scheme) */ +// short int mol_angle; /* for color scheme P_MOL_ANGLE */ int cluster; /* number of cluster */ - short int tested, cactive; /* for cluster search */ + int molecule; /* number of molecule */ + short int tested, cactive; /* for cluster search */ + short int coulomb; /* has value 1 if DNA-Coulomb interaction is attractive */ + short int added; /* has value 1 if particle has been added */ + short int reactive; /* has value 1 if particle can react */ + short int paired; /* has value 1 if belongs to base-paired molecule */ + int partner_molecule; /* number of partner molecule */ } t_particle; typedef struct @@ -434,6 +465,15 @@ typedef struct int color; /* color hue in case of different collisions */ } t_collision; +typedef struct +{ + int nparticles; /* number of particles */ + int particle[NPARTNERS+1]; /* list of particles */ + int npartners; /* number of partner molecules */ + int partner[NMAXPARTNERMOLECULES]; /* list of partner molecules */ + int connection_type[NMAXPARTNERMOLECULES]; /* types of particles in connection */ + short int added; /* has value 1 if molecule has been added */ +} t_molecule; typedef struct { @@ -456,5 +496,6 @@ typedef struct -int ncircles, nobstacles, nsegments, ngroups = 1, counter = 0; +int frame_time = 0, ncircles, nobstacles, nsegments, ngroups = 1, counter = 0, nmolecules = 0; +FILE *lj_log; diff --git a/global_pdes.c b/global_pdes.c index c803139..1286202 100644 --- a/global_pdes.c +++ b/global_pdes.c @@ -100,7 +100,9 @@ #define D_TREE 73 /* Christmas tree, to use with IOR_TREE */ #define D_MICHELSON 74 /* Michelson interferometer, to use with IOR_MICHELSON */ #define D_MICHELSON_MOVING 741 /* moving Michelson interferometer, to use with IOR_MICHELSON */ -#define D_RITCHEY_CHRETIEN 75 /* Ritchey-Chrétien telescope */ +#define D_RITCHEY_CHRETIEN_SPHERICAL 75 /* Ritchey-Chrétien telescope with spherical mirrors */ +#define D_RITCHEY_CHRETIEN_HYPERBOLIC 751 /* Ritchey-Chrétien telescope with hyperbolic mirrors */ +#define D_GRADIENT_INDEX_LENS 76 /* gradient index lens (only affects draw_billiard) */ /* for wave_sphere.c */ @@ -120,6 +122,7 @@ #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) */ +#define NMAXSOURCES 10 /* maximal number of sources */ #define C_SQUARE 0 /* square grid of circles */ #define C_HEX 1 /* hexagonal/triangular grid of circles */ @@ -185,6 +188,10 @@ #define IOR_WAVE_GUIDES_COUPLED_B 151 /* coupled wave guides, variant where only corners are reflecting */ #define IOR_WAVE_GUIDE_COATED 16 /* short coated S-shaped optical fiber, to use with D_WAVEGUIDE_S_SHORT */ #define IOR_MICHELSON 17 /* Michelson interferometer, to use with D_MICHELSON */ +#define IOR_GRADIENT_INDEX_LENS 18 /* gradient index lens (parabolic c(r)^2) */ +#define IOR_GRADIENT_INDEX_LENS_B 181 /* gradient index lens (parabolic c(r)) */ +#define IOR_LINEAR_X_A 19 /* IoR depending linearly on x */ +#define IOR_LINEAR_X_B 191 /* IoR depending linearly on x, with correct boundary values */ #define IOR_EARTH_DEM 20 /* digital elevation model (for waves on sphere) */ @@ -204,11 +211,14 @@ #define OSC_PERIODIC 0 /* periodic oscillation */ #define OSC_SLOWING 1 /* oscillation of slowing frequency (anti-chirp) */ #define OSC_WAVE_PACKET 2 /* Gaussian wave packet */ +#define OSC_WAVE_PACKETS 21 /* Gaussian wave packets */ #define OSC_CHIRP 3 /* chirp (linearly accelerating frequency) */ #define OSC_BEAM 4 /* periodic oscillation modulated by y cut-off */ #define OSC_BEAM_GAUSSIAN 41 /* periodic oscillation modulated by Gaussian in y */ #define OSC_BEAM_SINE 42 /* periodic oscillation modulated by sine in y */ #define OSC_BEAM_TWOPERIODS 5 /* sum of two periodic oscillations modulated by y cut-off */ +#define OSC_TWO_WAVES 6 /* two linear waves at an angle, separate */ +#define OSC_TWO_WAVES_ADDED 61 /* two linear waves at an angle, superimposed */ /* Wave packet types */ @@ -305,6 +315,8 @@ #define Z_NORM_GRADIENT_INTENSITY 24 /* gradient and intensity of polar angle */ #define Z_VORTICITY 25 /* curl of polar angle */ #define Z_VORTICITY_ABS 251 /* absolute value of curl of polar angle */ +#define Z_MAXTYPE_RPS 26 /* color or type with maximal density */ +// #define Z_ZERO 99 /* return zero */ /* for Schrodinger equation */ #define Z_MODULE 30 /* module squared of first two fields */ @@ -315,6 +327,7 @@ #define Z_MAXTYPE_RPSLZ 40 /* color of type with maximal density */ #define Z_THETA_RPSLZ 41 /* polar angle */ #define Z_NORM_GRADIENT_RPSLZ 42 /* gradient of polar angle */ +#define Z_NORM_GRADIENT_RPSLZ_ASYM 43 /* gradient of polar angle */ /* for Euler incompressible Euler equation */ #define Z_EULER_VORTICITY 50 /* vorticity of velocity */ @@ -415,7 +428,6 @@ typedef struct } t_wave_source; - int ncircles = NMAXCIRCLES; /* actual number of circles, can be decreased e.g. for random patterns */ int npolyline = NMAXPOLY; /* actual length of polyline */ int npolyrect = NMAXPOLY; /* actual number of polyrect */ @@ -438,8 +450,11 @@ t_circle circles_b[NMAXCIRCLES]; /* circular scatterers */ t_polygon polygons_b[NMAXCIRCLES]; /* polygonal scatterers */ t_vertex polyline_b[NMAXPOLY]; /* vertices of polygonal line */ +double courant2, courantb2; /* Courant parameters squared */ + // double julia_x = -0.5, julia_y = 0.5; /* parameters for Julia sets */ // double julia_x = 0.33267, julia_y = 0.06395; /* parameters for Julia sets */ double julia_x = 0.37468, julia_y = 0.21115; /* parameters for Julia sets */ -double wave_source_x, wave_source_y; /* position of wave source */ +double wave_source_x[NMAXSOURCES], wave_source_y[NMAXSOURCES]; /* position of wave source */ +double focus_x = XMAX; /* position of focus */ double michelson_position = 0.0; /* position of mirror in Michelson interferometer */ diff --git a/heat.c b/heat.c index f93ee6f..383b4ce 100644 --- a/heat.c +++ b/heat.c @@ -227,6 +227,15 @@ #define MESSAGE_LINTERLETTER 60 /* length of interval between letters for Morse code message */ #define MESSAGE_LSPACE 48 /* length of space for Morse code message */ #define MESSAGE_INITIAL_TIME 100 /* initial time before starting message for Morse code message */ +#define WAVE_PROFILE_X 2.1 /* value of x to sample wave profile */ +#define HRES 1 /* dummy, only used by rde.c */ +#define INITIAL_SHIFT 20.0 /* time shift of initial wave packet (in oscillation periods) */ +#define WAVE_PACKET_SHIFT 200.0 /* time shift between wave packets (in oscillation periods) */ +#define FADE_IN_OBSTACLE 0 /* set to 1 to fade color inside obstacles */ +#define SHADE_2D 0 /* set to 1 to add pseudo-3d shading effect */ +#define SHADE_SCALE_2D 0.05 /* lower value increases sensitivity of shading */ +#define N_SOURCES 1 /* number of sources, for option draw_sources */ +double light[2] = {0.40824829, 0.816496581}; /* location of light source for SHADE_2D option*/ /* end of constants only used by sub_wave and sub_maze */ #include "global_pdes.c" diff --git a/lennardjones.c b/lennardjones.c index 7653aec..9a512f2 100644 --- a/lennardjones.c +++ b/lennardjones.c @@ -36,8 +36,8 @@ #include #include -#define MOVIE 0 /* set to 1 to generate movie */ -#define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */ +#define MOVIE 1 /* 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 */ @@ -50,8 +50,6 @@ /* General geometrical parameters */ -// #define WINWIDTH 1440 /* window width */ -// #define WINHEIGHT 810 /* window height */ #define WINWIDTH 1760 /* window width */ #define WINHEIGHT 990 /* window height */ @@ -60,20 +58,20 @@ #define YMIN -1.125 #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ -#define INITXMIN -1.9 -#define INITXMAX 2.1 /* x interval for initial condition */ -#define INITYMIN -1.0 -#define INITYMAX 1.0 /* y interval for initial condition */ +#define INITXMIN -1.95 +#define INITXMAX 1.95 /* x interval for initial condition */ +#define INITYMIN -1.1 +#define INITYMAX 1.1 /* 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 -1.95 -#define ADDXMAX 1.95 /* x interval for adding particles */ -#define ADDYMIN 1.4 -#define ADDYMAX 3.7 /* y interval for adding particles */ +#define ADDXMIN -1.5 +#define ADDXMAX 1.5 /* x interval for adding particles */ +#define ADDYMIN -1.0 +#define ADDYMAX 1.0 /* y interval for adding particles */ #define ADDRMIN 4.75 #define ADDRMAX 6.0 /* r interval for adding particles */ @@ -85,13 +83,13 @@ #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 CIRCLE_PATTERN 8 /* 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 1 /* pattern of circles for additional 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 6 /* pattern of obstacles, see list in global_ljones.c */ +#define OBSTACLE_PATTERN 71 /* pattern of obstacles, see list in global_ljones.c */ #define ADD_FIXED_SEGMENTS 0 /* set to 1 to add fixed segments as obstacles */ #define SEGMENT_PATTERN 29 /* pattern of repelling segments, see list in global_ljones.c */ @@ -108,23 +106,21 @@ #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 INTERACTION 12 /* particle interaction, see list in global_ljones.c */ #define INTERACTION_B 12 /* 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 SPIN_INTER_FREQUENCY 2.0 /* angular frequency of spin-spin interaction */ +#define SPIN_INTER_FREQUENCY_B 2.0 /* angular frequency of spin-spin interaction for second particle type */ +#define MOL_ANGLE_FACTOR 1.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 9.0 /* minimal distance in Poisson disc process, controls density of particles */ +#define PDISC_DISTANCE 9.6 /* 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.75 /* parameter controlling the dimensions of domain */ -#define MU 0.009 /* parameter controlling radius of particles */ -#define MU_B 0.009 /* parameter controlling radius of particles of second type */ +#define LAMBDA 0.2 /* parameter controlling the dimensions of domain */ +#define MU 0.0087 /* parameter controlling radius of particles */ +#define MU_B 0.012 /* 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 */ @@ -133,8 +129,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 20 /* number of grid point for grid of disks */ -#define NGRIDY 10 /* number of grid point for grid of disks */ +#define NGRIDX 40 /* number of grid point for grid of disks */ +#define NGRIDY 20 /* 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 */ @@ -150,11 +146,11 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 1600 /* number of frames of movie */ -// #define NSTEPS 7275 /* number of frames of movie */ -#define NVID 65 /* number of iterations between images displayed on screen */ +#define NSTEPS 2500 /* number of frames of movie */ +// #define NSTEPS 200 /* number of frames of movie */ +#define NVID 50 /* 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 INITIAL_TIME 5 /* 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 */ @@ -164,7 +160,7 @@ #define PSLEEP 1 /* sleep time during pause */ #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 MID_FRAMES 100 /* number of still frames between parts of two-part movie */ // #define END_FRAMES 250 /* number of still frames at end of movie */ #define END_FRAMES 100 /* number of still frames at end of movie */ @@ -174,15 +170,14 @@ /* Plot type, see list in global_ljones.c */ -#define PLOT 19 -// #define PLOT_B 1 /* plot type for second movie */ -#define PLOT_B 18 /* plot type for second movie */ +#define PLOT 5 +#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 0 /* type of background coloring, see list in global_ljones.c */ -#define BG_COLOR_B 2 /* type of background coloring, see list in global_ljones.c */ +#define COLOR_BACKGROUND 1 /* 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 1 /* set to 1 to draw bonds between neighbours */ #define COLOR_BONDS 1 /* set to 1 to color bonds according to length */ @@ -198,13 +193,13 @@ #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_ANGLE 10 /* Color palette for angle representation */ #define COLOR_PALETTE_DIRECTION 0 /* 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 0 /* Color palette for cluster representation */ +#define COLOR_PALETTE_CLUSTER 11 /* Color palette for cluster representation */ #define BLACK 1 /* background */ @@ -240,25 +235,29 @@ #define PARTICLE_HUE_MIN 359.0 /* color of original particle */ #define PARTICLE_HUE_MAX 0.0 /* color of saturated particle */ // #define PARTICLE_EMAX 50000.0 /* energy of particle with hottest color */ -#define PARTICLE_EMIN 10.0 /* energy of particle with coolest color */ -#define PARTICLE_EMAX 50000.0 /* energy of particle with hottest color */ -#define HUE_TYPE0 280.0 /* hue of particles of type 0 */ -#define HUE_TYPE1 280.0 /* hue of particles of type 1 */ -#define HUE_TYPE2 70.0 /* hue of particles of type 2 */ -#define HUE_TYPE3 60.0 /* hue of particles of type 3 */ +#define PARTICLE_EMIN 10.0 /* energy of particle with coolest color */ +#define PARTICLE_EMAX 20000.0 /* energy of particle with hottest color */ +#define HUE_TYPE0 300.0 /* hue of particles of type 0 */ +#define HUE_TYPE1 00.0 /* hue of particles of type 1 */ +#define HUE_TYPE2 340.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 0 /* set to 1 for random circle radius */ #define DT_PARTICLE 3.0e-6 /* time step for particle displacement */ -#define KREPEL 150.0 /* constant in repelling force between particles */ +#define KREPEL 15.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 REPEL_RADIUS 25.0 /* radius in which repelling force acts (in units of particle radius) */ -#define DAMPING 20.0 /* damping coefficient of particles */ +#define DAMPING 500.0 /* damping coefficient of particles */ #define INITIAL_DAMPING 5000.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 2.0 /* mass of particle of radius MU_B */ +#define DAMPING_ROT 1.0e6 /* damping coefficient for rotation of particles */ +#define PARTICLE_MASS 2.0 /* mass of particle of radius MU */ +#define PARTICLE_MASS_B 16.0 /* mass of particle of radius MU_B */ #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 50.0 /* initial velocity range */ @@ -271,11 +270,11 @@ #define THERMOSTAT 1 /* 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.0005 /* initial inverse temperature */ +#define BETA 0.00007 /* 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 2.0e11 /* harmonic potential of obstacles */ -#define NBH_DIST_FACTOR 5.0 /* radius in which to count neighbours */ +#define NBH_DIST_FACTOR 5.0 /* radius in which to count neighbours */ #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 */ @@ -288,39 +287,40 @@ #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 50000.0 /* value of electric field */ +#define EFIELD 100000.0 /* value of electric field */ #define ADD_BFIELD 0 /* set to 1 to add a magnetic field */ #define BFIELD 2.666666667 /* value of magnetic field */ -#define CHARGE -0.0 /* charge of particles of first type */ -#define CHARGE_B 0.0 /* charge of particles of second type */ +#define CHARGE 1.0 /* charge of particles of first type */ +#define CHARGE_B -1.5 /* charge of particles of second type */ #define INCREASE_E 0 /* set to 1 to increase electric field */ // #define EFIELD_FACTOR 2500000.0 /* factor by which 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 20000.0 /* factor by which to increase magnetic field */ -#define CHARGE_OBSTACLES 1 /* set to 1 for obstacles to be charged */ +#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 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 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 50.0 /* force constant in angular dynamics */ +#define KTORQUE -100.0 /* force constant in angular dynamics */ #define KTORQUE_BOUNDARY 1.0e6 /* constant in torque from the boundary */ #define KTORQUE_B 10.0 /* force constant in angular dynamics */ -#define KTORQUE_DIFF 150.0 /* force constant in angular dynamics for different particles */ +#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_MINUS 1 /* 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 5.0 /* range of spin-spin interaction for second type of particle */ +#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 0 /* type of temperature schedule, see TS_* in global_ljones */ -// #define BETA_FACTOR 0.000001 /* factor by which to change BETA during simulation */ -#define BETA_FACTOR 1000.0 /* factor by which to change BETA during simulation */ +#define INCREASE_BETA 1 /* set to 1 to increase BETA during simulation */ +#define BETA_SCHEDULE 3 /* type of temperature schedule, see TS_* in global_ljones */ +#define BETA_FACTOR 10.0 /* factor by which to change BETA during simulation */ +// #define BETA_FACTOR 0.08 /* 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 */ @@ -356,7 +356,7 @@ #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 1000.0 /* factor by which to change 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 */ @@ -364,10 +364,10 @@ #define ADD_PARTICLES 0 /* 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 5 /* time interval between adding further particles */ -#define N_ADD_PARTICLES 5 /* number of particles to add */ -#define FINAL_NOADD_PERIOD 0 /* final period where no particles are added */ +#define ADD_TIME 50 /* time at which to add first particle */ +#define ADD_PERIOD 2 /* time interval between adding further particles */ +#define N_ADD_PARTICLES 8 /* number of particles to add */ +#define FINAL_NOADD_PERIOD 4700 /* 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 0 /* set to 1 to have a tracer particle */ @@ -417,20 +417,21 @@ #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 1 /* set to 1 to simulate a chemical reaction (particles may change type) */ -#define RD_REACTION 23 /* type of reaction, see list in global_ljones.c */ -#define RD_TYPES 2 /* number of types in reaction-diffusion equation */ -#define RD_INITIAL_COND 99 /* initial condition of particles */ -#define REACTION_DIST 3.5 /* maximal distance for reaction to occur */ +#define REACTION_DIFFUSION 0 /* set to 1 to simulate a chemical reaction (particles may change type) */ +#define RD_REACTION 256 /* type of reaction, see list in global_ljones.c */ +#define RD_TYPES 6 /* number of types in reaction-diffusion equation */ +#define RD_INITIAL_COND 10 /* initial condition of particles */ +#define REACTION_DIST 4.0 /* maximal distance for reaction to occur */ #define REACTION_PROB 1.0 /* probability controlling reaction term */ -#define DISSOCIATION_PROB 0.0001 /* probability controlling dissociation reaction */ +#define DISSOCIATION_PROB 0.0 /* probability controlling dissociation reaction */ +#define KILLING_PROB 0.0015 /* probability of enzymes being killed */ #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 DELTAVMAX 150.0 /* maximal deltav allowed for pairing molecules */ -#define AGREGMAX 6 /* maximal number of partners for CHEM_AGGREGATION reaction */ -#define AGREG_DECOUPLE 10 /* minimal number of partners to decouple from thermostat */ +#define COLLISION_TIME 25 /* time during which collisions are shown */ +#define DELTAVMAX 1000.0 /* maximal deltav allowed for pairing molecules */ +#define AGREGMAX 11 /* maximal number of partners for CHEM_AGGREGATION reaction */ +#define AGREG_DECOUPLE 12 /* minimal number of partners to decouple from thermostat */ #define CHANGE_RADIUS 0 /* set to 1 to change particle radius during simulation */ #define MU_RATIO 0.666666667 /* ratio by which to increase radius */ @@ -458,29 +459,34 @@ #define PROP_MAX 0.9 /* max proportion of type 1 particles */ #define PAIR_PARTICLES 1 /* set to 1 to form particles pairs */ -#define RANDOMIZE_ANGLE 1 /* set to 1 for random orientation */ -#define DEACIVATE_CLOSE_PAIRS 0 /* set to 1 to test for closeness to other particles */ +#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 KSPRING_PAIRS 2.0e10 /* spring constant for pair interaction */ -#define NPARTNERS 16 /* number of partners of particles */ -#define NARMS 4 /* number of "arms" for certain paring types */ -#define PAIRING_TYPE 61 /* type of pairing, see POLY_ in global_ljones.c */ +#define THIRD_TYPE_PROPORTION 1.0 /* proportion of third type pairings, for certain pairing types */ + +#define KSPRING_PAIRS 1.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 8 /* number of partners of particles */ +#define NARMS 1 /* number of "arms" for certain paring types */ +#define PAIRING_TYPE 42 /* 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.0 /* ratio between equilibrium distance and radius (default: 1.0) */ -#define MU_C 0.014 /* radius of partner particle */ -#define PARTICLE_MASS_C 2.0 /* mass or partner particle */ -#define CHARGE_C 1.5 /* charge of partner particle */ +#define PAIR_DRATIO 0.9 /* ratio between equilibrium distance and radius (default: 1.0) */ +#define MU_C 0.0125 /* radius of partner particle */ +#define PARTICLE_MASS_C 8.0 /* mass or partner particle */ +#define CHARGE_C -1.0 /* charge of partner particle */ #define CLUSTER_COLOR_FACTOR 400 /* 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 16 /* number of partners of particles */ -#define NARMS_B 4 /* number of "arms" for certain paring types */ -#define PAIRING_TYPE_B 61 /* type of pairing, see POLY_ in global_ljones.c */ -#define MU_D 0.014 /* radius of partner particle */ -#define PARTICLE_MASS_D 2.0 /* mass or partner particle */ -#define CHARGE_D -1.5 /* charge of partner particle */ -// #define PARTNER_ANGLE_B 104.45 /* angle (in degrees) between anions for POLY_WATER case */ +#define NPARTNERS_B 2 /* number of partners of particles */ +#define NARMS_B 1 /* number of "arms" for certain paring types */ +#define PAIRING_TYPE_B 2 /* type of pairing, see POLY_ in global_ljones.c */ +#define MU_D 0.008 /* radius of partner particle */ +#define PARTICLE_MASS_D 1.0 /* mass or partner particle */ +#define CHARGE_D 0.75 /* charge of partner particle */ #define NXMAZE 12 /* width of maze */ #define NYMAZE 12 /* height of maze */ @@ -494,9 +500,9 @@ #define FLOOR_OMEGA 0 /* set to 1 to limit particle momentum to PMAX */ #define PMAX 1000.0 /* maximal force */ -#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 HASHX 80 /* size of hashgrid in x direction */ +#define HASHY 40 /* 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 */ @@ -512,7 +518,8 @@ #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)) #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))) +#define PAIR_FORCE ((PAIR_PARTICLES)||((REACTION_DIFFUSION)&&((RD_REACTION == CHEM_AGGREGATION)||(RD_REACTION == CHEM_AGGREGATION_CHARGE)||(RD_REACTION == CHEM_AGGREGATION_NNEIGH)))) +#define COMPUTE_PAIR_TORQUE (KTORQUE_PAIR_ANGLE != 0.0) double xshift = 0.0; /* x shift of shown window */ double xspeed = 0.0; /* x speed of obstacle */ @@ -654,6 +661,16 @@ double temperature_schedule(int i) 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; @@ -705,6 +722,16 @@ double temperature_schedule(int i) 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; @@ -1010,6 +1037,7 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA px[j] *= exp(- DT_PARTICLE*damping); py[j] *= exp(- DT_PARTICLE*damping); pangle[j] *= exp(- DT_PARTICLE*DAMPING_ROT); +// printf("Damping particle angular velocity\n"); } if ((THERMOSTAT_ON)&&(COUPLE_ANGLE_TO_THERMOSTAT)&&(particle[j].thermostat)) pangle[j] *= exp(- 0.5*DT_PARTICLE*xi); @@ -1365,7 +1393,8 @@ void evolve_segment_groups(t_segment segment[NMAXSEGMENTS], int time, t_group_se void animation() { double time, scale, diss, rgb[3], dissip, gradient[2], x, y, dx, dy, dt, xleft, xright, - a, b, length, fx, fy, force[2], totalenergy = 0.0, 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; + 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; double *qx, *qy, *px, *py, *qangle, *pangle, *pressure, *obstacle_speeds; int i, j, k, n, m, s, ij[2], i0, iplus, iminus, j0, jplus, jminus, p, q, p1, q1, p2, q2, total_neighbours = 0, min_nb, max_nb, close, wrapx = 0, wrapy = 0, nadd_particle = 0, nmove = 0, nsuccess = 0, @@ -1382,6 +1411,7 @@ void animation() t_group_data *group_speeds; t_collision *collisions; t_hashgrid *hashgrid; + t_molecule *molecule; t_lj_parameters params; char message[100]; @@ -1407,6 +1437,9 @@ void animation() 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)); @@ -1429,6 +1462,8 @@ void animation() 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); @@ -1445,7 +1480,11 @@ void animation() init_particle_config(particle); /* add some particles, beta */ - if (ADD_INITIAL_PARTICLES) add_particle_config(particle, -0.6, 1.6, -1.0, 1.0, MU_B); + if (ADD_INITIAL_PARTICLES) + { + add_particle_config(particle, INITXMIN - 0.5, INITXMAX + 0.5, INITYMAX + 0.1, YMAX-0.1, NGRIDX*5/4, 3, 0, MU); + add_particle_config(particle, INITXMIN - 0.5, INITXMAX + 0.5, YMIN + 0.1, INITYMIN-0.1, NGRIDX*5/4, 3, 0, MU); + } init_hashgrid(hashgrid); @@ -1461,7 +1500,7 @@ void animation() printf("Initializing configuration\n"); - params.nactive = initialize_configuration(particle, hashgrid, obstacle, px, py, pangle, tracer_n, segment); + params.nactive = initialize_configuration(particle, hashgrid, obstacle, px, py, pangle, tracer_n, segment, molecule); printf("%i active particles\n", params.nactive); @@ -1486,8 +1525,12 @@ void animation() for (i=0; i<=INITIAL_TIME + NSTEPS; i++) { - printf("Computing frame %d\n",i); + /* TEST */ +// if (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); @@ -1615,10 +1658,44 @@ void animation() 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 INITIAL_TIME)&&(REACTION_DIFFUSION)) { - ncollisions = update_types(particle, collisions, ncollisions, particle_numbers, i - INITIAL_TIME - 1, &delta_energy); + ncollisions = update_types(particle, molecule, 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)) { - for (k=0; k #include -#define MOVIE 0 /* set to 1 to generate movie */ +#define MOVIE 1 /* set to 1 to generate movie */ #define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */ #define SAVE_MEMORY 1 /* set to 1 to save memory when writing tiff images */ #define NO_EXTRA_BUFFER_SWAP 1 /* some OS require one less buffer swap when recording images */ @@ -48,8 +48,13 @@ #define WINWIDTH 1920 /* window width */ #define WINHEIGHT 1150 /* window height */ -#define NX 1500 /* number of grid points on x axis */ -#define NY 750 /* number of grid points on y axis */ +// #define NX 240 /* number of grid points on x axis */ +// #define NY 120 /* number of grid points on y axis */ +// #define NX 960 /* number of grid points on x axis */ +// #define NY 500 /* number of grid points on y axis */ +#define NX 780 /* number of grid points on x axis */ +#define NY 400 /* number of grid points on y axis */ +#define HRES 2 /* factor for high resolution plots */ #define XMIN -2.0 #define XMAX 2.0 /* x interval */ @@ -58,28 +63,32 @@ /* Choice of simulated equation */ -#define RDE_EQUATION 41 /* choice of reaction term, see list in global_3d.c */ -#define NFIELDS 5 /* number of fields in reaction-diffusion equation */ -#define NLAPLACIANS 5 /* number of fields for which to compute Laplacian */ +#define RDE_EQUATION 7 /* choice of reaction term, see list in global_3d.c */ +#define NFIELDS 3 /* number of fields in reaction-diffusion equation */ +#define NLAPLACIANS 0 /* number of fields for which to compute Laplacian */ #define SPHERE 1 /* set to 1 to simulate equation on sphere */ #define DPOLE 1 /* safety distance to poles */ -#define DSMOOTH 10 /* size of neighbourhood of poles that are smoothed */ -#define SMOOTHPOLE 0.24 /* smoothing coefficient at poles */ +#define DSMOOTH 1 /* size of neighbourhood of poles that are smoothed */ +#define SMOOTHPOLE 0.01 /* smoothing coefficient at poles */ +#define SMOOTHCOTPOLE 0.01 /* smoothing coefficient of cotangent at poles */ #define PHISHIFT 0.0 /* shift of phi in 2D plot (in degrees) */ +#define SMOOTHBLOCKS 1 /* set to 1 to use blocks of points near the poles */ +#define BLOCKDIST 64 /* distance to poles where points are blocked */ +#define ZERO_MERIDIAN 190.0 /* choice of zero meridian (will be at left/right boundary of 2d plot) */ #define ADD_POTENTIAL 0 /* set to 1 to add a potential (for Schrodinger equation) */ #define ADD_MAGNETIC_FIELD 0 /* set to 1 to add a magnetic field (for Schrodinger equation) - then set POTENTIAL 1 */ -#define ADD_FORCE_FIELD 0 /* set to 1 to add a foce field (for compressible Euler equation) */ +#define ADD_FORCE_FIELD 1 /* set to 1 to add a foce field (for compressible Euler equation) */ #define POTENTIAL 7 /* type of potential or vector potential, see list in global_3d.c */ -#define FORCE_FIELD 5 /* type of force field, see list in global_3d.c */ -#define ADD_CORIOLIS_FORCE 0 /* set to 1 to add Coriolis force (quasigeostrophic Euler equations) */ +#define FORCE_FIELD 6 /* type of force field, see list in global_3d.c */ +#define ADD_CORIOLIS_FORCE 1 /* set to 1 to add Coriolis force (quasigeostrophic Euler equations) */ #define VARIABLE_DEPTH 0 /* set to 1 for variable depth in shallow water equation */ #define SWATER_DEPTH 4 /* variable depth in shallow water equation */ #define ANTISYMMETRIZE_WAVE_FCT 0 /* set tot 1 to make wave function antisymmetric */ -#define ADAPT_STATE_TO_BC 0 /* to smoothly adapt initial state to obstacles */ -#define OBSTACLE_GEOMETRY 1 /* geometry of obstacles, as in B_DOMAIN */ +#define ADAPT_STATE_TO_BC 1 /* to smoothly adapt initial state to obstacles */ +#define OBSTACLE_GEOMETRY 84 /* geometry of obstacles, as in B_DOMAIN */ #define BC_STIFFNESS 50.0 /* controls region of boundary condition control */ #define CHECK_INTEGRAL 1 /* set to 1 to check integral of first field */ @@ -98,10 +107,10 @@ #define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */ #define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */ -#define LAMBDA 0.9 /* parameter controlling the dimensions of domain */ -#define MU 0.06 /* parameter controlling the dimensions of domain */ +#define LAMBDA 1.0 /* parameter controlling the dimensions of domain */ +#define MU 1.0 /* parameter controlling the dimensions of domain */ #define NPOLY 5 /* number of sides of polygon */ -#define APOLY 2.0 /* angle by which to turn polygon, in units of Pi/2 */ +#define APOLY 2.0 /* angle by which to turn polygon, in units of Pi/2 */ #define MDEPTH 7 /* depth of computation of Menger gasket */ #define MRATIO 5 /* ratio defining Menger gasket */ #define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ @@ -128,14 +137,20 @@ /* Physical parameters of wave equation */ -#define DT 0.0000002 +// #define DT 0.0000001 +#define DT 0.00000015 +// #define DT 0.0000002 +// #define DT 0.0000012 +// #define DT 0.000003 +// #define DT 0.0000022 +// #define DT 0.0000012 -#define VISCOSITY 0.075 +#define VISCOSITY 0.02 #define POISSON_STIFFNESS 1.0 /* stiffness of Poisson equation solver for incompressible Euler */ #define DISSIPATION 0.0 #define RPSA 0.75 /* parameter in Rock-Paper-Scissors-type interaction */ -#define RPSLZB 0.75 /* second parameter in Rock-Paper-Scissors-Lizard-Spock type interaction */ +#define RPSLZB 0.0 /* second parameter in Rock-Paper-Scissors-Lizard-Spock type interaction */ #define K_AC 0.1 /* force constant in Allen-Cahn equation */ #define EPSILON 0.8 /* time scale separation */ @@ -148,23 +163,24 @@ #define BZQ 0.0008 /* parameter in BZ equation */ #define BZF 1.2 /* parameter in BZ equation */ #define B_FIELD 10.0 /* magnetic field */ -#define G_FIELD 0.75 /* gravity */ +#define G_FIELD 0.03 /* gravity/Coriolis force */ #define BC_FIELD 1.0e-5 /* constant in repulsive field from obstacles */ #define AB_RADIUS 0.2 /* radius of region with magnetic field for Aharonov-Bohm effect */ #define K_EULER 50.0 /* constant in stream function integration of Euler equation */ #define K_EULER_INC 0.5 /* constant in incompressible Euler equation */ +#define C_EULER_COMP 0.1 /* constant in compressible Euler equation */ #define SMOOTHEN_VORTICITY 0 /* set to 1 to smoothen vorticity field in Euler equation */ -#define SMOOTHEN_VELOCITY 0 /* set to 1 to smoothen velocity field in Euler equation */ +#define SMOOTHEN_VELOCITY 1 /* set to 1 to smoothen velocity field in Euler equation */ #define SMOOTHEN_PERIOD 10 /* period between smoothenings */ -#define SMOOTH_FACTOR 0.15 /* factor by which to smoothen */ +#define SMOOTH_FACTOR 0.05 /* factor by which to smoothen */ #define ADD_OSCILLATING_SOURCE 0 /* set to 1 to add an oscillating wave source */ #define OSCILLATING_SOURCE_PERIOD 1 /* period of oscillating source */ #define OSCILLATING_SOURCE_OMEGA 0.2 /* frequency of oscillating source */ -#define ADD_TRACERS 0 /* set to 1 to add tracer particles (for Euler equations) */ -#define N_TRACERS 1000 /* number of tracer particles */ +#define ADD_TRACERS 1 /* set to 1 to add tracer particles (for Euler equations) */ +#define N_TRACERS 2000 /* number of tracer particles */ #define TRACERS_STEP 0.005 /* step size in tracer evolution */ #define T_OUT 2.0 /* outside temperature */ @@ -181,10 +197,9 @@ #define ADJUST_INTSTEP 0 /* set to 1 to decrease integration step when viscosity increases */ #define VISCOSITY_INITIAL_TIME 10 /* initial time during which viscosity remains constant */ #define VISCOSITY_FACTOR 100.0 /* factor by which to change viscosity */ -#define VISCOSITY_MAX 2.0 /* max value of viscosity beyond which NVID is increased and integration step is decrase, - for numerical stability */ +#define VISCOSITY_MAX 2.0 /* max value of viscosity beyond which NVID is increased and integration step is decrase, for numerical stability */ -#define CHANGE_RPSLZB 1 /* set to 1 to change second parameter in Rock-Paper-Scissors-Lizard-Spock equation */ +#define CHANGE_RPSLZB 0 /* set to 1 to change second parameter in Rock-Paper-Scissors-Lizard-Spock equation */ #define RPSLZB_CHANGE 0.75 /* factor by which to rpslzb parameter */ #define RPSLZB_INITIAL_TIME 0 /* initial time during which rpslzb remains constant */ #define RPSLZB_FINAL_TIME 500 /* final time during which rpslzb remains constant */ @@ -215,11 +230,11 @@ #define B_COND_TOP 0 #define B_COND_BOTTOM 0 - /* Parameters for length and speed of simulation */ -#define NSTEPS 3500 /* number of frames of movie */ -#define NVID 8 /* number of iterations between images displayed on screen */ +#define NSTEPS 1900 /* number of frames of movie */ +// #define NSTEPS 500 /* number of frames of movie */ +#define NVID 100 /* number of iterations between images displayed on screen */ #define ACCELERATION_FACTOR 1.0 /* factor by which to increase NVID in course of simulation */ #define DT_ACCELERATION_FACTOR 1.0 /* factor by which to increase time step in course of simulation */ #define MAX_DT 0.024 /* maximal value of integration step */ @@ -242,23 +257,25 @@ #define ROTATE_VIEW 1 /* set to 1 to rotate position of observer */ #define ROTATE_ANGLE 360.0 /* total angle of rotation during simulation */ +// #define ROTATE_ANGLE 90.0 /* total angle of rotation during simulation */ #define SHADE_3D 1 /* set to 1 to change luminosity according to normal vector */ #define SHADE_2D 0 /* set to 1 to change luminosity according to normal vector */ +#define VIEWPOINT_TRAJ 0 /* type of viewpoint trajectory */ +#define MAX_LATITUDE 45.0 /* maximal latitude for viewpoint trajectory VP_ORBIT2 */ #define DRAW_PERIODICISED 0 /* set to 1 to repeat wave periodically in x and y directions */ /* Plot type - color scheme */ -#define CPLOT 40 -#define CPLOT_B 42 +#define CPLOT 62 +#define CPLOT_B 64 /* Plot type - height of 3D plot */ -#define ZPLOT 42 /* z coordinate in 3D plot */ -#define ZPLOT_B 42 /* z coordinate in second 3D plot */ +#define ZPLOT 62 /* z coordinate in 3D plot */ +#define ZPLOT_B 64 /* z coordinate in second 3D plot */ #define AMPLITUDE_HIGH_RES 1 /* set to 1 to increase resolution of P_3D_AMPLITUDE plot */ -#define SHADE_3D 1 /* set to 1 to change luminosity according to normal vector */ #define NON_DIRICHLET_BC 0 /* set to 1 to draw only facets in domain, if field is not zero on boundary */ #define WRAP_ANGLE 1 /* experimental: wrap angle to [0, 2Pi) for interpolation in angle schemes */ #define FADE_IN_OBSTACLE 0 /* set to 1 to fade color inside obstacles */ @@ -275,7 +292,7 @@ #define PRINT_TIME 0 /* set to 1 to print running time */ #define PRINT_VISCOSITY 0 /* set to 1 to print viscosity */ -#define PRINT_RPSLZB 1 /* set to 1 to print rpslzb parameter */ +#define PRINT_RPSLZB 0 /* set to 1 to print rpslzb parameter */ #define PRINT_PROBABILITIES 0 /* set to 1 to print probabilities (for Ehrenfest urn configuration) */ #define PRINT_NOISE 0 /* set to 1 to print noise intensity */ #define PRINT_FLOW_SPEED 0 /* set to 1 to print speed of flow */ @@ -299,8 +316,8 @@ /* Color schemes, see list in global_pdes.c */ -#define COLOR_PALETTE 0 /* Color palette, see list in global_pdes.c */ -#define COLOR_PALETTE_B 17 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE 10 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE_B 17 /* Color palette, see list in global_pdes.c */ #define BLACK 1 /* black background */ #define COLOR_OUT_R 1.0 /* color outside domain */ @@ -313,17 +330,17 @@ #define SCALE 0 /* set to 1 to adjust color scheme to variance of field */ #define SLOPE 1.0 /* sensitivity of color on wave amplitude */ -#define VSCALE_AMPLITUDE 0.5 /* additional scaling factor for color scheme P_3D_AMPLITUDE */ +#define VSCALE_AMPLITUDE 100.0 /* additional scaling factor for color scheme P_3D_AMPLITUDE */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ -#define CURL_SCALE 0.000015 /* scaling factor for curl representation */ +#define CURL_SCALE 1.0 /* scaling factor for curl representation */ #define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */ -#define SLOPE_SCHROD_LUM 400.0 /* sensitivity of luminosity on module, for color scheme Z_ARGUMENT */ -#define MIN_SCHROD_LUM 0.2 /* minimal luminosity in color scheme Z_ARGUMENT*/ +#define SLOPE_SCHROD_LUM 100.0 /* sensitivity of luminosity on module, for color scheme Z_ARGUMENT */ +#define MIN_SCHROD_LUM 0.1 /* minimal luminosity in color scheme Z_ARGUMENT*/ #define VSCALE_PRESSURE 2.0 /* additional scaling factor for color scheme Z_EULER_PRESSURE */ #define PRESSURE_SHIFT 10.0 /* shift for color scheme Z_EULER_PRESSURE */ #define PRESSURE_LOG_SHIFT -2.5 /* shift for color scheme Z_EULER_PRESSURE */ #define VSCALE_WATER_HEIGHT 0.4 /* vertical scaling of water height */ -#define SHADE_SCALE_2D 1.0 /* controls "depth" of 2D shading */ +#define SHADE_SCALE_2D 0.25 /* controls "depth" of 2D shading */ #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 */ @@ -336,13 +353,15 @@ #define LOG_SCALE 0.5 /* scaling factor for energy log representation */ #define LOG_SHIFT 1.0 #define LOG_MIN 1.0e-3 /* floor value for log vorticity plot */ -#define VSCALE_SPEED 200.0 /* additional scaling factor for color scheme Z_EULER_SPEED */ +#define VSCALE_SPEED 50.0 /* additional scaling factor for color scheme Z_EULER_SPEED */ #define VMEAN_SPEED 0.0 /* mean value around which to scale for color scheme Z_EULER_SPEED */ -#define SHIFT_DENSITY 8.5 /* shift for color scheme Z_EULER_DENSITY */ -#define VSCALE_DENSITY 3.0 /* additional scaling factor for color scheme Z_EULER_DENSITY */ -#define VSCALE_VORTICITY 20.0 /* additional scaling factor for color scheme Z_EULERC_VORTICITY */ +#define SHIFT_DENSITY 1.0 /* shift for color scheme Z_EULER_DENSITY */ +#define VSCALE_DENSITY 30.0 /* additional scaling factor for color scheme Z_EULER_DENSITY */ +#define VSCALE_VORTICITY 15.0 /* additional scaling factor for color scheme Z_EULERC_VORTICITY */ #define VORTICITY_SHIFT 0.0 /* vertical shift of vorticity */ -#define ZSCALE_SPEED 0.3 /* additional scaling factor for z-coord Z_EULER_SPEED and Z_SWATER_SPEED */ +#define ZSCALE_SPEED 0.5 /* additional scaling factor for z-coord Z_EULER_SPEED and Z_SWATER_SPEED */ +#define ZSHIFT_SPEED 0.0 /* additional shift of z-coord Z_EULER_SPEED and Z_SWATER_SPEED */ +#define ZSCALE_NORMGRADIENT -0.0001 /* vertical scaling for Z_NORM_GRADIENT */ #define VSCALE_SWATER 250.0 /* additional scaling factor for color scheme Z_EULER_DENSITY */ #define NXMAZE 7 /* width of maze */ @@ -352,11 +371,11 @@ #define MAZE_XSHIFT 0.0 /* horizontal shift of maze */ #define MAZE_WIDTH 0.04 /* half width of maze walls */ -#define DRAW_COLOR_SCHEME 0 /* set to 1 to plot the color scheme */ -#define COLORBAR_RANGE 2.5 /* scale of color scheme bar */ +#define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */ +#define COLORBAR_RANGE 3.0 /* scale of color scheme bar */ #define COLORBAR_RANGE_B 2.5 /* scale of color scheme bar for 2nd part */ #define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ -#define CIRC_COLORBAR 1 /* set to 1 to draw circular color scheme */ +#define CIRC_COLORBAR 0 /* set to 1 to draw circular color scheme */ #define CIRC_COLORBAR_B 1 /* set to 1 to draw circular color scheme */ /* only for compatibility with wave_common.c */ @@ -386,9 +405,12 @@ #define FLUX_WINDOW 20 /* averaging window for energy flux */ #define ADD_WAVE_PACKET_SOURCES 0 /* set to 1 to add several sources emitting wave packets */ #define WAVE_PACKET_SOURCE_TYPE 1 /* type of wave packet sources */ +#define N_SOURCES 1 /* number of wave sources */ #define N_WAVE_PACKETS 15 /* number of wave packets */ #define WAVE_PACKET_RADIUS 20 /* radius of wave packets */ #define OSCIL_LEFT_YSHIFT 25.0 /* y-dependence of left oscillation (for non-horizontal waves) */ +#define INITIAL_SHIFT 20.0 /* time shift of initial wave packet (in oscillation periods) */ +#define WAVE_PACKET_SHIFT 200.0 /* time shift between wave packets (in oscillation periods) */ #define DRAW_WAVE_PROFILE 0 /* set to 1 to draw a profile of the wave */ #define HORIZONTAL_WAVE_PROFILE 0 /* set to 1 to draw wave profile vertically */ #define VERTICAL_WAVE_PROFILE 0 /* set to 1 to draw wave profile vertically */ @@ -403,33 +425,48 @@ #define DRAW_WAVE_TIMESERIES 0 /* set to 1 to draw a time series of the wave */ #define TIMESERIES_NVALUES 400 /* number of values plotted in time series */ #define DRAW_WAVE_SOURCE 0 /* set to 1 to draw source of wave at (wave_source_x, wave_source_y) */ -#define MESSAGE_LDASH 14 /* length of dash for Morse code message */ -#define MESSAGE_LDOT 8 /* length of dot for Morse code message */ -#define MESSAGE_LINTERVAL 54 /* length of interval between dashes/dots for Morse code message */ -#define MESSAGE_LINTERLETTER 60 /* length of interval between letters for Morse code message */ -#define MESSAGE_LSPACE 48 /* length of space for Morse code message */ -#define MESSAGE_INITIAL_TIME 100 /* initial time before starting message for Morse code message */ +#define MESSAGE_LDASH 1 /* length of dash for Morse code message */ +#define MESSAGE_LDOT 1 /* length of dot for Morse code message */ +#define MESSAGE_LINTERVAL 1 /* length of interval between dashes/dots for Morse code message */ +#define MESSAGE_LINTERLETTER 1 /* length of interval between letters for Morse code message */ +#define MESSAGE_LSPACE 1 /* length of space for Morse code message */ +#define MESSAGE_INITIAL_TIME 1 /* initial time before starting message for Morse code message */ /* end of constants added only for compatibility with wave_common.c */ - double u_3d[2] = {0.75, -0.45}; /* projections of basis vectors for REP_AXO_3D representation */ double v_3d[2] = {-0.75, -0.45}; double w_3d[2] = {0.0, 0.015}; -double light[3] = {0.816496581, 0.40824829, 0.40824829}; /* vector of "light" direction for P_3D_ANGLE color scheme */ -double observer[3] = {8.0, 8.0, 7.0}; /* location of observer for REP_PROJ_3D representation */ +double light[3] = {-0.40824829, -0.816496581, 0.40824829}; /* vector of "light" direction for P_3D_ANGLE color scheme */ +double observer[3] = {-8.0, -4.0, 4.0}; /* location of observer for REP_PROJ_3D representation */ int reset_view = 0; /* switch to reset 3D view parameters (for option ROTATE_VIEW) */ +/* constants for simulations on planets */ +#define ADD_DEM 1 /* add DEM (digital elevation model) */ +#define ADD_NEGATIVE_DEM 0 /* add DEM with bathymetric data */ +#define RSCALE_DEM 0.1 /* scaling factor of radial component for DEM */ +#define SMOOTH_DEM 0 /* set to 1 to smoothen DEM (to make altitude less constant) */ +#define DEM_SMOOTH_STEPS 1 /* number of smoothening steps */ +#define DEM_SMOOTH_HEIGHT 2.0 /* relative height below which to smoothen */ +#define DEM_MAXHEIGHT 9000.0 /* max height of DEM (estimated from Everest/Olympus Mons) */ +#define DEM_MAXDEPTH -10000 /* max depth of DEM */ +#define PLANET_SEALEVEL 0.0 /* sea level for flooded planet */ +#define VENUS_NODATA_FACTOR 0.5 /* altitude to assign to DEM points without data (fraction of mean altitude) */ + #define Z_SCALING_FACTOR 0.8 /* overall scaling factor of z axis for REP_PROJ_3D representation */ #define XY_SCALING_FACTOR 2.0 /* overall scaling factor for on-screen (x,y) coordinates after projection */ #define ZMAX_FACTOR 1.0 /* max value of z coordinate for REP_PROJ_3D representation */ #define XSHIFT_3D 0.0 /* overall x shift for REP_PROJ_3D representation */ -#define YSHIFT_3D 0.0 /* overall y shift for REP_PROJ_3D representation */ +#define YSHIFT_3D 0.0 /* overall y shift for REP_PROJ_3D representation */ #define BORDER_PADDING 0 /* distance from boundary at which to plot points, to avoid boundary effects due to gradient */ +#define DRAW_ARROW 0 /* set to 1 to draw arrow above sphere */ + +#define RSCALE 0.01 /* scaling factor of radial component */ +#define RSHIFT -0.01 /* shift in radial component */ +#define RMAX 2.0 /* max value of radial component */ +#define RMIN 0.5 /* min value of radial component */ +// #define COS_VISIBLE -1.1 /* limit on cosine of normal to shown facets */ +#define COS_VISIBLE -0.3 /* limit on cosine of normal to shown facets */ -#define RSCALE -0.01 /* scaling factor of radial component */ -#define RMAX 1.005 /* max value of radial component */ -#define RMIN 0.995 /* min value of radial component */ -#define COS_VISIBLE -1.1 /* limit on cosine of normal to shown facets */ /* For debugging purposes only */ #define FLOOR 1 /* set to 1 to limit wave amplitude to VMAX */ @@ -443,6 +480,9 @@ int reset_view = 0; /* switch to reset 3D view parameters (for option RO #define ASYM_SPEED_COLOR (VMEAN_SPEED == 0.0) +int block_sizes[NY]; /* table of block sizes for blocking around poles */ +int block_numbers[NY]; /* table of block numbers for blocking around poles */ + #include "global_pdes.c" #include "global_3d.c" /* constants and global variables */ @@ -665,7 +705,7 @@ void initialize_vector_potential(double vpotential_field[2*NX*NY]) } } -void initialize_gfield(double gfield[2*NX*NY], double bc_field[NX*NY], double bc_field2[NX*NY]) +void initialize_gfield(double gfield[2*NX*NY], double bc_field[NX*NY], double bc_field2[NX*NY], t_wave_sphere *wsphere, t_wave_sphere *wsphere_hr) /* initialize the exterior field, e.g. for the compressible Euler equation */ { int i, j; @@ -701,7 +741,37 @@ void initialize_gfield(double gfield[2*NX*NY], double bc_field[NX*NY], double bc gfield[NX*NY+(NX-1)*NY+j] = 0.0; } } - + else if (FORCE_FIELD == GF_EARTH) + { + dx = (XMAX - XMIN)/(double)NX; + dy = (YMAX - YMIN)/(double)NY; + init_earth_map_rde(wsphere, 1); + init_earth_map_rde(wsphere_hr, HRES); + + #pragma omp parallel for private(i,j) + for (i=1; i 0.0) { @@ -1213,7 +1284,7 @@ double gfield[2*NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY]) vx = rde[i*NY+j].dxv; vy = rde[i*NY+j].dyv; - phi_out[0][i*NY+j] = rho - intstep*(u*rhox + v*rhoy + rho*(ux + vy)); + phi_out[0][i*NY+j] = rho - intstep*C_EULER_COMP*(u*rhox + v*rhoy + rho*(ux + vy)); phi_out[1][i*NY+j] = u - intstep*(u*ux + v*uy + K_EULER_INC*rhox/rho); phi_out[2][i*NY+j] = v - intstep*(u*vx + v*vy + K_EULER_INC*rhoy/rho); @@ -1224,8 +1295,22 @@ double gfield[2*NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY]) } if (ADD_CORIOLIS_FORCE) { - phi_out[1][i*NY+j] += intstep*G_FIELD*v; - phi_out[2][i*NY+j] -= intstep*G_FIELD*u; + if (SPHERE) + { + phi_out[1][i*NY+j] += intstep*G_FIELD*v*wsphere[i*NY+j].ctheta; + phi_out[2][i*NY+j] -= intstep*G_FIELD*u*wsphere[i*NY+j].reg_cottheta; +// phi_out[1][i*NY+j] += intstep*G_FIELD*v; +// phi_out[2][i*NY+j] -= intstep*G_FIELD*u; +// phi_out[1][i*NY+j] -= intstep*G_FIELD*v; +// phi_out[2][i*NY+j] += intstep*G_FIELD*u; +// phi_out[1][i*NY+j] -= intstep*G_FIELD*v*wsphere[i*NY+j].ctheta; +// phi_out[2][i*NY+j] += intstep*G_FIELD*u*wsphere[i*NY+j].ctheta; + } + else + { + phi_out[1][i*NY+j] += intstep*G_FIELD*v; + phi_out[2][i*NY+j] -= intstep*G_FIELD*u; + } } break; } @@ -1340,11 +1425,58 @@ void evolve_wave(double *phi[NFIELDS], double *phi_tmp[NFIELDS], short int xy_in evolve_wave_half(phi_tmp, phi, xy_in, potential_field, vector_potential_field, gfield, rde, wsphere); } +void update_tracer_table(double tracers[2*N_TRACERS*NSTEPS], t_rde rde[NX*NY], int time) +/* update tracer information in rde */ +{ + int tracer, t, t1, maxtime, i, j, n, ij[2], length = 50, cell, oldcell; + double x, y; + + #pragma omp parallel for private(cell) + for (cell=0; cell time) maxtime = time; + + #pragma omp parallel for private(tracer) + for (tracer = 0; tracer < N_TRACERS; tracer++) + { + for (t = 0; t < maxtime; t++) + { + t1 = time - t; + + x = tracers[t1*2*N_TRACERS + 2*tracer]; + y = tracers[t1*2*N_TRACERS + 2*tracer + 1]; + + xy_to_ij(x, y, ij); + cell = ij[0]*NY + ij[1]; + + n = rde[cell].n_tracer_pts; + if (n < NMAX_TRACER_PTS) + { + rde[cell].tracerx[n] = x; + rde[cell].tracery[n] = y; + rde[cell].n_tracer_pts++; + rde[cell].tracer = length - t; + } +// else printf("More than %i tracer points per cell\n", NMAX_TRACER_PTS); + + if ((cell != oldcell)&&(t > 0)) + rde[cell].prev_cell = oldcell; + + oldcell = cell; + } + } +} -void evolve_tracers(double *phi[NFIELDS], double tracers[2*N_TRACERS*NSTEPS], int time, int nsteps, double step) +void evolve_tracers(double *phi[NFIELDS], double tracers[2*N_TRACERS*NSTEPS], t_rde rde[NX*NY], int time, int nsteps, double step) /* time steps of tracer particle evolution (for Euler equation) */ { - int tracer, i, j, t, ij[2], iplus, jplus; + int tracer, i, j, n, t, ij[2], iplus, jplus, prev_cell, new_cell; double x, y, xy[2], vx, vy; step = TRACERS_STEP; @@ -1408,6 +1540,8 @@ void evolve_tracers(double *phi[NFIELDS], double tracers[2*N_TRACERS*NSTEPS], in tracers[(time+1)*2*N_TRACERS + 2*tracer + 1] = y; } } + + if ((PLOT_3D)&&(time+1 < NSTEPS)) update_tracer_table(tracers, rde, time); } @@ -1542,7 +1676,7 @@ void draw_color_bar_palette(int plot, double range, int palette, int circular, i if (circular) draw_circular_color_scheme_palette_fade(XMAX - 2.0*width, YMAX - 2.0*width, 1.0*width, plot, -range, range, palette, fade, fade_value); else if (ROTATE_COLOR_SCHEME) - draw_color_scheme_palette_fade(XMIN + 0.8, YMIN + 0.1, XMAX - 0.8, YMIN + 0.1 + width, plot, -range, range, palette, fade, fade_value); + draw_color_scheme_palette_fade(XMIN + 0.8, YMIN + 0.05, XMAX - 0.8, YMIN + 0.05 + width, plot, -range, range, palette, fade, fade_value); else draw_color_scheme_palette_fade(XMAX - 1.5*width, YMIN + 0.1, XMAX - 0.5*width, YMAX - 0.1, plot, -range, range, palette, fade, fade_value); } @@ -1599,35 +1733,73 @@ void viewpoint_schedule(int i) /* change position of observer */ { int j; - double angle, ca, sa; - static double observer_initial[3]; + double angle, ca, sa, r1, interpolate, rho; + static double observer_initial[3], r, ratio, rho0, zmax; static int first = 1; if (first) { for (j=0; j<3; j++) observer_initial[j] = observer[j]; + r1 = observer[0]*observer[0] + observer[1]*observer[1]; + r = sqrt(r1 + observer[2]*observer[2]); + ratio = r/sqrt(r1); + rho0 = module2(observer[0], observer[1]); + if (vabs(rho0) < 0.001) rho0 = 0.001; + zmax = r*sin(MAX_LATITUDE*PI/180.0); first = 0; } - angle = (ROTATE_ANGLE*DPI/360.0)*(double)i/(double)NSTEPS; + interpolate = (double)i/(double)NSTEPS; + angle = (ROTATE_ANGLE*DPI/360.0)*interpolate; +// printf("i = %i, interpolate = %.3lg, angle = %.3lg\n", i, interpolate, angle); ca = cos(angle); sa = sin(angle); - observer[0] = ca*observer_initial[0] - sa*observer_initial[1]; - observer[1] = sa*observer_initial[0] + ca*observer_initial[1]; + switch (VIEWPOINT_TRAJ) + { + case (VP_HORIZONTAL): + { + observer[0] = ca*observer_initial[0] - sa*observer_initial[1]; + observer[1] = sa*observer_initial[0] + ca*observer_initial[1]; + break; + } + case (VP_ORBIT): + { + observer[0] = ca*observer_initial[0] - sa*observer_initial[1]*ratio; + observer[1] = ca*observer_initial[1] + sa*observer_initial[0]*ratio; + observer[2] = ca*observer_initial[2]; + break; + } + case (VP_ORBIT2): + { + observer[0] = ca*observer_initial[0] - sa*observer_initial[1]*ratio; + observer[1] = ca*observer_initial[1] + sa*observer_initial[0]*ratio; + observer[2] = sa*zmax; + break; + } + case (VP_POLAR): + { + rho = -sa*observer_initial[2] + ca*rho0; + observer[0] = observer_initial[0]*rho/rho0; + observer[1] = observer_initial[1]*rho/rho0; + observer[2] = ca*observer_initial[2] + sa*rho0; + break; + } + } + printf("Angle %.3lg, Observer position (%.3lg, %.3lg, %.3lg)\n", angle, observer[0], observer[1], observer[2]); } void animation() { - double time = 0.0, scale, dx, var, jangle, cosj, sinj, sqrintstep, + double time = 0.0, scale, dx, var, jangle, cosj, sinj, sqrintstep, phishift, thetashift, amp, intstep0, viscosity_printed, fade_value, noise = NOISE_INTENSITY, x, y, sign, phase; double *phi[NFIELDS], *phi_tmp[NFIELDS], *potential_field, *vector_potential_field, *tracers, *gfield, *bc_field, *bc_field2; short int *xy_in; int i, j, k, s, nvid, field; static int counter = 0; t_rde *rde; - t_wave_sphere *wsphere; + t_wave_sphere *wsphere, *wsphere_hr; /* Since NX and NY are big, it seemed wiser to use some memory allocation here */ for (i=0; i= INITIAL_TIME)&&(DOUBLE_MOVIE)) { - draw_wave_rde(1, phi, xy_in, rde, wsphere, potential_field, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 0, 1.0, REFRESH_B); - if (ADD_TRACERS) draw_tracers(phi, tracers, i, 0, 1.0); + draw_wave_rde(1, phi, xy_in, rde, wsphere, wsphere_hr, potential_field, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 0, 1.0, REFRESH_B); + if ((ADD_TRACERS)&&(!PLOT_3D)) draw_tracers(phi, tracers, i, 0, 1.0); // draw_billiard(); if (PRINT_PARAMETERS) print_parameters(phi, rde, xy_in, time, PRINT_LEFT, viscosity_printed, noise); if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, CIRC_COLORBAR_B, 0, 1.0); @@ -1938,8 +2151,8 @@ void animation() { if (DOUBLE_MOVIE) { - draw_wave_rde(0, phi, xy_in, rde, wsphere, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1); - if (ADD_TRACERS) draw_tracers(phi, tracers, NSTEPS, 0, 1.0); + draw_wave_rde(0, phi, xy_in, rde, wsphere, wsphere_hr, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1); + if ((ADD_TRACERS)&&(!PLOT_3D)) draw_tracers(phi, tracers, NSTEPS, 0, 1.0); // draw_billiard(); if (PRINT_PARAMETERS) print_parameters(phi, rde, xy_in, time, PRINT_LEFT, viscosity_printed, noise); if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, 0, 1.0); @@ -1950,16 +2163,16 @@ void animation() else for (i=0; i YMAX) particles[n].yc += YMIN - YMAX; + if (particles[n].yc > INITYMAX) particles[n].yc += INITYMIN - INITYMAX; // else if (particles[n].yc < YMIN) particles[n].yc += YMAX - YMIN; particles[n].radius = MU; /* activate only circles that intersect the domain */ @@ -1033,9 +1038,15 @@ void init_particle_config(t_particle particles[NMAXCIRCLES]) printf("Function init_circle_config not defined for this pattern \n"); } } + + for (i=0; i type2) + { + max = type1; + min = type2; + } + else + { + max = type2; + min = type1; + } + + /* enzymes */ + if ((min >= 3)&&(max == 7)) return(-10); + + if (max - min >= 2) return(1); + if (min == 1) return(-4*coulomb); + if (min%2 == 1) return(-10*coulomb); + return(1); +} + int compute_particle_interaction(int i, int k, double force[2], double *torque, t_particle* particle, double distance, double krepel, double ca, double sa, double ca_rel, double sa_rel) /* compute repelling force and torque of particle #k on particle #i */ /* returns 1 if distance between particles is smaller than NBH_DIST_FACTOR*MU */ { - double x1, y1, x2, y2, r, f, angle, aniso, fx, fy, ff[2], dist_scaled, spin_f, ck, sk, ck_rel, sk_rel, alpha, amp, charge; + double x1, y1, x2, y2, r, f, angle, aniso, fx, fy, ff[2], dist_scaled, spin_f, ck, sk, ck_rel, sk_rel, alpha, amp, charge, f1; static double dxhalf = 0.5*(BCXMAX - BCXMIN), dyhalf = 0.5*(BCYMAX - BCYMIN); int wwrapx, wwrapy, twrapx, twrapy; @@ -3936,7 +4027,75 @@ int compute_particle_interaction(int i, int k, double force[2], double *torque, force[1] = f*sa; break; } - + case (I_COULOMB_PENTA): + { + charge = particle[i].charge*particle[k].charge; + if (charge < 0.0) + { + penta_lj_force(distance, ca, sa, ca_rel, sa_rel, ff, particle[k]); + force[0] = -charge*krepel*ff[0]; + force[1] = -charge*krepel*ff[1]; + } + else + { + f = krepel*lennard_jones_force(distance, particle[k]); + force[0] = f*ca; + force[1] = f*sa; + } + break; + } + case (I_COULOMB_IMAGINARY): + { + charge = particle[i].charge*particle[k].charge; + f = -100.0*krepel*charge/(1.0e-12 + distance*distance); + if (charge <= 0.0) + f1 = 0.01*krepel*lennard_jones_force(distance, particle[k]); + else f1 = 0.0; + force[0] = f1*ca - f*sa; + force[1] = f1*sa + f*ca; + break; + } + case (I_DNA_CHARGED): + { + f = 0.2*krepel*lennard_jones_force(distance, particle[k]); + charge = CHARGE*(double)dna_charge(particle[i].type, particle[k].type, particle[i].coulomb); + if ((RD_REACTION == CHEM_DNA_ENZYME)||(RD_REACTION == CHEM_DNA_ENZYME_REPAIR)) + { + /* make some interactions repulsive */ + if ((particle[i].reactive == 0)&&(particle[k].reactive == 0)&&(particle[i].added == particle[k].added)) + charge = vabs(charge); + + /* TEST */ + /* make interactions repulsive between base-paired and other molecules */ +// else if (particle[i].paired != particle[k].paired) +// charge = vabs(charge); + } + f -= krepel*charge/(1.0e-12 + distance*distance); + force[0] = f*ca; + force[1] = f*sa; + break; + } + case (I_DNA_CHARGED_B): + { + f = 0.2*krepel*lennard_jones_force(distance, particle[k]); + charge = CHARGE*(double)dna_charge(particle[i].type, particle[k].type, particle[i].coulomb); + if (particle[i].added != particle[k].added) charge *= 1.5; + if ((RD_REACTION == CHEM_DNA_ENZYME)||(RD_REACTION == CHEM_DNA_ENZYME_REPAIR)) + { + /* make some interactions repulsive */ + if ((particle[i].reactive == 0)&&(particle[k].reactive == 0)&&(particle[i].added == particle[k].added)) + charge = vabs(charge); + + /* TEST */ + /* make interactions repulsive between base-paired and other molecules */ +// else if (particle[i].paired != particle[k].paired) +// charge = vabs(charge); + } + f -= krepel*charge/(1.0e-12 + distance*distance); + force[0] = f*ca; + force[1] = f*sa; + break; + } } if (ROTATION) @@ -4166,34 +4325,50 @@ void print_partners(int i, t_particle particle[NMAXCIRCLES]) { int p; - printf("Particle %i, type %i, charge %.3lg\n", i, particle[i].type, particle[i].charge); - printf("%i partners: ", particle[i].npartners); + if (particle[i].active) + { +// printf("Particle %i, type %i, charge %.3lg\n", i, particle[i].type, particle[i].charge); + fprintf(lj_log, "Particle %i, type %i, charge %.3lg\n", i, particle[i].type, particle[i].charge); - for (p=0; p nmolecules) nmolecules = m + 1; + np = molecule[m].nparticles; + if (np < NPARTNERS+1) + { + molecule[m].particle[np] = i; + molecule[m].nparticles++; + } + molecule[m].added = particle[i].added; + } + printf("Found %i molecules\n", nmolecules); + + /* for debugging */ + for (m=0; m NMAXCIRCLES) @@ -4549,15 +5162,55 @@ void init_particle_pairs(t_particle particle[NMAXCIRCLES]) exit(1); } + nmolecules = ncircles; + for (i=0; i p) - { - particle[q].p0 = p; - particle[q].p1 = particle[p].partner[0]; - } - } - } +// for (i=0; i p) +// { +// particle[q].p0 = p; +// particle[q].p1 = particle[p].partner[0]; +// } +// } +// } printf("ncircles = %i\n", ncircles); -// for (i=0; i 0.5*(YMAX - YMIN)) y -= (YMAX - YMIN); else if (y < 0.5*(YMIN - YMAX)) y += (YMAX - YMIN); angle = argument(x, y)*MOL_ANGLE_FACTOR; +// printf("Particle p = %i, mol_angle = %i\n", p, particle.mol_angle); +// angle = argument(x, y)*(double)particle.mol_angle; +// angle = argument(x, y)*(double)(other_particle[particle.p0].npartners); while (angle > DPI) angle -= DPI; while (angle < 0.0) angle += DPI; hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*(angle)/DPI; @@ -5599,6 +6257,15 @@ int draw_special_particle(t_particle particle, double xc1, double yc1, double ra } break; } +// case (CHEM_DNA_ENZYME): +// { +// if (particle.type == 7) +// { +// /* TODO */ +// } +// return(0); +// break; +// } } return(1); } @@ -6091,6 +6758,17 @@ void draw_particles(t_particle particle[NMAXCIRCLES], int plot, double beta, t_c radius *= 0.75; break; } + case (I_COULOMB_PENTA): + { + nsides = 5; + break; + } + case (I_DNA_CHARGED): + { + if (particle[j].type == 7) nsides = 3; + else nsides = NSEG; + break; + } default: nsides = NSEG; } @@ -7702,11 +8380,99 @@ double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacl } } +void dissociate_particles_findp(int i, int j, t_particle particle[NMAXCIRCLES]) +/* dissociate partnered particles i and j */ +{ + int np, nq, r, p, q; + + np = particle[i].npartners; + nq = particle[j].npartners; + + /* find which partner of j is particle i */ + p = 0; + while ((particle[i].partner[p] != j)&&(p < np)) p++; + + /* remove partner p from partner list */ + for (q=p; q 0.5) r = 0.5; + if (r > 2.0*eq_distance) r = 2.0*eq_distance; +// if (r > 1.5*eq_distance) r = 1.5*eq_distance; ca = dx/r; sa = dy/r; - if (r < 0.5) force = KSPRING_PAIRS*(r - eq_distance); - else force = 0.0; + /* TODO: adjust max distance */ + if (r < 1.5*eq_distance) force = KSPRING_PAIRS*(r - eq_distance); + else + { +// printf("Dissociating partners %i and %i because max distance exceeded\n", j, n); + force = 0.0; +// dissociate_particles_findp(j, n, particle); +// dissociate_molecule(j, n, particle); + } f[0] = force*ca; f[1] = force*sa; + + if (ROTATION) + { + sangle = sin(SPIN_INTER_FREQUENCY*(particle[n].angle - particle[j].angle)); + if (sangle > 0.0) *torque = KTORQUE_PAIRS*(1.0 + sangle); + else *torque = KTORQUE_PAIRS*(-1.0 + sangle); + + if (COMPUTE_PAIR_TORQUE) + { + torque2 = KTORQUE_PAIR_ANGLE*sangle; + f[0] -= torque2*sa; + f[1] += torque2*ca; + } + } } @@ -7784,9 +8572,10 @@ void compute_particle_force(int j, double krepel, t_particle particle[NMAXCIRCLE if (PAIR_FORCE) for (l=0; l 0.0) + add_particle(x, y, 5.0*V_INITIAL*(double)rand()/RAND_MAX, -10.0*V_INITIAL, PARTICLE_MASS, type, particle); + else + add_particle(x, y, 5.0*V_INITIAL*(double)rand()/RAND_MAX, 10.0*V_INITIAL, PARTICLE_MASS, type, particle); particle[ncircles - 1].eq_dist = EQUILIBRIUM_DIST; particle[ncircles - 1].thermostat = 1; px[ncircles - 1] = particle[ncircles - 1].vx; py[ncircles - 1] = particle[ncircles - 1].vy; + +// init_molecule_data(particle, molecule); return (ncircles); @@ -9545,15 +10356,19 @@ int chem_merge_molecule(int i, int type2, int maxpartners, t_particle particle[N } -int chem_multi_glue_molecule(int i, int type2, int maxpartners, int require_charge, int equalize_charge, t_particle particle[NMAXCIRCLES], t_collision *collisions, int ncollisions, double reaction_prob) +int chem_multi_glue_molecule(int i, int type2, int maxpartners, int no_triangles, int no_cluster, int require_charge, int equalize_charge, t_particle particle[NMAXCIRCLES], t_molecule molecule[NMAXCIRCLES], t_collision *collisions, int ncollisions, double reaction_prob) /* glue molecule containing particle i with another molecule of type type 2 */ /* having at most nmaxpartners partners already */ +/* if no_triangles >= 1, particles that have a common partner are not paired */ +/* if no_cluster = 1, particles are not paired if they belong to the same cluster */ { - int k, p, q, p1, q1, np, nq, n, m, closeby = 0, reaction = 0, jp[NMAXPARTNERS], jq[NMAXPARTNERS], r, kk, different; - double distance, angle, move_factor, xnew, ynew, vxnew, vynew, m1, m2, mr1, mr2, deltav; - short int charge_condition; + int k, p, q, p1, q1, p2, np, nq, n, m, closeby = 0, reaction = 0, jp[NMAXPARTNERS], jq[NMAXPARTNERS], r, kk, different, np2, moli, molk, mp, newpartner, pp, type1; + double distance, angle, move_factor, xnew, ynew, vxnew, vynew, m1, m2, mr1, mr2, deltav, x, y; + short int charge_condition, triangle_condition; + type1 = particle[i].type; np = particle[i].npartners; + moli = particle[i].molecule; if (np > maxpartners) return(ncollisions); m1 = 1.0/particle[i].mass_inv; @@ -9568,8 +10383,98 @@ int chem_multi_glue_molecule(int i, int type2, int maxpartners, int require_char { k = particle[i].hashneighbour[p]; nq = particle[k].npartners; + molk = particle[k].molecule; charge_condition = (!require_charge)||(particle[i].charge*particle[k].charge < 0.0); - if ((particle[k].active)&&(charge_condition)&&(particle[k].type == type2)&&(nq <= maxpartners)) + triangle_condition = 1; + + /* exclude mergers between added particles */ + if (particle[i].added == particle[k].added) + { + triangle_condition = particle[i].reactive*particle[k].reactive; + /* allow mergers if molecules are paired to same DNA strand */ + if ((particle[i].paired)&&(particle[k].paired)&&(particle[i].partner_molecule % NGRIDY == particle[k].partner_molecule % NGRIDY)) triangle_condition = 1; + } + + if ((RD_REACTION == CHEM_DNA_ENZYME)||(RD_REACTION == CHEM_DNA_ENZYME_REPAIR)) + { + /* exclude mergers of added molecules to original ones via backbone */ + if ((particle[i].added != particle[k].added)&&(type2 <= 2)&&(type1 <= 2)) + triangle_condition = 0; + + /* NEW TEST */ +// if (molk == moli) triangle_condition = 0; + + /* NEW TEST */ + /* avoid original backbones to form a loop */ + if ((!particle[i].added)&&(!particle[k].added)&&((molk-moli > NGRIDX/2)||(moli-molk > NGRIDX/2))) triangle_condition = 0; + } + + if (no_triangles == 1) /* do not pair i and k if they have a common partner */ + { + for (q=0; q= 3) /* exclude connections of different type between the same molecules */ + { + if (molk != moli) for (p1=0; p1= 4) /* also exclude second-order partners */ + { + for (q=0; q 1) + { + n = particle[k].partner[p1]; + x = particle[i].xc - particle[n].xc; + y = particle[i].yc - particle[n].yc; + /* deal with periodic boundary conditions */ + if (x > 0.5*(XMAX - XMIN)) x -= (XMAX - XMIN); + else if (x < 0.5*(XMIN - XMAX)) x += (XMAX - XMIN); + if (y > 0.5*(YMAX - YMIN)) y -= (YMAX - YMIN); + else if (y < 0.5*(YMIN - YMAX)) y += (YMAX - YMIN); + distance = module2(x, y); + + if (distance > particle[i].radius) + { + np = particle[i].npartners; + particle[i].npartners++; + particle[i].partner[np] = n; + particle[i].partner_eqd[np] = distance; + + nq = particle[n].npartners; + particle[n].npartners++; + particle[n].partner[nq] = i; + particle[n].partner_eqd[nq] = distance; + } + } } } } @@ -9671,6 +10679,35 @@ int chem_multi_glue_molecule(int i, int type2, int maxpartners, int require_char if (reaction) { + /* TEST */ + if ((RD_REACTION == CHEM_DNA_ENZYME)||(RD_REACTION == CHEM_DNA_ENZYME_REPAIR)) + { + /* consider merged marticles as reactive */ + if ((particle[i].added == 1)&&(particle[i].type >= 3)&&(particle[i].type <= 6)) + { + particle[i].paired = 1; + particle[i].partner_molecule = molk; + for (p1 = 0; p1= 3)&&(type2 <= 6)) + { + particle[k].paired = 1; + particle[k].partner_molecule = moli; + for (p1 = 0; p1= 3)&&(particle[i].type <= 6)) + { + particle[i].paired = 0; + for (q1 = 0; q1 <= particle[i].npartners; q1++) + { + n = particle[i].partner[q1]; + particle[n].paired = 0; + } + + particle[k].paired = 0; + for (q1 = 0; q1 <= particle[k].npartners; q1++) + { + n = particle[k].partner[q1]; + particle[n].paired = 0; + } + } + } + } + + /* update molecule data */ + if ((reaction1)&&(moli != mol_diss)) + { + printf("Splitting molecules %i and %i\n", moli, mol_diss); + + mp = molecule[moli].npartners; + pp = 0; + while (molecule[moli].partner[pp] != mol_diss) pp++; + if (pp < mp) + { + while (pp < mp-1) + { + molecule[moli].partner[pp] = molecule[moli].partner[pp+1]; + molecule[moli].connection_type[pp] = molecule[moli].connection_type[pp+1]; + pp++; + } + molecule[moli].npartners--; + } + + mp = molecule[mol_diss].npartners; + pp = 0; + while (molecule[mol_diss].partner[pp] != moli) pp++; + if (pp < mp) + { + while (pp < mp-1) + { + molecule[mol_diss].partner[pp] = molecule[mol_diss].partner[pp+1]; + molecule[mol_diss].connection_type[pp] = molecule[mol_diss].connection_type[pp+1]; + pp++; + } + molecule[mol_diss].npartners--; + } +// sleep(3); + } + } + p++; + } + + if (reaction) + { + printf("dissociating molecule %i\n", moli); + collisions[ncollisions].x = particle[i].xc; + collisions[ncollisions].y = particle[i].yc; + collisions[ncollisions].time = COLLISION_TIME; + collisions[ncollisions].color = type_hue(particle[i].type); + + if (ncollisions < 2*NMAXCOLLISIONS - 1) ncollisions++; + else printf("Too many collisions\n"); + } + + return(ncollisions); +} void change_type_proportion(t_particle particle[NMAXCIRCLES], double prop) /* change proportion of particles of types 1 or 2 */ @@ -9820,13 +11101,378 @@ int update_cluster_color(t_particle particle[NMAXCIRCLES]) return(nclusters); } -int update_types(t_particle particle[NMAXCIRCLES], t_collision *collisions, int ncollisions, int *particle_numbers, int time, double *delta_e) + +int repair_dna(t_particle particle[NMAXCIRCLES], t_molecule molecule[NMAXCIRCLES], t_collision *collisions, int ncollisions) +/* repair unwanted connections in DNA, for CHEM_DNA_ENZYME_REPAIR reaction type */ +{ + int i, mol, moli, j, molj, p, q, p1, q1, molp, molpp, molq, molqq, r, type, type1, type2, type3, type4, np, nq, pp, qq, mp, split, common_partner, smol1, smol2, npartners; + int p_table[NMAXPARTNERMOLECULES], q_table[NMAXPARTNERMOLECULES]; + double dist; + + printf("Repairing DNA\n"); + /* Prevent single base from connecting to two different molecules */ + for (moli=0; moli= 3) + { + for (p=0; p= 3)&&(molp != molj)) + { + smol1 = molj; + smol2 = molp; + split = 1; + } + } + } + } + if (split) + { + printf("Splitting molecule %i from molecules %i and %i\n\n\n", moli, smol1, smol2); + + for (p=0; p= 3) + { + i = p1; + for (q=0; q= 3) + for (p=0; p= 3) + for (p=0; p= 3) + { + p_table[np] = molp; + np++; + } + } + + /* find base-neighbours for molj (normally there should be at most one) */ + nq = 0; + for (q=0; q= 3)) + if (type3 >= 3) + { + q_table[nq] = molq; + nq++; + } + } + + /* test for split, split occurs if each molecule has a base-neighbour, + and they are not neighbours */ + if ((np == 0)||(nq == 0)) split = 0; + else + { + split = 1; + + printf("Base partners of molecule %i(%i): ", moli, type1); + for (p=0; p= 1)) + if ((particle[qq].molecule == molj)&&(deltatype == 1)) + { + dist = module2(particle[pp].xc - particle[qq].xc, particle[pp].yc - particle[qq].yc); +// if (dist > 10.0*MU) +// if (dist > 15.0*MU) + if (dist > 20.0*MU) + { + split = 1; + i = pp; + fprintf(lj_log, "\n\n\n Time = %i\n", frame_time); + fprintf(lj_log, "Particles %i(%i) in molecule %i and %i(%i) in molecule %i are at distance %.3lg\n", pp, particle[pp].type, moli, qq, particle[qq].type, molj, dist); + } + } + } + } + } + + if (split) + { + printf("Splitting molecule %i from molecule %i because of wrong connection\n", moli, molj); + fprintf(lj_log, "Splitting molecule %i from molecule %i because of wrong connection\n", moli, molj); + + ncollisions = chem_split_molecule(i, particle, molecule, collisions, ncollisions); + } + } + } + + return(ncollisions); +} + + +int update_types(t_particle particle[NMAXCIRCLES], t_molecule molecule[NMAXCIRCLES], t_collision *collisions, int ncollisions, int *particle_numbers, int time, double *delta_e) /* update the types in case of reaction-diffusion equation */ { - int i, j, k, n, n3, n4, p, type, oldncollisions, delta_n; - double distance, rnd, p1, p2; + int i, j, k, n, n3, n4, p, type, atype, btype, oldncollisions, delta_n; + double distance, rnd, p1, p2, reac_dist; static double inv_masses[RD_TYPES+1], radii[RD_TYPES+1]; - static int first = 1; + static int first = 1, repair_counter = 0; if (first) /* compute total mass and mass ratios */ { @@ -10178,7 +11824,7 @@ int update_types(t_particle particle[NMAXCIRCLES], t_collision *collisions, int { if (particle[i].npartners < AGREGMAX) for (k=0; k<3; k++) { - ncollisions = chem_multi_glue_molecule(i, k, AGREGMAX, 0, 0, particle, collisions, ncollisions, REACTION_PROB); + ncollisions = chem_multi_glue_molecule(i, k, AGREGMAX, 0, 0, 0, 0, particle, molecule, collisions, ncollisions, REACTION_PROB); } /* decouple particles with several partners from thermostat */ @@ -10189,7 +11835,6 @@ int update_types(t_particle particle[NMAXCIRCLES], t_collision *collisions, int if ((PLOT == P_CLUSTER)||(PLOT_B == P_CLUSTER)) update_cluster_color(particle); - printf("%i collisions\n", ncollisions); delta_n = ncollisions - oldncollisions; printf("delta_n = %i\n", delta_n); @@ -10201,7 +11846,7 @@ int update_types(t_particle particle[NMAXCIRCLES], t_collision *collisions, int { if (particle[i].npartners < AGREGMAX) for (k=0; k<3; k++) { - ncollisions = chem_multi_glue_molecule(i, k, AGREGMAX, 1, 0, particle, collisions, ncollisions, REACTION_PROB); + ncollisions = chem_multi_glue_molecule(i, k, AGREGMAX, 0, 0, 1, 0, particle, molecule, collisions, ncollisions, REACTION_PROB); } /* decouple particles with several partners from thermostat */ @@ -10212,7 +11857,375 @@ int update_types(t_particle particle[NMAXCIRCLES], t_collision *collisions, int if ((PLOT == P_CLUSTER)||(PLOT_B == P_CLUSTER)) update_cluster_color(particle); + printf("%i collisions\n", ncollisions); + delta_n = ncollisions - oldncollisions; + printf("delta_n = %i\n", delta_n); + return(ncollisions); + } + case (CHEM_AGGREGATION_NNEIGH): + { + for (i=0; i= AGREG_DECOUPLE) particle[i].thermostat = 0; + } + /* update cluster color scheme */ + if ((PLOT == P_CLUSTER)||(PLOT_B == P_CLUSTER)) + update_cluster_color(particle); + + printf("%i collisions\n", ncollisions); + delta_n = ncollisions - oldncollisions; + printf("delta_n = %i\n", delta_n); + return(ncollisions); + } + case (CHEM_DNA): + { + for (i=0; i= AGREG_DECOUPLE) particle[i].thermostat = 0; + } + + /* update cluster color scheme */ + if ((PLOT == P_CLUSTER)||(PLOT_B == P_CLUSTER)) + update_cluster_color(particle); + + printf("%i collisions\n", ncollisions); + delta_n = ncollisions - oldncollisions; + printf("delta_n = %i\n", delta_n); + return(ncollisions); + } + case (CHEM_DNA_ALT): + { + for (i=0; i 0) + ncollisions = chem_multi_glue_molecule(i, btype, AGREGMAX, 0, 0, 0, 0, particle, molecule, collisions, ncollisions, REACTION_PROB); + } + + /* decouple particles with several partners from thermostat */ + if (particle[i].npartners >= AGREG_DECOUPLE) particle[i].thermostat = 0; + } + + /* update cluster color scheme */ +// if ((PLOT == P_CLUSTER)||(PLOT_B == P_CLUSTER)) + update_cluster_color(particle); + + printf("%i collisions\n", ncollisions); + delta_n = ncollisions - oldncollisions; + printf("delta_n = %i\n", delta_n); + return(ncollisions); + } + case (CHEM_DNA_DOUBLE): + { + for (i=0; i 0) + ncollisions = chem_multi_glue_molecule(i, btype, AGREGMAX, 0, 0, 0, 0, particle, molecule, collisions, ncollisions, REACTION_PROB); + } + + /* decouple particles with several partners from thermostat */ + if (particle[i].npartners >= AGREG_DECOUPLE) particle[i].thermostat = 0; + } + + /* update cluster color scheme */ + if ((PLOT == P_CLUSTER)||(PLOT_B == P_CLUSTER)) + update_cluster_color(particle); + + printf("%i collisions\n", ncollisions); + delta_n = ncollisions - oldncollisions; + printf("delta_n = %i\n", delta_n); + return(ncollisions); + } + case (CHEM_DNA_DSPLIT): + { + for (i=0; i 0) + ncollisions = chem_multi_glue_molecule(i, btype, AGREGMAX, 3, 0, 0, 0, particle, molecule, collisions, ncollisions, REACTION_PROB); + } + + /* split molecule with a certain probability */ + if ((double)rand()/RAND_MAX < DISSOCIATION_PROB) + ncollisions = chem_split_molecule(i, particle, molecule, collisions, ncollisions); + + /* decouple particles with several partners from thermostat */ + if (particle[i].npartners >= AGREG_DECOUPLE) particle[i].thermostat = 0; + } + + /* update cluster color scheme */ + if ((PLOT == P_CLUSTER)||(PLOT_B == P_CLUSTER)) + update_cluster_color(particle); + + printf("%i collisions\n", ncollisions); + delta_n = ncollisions - oldncollisions; + printf("delta_n = %i\n", delta_n); + return(ncollisions); + } + case (CHEM_DNA_BASE_SPLIT): + { + for (i=0; i 0) + ncollisions = chem_multi_glue_molecule(i, btype, AGREGMAX, 3, 0, 0, 0, particle, molecule, collisions, ncollisions, REACTION_PROB); + } + + /* split molecule with a certain probability */ + if ((double)rand()/RAND_MAX < DISSOCIATION_PROB) + ncollisions = chem_split_molecule(i, particle, molecule, collisions, ncollisions); + + reac_dist = 2.5*MU; + + switch (particle[i].type) { + case (3): + { + ncollisions = chem_split_molecule_if_nearby(i, 3, particle, molecule, collisions, ncollisions, reac_dist, REACTION_PROB); + ncollisions = chem_split_molecule_if_nearby(i, 5, particle, molecule, collisions, ncollisions, reac_dist, REACTION_PROB); + ncollisions = chem_split_molecule_if_nearby(i, 6, particle, molecule, collisions, ncollisions, reac_dist, REACTION_PROB); + break; + } + case (4): + { + ncollisions = chem_split_molecule_if_nearby(i, 4, particle, molecule, collisions, ncollisions, reac_dist, REACTION_PROB); + ncollisions = chem_split_molecule_if_nearby(i, 5, particle, molecule, collisions, ncollisions, reac_dist, REACTION_PROB); + ncollisions = chem_split_molecule_if_nearby(i, 6, particle, molecule, collisions, ncollisions, reac_dist, REACTION_PROB); + break; + } + case (5): + { + ncollisions = chem_split_molecule_if_nearby(i, 3, particle, molecule, collisions, ncollisions, reac_dist, REACTION_PROB); + ncollisions = chem_split_molecule_if_nearby(i, 3, particle, molecule, collisions, ncollisions, reac_dist, REACTION_PROB); + ncollisions = chem_split_molecule_if_nearby(i, 5, particle, molecule, collisions, ncollisions, reac_dist, REACTION_PROB); + break; + } + case (6): + { + ncollisions = chem_split_molecule_if_nearby(i, 3, particle, molecule, collisions, ncollisions, reac_dist, REACTION_PROB); + ncollisions = chem_split_molecule_if_nearby(i, 4, particle, molecule, collisions, ncollisions, reac_dist, REACTION_PROB); + ncollisions = chem_split_molecule_if_nearby(i, 6, particle, molecule, collisions, ncollisions, reac_dist, REACTION_PROB); + break; + } + } + + /* decouple particles with several partners from thermostat */ + if (particle[i].npartners >= AGREG_DECOUPLE) particle[i].thermostat = 0; + } + + /* update cluster color scheme */ + if ((PLOT == P_CLUSTER)||(PLOT_B == P_CLUSTER)) + update_cluster_color(particle); + + printf("%i collisions\n", ncollisions); + delta_n = ncollisions - oldncollisions; + printf("delta_n = %i\n", delta_n); + return(ncollisions); + } + case (CHEM_DNA_ENZYME): + { + for (i=0; i 0) + ncollisions = chem_multi_glue_molecule(i, btype, AGREGMAX, 3, 0, 0, 0, particle, molecule, collisions, ncollisions, REACTION_PROB); + } + + /* split molecule with a certain probability */ + /* TODO update molecules after split */ + if ((double)rand()/RAND_MAX < DISSOCIATION_PROB) + ncollisions = chem_split_molecule(i, particle, molecule, collisions, ncollisions); + + reac_dist = 2.5*MU; + + switch (particle[i].type) { + case (3): + { + ncollisions = chem_local_split_molecule_if_nearby(i, 4, 7, particle, molecule, collisions, ncollisions, reac_dist, REACTION_PROB); + break; + } + case (4): + { + ncollisions = chem_local_split_molecule_if_nearby(i, 3, 7, particle, molecule, collisions, ncollisions, reac_dist, REACTION_PROB); + break; + } + case (5): + { + ncollisions = chem_local_split_molecule_if_nearby(i, 6, 7, particle, molecule, collisions, ncollisions, reac_dist, REACTION_PROB); + break; + } + case (6): + { + ncollisions = chem_local_split_molecule_if_nearby(i, 5, 7, particle, molecule, collisions, ncollisions, reac_dist, REACTION_PROB); + break; + } + case (7): + { + if ((double)rand()/RAND_MAX < KILLING_PROB) particle[i].active = 0; + break; + } + } + + /* decouple particles with several partners from thermostat */ + if (particle[i].npartners >= AGREG_DECOUPLE) particle[i].thermostat = 0; + } + + /* update cluster color scheme */ + if ((PLOT == P_CLUSTER)||(PLOT_B == P_CLUSTER)) + update_cluster_color(particle); + + /* for debugging */ +// for (i=0; i 0) +// { +// printf("Molecule %i has %i partners: ", i, molecule[i].npartners); +// for (j=0; j 0) /* TEST */ + ncollisions = chem_multi_glue_molecule(i, btype, AGREGMAX, 3, 0, 0, 0, particle, molecule, collisions, ncollisions, REACTION_PROB); +// ncollisions = chem_multi_glue_molecule(i, btype, AGREGMAX, 3, 0, 0, 2, particle, molecule, collisions, ncollisions, REACTION_PROB); + } + + /* split molecule with a certain probability */ + /* TODO update molecules after split */ + if ((double)rand()/RAND_MAX < DISSOCIATION_PROB) + ncollisions = chem_split_molecule(i, particle, molecule, collisions, ncollisions); + + reac_dist = 2.5*MU; + + switch (particle[i].type) { + case (3): + { + ncollisions = chem_local_split_molecule_if_nearby(i, 4, 7, particle, molecule, collisions, ncollisions, reac_dist, REACTION_PROB); + break; + } + case (4): + { + ncollisions = chem_local_split_molecule_if_nearby(i, 3, 7, particle, molecule, collisions, ncollisions, reac_dist, REACTION_PROB); + break; + } + case (5): + { + ncollisions = chem_local_split_molecule_if_nearby(i, 6, 7, particle, molecule, collisions, ncollisions, reac_dist, REACTION_PROB); + break; + } + case (6): + { + ncollisions = chem_local_split_molecule_if_nearby(i, 5, 7, particle, molecule, collisions, ncollisions, reac_dist, REACTION_PROB); + break; + } + case (7): + { + if ((double)rand()/RAND_MAX < KILLING_PROB) particle[i].active = 0; + break; + } + } + + /* decouple particles with several partners from thermostat */ + if (particle[i].npartners >= AGREG_DECOUPLE) particle[i].thermostat = 0; + } + +// ncollisions = repair_dna(particle, molecule, collisions, ncollisions); + + repair_counter++; + if (repair_counter >= 2) + { + ncollisions = repair_dna(particle, molecule, collisions, ncollisions); +// ncollisions = check_dna_pairing(particle, molecule, collisions, ncollisions); + repair_counter = 0; + } + + /* update cluster color scheme */ + if ((PLOT == P_CLUSTER)||(PLOT_B == P_CLUSTER)) + update_cluster_color(particle); + + /* for debugging */ +// for (i=0; i 0) +// { +// printf("Molecule %i has %i partners: ", i, molecule[i].npartners); +// for (j=0; j 0.0) sign = 1.0; + else sign = -1.0; + + if (add) + { + phi[1][i*NY+j] += sign*module; + phi[0][i*NY+j] -= sign*module; /* approximate, stream function should solve Poisson equation */ + } + else + { + phi[1][i*NY+j] = sign*module; + phi[0][i*NY+j] = -sign*module; /* approximate, stream function should solve Poisson equation */ + } + } + else + { + phi[0][i*NY+j] = 0.0; + phi[1][i*NY+j] = 0.0; + } + } + break; + } + case (E_EULER_COMP): + { + for (i=0; i 0.0) sign = 1.0; + else sign = -1.0; + +// phi[0][i*NY+j] = 1.0; + /* nonconstant density to make things more interesting */ + /* TODO: correct initialization of phi[1], phi[2] */ + if (add) + { +// phi[0][i*NY+j] += 0.5 + density_mod*module/amp; + phi[0][i*NY+j] += density_mod*module/amp; + phi[1][i*NY+j] -= sign*module*(xy[1]-theta0)/vabs(scale); + phi[2][i*NY+j] += sign*module*pscal1/vabs(scale); + for (k=1; k<2; k++) phi[k][i*NY+j] *= wsphere[i*NY+j].stheta; + /* approximate, stream function should solve Poisson equation */ + } + else + { + phi[0][i*NY+j] = 1.0 + density_mod*module/amp; + phi[1][i*NY+j] = -sign*module*(xy[1]-theta0)/vabs(scale); + phi[2][i*NY+j] = sign*module*pscal1/vabs(scale); + for (k=1; k<2; k++) phi[k][i*NY+j] *= wsphere[i*NY+j].stheta; + } + + } + else + { + phi[0][i*NY+j] = 1.0; + phi[1][i*NY+j] = 0.0; + phi[2][i*NY+j] = 0.0; + } + } + break; + } + } +} + + +void init_vortex_state_sphere_mod(int add, double amp, double phi0, double theta0, double scale, double density_mod, double *phi[NFIELDS], short int xy_in[NX*NY], t_wave_sphere wsphere[NX*NY]) +/* initialise field with vortex at position (x,y) with amplitude and variance scale */ +/* for incompressible Euler, phi[0] is stream function, phi[1] is vorticity */ +/* for compressible Euler, phi[0] is the density, phi[1] and phi[2] are velocity components */ +{ + int i, j, k, ij[2]; + double xy[2], dist, dist2, module, phase, scale2, sign, pscal, pscal1, pscal2; + double x0, y0, z0, x, y, z, a, b, c, ephi, etheta; + +// scale2 = scale*scale; + switch (RDE_EQUATION) { + case (E_EULER_INCOMP): + { + for (i=0; i 0.0) sign = 1.0; + else sign = -1.0; + + if (add) + { + phi[1][i*NY+j] += sign*module; + phi[0][i*NY+j] -= sign*module; /* approximate, stream function should solve Poisson equation */ + } + else + { + phi[1][i*NY+j] = sign*module; + phi[0][i*NY+j] = -sign*module; /* approximate, stream function should solve Poisson equation */ + } + } + else + { + phi[0][i*NY+j] = 0.0; + phi[1][i*NY+j] = 0.0; + } + } + break; + } + case (E_EULER_COMP): + { + xy_to_ij(XMIN + (XMAX-XMIN)*phi0/DPI, YMIN + (YMAX-YMIN)*theta0/PI, ij); + x0 = wsphere[ij[0]*NY+ij[1]].x; + y0 = wsphere[ij[0]*NY+ij[1]].y; + z0 = wsphere[ij[0]*NY+ij[1]].z; + +// x0 = cos(phi0)*sin(theta0); +// y0 = sin(phi0)*sin(theta0); +// z0 = cos(theta0); + for (i=0; i 0.0) sign = 1.0; + else sign = -1.0; + +// phi[0][i*NY+j] = 1.0; + /* nonconstant density to make things more interesting */ + /* TODO: correct initialization of phi[1], phi[2] */ + if (add) + { +// phi[0][i*NY+j] += 0.5 + density_mod*module/amp; + phi[0][i*NY+j] += density_mod*module/amp; + phi[1][i*NY+j] += sign*module*dist*ephi/vabs(scale); + phi[2][i*NY+j] += sign*module*dist*etheta/vabs(scale); + for (k=1; k<2; k++) phi[k][i*NY+j] *= wsphere[i*NY+j].stheta; + /* approximate, stream function should solve Poisson equation */ + } + else + { + phi[0][i*NY+j] = 1.0 + density_mod*module/amp; + phi[1][i*NY+j] = sign*module*dist*ephi/vabs(scale); + phi[2][i*NY+j] = sign*module*dist*etheta/vabs(scale); + for (k=1; k<2; k++) phi[k][i*NY+j] *= wsphere[i*NY+j].stheta; + } + + } + else + { + phi[0][i*NY+j] = 1.0; + phi[1][i*NY+j] = 0.0; + phi[2][i*NY+j] = 0.0; + } + } + break; + } + } +} + void init_shear_flow(double amp, double delta, double rho, int nx, int ny, double yshift, double *phi[NFIELDS], short int xy_in[NX*NY]) /* initialise field with a shear flow */ /* phi[0] is stream function, phi[1] is vorticity */ @@ -790,6 +1017,64 @@ void init_laminar_flow(double amp, double xmodulation, double ymodulation, doubl } } +void init_laminar_flow_earth(double amp, double *phi[NFIELDS], short int xy_in[NX*NY]) +/* initialise field with a laminar flow in x direction */ +/* for Earth simulation with trades and westerlies */ +/* amp is global amplitude */ +{ + int i, j; + double xy[2], y1; + + switch (RDE_EQUATION) { + case (E_EULER_INCOMP): /* TODO */ + { + for (i=0; i max) + { + max = phi[k][i*NY+j]; + kmax = k; + } + } + angle = colors_rps(kmax); + rde[i*NY+j].theta = angle; + } + else rde[i*NY+j].theta = 0.0; + } + break; + } case (Z_MAXTYPE_RPSLZ): { #pragma omp parallel for private(i,j,max,kmax) @@ -1798,24 +2148,39 @@ void compute_gradient_euler_plane(double phi[NX*NY], double gradient[2*NX*NY], d void compute_gradient_euler_sphere(double phi[NX*NY], double gradient[2*NX*NY], t_wave_sphere wsphere[NX*NY], double yshift) /* compute the gradient of the field */ { - int i, j, iplus, iminus, jplus, jminus, padding = 0; - double deltaphi, maxgradient = 1.0e10, sintheta, invstheta; + int i, j, iplus, iminus, jplus, jminus, n, b, k, p, q, i1, b0; + double deltaphi, maxgradient = 1.0e10, sintheta, invstheta, sum1, sum2, factor; // double dx = (XMAX-XMIN)/((double)NX); - static short int first = 1; - static double dphi, dtheta, cphiphi, ctheta, dt, vdrift; + static short int first = 1, jsouth, jnorth; + static double dphi, dtheta, cphiphi, ctheta, dt, vdrift, xyfactor; if (first) { dphi = DPI/(double)NX; dtheta = PI/(double)NY; + xyfactor = 1.0/(dtheta*(double)NX); + + if (SMOOTHBLOCKS) + { + jsouth = BLOCKDIST; + jnorth = NY - BLOCKDIST; + } + else + { + jsouth = 1; + jnorth = NY-1; + } first = 0; } -// dx = (XMAX-XMIN)/((double)NX); - +// #pragma omp parallel for private(i) +// for (i=0; i<2*NX*NY; i++) gradient[i] = 0.0; + +// printf("Computing gradient Euler sphere\n"); + #pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus) - for (j=1; j n*b) + { + factor = 1.0/(NX - n*b); + sum1 = 0.0; + sum2 = 0.0; + for (i=n*b; iNY-BLOCKDIST-1; j--) + { + b = block_sizes[j]; + n = block_numbers[j]; + factor = 1.0/(double)b; + + for (k=0; k n*b) + { + factor = 1.0/(NX - n*b); + sum1 = 0.0; + sum2 = 0.0; + for (i=n*b; iNY-BLOCKDIST-1-b; j--) + { + for (i=0; i n*b) + { + factor = 1.0/(NX - n*b); + sum1 = 0.0; + sum2 = 0.0; + sum3 = 0.0; + sum4 = 0.0; + for (i=n*b; i=jnorth; j--) + { + b = block_sizes[j]; + n = block_numbers[j]; + factor = 1.0/(double)b; + + for (k=0; k n*b) + { + factor = 1.0/(NX - n*b); + sum1 = 0.0; + sum2 = 0.0; + sum3 = 0.0; + sum4 = 0.0; + for (i=n*b; i n*b) + { + factor = 1.0/(NX - n*b); + sum = 0.0; + for (i=n*b; iNY-BLOCKDIST-1; j--) + { + b = block_sizes[j]; + n = block_numbers[j]; + factor = 1.0/(double)b; + + for (k=0; k n*b) + { + factor = 1.0/(NX - n*b); + sum = 0.0; + for (i=n*b; i RMAX) r = RMAX; if (r < RMIN) r = RMIN; wsphere[i*NY+j].radius = r; @@ -3305,7 +4415,8 @@ void compute_light_angle_sphere_rde(short int xy_in[NX*NY], t_rde rde[NX*NY], t_ /* TODO ? Avoid artifacts due to singularity at north pole */ for (j=NY - DPOLE; j RMAX) r = RMAX; // if (r < RMIN) r = RMIN; if (r < 1.0) r = 1.0; @@ -3313,7 +4424,8 @@ void compute_light_angle_sphere_rde(short int xy_in[NX*NY], t_rde rde[NX*NY], t_ } for (j=0; j RMAX) r = RMAX; // if (r < RMIN) r = RMIN; if (r < 1.0) r = 1.0; @@ -3355,7 +4467,19 @@ void compute_light_angle_sphere_rde(short int xy_in[NX*NY], t_rde rde[NX*NY], t_ } else { - deltar = (wsphere[(i+1)*NY+j].radius - wsphere[i*NY+j].radius)/dphi; + if (SMOOTHBLOCKS) + { + if ((j > BLOCKDIST)&&(j < NY-1-BLOCKDIST)) + deltar = (wsphere[(i+1)*NY+j].radius - wsphere[i*NY+j].radius)/dphi; + else + { + b = block_sizes[j]; + iplus = i + b; + if (iplus > NX) iplus -= NX; + deltar = (wsphere[iplus*NY+j].radius - wsphere[i*NY+j].radius)/((double)b*dphi); + } + } + else deltar = (wsphere[(i+1)*NY+j].radius - wsphere[i*NY+j].radius)/dphi; deltai[0] = -wsphere[i*NY+j].radius*wsphere[i*NY+j].sphi; deltai[0] += deltar*wsphere[i*NY+j].cphi; @@ -3493,7 +4617,7 @@ void compute_light_angle_sphere_rde(short int xy_in[NX*NY], t_rde rde[NX*NY], t_ void compute_light_angle_sphere_rde_2d(short int xy_in[NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY], double potential[NX*NY], int movie, int transparent) /* computes cosine of angle between normal vector and vector light */ { - int i, j; + int i, j, b, iplus, iminus; short int draw; double gradx, grady, norm, pscal; static double dx, dy, vscale2, vshift; @@ -3517,10 +4641,25 @@ void compute_light_angle_sphere_rde_2d(short int xy_in[NX*NY], t_rde rde[NX*NY], { if (xy_in[i*NY+j]) { + if (SMOOTHBLOCKS) { - gradx = (*rde[(i+1)*NY+j].p_zfield[movie] - *rde[(i-1)*NY+j].p_zfield[movie])/dx; - grady = (*rde[i*NY+j+1].p_zfield[movie] - *rde[i*NY+j-1].p_zfield[movie])/dy; + if ((j > BLOCKDIST)&&(j < NY-1-BLOCKDIST)) + { + gradx = (*rde[(i+1)*NY+j].p_zfield[movie] - *rde[(i-1)*NY+j].p_zfield[movie])/dx; + } + else + { + b = block_sizes[j]; + iplus = i + b; + if (iplus >= NX) iplus -= NX; + iminus = i - b; + if (iminus < 0) iminus += NX; + gradx = (*rde[iplus*NY+j].p_zfield[movie] - *rde[iminus*NY+j].p_zfield[movie])/dx; + } } + else gradx = (*rde[(i+1)*NY+j].p_zfield[movie] - *rde[(i-1)*NY+j].p_zfield[movie])/dx; + + grady = (*rde[i*NY+j+1].p_zfield[movie] - *rde[i*NY+j-1].p_zfield[movie])/dy; norm = sqrt(vscale2 + gradx*gradx + grady*grady); pscal = -gradx*light[0] - grady*light[1] + SHADE_SCALE_2D; @@ -3535,7 +4674,22 @@ void compute_light_angle_sphere_rde_2d(short int xy_in[NX*NY], t_rde rde[NX*NY], { if (xy_in[j]) { - gradx = (*rde[NY+j].p_zfield[movie] - *rde[(NY-1)*NY+j].p_zfield[movie])/dx; + if (SMOOTHBLOCKS) + { + if ((j > BLOCKDIST)&&(j < NY-1-BLOCKDIST)) + { + gradx = (*rde[NY+j].p_zfield[movie] - *rde[(NX-1)*NY+j].p_zfield[movie])/dx; + } + else + { + b = block_sizes[j]; + iplus = b; + iminus = NX - b; + gradx = (*rde[iplus*NY+j].p_zfield[movie] - *rde[iminus*NY+j].p_zfield[movie])/dx; + } + } + else gradx = (*rde[NY+j].p_zfield[movie] - *rde[(NY-1)*NY+j].p_zfield[movie])/dx; + grady = (*rde[j+1].p_zfield[movie] - *rde[j-1].p_zfield[movie])/dy; norm = sqrt(vscale2 + gradx*gradx + grady*grady); @@ -3549,10 +4703,24 @@ void compute_light_angle_sphere_rde_2d(short int xy_in[NX*NY], t_rde rde[NX*NY], /* i=NX-1 */ for (j=1; j BLOCKDIST)&&(j < NY-1-BLOCKDIST)) + { + gradx = (*rde[NY+j].p_zfield[movie] - *rde[(NY-1)*NY+j].p_zfield[movie])/dx; + } + else + { + b = block_sizes[j]; + iplus = b - 1; + iminus = NX - 1 - b; + gradx = (*rde[iplus*NY+j].p_zfield[movie] - *rde[iminus*NY+j].p_zfield[movie])/dx; + } + } + else gradx = (*rde[NY+j].p_zfield[movie] - *rde[(NY-1)*NY+j].p_zfield[movie])/dx; + grady = (*rde[j+1].p_zfield[movie] - *rde[j-1].p_zfield[movie])/dy; norm = sqrt(vscale2 + gradx*gradx + grady*grady); @@ -3725,6 +4893,11 @@ void compute_field_color_rde(double value, int cplot, int palette, double rgb[3] hsl_to_rgb_palette(360.0*(1.0 - vabs(color_amplitude(value, 1.0, 0))), 0.9, 0.5, rgb, palette); break; } + case (Z_MAXTYPE_RPS): + { + hsl_to_rgb_palette(value, 0.9, 0.5, rgb, palette); + break; + } case (Z_MAXTYPE_RPSLZ): { hsl_to_rgb_palette(value, 0.9, 0.5, rgb, palette); @@ -3740,6 +4913,11 @@ void compute_field_color_rde(double value, int cplot, int palette, double rgb[3] color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 0, rgb); break; } + case (Z_NORM_GRADIENT_RPSLZ_ASYM): + { + color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 0, rgb); + break; + } case (Z_MODULE): { color_scheme_asym_palette(COLOR_SCHEME, palette, VSCALE_AMPLITUDE*value, 1.0, 0, rgb); @@ -3807,12 +4985,14 @@ void compute_field_color_rde(double value, int cplot, int palette, double rgb[3] } case (Z_EULER_DIRECTION): { - hsl_to_rgb_palette(360.0*value/DPI, 0.9, 0.5, rgb, palette); + if (SPHERE) hsl_to_rgb_palette(360.0*(1.0 - value/DPI), 0.9, 0.5, rgb, palette); + else hsl_to_rgb_palette(360.0*value/DPI, 0.9, 0.5, rgb, palette); break; } case (Z_EULER_DIRECTION_SPEED): { - hsl_to_rgb_palette(360.0*value/DPI, 0.9, 0.5, rgb, palette); + if (SPHERE) hsl_to_rgb_palette(360.0*(1.0 - value/DPI), 0.9, 0.5, rgb, palette); + else hsl_to_rgb_palette(360.0*value/DPI, 0.9, 0.5, rgb, palette); break; } case (Z_SWATER_HEIGHT): @@ -3978,7 +5158,7 @@ double compute_depth_colors_rde(int i, int j, t_rde rde[NX*NY], double potential } -void compute_rde_fields(double *phi[NFIELDS], short int xy_in[NX*NY], int zplot, int cplot, t_rde rde[NX*NY]) +void compute_rde_fields(double *phi[NFIELDS], short int xy_in[NX*NY], int zplot, int cplot, t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY]) /* compute the necessary auxiliary fields */ { int i, j; @@ -3988,18 +5168,19 @@ void compute_rde_fields(double *phi[NFIELDS], short int xy_in[NX*NY], int zplot, { if ((COMPUTE_THETA)||(COMPUTE_THETAZ)) compute_theta(phi, xy_in, rde); - + if ((zplot == Z_NORM_GRADIENT)||(cplot == Z_NORM_GRADIENT)) { compute_gradient_theta(rde); // compute_gradient_polar(rde, 1.0); - compute_gradient_polar(rde, 0.003); +// compute_gradient_polar(rde, 0.001); + compute_gradient_polar(rde, ZSCALE_NORMGRADIENT); } if ((zplot == Z_NORM_GRADIENTX)||(cplot == Z_NORM_GRADIENTX)||(zplot == Z_ANGLE_GRADIENTX)||(cplot == Z_ANGLE_GRADIENTX)) { compute_gradient_rde(phi[0], rde); - compute_gradient_polar(rde, 0.03); + compute_gradient_polar(rde, ZSCALE_NORMGRADIENT); } if ((zplot == Z_VORTICITY)||(cplot == Z_VORTICITY)||(zplot == Z_VORTICITY_ABS)||(cplot == Z_VORTICITY_ABS)) @@ -4007,12 +5188,16 @@ void compute_rde_fields(double *phi[NFIELDS], short int xy_in[NX*NY], int zplot, compute_gradient_theta(rde); compute_curl(rde); } + + if ((zplot == Z_MAXTYPE_RPS)||(cplot == Z_MAXTYPE_RPS)) + compute_theta_rpslz(phi, xy_in, rde, cplot); + break; } case (E_RPSLZ): { compute_theta_rpslz(phi, xy_in, rde, cplot); - if ((zplot == Z_NORM_GRADIENT_RPSLZ)||(cplot == Z_NORM_GRADIENT_RPSLZ)) + if ((zplot == Z_NORM_GRADIENT_RPSLZ)||(cplot == Z_NORM_GRADIENT_RPSLZ)||(zplot == Z_NORM_GRADIENT_RPSLZ_ASYM)||(cplot == Z_NORM_GRADIENT_RPSLZ_ASYM)) { compute_gradient_theta(rde); compute_gradient_polar(rde, 0.005); @@ -4039,7 +5224,7 @@ void compute_rde_fields(double *phi[NFIELDS], short int xy_in[NX*NY], int zplot, case (E_EULER_COMP): { if ((zplot == Z_EULER_SPEED)||(cplot == Z_EULER_SPEED)||(zplot == Z_EULER_DIRECTION_SPEED)||(cplot == Z_EULER_DIRECTION_SPEED)) - compute_speed(phi, rde); + compute_speed(phi, rde, wsphere); if ((zplot == Z_EULER_DIRECTION)||(cplot == Z_EULER_DIRECTION)||(zplot == Z_EULER_DIRECTION_SPEED)||(cplot == Z_EULER_DIRECTION_SPEED)) compute_direction(phi, rde); if ((zplot == Z_EULERC_VORTICITY)||(cplot == Z_EULERC_VORTICITY)) @@ -4050,7 +5235,7 @@ void compute_rde_fields(double *phi[NFIELDS], short int xy_in[NX*NY], int zplot, { if (zplot == Z_SWATER_HEIGHT) adjust_height(phi, rde); if ((zplot == Z_SWATER_SPEED)||(cplot == Z_SWATER_SPEED)||(zplot == Z_SWATER_DIRECTION_SPEED)||(cplot == Z_SWATER_DIRECTION_SPEED)) - compute_speed(phi, rde); + compute_speed(phi, rde, wsphere); if ((zplot == Z_SWATER_DIRECTION_SPEED)||(cplot == Z_SWATER_DIRECTION_SPEED)) compute_direction(phi, rde); } @@ -4118,6 +5303,12 @@ void init_zfield_rde(double *phi[NFIELDS], short int xy_in[NX*NY], int zplot, t_ for (i=0; i NX) ii -= NX; + + glVertex2i(HRES*ii, HRES*j); + glVertex2i(HRES*(ii+1), HRES*j); + glVertex2i(HRES*(ii+1), HRES*(j+1)); + glVertex2i(HRES*ii, HRES*(j+1)); } glEnd (); + + /* draw the continents */ + if (RDE_PLANET) + { + glBegin(GL_QUADS); + for (i=0; i HRES*NX) ii -= HRES*NX; + + glVertex2i(ii, j); + glVertex2i(ii+1, j); + glVertex2i(ii+1, j+1); + glVertex2i(ii, j+1); + } + } + glEnd (); + } + if (DRAW_BILLIARD) draw_billiard(0, 1.0); } -void draw_wave_sphere_rde_ij(int i, int iplus, int j, int jplus, int jcolor, int movie, double *phi[NFIELDS], short int xy_in[NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY], int zplot, int cplot, int palette, int fade, double fade_value) +void draw_wave_sphere_rde_ij(int i, int iplus, int j, int jplus, int jcolor, int movie, double *phi[NFIELDS], short int xy_in[NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY], t_wave_sphere wsphere_hr[HRES*HRES*NX*NY], int zplot, int cplot, int palette, int fade, double fade_value) /* draw wave at simulation grid point (i,j) */ { - int k, l; - short int draw, notdraw, draw_bc=1; - double xyz[3], ca; + int k, l, n, m, s, p, q, i1, j1, prev_cell, cell, iplus1, jplus1; + short int draw, drawij, notdraw, draw_bc=1, draw_hr[HRES*HRES]; + double xyz[3], ca, rgb[3], r_hr[HRES*HRES], g_hr[HRES*HRES], b_hr[HRES*HRES], lfactor, ratio, xt[NMAX_TRACER_PTS+1], yt[NMAX_TRACER_PTS+1], zt[NMAX_TRACER_PTS+1]; if (NON_DIRICHLET_BC) draw_bc = (xy_in[i*NY+j])&&(xy_in[iplus*NY+j])&&(xy_in[i*NY+jplus])&&(xy_in[iplus*NY+jplus]); -// if (FLOODING) draw = wsphere[i*NY+j].draw_wave; -// else draw = ((TWOSPEEDS)||(xy_in[i*NY+j])); -// else - draw = (xy_in[i*NY+j]); - if (draw) glColor3f(rde[i*NY+jcolor].rgb[0], rde[i*NY+jcolor].rgb[1], rde[i*NY+jcolor].rgb[2]); - else - { - ca = wsphere[i*NY+j].cos_angle; - ca = (ca + 1.0)*0.4 + 0.2; - if (fade) ca *= fade_value; -// 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); - } -// if (FLOODING) -// { -// glBegin(GL_TRIANGLE_FAN); -// if (ij_to_sphere(i, j, *wave[i*NY+j].p_zfield[movie], wsphere, xyz, wsphere[i*NY+j].draw_wave)) -// draw_vertex_sphere(xyz); -// if (ij_to_sphere(iplus, j, *wave[iplus*NY+j].p_zfield[movie], wsphere, xyz, wsphere[iplus*NY+j].draw_wave)) -// draw_vertex_sphere(xyz); -// if (ij_to_sphere(iplus, jplus, *wave[iplus*NY+j+1].p_zfield[movie], wsphere, xyz, wsphere[iplus*NY+j+1].draw_wave)) -// draw_vertex_sphere(xyz); -// if (ij_to_sphere(i, jplus, *wave[i*NY+j+1].p_zfield[movie], wsphere, xyz, wsphere[i*NY+j+1].draw_wave)) -// draw_vertex_sphere(xyz); -// glEnd (); -// } -// else + draw = wsphere[i*NY+j].indomain; + for (p=0; p 0)) + { + prev_cell = rde[cell].prev_cell; + lfactor = (double)rde[cell].tracer/50.0; + for (k=0; k<3; k++) rgb[k] *= 1.0 + lfactor; + glColor3f(rgb[0], rgb[1], rgb[2]); + n = rde[cell].n_tracer_pts; + + #pragma omp parallel for private(k) + for (k=0; k= 0)&&(prev_cell != cell)) + { + xt[n] = rde[prev_cell].tracerx[m]; + yt[n] = rde[prev_cell].tracery[m]; + zt[n] = wsphere[prev_cell].radius; + n++; + } + + /* do some smoothing of the curve */ + #pragma omp parallel for private(k) + for (k=0; kimid; i--) { for (j=0; j<=jmax; j++) - draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, zplot, cplot, palette, fade, fade_value); + draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); } for (i=imid; i>imin; i--) { for (j=0; j<=jmid; j++) - draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, zplot, cplot, palette, fade, fade_value); + draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); for (j=jmax; j>=jmid; j--) - draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, zplot, cplot, palette, fade, fade_value); + draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); } for (i=imax+1; i=jmid; j--) - draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, zplot, cplot, palette, fade, fade_value); + draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); + } + + if (imin >= 1) + { +// for (j=0; j<=jmid; j++) +// draw_wave_sphere_rde_ij(imin-1, imin, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); + for (j=jmax; j>=jmid; j--) + draw_wave_sphere_rde_ij(imin-1, imin, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); } } else @@ -4604,74 +5951,88 @@ void draw_wave_sphere_3d_rde(int movie, double *phi[NFIELDS], short int xy_in[NX for (i=imax; i=jmid; j--) - draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, zplot, cplot, palette, fade, fade_value); + draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); } for (i=imax-1; i>=0; i--) { for (j=0; j<=jmax; j++) - draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, zplot, cplot, palette, fade, fade_value); + draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); } for (j=0; j<=jmax; j++) - draw_wave_sphere_rde_ij(NX-1, 0, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, zplot, cplot, palette, fade, fade_value); + draw_wave_sphere_rde_ij(NX-1, 0, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); for (i=NX-2; i>=imin; i--) { for (j=0; j<=jmid; j++) - draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, zplot, cplot, palette, fade, fade_value); + draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); for (j=jmax; j>=jmid; j--) - draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, zplot, cplot, palette, fade, fade_value); + draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); } + + if (imin >= 1) + { + /* experimental */ + for (j=jmid/3; j<=jmid; j++) + draw_wave_sphere_rde_ij(imin-1, imin, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); + for (j=jmax; j>=jmid; j--) + draw_wave_sphere_rde_ij(imin-1, imin, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); + } + } /* North pole */ for (i=0; iimid; i--) { for (j=jmax; j>=0; j--) - draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, zplot, cplot, palette, fade, fade_value); + draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); } for (i=imid; i>imin; i--) { for (j=jmax; j>=jmid; j--) - draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, zplot, cplot, palette, fade, fade_value); + draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); for (j=0; j<=jmid; j++) - draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, zplot, cplot, palette, fade, fade_value); + draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); } for (i=imax+1; i=0; j--) - draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, zplot, cplot, palette, fade, fade_value); + draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); } for (j=jmax; j>=0; j--) - draw_wave_sphere_rde_ij(NX-1, 0, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, zplot, cplot, palette, fade, fade_value); + draw_wave_sphere_rde_ij(NX-1, 0, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); for (i=0; i<=imin; i++) { for (j=jmax; j>=jmid; j--) - draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, zplot, cplot, palette, fade, fade_value); + draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); for (j=0; j<=jmid; j++) - draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, zplot, cplot, palette, fade, fade_value); + draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); } } else @@ -4679,42 +6040,44 @@ void draw_wave_sphere_3d_rde(int movie, double *phi[NFIELDS], short int xy_in[NX for (i=imax; i=0; j--) - draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, zplot, cplot, palette, fade, fade_value); + draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); } for (i=imid; i=jmid; j--) - draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, zplot, cplot, palette, fade, fade_value); + draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); for (j=0; j<=jmid; j++) - draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, zplot, cplot, palette, fade, fade_value); + draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); } for (i=imax-1; i>=0; i--) { for (j=jmax; j>=0; j--) - draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, zplot, cplot, palette, fade, fade_value); + draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); } for (j=jmax; j>=0; j--) - draw_wave_sphere_rde_ij(NX-1, 0, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, zplot, cplot, palette, fade, fade_value); + draw_wave_sphere_rde_ij(NX-1, 0, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); for (i=NX-2; i>=imin-1; i--) { for (j=jmax; j>=jmid; j--) - draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, zplot, cplot, palette, fade, fade_value); + draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); for (j=0; j<=jmid; j++) - draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, zplot, cplot, palette, fade, fade_value); + draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); } } /* South pole */ for (i=0; i0; j--) - draw_wave_sphere_rde_ij(i, i+1, j-1, j, j, movie, phi, xy_in, rde, wsphere, zplot, cplot, palette, fade, fade_value); + draw_wave_sphere_rde_ij(i, i+1, j-1, j, j, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); for (j=2; j>0; j--) - draw_wave_sphere_rde_ij(NX-1, 0, j-1, j, DPOLE, movie, phi, xy_in, rde, wsphere, zplot, cplot, palette, fade, fade_value); + draw_wave_sphere_rde_ij(NX-1, 0, j-1, j, DPOLE, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); } + + } void draw_wave_3d_ij_rde(int i, int j, int movie, double *phi[NFIELDS], short int xy_in[NX*NY], t_rde rde[NX*NY], double potential[NX*NY], int zplot, int cplot, int palette, int fade, double fade_value) @@ -4954,8 +6317,9 @@ void draw_periodicised_wave_3d(int movie, double *phi[NFIELDS], short int xy_in[ void draw_wave_rde(int movie, double *phi[NFIELDS], short int xy_in[NX*NY], t_rde rde[NX*NY], - t_wave_sphere wsphere[NX*NY], double potential[NX*NY], int zplot, - int cplot, int palette, int fade, double fade_value, int refresh) + t_wave_sphere wsphere[NX*NY], t_wave_sphere wsphere_hr[HRES*HRES*NX*NY], + double potential[NX*NY], int zplot, int cplot, int palette, int fade, + double fade_value, int refresh) { int i, j, k, l, draw = 1; double xy[2], xy_screen[2], rgb[3], pos[2], ca, rgb_e[3], rgb_w[3], rgb_n[3], rgb_s[3]; @@ -4966,7 +6330,7 @@ void draw_wave_rde(int movie, double *phi[NFIELDS], short int xy_in[NX*NY], t_rd if (refresh) { // printf("Computing fields\n"); - compute_rde_fields(phi, xy_in, zplot, cplot, rde); + compute_rde_fields(phi, xy_in, zplot, cplot, rde, wsphere); // printf("Computed fields\n"); if ((PLOT_3D)&&(SHADE_3D)) { @@ -4982,46 +6346,50 @@ void draw_wave_rde(int movie, double *phi[NFIELDS], short int xy_in[NX*NY], t_rd if (PLOT_3D) { if (PLOT_SPHERE) - draw_wave_sphere_3d_rde(movie, phi, xy_in, rde, wsphere, potential, zplot, cplot, palette, fade, fade_value, 0); + draw_wave_sphere_3d_rde(movie, phi, xy_in, rde, wsphere, wsphere_hr, potential, zplot, cplot, palette, fade, fade_value, 0); else if (DRAW_PERIODICISED) draw_periodicised_wave_3d(movie, phi, xy_in, rde, potential, zplot, cplot, palette, fade, fade_value); else draw_wave_3d_rde(movie, phi, xy_in, rde, potential, zplot, cplot, palette, fade, fade_value); } - else draw_wave_2d_rde(xy_in, rde); + else draw_wave_2d_rde(xy_in, rde, wsphere, wsphere_hr, fade, fade_value); } void draw_tracers(double *phi[NFIELDS], double tracers[2*N_TRACERS*NSTEPS], int time, int fade, double fade_value) /* draw trajectories of tracers */ { - int tracer, t, t1, length = 50; + int tracer, t, t1, length = 50, ij[2]; double x1, y1, x2, y2, lum; glColor3f(1.0, 1.0, 1.0); glLineWidth(1); + printf("Drawing tracers\n"); + t1 = time - length; if (t1 < 1) t1 = 1; for (t = t1 + 1; t < time; t++) for (tracer = 0; tracer < N_TRACERS; tracer++) { - x1 = tracers[2*(t-1)*N_TRACERS + 2*tracer]; + x1 = -tracers[2*(t-1)*N_TRACERS + 2*tracer]; y1 = tracers[2*(t-1)*N_TRACERS + 2*tracer + 1]; - x2 = tracers[2*t*N_TRACERS + 2*tracer]; + x2 = -tracers[2*t*N_TRACERS + 2*tracer]; y2 = tracers[2*t*N_TRACERS + 2*tracer + 1]; lum = 1.0 - 0.75*(double)(time - t)/(double)length; if (fade) lum *= fade_value; glColor3f(lum, lum, lum); - if (module2(x2 - x1, y2 - y1) < 0.2) draw_line(x1, y1, x2, y2); + + if (module2(x2 - x1, y2 - y1) < 0.2) draw_line_hres(x1, y1, x2, y2); // printf("time = %i, tracer = %i, coord = %i, x1 = %.2lg, y1 = %.2lg, x2 = %.2lg, y2 = %.2lg\n", t, tracer,2*t*N_TRACERS + 2*tracer, x1, y1, x2, y2); } } + void compute_average_speeds(double *phi[NFIELDS], t_rde rde[NX*NY], double *speed1, double *speed2) { int i, j; @@ -5220,12 +6588,6 @@ void smooth_poles_half(double phi_in[NX*NY], double phi_out[NX*NY]) } smooth_row(j, phi_in, phi_out); } - -// #pragma omp parallel for private(i) -// for (i=0; i= 0) + { + block_sizes[j] = b; + block_numbers[j] = NX/b; + if (j*b < BLOCKDIST) b*=2; + j--; + } + + /* poles */ + block_sizes[0] = NX; + block_numbers[0] = 1; + + for (j=0; j n*b) + { + factor = 1.0/(NX - n*b); + value = 0.0; + for (i=n*b; iNY-BLOCKDIST-1; j--) + { + b = block_sizes[j]; + n = block_numbers[j]; + factor = 1.0/(double)b; + + for (k=0; k n*b) + { + factor = 1.0/(NX - n*b); + value = 0.0; + for (i=n*b; i 0.0) /* primary mirror */ + { + if (y1 < MU) return(1); + if (x > 0.5*d + WALL_WIDTH) return(1); + return((x-x01)*(x-x01) > a1*a1*(1.0 + y*y/(b1*b1))); + } + else /* secondary mirror */ + { + if (y1 > 0.3) return(1); +// if (x < -0.5*d - WALL_WIDTH) return(1); + return((x > x02)||((x-x02)*(x-x02) < a2*a2*(1.0 + y*y/(b2*b2)))); + } +} + int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_circle *circles) /* returns 1 if (x,y) represents a point in the billiard */ { @@ -3420,7 +3515,7 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_ if ((vabs(x) < LAMBDA)&&(vabs(y - h) < 0.5*WALL_WIDTH)) return(2); return(0); } - case (D_RITCHEY_CHRETIEN): + case (D_RITCHEY_CHRETIEN_SPHERICAL): { /* LAMBDA is magnification M */ d = 2.5; @@ -3443,6 +3538,15 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_ return(module2(x-x1, y1) < r1); } } + case (D_RITCHEY_CHRETIEN_HYPERBOLIC): + { + return(rc_hyp(x, y)); + } + case (D_GRADIENT_INDEX_LENS): + { + /* use IOR_GRADIENT_INDEX_LENS for index of refraction control */ + return(1); + } case (D_MENGER): { x1 = 0.5*(x+1.0); @@ -3650,6 +3754,80 @@ void hex_transfo(double u, double v, double *x, double *y) *y = rb*v; } +void draw_rc_hyp() +/* draw D_RITCHEY_CHRETIEN_HYPERBOLIC domain */ +{ + static int first = 1; + static double m, d, dprime, b, r1, r2, e1, e2, a1, b1, a2, b2, x01, x02, dy; + double x, y, f, f1, f2; + + if (first) + { + m = LAMBDA; /* secondary magnification */ + d = 2.5; /* distance between mirrors */ + b = 3.0; /* distance between secondary mirror and effective focal point */ + + f = m*d + b; + f1 = f/m; /* focal distance of primary mirror */ + dprime = f1 - d; /* distance between secondary mirror and primary focal point */ + f2 = m*dprime/(m+1); /* focal distance of secondary mirror */ + + r1 = 2.0*f/m; /* radius of curvature of primary mirror */ + r2 = 2.0*b/(m - 1.0); /* radius of curvature of secondary mirror */ + e1 = sqrt(1.0 + 2.0*b/(m*m*m*d)); /* eccentricity of primary mirror */ + e2 = sqrt(1.0 + 2.0*(m*(2.0*m - 1.0) + b/d)/((m-1)*(m-1)*(m-1))); + + a1 = f1/e1; + b1 = sqrt(a1*r1); /* semi-axes of primary mirror */ + x01 = 0.5*d + a1; /* center of primary hyperbola */ + + a2 = f2/e2; + b2 = sqrt(a2*r2); /* semi-axes of secondary mirror */ + x02 = -0.5*d + a2; /* center of secondary hyperbola */ + + dy = (YMAX - YMIN)/(double)NSEG; + first = 0; + } + + /* primary mirror */ + glBegin(GL_LINE_STRIP); + draw_vertex(0.5*d + WALL_WIDTH, YMAX); + draw_vertex(0.5*d + WALL_WIDTH, MU); + for (y = MU; y < YMAX; y += dy) + { + x = x01 - a1*sqrt(1.0 + y*y/(b1*b1)); + draw_vertex(x, y); + } + glEnd(); + + glBegin(GL_LINE_STRIP); + draw_vertex(0.5*d + WALL_WIDTH, YMIN); + draw_vertex(0.5*d + WALL_WIDTH, -MU); + for (y = -MU; y > YMIN; y -= dy) + { + x = x01 - a1*sqrt(1.0 + y*y/(b1*b1)); + draw_vertex(x, y); + } + glEnd(); + + /* secondary mirror */ + glBegin(GL_LINE_STRIP); + draw_vertex(XMIN, -0.3); + for (y = -0.3; y < 0.3; y += dy) + { + x = x02 - a2*sqrt(1.0 + y*y/(b2*b2)); + draw_vertex(x, y); + } + draw_vertex(XMIN, 0.3); + glEnd(); + + /* focal plane */ + glBegin(GL_LINE_STRIP); + draw_vertex(b - 0.5*d, YMIN); + draw_vertex(b - 0.5*d, YMAX); + glEnd(); +} + void draw_billiard(int fade, double fade_value) /* draws the billiard boundary */ { @@ -5411,7 +5589,7 @@ void draw_billiard(int fade, double fade_value) /* draws the billiard bound write_text(pos[0], pos[1], message); break; } - case (D_RITCHEY_CHRETIEN): + case (D_RITCHEY_CHRETIEN_SPHERICAL): { d = 2.5; b = 3.5; @@ -5437,6 +5615,21 @@ void draw_billiard(int fade, double fade_value) /* draws the billiard bound draw_line(0.5*d + WALL_WIDTH, -MU, 0.5*d + WALL_WIDTH, YMIN); break; } + case (D_RITCHEY_CHRETIEN_HYPERBOLIC): + { + draw_rc_hyp(); + break; + } + case (D_GRADIENT_INDEX_LENS): + { + draw_line(-LAMBDA, YMIN, -LAMBDA, YMAX); + draw_line(LAMBDA, YMIN, LAMBDA, YMAX); + draw_line(-LAMBDA, -1.0, LAMBDA, -1.0); + draw_line(-LAMBDA, 1.0, LAMBDA, 1.0); + /* focal plane */ + draw_line(WAVE_PROFILE_X, YMIN, WAVE_PROFILE_X, YMAX); + break; + } case (D_MENGER): { glLineWidth(3); @@ -6769,7 +6962,7 @@ void compute_laplacian(double phi[NX*NY], t_laplacian laplace[NX*NY], double del double oscillating_bc(int time, int j) { int ij[2], jmin, jmax; - double t, phase, a, envelope, omega, amp, dist2; + double t, t1, phase, a, envelope, omega, amp, dist2; switch (OSCILLATION_SCHEDULE) { @@ -6793,11 +6986,18 @@ double oscillating_bc(int time, int j) case (OSC_WAVE_PACKET): { t = (double)time*OMEGA; -// a = 10.0; -// a = 0.02/OMEGA; a = sqrt(INITIAL_VARIANCE)/OMEGA; phase = AMPLITUDE*cos(t); - envelope = exp(-(t-0.2)*(t-0.2)/(a*a))*sqrt(DPI/a); + envelope = exp(-(t-INITIAL_SHIFT)*(t-INITIAL_SHIFT)/(a*a))*sqrt(DPI/a); + return(phase*envelope); + } + case (OSC_WAVE_PACKETS): + { + t = (double)time*OMEGA; + t1 = t - (double)((int)(t/WAVE_PACKET_SHIFT + 0.5))*WAVE_PACKET_SHIFT; + a = sqrt(INITIAL_VARIANCE)/OMEGA; + phase = AMPLITUDE*cos(t); + envelope = exp(-(t1-INITIAL_SHIFT)*(t1-INITIAL_SHIFT)/(a*a))*sqrt(DPI/a); return(phase*envelope); } case (OSC_CHIRP): @@ -6857,6 +7057,23 @@ double oscillating_bc(int time, int j) amp *= exp(-dist2/INITIAL_VARIANCE); return(amp); } + case (OSC_TWO_WAVES): + { + t = (double)time*OMEGA + (double)(NY-j)*OSCIL_LEFT_YSHIFT/(double)NY; + t1 = (double)time*OMEGA + (double)(j)*OSCIL_LEFT_YSHIFT/(double)NY; + if (j < NY/2) amp = cos(t); + else amp = cos(t1); + return(AMPLITUDE*amp); + } + case (OSC_TWO_WAVES_ADDED): + { + t = (double)time*OMEGA + (double)(NY-j)*OSCIL_LEFT_YSHIFT/(double)NY; + t1 = (double)time*OMEGA + (double)(j)*OSCIL_LEFT_YSHIFT/(double)NY; + amp = 0.0; + if (t > 0.0) amp += cos(t); + if (t1 > 0.0) amp += cos(t1); + return(AMPLITUDE*amp); + } } } @@ -7256,6 +7473,7 @@ void init_ior_2d(short int *xy_in[NX], double *tcc_table[NX], double *tgamma_tab { printf("Initializing IOR_LENS_WALL\n"); for (i=0; i LAMBDA)) + { + tcc_table[i][j] = courant2; + tgamma_table[i][j] = GAMMA; + } + else if (y > 1.0) + { + tcc_table[i][j] = 0.0; + tgamma_table[i][j] = 0.0; + } + else + { + tcc_table[i][j] = courantb2/(1.0 - MU*y*y); + tgamma_table[i][j] = GAMMAB; + } + } + } + break; + } + case (IOR_GRADIENT_INDEX_LENS_B): + { + /* focal distance is f = 1/(4*LAMBDA*n0*MU) */ + /* with n0 = COURANT/COURANTB */ + for (i=0; i LAMBDA)) + { + tcc_table[i][j] = courant2; + tgamma_table[i][j] = GAMMA; + } + else if (y > 1.0) + { + tcc_table[i][j] = 0.0; + tgamma_table[i][j] = 0.0; + } + else + { + speed = 1.0/(1.0 - MU*y*y); + speed *= speed; + tcc_table[i][j] = courantb2*speed; + tgamma_table[i][j] = GAMMAB; + } + } + } + break; + } + case (IOR_LINEAR_X_A): + { + /* Warning: Depending on COURANT and COURANTB */ + /* this may generate a wave speed of the form |a - bx| */ + /* Use IOR_LINEAR_X_B instead to avoid this */ + for (i=0; i LAMBDA)) + { + tcc[i*NY+j] = courant2; + tgamma[i*NY+j] = GAMMA; + } + else if (y > 1.0) + { + tcc[i*NY+j] = 0.0; + tgamma[i*NY+j] = 0.0; + } + else + { + tcc[i*NY+j] = courantb2/(1.0 - MU*y*y); + tgamma[i*NY+j] = GAMMAB; + } + } + } + break; + } + case (IOR_GRADIENT_INDEX_LENS_B): + { + /* focal distance if f = 1/(4*LAMBDA*n0*MU) */ + /* with n0 = COURANT/COURANTB */ + for (i=0; i LAMBDA)) + { + tcc[i*NY+j] = courant2; + tgamma[i*NY+j] = GAMMA; + } + else if (y > 1.0) + { + tcc[i*NY+j] = 0.0; + tgamma[i*NY+j] = 0.0; + } + else + { + speed = 1.0/(1.0 - MU*y*y); + speed *= speed; + tcc[i*NY+j] = courantb2*speed; + tgamma[i*NY+j] = GAMMAB; + } + } + } + break; + } default: { 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 0) hmin = rgbval; + } + else hmin = rgbval; + } + if (rgbval > hmax) hmax = rgbval; + } + + printf("hmin = %i, hmax = %i, hsea = %i\n", hmin, hmax, (int)hsea); + + if (B_DOMAIN == D_SPHERE_VENUS) + { + hsum = 0; + for (j=0; j nx-1) ii -= nx; + if (ii < 0) ii = 0; + jj = (int)(ry*(double)(rny-nshift - j)); + if (jj > ny-1) jj = ny-1; + if (jj < 0) jj = 0; + height_values[i*rny+j] = ((double)rgb_values[3*(jj*nx+ii)]-hsea)*cratio; + wsphere[i*rny+j].altitude = ((double)rgb_values[3*(jj*nx+ii)]-hsea)*cratio; + + /* take care of black areas (missing data) on venus */ + if ((B_DOMAIN == D_SPHERE_VENUS)&&(rgb_values[3*(jj*nx+ii)] == 0)) + { + height_values[i*rny+j] = VENUS_NODATA_FACTOR*hmean*cratio; + wsphere[i*rny+j].altitude = VENUS_NODATA_FACTOR*hmean*cratio; + } + + if (OTHER_PLANET) + wsphere[i*rny+j].indomain = (wsphere[i*rny+j].altitude < vshift); + +// if (wsphere[i*rny+j].indomain) printf("rgb = %i, altitude = %.3lg\n", rgb_values[3*(jj*nx+ii)], height_values[i*rny+j]); + } + + /* smooth values in case of high resolution */ + if ((rnx > nx)||(rny > ny)) + { + for (i=1; i nmaxpixels) + { + printf("Image too large, increase nmaxpixels in init_earth_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)(res*NY-nshift - j)); + if (jj > ny-1) jj = ny-1; + if (jj < 0) jj = 0; + wsphere[i*res*NY+j].r = (double)rgb_values[3*(jj*nx+ii)]*cratio; + wsphere[i*res*NY+j].g = (double)rgb_values[3*(jj*nx+ii)+1]*cratio; + wsphere[i*res*NY+j].b = (double)rgb_values[3*(jj*nx+ii)+2]*cratio; + +// printf("RGB at (%i, %i) = (%.3lg, %3.lg, %.3lg)\n", i, j, wsphere[i*NY+j].r, wsphere[i*NY+j].g, wsphere[i*NY+j].b); + + /* decide which points are in the Sea */ + diff = iabs(rgb_values[3*(jj*nx+ii)] - 10); + diff += iabs(rgb_values[3*(jj*nx+ii)+1] - 10); + diff += iabs(rgb_values[3*(jj*nx+ii)+2] - 51); + wsphere[i*res*NY+j].indomain = (diff < 15); + } + + /* smooth colors in case of high resolution */ + if ((res*NX > nx)||(res*NY > ny)) + for (i=1; i COS_VISIBLE); } +int ij_to_sphere_hres(int i, int j, double r, t_wave_sphere wsphere[HRES*HRES*NX*NY], double xyz[3], int use_wave_radius) +/* convert spherical to rectangular coordinates */ +{ + double pscal, newr; + static double norm_observer; + static int first = 1; + + if (first) + { + norm_observer = sqrt(observer[0]*observer[0] + observer[1]*observer[1] + observer[2]*observer[2]); + first = 0; + } + + xyz[0] = wsphere[i*HRES*NY+j].x; + xyz[1] = wsphere[i*HRES*NY+j].y; + xyz[2] = wsphere[i*HRES*NY+j].z; + + pscal = xyz[0]*observer[0] + xyz[1]*observer[1] + xyz[2]*observer[2]; + + if (use_wave_radius) + { + newr = wsphere[i*HRES*NY+j].radius; + xyz[0] *= newr; + xyz[1] *= newr; + xyz[2] *= newr; + } + else + { + newr = wsphere[i*HRES*NY+j].radius_dem; + xyz[0] *= newr; + xyz[1] *= newr; + xyz[2] *= newr; + } + + return(pscal/norm_observer > COS_VISIBLE); +} + void xyz_to_xy(double x, double y, double z, double xy_out[2]) { int i; @@ -165,6 +727,36 @@ void xyz_to_xy(double x, double y, double z, double xy_out[2]) } } +void draw_vertex_in_spherical_coords(double x, double y, double r, int i) +{ + double phi, theta, x1, y1, z1, xy_screen[2]; + static double phi_ratio, theta_ratio, phi_offset, theta_offset; + static int first = 1; + + if (first) + { + phi_ratio = DPI/(XMAX - XMIN); + theta_ratio = PI/(YMAX - YMIN); + phi_offset = phi_ratio*XMIN; + theta_offset = theta_ratio*YMIN; + first = 0; + } + + phi = phi_ratio*x - phi_offset; + theta = theta_ratio*y - theta_offset; +// phi = DPI*(x - XMIN)/(XMAX - XMIN); +// theta = PI*(y - YMIN)/(YMAX - YMIN); + x1 = r*cos(phi)*sin(theta); + y1 = r*sin(phi)*sin(theta); + z1 = -r*cos(theta); + +// printf("(phi, theta) = (%.5lg, %.5lg)\n", phi, theta); +// printf("(x1, y1, z1) = (%.5lg, %.5lg, %.5lg)\n", x1, y1, z1); + + xyz_to_xy(x1, y1, z1, xy_screen); + glVertex2d(xy_screen[0], xy_screen[1]); +} + int xy_in_billiard_sphere(int i, int j, t_wave_sphere wsphere[NX*NY]) /* returns 1 if (x,y) represents a point in the billiard */ { @@ -1871,6 +2463,18 @@ void draw_circular_color_scheme_palette_3d(double x1, double y1, double radius, color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb); break; } + case (Z_EULER_DIRECTION): + { + value = dy_phase*(double)(j); + color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb); + break; + } + case (Z_EULER_DIRECTION_SPEED): + { + value = 0.5 - dy_phase*(double)(j); + color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb); + break; + } case (Z_EULER_VORTICITY): { value = min + 1.0*dy*(double)(j); diff --git a/wave_3d.c b/wave_3d.c index 5848240..7a8ba4c 100644 --- a/wave_3d.c +++ b/wave_3d.c @@ -44,7 +44,7 @@ #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 DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */ #define SAVE_MEMORY 1 /* set to 1 to save memory when writing tiff images */ #define NO_EXTRA_BUFFER_SWAP 1 /* some OS require one less buffer swap when recording images */ @@ -52,47 +52,23 @@ #define WINWIDTH 1920 /* window width */ #define WINHEIGHT 1150 /* window height */ -// // #define NX 1920 /* number of grid points on x axis */ +// #define NX 1920 /* number of grid points on x axis */ // #define NY 1150 /* number of grid points on y axis */ #define NX 3000 /* number of grid points on x axis */ #define NY 1600 /* number of grid points on y axis */ -// #define NX 3840 /* number of grid points on x axis */ -// #define NY 2300 /* number of grid points on y axis */ -// #define XMIN -2.0 -// #define XMAX 2.0 /* x interval */ -// #define YMIN -1.197916667 -// #define YMAX 1.197916667 /* y interval for 9/16 aspect ratio */ -#define XMIN -1.669565217 -#define XMAX 1.669565217 /* x interval */ -#define YMIN -1.0 -#define YMAX 1.0 /* y interval for 9/16 aspect ratio */ +#define XMIN -2.0 +#define XMAX 2.0 /* x interval */ +#define YMIN -1.197916667 +#define YMAX 1.197916667 /* y interval for 9/16 aspect ratio */ -#define HIGHRES 0 /* set to 1 if resolution of grid is double that of displayed image */ - -// #define WINWIDTH 1280 /* window width */ -// #define WINHEIGHT 720 /* window height */ -// -// // #define NX 1280 /* number of grid points on x axis */ -// // #define NX 720 /* number of grid points on x axis */ -// // #define NY 720 /* number of grid points on y axis */ -// #define NX 2560 /* number of grid points on x axis */ -// // #define NX 1440 /* number of grid points on x axis */ -// #define NY 1440 /* number of grid points on y axis */ -// -// // #define NX 360 /* number of grid points on x axis */ -// // #define NY 360 /* number of grid points on y axis */ -// -// #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 HIGHRES 1 /* set to 1 if resolution of grid is double that of displayed image */ #define JULIA_SCALE 0.8 /* scaling for Julia sets */ /* Choice of the billiard table */ -#define B_DOMAIN 17 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN 76 /* choice of domain shape, see list in global_pdes.c */ #define CIRCLE_PATTERN 2 /* pattern of circles or polygons, see list in global_pdes.c */ @@ -100,8 +76,8 @@ #define B_DOMAIN_B 20 /* second domain shape, for comparisons */ #define CIRCLE_PATTERN_B 0 /* second pattern of circles or polygons */ -#define VARIABLE_IOR 0 /* set to 1 for a variable index of refraction */ -#define IOR 9 /* choice of index of refraction, see list in global_pdes.c */ +#define VARIABLE_IOR 1 /* set to 1 for a variable index of refraction */ +#define IOR 181 /* choice of index of refraction, see list in global_pdes.c */ #define IOR_TOTAL_TURNS 1.0 /* total angle of rotation for IOR_PERIODIC_WELLS_ROTATING */ #define MANDEL_IOR_SCALE -0.05 /* parameter controlling dependence of IoR on Mandelbrot escape speed */ @@ -109,11 +85,11 @@ #define NPOISSON 1000 /* number of points for Poisson C_RAND_POISSON arrangement */ #define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */ -#define LAMBDA 3.0 /* parameter controlling the dimensions of domain */ -#define MU 0.14 /* parameter controlling the dimensions of domain */ +#define LAMBDA 1.0 /* parameter controlling the dimensions of domain */ +#define MU 0.35 /* parameter controlling the dimensions of domain */ #define NPOLY 6 /* number of sides of polygon */ #define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */ -#define MDEPTH 2 /* depth of computation of Menger gasket */ +#define MDEPTH 7 /* depth of computation of Menger gasket */ #define MRATIO 3 /* ratio defining Menger gasket */ #define MANDELLEVEL 2000 /* iteration level for Mandelbrot set */ #define MANDELLIMIT 20.0 /* limit value for approximation of Mandelbrot set */ @@ -140,17 +116,19 @@ /* Physical parameters of wave equation */ #define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */ -#define OSCILLATE_LEFT 1 /* set to 1 to add oscilating boundary condition on the left */ +#define OSCILLATE_LEFT 0 /* set to 1 to add oscilating boundary condition on the left */ #define OSCILLATE_TOPBOT 0 /* set to 1 to enforce a planar wave on top and bottom boundary */ #define OSCILLATION_SCHEDULE 3 /* oscillation schedule, see list in global_pdes.c */ #define OSCIL_YMAX 0.35 /* defines oscillation range */ +#define INITIAL_SHIFT 20.0 /* time shift of initial wave packet (in oscillation periods) */ +#define WAVE_PACKET_SHIFT 200.0 /* time shift between wave packets (in oscillation periods) */ -#define OMEGA 0.015 /* frequency of periodic excitation */ -#define AMPLITUDE 1.0 /* amplitude of periodic excitation */ -#define ACHIRP 0.2 /* acceleration coefficient in chirp */ -#define DAMPING 0.0 /* damping of periodic excitation */ -#define COURANT 0.1 /* Courant number */ -#define COURANTB 0.01 /* Courant number in medium B */ +#define OMEGA 0.025 /* frequency of periodic excitation */ +#define AMPLITUDE 0.5 /* amplitude of periodic excitation */ +#define ACHIRP 0.2 /* acceleration coefficient in chirp */ +#define DAMPING 0.0 /* damping of periodic excitation */ +#define COURANT 0.1 /* Courant number */ +#define COURANTB 0.08 /* Courant number in medium B */ #define GAMMA 0.0 /* damping factor in wave equation */ #define GAMMAB 1.0e-6 /* damping factor in wave equation */ #define GAMMA_SIDES 1.0e-4 /* damping factor on boundary */ @@ -164,8 +142,8 @@ /* Increasing COURANT speeds up the simulation, but decreases accuracy */ /* For similar wave forms, COURANT^2*GAMMA should be kept constant */ -#define ADD_OSCILLATING_SOURCE 0 /* set to 1 to add an oscillating wave source */ -#define OSCILLATING_SOURCE_PERIOD 30 /* period of oscillating source */ +#define ADD_OSCILLATING_SOURCE 1 /* set to 1 to add an oscillating wave source */ +#define OSCILLATING_SOURCE_PERIOD 9 /* period of oscillating source */ #define ALTERNATE_OSCILLATING_SOURCE 1 /* set to 1 to alternate sign of oscillating source */ #define ADD_WAVE_PACKET_SOURCES 0 /* set to 1 to add several sources emitting wave packets */ @@ -175,17 +153,14 @@ /* Boundary conditions, see list in global_pdes.c */ -// #define B_COND 1 -#define B_COND 3 +#define B_COND 2 #define PRECOMPUTE_BC 0 /* set to 1 to compute neighbours for Laplacian in advance */ /* Parameters for length and speed of simulation */ -#define NSTEPS 1000 /* number of frames of movie */ -// #define NSTEPS 300 /* number of frames of movie */ -#define NVID 20 /* number of iterations between images displayed on screen */ -// #define NVID 10 /* number of iterations between images displayed on screen */ +#define NSTEPS 1800 /* number of frames of movie */ +#define NVID 15 /* number of iterations between images displayed on screen */ #define NSEG 1000 /* number of segments of boundary */ #define INITIAL_TIME 0 /* time after which to start saving frames */ #define BOUNDARY_WIDTH 2 /* width of billiard boundary */ @@ -202,19 +177,18 @@ /* Parameters of initial condition */ -#define INITIAL_AMP 0.75 /* amplitude of initial condition */ -// #define INITIAL_VARIANCE 0.000025 /* variance of initial condition */ -#define INITIAL_VARIANCE 0.00005 /* variance of initial condition */ -#define INITIAL_WAVELENGTH 0.05 /* wavelength of initial condition */ +#define INITIAL_AMP 1.0 /* amplitude of initial condition */ +#define INITIAL_VARIANCE 0.00001 /* variance of initial condition */ +#define INITIAL_WAVELENGTH 0.025 /* wavelength of initial condition */ /* Plot type, see list in global_pdes.c */ -#define ZPLOT 108 /* wave height */ -#define CPLOT 108 /* color scheme */ +#define ZPLOT 103 /* wave height */ +#define CPLOT 103 /* color scheme */ -#define ZPLOT_B 103 -#define CPLOT_B 103 /* plot type for second movie */ +#define ZPLOT_B 108 +#define CPLOT_B 108 /* plot type for second movie */ #define CHANGE_LUMINOSITY 1 /* set to 1 to let luminosity depend on energy flux intensity */ #define FLUX_WINDOW 30 /* size of averaging window of flux intensity */ @@ -240,12 +214,11 @@ #define ROTATE_VIEW 0 /* set to 1 to rotate position of observer */ #define ROTATE_ANGLE 360.0 /* total angle of rotation during simulation */ -// #define ROTATE_ANGLE 45.0 /* total angle of rotation during simulation */ /* Color schemes */ -#define COLOR_PALETTE 13 /* Color palette, see list in global_pdes.c */ -#define COLOR_PALETTE_B 17 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE 17 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE_B 14 /* Color palette, see list in global_pdes.c */ #define BLACK 1 /* background */ @@ -254,11 +227,11 @@ #define SCALE 0 /* set to 1 to adjust color scheme to variance of field */ #define SLOPE 1.0 /* sensitivity of color on wave amplitude */ #define VSCALE_AMPLITUDE 2.0 /* additional scaling factor for color scheme P_3D_AMPLITUDE */ -#define VSCALE_ENERGY 1.5 /* additional scaling factor for color scheme P_3D_ENERGY */ +#define VSCALE_ENERGY 3.0 /* additional scaling factor for color scheme P_3D_ENERGY */ #define PHASE_FACTOR 20.0 /* factor in computation of phase in color scheme P_3D_PHASE */ #define PHASE_SHIFT 0.0 /* shift of phase in color scheme P_3D_PHASE */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ -#define E_SCALE 50.0 /* scaling factor for energy representation */ +#define E_SCALE 100.0 /* scaling factor for energy representation */ #define LOG_SCALE 0.75 /* scaling factor for energy log representation */ #define LOG_SHIFT 0.5 /* shift of colors on log scale */ #define LOG_ENERGY_FLOOR -10.0 /* floor value for log of (total) energy */ @@ -290,16 +263,15 @@ #define MAZE_WIDTH 0.02 /* half width of maze walls */ #define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */ -#define COLORBAR_RANGE 1.0 /* scale of color scheme bar */ -#define COLORBAR_RANGE_B 1.5 /* scale of color scheme bar for 2nd part */ +#define COLORBAR_RANGE 3.0 /* scale of color scheme bar */ +#define COLORBAR_RANGE_B 2.0 /* scale of color scheme bar for 2nd part */ #define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ #define SAVE_TIME_SERIES 0 /* set to 1 to save wave time series at a point */ #define ADD_POTENTIAL 1 /* set to 1 to add potential to z coordinate */ -// #define POT_MAZE 7 #define POTENTIAL 10 -#define POT_FACT 20.0 +#define POT_FACT 200.0 #define DRAW_WAVE_PROFILE 0 /* set to 1 to draw a profile of the wave */ #define MU_B 1.0 /* parameter controlling the dimensions of domain */ #define HORIZONTAL_WAVE_PROFILE 0 /* set to 1 to draw wave profile vertically */ @@ -311,6 +283,10 @@ #define AVERAGE_WAVE_PROFILE 1 /* set to 1 to draw time-average of wave profile squared*/ #define TIMESERIES_NVALUES 400 /* number of values plotted in time series */ #define DRAW_WAVE_SOURCE 0 /* set to 1 to draw source of wave at (wave_source_x, wave_source_y) */ +#define HRES 1 /* dummy, only used by rde.c */ +#define SHADE_2D 0 /* set to 1 to add pseudo-3d shading effect */ +#define SHADE_SCALE_2D 0.05 /* lower value increases sensitivity of shading */ +#define N_SOURCES 1 /* number of sources, for option draw_sources */ /* end of constants only used by sub_wave and sub_maze */ @@ -330,8 +306,8 @@ int reset_view = 0; /* switch to reset 3D view parameters (for option RO #define Z_SCALING_FACTOR 0.15 /* overall scaling factor of z axis for REP_PROJ_3D representation */ #define XY_SCALING_FACTOR 1.7 /* overall scaling factor for on-screen (x,y) coordinates after projection */ #define ZMAX_FACTOR 1.0 /* max value of z coordinate for REP_PROJ_3D representation */ -#define XSHIFT_3D 0.1 /* overall x shift for REP_PROJ_3D representation */ -#define YSHIFT_3D 0.1 /* overall y shift for REP_PROJ_3D representation */ +#define XSHIFT_3D -0.1 /* overall x shift for REP_PROJ_3D representation */ +#define YSHIFT_3D 0.3 /* overall y shift for REP_PROJ_3D representation */ #include "global_pdes.c" /* constants and global variables */ @@ -1174,73 +1150,8 @@ void animation() if ((ADD_OSCILLATING_SOURCE)&&(i%OSCILLATING_SOURCE_PERIOD == 1)) { if (ALTERNATE_OSCILLATING_SOURCE) sign = -sign; - add_circular_wave_mod(sign, -0.5, 0.0, phi, psi, xy_in); - -// phase_shift = 1/30.0; -// // phase_shift = 1/60.0; -// -// if (first_source) for (k=0; k<25; k++) -// { -// angle = 0.0; -// omega = DPI/25.0; -// wave_source[k].xc = 0.05*cos((double)k*omega);; -// wave_source[k].yc = 0.05*sin((double)k*omega); -// wave_source[k].phase = 0.99 - 1.4*sin(0.7*(1.0 + wave_source[k].xc/0.05)); -// wave_source[k].amp = 1.0; -// // if (wave_source[k].phase > 0.0) wave_source[k].sign = 1; -// // else wave_source[k].sign = -1; -// wave_source[k].sign = 1; -// } -// first_source = 0; -// -// for (k=0; k<25; k++) -// wave_source[k].phase += 1.4*sin(0.7*(1.0 + wave_source[k].xc*cos(angle)/0.05 + wave_source[k].yc*sin(angle)/0.05)); -// -// angle = DPI*(double)i/(double)NSTEPS; -// // phase_shift = 0.02 + 0.08*(double)i/(double)NSTEPS; -// // printf("Phase shift = %.3lg\n", phase_shift); -// -// for (k=0; k<25; k++) -// { -// // wave_source[k].phase += 0.07; -// wave_source[k].phase += phase_shift; -// wave_source[k].phase -= 1.4*sin(0.7*(1.0 + wave_source[k].xc*cos(angle)/0.05 + wave_source[k].yc*sin(angle)/0.05)); -// -// if (wave_source[k].phase > 1.0) -// { -// add_circular_wave_mod((double)wave_source[k].sign*wave_source[k].amp, wave_source[k].xc, wave_source[k].yc, phi, psi, xy_in); -// printf("Adding wave at (%.2lg, %.2lg)\n", wave_source[k].xc, wave_source[k].yc); -// wave_source[k].phase -= 1.0; -// wave_source[k].sign *= -1; -// } -// } - -// for (j=0; j 0.083333333*XMIN)&&(x1 < 0.083333333*XMAX)) -// { -// add_circular_wave_mod(sign1, x1, y, phi, psi, xy_in); -// printf("Adding wave at (%.2lg, %.2lg)\n", x1, y); -// } -// sign1 = -sign1; -// } -// source_counter++; -// if (p > 0) q = p; -// else q = -p; -// if (source_counter >= q) -// { -// source_counter = 0; -// sign = -sign; -// } + add_circular_wave_mod(sign, -1.75, 0.3, phi, psi, xy_in); + add_circular_wave_mod(sign, -1.75, -0.3, phi, psi, xy_in); } if (PRINT_SPEED) print_speed_3d(speed, 0, 1.0); diff --git a/wave_billiard.c b/wave_billiard.c index d8f321b..475ee82 100644 --- a/wave_billiard.c +++ b/wave_billiard.c @@ -49,7 +49,7 @@ #define NO_EXTRA_BUFFER_SWAP 1 /* some OS require one less buffer swap when recording images */ #define VARIABLE_IOR 0 /* set to 1 for a variable index of refraction */ -#define IOR 17 /* choice of index of refraction, see list in global_pdes.c */ +#define IOR 191 /* 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 */ @@ -57,23 +57,26 @@ #define WINWIDTH 1920 /* window width */ #define WINHEIGHT 1150 /* window height */ +// #define NX 1920 /* number of grid points on x axis */ +// #define NY 1150 /* number of grid points on y axis */ #define NX 3840 /* number of grid points on x axis */ #define NY 2300 /* number of grid points on y axis */ -#define XMIN -1.6 -#define XMAX 2.4 /* x interval */ +#define XMIN -2.0 +#define XMAX 2.0 /* x interval */ #define YMIN -1.197916667 #define YMAX 1.197916667 /* y interval for 9/16 aspect ratio */ #define HIGHRES 1 /* set to 1 if resolution of grid is double that of displayed image */ +#define HRES 1 /* dummy, only used by rde.c */ #define JULIA_SCALE 1.0 /* scaling for Julia sets */ /* Choice of the billiard table */ -#define B_DOMAIN 75 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN 20 /* choice of domain shape, see list in global_pdes.c */ -#define CIRCLE_PATTERN 103 /* pattern of circles or polygons, see list in global_pdes.c */ +#define CIRCLE_PATTERN 1 /* pattern of circles or polygons, see list in global_pdes.c */ #define COMPARISON 0 /* set to 1 to compare two different patterns (beta) */ #define B_DOMAIN_B 20 /* second domain shape, for comparisons */ @@ -83,9 +86,9 @@ #define NPOISSON 1000 /* number of points for Poisson C_RAND_POISSON arrangement */ #define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */ -#define LAMBDA 3.0 /* parameter controlling the dimensions of domain */ -#define MU 0.14 /* parameter controlling the dimensions of domain */ -#define MU_B 0.42 /* parameter controlling the dimensions of domain */ +#define LAMBDA 0.2 /* parameter controlling the dimensions of domain */ +#define MU 0.012 /* parameter controlling the dimensions of domain */ +#define MU_B 1.0 /* parameter controlling the dimensions of domain */ #define NPOLY 6 /* number of sides of polygon */ #define APOLY -0.666666666666 /* angle by which to turn polygon, in units of Pi/2 */ #define MDEPTH 6 /* depth of computation of Menger gasket */ @@ -93,8 +96,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 60 /* number of grid point for grid of disks */ -#define NGRIDY 10 /* number of grid point for grid of disks */ +#define NGRIDX 24 /* number of grid point for grid of disks */ +#define NGRIDY 40 /* number of grid point for grid of disks */ #define WALL_WIDTH 0.1 /* width of wall separating lenses */ #define X_SHOOTER -0.2 @@ -114,17 +117,19 @@ /* Physical parameters of wave equation */ #define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */ -#define OSCILLATE_LEFT 1 /* set to 1 to add oscilating boundary condition on the left */ +#define OSCILLATE_LEFT 0 /* set to 1 to add oscilating boundary condition on the left */ #define OSCILLATE_TOPBOT 0 /* set to 1 to enforce a planar wave on top and bottom boundary */ -#define OSCILLATION_SCHEDULE 0 /* oscillation schedule, see list in global_pdes.c */ +#define OSCILLATION_SCHEDULE 61 /* oscillation schedule, see list in global_pdes.c */ #define OSCIL_YMAX 0.35 /* defines oscillation range */ +#define INITIAL_SHIFT 20.0 /* time shift of initial wave packet (in oscillation periods) */ +#define WAVE_PACKET_SHIFT 200.0 /* time shift between wave packets (in oscillation periods) */ -#define OMEGA 0.015 /* frequency of periodic excitation */ -#define AMPLITUDE 1.0 /* amplitude of periodic excitation */ +#define OMEGA 0.025 /* frequency of periodic excitation */ +#define AMPLITUDE 0.5 /* amplitude of periodic excitation */ #define ACHIRP 0.25 /* acceleration coefficient in chirp */ #define DAMPING 0.0 /* damping of periodic excitation */ -#define COURANT 0.1 /* Courant number */ -#define COURANTB 0.03 /* Courant number in medium B */ +#define COURANT 0.2 /* Courant number */ +#define COURANTB 0.0 /* Courant number in medium B */ #define GAMMA 0.0 /* damping factor in wave equation */ #define GAMMAB 0.0 /* damping factor in wave equation */ #define GAMMA_SIDES 1.0e-4 /* damping factor on boundary */ @@ -132,15 +137,16 @@ #define KAPPA 0.0 /* "elasticity" term enforcing oscillations */ #define KAPPA_SIDES 5.0e-4 /* "elasticity" term on absorbing boundary */ #define KAPPA_TOPBOT 0.0 /* "elasticity" term on absorbing boundary */ -#define OSCIL_LEFT_YSHIFT 0.0 /* y-dependence of left oscillation (for non-horizontal waves) */ +#define OSCIL_LEFT_YSHIFT 40.0 /* y-dependence of left oscillation (for non-horizontal waves) */ /* The Courant number is given by c*DT/DX, where DT is the time step and DX the lattice spacing */ /* The physical damping coefficient is given by GAMMA/(DT)^2 */ /* Increasing COURANT speeds up the simulation, but decreases accuracy */ /* For similar wave forms, COURANT^2*GAMMA should be kept constant */ -#define ADD_OSCILLATING_SOURCE 0 /* set to 1 to add an oscillating wave source */ -#define OSCILLATING_SOURCE_PERIOD 4 /* period of oscillating source */ +#define ADD_OSCILLATING_SOURCE 1 /* set to 1 to add an oscillating wave source */ +#define OSCILLATING_SOURCE_PERIOD 8 /* period of oscillating source */ #define ALTERNATE_OSCILLATING_SOURCE 1 /* set to 1 to alternate sign of oscillating source */ +#define N_SOURCES 2 /* number of sources, for option draw_sources */ #define ADD_WAVE_PACKET_SOURCES 0 /* set to 1 to add several sources emitting wave packets */ #define WAVE_PACKET_SOURCE_TYPE 3 /* type of wave packet sources */ @@ -151,14 +157,14 @@ /* Boundary conditions, see list in global_pdes.c */ -#define B_COND 3 +#define B_COND 2 /* Parameters for length and speed of simulation */ -#define NSTEPS 2200 /* number of frames of movie */ -#define NVID 15 /* number of iterations between images displayed on screen */ +#define NSTEPS 1600 /* number of frames of movie */ +#define NVID 8 /* number of iterations between images displayed on screen */ #define NSEG 1000 /* number of segments of boundary */ -#define INITIAL_TIME 50 /* time after which to start saving frames */ +#define INITIAL_TIME 0 /* time after which to start saving frames */ #define BOUNDARY_WIDTH 2 /* width of billiard boundary */ #define PRINT_SPEED 0 /* print speed of moving source */ #define PRINT_FREQUENCY 0 /* print frequency (for phased array) */ @@ -173,9 +179,9 @@ /* Parameters of initial condition */ -#define INITIAL_AMP 0.75 /* amplitude of initial condition */ -#define INITIAL_VARIANCE 0.01 /* variance of initial condition */ -#define INITIAL_WAVELENGTH 0.025 /* wavelength of initial condition */ +#define INITIAL_AMP 1.5 /* amplitude of initial condition */ +#define INITIAL_VARIANCE 0.00001 /* variance of initial condition */ +#define INITIAL_WAVELENGTH 0.025 /* wavelength of initial condition */ /* Plot type, see list in global_pdes.c */ @@ -185,8 +191,8 @@ /* Color schemes */ -#define COLOR_PALETTE 17 /* Color palette, see list in global_pdes.c */ -#define COLOR_PALETTE_B 13 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE 11 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE_B 14 /* Color palette, see list in global_pdes.c */ #define BLACK 1 /* background */ @@ -197,12 +203,15 @@ #define PHASE_FACTOR 1.0 /* factor in computation of phase in color scheme P_3D_PHASE */ #define PHASE_SHIFT 0.0 /* shift of phase in color scheme P_3D_PHASE */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ -#define E_SCALE 75.0 /* scaling factor for energy representation */ -#define LOG_SCALE 0.75 /* scaling factor for energy log representation */ -#define LOG_SHIFT 1.5 /* shift of colors on log scale */ +#define E_SCALE 55.0 /* scaling factor for energy representation */ +#define LOG_SCALE 1.0 /* scaling factor for energy log representation */ +#define LOG_SHIFT 0.5 /* shift of colors on log scale */ #define FLUX_SCALE 5.0e3 /* scaling factor for energy flux represtnation */ #define AVRG_E_FACTOR 0.95 /* controls time window size in P_AVERAGE_ENERGY scheme */ #define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */ +#define FADE_IN_OBSTACLE 0 /* set to 1 to fade color inside obstacles */ +#define SHADE_2D 1 /* set to 1 to add pseudo-3d shading effect */ +#define SHADE_SCALE_2D 0.05 /* lower value increases sensitivity of shading */ #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 */ @@ -212,23 +221,23 @@ #define HUEAMP -180.0 /* amplitude of variation of hue for color scheme C_HUE */ #define DRAW_COLOR_SCHEME 0 /* set to 1 to plot the color scheme */ -#define COLORBAR_RANGE 3.0 /* scale of color scheme bar */ -#define COLORBAR_RANGE_B 1.5 /* scale of color scheme bar for 2nd part */ -#define ROTATE_COLOR_SCHEME 1 /* set to 1 to draw color scheme horizontally */ +#define COLORBAR_RANGE 2.0 /* scale of color scheme bar */ +#define COLORBAR_RANGE_B 0.8 /* scale of color scheme bar for 2nd part */ +#define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ #define CIRC_COLORBAR 0 /* set to 1 to draw circular color scheme */ #define CIRC_COLORBAR_B 0 /* set to 1 to draw circular color scheme */ #define DRAW_WAVE_PROFILE 1 /* set to 1 to draw a profile of the wave */ #define HORIZONTAL_WAVE_PROFILE 0 /* set to 1 to draw wave profile vertically */ #define VERTICAL_WAVE_PROFILE 1 /* set to 1 to draw wave profile vertically */ -#define WAVE_PROFILE_X 2.1 /* value of x to sample wave profile */ -#define WAVE_PROFILE_Y -1.0 /* value of y to sample wave profile */ +#define WAVE_PROFILE_X 1.9 /* value of x to sample wave profile */ +#define WAVE_PROFILE_Y 0.075 /* value of y to sample wave profile */ #define PROFILE_AT_BOTTOM 1 /* draw wave profile at bottom instead of top */ #define AVERAGE_WAVE_PROFILE 1 /* set to 1 to draw time-average of wave profile squared*/ #define DRAW_WAVE_TIMESERIES 0 /* set to 1 to draw a time series of the wave, 2 to also draw it at the top */ #define TIMESERIES_NVALUES 400 /* number of values plotted in time series */ #define SAVE_TIME_SERIES 0 /* set to 1 to save wave time series at a point */ -#define DRAW_WAVE_SOURCE 0 /* set to 1 to draw source of wave at (wave_source_x, wave_source_y) */ +#define DRAW_WAVE_SOURCE 1 /* set to 1 to draw source of wave at (wave_source_x, wave_source_y), set to 2 to also draw focus */ #define MESSAGE_LDASH 14 /* length of dash for Morse code message */ #define MESSAGE_LDOT 8 /* length of dot for Morse code message */ @@ -257,6 +266,8 @@ #define MEAN_FLUX (PLOT == P_TOTAL_ENERGY_FLUX)||(PLOT_B == P_TOTAL_ENERGY_FLUX) #define REFRESH_IOR ((IOR == IOR_PERIODIC_WELLS_ROTATING)||(IOR == IOR_PERIODIC_WELLS_ROTATING_LARGE)) +double light[2] = {0.40824829, 0.816496581}; /* location of light source for SHADE_2D option*/ + #include "global_pdes.c" /* constants and global variables */ #include "sub_maze.c" /* support for generating mazes */ #include "sub_wave.c" /* common functions for wave_billiard, heat and schrodinger */ @@ -264,8 +275,6 @@ FILE *time_series_left, *time_series_right; -double courant2, courantb2; /* Courant parameters squared */ - /*********************/ /* animation part */ /*********************/ @@ -296,22 +305,11 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX for (j=0; j vmax[pnb]) vmax[pnb] = average_vals[i]; } - else + else /*if (!fade)*/ { for (i=imin; i NX-1) ij[0] = NX-1; + ival = ij[0]; + ival -= (ival%size); + } + + if ((average)&&(!fade)) { for (j=jmin; j vmax[pnb]) vmax[pnb] = average_vals[j]; } - else + + + else if (!fade) { for (j=jmin; j 0.0)) +// { +// if (vmax[pnb] > -vmin[pnb]) vmin[pnb] = -vmax[pnb]; +// else if (vmax[pnb] < -vmin[pnb]) vmax[pnb] = -vmin[pnb]; +// } + deltav = vmax[pnb]-vmin[pnb]; if (deltav <= 0.0) deltav = 0.01; a = deltai/deltav; @@ -1340,8 +1360,11 @@ void draw_wave_profile_vertical(double *values, int size, int average, int pnb, glEnd(); /* slightly decrease extramal values in case signal gets weaker */ - vmax[pnb] *= 0.99; - vmin[pnb] *= 0.99; + if (!fade) + { + vmax[pnb] *= 0.99; + vmin[pnb] *= 0.99; + } } @@ -1692,17 +1715,57 @@ void draw_wave_timeseries_old(double *values, int size, int fade, double fade_va draw_entrance_timeseries(1.0, 1.0, values, size, fade, fade_value); } -void draw_wave_highres_palette(int size, double *phi[NX], double *psi[NX], double *total_energy[NX], double *average_energy[NX], double *total_flux, short int *xy_in[NX], double scale, int time, int plot, int palette, int pnumber, int fade, double fade_value) + +void init_fade_table(double *tcc_table[NX], double *fade_table) +{ + int i, j; + + for (i=0; i= 1) { glColor3f(0.0, 0.0, 0.0); - draw_circle(wave_source_x, wave_source_y, 0.005, 100); + for (k=0; k= 2) + draw_circle(focus_x, 0.0, 0.005, 100); } + free(values); + free(rgbvals); } /* modified function for "flattened" wave tables */ diff --git a/wave_comparison.c b/wave_comparison.c index 1125c23..3f455df 100644 --- a/wave_comparison.c +++ b/wave_comparison.c @@ -270,6 +270,14 @@ #define AVERAGE_WAVE_PROFILE 1 /* set to 1 to draw time-average of wave profile squared*/ #define TIMESERIES_NVALUES 400 /* number of values plotted in time series */ #define DRAW_WAVE_SOURCE 0 /* set to 1 to draw source of wave at (wave_source_x, wave_source_y) */ +#define HRES 1 /* dummy, only used by rde.c */ +#define INITIAL_SHIFT 20.0 /* time shift of initial wave packet (in oscillation periods) */ +#define WAVE_PACKET_SHIFT 200.0 /* time shift between wave packets (in oscillation periods) */ +#define FADE_IN_OBSTACLE 0 /* set to 1 to fade color inside obstacles */ +#define SHADE_2D 0 /* set to 1 to add pseudo-3d shading effect */ +#define SHADE_SCALE_2D 0.05 /* lower value increases sensitivity of shading */ +#define N_SOURCES 1 /* number of sources, for option draw_sources */ +double light[2] = {0.40824829, 0.816496581}; /* location of light source for SHADE_2D option*/ /* end of constants only used by sub_wave and sub_maze */ diff --git a/wave_energy.c b/wave_energy.c index e3c17b2..49322fd 100644 --- a/wave_energy.c +++ b/wave_energy.c @@ -244,6 +244,14 @@ #define AVERAGE_WAVE_PROFILE 1 /* set to 1 to draw time-average of wave profile squared*/ #define TIMESERIES_NVALUES 400 /* number of values plotted in time series */ #define DRAW_WAVE_SOURCE 0 /* set to 1 to draw source of wave at (wave_source_x, wave_source_y) */ +#define HRES 1 /* dummy, only used by rde.c */ +#define INITIAL_SHIFT 20.0 /* time shift of initial wave packet (in oscillation periods) */ +#define WAVE_PACKET_SHIFT 200.0 /* time shift between wave packets (in oscillation periods) */ +#define FADE_IN_OBSTACLE 0 /* set to 1 to fade color inside obstacles */ +#define SHADE_2D 0 /* set to 1 to add pseudo-3d shading effect */ +#define SHADE_SCALE_2D 0.05 /* lower value increases sensitivity of shading */ +#define N_SOURCES 1 /* number of sources, for option draw_sources */ +double light[2] = {0.40824829, 0.816496581}; /* location of light source for SHADE_2D option*/ /* end of constants only used by sub_wave and sub_maze */ #include "global_pdes.c" /* constants and global variables */ diff --git a/wave_sphere.c b/wave_sphere.c index 1a706d5..6434467 100644 --- a/wave_sphere.c +++ b/wave_sphere.c @@ -320,6 +320,11 @@ #define ADD_POTENTIAL 0 /* set to 1 to add potential to z coordinate */ #define POTENTIAL 10 #define POT_FACT 20.0 +#define HRES 1 /* dummy, only used by rde.c */ +#define INITIAL_SHIFT 20.0 /* time shift of initial wave packet (in oscillation periods) */ +#define WAVE_PACKET_SHIFT 200.0 /* time shift between wave packets (in oscillation periods) */ +#define FADE_IN_OBSTACLE 0 /* set to 1 to fade color inside obstacles */ +#define N_SOURCES 1 /* number of sources, for option draw_sources */ /* end of constants only used by sub_wave and sub_maze */ /* For debugging purposes only */