/*********************/ /* Graphics routines */ /*********************/ #include "colors_waves.c" // #define HUE_TYPE0 260.0 /* hue of particles of type 0 */ // #define HUE_TYPE0 300.0 /* hue of particles of type 0 */ // #define HUE_TYPE1 90.0 /* hue of particles of type 1 */ int writetiff(char *filename, char *description, int x, int y, int width, int height, int compression) { TIFF *file; GLubyte *image, *p; int i; file = TIFFOpen(filename, "w"); if (file == NULL) { return 1; } image = (GLubyte *) malloc(width * height * sizeof(GLubyte) * 3); /* OpenGL's default 4 byte pack alignment would leave extra bytes at the end of each image row so that each full row contained a number of bytes divisible by 4. Ie, an RGB row with 3 pixels and 8-bit componets would be laid out like "RGBRGBRGBxxx" where the last three "xxx" bytes exist just to pad the row out to 12 bytes (12 is divisible by 4). To make sure the rows are packed as tight as possible (no row padding), set the pack alignment to 1. */ glPixelStorei(GL_PACK_ALIGNMENT, 1); glReadPixels(x, y, width, height, GL_RGB, GL_UNSIGNED_BYTE, image); TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32_t) width); TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32_t) height); TIFFSetField(file, TIFFTAG_BITSPERSAMPLE, 8); TIFFSetField(file, TIFFTAG_COMPRESSION, compression); TIFFSetField(file, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB); TIFFSetField(file, TIFFTAG_ORIENTATION, ORIENTATION_BOTLEFT); TIFFSetField(file, TIFFTAG_SAMPLESPERPIXEL, 3); TIFFSetField(file, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG); TIFFSetField(file, TIFFTAG_ROWSPERSTRIP, 1); TIFFSetField(file, TIFFTAG_IMAGEDESCRIPTION, description); p = image; for (i = height - 1; i >= 0; i--) { // if (TIFFWriteScanline(file, p, height - i - 1, 0) < 0) if (TIFFWriteScanline(file, p, i, 0) < 0) { free(image); TIFFClose(file); return 1; } p += width * sizeof(GLubyte) * 3; } TIFFClose(file); if (SAVE_MEMORY) free(image); return 0; } void init() /* initialisation of window */ { glLineWidth(3); glClearColor(0.0, 0.0, 0.0, 1.0); glClear(GL_COLOR_BUFFER_BIT); glOrtho(XMIN, XMAX, YMIN, YMAX , -1.0, 1.0); // glOrtho(0.0, NX, 0.0, NY, -1.0, 1.0); } void blank() { if (BLACK) glClearColor(0.0, 0.0, 0.0, 1.0); else glClearColor(1.0, 1.0, 1.0, 1.0); glClear(GL_COLOR_BUFFER_BIT); } void save_frame_lj() { static int counter = 0; char *name="lj.", n2[100]; char format[6]=".%05i"; counter++; // printf (" p2 counter = %d \n",counter); strcpy(n2, name); sprintf(strstr(n2,"."), format, counter); strcat(n2, ".tif"); printf(" saving frame %s \n",n2); writetiff(n2, "Particles with Lennard-Jones interaction in a planar domain", 0, 0, WINWIDTH, WINHEIGHT, COMPRESSION_LZW); } void save_frame_lj_counter(int counter) { char *name="lj.", n2[100]; char format[6]=".%05i"; strcpy(n2, name); sprintf(strstr(n2,"."), format, counter); strcat(n2, ".tif"); printf(" saving frame %s \n",n2); writetiff(n2, "Particles with Lennard-Jones interaction in a planar domain", 0, 0, WINWIDTH, WINHEIGHT, COMPRESSION_LZW); } void write_text_fixedwidth( double x, double y, char *st) { int l, i; l=strlen( st ); // see how many characters are in text string. glRasterPos2d( x, y); // location to start printing text for( i=0; i < l; i++) // loop until i is greater then l { // glutBitmapCharacter(GLUT_BITMAP_TIMES_ROMAN_24, st[i]); // Print a character on the screen // glutBitmapCharacter(GLUT_BITMAP_8_BY_13, st[i]); // Print a character on the screen glutBitmapCharacter(GLUT_BITMAP_9_BY_15, st[i]); // Print a character on the screen } } void write_text( double x, double y, char *st) { int l,i; l=strlen( st ); // see how many characters are in text string. glRasterPos2d( x, y); // location to start printing text for( i=0; i < l; i++) // loop until i is greater then l { glutBitmapCharacter(GLUT_BITMAP_TIMES_ROMAN_24, st[i]); // Print a character on the screen // glutBitmapCharacter(GLUT_BITMAP_8_BY_13, st[i]); // Print a character on the screen } } /*********************/ /* some basic math */ /*********************/ double vabs(double x) /* absolute value */ { double res; if (x<0.0) res = -x; else res = x; return(res); } double argument(double x, double y) { double alph; if (x!=0.0) { alph = atan(y/x); if (x<0.0) alph += PI; } else { alph = PID; if (y<0.0) alph = PI*1.5; } return(alph); } double ipow(double x, int n) { double y; int i; y = x; for (i=1; i= 1 || S == 0); X = V1 * sqrt(-2 * log(S) / S); } else X = V2 * sqrt(-2 * log(S) / S); phase = 1 - phase; return X; } /*********************/ /* drawing routines */ /*********************/ void erase_area(double x, double y, double dx, double dy) { double pos[2], rgb[3]; hsl_to_rgb(220.0, 0.8, 0.7, rgb); glColor3f(rgb[0], rgb[1], rgb[2]); glBegin(GL_QUADS); glVertex2d(x - dx, y - dy); glVertex2d(x + dx, y - dy); glVertex2d(x + dx, y + dy); glVertex2d(x - dx, y + dy); glEnd(); } void erase_area_rgb(double x, double y, double dx, double dy, double rgb[3]) { double pos[2]; glColor3f(rgb[0], rgb[1], rgb[2]); glBegin(GL_QUADS); glVertex2d(x - dx, y - dy); glVertex2d(x + dx, y - dy); glVertex2d(x + dx, y + dy); glVertex2d(x - dx, y + dy); glEnd(); } void erase_rect_rgb(double xmin, double ymin, double xmax, double ymax, double rgb[3]) { double pos[2]; glColor3f(rgb[0], rgb[1], rgb[2]); glBegin(GL_QUADS); glVertex2d(xmin, ymin); glVertex2d(xmax, ymin); glVertex2d(xmax, ymax); glVertex2d(xmin, ymax); glEnd(); } void erase_area_hsl(double x, double y, double dx, double dy, double h, double s, double l) { double pos[2], rgb[3]; hsl_to_rgb(h, s, l, rgb); erase_area_rgb(x, y, dx, dy, rgb); } void erase_area_hsl_turbo(double x, double y, double dx, double dy, double h, double s, double l) { double pos[2], rgb[3]; hsl_to_rgb_turbo(h, s, l, rgb); erase_area_rgb(x, y, dx, dy, rgb); } void erase_rectangle_hsl_turbo(double xmin, double ymin, double xmax, double ymax, double h, double s, double l) { double pos[2], rgb[3]; hsl_to_rgb_turbo(h, s, l, rgb); erase_rect_rgb(xmin, ymin, xmax, ymax, rgb); } void draw_line(double x1, double y1, double x2, double y2) { glBegin(GL_LINE_STRIP); glVertex2d(x1, y1); glVertex2d(x2, y2); glEnd(); } void draw_arrow(double x1, double y1, double x2, double y2, double angle, double length) { double alpha, beta, x3, y3, x4, y4, x5, y5; alpha = argument(x2 - x1, y2 - y1); beta = angle*PI/180.0; x3 = x2 - length*cos(alpha - beta); y3 = y2 - length*sin(alpha - beta); x4 = x2 - length*cos(alpha + beta); y4 = y2 - length*sin(alpha + beta); x5 = x2 - 0.5*length*cos(alpha); y5 = y2 - 0.5*length*sin(alpha); glBegin(GL_LINE_STRIP); glVertex2d(x1, y1); glVertex2d(x5, y5); glEnd(); glBegin(GL_TRIANGLE_FAN); glVertex2d(x2, y2); glVertex2d(x3, y3); glVertex2d(x4, y4); glEnd(); } void draw_rectangle(double x1, double y1, double x2, double y2) { glBegin(GL_LINE_LOOP); glVertex2d(x1, y1); glVertex2d(x2, y1); glVertex2d(x2, y2); glVertex2d(x1, y2); glEnd(); } void draw_colored_rectangle(double x1, double y1, double x2, double y2, double rgb[3]) { glColor3f(rgb[0], rgb[1], rgb[2]); glBegin(GL_TRIANGLE_FAN); glVertex2d(x1, y1); glVertex2d(x2, y1); glVertex2d(x2, y2); glVertex2d(x1, y2); glEnd(); } void draw_triangle(double x1, double y1, double x2, double y2, double x3, double y3) { glBegin(GL_LINE_LOOP); glVertex2d(x1, y1); glVertex2d(x2, y2); glVertex2d(x3, y3); glEnd(); } void draw_colored_triangle(double x1, double y1, double x2, double y2, double x3, double y3, double rgb[3]) { glColor3f(rgb[0], rgb[1], rgb[2]); glBegin(GL_TRIANGLE_FAN); glVertex2d(x1, y1); glVertex2d(x2, y2); glVertex2d(x3, y3); glEnd(); } double triangle_area(double x1, double y1, double x2, double y2, double x3, double y3) { return((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)); } void draw_colored_triangle_turbo(double x1, double y1, double x2, double y2, double x3, double y3, double h, double s, double l) { double rgb[3]; hsl_to_rgb_turbo(h, s, l, rgb); draw_colored_triangle(x1, y1, x2, y2, x3, y3, rgb); } void draw_circle(double x, double y, double r, int nseg) { int i; double pos[2], alpha, dalpha; dalpha = DPI/(double)nseg; glBegin(GL_LINE_LOOP); for (i=0; i<=nseg; i++) { alpha = (double)i*dalpha; glVertex2d(x + r*cos(alpha), y + r*sin(alpha)); } glEnd(); } void draw_colored_circle(double x, double y, double r, int nseg, double rgb[3]) { int i; double pos[2], alpha, dalpha; dalpha = DPI/(double)nseg; glColor3f(rgb[0], rgb[1], rgb[2]); glBegin(GL_TRIANGLE_FAN); glVertex2d(x, y); for (i=0; i<=nseg; i++) { alpha = (double)i*dalpha; glVertex2d(x + r*cos(alpha), y + r*sin(alpha)); } glEnd(); } void draw_circle_precomp(double x, double y, double r) { int i; glBegin(GL_LINE_LOOP); for (i=0; i<=NSEG; i++) glVertex2d(x + r*cosangle[i], y + r*sinangle[i]); glEnd(); } void draw_colored_circle_precomp(double x, double y, double r, double rgb[3]) { int i; glColor3f(rgb[0], rgb[1], rgb[2]); glBegin(GL_TRIANGLE_FAN); glVertex2d(x, y); for (i=0; i<=NSEG; i++) glVertex2d(x + r*cosangle[i], y + r*sinangle[i]); glEnd(); } void draw_polygon(double x, double y, double r, int nsides, double angle) { int i; double pos[2], alpha, dalpha; dalpha = DPI/(double)nsides; glColor3f(1.0, 1.0, 1.0); glBegin(GL_LINE_LOOP); for (i=0; i<=nsides; i++) { alpha = angle + (double)i*dalpha; glVertex2d(x + r*cos(alpha), y + r*sin(alpha)); } glEnd(); } void draw_colored_polygon(double x, double y, double r, int nsides, double angle, double rgb[3]) { int i; double pos[2], alpha, dalpha; dalpha = DPI/(double)nsides; glColor3f(rgb[0], rgb[1], rgb[2]); glBegin(GL_TRIANGLE_FAN); glVertex2d(x, y); for (i=0; i<=nsides; i++) { alpha = angle + (double)i*dalpha; glVertex2d(x + r*cos(alpha), y + r*sin(alpha)); } glEnd(); } void draw_rhombus(double x, double y, double r, double angle) { int i; static int first = 1; static double ratio; if (first) { ratio = tan(0.1*PI); first = 0; } glBegin(GL_LINE_LOOP); glVertex2d(x + r*cos(angle), y + r*sin(angle)); glVertex2d(x - ratio*r*sin(angle), y + ratio*r*cos(angle)); glVertex2d(x - r*cos(angle), y - r*sin(angle)); glVertex2d(x + ratio*r*sin(angle), y - ratio*r*cos(angle)); glEnd(); } void draw_colored_rhombus(double x, double y, double r, double angle, double rgb[3]) { int i; static int first = 1; static double ratio; if (first) { ratio = tan(0.1*PI); first = 0; } glColor3f(rgb[0], rgb[1], rgb[2]); glBegin(GL_TRIANGLE_FAN); glVertex2d(x, y); glVertex2d(x + r*cos(angle), y + r*sin(angle)); glVertex2d(x - ratio*r*sin(angle), y + ratio*r*cos(angle)); glVertex2d(x - r*cos(angle), y - r*sin(angle)); glVertex2d(x + ratio*r*sin(angle), y - ratio*r*cos(angle)); glVertex2d(x + r*cos(angle), y + r*sin(angle)); glEnd(); } void draw_rotated_rect(double x, double y, double l, double w, double angle) { int i; double ca, sa, lca, lsa, wca, wsa; ca = cos(angle); sa = sin(angle); lca = l*ca; lsa = l*sa; wca = w*ca; wsa = w*sa; glBegin(GL_LINE_LOOP); glVertex2d(x + lca - wsa, y + lsa + wca); glVertex2d(x - lca - wsa, y - lsa + wca); glVertex2d(x - lca + wsa, y - lsa - wca); glVertex2d(x + lca + wsa, y + lsa - wca); glEnd(); } void draw_colored_rotated_rect(double x, double y, double l, double w, double angle, double rgb[3]) { int i; double ca, sa, lca, lsa, wca, wsa; ca = cos(angle); sa = sin(angle); lca = l*ca; lsa = l*sa; wca = w*ca; wsa = w*sa; glColor3f(rgb[0], rgb[1], rgb[2]); glBegin(GL_TRIANGLE_FAN); glVertex2d(x, y); glVertex2d(x + lca - wsa, y + lsa + wca); glVertex2d(x - lca - wsa, y - lsa + wca); glVertex2d(x - lca + wsa, y - lsa - wca); glVertex2d(x + lca + wsa, y + lsa - wca); glVertex2d(x + lca - wsa, y + lsa + wca); glEnd(); } void draw_colored_sector(double xc, double yc, double r1, double r2, double angle1, double angle2, double rgb[3], int nsides) { int i; double angle, dangle; dangle = (angle2 - angle1)/(double)nsides; glColor3f(rgb[0], rgb[1], rgb[2]); glBegin(GL_TRIANGLE_FAN); glVertex2d(xc + r1*cos(angle1), yc + r1*sin(angle1)); for (i = 0; i < nsides+1; i++) { angle = angle1 + dangle*(double)i; glVertex2d(xc + r2*cos(angle), yc + r2*sin(angle)); } glEnd(); glBegin(GL_TRIANGLE_FAN); glVertex2d(xc + r2*cos(angle2), yc + r2*sin(angle2)); for (i = 0; i < nsides+1; i++) { angle = angle1 + dangle*(double)i; glVertex2d(xc + r1*cos(angle), yc + r1*sin(angle)); } glEnd(); } double neighbour_color(int nnbg) { if (nnbg > 7) nnbg = 7; switch(nnbg){ case (7): return(340.0); case (6): return(310.0); case (5): return(260.0); case (4): return(200.0); case (3): return(140.0); case (2): return(100.0); case (1): return(70.0); default: return(30.0); } } double partner_color(int np) { switch(np){ case (0): return(340.0); case (1): return(260.0); case (2): return(210.0); case (3): return(140.0); case (4): return(70.0); default: return(20.0); // case (0): return(70.0); // case (1): return(200.0); // case (2): return(280.0); // case (3): return(140.0); // case (4): return(320.0); // default: return(20.0); } } double type_hue(int type) { int hue; double t2; static double b, hmax; static int first = 1; if (first) { hmax = 360.0; b = 16.0*(hmax - HUE_TYPE3); first = 0; } if ((RD_REACTION == CHEM_CATALYTIC_A2D)&&(type == 4)) return(HUE_TYPE3); if ((RD_REACTION == CHEM_ABDACBE)&&(type == 4)) return(HUE_TYPE3); if ((RD_REACTION == CHEM_ABDACBE)&&(type == 5)) return(280.0); switch (type) { case (0): return(HUE_TYPE0); case (1): return(HUE_TYPE1); case (2): return(HUE_TYPE2); case (3): return(HUE_TYPE3); case (4): return(HUE_TYPE4); case (5): return(HUE_TYPE5); case (6): return(HUE_TYPE6); case (7): { if (RD_REACTION == CHEM_BZ) return(HUE_TYPE2); else return(HUE_TYPE7); } default: { if (RD_REACTION == CHEM_BZ) { if (type == 7) return(HUE_TYPE2); if (type == 8) type = 5; } else if ((RD_REACTION == CHEM_BRUSSELATOR)&&(type >= 5)) return(70.0); t2 = (double)(type*type); hue = (hmax*t2 - b)/t2; return(hue); } } } void set_type_color(int type, double lum, double rgb[3]) { int hue; if (COUNT_PARTNER_TYPE) hue = partner_color(type-1); else hue = type_hue(type); hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE); glColor3f(lum*rgb[0], lum*rgb[1], lum*rgb[2]); } double distance_to_segment(double x, double y, double x1, double y1, double x2, double y2) /* distance of (x,y) to segment from (x1,y1) to (x2,y2) */ { double xp, yp, angle, length, ca, sa; angle = argument(x2 - x1, y2 - y1); length = module2(x2 - x1, y2 - y1); ca = cos(angle); sa = sin(angle); xp = ca*(x - x1) + sa*(y - y1); yp = -sa*(x - x1) + ca*(y - y1); if ((xp >= 0)&&(xp <= length)) return(vabs(yp)); else if (xp < 0) return(module2(xp, yp)); else return(module2(xp-length, yp)); } int in_polygon(double x, double y, double r, int npoly, double apoly) /* test whether (x,y) is in regular polygon of npoly sides inscribed in circle of radius r, turned by apoly Pi/2 */ { int condition = 1, k; double omega, cw, angle; omega = DPI/((double)npoly); cw = cos(omega*0.5); for (k=0; k 0)&&(npoints < nmax)) { /* randomly select an active circle */ i = rand()%(npoints); while (!active_poisson[i]) i = rand()%(npoints); /* generate new candidates */ naccepted = 0; for (j=0; j= dpoisson*dpoisson); /* new circle is in domain */ far = far*(x < xmax)*(x > xmin)*(y < ymax)*(y > ymin); } if (far) /* accept new circle */ { printf("New npoint at (%.3f,%.3f) accepted\n", x, y); point[npoints].xc = x; point[npoints].yc = y; active_poisson[npoints] = 1; npoints++; n_p_active++; naccepted++; } } if (naccepted == 0) /* inactivate circle i */ { active_poisson[i] = 0; n_p_active--; } printf("%i active npoints\n", n_p_active); } printf("Generated %i points\n", npoints); return(npoints); } void init_particle_config(t_particle particles[NMAXCIRCLES]) /* initialise particle configuration */ { int i, j, k, n, ncirc0, n_p_active, ncandidates = PDISC_CANDIDATES, naccepted; double dx, dy, p, phi, r, r0, ra[5], sa[5], height, x, y = 0.0, gamma, dpoisson = PDISC_DISTANCE*MU, xx[4], yy[4]; short int active_poisson[NMAXCIRCLES], far; switch (CIRCLE_PATTERN) { case (C_SQUARE): { ncircles = NGRIDX*NGRIDY; dx = (INITXMAX - INITXMIN)/((double)NGRIDX); dy = (INITYMAX - INITYMIN)/((double)NGRIDY); for (i = 0; i < NGRIDX; i++) for (j = 0; j < NGRIDY; j++) { n = NGRIDY*i + j; particles[n].xc = INITXMIN + ((double)i - 0.5)*dx; particles[n].yc = INITYMIN + ((double)j - 0.5)*dy; particles[n].radius = MU; /* activate only circles that intersect the domain */ if ((particles[n].yc < INITYMAX + MU)&&(particles[n].yc > INITYMIN - MU)&&(particles[n].xc < INITXMAX + MU)&&(particles[n].xc > INITXMIN - MU)) particles[n].active = 1; else particles[n].active = 0; } break; } case (C_HEX): { ncircles = NGRIDX*(NGRIDY+1); dx = (INITXMAX - INITXMIN)/((double)NGRIDX); dy = (INITYMAX - INITYMIN)/((double)NGRIDY); // dx = dy*0.5*sqrt(3.0); for (i = 0; i < NGRIDX; i++) for (j = 0; j < NGRIDY+1; j++) { n = (NGRIDY+1)*i + j; // particles[n].xc = ((double)(i-NGRIDX/2) + 0.5)*dx; /* is +0.5 needed? */ particles[n].xc = INITXMIN + ((double)i - 0.5)*dx; particles[n].yc = INITYMIN + ((double)j - 0.5)*dy; if ((i+NGRIDX)%2 == 1) particles[n].yc += 0.5*dy; 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 */ if ((particles[n].yc < INITYMAX + MU)&&(particles[n].yc > INITYMIN - MU)&&(particles[n].xc < INITXMAX + MU)&&(particles[n].xc > INITXMIN - MU)) particles[n].active = 1; else particles[n].active = 0; } break; } case (C_RAND_DISPLACED): { ncircles = NGRIDX*NGRIDY; dy = (YMAX - YMIN)/((double)NGRIDY); for (i = 0; i < NGRIDX; i++) for (j = 0; j < NGRIDY; j++) { n = NGRIDY*i + j; particles[n].xc = ((double)(i-NGRIDX/2) + 0.5*((double)rand()/RAND_MAX - 0.5))*dy; particles[n].yc = YMIN + ((double)j + 0.5 + 0.5*((double)rand()/RAND_MAX - 0.5))*dy; particles[n].radius = MU; // particles[n].radius = MU*sqrt(1.0 + 0.8*((double)rand()/RAND_MAX - 0.5)); particles[n].active = 1; } break; } case (C_RAND_PERCOL): { ncircles = NGRIDX*NGRIDY; dy = (YMAX - YMIN)/((double)NGRIDY); for (i = 0; i < NGRIDX; i++) for (j = 0; j < NGRIDY; j++) { n = NGRIDY*i + j; particles[n].xc = ((double)(i-NGRIDX/2) + 0.5)*dy; particles[n].yc = YMIN + ((double)j + 0.5)*dy; particles[n].radius = MU; p = (double)rand()/RAND_MAX; if (p < P_PERCOL) particles[n].active = 1; else particles[n].active = 0; } break; } case (C_RAND_POISSON): { ncircles = NPOISSON; for (n = 0; n < NPOISSON; n++) { // particles[n].xc = LAMBDA*(2.0*(double)rand()/RAND_MAX - 1.0); particles[n].xc = (XMAX - XMIN)*(double)rand()/RAND_MAX + XMIN; particles[n].yc = (YMAX - YMIN)*(double)rand()/RAND_MAX + YMIN; particles[n].radius = MU; particles[n].active = 1; } break; } case (C_CLOAK): { ncircles = 200; for (i = 0; i < 40; i++) for (j = 0; j < 5; j++) { n = 5*i + j; phi = (double)i*DPI/40.0; r = LAMBDA*0.5*(1.0 + (double)j/5.0); particles[n].xc = r*cos(phi); particles[n].yc = r*sin(phi); particles[n].radius = MU; particles[n].active = 1; } break; } case (C_CLOAK_A): /* optimized model A1 by C. Jo et al */ { ncircles = 200; ra[0] = 0.0731; sa[0] = 1.115; ra[1] = 0.0768; sa[1] = 1.292; ra[2] = 0.0652; sa[2] = 1.464; ra[3] = 0.056; sa[3] = 1.633; ra[4] = 0.0375; sa[4] = 1.794; for (i = 0; i < 40; i++) for (j = 0; j < 5; j++) { n = 5*i + j; phi = (double)i*DPI/40.0; r = LAMBDA*sa[j]; particles[n].xc = r*cos(phi); particles[n].yc = r*sin(phi); particles[n].radius = LAMBDA*ra[j]; particles[n].active = 1; } break; } case (C_LASER): { ncircles = 17; xx[0] = 0.5*(X_SHOOTER + X_TARGET); xx[1] = LAMBDA - 0.5*(X_TARGET - X_SHOOTER); xx[2] = -xx[0]; xx[3] = -xx[1]; yy[0] = 0.5*(Y_SHOOTER + Y_TARGET); yy[1] = 1.0 - 0.5*(Y_TARGET - Y_SHOOTER); yy[2] = -yy[0]; yy[3] = -yy[1]; for (i = 0; i < 4; i++) for (j = 0; j < 4; j++) { particles[4*i + j].xc = xx[i]; particles[4*i + j].yc = yy[j]; } particles[ncircles - 1].xc = X_TARGET; particles[ncircles - 1].yc = Y_TARGET; for (i=0; i 0)&&(ncircles < NMAXCIRCLES)) { /* randomly select an active circle */ i = rand()%(ncircles); while (!active_poisson[i]) i = rand()%(ncircles); // printf("Starting from circle %i at (%.3f,%.3f)\n", i, particles[i].xc, particles[i].yc); /* generate new candidates */ naccepted = 0; for (j=0; j= dpoisson*dpoisson); /* new circle is in domain */ far = far*(x < INITXMAX)*(x > INITXMIN)*(y < INITYMAX)*(y > INITYMIN); // far = far*(vabs(x) < LAMBDA)*(y < INITYMAX)*(y > INITYMIN); } if (far) /* accept new circle */ { printf("New circle at (%.3f,%.3f) accepted\n", x, y); particles[ncircles].xc = x; particles[ncircles].yc = y; particles[ncircles].radius = MU; particles[ncircles].active = 1; active_poisson[ncircles] = 1; ncircles++; n_p_active++; naccepted++; } // else printf("Rejected\n"); } if (naccepted == 0) /* inactivate circle i */ { // printf("No candidates work, inactivate circle %i\n", i); active_poisson[i] = 0; n_p_active--; } printf("%i active circles\n", n_p_active); } printf("Generated %i circles\n", ncircles); break; } case (C_GOLDEN_MEAN): { ncircles = 300; gamma = (sqrt(5.0) - 1.0)*0.5; /* golden mean */ height = YMAX - YMIN; dx = 2.0*LAMBDA/((double)ncircles); for (n = 0; n < ncircles; n++) { particles[n].xc = -LAMBDA + n*dx; particles[n].yc = y; y += height*gamma; if (y > YMAX) y -= height; particles[n].radius = MU; particles[n].active = 1; } /* test for circles that overlap top or bottom boundary */ ncirc0 = ncircles; for (n=0; n < ncirc0; n++) { if (particles[n].yc + particles[n].radius > YMAX) { particles[ncircles].xc = particles[n].xc; particles[ncircles].yc = particles[n].yc - height; particles[ncircles].radius = MU; particles[ncircles].active = 1; ncircles ++; } else if (particles[n].yc - particles[n].radius < YMIN) { particles[ncircles].xc = particles[n].xc; particles[ncircles].yc = particles[n].yc + height; particles[ncircles].radius = MU; particles[ncircles].active = 1; ncircles ++; } } break; } case (C_GOLDEN_SPIRAL): { ncircles = 1; particles[0].xc = 0.0; particles[0].yc = 0.0; gamma = (sqrt(5.0) - 1.0)*PI; /* golden mean times 2Pi */ phi = 0.0; r0 = 2.0*MU; r = r0 + MU; for (i=0; i<1000; i++) { x = r*cos(phi); y = r*sin(phi); phi += gamma; r += MU*r0/r; if ((vabs(x) < LAMBDA)&&(vabs(y) < YMAX + MU)) { particles[ncircles].xc = x; particles[ncircles].yc = y; ncircles++; } } for (i=0; i YMIN - MU)) particles[i].active = 1; // printf("i = %i, circlex = %.3lg, circley = %.3lg\n", i, particles[i].xc, particles[i].yc); } break; } case (C_SQUARE_HEX): { ncircles = NGRIDX*(NGRIDY+1); dy = (YMAX - YMIN)/((double)NGRIDY); dx = dy*0.5*sqrt(3.0); for (i = 0; i < NGRIDX; i++) for (j = 0; j < NGRIDY+1; j++) { n = (NGRIDY+1)*i + j; particles[n].xc = ((double)(i-NGRIDX/2) + 0.5)*dy; /* is +0.5 needed? */ particles[n].yc = YMIN + ((double)j - 0.5)*dy; if (((i+NGRIDX)%4 == 2)||((i+NGRIDX)%4 == 3)) particles[n].yc += 0.5*dy; particles[n].radius = MU; /* activate only circles that intersect the domain */ if ((particles[n].yc < YMAX + MU)&&(particles[n].yc > YMIN - MU)) particles[n].active = 1; else particles[n].active = 0; } break; } case (C_POOL_TABLE): { for (i=1; i<6; i++) for (j=0; j ymin - radius)&&(particles[n].xc < xmax + radius)&&(particles[n].xc > xmin - radius)) particles[n].active = 1; else particles[n].active = 0; } ncircles += NGRIDX*(NGRIDY+1); break; } case (C_POISSON_DISC): { ncirc0 = ncircles; printf("Generating new Poisson disc sample\n"); /* generate first circle */ particles[ncirc0].xc = (xmax - xmin)*(double)rand()/RAND_MAX + xmin; particles[ncirc0].yc = (ymax - ymin)*(double)rand()/RAND_MAX + ymin; active_poisson[0] = 1; // // particles[0].active = 1; n_p_active = 1; newcircles = 1; while ((n_p_active > 0)&&(ncircles < NMAXCIRCLES)) { /* randomly select an active circle */ i = rand()%(newcircles); while (!active_poisson[i]) i = rand()%(ncircles); // printf("Starting from circle %i at (%.3f,%.3f)\n", i, particles[i].xc, particles[i].yc); /* generate new candidates */ naccepted = 0; for (j=0; j= dpoisson*dpoisson); /* new circle is in domain */ far = far*(x < xmax)*(x > xmin)*(y < ymax)*(y > ymin); // far = far*(vabs(x) < LAMBDA)*(y < INITYMAX)*(y > INITYMIN); } if (far) /* accept new circle */ { printf("New circle at (%.3f,%.3f) accepted\n", x, y); particles[ncircles].xc = x; particles[ncircles].yc = y; particles[ncircles].radius = radius; particles[ncircles].active = 1; particles[ncircles].added = 1; particles[ncircles].coulomb = 1; particles[ncircles].reactive = reactive; active_poisson[ncircles] = 1; ncircles++; ncirc0++; n_p_active++; naccepted++; } // else printf("Rejected\n"); } if (naccepted == 0) /* inactivate circle i */ { // printf("No candidates work, inactivate circle %i\n", i); active_poisson[i] = 0; n_p_active--; } printf("%i active circles\n", n_p_active); } printf("Generated %i circles\n", ncircles); break; } default: { printf("Function init_circle_config not defined for this pattern \n"); } } } void init_people_config(t_person people[NMAXCIRCLES]) /* initialise particle configuration */ { t_particle particles[NMAXCIRCLES]; int n; init_particle_config(particles); for (n=0; n angle[table[j+1]]) { temp = table[j]; table[j] = table[j+1]; table[j+1] = temp; } } } /* TEST */ for (i=0; i 0.0); } int old_triangle_overlap(t_otriangle triangle1, t_otriangle triangle2, t_obstacle obstacle[NMAXOBSTACLES]) /* returns 1 if triangles overlap */ { int j, jplus, jminus; /* find first vertex equal to triangle1.i[0] */ j = 0; while (triangle2.i[j] != triangle1.i[0]) j++; if (j < 3) { jplus = (j+1)%3; jminus = (j+5)%3; if (triangle1.i[1] == triangle2.i[jplus]) return(triangle_vertex_overlap(0, 1, 2, j, jplus, jminus, triangle1, triangle2, obstacle)); else if (triangle1.i[2] == triangle2.i[jminus]) return(triangle_vertex_overlap(2, 0, 1, j, jminus, jplus, triangle1, triangle2, obstacle)); else return(0); } else { if ((triangle1.i[1] == triangle2.i[1])&&(triangle1.i[2] == triangle2.i[2])) return(triangle_vertex_overlap(1, 2, 0, 1, 2, 0, triangle1, triangle2, obstacle)); else return(0); } } int triangle_overlap(t_otriangle triangle1, t_otriangle triangle2) /* returns 1 if triangles overlap */ { int j, jplus, jminus; // printf("Testing triangles (%i, %i, %i) and (%i, %i, %i)\n", triangle1.i[0], triangle1.i[1], triangle1.i[2], triangle2.i[0], triangle2.i[1], triangle2.i[2]); /* find first vertex equal to triangle1.i[0] */ j = 0; while ((j<3)&&(triangle2.i[j] != triangle1.i[0])) j++; if (j < 3) { jplus = (j+1)%3; jminus = (j+5)%3; if (triangle1.i[1] == triangle2.i[jplus]) return(1); else if (triangle1.i[2] == triangle2.i[jminus]) return(1); else return(0); } else { if ((triangle1.i[1] == triangle2.i[1])&&(triangle1.i[2] == triangle2.i[2])) return(1); else return(0); } } int init_obstacle_facets(t_otriangle otriangle[NMAX_TRIANGLES_PER_OBSTACLE*NMAXOBSTACLES], t_ofacet ofacet[NMAXOBSTACLES]) /* group triangles into facets, for option COUPLE_OBSTACLES */ { int i, j, k, n, ik, f, nt; printf("Initializing obstacle facets from %i triangles\n", n_otriangles); // printf("1\n"); // printf("2\n"); for (i=0; i NMAXOBSTACLES) { printf("NMAXOBSTACLES should be at least %i\n", n_ofacets); exit(1); } return(n_ofacets); } double otriangle_area(t_obstacle obstacle[NMAXOBSTACLES], t_otriangle triangle) { int s1, s2, s3; s1 = triangle.i[0]; s2 = triangle.i[1]; s3 = triangle.i[2]; return(vabs(triangle_area(obstacle[s1].xc, obstacle[s1].yc, obstacle[s2].xc, obstacle[s2].yc, obstacle[s3].xc, obstacle[s3].yc))); } int init_obstacle_coupling(t_obstacle obstacle[NMAXOBSTACLES], t_otriangle otriangle[NMAX_TRIANGLES_PER_OBSTACLE*NMAXOBSTACLES], t_ofacet ofacet[NMAXOBSTACLES]) /* initialize list of neighbours, for option COUPLE_OBSTACLES */ { int i, j, n, m, table[NMAX_OBSTACLE_NEIGHBOURS], neigh_table[NMAX_OBSTACLE_NEIGHBOURS]; double dist, bdistx, bdisty, angle[NMAX_OBSTACLE_NEIGHBOURS], dist_table[NMAX_OBSTACLE_NEIGHBOURS]; short int near_boundary; printf("Initializing obstacle neighbours\n"); n_otriangles = 0; for (i=0; i= BCXMAX - bdistx)||(obstacle[i].yc <= BCYMIN + bdisty)||(obstacle[i].yc >= BCYMAX - bdisty)); switch (OBSTACLE_PINNING_TYPE) { case (OP_NPINNED): { obstacle[i].pinned = (obstacle[i].nneighb <= NMAX_OBSTACLE_PINNED); break; } case (OP_LEFT): { obstacle[i].pinned = (obstacle[i].xc <= BCXMIN + bdistx); break; } case (OP_LEFTTOPBOT): { obstacle[i].pinned = near_boundary; break; } case (OP_CORNERS): { obstacle[i].pinned = (((obstacle[i].xc <= BCXMIN + bdistx)||(obstacle[i].xc >= BCXMAX - bdistx))&&((obstacle[i].yc <= BCYMIN + bdisty)||(obstacle[i].yc >= BCYMAX - bdisty))); break; } case (OP_LEFTCORNERS): { obstacle[i].pinned = ((near_boundary)&&((obstacle[i].xc < 0.0)||(obstacle[i].xc > BCXMAX - bdistx))); break; } case (OP_BNDRY_STEP): { obstacle[i].pinned = ((near_boundary)&&(obstacle[i].chessboard == 0)); } } } printf("Found %i facets\n", n_ofacets); // for (i=0; i OBSXMAX) obstacle[n].xc += (OBSXMAX - OBSXMIN); obstacle[n].radius = OBSTACLE_RADIUS; obstacle[n].active = 1; n++; } nobstacles = n; break; } case (O_SQUARE): { dx = (OBSXMAX - OBSXMIN)/(double)NOBSX; dy = (OBSYMAX - OBSYMIN)/(double)NOBSY; n = 0; if ((ADD_FIXED_SEGMENTS)&&(SEGMENT_PATTERN == S_CYLINDER)) /* pseudo-cylindrical boundary conditions (e.g. for Hall effect) */ { jmin = 1; jmax = NOBSY-1; } else { jmin = 0; jmax = NOBSY; } for (i=0; i OBSXMAX) obstacle[n].xc += (OBSXMAX - OBSXMIN); obstacle[n].radius = OBSTACLE_RADIUS; obstacle[n].active = 1; obstacle[n].chessboard = (i+j)%BDRY_PINNING_STEP; n++; } nobstacles = n; break; } case (O_SQUARE_TWOMASSES): { dx = (XMAX - XMIN)/(double)NOBSX; dy = (YMAX - YMIN)/(double)NOBSY; n = 0; if ((ADD_FIXED_SEGMENTS)&&(SEGMENT_PATTERN == S_CYLINDER)) /* pseudo-cylindrical boundary conditions (e.g. for Hall effect) */ { jmin = 1; jmax = NOBSY-1; } else { jmin = 0; jmax = NOBSY; } for (i=0; i XMAX) obstacle[n].xc += (XMAX - XMIN); obstacle[n].radius = OBSTACLE_RADIUS; if ((i+j)%2 == 0) { obstacle[n].radius *= sqrt(2.0); obstacle[n].mass *= 2.0; } obstacle[n].active = 1; obstacle[n].chessboard = (i+j)%BDRY_PINNING_STEP; n++; } nobstacles = n; break; } case (O_SIDES): { n = 0; for (i = 0; i < 9; i++) for (j = 0; j < 5; j++) if ((i == 0)||(i == 8)||(j == 0)||(j == 4)) { obstacle[n].xc = BCXMIN + 0.125*((double)i)*(BCXMAX - BCXMIN); obstacle[n].yc = BCYMIN + 0.25*((double)j)*(BCYMAX - BCYMIN); obstacle[n].radius = OBSTACLE_RADIUS; obstacle[n].active = 1; n++; } nobstacles = n; break; } case (O_SIDES_B): { n = 0; nx = 16; ny = 8; dx = (BCXMAX - BCXMIN)/(double)nx; dy = (BCYMAX - BCYMIN)/(double)ny; for (i = 0; i < nx+1; i++) for (j = 0; j < ny+1; j++) if ((i == 0)||(i == nx)||(j == 0)||(j == ny)) { obstacle[n].xc = BCXMIN + (double)i*dx; obstacle[n].yc = BCYMIN + (double)j*dy; obstacle[n].radius = OBSTACLE_RADIUS; obstacle[n].active = 1; n++; } nobstacles = n; break; } case (O_SIEVE): { n = 0; width = 1.0; dx = width/12.0; dy = dx*0.6; for (i = 0; i < 13; i++) { obstacle[n].xc = -1.0 + (double)i*dx; obstacle[n].yc = 0.7 - (double)i*dy; obstacle[n].radius = 0.95*OBSTACLE_RADIUS; obstacle[n].omega = 1.25*OBSTACLE_OMEGA; obstacle[n].active = 1; n++; } for (i = 0; i < 13; i++) { obstacle[n].xc = -1.0 + (double)i*dx; obstacle[n].yc = 0.2 - (double)i*dy; obstacle[n].radius = 1.2*OBSTACLE_RADIUS; obstacle[n].omega = OBSTACLE_OMEGA; obstacle[n].active = 1; n++; } for (i = 0; i < 13; i++) { obstacle[n].xc = -1.0 + (double)i*dx; obstacle[n].yc = -0.3 - (double)i*dy; obstacle[n].radius = 1.6*OBSTACLE_RADIUS; obstacle[n].omega = 0.75*OBSTACLE_OMEGA; obstacle[n].active = 1; n++; } nobstacles = n; break; } case (O_SIEVE_B): { n = 0; width = 1.2; dx = width/14.0; dy = dx*0.6; for (i = 0; i < 15; i++) { obstacle[n].xc = -1.2 + (double)i*dx; obstacle[n].yc = 0.85 - (double)i*dy; obstacle[n].radius = 0.9*OBSTACLE_RADIUS; obstacle[n].omega = 1.25*OBSTACLE_OMEGA; obstacle[n].active = 1; n++; } for (i = 0; i < 15; i++) { obstacle[n].xc = -1.2 + (double)i*dx; obstacle[n].yc = 0.35 - (double)i*dy; obstacle[n].radius = 1.2*OBSTACLE_RADIUS; obstacle[n].omega = OBSTACLE_OMEGA; obstacle[n].active = 1; n++; } for (i = 0; i < 15; i++) { obstacle[n].xc = -1.2 + (double)i*dx; obstacle[n].yc = -0.15 - (double)i*dy; obstacle[n].radius = 1.6*OBSTACLE_RADIUS; obstacle[n].omega = 0.75*OBSTACLE_OMEGA; obstacle[n].active = 1; n++; } if (RATTLE_OBSTACLES) for (i = 0; i < n; i++) { obstacle[i].oscillate = 1; obstacle[i].period = 8; obstacle[i].amplitude = 0.0015; obstacle[i].phase = (double)i*DPI/10.0; } nobstacles = n; break; } case (O_SIEVE_LONG): { n = 0; ntot = 36; width = 1.2; dx = (XMAX - XMIN)/(double)ntot; dy = dx*0.3; for (i = 0; i < ntot + 1; i++) { obstacle[n].xc = XMIN + (double)i*dx; obstacle[n].yc = 0.6 - (double)i*dy; obstacle[n].radius = OBSTACLE_RADIUS*(2.4 - 1.6*(double)i/(double)ntot); obstacle[n].omega = OBSTACLE_OMEGA; obstacle[n].active = 1; n++; } if (RATTLE_OBSTACLES) for (i = 0; i < n; i++) { obstacle[i].oscillate = 1; obstacle[i].period = 8; obstacle[i].amplitude = 0.0018; obstacle[i].phase = (double)i*DPI/10.0; } nobstacles = n; break; } case (O_SIEVE_LONG_B): { n = 0; ntot = 45; width = 1.2; dx = (XMAX - XMIN)/(double)ntot; dy = dx*0.2; x = XMIN; y = 0.4; c = 0.019; while (x < XMAX) { obstacle[n].xc = x; obstacle[n].yc = y; obstacle[n].radius = OBSTACLE_RADIUS; obstacle[n].omega = OBSTACLE_OMEGA; obstacle[n].active = 1; x += dx; y -= dy; dx += c*dx; dy += c*dy; n++; } if (RATTLE_OBSTACLES) for (i = 0; i < n; i++) { obstacle[i].oscillate = 1; obstacle[i].period = 8; obstacle[i].amplitude = 0.0018; obstacle[i].phase = (double)i*DPI/10.0; } nobstacles = n; break; } case (O_SIEVE_VIDEO1400): { n = 0; ntot = 28; width = 1.2; dx = 2.9/(double)ntot; dy = dx*0.3; for (i = 0; i < ntot + 1; i++) { obstacle[n].xc = 2.5 + (double)i*dx; obstacle[n].yc = 0.6 - (double)i*dy; obstacle[n].radius = OBSTACLE_RADIUS*(2.4 - 1.6*(double)i/(double)ntot); obstacle[n].omega = OBSTACLE_OMEGA; obstacle[n].active = 1; n++; } if (RATTLE_OBSTACLES) for (i = 0; i < n; i++) { obstacle[i].oscillate = 1; obstacle[i].period = 8; obstacle[i].amplitude = 0.0018; obstacle[i].phase = (double)i*DPI/10.0; } nobstacles = n; break; } case (O_POISSON_DISC): { point = (t_point *)malloc(NMAXOBSTACLES*sizeof(t_point)); nobstacles = generate_poisson_discs(point, NMAXOBSTACLES, OBSXMIN, OBSXMAX, OBSYMIN, OBSYMAX, OBSTACLE_PISC_DISTANCE); printf("Generated %i obstacles\n", nobstacles); for (n=0; n XMAX) obstacle[n].xc += (XMAX - XMIN); obstacle[n].radius = OBSTACLE_RADIUS; obstacle[n].active = 1; } free(point); break; } default: { printf("Function init_obstacle_config not defined for this pattern \n"); } } if (CHARGE_OBSTACLES) for (n=0; n= 0.0) return(width); else switch (nozzle_shape) { case (NZ_STRAIGHT): return(width); case (NZ_BELL): return(sqrt(width*width - 0.5*x)); case (NZ_GLAS): return(sqrt(width*width - 1.2*x) + 1.0*x); case (NZ_CONE): return(width - (sqrt(width*width + 0.5) - width)*x); case (NZ_TRUMPET): return(width + (sqrt(width*width + LAMBDA)-width)*x*x); case (NZ_BROAD): { if (-x < 0.1) return(width - (0.5 - width)*x/0.1); else return(0.5); } case (NZ_DELAVAL): { a = 1.5; b = 0.05; return(sqrt(width*width - 0.5*x) + a*b*x*(1.0 + x)/(b + x*x)); // a = (sqrt(width*width+0.5) - width)/sqrt(0.5); // c = (a*sqrt(0.5*h) - width)/(h*h); // if (-x < h) return(width + c*x*x); // else return(width + a * sqrt(-0.5*x)); } default: return(0.0); } } void add_rocket_to_segments(t_segment segment[NMAXSEGMENTS], double x0, double y0, int rocket_shape, int nozzle_shape, int nsides, int group) /* add one or several rocket_shaped set of segments */ { int i, j, cycle = 0, nsegments0; double angle, dx, x1, y1, x2, y2, nozx, nozy, a, b, c, w; nsegments0 = nsegments; /* compute intersection point of nozzle and ellipse */ angle = -PID + DPI/(double)NPOLY; nozx = 0.7*LAMBDA*cos(angle); nozy = y0 + YMIN + LAMBDA*(1.7 + 0.7*sin(angle)); /* form of combustion chamber */ switch (rocket_shape) { case (RCK_DISC): /* circular chamber */ { for (i=1; i0)||(j0)||(j!=NYMAZE/2))&&(maze[n].west)) add_rectangle_to_segments(x1, y1, x1 - width, y1 + dy, segment, 0); if (maze[n].south) add_rectangle_to_segments(x1, y1, x1 + dx, y1 - width, segment, 0); } /* top side of maze */ add_rectangle_to_segments(YMIN + padding + MAZE_XSHIFT, YMAX - padding, YMAX - padding + MAZE_XSHIFT, YMAX - padding - width, segment, 0); /* right side of maze */ x1 = YMAX - padding + MAZE_XSHIFT; if (diag) { y1 = YMIN + padding + dy; add_rectangle_to_segments(x1, y1, x1 - width, YMAX + 10.0, segment, 0); } else { y1 = YMIN + padding + dy*((double)NYMAZE/2); add_rectangle_to_segments(x1, YMIN - 1.0, x1 - width, y1 - dy, segment, 0); add_rectangle_to_segments(x1, y1, x1 - width, YMAX + 1.0, segment, 0); } /* left side of maze */ x1 = YMIN + padding + MAZE_XSHIFT; add_rectangle_to_segments(x1, YMIN - 1.0, x1 - width, YMIN + padding, segment, 0); add_rectangle_to_segments(x1, YMAX - padding, x1 - width, YMAX + 10.0, segment, 0); if (diag) { add_rotated_angle_to_segments(XMIN, YMAX - 0.5*dy, x1, YMAX - dy - 2.0*width, width, 0, segment, 0); add_rectangle_to_segments(XMIN, YMAX - 0.5*dy, XMIN - width, YMAX + 10.0, segment, 0); } free(maze); } void translate_segments(t_segment segment[NMAXSEGMENTS], double deltax[2], double deltay[2]) /* rotates the repelling segments by given angle */ { int i, group; for (i=0; i DPI) { segment[i].angle1 -= DPI; segment[i].angle2 -= DPI; } while (segment[i].angle2 < 0.0) { segment[i].angle1 += DPI; segment[i].angle2 += DPI; } } } int add_shovel_to_belt(double x, double y, double angle, double size, double width, t_segment segment[NMAXSEGMENTS], t_belt belt) /* add showver to conveyor belt */ { int i, n, iminus, iplus; double sq2, angle1, angle2, cg, sg, t; if (nsegments + 6 > NMAXSEGMENTS) { printf("NMAXSEGMENTS too small"); return(0); } n = nsegments; // sq2 = sqrt(2.0)/2.0; cg = cos(PI/6.0); sg = sin(PI/6.0); t = size + width*(cg - 1.0)/sg; segment[n].x1 = 0.0; segment[n].y1 = 0.0; segment[n+1].x1 = 0.0; segment[n+1].y1 = size; segment[n+2].x1 = size*sg; segment[n+2].y1 = size*(1.0+cg); segment[n+3].x1 = size*sg + width*cg; segment[n+3].y1 = size*(1.0+cg) - width*sg; segment[n+4].x1 = width; segment[n+4].y1 = size*(1.0 + cg) - width*sg - t*cg; // segment[n+2].x1 = size*sq2; // segment[n+2].y1 = size*(1.0+sq2); // // segment[n+3].x1 = size*sq2 + width*sq2; // segment[n+3].y1 = size*(1.0+sq2) - width*sq2; // // segment[n+4].x1 = width; // segment[n+4].y1 = size + width*(1.0 - sqrt(2.0)); segment[n+5].x1 = width; segment[n+5].y1 = 0.0; for (i=0; i<6; i++) { segment[n+i].x2 = segment[n+i+1].x1; segment[n+i].y2 = segment[n+i+1].y1; } segment[n+6].x2 = segment[n].x1; segment[n+6].y2 = segment[n].y1; for (i=n; i n + 5) iplus = n; angle1 = argument(segment[iplus].x1 - segment[i].x1, segment[iplus].y1 - segment[i].y1) + PID; angle2 = argument(segment[i].x1 - segment[iminus].x1, segment[i].y1 - segment[iminus].y1) + PID; if (angle2 < angle1) angle2 += DPI; segment[i].angle1 = angle1; segment[i].angle2 = angle2; } for (i=n; i NMAXSEGMENTS) { printf("NMAXSEGMENTS too small"); return(0); } if (nbelts + 1 > NMAXBELTS) { printf("NMAXBELTS too small"); return(0); } segment[n].x1 = x1 - ty*width; segment[n].y1 = y1 + tx*width; for (i=0; i<7; i++) { segment[n+1+i].x1 = x2 + width*sin(-angle + (double)i*PI/6.0); segment[n+1+i].y1 = y2 + width*cos(-angle + (double)i*PI/6.0); } for (i=7; i<14; i++) { segment[n+1+i].x1 = x1 + width*sin(-angle + (double)(i-1)*PI/6.0); segment[n+1+i].y1 = y1 + width*cos(-angle + (double)(i-1)*PI/6.0); } for (i=0; i<13; i++) { segment[n+i].x2 = segment[n+i+1].x1; segment[n+i].y2 = segment[n+i+1].y1; } segment[n+13].x2 = segment[n].x1; segment[n+13].y2 = segment[n].y1; for (i=n; i n + 13) iplus = n; angle1 = argument(segment[iplus].x1 - segment[i].x1, segment[iplus].y1 - segment[i].y1) + PID; angle2 = argument(segment[i].x1 - segment[iminus].x1, segment[i].y1 - segment[iminus].y1) + PID; if (angle2 < angle1) angle2 += DPI; segment[i].angle1 = angle1; segment[i].angle2 = angle2; } nsegments += 14; conveyor_belt[nbelts].x1 = x1; conveyor_belt[nbelts].x2 = x2; conveyor_belt[nbelts].y1 = y1; conveyor_belt[nbelts].y2 = y2; conveyor_belt[nbelts].speed = speed; conveyor_belt[nbelts].width = width; conveyor_belt[nbelts].position = 0.0; conveyor_belt[nbelts].length = length; conveyor_belt[nbelts].angle = angle; conveyor_belt[nbelts].tx = tx; conveyor_belt[nbelts].ty = ty; conveyor_belt[nbelts].nshovels = nshovels; nbelts++; if (nshovels > NMAXSHOVELS) { printf("NMAXSHOVELS too small"); conveyor_belt[nbelts].nshovels = 0; return(0); } if (nshovels > 0) { beltlength = 2.0*length + DPI*width; shovel_dist = beltlength/(double)nshovels; shovel_pos = 0.0; shovel = 0; while (shovel_pos < length) { x = x1 - width*ty + shovel_pos*tx; y = y1 + width*tx + shovel_pos*ty; conveyor_belt[nbelts-1].shovel_segment[shovel] = nsegments; add_shovel_to_belt(x, y, angle, width, 0.3*width, segment, conveyor_belt[nbelts-1]); conveyor_belt[nbelts-1].shovel_pos[shovel] = shovel_pos; shovel_pos += shovel_dist; shovel++; } while (shovel_pos < length + width*PI) { beta = (shovel_pos - length)/width; x = x2 + width*cos(angle + PID - beta); y = y2 + width*sin(angle + PID - beta); conveyor_belt[nbelts-1].shovel_segment[shovel] = nsegments; add_shovel_to_belt(x, y, angle - beta, width, 0.3*width, segment, conveyor_belt[nbelts-1]); conveyor_belt[nbelts-1].shovel_pos[shovel] = shovel_pos; shovel_pos += shovel_dist; shovel++; } while (shovel_pos < 2.0*length + width*PI) { x = x2 + width*ty - (shovel_pos - length - width*PI)*tx; y = y2 - width*tx - (shovel_pos - length - width*PI)*ty; conveyor_belt[nbelts-1].shovel_segment[shovel] = nsegments; add_shovel_to_belt(x, y, angle - PI, width, 0.3*width, segment, conveyor_belt[nbelts-1]); conveyor_belt[nbelts-1].shovel_pos[shovel] = shovel_pos; shovel_pos += shovel_dist; shovel++; } while (shovel_pos < beltlength) { beta = (shovel_pos - 2.0*length - width*PI)/width; x = x1 + width*cos(angle - PID - beta); y = y1 + width*sin(angle - PID - beta); conveyor_belt[nbelts-1].shovel_segment[shovel] = nsegments; add_shovel_to_belt(x, y, angle - PI - beta, width, 0.3*width, segment, conveyor_belt[nbelts-1]); conveyor_belt[nbelts-1].shovel_pos[shovel] = shovel_pos; shovel_pos += shovel_dist; shovel++; } } return(1); } void init_segment_config(t_segment segment[NMAXSEGMENTS], t_belt conveyor_belt[NMAXBELTS]) /* initialise linear obstacle configuration */ { int i, j, cycle = 0, iminus, iplus, nsides, n, concave = 1; double angle, angle2, dangle, dx, width, height, a, b, length, xmid = 0.5*(BCXMIN + BCXMAX), lpocket, r, x, x1, y1, x2, y2, nozx, nozy, y, yy, dy, ca, sa, padding; double lfoot[NTREES][2], rfoot[NTREES][2]; /* set default to no conveyor, no torque */ for (i=0; i NMAXSEGMENTS) { printf("Error: NMAXSEGMENTS is too small\n"); exit(1); } for (i=0; i NMAXSEGMENTS) { printf("Error: NMAXSEGMENTS is too small\n"); exit(1); } for (i=0; i 7) segment[i].inactivate = 1; else segment[i].inactivate = 0; } break; } case (S_DAM_WITH_HOLE_AND_RAMP): { add_rectangle_to_segments(DAM_WIDTH, BCYMIN - 0.5, -DAM_WIDTH, BCYMIN + 0.2, segment, 0); add_rectangle_to_segments(DAM_WIDTH, BCYMIN + 0.3, -DAM_WIDTH, LAMBDA, segment, 0); add_rectangle_to_segments(DAM_WIDTH, BCYMIN + 0.2, -DAM_WIDTH, BCYMIN + 0.3, segment, 0); r = 1.0; for (i=0; i<10; i++) { angle = 0.1*PID*(double)i; dangle = 0.1*PID; x1 = XMAX - r + (r + MU)*cos(angle); y1 = YMIN + r - (r + MU)*sin(angle); x2 = XMAX - r + (r + MU)*cos(angle + dangle); y2 = YMIN + r - (r + MU)*sin(angle + dangle); add_rotated_angle_to_segments(x1, y1, x2, y2, MU, 0, segment, 0); } cycle = 0; concave = 0; /* add_rectangle_to_segments already deals with concave corners */ for (i=0; i 7)&&(i < 12)) segment[i].inactivate = 1; else segment[i].inactivate = 0; } break; } case (S_MAZE): { init_maze_segments(segment, 0); cycle = 0; for (i=0; i nsegments - 1) iplus = 0; } angle = argument(segment[iplus].x1 - segment[i].x1, segment[iplus].y1 - segment[i].y1) + PID; angle2 = argument(segment[i].x1 - segment[iminus].x1, segment[i].y1 - segment[iminus].y1) + PID; if (angle2 < angle) angle2 += DPI; segment[i].angle1 = angle; segment[i].angle2 = angle2; printf("i = %i, iplus = %i, iminus = %i, angle1 = %.0f, angle2 = %.0f\n", i, iplus, iminus, angle*360.0/DPI, angle2*360.0/DPI); } /* make copy of initial values in case of rotation/translation */ if ((ROTATE_BOUNDARY)||(MOVE_BOUNDARY)||(MOVE_SEGMENT_GROUPS)) for (i=0; i 0.05)&&(y1 < b - 0.05)); } case (RCK_RECT_HAT) : { // printf("(%.2lg,%.2lg) in rocket?\n", x, y); a = 0.5*LAMBDA; b = (0.49*PI-0.25)*LAMBDA; y1 = y - YMIN - LAMBDA; if (vabs(x) > 0.95*a) return(0); if (y1 < 0.05) return(0); if (y1 < b - 0.05) return(1); return(y1 < a + b - 0.05 - vabs(x)); // return(1); } case (RCK_RECT_BAR) : { a = 0.5*LAMBDA; b = (0.49*PI-0.25)*LAMBDA; c = 0.5*a; y1 = y - YMIN - LAMBDA; if (vabs(x) > 0.95*a) return(0); if (vabs(x) < 0.1*a) return(0); if (y1 < 0.05 + c) return(0); if (y1 < b - 0.05 + c) return(1); return(y1 < a + b - 0.05 + c - vabs(x)); } } } int in_segment_region(double x, double y, t_segment segment[NMAXSEGMENTS]) /* returns 1 if (x,y) is inside region delimited by obstacle segments */ { int i; double angle, dx, height, width, length, theta, lx, ly, x1, y1, x2, y2, padding, ca, sa, r; if (x >= BCXMAX) return(0); if (x <= BCXMIN) return(0); if (y >= BCYMAX) return(0); if (y <= BCYMIN) return(0); switch (SEGMENT_PATTERN) { case (S_CUP): { angle = APOLY*PID; dx = (BCYMAX - BCYMIN)/tan(angle); if (y < BCYMAX - (BCYMAX - BCYMIN)*(x - BCXMIN)/dx) return(0); if (y < BCYMAX - (BCYMAX - BCYMIN)*(BCXMAX - x)/dx) return(0); } case (S_HOURGLASS): { angle = APOLY*PID; width = 2.5*MU; height = 2.5*MU; x = vabs(x); y = vabs(y); if ((x >= width)&&(x - width >= (y - height)*(BCXMAX - width)/(BCYMAX - height))) return(0); return(1); } case (S_PENTA): { height = 0.5*(BCYMAX - BCYMIN); width = height/sqrt(3.0); if (y < BCYMIN + height*(1.0 - (x - BCXMIN)/width)) return(0); if (y > BCYMAX - height*(1.0 - (x - BCXMIN)/width)) return(0); return(1); } case (S_CENTRIFUGE): { angle = argument(x,y); theta = DPI/(double)NPOLY; while (angle > theta) angle -= theta; while (angle < 0.0) angle += theta; if (angle < 0.1) return(0); if (angle > 0.9) return(0); return(1); } case (S_POLY_ELLIPSE): { if (x*x + y*y/(LAMBDA*LAMBDA) < 0.95) return(1); else return(0); } case (S_CENTRIFUGE_RND): { if (module2(x,y) > 0.75) return(0); angle = argument(x,y); theta = DPI/(double)NPOLY; while (angle > theta) angle -= theta; while (angle < 0.0) angle += theta; if (angle < 0.1) return(0); if (angle > 0.9) return(0); return(1); } case (S_CENTRIFUGE_LEAKY): { if (module2(x,y) > 0.75) return(0); angle = argument(x,y); theta = DPI/(double)NPOLY; while (angle > theta) angle -= theta; while (angle < 0.0) angle += theta; if (angle < 0.1) return(0); if (angle > 0.9) return(0); return(1); } case (S_CIRCLE_EXT): { if (module2(x - SEGMENTS_X0, y - SEGMENTS_Y0) > LAMBDA + 2.0*MU) return(1); else return(0); } case (S_TWO_CIRCLES_EXT): { if ((module2(x - SEGMENTS_X0, y - SEGMENTS_Y0) > LAMBDA + 2.0*MU)&&(module2(x + SEGMENTS_X0, y - SEGMENTS_Y0) > TWO_CIRCLES_RADIUS_RATIO*LAMBDA + 2.0*MU)) return(1); else return(0); } case (S_ROCKET_NOZZLE): { if (x < 0.0) return(0); else if (x > 1.4*LAMBDA) return(0); else { lx = 0.7*LAMBDA; ly = 0.7*LAMBDA; x -= lx; if (x*x/(lx*lx) + y*y/(ly*ly) < 0.95) return(1); } return(0); } case (S_ROCKET_NOZZLE_ROTATED): { y1 = y - ysegments[0]; if (y1 < YMIN + LAMBDA) return(0); else if (y1 > YMIN + 2.4*LAMBDA) return(0); else if (in_rocket(x - xsegments[0], y1, ROCKET_SHAPE)) return(1); // { // ly = 0.7*LAMBDA; // lx = 0.7*LAMBDA; // x1 = x - xsegments[0]; // y1 -= YMIN + 1.7*LAMBDA; // if (x1*x1/(lx*lx) + y1*y1/(ly*ly) + MU*MU < 0.925) return(1); // } return(0); } case (S_TWO_ROCKETS): { y1 = y - ysegments[0]; y2 = y - ysegments[1]; if ((y1 < YMIN + LAMBDA)&&(y2 < YMIN + LAMBDA)) return(0); else if ((y1 > YMIN + 3.5*LAMBDA)&&(y2 > YMIN + 3.5*LAMBDA)) return(0); else { if (in_rocket(x - xsegments[0], y1, ROCKET_SHAPE_B)) return(1); if (in_rocket(x - xsegments[1], y2, ROCKET_SHAPE)) return(1); } return(0); } case (S_DAM): { if (vabs(x) > DAM_WIDTH) return(1); else if (y > LAMBDA) return(1); else return(0); } case (S_EXT_RECTANGLE): { padding = 0.1; if (vabs(x) > LAMBDA + padding) return(1); else if (vabs(y) > 0.1*LAMBDA + padding) return(1); else return(0); } case (S_EXT_CIRCLE_RECT): { padding = 0.1; if ((vabs(x - SEGMENTS_X0) < LAMBDA + padding)&&(vabs(y - SEGMENTS_Y0) < 0.1*LAMBDA + padding)) return(0); else if (module2(x + SEGMENTS_X0, y - SEGMENTS_Y0) < 0.5*LAMBDA + padding) return(0); else return(1); } case (S_BIN_OPENING): { padding = 3.0*MU; if (y < 0.0) return(1); if ((y > 1.0 - LAMBDA + padding)&&(vabs(x) < LAMBDA - padding)) return(1); return(0); } case (S_POLYGON_EXT): { padding = 3.0*MU; if (in_polygon(x, y, LAMBDA + padding, NPOLY, APOLY)) return(0); else return(1); } case (S_WEDGE_EXT): { padding = 3.0*MU; angle = AWEDGE*PID; ca = cos(APOLY*PID); sa = sin(APOLY*PID); x1 = x*ca - y*sa; y1 = x*sa + y*ca; if (vabs(y1) - padding > (LAMBDA-x1)*sin(angle)/(1.0+cos(angle))) return(1); if (vabs(y1) + padding < -x1*tan(angle)) return(1); return(0); } case (S_MIXER): { padding = 1.5*MU; r = module2(x,y); if (r > LAMBDA + padding) return(1); if (r < 2.0*padding) return(0); for (i=0; i 0.0)&&(vabs(y1) < 0.025 + padding)) return(0); } return(1); } case (S_AIRFOIL): { if (vabs(x) > LAMBDA + 4.0*MU) return(1); padding = 0.35; angle = APOLY*PID; ca = cos(angle); sa = sin(angle); x1 = x*ca - y*sa; y1 = x*sa + y*ca; y1 += 0.5*x1*x1; return(x1*x1 + 25.0*y1*y1 > LAMBDA*LAMBDA*(1.0 + padding)); } case (S_COANDA): { return(vabs(y - SEGMENTS_Y0 - 0.25*cos(PI*x/XMAX)) > 0.1); } case (S_COANDA_SHORT): { x1 = XMIN + 0.1; length = 0.85*(XMAX - XMIN - 0.2); if (x > x1 + length + 0.1) return(1); return(vabs(y - SEGMENTS_Y0 + 0.2*cos(DPI*(x-x1)/length)) > 0.05); } case (S_CYLINDER): { if (y > YMAX - 0.05 - MU) return(0); if (y < YMIN + 0.05 + MU) return(0); return(1); } case (S_MAZE): { for (i=0; i DPI) { segment[i].angle1 -= DPI; segment[i].angle2 -= DPI; } } } } /* Computation of interaction force */ double lennard_jones_force(double r, t_particle particle) { int i; double rmin = 0.01, rplus, ratio = 1.0; if (r > REPEL_RADIUS*particle.radius) return(0.0); else { if (r > rmin) rplus = r; else rplus = rmin; // ratio = ipow(particle.eq_dist*particle.radius/rplus, 6); ratio = ipow(particle.eq_dist*particle.radius/rplus, 3); ratio = ratio*ratio; return((ratio - 2.0*ratio*ratio)/rplus); } } double harmonic_force(double r, t_particle particle) { int i; double rmin = 0.01, rplus, ratio = 1.0; if (r > REPEL_RADIUS*particle.radius) return(0.0); else { // if (r > rmin) rplus = r; // else rplus = rmin; // ratio = rplus/particle.eq_dist*particle.radius; ratio = r/particle.eq_dist*particle.radius; if (ratio < 1.0) return(ratio - 1.0); else return(0.0); } } double coulomb_force(double r, t_particle particle) { int i; double rmin = 0.01, rplus, ratio = 1.0; if (r > REPEL_RADIUS*particle.radius) return(0.0); else { // if (r > rmin) rplus = r; // else rplus = rmin; // ratio = rplus/particle.eq_dist*particle.radius; ratio = r/particle.eq_dist*particle.radius; return(-1.0/(ratio*ratio + rmin*rmin)); // if (ratio < 1.0) return(ratio - 1.0); // else return(0.0); } } void aniso_lj_force(double r, double ca, double sa, double ca_rel, double sa_rel, double force[2], t_particle particle) { int i; double rmin = 0.01, rplus, ratio = 1.0, c2, s2, c4, s4, a, aprime, f1, f2; if (r > REPEL_RADIUS*particle.radius) { force[0] = 0.0; force[1] = 0.0; } else { if (r > rmin) rplus = r; else rplus = rmin; // ratio = ipow(particle.eq_dist*particle.radius/rplus, 6); ratio = ipow(particle.eq_dist*particle.radius/rplus, 3); ratio = ratio*ratio; /* cos(2phi) and sin(2phi) */ c2 = ca_rel*ca_rel - sa_rel*sa_rel; s2 = 2.0*ca_rel*sa_rel; /* cos(4phi) and sin(4phi) */ c4 = c2*c2 - s2*s2; s4 = 2.0*c2*s2; a = 0.5*(9.0 - 7.0*c4); aprime = 14.0*s4; f1 = ratio*(a - ratio)/rplus; f2 = ratio*aprime/rplus; force[0] = f1*ca - f2*sa; force[1] = f1*sa + f2*ca; } } void penta_lj_force(double r, double ca, double sa, double ca_rel, double sa_rel, double force[2], t_particle particle) { int i; double rmin = 0.01, rplus, ratio = 1.0, c2, s2, c4, s4, c5, s5, a, aprime, f1, f2; static double a0, b0; static int first = 1; if (first) { a0 = cos(0.1*PI) + 0.5; b0 = a0 - 1.0; first = 0; } if (r > REPEL_RADIUS*particle.radius) { force[0] = 0.0; force[1] = 0.0; } else { if (r > rmin) rplus = r; else rplus = rmin; // ratio = ipow(particle.eq_dist*particle.radius/rplus, 6); ratio = ipow(particle.eq_dist*particle.radius/rplus, 3); ratio = ratio*ratio; /* cos(2phi) and sin(2phi) */ c2 = ca_rel*ca_rel - sa_rel*sa_rel; s2 = 2.0*ca_rel*sa_rel; /* cos(4phi) and sin(4phi) */ c4 = c2*c2 - s2*s2; s4 = 2.0*c2*s2; /* cos(5phi) and sin(5phi) */ c5 = ca_rel*c4 - sa_rel*s4; s5 = sa_rel*c4 + ca_rel*s4; a = a0 - b0*c5; aprime = 5.0*b0*s5; f1 = ratio*(a - ratio)/rplus; f2 = ratio*aprime/rplus; force[0] = f1*ca - f2*sa; force[1] = f1*sa + f2*ca; } } double old_golden_ratio_force(double r, t_particle particle) /* potential with two minima at distances whose ratio is the golden ratio Phi */ /* old version that does not work very well */ { int i; double x, y, z, rplus, ratio = 1.0, phi, a, phi3; static int first = 1; static double rmin, b, c, d; if (first) { rmin = 0.5*particle.radius; phi = 0.5*(1.0 + sqrt(5.0)); phi3 = 1.0/(phi*phi*phi); a = 0.66; b = 1.0 + phi3 + a; d = phi3*a; c = phi3 + a + d; // b = 7.04; // c = 13.66; // d = 6.7; first = 0; printf("a = %.4lg, b = %.4lg, c = %.4lg, d = %.4lg\n", a, b, c, d); } if (r > REPEL_RADIUS*particle.radius) return(0.0); else { if (r > rmin) rplus = r; else rplus = rmin; x = particle.eq_dist*particle.radius/rplus; y = x*x*x; z = d - c*y + b*y*y - y*y*y; return(x*z/rplus); } } double golden_ratio_force(double r, t_particle particle) /* potential with two minima at distances whose ratio is the golden ratio Phi */ /* piecewise polynomial/LJ version */ { int i; double x, rplus, xm6, y1; static int first = 1; static double rmin, phi, a, h1, h2, phi6; if (first) { rmin = 0.5*particle.radius; phi = 0.5*(1.0 + sqrt(5.0)); a = 1.2; h1 = 1.0; /* inner potential well depth */ h2 = 10.0; /* outer potential well depth */ phi6 = ipow(phi, 6); first = 0; } if (r > REPEL_RADIUS*particle.radius) return(0.0); else { if (r > rmin) rplus = r; else rplus = rmin; x = rplus/(particle.eq_dist*particle.radius); // xm6 = 1.0/ipow(x, 6); xm6 = 1.0/ipow(x, 3); xm6 = xm6*xm6; if (x <= 1.0) return(12.0*h1*xm6*(1.0 - xm6)/x); else if (x <= a) { y1 = ipow(a - 1.0, 3); return(6.0*h1*(x - 1.0)*(a - x)/y1); } else if (x <= phi) { y1 = ipow(phi - a, 3); return(6.0*h2*(x - a)*(x - phi)/y1); } else return(12.0*h2*phi6*(1.0 - phi6*xm6)*xm6/x); } } void dipole_lj_force(double r, double ca, double sa, double ca_rel, double sa_rel, double force[2], t_particle particle) { int i; double rmin = 0.01, rplus, ratio = 1.0, a, aprime, f1, f2; if (r > REPEL_RADIUS*MU) { force[0] = 0.0; force[1] = 0.0; } else { if (r > rmin) rplus = r; else rplus = rmin; // ratio = ipow(particle.eq_dist*particle.radius/rplus, 6); ratio = ipow(particle.eq_dist*particle.radius/rplus, 3); ratio = ratio*ratio; a = 1.0 + 0.25*ca_rel; aprime = -0.25*sa_rel; f1 = ratio*(a - ratio)/rplus; f2 = ratio*aprime/rplus; force[0] = f1*ca - f2*sa; force[1] = f1*sa + f2*ca; } } void quadrupole_lj_force(double r, double ca, double sa, double ca_rel, double sa_rel, double force[2], t_particle particle) { int i; double rmin = 0.01, rplus, ratio = 1.0, a, aprime, f1, f2, ca2, sa2, x, y, dplus, dminus; static int first = 1; static double a0, b0, aplus, aminus; if (first) { dplus = cos(0.2*PI)*cos(0.1*PI); // dminus = 0.8*dplus; dminus = QUADRUPOLE_RATIO*dplus; aplus = ipow(1.0/dplus, 6); aminus = ipow(1.0/dminus, 6); // aminus = ipow(cos(0.2*PI)*(0.25 + 0.5*sin(0.1*PI)), 6); a0 = 0.5*(aplus + aminus); b0 = 0.5*(aplus - aminus); first = 0; } if (r > REPEL_RADIUS*particle.radius) { force[0] = 0.0; force[1] = 0.0; } else { if (r > rmin) rplus = r; else rplus = rmin; // ratio = ipow(particle.eq_dist*particle.radius/rplus, 6); ratio = ipow(particle.eq_dist*particle.radius/rplus, 3); ratio = ratio*ratio; /* cos(2*phi) and sin(2*phi) */ ca2 = ca_rel*ca_rel - sa_rel*sa_rel; sa2 = 2.0*ca_rel*sa_rel; a = a0 + b0*ca2; // if (a == 0.0) a = 1.0e-10; aprime = -2.0*b0*sa2; f1 = ratio*(a - ratio)/rplus; f2 = ratio*aprime/rplus; force[0] = f1*ca - f2*sa; force[1] = f1*sa + f2*ca; } } void quadrupole_lj_force2(double r, double ca, double sa, double ca_rel, double sa_rel, double force[2], t_particle particle) { int i; double rmin = 0.01, rplus, ratio = 1.0, a, aprime, f1, f2, ca2, sa2, x, y, eqdist; static int first = 1; static double aplus, aminus, a0, b0; if (first) { aplus = ipow(cos(0.2*PI)*cos(0.1*PI), 6); aminus = 0.1*aplus; // aminus = 0.0; // aminus = -2.0*ipow(cos(0.2*PI)*(0.5*sin(0.1*PI)), 6); // aminus = ipow(cos(0.2*PI)*(0.25 + 0.5*sin(0.1*PI)), 6); a0 = 0.5*(aplus + aminus); b0 = 0.5*(aplus - aminus); first = 0; } if (r > REPEL_RADIUS*particle.radius) { force[0] = 0.0; force[1] = 0.0; } else { if (r > rmin) rplus = r; else rplus = rmin; /* correct distance */ // ratio = ipow(particle.eq_dist*particle.radius/rplus, 6); ratio = ipow(particle.eq_dist*particle.radius/rplus, 3); ratio = ratio*ratio; /* cos(2*phi) and sin(2*phi) */ ca2 = ca_rel*ca_rel - sa_rel*sa_rel; sa2 = 2.0*ca_rel*sa_rel; a = a0 + b0*ca2; if (a == 0.0) a = 1.0e-10; aprime = -2.0*b0*sa2; // f1 = ratio*(a - ratio)/rplus; // f2 = ratio*aprime/rplus; f1 = ratio*(aplus - ratio)/(rplus); f2 = ratio*(aminus - ratio)/(rplus); // force[0] = f1*ca_rel - f2*sa_rel; // force[1] = f1*sa_rel + f2*ca_rel; force[0] = f1*ca - f2*sa; force[1] = f1*sa + f2*ca; } } double water_torque(double r, double ca, double sa, double ca_rel, double sa_rel, double ck_rel, double sk_rel) /* compute torque of water molecule #k on water molecule #j (for interaction I_LJ_WATER) - OLD VERSION */ { double c1p, c1m, c2p, c2m, s2p, s2m, s21, s21p, s21m, c21, c21p, c21m, torque; double r2, rd, rd2, rr[3][3]; static double cw = -0.5, sw = 0.866025404, delta = 1.5*MU, d2 = 2.25*MU*MU; int i, j; c1p = ck_rel*cw - sk_rel*sw; c1m = ck_rel*cw + sk_rel*sw; c2p = ca_rel*cw - sa_rel*sw; c2m = ca_rel*cw + sa_rel*sw; s2p = sa_rel*cw + ca_rel*sw; s2m = sa_rel*cw - ca_rel*sw; s21 = sa_rel*ck_rel - ca_rel*sk_rel; c21 = ca_rel*ck_rel + sa_rel*sk_rel; s21p = s21*cw - c21*sw; s21m = s21*cw + c21*sw; c21p = c21*cw + s21*sw; c21m = c21*cw - s21*sw; r2 = r*r; rd = 2.0*r*delta; rd2 = r2 + d2; rr[0][0] = r; rr[0][1] = sqrt(rd2 + rd*c2p); rr[0][2] = sqrt(rd2 + rd*c2m); rr[1][0] = sqrt(rd2 - rd*c1p); rr[2][0] = sqrt(rd2 - rd*c1m); rr[1][1] = sqrt(r2 + rd*(c2p - c1p) + 2.0*d2*(1.0 - c21)); rr[1][2] = sqrt(r2 + rd*(c2m - c1p) + 2.0*d2*(1.0 - c21m)); rr[2][1] = sqrt(r2 + rd*(c2p - c1m) + 2.0*d2*(1.0 - c21p)); rr[2][2] = sqrt(r2 + rd*(c2m - c1m) + 2.0*d2*(1.0 - c21)); for (i=0; i<3; i++) for (j=0; j<3; j++) { if (rr[i][j] < 1.0e-4) rr[i][j] = 1.0e-4; rr[i][j] = rr[i][j]*rr[i][j]*rr[i][j]; } torque = rd*(s2p/rr[0][1] + s2m/rr[0][2]); torque += -0.5*rd*(s2p/rr[1][1] + s2p/rr[2][1] + s2m/rr[1][2] + s2m/rr[2][2]); torque += d2*(s21/rr[1][1] + s21/rr[2][2] + s21m/rr[1][2] + s21p/rr[2][1]); return(torque); } double water_force(double r, double ca, double sa, double ca_rel, double sa_rel, double ck_rel, double sk_rel, double f[2]) /* compute force and torque of water molecule #k on water molecule #j (for interaction I_LJ_WATER) */ { double x1[3], y1[3], x2[3], y2[3], rr[3][3], dx[3][3], dy[3][3], fx[3][3], fy[3][3], m[3][3], torque = 0.0; static double cw[3], sw[3], q[3], d[3], delta = 1.25*MU, dmin = 0.5*MU, fscale = 1.0; int i, j; static int first = 1; if (first) { cw[0] = 1.0; cw[1] = -0.5; cw[2] = -0.5; sw[0] = 0.0; sw[1] = 866025404; sw[2] = -866025404; /* sines and cosines of angles */ q[0] = -2.0; q[1] = 1.0; q[2] = 1.0; /* charges */ d[0] = 0.5*delta; d[1] = delta; d[2] = delta; /* distances to center */ first = 0; } /* positions of O and H atoms */ for (i=0; i<3; i++) { x1[i] = d[i]*(ca_rel*cw[i] - sa_rel*sw[i]); y1[i] = d[i]*(ca_rel*sw[i] + sa_rel*cw[i]); x2[i] = r + d[i]*(ck_rel*cw[i] - sk_rel*sw[i]); y2[i] = d[i]*(ck_rel*sw[i] + sk_rel*cw[i]); } /* relative positions */ for (i=0; i<3; i++) for (j=0; j<3; j++) { dx[i][j] = x2[j] - x1[i]; dy[i][j] = y2[j] - y1[i]; rr[i][j] = module2(dx[i][j], dy[i][j]); if (rr[i][j] < dmin) rr[i][j] = dmin; rr[i][j] = ipow(rr[i][j],3); // rr[i][j] = rr[i][j]*rr[i][j]*rr[i][j]; } /* forces between particles */ for (i=0; i<3; i++) for (j=0; j<3; j++) { fx[i][j] = -q[i]*q[j]*dx[i][j]/rr[i][j]; fy[i][j] = -q[i]*q[j]*dy[i][j]/rr[i][j]; } /* torques between particles */ for (i=0; i<3; i++) for (j=0; j<3; j++) { m[i][j] = x1[i]*fy[i][j] - y1[i]*fx[i][j]; } /* total force */ f[0] = 0.0; f[1] = 0.0; for (i=0; i<3; i++) for (j=0; j<3; j++) { f[0] += fscale*fx[i][j]; f[1] += fscale*fy[i][j]; torque += fscale*m[i][j]; } return(torque); } void segment_interaction(t_particle part_i, t_particle part_k, double ca, double sa, double distance, double ffm[3], int charge) /* compute interaction for I_SEGMENT interaction */ /* computes force and torque of particle k on particle i */ { int s; double x0, y0, x, y, fx, fy, cb, sb, cc, sc, cab, sab, ccb, scb, r; double width, torque, mu_i, mu_k, r3, endfact, cfact; fx = 0.0; fy = 0.0; torque = 0.0; mu_i = part_i.radius; mu_k = part_k.radius; width = 0.2*mu_i; cb = cos(part_i.angle); sb = sin(part_i.angle); cc = cos(part_k.angle); sc = sin(part_k.angle); cab = ca*cb + sa*sb; sab = sa*cb - ca*sb; ccb = cc*cb + sc*sb; scb = sc*cb - cc*sb; /* relative coordinates of particle k in basis where particle i is horizontal */ x0 = distance*cab; y0 = distance*sab; for (s=-1; s<=1; s+=2) { x = x0 + (double)s*mu_k*ccb; y = y0 + (double)s*mu_k*scb; if (vabs(y) < width) { if (vabs(x) < mu_i) { // if (((double)s*scb < 0.0)&&(y > 0.0)) fy += y - width; if (y > 0.0) fy += y - width; else fy += y + width; } else if ((x >= mu_i)&&(x < mu_i+width)) { r = module2(x-mu_i, y) + 1.0e-10; if (r < width) { fx -= (width - r)*(x-mu_i)/r; fy -= (width - r)*y/r; } } else if ((x <= -mu_i)&&(x > -mu_i-width)) { r = module2(x+mu_i, y) + 1.0e-10; if (r < width) { fx -= (width - r)*(x+mu_i)/r; fy -= (width - r)*y/r; } } torque += x*fy - y*fx; } if (charge) /* add Coulomb force between ends */ { cfact = CHARGE*KREPEL/KSPRING_OBSTACLE; r = module2(x-mu_i,y) + 1.0e-2; r3 = r*r*r; fx += cfact*(double)s*(x+mu_i)/r3; fy += cfact*(double)s*y/r3; torque += mu_i*fx; r = module2(x+mu_i,y) + 1.0e-2; r3 = r*r*r; fx -= cfact*(double)s*(x+mu_i)/r3; fy -= cfact*(double)s*y/r3; torque -= mu_i*fy; } } ffm[0] = cb*fx - sb*fy; ffm[1] = sb*fx + cb*fy; ffm[2] = torque; } void polygon_interaction(t_particle part_i, t_particle part_k, double ca, double sa, double distance, double ffm[3], int charge) /* compute interaction for I_POLYGON interaction */ /* computes force and torque of particle k on particle i */ { int s, k, s0, s1; double x0, y0, x, y, f, fx, fy, cb, sb, cab, sab, r, r0, z2, phi, ckt, skt, d, fx1, fy1, dmax2; double torque, mu_i, mu_k, r3, width, gamma_beta, rot, cr, sr, xx[NPOLY], yy[NPOLY], dd[NPOLY], z, z0, z02, d1, dmin; static double theta, twotheta, ctheta, stheta, cornerx[NPOLY+1], cornery[NPOLY+1], nx[NPOLY], ny[NPOLY]; static int first = 1; if (first) { theta = PI/(double)NPOLY; twotheta = 2.0*theta; ctheta = cos(theta); stheta = sin(theta); for (s=0; s width) f = width; fx1 = f*nx[k]; fy1 = f*ny[k]; fx -= fx1; fy -= fy1; torque += x*fy1 - y*fx1; // fx -= f*nx[k]; // fy -= f*ny[k]; // torque += x*fy - y*fx; } else for (s1=0; s1 width) f = width; fx1 -= f*(x - mu_i*cornerx[s1])/d1; fy1 -= f*(y - mu_i*cornery[s1])/d1; } fx += fx1; fy += fy1; torque += x*fy1 - y*fx1; } } } ffm[0] = cb*fx - sb*fy; ffm[1] = sb*fx + cb*fy; ffm[2] = torque; } int dna_charge(int type1, int type2, int coulomb) /* returns -1 for matching nucleotide ends, 1 for non-matching, and 0 else */ { int min, max; // if (coulomb == 0) return(0); if (type1 == 0) return(0); if (type2 == 0) return(0); if ((type1 == 7)&&(type2 == 7)) return(2); if (type1 == type2) return(1); if (type1 > 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, f1, ffm[3]; static double dxhalf = 0.5*(BCXMAX - BCXMIN), dyhalf = 0.5*(BCYMAX - BCYMIN); int wwrapx, wwrapy, twrapx, twrapy; if (BOUNDARY_COND == BC_GENUS_TWO) { dxhalf *= 0.75; dyhalf *= 0.75; } x1 = particle[i].xc; y1 = particle[i].yc; x2 = particle[k].xc; y2 = particle[k].yc; wwrapx = ((BOUNDARY_COND == BC_KLEIN)||(BOUNDARY_COND == BC_BOY)||(BOUNDARY_COND == BC_GENUS_TWO))&&(vabs(x2 - x1) > dxhalf); wwrapy = ((BOUNDARY_COND == BC_BOY)||(BOUNDARY_COND == BC_GENUS_TWO))&&(vabs(y2 - y1) > dyhalf); twrapx = ((BOUNDARY_COND == BC_KLEIN)||(BOUNDARY_COND == BC_BOY))&&(vabs(x2 - x1) > dxhalf); twrapy = (BOUNDARY_COND == BC_BOY)&&(vabs(y2 - y1) > dyhalf); switch (particle[k].interaction) { case (I_COULOMB): { f = -krepel/(1.0e-8 + distance*distance); force[0] = f*ca; force[1] = f*sa; break; } case (I_LENNARD_JONES): { f = krepel*lennard_jones_force(distance, particle[k]); force[0] = f*ca; force[1] = f*sa; break; } case (I_LJ_DIRECTIONAL): { aniso_lj_force(distance, ca, sa, ca_rel, sa_rel, ff, particle[k]); force[0] = krepel*ff[0]; force[1] = krepel*ff[1]; break; } case (I_LJ_PENTA): { penta_lj_force(distance, ca, sa, ca_rel, sa_rel, ff, particle[k]); force[0] = krepel*ff[0]; force[1] = krepel*ff[1]; break; } case (I_GOLDENRATIO): { f = krepel*golden_ratio_force(distance, particle[k]); force[0] = f*ca; force[1] = f*sa; break; } case (I_LJ_DIPOLE): { dipole_lj_force(distance, ca, sa, ca_rel, sa_rel, ff, particle[k]); force[0] = krepel*ff[0]; force[1] = krepel*ff[1]; break; } case (I_LJ_QUADRUPOLE): { quadrupole_lj_force(distance, ca, sa, ca_rel, sa_rel, ff, particle[k]); force[0] = krepel*ff[0]; force[1] = krepel*ff[1]; break; } case (I_LJ_WATER): { f = krepel*lennard_jones_force(distance, particle[k]); force[0] = f*ca; force[1] = f*sa; break; } case (I_VICSEK): { force[0] = 0.0; force[1] = 0.0; break; } case (I_VICSEK_REPULSIVE): { f = krepel*coulomb_force(distance, particle[k]); // f = krepel*harmonic_force(distance, particle[k]); // f = krepel*lennard_jones_force(distance, particle[k]); force[0] = f*ca; force[1] = f*sa; break; } case (I_VICSEK_SPEED): { f = cos(0.5*(particle[k].angle - particle[i].angle)); force[0] = f*KSPRING_VICSEK*(particle[k].vx - particle[i].vx); force[1] = f*KSPRING_VICSEK*(particle[k].vy - particle[i].vy); break; } case (I_VICSEK_SHARK): { if (particle[k].type == particle[i].type) { f = cos(0.5*(particle[k].angle - particle[i].angle)); force[0] = f*KSPRING_VICSEK*(particle[k].vx - particle[i].vx); force[1] = f*KSPRING_VICSEK*(particle[k].vy - particle[i].vy); } else if (particle[i].type != 2) { f = krepel*coulomb_force(distance, particle[k]); force[0] = f*ca; force[1] = f*sa; } else { if (VICSEK_REPULSION > 0.0) { // f = VICSEK_REPULSION*harmonic_force(distance, particle[k]); f = VICSEK_REPULSION*coulomb_force(distance, particle[k]); force[0] = f*ca; force[1] = f*sa; } else { force[0] = 0.0; force[1] = 0.0; } } break; } case (I_COULOMB_LJ): { charge = particle[i].charge*particle[k].charge; f = -100.0*krepel*charge/(1.0e-12 + distance*distance); if (charge <= 0.0) f += COULOMB_LJ_FACTOR*krepel*lennard_jones_force(distance, particle[k]); force[0] = f*ca; 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; } case (I_SEGMENT): { segment_interaction(particle[i], particle[k], ca, sa, distance, ffm, 0); force[0] = KSPRING_OBSTACLE*ffm[0]; force[1] = KSPRING_OBSTACLE*ffm[1]; *torque = KTORQUE*ffm[2]; f = 0.1*krepel*lennard_jones_force(distance, particle[k]); force[0] += f*ca; force[1] += f*sa; break; } case (I_SEGMENT_CHARGED): { segment_interaction(particle[i], particle[k], ca, sa, distance, ffm, 1); force[0] = KSPRING_OBSTACLE*ffm[0]; force[1] = KSPRING_OBSTACLE*ffm[1]; *torque = KTORQUE*ffm[2]; // printf("Force between particles %i and %i: (%.3lg, %.3lg)\n", i, k, force[0], force[1]); // f = 0.1*krepel*lennard_jones_force(distance, particle[k]); // force[0] += f*ca; // force[1] += f*sa; break; } case (I_POLYGON): { polygon_interaction(particle[i], particle[k], ca, sa, distance, ffm, 0); force[0] = KSPRING_OBSTACLE*ffm[0]; force[1] = KSPRING_OBSTACLE*ffm[1]; *torque = KTORQUE*ffm[2]; f = 0.1*krepel*lennard_jones_force(distance, particle[k]); force[0] += f*ca; force[1] += f*sa; break; } case (I_POLYGON_ALIGN): { polygon_interaction(particle[i], particle[k], ca, sa, distance, ffm, 0); force[0] = KSPRING_OBSTACLE*ffm[0]; force[1] = KSPRING_OBSTACLE*ffm[1]; *torque = KTORQUE*ffm[2]; f = 0.1*krepel*lennard_jones_force(distance, particle[k]); force[0] += f*ca; force[1] += f*sa; break; } } if (ROTATION) { dist_scaled = distance/(particle[i].spin_range*particle[i].radius); switch (particle[k].interaction) { case (I_LJ_WATER): { ck = cos(particle[k].angle); sk = sin(particle[k].angle); ck_rel = ca*ck + sa*sk; sk_rel = sa*ck - ca*sk; // *torque = (-3.0*ca_rel*sk_rel + 2.0*sa_rel*ck_rel)/(1.0e-12 + dist_scaled*dist_scaled*dist_scaled); // *torque = water_torque(distance, ca, sa, ca_rel, sa_rel, ck_rel, sk_rel); // *torque = (0.5*sin(angle) + 0.5*sin(2.0*angle) - 0.45*sin(3.0*angle))/(1.0e-12 + dist_scaled*dist_scaled*dist_scaled); *torque = water_force(distance, ca, sa, ca_rel, sa_rel, ck_rel, sk_rel, ff); force[0] += ff[0]; force[1] += ff[1]; // printf("force = (%.3lg, %.3lg)\n", ff[0], ff[1]); break; } case (I_VICSEK): { if (dist_scaled > 1.0) *torque = 0.0; else if (twrapx||twrapy) *torque = sin(-particle[k].angle - particle[i].angle); else *torque = sin(particle[k].angle - particle[i].angle); break; } case (I_VICSEK_REPULSIVE): { if (dist_scaled > 1.0) *torque = 0.0; else if (twrapx||twrapy) *torque = sin(-particle[k].angle - particle[i].angle); else *torque = sin(particle[k].angle - particle[i].angle); break; } case (I_VICSEK_SPEED): { if (dist_scaled > 1.0) *torque = 0.0; else if (twrapx||twrapy) *torque = sin(-particle[k].angle - particle[i].angle); else *torque = sin(particle[k].angle - particle[i].angle); break; } case (I_VICSEK_SHARK): { if (dist_scaled > 10.0) *torque = 0.0; else if (particle[k].type == particle[i].type) /* fish adjusting direction */ { if (twrapx||twrapy) *torque = sin(-particle[k].angle - particle[i].angle); else *torque = sin(particle[k].angle - particle[i].angle); } else if (particle[k].type == 2) /* fish fleeing a shark */ { alpha = argument(ca,sa); particle[i].angle = alpha + PI; *torque = 0.0; } else /* shark tracking fish */ { *torque = cos(particle[k].angle)*sa - sin(particle[k].angle)*ca; } break; } case (I_SEGMENT): { /* Do nothing, torque already computed */ break; } case (I_SEGMENT_CHARGED): { /* Do nothing, torque already computed */ break; } case (I_POLYGON): { /* Do nothing, torque already computed */ break; } case (I_POLYGON_ALIGN): { spin_f = particle[i].spin_freq; if (twrapx||twrapy) *torque += sin(spin_f*(-particle[k].angle - particle[i].angle))/(1.0e-8 + dist_scaled*dist_scaled); else { if (NPOLY%2 == 0) *torque += sin(spin_f*(particle[k].angle - particle[i].angle))/(1.0e-8 + dist_scaled*dist_scaled); else *torque -= sin(spin_f*(particle[k].angle - particle[i].angle))/(1.0e-8 + dist_scaled*dist_scaled); } break; } default: { spin_f = particle[i].spin_freq; if (twrapx||twrapy) *torque = sin(spin_f*(-particle[k].angle - particle[i].angle))/(1.0e-8 + dist_scaled*dist_scaled); else *torque = sin(spin_f*(particle[k].angle - particle[i].angle))/(1.0e-8 + dist_scaled*dist_scaled); } } if (particle[i].type == particle[k].type) { if (particle[i].type == 0) *torque *= KTORQUE; else *torque *= KTORQUE_B; } else *torque *= KTORQUE_DIFF; } else *torque = 0.0; if ((distance < NBH_DIST_FACTOR*particle[i].radius)&&(k != i)) return(1); // if ((distance < NBH_DIST_FACTOR*particle[i].radius)) return(1); else return(0); } int compute_repelling_force(int i, int j, double force[2], double *torque, t_particle* particle, double krepel) /* compute repelling force of neighbour #j on particle #i */ /* returns 1 if distance between particles is smaller than NBH_DIST_FACTOR*MU */ { double distance, ca, sa, cj, sj, ca_rel, sa_rel, f[2], ff[2], torque1, ck, sk, ck_rel, sk_rel; static double distmin = 10.0*((XMAX - XMIN)/HASHX + (YMAX - YMIN)/HASHY); int interact, k; if (BOUNDARY_COND == BC_GENUS_TWO) distmin *= 2.0; k = particle[i].hashneighbour[j]; distance = module2(particle[i].deltax[j], particle[i].deltay[j]); /* for monitoring purposes */ // if (distance > distmin) // { // printf("i = %i, hashcell %i, j = %i, hashcell %i\n", i, particle[i].hashcell, k, particle[k].hashcell); // printf("X = (%.3lg, %.3lg)\n", particle[i].xc, particle[i].yc); // printf("Y = (%.3lg, %.3lg) d = %.3lg\n", particle[k].xc, particle[k].yc, distance); // } if ((distance == 0.0)||(i == k)) { force[0] = 0.0; force[1] = 0.0; *torque = 0.0; return(1); } else if (distance > REPEL_RADIUS*particle[i].radius) { force[0] = 0.0; force[1] = 0.0; *torque = 0.0; return(0); } else { /* to avoid numerical problems, assign minimal value to distance */ if (distance < 0.1*particle[i].radius) distance = 0.1*particle[i].radius; ca = (particle[i].deltax[j])/distance; sa = (particle[i].deltay[j])/distance; /* compute relative angle in case particles can rotate */ if (ROTATION) { cj = cos(particle[j].angle); sj = sin(particle[j].angle); ca_rel = ca*cj + sa*sj; sa_rel = sa*cj - ca*sj; } else { ca_rel = ca; sa_rel = sa; } interact = compute_particle_interaction(i, k, f, torque, particle, distance, krepel, ca, sa, ca_rel, sa_rel); if (SYMMETRIZE_FORCE) { torque1 = *torque; // compute_particle_interaction(k, i, ff, torque, particle, distance, krepel, ca, sa, ca_rel, sa_rel); ck = cos(particle[j].angle); sk = sin(particle[j].angle); ck_rel = ca*ck + sa*sk; sk_rel = sa*ck - ca*sk; compute_particle_interaction(k, i, ff, torque, particle, distance, krepel, -ca, -sa, -ck_rel, -sk_rel); force[0] = 0.5*(f[0] - ff[0]); force[1] = 0.5*(f[1] - ff[1]); *torque = 0.5*(torque1 - *torque); // *torque = 0.5*(*torque + torque1); } else { force[0] = f[0]; force[1] = f[1]; } // printf("force = (%.3lg, %.3lg), torque = %.3lg\n", f[0], f[1], *torque); return(interact); } } int resample_particle(int n, int maxtrials, t_particle particle[NMAXCIRCLES]) /* resample y coordinate of particle n, returns 1 if no collision is created */ { double x, y, dist, dmin = 10.0; int i, j, closeby = 0, success = 0, trials = 0; while ((!success)&&(trials < maxtrials)) { success = 1; x = particle[n].xc - MU*(double)rand()/RAND_MAX; y = 0.95*(BCYMIN + (BCYMAX - BCYMIN)*(double)rand()/RAND_MAX); i = 0; while ((success)&&(i 0) { printf("Deactivating particle cluster %i\n", i); particle[i].active = 0; for (l=0; l nmolecules) nmolecules = m + 1; np = molecule[m].nparticles; // printf("np = %i\n", np); if (np < NPARTNERS+1) { molecule[m].particle[np] = i; molecule[m].nparticles++; } molecule[m].added = particle[i].added; } printf("Found %i molecules\n", nmolecules); fprintf(lj_log, "Found %i molecules\n", nmolecules); /* for debugging */ for (m=0; m NMAXCIRCLES) { printf("Error: NMAXCIRCLES is too small\n"); exit(1); } nmolecules = ncircles; 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= NMAXCIRCLES)) { printf("Cannot add particle at (%.3lg, %.3lg)\n", x, y); return(0); } else { i = ncircles; particle[i].type = type; particle[i].xc = x; particle[i].yc = y; particle[i].radius = MU; // particle[i].radius = MU*sqrt(mass); particle[i].active = 1; particle[i].neighb = 0; particle[i].diff_neighb = 0; particle[i].thermostat = 1; particle[i].charge = CHARGE; particle[i].energy = 0.0; particle[i].emean = 0.0; particle[i].dirmean = 0.0; if (RANDOM_RADIUS) particle[i].radius = particle[i].radius*(RANDOM_RADIUS_MIN + RANDOM_RADIUS_RANGE*((double)rand()/RAND_MAX)); particle[i].mass_inv = 1.0/mass; if (particle[i].type == 0) particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT; else particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT; if ((RANDOM_RADIUS)&&(ADAPT_MASS_TO_RADIUS > 0.0)) particle[i].mass_inv *= 1.0/(1.0 + pow(particle[i].radius/MU, ADAPT_MASS_TO_RADIUS)); if ((RANDOM_RADIUS)&&(ADAPT_DAMPING_TO_RADIUS > 0.0)) particle[i].damping = 1.0 + ADAPT_DAMPING_FACTOR*pow(particle[i].radius/MU, ADAPT_DAMPING_TO_RADIUS); particle[i].vx = vx; particle[i].vy = vy; particle[i].energy = (particle[i].vx*particle[i].vx + particle[i].vy*particle[i].vy)*particle[i].mass_inv; particle[i].angle = DPI*(double)rand()/RAND_MAX; particle[i].omega = OMEGA_INITIAL*gaussian(); // printf("Particle[%i].omega = %.4lg\n", i, particle[i].omega); // if (particle[i].type == 1) // { // particle[i].interaction = INTERACTION_B; // particle[i].eq_dist = EQUILIBRIUM_DIST_B; // particle[i].spin_range = SPIN_RANGE_B; // particle[i].spin_freq = SPIN_INTER_FREQUENCY_B; // } if ((PLOT == P_INITIAL_POS)||(PLOT_B == P_INITIAL_POS)) { switch (INITIAL_POS_TYPE) { case (IP_X): { particle[i].color_hue = 360.0*(particle[i].xc - INITXMIN)/(INITXMAX - INITXMIN); break; } case (IP_Y): { particle[i].color_hue = 360.0*(particle[i].yc - INITYMIN)/(INITYMAX - INITYMIN); break; } } } else if ((PLOT == P_NUMBER)||(PLOT_B == P_NUMBER)) particle[i].color_hue = 360.0*(double)(i%N_PARTICLE_COLORS)/(double)N_PARTICLE_COLORS; if ((PAIR_PARTICLES)&&(type == 0)) { angle = DPI*(double)rand()/(double)RAND_MAX; for (k=0; k 0.0) { hue = ENERGY_HUE_MIN + (ENERGY_HUE_MAX - ENERGY_HUE_MIN)*ej/PARTICLE_EMAX; if (hue > ENERGY_HUE_MIN) hue = ENERGY_HUE_MIN; if (hue < ENERGY_HUE_MAX) hue = ENERGY_HUE_MAX; } *radius = particle.radius; *width = BOUNDARY_WIDTH; break; } case (P_NEIGHBOURS): { hue = neighbour_color(particle.neighb); *radius = particle.radius; *width = BOUNDARY_WIDTH; break; } case (P_BONDS): { hue = neighbour_color(particle.neighb); *radius = particle.radius; *width = 1; break; } case (P_ANGLE): { angle = particle.angle; hue = angle*particle.spin_freq/DPI; hue -= (double)((int)hue); huex = (DPI - angle)*particle.spin_freq/DPI; huex -= (double)((int)huex); angle = PI - angle; if (angle < 0.0) angle += DPI; huey = angle*particle.spin_freq/DPI; huey -= (double)((int)huey); hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*(hue); huex = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*(huex); huey = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*(huey); *radius = particle.radius; *width = BOUNDARY_WIDTH; break; } case (P_TYPE): { hue = type_hue(particle.type); *radius = particle.radius; *width = BOUNDARY_WIDTH; break; } case (P_DIRECTION): { hue = argument(particle.vx, particle.vy); if (hue > DPI) hue -= DPI; if (hue < 0.0) hue += DPI; hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*(hue)/DPI; *radius = particle.radius; *width = BOUNDARY_WIDTH; break; } case (P_DIRECT_ENERGY): { hue = argument(particle.vx, particle.vy); if (hue > DPI) hue -= DPI; if (hue < 0.0) hue += DPI; hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*(hue)/DPI; if (particle.energy < 0.1*PARTICLE_EMAX) lum = 10.0*particle.energy/PARTICLE_EMAX; else lum = 1.0; *radius = particle.radius; *width = BOUNDARY_WIDTH; break; } case (P_ANGULAR_SPEED): { hue = 160.0*(1.0 + tanh(SLOPE*particle.omega)); // printf("omega = %.3lg, hue = %.3lg\n", particle[j].omega, hue); *radius = particle.radius; *width = BOUNDARY_WIDTH; break; } case (P_DIFF_NEIGHB): { hue = (double)(particle.diff_neighb+1)/(double)(particle.neighb+1); // hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*hue; // hue = 180.0*(1.0 + hue); hue = 20.0 + 320.0*hue; *radius = particle.radius; *width = BOUNDARY_WIDTH; break; } case (P_THERMOSTAT): { if (particle.thermostat) hue = 30.0; else hue = 270.0; *radius = particle.radius; *width = BOUNDARY_WIDTH; break; } case (P_INITIAL_POS): { hue = particle.color_hue; *radius = particle.radius; *width = BOUNDARY_WIDTH; break; } case (P_NUMBER): { hue = particle.color_hue; *radius = particle.radius; *width = BOUNDARY_WIDTH; break; } case (P_EMEAN): { ej = particle.emean; if (ej > 0.0) { hue = ENERGY_HUE_MIN + (ENERGY_HUE_MAX - ENERGY_HUE_MIN)*ej/PARTICLE_EMAX; if (hue > ENERGY_HUE_MIN) hue = ENERGY_HUE_MIN; if (hue < ENERGY_HUE_MAX) hue = ENERGY_HUE_MAX; } *radius = particle.radius; *width = BOUNDARY_WIDTH; break; } case (P_EMEAN_DENSITY): { ej = particle.emean; cl = particle.cluster; ej *= PARTICLE_MASS/cluster[cl].mass; if (ej > 0.0) { hue = ENERGY_HUE_MIN + (ENERGY_HUE_MAX - ENERGY_HUE_MIN)*ej/PARTICLE_EMAX; if (hue > ENERGY_HUE_MIN) hue = ENERGY_HUE_MIN; if (hue < ENERGY_HUE_MAX) hue = ENERGY_HUE_MAX; } *radius = particle.radius; *width = BOUNDARY_WIDTH; break; } case (P_LOG_EMEAN): { ej = particle.emean; if (ej > 0.0) { ej = log(ej/PARTICLE_EMIN)/log(PARTICLE_EMAX/PARTICLE_EMIN); if (ej < 0.0) ej = 0.0; else if (ej > 1.0) ej = 1.0; hue = ENERGY_HUE_MIN + (ENERGY_HUE_MAX - ENERGY_HUE_MIN)*ej; // if (hue > ENERGY_HUE_MIN) hue = ENERGY_HUE_MIN; // if (hue < ENERGY_HUE_MAX) hue = ENERGY_HUE_MAX; } *radius = particle.radius; *width = BOUNDARY_WIDTH; break; } case (P_DIRECT_EMEAN): { hue = particle.dirmean + COLOR_HUESHIFT*PI; // printf("dirmean = %.3lg\n", particle.dirmean); if (hue > DPI) hue -= DPI; if (hue < 0.0) hue += DPI; hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*(hue)/DPI; ej = particle.emean; if (ej < 0.5*PARTICLE_EMAX) lum = 2.0*ej/PARTICLE_EMAX; else lum = 1.0; *radius = particle.radius; *width = BOUNDARY_WIDTH; break; } case (P_NOPARTICLE): { hue = 0.0; lum = 1.0; *radius = particle.radius; *width = BOUNDARY_WIDTH; break; } case (P_NPARTNERS): { hue = partner_color(particle.npartners); *radius = particle.radius; *width = BOUNDARY_WIDTH; break; } case (P_CHARGE): { hue = (-tanh(0.5*particle.charge)+1.0)*180.0; *radius = particle.radius; *width = BOUNDARY_WIDTH; break; } case (P_MOL_ANGLE): { p = particle.p0; q = particle.p1; x = other_particle[q].xc - other_particle[p].xc; y = other_particle[q].yc - other_particle[p].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); 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; *radius = particle.radius; *width = BOUNDARY_WIDTH; break; } case (P_CLUSTER): { // cluster = (double)(particle.cluster)/(double)(ncircles); ccluster = (double)(particle.cluster_color)/(double)(ncircles); ccluster -= (double)((int)ccluster); hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*ccluster; *radius = particle.radius; *width = BOUNDARY_WIDTH; break; } case (P_CLUSTER_SIZE): { // cluster = (double)(particle.cluster)/(double)(ncircles); ccluster = 1.0 - 5.0/((double)particle.cluster_size + 4.0); hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*ccluster; *radius = particle.radius; *width = BOUNDARY_WIDTH; break; } case (P_CLUSTER_SELECTED): { cl = particle.cluster; if (cluster[cl].selected) hue = COLOR_HUE_CLUSTER_SELECTED; else hue = COLOR_HUE_CLUSTER_NOT_SELECTED; *radius = particle.radius; *width = BOUNDARY_WIDTH; break; } case (P_COLLISION): { hue = (double)particle.collision; if (hue > 0.0) hue = atan(0.25*(0.03*hue + 1.0))/PID; // { // hue += 10.0; // hue *= 1.0/(40.0 + hue); // } hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*hue; *radius = particle.radius; *width = BOUNDARY_WIDTH; break; } case (P_RADIUS): { // hue = atan(5.0*(particle.radius/MU - 0.75))/PID; hue = (particle.radius/MU - RANDOM_RADIUS_MIN)/RANDOM_RADIUS_RANGE; // hue = 0.5*(hue + 1.0); hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*hue; *radius = particle.radius; *width = BOUNDARY_WIDTH; break; } } switch (plot) { case (P_KINETIC): { // hsl_to_rgb_turbo(hue, 0.9, 0.5, rgb); // hsl_to_rgb_turbo(hue, 0.9, 0.5, rgbx); // hsl_to_rgb_turbo(hue, 0.9, 0.5, rgby); hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_EKIN); hsl_to_rgb_palette(hue, 0.9, 0.5, rgbx, COLOR_PALETTE_EKIN); hsl_to_rgb_palette(hue, 0.9, 0.5, rgby, COLOR_PALETTE_EKIN); break; } case (P_BONDS): { hsl_to_rgb_turbo(hue, 0.9, 0.5, rgb); hsl_to_rgb_turbo(hue, 0.9, 0.5, rgbx); hsl_to_rgb_turbo(hue, 0.9, 0.5, rgby); break; } case (P_DIRECTION): { hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_DIRECTION); hsl_to_rgb_palette(hue, 0.9, 0.5, rgbx, COLOR_PALETTE_DIRECTION); hsl_to_rgb_palette(hue, 0.9, 0.5, rgby, COLOR_PALETTE_DIRECTION); break; } case (P_ANGLE): { hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_ANGLE); hsl_to_rgb_palette(hue, 0.9, 0.5, rgbx, COLOR_PALETTE_ANGLE); hsl_to_rgb_palette(hue, 0.9, 0.5, rgby, COLOR_PALETTE_ANGLE); break; } case (P_DIRECT_ENERGY): { hsl_to_rgb_twilight(hue, 0.9, 0.5, rgb); hsl_to_rgb_twilight(hue, 0.9, 0.5, rgbx); hsl_to_rgb_twilight(hue, 0.9, 0.5, rgby); for (i=0; i<3; i++) { rgb[i] *= lum; rgbx[i] *= lum; rgby[i] *= lum; } break; } case (P_DIFF_NEIGHB): { hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_DIFFNEIGH); hsl_to_rgb_palette(hue, 0.9, 0.5, rgbx, COLOR_PALETTE_DIFFNEIGH); hsl_to_rgb_palette(hue, 0.9, 0.5, rgby, COLOR_PALETTE_DIFFNEIGH); // hsl_to_rgb_twilight(hue, 0.9, 0.5, rgb); // hsl_to_rgb_twilight(hue, 0.9, 0.5, rgbx); // hsl_to_rgb_twilight(hue, 0.9, 0.5, rgby); break; } case (P_INITIAL_POS): { hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_INITIAL_POS); hsl_to_rgb_palette(hue, 0.9, 0.5, rgbx, COLOR_PALETTE_INITIAL_POS); hsl_to_rgb_palette(hue, 0.9, 0.5, rgby, COLOR_PALETTE_INITIAL_POS); break; } case (P_EMEAN): { hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_EKIN); hsl_to_rgb_palette(hue, 0.9, 0.5, rgbx, COLOR_PALETTE_EKIN); hsl_to_rgb_palette(hue, 0.9, 0.5, rgby, COLOR_PALETTE_EKIN); break; } case (P_LOG_EMEAN): { hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_EKIN); hsl_to_rgb_palette(hue, 0.9, 0.5, rgbx, COLOR_PALETTE_EKIN); hsl_to_rgb_palette(hue, 0.9, 0.5, rgby, COLOR_PALETTE_EKIN); break; } case (P_DIRECT_EMEAN): { hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_DIRECTION); hsl_to_rgb_palette(hue, 0.9, 0.5, rgbx, COLOR_PALETTE_DIRECTION); hsl_to_rgb_palette(hue, 0.9, 0.5, rgby, COLOR_PALETTE_DIRECTION); for (i=0; i<3; i++) { rgb[i] *= lum; rgbx[i] *= lum; rgby[i] *= lum; } break; } case (P_CHARGE): { hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_CHARGE); hsl_to_rgb_palette(hue, 0.9, 0.5, rgbx, COLOR_PALETTE_CHARGE); hsl_to_rgb_palette(hue, 0.9, 0.5, rgby, COLOR_PALETTE_CHARGE); break; } case (P_MOL_ANGLE): { hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_ANGLE); hsl_to_rgb_palette(hue, 0.9, 0.5, rgbx, COLOR_PALETTE_ANGLE); hsl_to_rgb_palette(hue, 0.9, 0.5, rgby, COLOR_PALETTE_ANGLE); break; } case (P_CLUSTER): { hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_CLUSTER); hsl_to_rgb_palette(hue, 0.9, 0.5, rgbx, COLOR_PALETTE_CLUSTER); hsl_to_rgb_palette(hue, 0.9, 0.5, rgby, COLOR_PALETTE_CLUSTER); break; } case (P_CLUSTER_SIZE): { hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_CLUSTER_SIZE); hsl_to_rgb_palette(hue, 0.9, 0.5, rgbx, COLOR_PALETTE_CLUSTER_SIZE); hsl_to_rgb_palette(hue, 0.9, 0.5, rgby, COLOR_PALETTE_CLUSTER_SIZE); break; } case (P_CLUSTER_SELECTED): { hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_CLUSTER_SELECTED); hsl_to_rgb_palette(hue, 0.9, 0.5, rgbx, COLOR_PALETTE_CLUSTER_SELECTED); hsl_to_rgb_palette(hue, 0.9, 0.5, rgby, COLOR_PALETTE_CLUSTER_SELECTED); break; } default: { hsl_to_rgb(hue, 0.9, 0.5, rgb); hsl_to_rgb(hue, 0.9, 0.5, rgbx); hsl_to_rgb(hue, 0.9, 0.5, rgby); } } } void set_segment_group_color(int group, double lum, double rgb[3]) { switch (group) { case (1): { hsl_to_rgb_palette(270.0, 0.9, 0.5, rgb, COLOR_PALETTE); break; } case (2): { hsl_to_rgb_palette(90.0, 0.9, 0.5, rgb, COLOR_PALETTE); break; } default: { rgb[0] = 1.0; rgb[1] = 1.0; rgb[2] = 1.0; } } glColor3f(lum*rgb[0], lum*rgb[1], lum*rgb[2]); } void set_segment_pressure_color(double pressure, double lum, double rgb[3]) { double hue; // if (pressure < 0.0) pressure = 0.0; hue = SEGMENT_HUE_MIN + (SEGMENT_HUE_MAX - SEGMENT_HUE_MIN)*pressure/SEGMENT_PMAX; if (hue > SEGMENT_HUE_MIN) hue = SEGMENT_HUE_MIN; else if (hue < SEGMENT_HUE_MAX) hue = SEGMENT_HUE_MAX; // hsl_to_rgb_turbo(hue, 0.9, lum, rgb); hsl_to_rgb_palette(hue, 0.9, lum, rgb, COLOR_PALETTE_PRESSURE); glColor3f(rgb[0], rgb[1], rgb[2]); } void draw_altitude_lines() { int i, i1, i2; double x, y; glColor3f(0.5, 0.5, 0.5); glLineWidth(1); i1 = (int)(YMIN + ytrack) - 1.0; i2 = (int)(YMAX + ytrack) + 1.0; for (i = i1; i < i2; i++) { y = (double)i - ytrack; draw_line(BCXMIN, y, XMAX - 1.8, y); } i1 = (int)(XMIN + xtrack) - 1.0; i2 = (int)(XMAX - 1.8 + xtrack) + 1.0; for (i = i1; i < i2; i++) { x = (double)i - xtrack; if (x < XMAX - 1.8) draw_line(x, YMIN, x, YMAX); } } void draw_one_triangle(t_particle particle, int same_table[9*HASHMAX], int p, int q, int nsame) { double x, y, dx, dy; int k; x = particle.xc + (double)p*(BCXMAX - BCXMIN); y = particle.yc + (double)q*(BCYMAX - BCYMIN); if (TRACK) { x -= xtrack; y -= ytrack; } glBegin(GL_TRIANGLE_FAN); glVertex2d(x, y); for (k=0; k 1) { if (plot == P_TYPE) { hue = type_hue(t0); hsl_to_rgb(hue, 0.9, 0.3, rgb); } else { compute_particle_colors(particle[i], cluster, plot, rgb, rgbx, rgby, &radius, &width, particle); } glColor3f(rgb[0], rgb[1], rgb[2]); if (PERIODIC_BC) for (p=-1; p<2; p++) for (q=-1; q<2; q++) draw_one_triangle(particle[i], same_table, p, q, nsame); else draw_one_triangle(particle[i], same_table, 0, 0, nsame); } } printf("\n"); } void draw_one_particle_links(t_particle particle) /* draw links of one particle */ { int i, j, k; double x1, x2, y1, y2, length, linkcolor,periodx, periody, xt1, yt1, xt2, yt2; glLineWidth(LINK_WIDTH); // if (particle.active) // { // radius = particle[j].radius; for (k = 0; k < particle.hash_nneighb; k++) { x1 = particle.xc; if (CENTER_VIEW_ON_OBSTACLE) x1 -= xshift; y1 = particle.yc; x2 = x1 + particle.deltax[k]; y2 = y1 + particle.deltay[k]; length = module2(particle.deltax[k], particle.deltay[k])/particle.radius; if (COLOR_BONDS) { if (length < 1.5) linkcolor = 1.0; else linkcolor = 1.0 - 0.75*(length - 1.5)/(NBH_DIST_FACTOR - 1.5); glColor3f(linkcolor, linkcolor, linkcolor); } if (length < 1.0*NBH_DIST_FACTOR) draw_line(x1, y1, x2, y2); } } int draw_special_particle(t_particle particle, double xc1, double yc1, double radius, double angle, int nsides, double rgb[3], short int fill) /* draw special particles shapes, for chemical reactions */ { double x1, y1, x2, y2, omega, r; int wsign, i, j; switch(RD_REACTION) { case (CHEM_AAB): { if (particle.type == 2) { for (wsign = -1; wsign <= 1; wsign+=2) { x1 = xc1 + (double)wsign*0.7*radius*cos(particle.angle); y1 = yc1 + (double)wsign*0.7*radius*sin(particle.angle); if (fill) draw_colored_polygon(x1, y1, 0.7*radius, nsides, angle + APOLY*PID, rgb); else draw_polygon(x1, y1, 0.7*radius, nsides, angle + APOLY*PID); } return(0); } break; } case (CHEM_ABC): { if (particle.type == 3) { for (wsign = -1; wsign <= 1; wsign+=2) { x1 = xc1 + (double)wsign*0.7*radius*cos(particle.angle); y1 = yc1 + (double)wsign*0.7*radius*sin(particle.angle); if (wsign == 1) { if (fill) draw_colored_polygon(x1, y1, 1.2*MU_B, nsides, angle + APOLY*PID, rgb); else draw_polygon(x1, y1, 1.2*MU_B, nsides, angle + APOLY*PID); } else { if (fill) draw_colored_polygon(x1, y1, 1.2*MU, nsides, angle + APOLY*PID, rgb); else draw_polygon(x1, y1, 1.2*MU, nsides, angle + APOLY*PID); } } return(0); } break; } case (CHEM_A2BC): { if (particle.type == 3) { draw_colored_polygon(xc1, yc1, 1.2*MU_B, nsides, angle + APOLY*PID, rgb); for (wsign = -1; wsign <= 1; wsign+=2) { x1 = xc1 + 1.5*radius*cos(particle.angle + 0.6*(double)wsign*PI); y1 = yc1 + 1.5*radius*sin(particle.angle + 0.6*(double)wsign*PI); if (fill) draw_colored_polygon(x1, y1, 1.2*MU, nsides, angle + APOLY*PID, rgb); else draw_polygon(x1, y1, 1.2*MU, nsides, angle + APOLY*PID); } return(0); } break; } case (CHEM_CATALYSIS): { if (particle.type == 3) { for (wsign = -1; wsign <= 1; wsign+=2) { x1 = xc1 + (double)wsign*0.7*radius*cos(particle.angle); y1 = yc1 + (double)wsign*0.7*radius*sin(particle.angle); if (fill) draw_colored_polygon(x1, y1, 1.2*MU, nsides, angle + APOLY*PID, rgb); else draw_polygon(x1, y1, 1.2*MU, nsides, angle + APOLY*PID); } return(0); } break; } case (CHEM_AUTOCATALYSIS): { if (particle.type == 2) { for (wsign = -1; wsign <= 1; wsign+=2) { x1 = xc1 + (double)wsign*MU*0.7*cos(particle.angle); y1 = yc1 + (double)wsign*MU*0.7*sin(particle.angle); if (fill) draw_colored_polygon(x1, y1, MU, nsides, angle + APOLY*PID, rgb); else draw_polygon(x1, y1, MU, nsides, angle + APOLY*PID); } return(0); } break; } case (CHEM_BAA): { if (particle.type == 2) { for (wsign = -1; wsign <= 1; wsign+=2) { x1 = xc1 + (double)wsign*1.2*radius*cos(particle.angle); y1 = yc1 + (double)wsign*1.2*radius*sin(particle.angle); if (fill) draw_colored_polygon(x1, y1, 0.9*radius, nsides, angle + APOLY*PID, rgb); else draw_polygon(x1, y1, 0.9*radius, nsides, angle + APOLY*PID); } return(0); } break; } case (CHEM_AABAA): { if (particle.type == 2) { for (wsign = -1; wsign <= 1; wsign+=2) { x1 = xc1 + (double)wsign*0.7*radius*cos(particle.angle); y1 = yc1 + (double)wsign*0.7*radius*sin(particle.angle); if (fill) draw_colored_polygon(x1, y1, 0.9*radius, nsides, angle + APOLY*PID, rgb); else draw_polygon(x1, y1, 0.9*radius, nsides, angle + APOLY*PID); } return(0); } break; } case (CHEM_POLYMER): { if (particle.type >= 3) { omega = DPI/(double)(particle.type - 2); draw_colored_polygon(xc1, yc1, 1.2*MU_B, nsides, angle + APOLY*PID, rgb); for (i=0; i= 3) { omega = DPI/(double)(particle.type - 2); draw_colored_polygon(xc1, yc1, 1.2*MU_B, nsides, angle + APOLY*PID, rgb); for (i=0; i= 2) { omega = DPI/(double)(particle.type - 1); draw_colored_polygon(xc1, yc1, MU, nsides, angle + APOLY*PID, rgb); for (i=0; i= 2)&&(particle.type <= 4)) { omega = DPI/(double)(particle.type - 1); if (fill) draw_colored_polygon(xc1, yc1, MU, nsides, angle + APOLY*PID, rgb); else draw_polygon(xc1, yc1, MU, nsides, angle + APOLY*PID); for (i=0; i 0.0)) { pradius = particle.radius; if (ROTATION) { angle = angle + APOLY*PID; ca = cos(angle); sa = sin(angle); } glLineWidth(2); glColor3f(1.0, 1.0, 1.0); x1 = xc1 - pradius*ca; y1 = yc1 - pradius*sa; x2 = xc1 + pradius*ca; y2 = yc1 + pradius*sa; draw_line(x1, y1, x2, y2); x1 = xc1 - pradius*sa; y1 = yc1 + pradius*ca; x2 = xc1 + pradius*sa; y2 = yc1 - pradius*ca; draw_line(x1, y1, x2, y2); } if ((DRAW_MINUS)&&(particle.charge < 0.0)) { pradius = particle.radius; if (ROTATION) { angle = angle + APOLY*PID; ca = cos(angle); sa = sin(angle); } glLineWidth(2); glColor3f(1.0, 1.0, 1.0); x1 = xc1 - pradius*ca; y1 = yc1 - pradius*sa; x2 = xc1 + pradius*ca; y2 = yc1 + pradius*sa; draw_line(x1, y1, x2, y2); } } glLineWidth(width); glColor3f(1.0, 1.0, 1.0); if (REACTION_DIFFUSION) cont = draw_special_particle(particle, xc1, yc1, radius, angle, nsides, rgb, 0); if ((particle.interaction == I_LJ_QUADRUPOLE)||(particle.interaction == I_LJ_DIPOLE)||(particle.interaction == I_VICSEK)||(particle.interaction == I_VICSEK_REPULSIVE)||(particle.interaction == I_VICSEK_SPEED)||(particle.interaction == I_VICSEK_SHARK)) draw_rhombus(xc1, yc1, radius, angle + APOLY*PID); else if (particle.interaction == I_SEGMENT) draw_rotated_rect(xc1, yc1, particle.radius, 0.1*particle.radius, angle + APOLY*PID); else if (particle.interaction == I_SEGMENT_CHARGED) { draw_rotated_rect(xc1, yc1, particle.radius, 0.1*particle.radius, angle + APOLY*PID); /* TODO: draw ends */ } else if (cont) { if (nsides == NSEG) draw_circle_precomp(xc1, yc1, radius); else draw_polygon(xc1, yc1, radius, nsides, angle + APOLY*PID); } if (particle.interaction == I_LJ_WATER) for (wsign = -1; wsign <= 1; wsign+=2) { wangle = particle.angle + (double)wsign*DPI/3.0; x1 = xc1 + particle.radius*cos(wangle); y1 = yc1 + particle.radius*sin(wangle); draw_colored_polygon(x1, y1, 0.5*radius, nsides, angle + APOLY*PID, rgb); glColor3f(1.0, 1.0, 1.0); draw_polygon(x1, y1, 0.5*radius, nsides, angle + APOLY*PID); } // printf("Particle radius %.3f, pradius %.3f, mass %.3f\n", radius, particle.radius, 1.0/particle.mass_inv); // sleep(1); } void draw_collisions(t_collision *collisions, int ncollisions) /* draw discs where collisions happen */ { int i, j; double rgb[3], lum, x, y, x1, y1; for (i=0; i 0) { lum = (double)collisions[i].time/(double)COLLISION_TIME; if (collisions[i].color == 0.0) for (j=0; j<3; j++) rgb[j] = lum; else hsl_to_rgb_palette(collisions[i].color, 0.9, lum, rgb, COLOR_PALETTE); x = collisions[i].x; y = collisions[i].y; if (CENTER_VIEW_ON_OBSTACLE) x1 = x - xshift; else x1 = x; if (TRACK) { x1 -= xtrack; y1 = y - ytrack; } else y1 = y; draw_colored_polygon(x1, y1, COLLISION_RADIUS*MU, NSEG, 0.0, rgb); collisions[i].time--; } } void draw_trajectory(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTICLES], int traj_position, int traj_length, t_particle *particle, t_cluster *cluster, int *tracer_n, int plot) /* draw tracer particle trajectory */ { int i, j, time, p, width; double x1, x2, y1, y2, rgb[3], rgbx[3], rgby[3], radius, lum; // blank(); glLineWidth(TRAJECTORY_WIDTH); printf("drawing trajectory\n"); if (traj_length < TRAJECTORY_LENGTH*TRACER_STEPS) for (i=0; i < traj_length-1; i++) // for (i=traj_length-2; i >= 0; i--) for (j=0; j 0.1)&&(module2(x2 - x1, y2 - y1) < 0.25*(YMAX - YMIN))) draw_line(x1, y1, x2, y2); } else { // for (i = traj_length-2; i >= traj_position + 1; i--) for (i = traj_position + 1; i < traj_length-1; i++) for (j=0; j 0.1)&&(module2(x2 - x1, y2 - y1) < 0.25*(YMAX - YMIN))) draw_line(x1, y1, x2, y2); } for (i=0; i < traj_position-1; i++) // for (i = traj_position-2; i >= 0; i--) for (j=0; j 0.1)&&(module2(x2 - x1, y2 - y1) < 0.25*(YMAX - YMIN))) draw_line(x1, y1, x2, y2); } } } void old_draw_trajectory(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTICLES], int traj_position, int traj_length, t_particle *particle, t_cluster *cluster, int *tracer_n, int plot) /* draw tracer particle trajectory */ { int i, j, time, p, width; double x1, x2, y1, y2, rgb[3], rgbx[3], rgby[3], radius, lum; // blank(); glLineWidth(TRAJECTORY_WIDTH); printf("drawing trajectory\n"); for (j=0; j= 0; i--) { x1 = trajectory[j*TRAJECTORY_LENGTH*TRACER_STEPS + i].xc; x2 = trajectory[j*TRAJECTORY_LENGTH*TRACER_STEPS + i+1].xc; y1 = trajectory[j*TRAJECTORY_LENGTH*TRACER_STEPS + i].yc; y2 = trajectory[j*TRAJECTORY_LENGTH*TRACER_STEPS + i+1].yc; time = traj_length - i; lum = 0.9 - TRACER_LUM_FACTOR*(double)time/(double)(TRAJECTORY_LENGTH*TRACER_STEPS); if (lum < 0.0) lum = 0.0; glColor3f(lum*rgb[0], lum*rgb[1], lum*rgb[2]); if ((lum > 0.1)&&(module2(x2 - x1, y2 - y1) < 0.25*(YMAX - YMIN))) draw_line(x1, y1, x2, y2); // printf("tracer = %i, (x1, y1) = (%.3lg, %.3lg), (x2, y2) = (%.3lg, %.3lg)\n", j, x1, y1, x2, y2); } else { // for (i = traj_position + 1; i < traj_length-1; i++) for (i = traj_length-2; i >= traj_position + 1; i--) { x1 = trajectory[j*TRAJECTORY_LENGTH*TRACER_STEPS + i].xc; x2 = trajectory[j*TRAJECTORY_LENGTH*TRACER_STEPS + i+1].xc; y1 = trajectory[j*TRAJECTORY_LENGTH*TRACER_STEPS + i].yc; y2 = trajectory[j*TRAJECTORY_LENGTH*TRACER_STEPS + i+1].yc; time = traj_position + traj_length - i; lum = 0.9 - TRACER_LUM_FACTOR*(double)time/(double)(TRAJECTORY_LENGTH*TRACER_STEPS); if (lum < 0.0) lum = 0.0; glColor3f(lum*rgb[0], lum*rgb[1], lum*rgb[2]); if ((lum > 0.1)&&(module2(x2 - x1, y2 - y1) < 0.25*(YMAX - YMIN))) draw_line(x1, y1, x2, y2); } // for (i=0; i < traj_position-1; i++) for (i = traj_position-2; i >= 0; i--) { x1 = trajectory[j*TRAJECTORY_LENGTH*TRACER_STEPS + i].xc; x2 = trajectory[j*TRAJECTORY_LENGTH*TRACER_STEPS + i+1].xc; y1 = trajectory[j*TRAJECTORY_LENGTH*TRACER_STEPS + i].yc; y2 = trajectory[j*TRAJECTORY_LENGTH*TRACER_STEPS + i+1].yc; time = traj_position - i; lum = 0.9 - TRACER_LUM_FACTOR*(double)time/(double)(TRAJECTORY_LENGTH*TRACER_STEPS); if (lum < 0.0) lum = 0.0; glColor3f(lum*rgb[0], lum*rgb[1], lum*rgb[2]); if ((lum > 0.1)&&(module2(x2 - x1, y2 - y1) < 0.25*(YMAX - YMIN))) draw_line(x1, y1, x2, y2); } } } } void color_background(t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NMAXOBSTACLES], int bg_color, t_hashgrid hashgrid[HASHX*HASHY]) /* color background according to particle properties */ { int i, j, k, n, p, q, m, nnb, number, avrg_fact, obs; double rgb[3], hue, value, p1, p2, pp1, pp2, oldhue, valx, valy, lum; static int first = 1, counter = 0; p1 = 0.75; p2 = 1.0 - p1; pp1 = 0.95; pp2 = 1.0 - pp1; glBegin(GL_QUADS); for (i=0; i ENERGY_HUE_MIN) hue = ENERGY_HUE_MIN; if (hue < ENERGY_HUE_MAX) hue = ENERGY_HUE_MAX; hue = p1*oldhue + p2*hue; hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_EKIN); break; } case (BG_EOBSTACLES): { value = 0.0; nnb = hashgrid[n].nneighb; for (q=0; q ENERGY_HUE_MIN) hue = ENERGY_HUE_MIN; if (hue < ENERGY_HUE_MAX) hue = ENERGY_HUE_MAX; hue = p1*oldhue + p2*hue; hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_EKIN); break; } case (BG_EKIN_OBSTACLES): { value = 0.0; nnb = hashgrid[n].nneighb; for (q=0; q ENERGY_HUE_MIN) hue = ENERGY_HUE_MIN; if (hue < ENERGY_HUE_MAX) hue = ENERGY_HUE_MAX; hue = p1*oldhue + p2*hue; hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_EKIN); break; } case (BG_DIR_OBSTACLES): { valx = 0.0; valy = 0.0; nnb = hashgrid[n].nneighb; for (q=0; q DPI) hue -= DPI; if (hue < 0.0) hue += DPI; hue = pp1*oldhue + pp2*hue; lum = module2(valx, valy)/OBSTACLE_VMAX; if (lum > 1.0) lum = 1.0; hsl_to_rgb_palette(hue*360.0/DPI, 0.9, 0.5*lum, rgb, COLOR_PALETTE_DIRECTION); break; } case (BG_POS_OBSTACLES): { valx = 0.0; valy = 0.0; nnb = hashgrid[n].nneighb; // for (q=0; q DPI) hue -= DPI; if (hue < 0.0) hue += DPI; hue = pp1*oldhue + pp2*hue; lum = module2(valx, valy); if (lum > 1.0) lum = 1.0; // lum = 1.0; hsl_to_rgb_palette(hue*360.0/DPI, 0.9, 0.5*lum, rgb, COLOR_PALETTE_DIRECTION); break; } case (BG_FORCE): { value = 0.0; nnb = hashgrid[n].nneighb; for (q=0; q 0)) color_background(particle, bg_color, hashgrid); // else if ((!TRACER_PARTICLE)&&(!COLOR_BACKGROUND)&&(!DRAW_OBSTACLE_LINKS)&&(!FILL_OBSTACLE_TRIANGLES)) blank(); // printf("Drawing particles 1\n"); glColor3f(1.0, 1.0, 1.0); /* show region of partial thermostat */ if (PARTIAL_THERMO_COUPLING) { switch (PARTIAL_THERMO_REGION) { case (TH_INBOX): { if (INCREASE_BETA) { logratio = log(beta/BETA)/log(0.5*BETA_FACTOR); if (logratio > 1.0) logratio = 1.0; else if (logratio < 0.0) logratio = 0.0; if (BETA_FACTOR > 1.0) hue = PARTICLE_HUE_MAX - (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*logratio; else hue = PARTICLE_HUE_MIN - (PARTICLE_HUE_MIN - PARTICLE_HUE_MAX)*logratio; } else hue = 0.25*PARTICLE_HUE_MIN + 0.75*PARTICLE_HUE_MAX; erase_area_hsl_turbo(0.0, YMIN, 2.0*PARTIAL_THERMO_WIDTH, PARTIAL_THERMO_HEIGHT*(YMAX - YMIN), hue, 0.9, 0.15); break; } case(TH_LAYER): { hue = 0.75*PARTICLE_HUE_MIN + 0.25*PARTICLE_HUE_MAX; erase_area_hsl_turbo(0.0, YMIN, XMAX, PARTIAL_THERMO_HEIGHT - YMIN, hue, 0.9, 0.15); break; } case (TH_RING): { hue = 0.75*PARTICLE_HUE_MIN + 0.25*PARTICLE_HUE_MAX; hsl_to_rgb_turbo(hue, 0.9, 0.15, rgb); draw_colored_circle(0.0, 0.0, PARTIAL_THERMO_WIDTH, 180, rgb); break; } case (TH_RING_EXPAND): { hue = 0.75*PARTICLE_HUE_MIN + 0.25*PARTICLE_HUE_MAX; hsl_to_rgb_turbo(hue, 0.9, 0.15, rgb); draw_colored_circle(0.0, 0.0, params.thermo_radius, 180, rgb); break; } case(TH_INIT): { hue = 0.75*PARTICLE_HUE_MIN + 0.25*PARTICLE_HUE_MAX; erase_rectangle_hsl_turbo(INITXMIN, INITYMIN, INITXMAX, INITYMAX, hue, 0.9, 0.15); break; } case(TH_THERMO): { if (INCREASE_BETA) { logratio = log(beta/BETA)/log(0.5*BETA_FACTOR); if (logratio > 1.0) logratio = 1.0; else if (logratio < 0.0) logratio = 0.0; if (BETA_FACTOR > 1.0) hue = PARTICLE_HUE_MAX - (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*logratio; else hue = PARTICLE_HUE_MIN - (PARTICLE_HUE_MIN - PARTICLE_HUE_MAX)*logratio; } else hue = 0.25*PARTICLE_HUE_MIN + 0.75*PARTICLE_HUE_MAX; erase_rectangle_hsl_turbo(THERMOXMIN, THERMOYMIN, THERMOXMAX, THERMOYMAX, hue, 0.9, 0.15); break; } case (TH_CONE): { if (INCREASE_BETA) { logratio = log(beta/BETA)/log(0.5*BETA_FACTOR); if (logratio > 1.0) logratio = 1.0; else if (logratio < 0.0) logratio = 0.0; if (BETA_FACTOR > 1.0) hue = PARTICLE_HUE_MAX - (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*logratio; else hue = PARTICLE_HUE_MIN - (PARTICLE_HUE_MIN - PARTICLE_HUE_MAX)*logratio; } else hue = 0.25*PARTICLE_HUE_MIN + 0.75*PARTICLE_HUE_MAX; draw_colored_triangle_turbo(0.05, -0.2, LAMBDA, LAMBDA, -LAMBDA, LAMBDA, hue, 0.9, 0.15); draw_colored_triangle_turbo(0.05, -0.2, -LAMBDA, LAMBDA, -0.04, -0.2, hue, 0.9, 0.15); break; } } } // printf("Drawing particles 2\n"); /* draw "altitude lines" */ if (ALTITUDE_LINES) draw_altitude_lines(); /* draw the bonds first */ if ((DRAW_BONDS)||(plot == P_BONDS)) { glLineWidth(LINK_WIDTH); for (j=0; j 0)) draw_collisions(collisions, ncollisions); /* determine particle color and size */ for (j=0; j= 0.0)) for (i=-1; i<2; i++) for (k=-1; k<2; k++) { x = x1 + (double)i*periodx; y = y1 + (double)k*periody; if (x < 1.2*particle[j].radius) draw_one_particle(particle[j], x, y, radius, angle, nsides, width, rgb); } else if ((x1 >= 0.0)&&(y1 < 0.0)) for (i=-1; i<2; i++) for (k=-1; k<2; k++) { x = x1 + (double)i*periodx; y = y1 + (double)k*periody; if (y < 1.2*particle[j].radius) draw_one_particle(particle[j], x, y, radius, angle, nsides, width, rgb); } } else if (plot != P_NOPARTICLE) draw_one_particle(particle[j], particle[j].xc, particle[j].yc, radius, angle, nsides, width, rgb); } // /* draw spin vectors */ if ((DRAW_SPIN)||(DRAW_SPIN_B)) { glLineWidth(width); for (j=0; j ENERGY_HUE_MIN) hue = ENERGY_HUE_MIN; if (hue < ENERGY_HUE_MAX) hue = ENERGY_HUE_MAX; hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_EKIN); break; } case (OC_DIRECTION): { hue = argument(obstacle.vx, obstacle.vy); if (hue < 0.0) hue += DPI; hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*(hue)/DPI; hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_DIRECTION); break; } case (OC_DIRECTION_ENERGY): { hue = argument(obstacle.vx, obstacle.vy); etot = obstacle.energy; if (hue < 0.0) hue += DPI; hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*(hue)/DPI; lum = etot/OBSTACLE_EMAX; if (lum > 0.5) lum = 0.5; hsl_to_rgb_palette(hue, 0.9, lum, rgb, COLOR_PALETTE_DIRECTION); break; } } } int area_to_rgb(t_obstacle obstacle[NMAXOBSTACLES], int i, int n1, int n2, double rgb[3]) /* computes rgb code of triangle according to area */ /* returns 1 if angle of triangle is less than PI/2 */ /* OLD VERSION */ { int p, q, k, n, counter; double area, angle1, angle2; static int first; static double mean_area; if (first) { /* compute mean area of triangles */ if (FILL_OBSTACLE_TRIANGLES) { counter = 0; mean_area = 0.0; for (p = 0; p < nobstacles; p++) { n = obstacle[p].nneighb; for (q = 0; q < n - 1; q++) { n1 = obstacle[p].neighb[q]; n2 = obstacle[p].neighb[q+1]; mean_area += vabs(triangle_area(obstacle[p].xc, obstacle[p].yc, obstacle[n1].xc, obstacle[n1].yc, obstacle[n2].xc, obstacle[n2].yc)); counter++; } n1 = obstacle[p].neighb[n-1]; n2 = obstacle[p].neighb[0]; mean_area += vabs(triangle_area(obstacle[p].xc, obstacle[p].yc, obstacle[n1].xc, obstacle[n1].yc, obstacle[n2].xc, obstacle[n2].yc)); counter++; } mean_area *= 1.0/(double)counter; } first = 0; } area = vabs(triangle_area(obstacle[i].xc, obstacle[i].yc, obstacle[n1].xc, obstacle[n1].yc, obstacle[n2].xc, obstacle[n2].yc)); area = OBSTACLE_AREA_SHADE_FACTOR*(area - mean_area); if (area > 1.0) area = 1.0; if (area < 0.0) area = 0.0; for (k=0; k<3; k++) rgb[k] = area; angle1 = argument(obstacle[n1].xc - obstacle[i].xc, obstacle[n1].yc - obstacle[i].yc); angle2 = argument(obstacle[n2].xc - obstacle[i].xc, obstacle[n2].yc - obstacle[i].yc); if (angle2 < angle1) angle2 += DPI; return(angle2 - angle1 < 1.5*PID); } void old_draw_obstacle_triangles(t_obstacle obstacle[NMAXOBSTACLES]) /* draw the triangles between obstacles, for option FILL_OBSTACLE_TRIANGLES */ { int i, j, k, n, n1, n2, neigh; short int acute; double rgb[3], angle1, angle2; if (FILL_OBSTACLE_TRIANGLES) for (i = 0; i < nobstacles; i++) { n = obstacle[i].nneighb; for (j = 0; j < n - 1; j++) { n1 = obstacle[i].neighb[j]; n2 = obstacle[i].neighb[j+1]; acute = area_to_rgb(obstacle, i, n1, n2, rgb); if (acute) draw_colored_triangle(obstacle[i].xc, obstacle[i].yc, obstacle[n1].xc, obstacle[n1].yc, obstacle[n2].xc, obstacle[n2].yc, rgb); } n1 = obstacle[i].neighb[n-1]; n2 = obstacle[i].neighb[0]; acute = area_to_rgb(obstacle, i, n1, n2, rgb); if (acute) draw_colored_triangle(obstacle[i].xc, obstacle[i].yc, obstacle[n1].xc, obstacle[n1].yc, obstacle[n2].xc, obstacle[n2].yc, rgb); } glColor3f(1.0, 1.0, 1.0); for (i = 0; i < nobstacles; i++) for (j = 0; j < obstacle[i].nneighb; j++) { neigh = obstacle[i].neighb[j]; draw_line(obstacle[i].xc, obstacle[i].yc, obstacle[neigh].xc, obstacle[neigh].yc); } } void facet_area_to_rgb(t_obstacle obstacle[NMAXOBSTACLES], t_otriangle triangle[NMAX_TRIANGLES_PER_OBSTACLE*NMAXOBSTACLES], t_ofacet facet[NMAXOBSTACLES], int nfacet, double rgb[3]) /* computes rgb code of triangle according to area */ /* returns 1 if angle of triangle is less than PI/2 */ { int p, q, k, n, counter; double area, area0; static int first; static double mean_area; if (first) { /* compute mean area of triangles */ if (FILL_OBSTACLE_TRIANGLES) { counter = 0; mean_area = 0.0; for (p = 0; p < n_otriangles; p++) { mean_area += otriangle_area(obstacle, triangle[p]); counter++; } mean_area *= 1.0/(double)counter; for (p = 0; p < n_ofacets; p++) { area0 = 0.0; for (k = 0; k < facet[p].ntriangles; k++) area0 += triangle[facet[p].triangle[k]].area0; facet[p].area0 = area0/(double)facet[p].ntriangles; } } first = 0; } area = 0.0; for (k = 0; k < facet[nfacet].ntriangles; k++) { p = facet[nfacet].triangle[k]; area += otriangle_area(obstacle, triangle[p]); } area *= 1.0/(double)facet[nfacet].ntriangles; area = OBSTACLE_AREA_SHADE_FACTOR*(area - facet[nfacet].area0); if (area > 1.0) area = 1.0; if (area < 0.0) area = 0.0; for (k=0; k<3; k++) rgb[k] = area; } double triangle_area_to_rgb(t_obstacle obstacle[NMAXOBSTACLES], t_otriangle otriangle[NMAX_TRIANGLES_PER_OBSTACLE*NMAXOBSTACLES], int triangle, double rgb[3]) { double area; int k, p; area = otriangle_area(obstacle, otriangle[triangle]); area = 0.6 + OBSTACLE_AREA_SHADE_FACTOR*(area - otriangle[triangle].area0); if (area > 0.9) area = 0.9; if (area < 0.1) area = 0.1; for (k=0; k<3; k++) rgb[k] = area; return(area); } void draw_obstacle_triangles(t_obstacle obstacle[NMAXOBSTACLES], t_otriangle otriangle[NMAX_TRIANGLES_PER_OBSTACLE*NMAXOBSTACLES], t_ofacet ofacet[NMAXOBSTACLES]) /* draw the triangles between obstacles, for option FILL_OBSTACLE_TRIANGLES */ { int i, j, k, n, n1, n2, neigh, nf, f, t, s1, s2, s3; short int acute; double rgb[3], pvect; printf("drawing %i facets with %i triangles\n", n_ofacets, n_otriangles); if (FILL_OBSTACLE_TRIANGLES) { /* determine acute angles */ for (t = 0; t < n_otriangles; t++) { s1 = otriangle[t].i[0]; s2 = otriangle[t].i[1]; s3 = otriangle[t].i[2]; pvect = (obstacle[s2].xc - obstacle[s1].xc)*(obstacle[s3].yc - obstacle[s1].yc); pvect -= (obstacle[s2].yc - obstacle[s1].yc)*(obstacle[s3].xc - obstacle[s1].xc); otriangle[t].acute = (pvect > 0.0); } if (SHADE_OBSTACLE_FACETS) for (i = 0; i < n_ofacets; i++) { facet_area_to_rgb(obstacle, otriangle, ofacet, i, rgb); nf = ofacet[i].ntriangles; for (f = 0; f < nf; f++) { t = ofacet[i].triangle[f]; if (otriangle[t].acute) { s1 = otriangle[t].i[0]; s2 = otriangle[t].i[1]; s3 = otriangle[t].i[2]; draw_colored_triangle(obstacle[s1].xc, obstacle[s1].yc, obstacle[s2].xc, obstacle[s2].yc, obstacle[s3].xc, obstacle[s3].yc, rgb); } } } else for (t = 0; t < n_otriangles; t++) if (otriangle[t].acute) { triangle_area_to_rgb(obstacle, otriangle, t, rgb); s1 = otriangle[t].i[0]; s2 = otriangle[t].i[1]; s3 = otriangle[t].i[2]; draw_colored_triangle(obstacle[s1].xc, obstacle[s1].yc, obstacle[s2].xc, obstacle[s2].yc, obstacle[s3].xc, obstacle[s3].yc, rgb); } } glColor3f(1.0, 1.0, 1.0); for (i = 0; i < nobstacles; i++) for (j = 0; j < obstacle[i].nneighb; j++) { neigh = obstacle[i].neighb[j]; draw_line(obstacle[i].xc, obstacle[i].yc, obstacle[neigh].xc, obstacle[neigh].yc); } } void draw_container(double xmin, double xmax, t_obstacle obstacle[NMAXOBSTACLES], t_segment segment[NMAXSEGMENTS], t_belt *conveyor_belt, int wall) /* draw the container, for certain boundary conditions */ { int i, j, k, n, n1, n2, neigh; double rgb[3], x, y, r, phi, angle, dx, dy, ybin, x1, x2, y1, y2, h, bangle, blength, pos0, pos, bsegment, tx, ty, bwidth, ca, sa, ekin, epot, etot, hue, area; char message[100]; /* draw fixed obstacles */ if (ADD_FIXED_OBSTACLES) { glLineWidth(CONTAINER_WIDTH); if (CHARGE_OBSTACLES) hsl_to_rgb(30.0, 0.1, 0.5, rgb); else if (!OSCILLATE_OBSTACLES) hsl_to_rgb(300.0, 0.1, 0.5, rgb); for (i = 0; i < nobstacles; i++) { if (OSCILLATE_OBSTACLES) { obstacle_hue(obstacle[i], rgb); } draw_colored_circle_precomp(obstacle[i].xc - xtrack, obstacle[i].yc - ytrack, obstacle[i].radius, rgb); } glColor3f(1.0, 1.0, 1.0); for (i = 0; i < nobstacles; i++) { x = obstacle[i].xc - xtrack; y = obstacle[i].yc - ytrack; r = obstacle[i].radius; draw_circle_precomp(x, y, r); if (CHARGE_OBSTACLES) { draw_line(x - r, y, x + r, y); draw_line(x, y - r, x, y + r); } if (ROTATE_OBSTACLES) { ca = cos(obstacle[i].angle); sa = sin(obstacle[i].angle); // printf("Obstacle %i angle = %.5lg\n", i, obstacle[i].angle); draw_line(x - r*ca, y - r*sa, x + r*ca, y + r*sa); draw_line(x + r*sa, y - r*ca, x - r*sa, y + r*ca); } } } if (ADD_FIXED_SEGMENTS) { glLineWidth(CONTAINER_WIDTH); glColor3f(1.0, 1.0, 1.0); for (i = 0; i < nsegments; i++) if (segment[i].active) { if (COLOR_SEG_GROUPS) set_segment_group_color(segment[i].group, 1.0, rgb); else if (SHOW_SEGMENTS_PRESSURE) set_segment_pressure_color(segment[i].avrg_pressure, 1.0, rgb); draw_line(segment[i].x1 - xtrack, segment[i].y1 - ytrack, segment[i].x2 - xtrack, segment[i].y2 - ytrack); } /* draw conveyor belt */ if (ADD_CONVEYOR_FORCE) for (i=0; i 0.0) pos0 = conveyor_belt[i].position - bsegment*(double)((int)(conveyor_belt[i].position/bsegment)); else { pos = -conveyor_belt[i].position; pos0 = pos - bsegment*(double)((int)(pos/bsegment)); pos0 = bsegment - pos0; } while (pos0 < conveyor_belt[i].length) { x = x1 + tx*pos0; y = y1 + ty*pos0; draw_line(x - r*ty, y + r*tx, x - bwidth*ty, y + bwidth*tx); pos0 += bsegment; } while (pos0 < conveyor_belt[i].length + PI*bwidth) { pos = (pos0 - conveyor_belt[i].length)/bwidth; angle = PID + conveyor_belt[i].angle - pos; ca = cos(angle); sa = sin(angle); draw_line(x2 + r*ca, y2 + r*sa, x2 + bwidth*ca, y2 + bwidth*sa); pos0 += bsegment; } while (pos0 < 2.0*conveyor_belt[i].length + PI*bwidth) { pos = pos0 - conveyor_belt[i].length - PI*bwidth; x = x2 - tx*pos; y = y2 - ty*pos; draw_line(x + r*ty, y - r*tx, x + bwidth*ty, y - bwidth*tx); pos0 += bsegment; } while (pos0 < 2.0*conveyor_belt[i].length + DPI*bwidth) { pos = (pos0 - 2.0*conveyor_belt[i].length - PI*bwidth)/bwidth; angle = -PID + conveyor_belt[i].angle - pos; ca = cos(angle); sa = sin(angle); draw_line(x1 + r*ca, y1 + r*sa, x1 + bwidth*ca, y1 + bwidth*sa); pos0 += bsegment; } } } switch (BOUNDARY_COND) { case (BC_SCREEN): { /* do nothing */ break; } case (BC_RECTANGLE): { glColor3f(1.0, 1.0, 1.0); glLineWidth(CONTAINER_WIDTH); draw_line(BCXMIN, BCYMIN, BCXMAX, BCYMIN); draw_line(BCXMIN, BCYMAX, BCXMAX, BCYMAX); if (!SYMMETRIC_DECREASE) draw_line(BCXMAX, BCYMIN, BCXMAX, BCYMAX); draw_line(xmin, BCYMIN, xmin, BCYMAX); // draw_line(XMIN, 0.5*(BCYMIN + BCYMAX), xmin, 0.5*(BCYMIN + BCYMAX)); if (SYMMETRIC_DECREASE) { draw_line(xmax, BCYMIN, xmax, BCYMAX); draw_line(XMAX, 0.5*(BCYMIN + BCYMAX), xmax, 0.5*(BCYMIN + BCYMAX)); } break; } case (BC_CIRCLE): { glLineWidth(CONTAINER_WIDTH); hsl_to_rgb(300.0, 0.1, 0.5, rgb); for (i=-1; i<2; i++) { if (CENTER_VIEW_ON_OBSTACLE) x = 0.0; else x = xmin + (double)i*(OBSXMAX - OBSXMIN); draw_colored_circle_precomp(x, 0.0, OBSTACLE_RADIUS, rgb); glColor3f(1.0, 1.0, 1.0); draw_circle_precomp(x, 0.0, OBSTACLE_RADIUS); glColor3f(0.0, 0.0, 0.0); sprintf(message, "Mach %.3f", xspeed/20.0); // sprintf(message, "Speed %.2f", xspeed); write_text(x-0.17, -0.025, message); } break; } case (BC_PERIODIC_CIRCLE): { glLineWidth(CONTAINER_WIDTH); hsl_to_rgb(300.0, 0.1, 0.5, rgb); for (i=-1; i<2; i++) { if (CENTER_VIEW_ON_OBSTACLE) x = 0.0; else x = xmin + (double)i*(OBSXMAX - OBSXMIN); draw_colored_circle_precomp(x, 0.0, OBSTACLE_RADIUS, rgb); glColor3f(1.0, 1.0, 1.0); draw_circle_precomp(x, 0.0, OBSTACLE_RADIUS); glColor3f(0.0, 0.0, 0.0); sprintf(message, "Mach %.2f", xspeed/20.0); // sprintf(message, "Speed %.2f", xspeed); write_text(x-0.17, -0.025, message); } break; } case (BC_PERIODIC_TRIANGLE): { glLineWidth(CONTAINER_WIDTH); hsl_to_rgb(300.0, 0.1, 0.5, rgb); for (i=-1; i<2; i++) { if (CENTER_VIEW_ON_OBSTACLE) x = 0.0; else x = xmin + (double)i*(OBSXMAX - OBSXMIN); x1 = x + OBSTACLE_RADIUS; x2 = x - OBSTACLE_RADIUS; h = 2.0*OBSTACLE_RADIUS*tan(APOLY*PID); draw_colored_triangle(x1, 0.0, x2, h, x2, -h, rgb); glColor3f(1.0, 1.0, 1.0); draw_triangle(x1, 0.0, x2, h, x2, -h); glColor3f(0.0, 0.0, 0.0); sprintf(message, "Mach %.2f", xspeed*3.0/40.0); write_text(x-0.25, -0.025, message); } break; } case (BC_PERIODIC_FUNNEL): { glLineWidth(CONTAINER_WIDTH); hsl_to_rgb(300.0, 0.1, 0.5, rgb); for (i=-1; i<2; i++) { if (CENTER_VIEW_ON_OBSTACLE) x = 0.0; else x = xmin + (double)i*(OBSXMAX - OBSXMIN); for (j=-1; j<2; j+=2) { draw_colored_circle_precomp(x, (double)j*(FUNNEL_WIDTH + OBSTACLE_RADIUS), OBSTACLE_RADIUS, rgb); glColor3f(1.0, 1.0, 1.0); draw_circle_precomp(x, (double)j*(FUNNEL_WIDTH + OBSTACLE_RADIUS), OBSTACLE_RADIUS); } glColor3f(0.0, 0.0, 0.0); sprintf(message, "Mach %.2f", xspeed/20.0); write_text(x-0.17, 0.75, message); } break; } case (BC_RECTANGLE_LID): { glColor3f(1.0, 1.0, 1.0); glLineWidth(CONTAINER_WIDTH); draw_line(BCXMIN, BCYMIN, BCXMAX, BCYMIN); draw_line(BCXMIN, BCYMIN, BCXMIN, BCYMAX); draw_line(BCXMAX, BCYMIN, BCXMAX, BCYMAX); hsl_to_rgb(300.0, 0.1, 0.5, rgb); draw_colored_rectangle(BCXMIN + 0.05, ylid, BCXMAX - 0.05, ylid + LID_WIDTH, rgb); glColor3f(1.0, 1.0, 1.0); draw_rectangle(BCXMIN + 0.05, ylid, BCXMAX - 0.05, ylid + LID_WIDTH); break; } case (BC_RECTANGLE_WALL): { glColor3f(1.0, 1.0, 1.0); glLineWidth(CONTAINER_WIDTH); draw_rectangle(BCXMIN, BCYMIN, BCXMAX, BCYMAX); draw_line(0.5*(BCXMIN+BCXMAX), BCYMAX, 0.5*(BCXMIN+BCXMAX), BCYMAX + 0.5*WALL_WIDTH); draw_line(0.5*(BCXMIN+BCXMAX), BCYMIN, 0.5*(BCXMIN+BCXMAX), BCYMIN - 0.5*WALL_WIDTH); if (wall) { hsl_to_rgb(300.0, 0.1, 0.5, rgb); draw_colored_rectangle(xwall - 0.5*WALL_WIDTH, BCYMIN + 0.025, xwall + 0.5*WALL_WIDTH, BCYMAX - 0.025, rgb); glColor3f(1.0, 1.0, 1.0); draw_rectangle(xwall - 0.5*WALL_WIDTH, BCYMIN + 0.025, xwall + 0.5*WALL_WIDTH, BCYMAX - 0.025); } break; } case (BC_EHRENFEST): { glLineWidth(CONTAINER_WIDTH); glColor3f(1.0, 1.0, 1.0); phi = asin(EHRENFEST_WIDTH/EHRENFEST_RADIUS); glBegin(GL_LINE_LOOP); for (i=0; i<=NSEG; i++) { angle = -PI + phi + (double)i*2.0*(PI - phi)/(double)NSEG; glVertex2d(1.0 + EHRENFEST_RADIUS*cos(angle), EHRENFEST_RADIUS*sin(angle)); } for (i=0; i<=NSEG; i++) { angle = phi + (double)i*2.0*(PI - phi)/(double)NSEG; glVertex2d(-1.0 + EHRENFEST_RADIUS*cos(angle), EHRENFEST_RADIUS*sin(angle)); } glEnd(); break; } case (BC_SCREEN_BINS): { glLineWidth(CONTAINER_WIDTH); glColor3f(1.0, 1.0, 1.0); dy = (YMAX - YMIN)/((double)NGRIDX + 3); dx = dy/cos(PI/6.0); ybin = 2.75*dy; for (i=-1; i<=NGRIDX; i++) { x = ((double)i - 0.5*(double)NGRIDX + 0.5)*dx; draw_line(x, YMIN, x, YMIN + ybin); } break; } case (BC_GENUS_TWO): { hsl_to_rgb(300.0, 0.1, 0.5, rgb); draw_colored_rectangle(0.0, 0.0, BCXMAX, BCYMAX, rgb); break; } default: { /* do nothing */ } } } void print_omega(double angle, double angular_speed, double fx, double fy) { char message[100]; double rgb[3], y1, frac, absa; static double xleftbox, xlefttext, xrightbox, xrighttext, y = YMAX - 0.1, ymin = YMIN + 0.05; static int first = 1; if (first) { xrightbox = XMAX - 0.54; xrighttext = xrightbox - 0.48; first = 0; } y1 = y; if (PRINT_ANGLE) { erase_area_hsl(xrightbox, y + 0.025, 0.42, 0.05, 0.0, 0.9, 0.0); glColor3f(1.0, 1.0, 1.0); angle = angle*360.0/DPI; absa = vabs(angle); frac = absa - (double)((int)absa); sprintf(message, "Angle = %3d.%d degrees", (int)absa, (int)(10.0*frac)); write_text(xrighttext + 0.1, y, message); y1 -= 0.1; } if (PRINT_OMEGA) { erase_area_hsl(xrightbox, y1 + 0.025, 0.42, 0.05, 0.0, 0.9, 0.0); glColor3f(1.0, 1.0, 1.0); sprintf(message, "Angular speed = %.4f", angular_speed); write_text(xrighttext + 0.1, y1, message); y1 -= 0.1; } if (PRINT_SEGMENTS_FORCE) { erase_area_hsl(xrightbox, y1 + 0.025, 0.42, 0.05, 0.0, 0.9, 0.0); glColor3f(1.0, 1.0, 1.0); sprintf(message, "Fx = %.4f", fx); write_text(xrighttext + 0.1, y1, message); y1 -= 0.1; erase_area_hsl(xrightbox, y1 + 0.025, 0.42, 0.05, 0.0, 0.9, 0.0); glColor3f(1.0, 1.0, 1.0); sprintf(message, "Fy = %.4f", fy); write_text(xrighttext + 0.1, y1, message); y1 -= 0.1; erase_area_hsl(xrightbox, y1 + 0.025, 0.42, 0.05, 0.0, 0.9, 0.0); glColor3f(1.0, 1.0, 1.0); sprintf(message, "Fy/Fx = %.4f", fy/fx); write_text(xrighttext + 0.1, y1, message); } } void compute_segments_force(t_lj_parameters *params, t_segment segment[NMAXSEGMENTS]) { int i; double fx = 0.0, fy = 0.0; for (i=0; ibdry_fx = fx; params->bdry_fy = fy; } void print_parameters(t_lj_parameters params, short int left, double pressure[N_PRESSURES], short int refresh) { char message[100]; int i, j, k; double density, hue, rgb[3], logratio, x, y, meanpress[N_PRESSURES], phi, sphi, dphi, pprint, mean_temp, lengthcontainer, boundary_force, fx, fy, r1, r2; static double xbox, xtext, xmid, xmidtext, xxbox, xxtext, pressures[N_P_AVERAGE], meanpressure = 0.0, maxpressure = 0.0, mean_fx, mean_fy; static double press[N_PRESSURES][N_P_AVERAGE], temp[N_T_AVERAGE], scale; static int first = 1, i_pressure, i_temp; if (first) { // scale = (XMAX - XMIN)/4.0; scale = (YMAX - YMIN)/2.5; if (left) { xbox = XMIN + 0.45*scale; xtext = XMIN + 0.12*scale; xxbox = XMAX - 0.39*scale; xxtext = XMAX - 0.73*scale; } else { xbox = XMAX - 0.41*scale; xtext = XMAX - 0.73*scale; xxbox = XMIN + 0.4*scale; xxtext = XMIN + 0.08*scale; } xmid = 0.5*(XMIN + XMAX) - 0.1*scale; // xmid = 0.5*(XMIN + XMAX) + 0.05*scale; xmidtext = xmid - 0.2*scale; for (i=0; i maxpressure) maxpressure = meanpress[j]; printf("Max pressure = %.5lg\n\n", maxpressure); } y = YMAX - 0.1*scale; if ((INCREASE_BETA)||(PRINT_TEMPERATURE)) /* print temperature */ { logratio = log(params.beta/BETA)/log(0.5*BETA_FACTOR); if (logratio > 1.0) logratio = 1.0; else if (logratio < 0.0) logratio = 0.0; if (BETA_FACTOR > 1.0) hue = PARTICLE_HUE_MAX - (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*logratio; else hue = PARTICLE_HUE_MIN - (PARTICLE_HUE_MIN - PARTICLE_HUE_MAX)*logratio; if (PRINT_LEFT) erase_area_hsl_turbo(xbox, y + 0.025*scale, 0.5*scale, 0.05*scale, hue, 0.9, 0.5); else erase_area_hsl_turbo(xmid + 0.1, y + 0.025*scale, 0.45*scale, 0.05*scale, hue, 0.9, 0.5); if ((hue < 90)||(hue > 270)) glColor3f(1.0, 1.0, 1.0); else glColor3f(0.0, 0.0, 0.0); sprintf(message, "Temperature %.2f", 1.0/params.beta); if (PRINT_LEFT) write_text(xtext, y, message); else write_text(xmidtext, y, message); // y -= 0.1; // erase_area_hsl(xxbox, y + 0.025, 0.37, 0.05, 0.0, 0.9, 0.0); // glColor3f(1.0, 1.0, 1.0); // sprintf(message, "Pressure %.3f", meanpressure); // write_text(xxtext, y, message); } if (DECREASE_CONTAINER_SIZE) /* print density */ { density = (double)ncircles/((lengthcontainer)*(INITYMAX - INITYMIN)); erase_area_hsl(xbox, y + 0.025*scale, 0.37*scale, 0.05*scale, 0.0, 0.9, 0.0); glColor3f(1.0, 1.0, 1.0); sprintf(message, "Density %.3f", density); write_text(xtext + 0.1, y, message); erase_area_hsl(xmid, y + 0.025*scale, 0.37*scale, 0.05*scale, 0.0, 0.9, 0.0); glColor3f(1.0, 1.0, 1.0); sprintf(message, "Temperature %.2f", params.mean_energy); write_text(xmidtext - 0.1, y, message); erase_area_hsl(xxbox, y + 0.025*scale, 0.37*scale, 0.05*scale, 0.0, 0.9, 0.0); glColor3f(1.0, 1.0, 1.0); sprintf(message, "Pressure %.3f", meanpressure); write_text(xxtext, y, message); } else if (INCREASE_KREPEL) /* print force constant */ { erase_area_hsl(xbox, y + 0.025*scale, 0.35*scale, 0.05*scale, 0.0, 0.9, 0.0); glColor3f(1.0, 1.0, 1.0); sprintf(message, "Force constant %.2f", params.krepel); write_text(xtext + 0.03, y, message); } else if (INCREASE_E) /* print electric field */ { erase_area_hsl(xbox, y + 0.025*scale, 0.27*scale, 0.05*scale, 0.0, 0.9, 0.0); glColor3f(1.0, 1.0, 1.0); sprintf(message, "E field = %.2f", 25.0*NVID*DT_PARTICLE*params.efield); write_text(xtext + 0.08, y, message); } else if (INCREASE_B) /* print magnetic field */ { erase_area_hsl(xbox, y + 0.025*scale, 0.27*scale, 0.05*scale, 0.0, 0.9, 0.0); glColor3f(1.0, 1.0, 1.0); sprintf(message, "B field = %.2f", 25.0*NVID*DT_PARTICLE*params.bfield); write_text(xtext + 0.08, y, message); } if (PRINT_NPARTICLES) /* print number of particles */ { erase_area_hsl(xbox, y + 0.025*scale, 0.27*scale, 0.05*scale, 0.0, 0.9, 0.0); glColor3f(1.0, 1.0, 1.0); sprintf(message, "%i electrons", params.nactive); write_text(xtext + 0.08, y, message); } if (PRINT_TYPE_PROP) /* print proportion of types */ { erase_area_hsl(xbox, y + 0.025*scale, 0.27*scale, 0.05*scale, 0.0, 0.9, 0.0); glColor3f(1.0, 1.0, 1.0); sprintf(message, "%.1f %% anions", 100.0*params.prop); write_text(xtext + 0.08, y, message); } if (RECORD_PRESSURES) { y = FUNNEL_WIDTH + OBSTACLE_RADIUS; for (i=0; i PARTICLE_HUE_MIN) hue = PARTICLE_HUE_MIN; if (hue < PARTICLE_HUE_MAX) hue = PARTICLE_HUE_MAX; hsl_to_rgb_turbo(hue, 0.9, 0.5, rgb); dphi = DPI/(double)N_PRESSURES; // x = 0.95*OBSTACLE_RADIUS*cos(phi); sphi = sin(phi); if (sphi < 0.0) draw_colored_sector(0.0, y, 0.95*OBSTACLE_RADIUS, OBSTACLE_RADIUS, phi, phi + dphi, rgb, 10); else draw_colored_sector(0.0, -y, 0.95*OBSTACLE_RADIUS, OBSTACLE_RADIUS, phi, phi + dphi, rgb, 10); } glColor3f(1.0, 1.0, 1.0); for (i=-1; i<2; i++) { k = N_PRESSURES/4 + i*N_PRESSURES/9; phi = DPI*(double)k/(double)N_PRESSURES; pprint = 0.0; for (j=-2; j<3; j++) pprint += meanpress[k + j]; sprintf(message, "p = %.0f", pprint*200.0/MAX_PRESSURE); write_text(0.85*OBSTACLE_RADIUS*cos(phi) - 0.1, -y + 0.85*OBSTACLE_RADIUS*sin(phi), message); } } if ((PARTIAL_THERMO_COUPLING)&&(!INCREASE_BETA)&&(!EXOTHERMIC)) { printf("Temperature %i in average: %.3lg\n", i_temp, params.mean_energy); temp[i_temp] = params.mean_energy; i_temp++; if (i_temp >= N_T_AVERAGE) i_temp = 0; mean_temp = 0.0; for (i=0; i 270)) glColor3f(1.0, 1.0, 1.0); else glColor3f(0.0, 0.0, 0.0); sprintf(message, "Temperature %.2f", mean_temp); write_text(xtext, y, message); } if (INCREASE_GRAVITY) { erase_area_hsl(xmid, y + 0.025*scale, 0.22*scale, 0.05*scale, 0.0, 0.9, 0.0); glColor3f(1.0, 1.0, 1.0); sprintf(message, "Gravity %.2f", params.gravity/GRAVITY); // write_text(xmidtext + 0.1, y, message); write_text(xmidtext, y, message); } if (CHANGE_RADIUS) { erase_area_hsl(xmid, y + 0.025*scale, 0.3*scale, 0.05*scale, 0.0, 0.9, 0.0); glColor3f(1.0, 1.0, 1.0); sprintf(message, "Radius %.4f", params.radius); write_text(xmidtext + 0.05, y, message); } if (PRINT_SEGMENTS_FORCE) { glColor3f(1.0, 1.0, 1.0); if (refresh) { fx = 0.01*params.bdry_fx/(double)params.nactive; fy = 0.01*params.bdry_fy/(double)params.nactive; /* average boundary force */ mean_fx = r2*mean_fx + r1*fx; mean_fy = r2*mean_fy + r1*fy; } if ((PRINT_ANGLE)||(PRINT_OMEGA)) draw_arrow(0.0, 0.0, FORCE_FACTOR*mean_fx, FORCE_FACTOR*mean_fy, 15.0, 0.1); else { erase_area_hsl(xmid, 0.0, 0.2*scale, 0.05*scale, 0.0, 0.9, 0.0); glColor3f(1.0, 1.0, 1.0); sprintf(message, "Fy = %.2f", mean_fy); write_text(xmidtext + 0.15*scale, -0.02*scale, message); if (mean_fx*mean_fx + mean_fy*mean_fy > 5.0*FORCE_FACTOR) draw_arrow(0.0, 0.0, FORCE_FACTOR*mean_fx, FORCE_FACTOR*mean_fy, 15.0, 0.1); } } if ((PRINT_ANGLE)||(PRINT_OMEGA)) print_omega(params.angle, params.omega, mean_fx, mean_fy); } void print_ehrenfest_parameters(t_particle particle[NMAXCIRCLES], double pleft, double pright) { char message[100]; int i, j, nleft1 = 0, nleft2 = 0, nright1 = 0, nright2 = 0; double density, hue, rgb[3], logratio, y, shiftx = 0.3, xmidplus, xmidminus; static double xleftbox, xlefttext, xmidbox, xmidtext, xrightbox, xrighttext, pressures[500][2], meanpressure[2]; static int first = 1, i_pressure, naverage = 500, n_pressure; if (first) { xleftbox = -0.85; xlefttext = xleftbox - 0.5; xrightbox = 1.0; xrighttext = xrightbox - 0.45; // xmid = 0.5*(XMIN + XMAX) - 0.1; // xmidtext = xmid - 0.24; meanpressure[0] = 0.0; meanpressure[1] = 0.0; for (i=0; i xmidplus) { if (particle[i].type == 0) nright1++; else nright2++; } } y = YMIN + 0.05; erase_area_hsl(xleftbox - shiftx, y + 0.025, 0.22, 0.05, 0.0, 0.9, 0.0); hsl_to_rgb(HUE_TYPE0, 0.9, 0.5, rgb); glColor3f(rgb[0], rgb[1], rgb[2]); sprintf(message, "%i particles", nleft1); write_text(xlefttext + 0.28 - shiftx, y, message); erase_area_hsl(xleftbox + shiftx, y + 0.025, 0.22, 0.05, 0.0, 0.9, 0.0); hsl_to_rgb(HUE_TYPE1, 0.9, 0.5, rgb); glColor3f(rgb[0], rgb[1], rgb[2]); sprintf(message, "%i particles", nleft2); write_text(xlefttext + 0.28 + shiftx, y, message); erase_area_hsl(xrightbox - shiftx, y + 0.025, 0.22, 0.05, 0.0, 0.9, 0.0); hsl_to_rgb(HUE_TYPE0, 0.9, 0.5, rgb); glColor3f(rgb[0], rgb[1], rgb[2]); sprintf(message, "%i particles", nright1); write_text(xrighttext + 0.28 - shiftx, y, message); erase_area_hsl(xrightbox + shiftx, y + 0.025, 0.22, 0.05, 0.0, 0.9, 0.0); hsl_to_rgb(HUE_TYPE1, 0.9, 0.5, rgb); glColor3f(rgb[0], rgb[1], rgb[2]); sprintf(message, "%i particles", nright2); write_text(xrighttext + 0.28 + shiftx, y, message); y = YMAX - 0.1; erase_area_hsl(xleftbox - 0.1, y + 0.025, 0.22, 0.05, 0.0, 0.9, 0.0); hsl_to_rgb_turbo(HUE_TYPE1, 0.9, 0.5, rgb); glColor3f(rgb[0], rgb[1], rgb[2]); sprintf(message, "Pressure %.2f", 0.001*meanpressure[0]/(double)ncircles); write_text(xlefttext + 0.25, y, message); erase_area_hsl(xrightbox - 0.1, y + 0.025, 0.22, 0.05, 0.0, 0.9, 0.0); hsl_to_rgb_turbo(HUE_TYPE0, 0.9, 0.5, rgb); glColor3f(rgb[0], rgb[1], rgb[2]); sprintf(message, "Pressure %.2f", 0.001*meanpressure[1]/(double)ncircles); write_text(xrighttext + 0.2, y, message); } void count_particle_number(t_particle *particle, int *particle_numbers, int time) { int type, i, n; n = time*(RD_TYPES+1); for (type = 0; type < RD_TYPES+1; type++) particle_numbers[n + type] = 0; if (COUNT_PARTNER_TYPE) { for (i=0; i 1) sprintf(message, "%i particles", npart); else sprintf(message, "%i particle", npart); write_text(xlefttext + 0.28, y, message); } void print_particle_types_number(t_particle *particle, int ntypes) { int i, ntype[10]; char message[100]; double y = YMAX - 0.1, rgb[3]; static double xleftbox, xlefttext; static int first = 1; if (first) { xleftbox = XMAX - 0.5; xlefttext = xleftbox - 0.45; first = 0; } for (i=0; i<10; i++) ntype[i] = 0; for (i=0; i=t_window) position[group] = 0; for (i=0; i 2) { xbox = xleftbox + (xrightbox - xleftbox)*(group-1)/(ngroups-2); xtext = xlefttext + (xrighttext - xlefttext)*(group-1)/(ngroups-2); } else { xbox = xrightbox - 0.2; xtext = xrighttext - 0.2; } // printf("xbox = %.2f, xtext = %.2f, av_vy = %.2f\n", xbox, xtext, av_vy); y = YMAX - 0.1*scale; erase_area_hsl(xbox, y + 0.025*scale, 0.25*scale, 0.05*scale, 0.0, 0.9, 0.0); set_segment_group_color(group, 1.0, rgb); sprintf(message, "Vx = %.4f", av_vx); write_text(xtext + 0.1, y, message); y -= 0.1*scale; erase_area_hsl(xbox, y + 0.025*scale, 0.25*scale, 0.05*scale, 0.0, 0.9, 0.0); set_segment_group_color(group, 1.0, rgb); sprintf(message, "Vy = %.4f", av_vy); write_text(xtext + 0.1, y, message); y -= 0.1*scale; erase_area_hsl(xbox, y + 0.025*scale, 0.3*scale, 0.05*scale, 0.0, 0.9, 0.0); set_segment_group_color(group, 1.0, rgb); sprintf(message, "V = %.4f", module2(av_vx, av_vy)); write_text(xtext + 0.1, y, message); y -= 0.1*scale; erase_area_hsl(xbox, y + 0.025*scale, 0.3*scale, 0.05*scale, 0.0, 0.9, 0.0); set_segment_group_color(group, 1.0, rgb); sprintf(message, "Omega = %.4f", av_omega); write_text(xtext + 0.1, y, message); } } void print_particles_speeds(t_particle particle[NMAXCIRCLES]) { char message[100]; double y = YMAX - 0.1, vx = 0.0, vy = 0.0; int i, nactive = 0; static double xleftbox, xlefttext, xrightbox, xrighttext, ymin = YMIN + 0.05; static int first = 1; if (first) { xrightbox = XMAX - 0.39; xrighttext = xrightbox - 0.55; first = 0; } for (i=0; i 0.0)&&(proj < 1.0)) { // distance = segment[i].nx*x + segment[i].ny*y - segment[i].c; distance = segment[i].nx*x + segment[i].ny*y - segment[i].c; /* case of interacting polygons */ // if (POLYGON_INTERACTION) // { // // distance = segment[i].nx*x + segment[i].ny*y - segment[i].c; // pangle = DPI/(double)NPOLY; // dx = 0.0; // dy = 0.0; // for (corner=0; corner 0.0)&&(pvect < pvmin)) pvect = pvmin; else if ((pvect < 0.0)&&(pvect > -pvmin)) pvect = -pvmin; // particle[j].torque += KTORQUE_BOUNDARY*pvect; particle[j].torque += KTORQUE_BOUNDARY*pvect*(SEGMENT_FORCE_EQR*r - vabs(distance)); } } /* compute force from concave corners */ if (segment[i].concave) { distance = module2(x - segment[i].x1, y - segment[i].y1); angle = argument(x - segment[i].x1, y - segment[i].y1); if (angle < segment[i].angle1) angle += DPI; /* added 24/9/22 */ else if (angle > segment[i].angle2) angle -= DPI; r = SEGMENT_FORCE_EQR*particle[j].radius; if ((distance < r)&&(angle > segment[i].angle1)&&(angle < segment[i].angle2)) { f = KSPRING_OBSTACLE*(r - distance); particle[j].fx += f*cos(angle); particle[j].fy += f*sin(angle); if ((MOVE_BOUNDARY)||(MOVE_SEGMENT_GROUPS)) { segment[i].fx -= f*cos(angle); segment[i].fy -= f*sin(angle); segment[i].torque -= (x - segment[i].xc)*f*sin(angle) - (y - segment[i].yc)*f*cos(angle); } } } } if (INACTIVATE_SEGMENTS_UNDER_PRESSURE) { ngroups = 0; for (group=0; group ngroups) ngroups = group; } for (group=1; group 0) { // group_pressure[group] *= 1.0/(double)segments_per_group[group]; if (group_pressure[group] > SEGMENT_P_INACTIVATE*(double)segments_per_group[group]) inactivate_group[group] = 1; else inactivate_group[group] = 0; } for (i=0; i 0)&&(inactivate_group[group])) segment[i].active = 0; } } switch(BOUNDARY_COND){ case (BC_SCREEN): { /* add harmonic force outside screen */ if (particle[j].xc > XMAX) particle[j].fx -= KSPRING_BOUNDARY*(particle[j].xc - XMAX); else if (particle[j].xc < XMIN) particle[j].fx += KSPRING_BOUNDARY*(XMIN - particle[j].xc); if (particle[j].yc > YMAX) particle[j].fy -= KSPRING_BOUNDARY*(particle[j].yc - YMAX); else if (particle[j].yc < YMIN) particle[j].fy += KSPRING_BOUNDARY*(YMIN - particle[j].yc); // if (particle[j].xc > BCXMAX) particle[j].fx -= KSPRING_BOUNDARY*(particle[j].xc - BCXMAX); // else if (particle[j].xc < BCXMIN) particle[j].fx += KSPRING_BOUNDARY*(BCXMIN - particle[j].xc); // if (particle[j].yc > BCYMAX) particle[j].fy -= KSPRING_BOUNDARY*(particle[j].yc - BCYMAX); // else if (particle[j].yc < BCYMIN) particle[j].fy += KSPRING_BOUNDARY*(BCYMIN - particle[j].yc); return(fperp); } case (BC_RECTANGLE): { /* add harmonic force outside rectangular box */ padding = particle[j].radius + 0.01; xmin = xleft + padding; xmax = xright - padding; ymin = BCYMIN + padding; ymax = BCYMAX - padding; if (particle[j].xc > xmax) { fperp = KSPRING_BOUNDARY*(particle[j].xc - xmax); particle[j].fx -= fperp; } else if (particle[j].xc < xmin) { fperp = KSPRING_BOUNDARY*(xmin - particle[j].xc); particle[j].fx += fperp; } if (particle[j].yc > ymax) { fperp = KSPRING_BOUNDARY*(particle[j].yc - ymax); particle[j].fy -= fperp; } else if (particle[j].yc < ymin) { fperp = KSPRING_BOUNDARY*(ymin - particle[j].yc); particle[j].fy += fperp; } // if (particle[j].xc > xmax) particle[j].fx -= KSPRING_BOUNDARY*(particle[j].xc - xmax); // else if (particle[j].xc < xmin) particle[j].fx += KSPRING_BOUNDARY*(xmin - particle[j].xc); // if (particle[j].yc > ymax) particle[j].fy -= KSPRING_BOUNDARY*(particle[j].yc - ymax); // else if (particle[j].yc < ymin) particle[j].fy += KSPRING_BOUNDARY*(ymin - particle[j].yc); return(fperp); } case (BC_CIRCLE): { /* add harmonic force outside screen */ if (particle[j].xc > BCXMAX) particle[j].fx -= KSPRING_BOUNDARY*(particle[j].xc - BCXMAX); else if (particle[j].xc < BCXMIN) particle[j].fx += KSPRING_BOUNDARY*(BCXMIN - particle[j].xc); if (particle[j].yc > BCYMAX) particle[j].fy -= KSPRING_BOUNDARY*(particle[j].yc - BCYMAX); else if (particle[j].yc < BCYMIN) particle[j].fy += KSPRING_BOUNDARY*(BCYMIN - particle[j].yc); /* add harmonic force from obstacle */ for (i=-2; i<2; i++) { x = xleft + (double)i*(OBSXMAX - OBSXMIN); if (vabs(particle[j].xc - x) < 1.1*OBSTACLE_RADIUS) { r = module2(particle[j].xc - x, particle[j].yc); if (r < 1.0e-5) r = 1.0e-05; cphi = (particle[j].xc - x)/r; sphi = particle[j].yc/r; padding = MU + 0.03; if (r < OBSTACLE_RADIUS + padding) { f = KSPRING_OBSTACLE*(OBSTACLE_RADIUS + padding - r); particle[j].fx += f*cphi; particle[j].fy += f*sphi; } } } return(fperp); } case (BC_PERIODIC_CIRCLE): { x = xleft; if (vabs(particle[j].xc - x) < 1.1*OBSTACLE_RADIUS) { r = module2(particle[j].xc - x, particle[j].yc); if (r < 1.0e-5) r = 1.0e-05; cphi = (particle[j].xc - x)/r; sphi = particle[j].yc/r; padding = MU + 0.03; if (r < OBSTACLE_RADIUS + padding) { f = KSPRING_OBSTACLE*(OBSTACLE_RADIUS + padding - r); particle[j].fx += f*cphi; particle[j].fy += f*sphi; } } return(f); } case (BC_PERIODIC_TRIANGLE): { x = xleft; x1 = x + OBSTACLE_RADIUS; x2 = x - OBSTACLE_RADIUS; h = 2.0*OBSTACLE_RADIUS*tanh(APOLY*PID); padding = MU + 0.03; // ytop = 0.5*h*(1.0 - (particle[j].xc - x)/OBSTACLE_RADIUS); if ((vabs(particle[j].xc - x) < OBSTACLE_RADIUS + padding)&&(vabs(particle[j].yc < h + padding))) { /* signed distances to side of triangle */ dleft = x2 - particle[j].xc; norm = module2(h, 2.0*OBSTACLE_RADIUS); if (particle[j].yc >= 0.0) { dplus = (h*particle[j].xc + 2.0*OBSTACLE_RADIUS*particle[j].yc - h*(x+OBSTACLE_RADIUS))/norm; if ((dleft < padding)&&(dleft > dplus)) /* left side is closer */ { f = KSPRING_OBSTACLE*(padding - dleft); particle[j].fx -= f; } else if (dplus < padding) /* top side is closer */ { f = KSPRING_OBSTACLE*(padding - dplus); particle[j].fx += f*h/norm; particle[j].fy += 2.0*f*OBSTACLE_RADIUS/norm; } } else { dminus = (h*particle[j].xc - 2.0*OBSTACLE_RADIUS*particle[j].yc - h*(x+OBSTACLE_RADIUS))/norm; if ((dleft < padding)&&(dleft > dminus)) /* left side is closer */ { f = KSPRING_OBSTACLE*(padding - dleft); particle[j].fx -= f; } else if (dminus < padding) /* bottom side is closer */ { f = KSPRING_OBSTACLE*(padding - dminus); particle[j].fx += f*h/norm; particle[j].fy += -2.0*f*OBSTACLE_RADIUS/norm; } } /* force from tip of triangle */ r = module2(particle[j].xc - x1, particle[j].yc); if (r < 0.5*padding) { if (r < 1.0e-5) r = 1.0e-05; cphi = (particle[j].xc - x1)/r; sphi = particle[j].yc/r; f = KSPRING_OBSTACLE*(0.5*padding - r); particle[j].fx += f*cphi; particle[j].fy += f*sphi; } } return(f); } case (BC_PERIODIC_FUNNEL): { x = xleft; padding = MU + 0.02; if (vabs(particle[j].yc) > FUNNEL_WIDTH - padding) for (i=-1; i<2; i+=2) { r = module2(particle[j].xc - x, particle[j].yc - (double)i*(FUNNEL_WIDTH + OBSTACLE_RADIUS)); if (r < 1.0e-5) r = 1.0e-05; cphi = (particle[j].xc - x)/r; sphi = (particle[j].yc - (double)i*(FUNNEL_WIDTH + OBSTACLE_RADIUS))/r; if (r < OBSTACLE_RADIUS + padding) { f = KSPRING_OBSTACLE*(OBSTACLE_RADIUS + padding - r); particle[j].fx += f*cphi; particle[j].fy += f*sphi; if (RECORD_PRESSURES) { angle = argument(cphi, sphi); if (angle < 0.0) angle += DPI; k = (int)((double)N_PRESSURES*angle/DPI); if (k >= N_PRESSURES) k = N_PRESSURES - 1; pressure[k] += f; } } } return(f); } case (BC_RECTANGLE_LID): { r = particle[j].radius; if (particle[j].yc < BCYMIN + r) particle[j].fy += KSPRING_BOUNDARY*(BCYMIN + r - particle[j].yc); else if (particle[j].yc > ylid - r) { fperp = KSPRING_BOUNDARY*(particle[j].yc - ylid + r); particle[j].fy -= fperp; } if (particle[j].yc < BCYMAX + r) { if (particle[j].xc > BCXMAX - r) particle[j].fx -= KSPRING_BOUNDARY*(particle[j].xc - BCXMAX + r); else if (particle[j].xc < BCXMIN + r) particle[j].fx += KSPRING_BOUNDARY*(BCXMIN + r - particle[j].xc); } return(fperp); } case (BC_RECTANGLE_WALL): { padding = particle[j].radius + 0.01; xmin = BCXMIN + padding; xmax = BCXMAX - padding; ymin = BCYMIN + padding; ymax = BCYMAX - padding; if (particle[j].xc > xmax) { fperp = KSPRING_BOUNDARY*(particle[j].xc - xmax); particle[j].fx -= fperp; tmp_pright += fperp; } else if (particle[j].xc < xmin) { fperp = KSPRING_BOUNDARY*(xmin - particle[j].xc); particle[j].fx += fperp; tmp_pleft += fperp; } if (particle[j].yc > ymax) { fperp = KSPRING_BOUNDARY*(particle[j].yc - ymax); particle[j].fy -= fperp; if (particle[j].xc > xwall) tmp_pright += fperp; else tmp_pleft += fperp; } else if (particle[j].yc < ymin) { fperp = KSPRING_BOUNDARY*(ymin - particle[j].yc); particle[j].fy += fperp; if (particle[j].xc > xwall) tmp_pright += fperp; else tmp_pleft += fperp; } if (wall) { *pleft += tmp_pleft/(2.0*(BCYMAX - BCYMIN) + 2.0*(xwall - BCXMIN)); *pright += tmp_pright/(2.0*(BCYMAX - BCYMIN) + 2.0*(BCXMAX - xwall)); } else { *pleft += tmp_pleft/(2.0*(BCYMAX - BCYMIN + BCXMAX - BCXMIN)); *pright += tmp_pright/(2.0*(BCYMAX - BCYMIN + BCXMAX - BCXMIN)); } if ((wall)&&(vabs(particle[j].xc - xwall) < 0.5*WALL_WIDTH + padding)) { if (particle[j].xc > xwall) { fperp = -KSPRING_BOUNDARY*(xwall + 0.5*WALL_WIDTH + padding - particle[j].xc); particle[j].fx -= fperp; *pright -= fperp/(BCYMAX - BCYMIN); } else { fperp = KSPRING_BOUNDARY*(particle[j].xc - xwall + 0.5*WALL_WIDTH + padding); particle[j].fx -= fperp; *pleft += fperp/(BCYMAX - BCYMIN); } return(fperp); } return(0.0); } case (BC_EHRENFEST): { rp = particle[j].radius; xtube = 1.0 - sqrt(EHRENFEST_RADIUS*EHRENFEST_RADIUS - EHRENFEST_WIDTH*EHRENFEST_WIDTH); distance = 0.0; /* middle tube */ if (vabs(particle[j].xc) <= xtube) { if (particle[j].yc > EHRENFEST_WIDTH - rp) { distance = particle[j].yc - EHRENFEST_WIDTH; particle[j].fy -= KSPRING_BOUNDARY*(distance + rp); } else if (particle[j].yc < -EHRENFEST_WIDTH + rp) { distance = - EHRENFEST_WIDTH - particle[j].yc; particle[j].fy += KSPRING_BOUNDARY*(distance + rp); } } /* right container */ else if (particle[j].xc > 0.0) { r = module2(particle[j].xc - 1.0, particle[j].yc); if ((r > EHRENFEST_RADIUS - rp)&&((particle[j].xc > 1.0)||(vabs(particle[j].yc) > EHRENFEST_WIDTH))) { cphi = (particle[j].xc - 1.0)/r; sphi = particle[j].yc/r; f = KSPRING_BOUNDARY*(EHRENFEST_RADIUS - r - rp); particle[j].fx += f*cphi; particle[j].fy += f*sphi; *pright -= f; } } /* left container */ else { r = module2(particle[j].xc + 1.0, particle[j].yc); if ((r > EHRENFEST_RADIUS - rp)&&((particle[j].xc < -1.0)||(vabs(particle[j].yc) > EHRENFEST_WIDTH))) { cphi = (particle[j].xc + 1.0)/r; sphi = particle[j].yc/r; f = KSPRING_BOUNDARY*(EHRENFEST_RADIUS - r - rp); particle[j].fx += f*cphi; particle[j].fy += f*sphi; *pleft -= f; } } /* add force from "corners" */ if ((vabs(particle[j].xc) - xtube < rp)&&(vabs(particle[j].yc) - EHRENFEST_WIDTH < rp)) { for (i=-1; i<=1; i+=2) for (k=-1; k<=1; k+=2) { distance = module2(particle[j].xc - (double)i*xtube, particle[j].yc - (double)k*EHRENFEST_WIDTH); if (distance < rp) { cphi = (particle[j].xc - (double)i*xtube)/distance; sphi = (particle[j].yc - (double)k*EHRENFEST_WIDTH)/distance; f = KSPRING_BOUNDARY*(rp - distance); particle[j].fx += f*cphi; particle[j].fy += f*sphi; } } } return(fperp); } case (BC_SCREEN_BINS): { /* add harmonic force outside screen */ if (particle[j].xc > XMAX) particle[j].fx -= KSPRING_BOUNDARY*(particle[j].xc - XMAX); else if (particle[j].xc < XMIN) particle[j].fx += KSPRING_BOUNDARY*(XMIN - particle[j].xc); if (particle[j].yc > YMAX + 10.0*MU) particle[j].fy -= KSPRING_BOUNDARY*(particle[j].yc - YMAX - 10.0*MU); else if (particle[j].yc < YMIN) particle[j].fy += KSPRING_BOUNDARY*(YMIN - particle[j].yc); /* force from the bins */ dy = (YMAX - YMIN)/((double)NGRIDX + 3); dx = dy/cos(PI/6.0); rp = particle[j].radius; width = rp + 0.05*dx; ybin = 2.75*dy; if (particle[j].yc < YMIN + ybin) for (i=-1; i<=NGRIDX; i++) { x = ((double)i - 0.5*(double)NGRIDX + 0.5)*dx; distance = vabs(particle[j].xc - x); if (distance < width) { if (particle[j].xc > x) particle[j].fx += KSPRING_BOUNDARY*(width - distance); else particle[j].fx -= KSPRING_BOUNDARY*(width - distance); } } else if (particle[j].yc < YMIN + ybin + particle[j].radius) for (i=-1; i<=NGRIDX; i++) { x = ((double)i - 0.5*(double)NGRIDX + 0.5)*dx; distance = module2(particle[j].xc - x, particle[j].yc - YMIN - ybin); if (distance < rp) { if (distance < 1.0e-8) distance = 1.0e-8; cphi = (particle[j].xc - x)/distance; sphi = (particle[j].yc - YMIN - ybin)/distance; f = KSPRING_BOUNDARY*(rp - distance); particle[j].fx += f*cphi; particle[j].fy += f*sphi; } } return(fperp); } case (BC_ABSORBING): { /* add harmonic force outside screen */ padding = 0.1; x = particle[j].xc; y = particle[j].yc; if ((x > BCXMAX + padding)||(x < BCXMIN - padding)||(y > BCYMAX + padding)||(y < BCYMIN - padding)) { particle[j].active = 0; particle[j].vx = 0.0; particle[j].vy = 0.0; particle[j].xc = BCXMAX + 2.0*padding; particle[j].yc = BCYMAX + 2.0*padding; } return(fperp); } case (BC_REFLECT_ABS): { /* add harmonic force outside screen */ padding = 0.1; x = particle[j].xc; y = particle[j].yc; if ((x > BCXMAX + padding)||(y > BCYMAX + padding)||(y < BCYMIN - padding)) { particle[j].active = 0; particle[j].vx = 0.0; particle[j].vy = 0.0; particle[j].xc = BCXMAX + 2.0*padding; particle[j].yc = BCYMAX + 2.0*padding; } else if (particle[j].yc < BCYMIN) particle[j].fy += KSPRING_BOUNDARY*(BCYMIN - particle[j].yc); return(fperp); } case (BC_REFLECT_ABS_BOTTOM): { /* add harmonic force outside screen */ padding = MU; x = particle[j].xc; y = particle[j].yc; if (y < BCYMIN) { particle[j].active = 0; particle[j].vx = 0.0; particle[j].vy = 0.0; particle[j].xc = BCXMAX + 2.0*padding; particle[j].yc = BCYMIN - 2.0*padding; } else { if (y > BCYMAX + padding) particle[j].fy -= KSPRING_BOUNDARY*(y - BCYMAX - padding); if (x > BCXMAX - padding) particle[j].fx -= KSPRING_BOUNDARY*(x - BCXMAX + padding); else if (x < BCXMIN + padding) particle[j].fx += KSPRING_BOUNDARY*(BCXMIN - x + padding); } return(0.0); } case (BC_REFLECT_ABS_RIGHT): { /* add harmonic force outside screen */ padding = MU; x = particle[j].xc; y = particle[j].yc; if (x > BCXMAX + padding) { particle[j].active = 0; particle[j].vx = 0.0; particle[j].vy = 0.0; particle[j].xc = BCXMAX + 2.0*padding; particle[j].yc = BCYMIN - 2.0*padding; } else { if (y > BCYMAX - padding) particle[j].fy -= KSPRING_BOUNDARY*(y - BCYMAX + padding); else if (y < BCYMIN + padding) particle[j].fy += KSPRING_BOUNDARY*(BCYMIN - y + padding); if (x < BCXMIN + padding) particle[j].fx += KSPRING_BOUNDARY*(BCXMIN - x + padding); } return(0.0); } } } 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*(BCXMAX - BCXMIN)) dx -= (BCXMAX-BCXMIN); else if (dx < -0.5*(BCXMAX - BCXMIN)) dx += (BCXMAX-BCXMIN); if (dy > 0.5*(BCYMAX - BCYMIN)) dy -= (BCYMAX-BCYMIN); else if (dy < -0.5*(BCYMAX - BCYMIN)) dy += (BCYMAX-BCYMIN); if ((PAIR_PARTICLES)&&(PAIRING_TYPE == POLY_SEG_POLYGON)) { rmax = 1.0*particle[j].radius; for (p=-1; p<=1; p+=2) for (q=-1; q<=1; q+=2) { x1 = particle[n].xc + (double)p*particle[n].radius*cos(particle[n].angle); y1 = particle[n].yc + (double)p*particle[n].radius*sin(particle[n].angle); x2 = particle[j].xc + (double)q*particle[n].radius*cos(particle[j].angle); y2 = particle[j].yc + (double)q*particle[n].radius*sin(particle[j].angle); r = module2(x2-x1, y2-y1); if ((r>0.0)&&(r < rmax)) { // f[0] += KSPRING_PAIRS*(x1-x2)/r; // f[1] += KSPRING_PAIRS*(y1-y2)/r; f[0] = 1.0e8*(x1-x2)/r; f[1] = 1.0e8*(y1-y2)/r; } } } // else { r = module2(dx, dy); if (r < 1.0e-10) r = 1.0e-10; 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; /* 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; /* TEST */ f[0] -= particle[j].vx*DAMPING; f[1] -= particle[j].vy*DAMPING; } if (ROTATION) { *torque = 0.0; if ((REACTION_DIFFUSION)&&(RD_REACTION == CHEM_POLYGON_AGGREGATION)) { // r = module2(dx, dy); // if (r < 2.0*MU) // { // if (NPOLY%2 == 0) // sangle = sin((double)NPOLY*(particle[n].angle - particle[j].angle)); // else sangle = sin((double)NPOLY*(particle[n].angle - particle[j].angle) - PI); // sangle = sin((double)NPOLY*(particle[n].angle - particle[j].angle - eq_angle)); // *torque += KTORQUE_PAIRS*sangle; /* stabilize relative angle */ q = particle[j].partner[0]; if ((particle[j].npartners == 1)&&(particle[q].npartners == 1)) { sangle = sin((double)NPOLY*(particle[n].angle - particle[j].angle - eq_angle)); *torque += KTORQUE_PAIRS*sangle; phi = argument(dx, dy); sangle = sin((double)NPOLY*(phi - particle[j].angle) - PI); *torque += KTORQUE_PAIRS*sangle; } /* TEST: add some damping to rotation */ *torque -= particle[j].omega*DAMPING_ROT; // sangle = sin((double)NPOLY*(phi - particle[j].angle) - PI); // *torque += KTORQUE_PAIRS*sangle; // } } else if ((PAIR_PARTICLES)&&(PAIRING_TYPE == POLY_SEG_POLYGON)) { /* TODO */ sangle = sin(particle[n].angle - particle[j].angle - eq_angle); // if (sangle > 0.0) *torque = KTORQUE_PAIRS*(1.0 + sangle); // else *torque = KTORQUE_PAIRS*(-1.0 + sangle); *torque = KTORQUE_PAIRS*sangle; // *torque = 0.0; // if (COMPUTE_PAIR_TORQUE) // { // torque2 = KTORQUE_PAIR_ANGLE*sangle; // f[0] -= torque2*sa; // f[1] += torque2*ca; // } } else { 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; } } } } void compute_particle_force(int j, double krepel, t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HASHX*HASHY]) /* compute force from other particles on particle j */ { int i0, j0, m0, k, l, m, q, close, n, test, different_cluster = 1; double fx = 0.0, fy = 0.0, force[2], torque = 0.0, torque_ij, x, y; particle[j].neighb = 0; if (REACTION_DIFFUSION) particle[j].diff_neighb = 0; for (k=0; k TYPE_PROPORTION); case (TTC_CHESSBOARD): { switch (CIRCLE_PATTERN) { case (C_SQUARE): return((i%NGRIDY + i/NGRIDY)%2 == 0); case (C_HEX): return(i%2 == 0); default: return(i%2 == 0); } } case (TTC_COANDA): { if (vabs(particle[i].yc) < LAMBDA) return(1); if (vabs(particle[i].yc) < LAMBDA + MU) return(-1); return(0); } default: return(0); } } int initialize_configuration(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HASHX*HASHY], t_obstacle obstacle[NMAXOBSTACLES], double px[NMAXCIRCLES], double py[NMAXCIRCLES], double pangle[NMAXCIRCLES], int tracer_n[N_TRACER_PARTICLES], t_segment segment[NMAXSEGMENTS], t_molecule molecule[NMAXCIRCLES]) /* initialize all particles, obstacles, and the hashgrid */ { int i, j, k, n, type, nactive = 0, hashcell; double x, y, h, xx, yy, rnd, angle; for (i=0; i < ncircles; i++) { /* set particle type */ particle[i].type = 0; // if ((TWO_TYPES)&&((double)rand()/RAND_MAX > TYPE_PROPORTION)) // if ((TWO_TYPES)&&(twotype_config(i, particle))) if (TWO_TYPES) { type = twotype_config(i, particle); if (type == 1) { particle[i].type = 1; // particle[i].type = 2; particle[i].radius = MU_B; } else if (type == -1) particle[i].active = 0; } if ((INTERACTION == I_VICSEK_SHARK)&&(i==1)) { particle[i].type = 2; particle[i].radius = MU_B; } particle[i].neighb = 0; particle[i].diff_neighb = 0; particle[i].thermostat = 1; particle[i].close_to_boundary = 0; particle[i].emean = 0.0; particle[i].damping = 1.0; particle[i].dirmean = 0.0; particle[i].paired = 0; particle[i].collision = 0; // particle[i].added = 0; // particle[i].coulomb = 1; // particle[i].energy = 0.0; // y = particle[i].yc; // if (y >= YMAX) y -= particle[i].radius; // if (y <= YMIN) y += particle[i].radius; if (RANDOM_RADIUS) particle[i].radius = particle[i].radius*(RANDOM_RADIUS_MIN + RANDOM_RADIUS_RANGE*((double)rand()/RAND_MAX)); if (particle[i].type == 0) { particle[i].interaction = INTERACTION; particle[i].eq_dist = EQUILIBRIUM_DIST; particle[i].spin_range = SPIN_RANGE; particle[i].spin_freq = SPIN_INTER_FREQUENCY; particle[i].mass_inv = 1.0/PARTICLE_MASS; if ((RANDOM_RADIUS)&&(ADAPT_MASS_TO_RADIUS > 0.0)) particle[i].mass_inv *= 1.0/(1.0 + pow(particle[i].radius/MU, ADAPT_MASS_TO_RADIUS)); if ((RANDOM_RADIUS)&&(ADAPT_DAMPING_TO_RADIUS > 0.0)) particle[i].damping = 1.0 + ADAPT_DAMPING_FACTOR*pow(particle[i].radius/MU, ADAPT_DAMPING_TO_RADIUS); particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT; particle[i].charge = CHARGE; } else { particle[i].interaction = INTERACTION_B; particle[i].eq_dist = EQUILIBRIUM_DIST_B; particle[i].spin_range = SPIN_RANGE_B; particle[i].spin_freq = SPIN_INTER_FREQUENCY_B; particle[i].mass_inv = 1.0/PARTICLE_MASS_B; particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT_B; particle[i].charge = CHARGE_B; } switch (V_INITIAL_TYPE) { case (VI_RANDOM): { particle[i].vx = V_INITIAL*gaussian(); particle[i].vy = V_INITIAL*gaussian(); break; } case (VI_COANDA): { if (vabs(particle[i].yc) < LAMBDA) particle[i].vx = V_INITIAL; else particle[i].vx = 0.0; particle[i].vy = 0.0; break; } } if ((INTERACTION == I_VICSEK_SHARK)&&(i==1)) { particle[i].vx *= 1000.0; particle[i].vy *= 1000.0; } particle[i].energy = (particle[i].vx*particle[i].vx + particle[i].vy*particle[i].vy)*particle[i].mass_inv; particle[i].emean = particle[i].energy; particle[i].dirmean = 0.0; px[i] = particle[i].vx; py[i] = particle[i].vy; if (ROTATION) { particle[i].angle = DPI*(double)rand()/RAND_MAX; particle[i].omega = OMEGA_INITIAL*gaussian(); if (COUPLE_ANGLE_TO_THERMOSTAT) particle[i].energy += particle[i].omega*particle[i].omega*particle[i].inertia_moment_inv; } else { particle[i].angle = 0.0; particle[i].omega = 0.0; } pangle[i] = particle[i].omega; if ((PLOT == P_INITIAL_POS)||(PLOT_B == P_INITIAL_POS)) { switch (INITIAL_POS_TYPE) { case (IP_X): { particle[i].color_hue = 360.0*(particle[i].xc - INITXMIN)/(INITXMAX - INITXMIN); break; } case (IP_Y): { particle[i].color_hue = 360.0*(particle[i].yc - INITYMIN)/(INITYMAX - INITYMIN); break; } } } if ((PLOT == P_NUMBER)||(PLOT_B == P_NUMBER)) particle[i].color_hue = 360.0*(double)(i%N_PARTICLE_COLORS)/(double)N_PARTICLE_COLORS; // if ((PLOT == P_CLUSTER)||(PLOT_B == P_CLUSTER)) { if (PAIR_PARTICLES) particle[i].cluster = rand()%(ncircles*CLUSTER_COLOR_FACTOR); else particle[i].cluster = rand()%ncircles; } } /* initialize dummy values in case particles are added */ for (i=ncircles; i < NMAXCIRCLES; i++) { particle[i].type = 0; particle[i].active = 0; particle[i].neighb = 0; particle[i].thermostat = 0; particle[i].energy = 0.0; particle[i].emean = 0.0; particle[i].damping = 1.0; particle[i].dirmean = 0.0; particle[i].mass_inv = 1.0/PARTICLE_MASS; particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT; particle[i].charge = 0.0; particle[i].vx = 0.0; particle[i].vy = 0.0; px[i] = 0.0; py[i] = 0.0; particle[i].angle = DPI*(double)rand()/RAND_MAX; particle[i].omega = 0.0; pangle[i] = 0.0; particle[i].interaction = INTERACTION; particle[i].eq_dist = EQUILIBRIUM_DIST; particle[i].spin_range = SPIN_RANGE; particle[i].spin_freq = SPIN_INTER_FREQUENCY; particle[i].close_to_boundary = 0; particle[i].coulomb = 1; particle[i].added = 1; particle[i].reactive = 1; particle[i].paired = 0; particle[i].collision = 0; } /* add particles at the bottom as seed */ if (PART_AT_BOTTOM) for (i=0; i<=NPART_BOTTOM; i++) { x = XMIN + (double)i*(XMAX - XMIN)/(double)NPART_BOTTOM; y = YMIN + 2.0*MU; add_particle(x, y, 0.0, 0.0, MASS_PART_BOTTOM, 0, particle); } if (PART_AT_BOTTOM) for (i=0; i<=NPART_BOTTOM; i++) { x = XMIN + (double)i*(XMAX - XMIN)/(double)NPART_BOTTOM; y = YMIN + 4.0*MU; add_particle(x, y, 0.0, 0.0, MASS_PART_BOTTOM, 0, particle); } /* add larger copies of particles (for Ehrenfest model)*/ if (EHRENFEST_COPY) { for (i=0; i < ncircles; i++) { n = ncircles + i; particle[n].xc = -particle[i].xc; particle[n].yc = particle[i].yc; particle[n].vx = -0.5*particle[i].vx; particle[n].vy = 0.5*particle[i].vy; px[n] = -0.5*px[i]; py[n] = 0.5*py[i]; particle[n].energy = particle[i].energy; particle[n].radius = 2.0*particle[i].radius; particle[n].type = 2; particle[n].mass_inv = 1.25*particle[i].mass_inv; particle[n].thermostat = 1; particle[n].interaction = particle[i].interaction; particle[n].eq_dist = 0.45*particle[i].eq_dist; if ((double)rand()/RAND_MAX > 0.6) particle[n].active = 1; } ncircles *= 2; } /* change type of tracer particle */ if (TRACER_PARTICLE) for (j=0; j 0.5)) i++; while ((!particle[i].active)||(module2(particle[i].xc + xx, particle[i].yc - yy) > 0.4)) i++; tracer_n[j] = i; particle[i].type = 2 + j; particle[i].radius *= 1.5; particle[i].mass_inv *= 1.0/TRACER_PARTICLE_MASS; particle[i].vx *= 0.1; particle[i].vy *= 0.1; particle[i].thermostat = 0; px[i] *= 0.1; py[i] *= 0.1; n_tracers++; } /* position-dependent particle type */ if (POSITION_DEPENDENT_TYPE) for (i=0; i EHRENFEST_RADIUS) particle[i].active = 0; } else if (BOUNDARY_COND == BC_RECTANGLE_WALL) { for (i=0; i< ncircles; i++) if (vabs(particle[i].xc - xwall) < WALL_WIDTH) particle[i].active = 0; } else if (BOUNDARY_COND == BC_GENUS_TWO) { for (i=0; i< ncircles; i++) if ((particle[i].xc > 0.0)&&(particle[i].yc > 0.0)) particle[i].active = 0; } if (ADD_FIXED_OBSTACLES) { for (i=0; i< ncircles; i++) for (j=0; j < nobstacles; j++) { if (module2(particle[i].xc - obstacle[j].xc, particle[i].yc - obstacle[j].yc) < OBSTACLE_RADIUS + particle[i].radius) particle[i].active = 0; hashcell = hash_cell(obstacle[j].xc, obstacle[j].yc); hashgrid[hashcell].charge = OBSTACLE_CHARGE; } } /* case of segment obstacles */ if (ADD_FIXED_SEGMENTS) for (i=0; i< ncircles; i++) if (!in_segment_region(particle[i].xc, particle[i].yc, segment)) particle[i].active = 0; /* case of reaction-diffusion equation/chemical reactions */ if (REACTION_DIFFUSION) for (i=0; i< ncircles; i++) { switch (RD_INITIAL_COND) { case (IC_UNIFORM): { particle[i].type = 1; break; } case (IC_UNIFORM2): { particle[i].type = 2; particle[i].radius = MU_B; particle[i].omega = OMEGA_INITIAL*gaussian(); break; } case (IC_RANDOM_UNIF): { particle[i].type = 1 + (int)(RD_TYPES*(double)rand()/(double)RAND_MAX); break; } case (IC_RANDOM_TWO): { if ((double)rand()/(double)RAND_MAX < TYPE_PROPORTION) particle[i].type = 1; else { particle[i].type = 2; particle[i].radius = MU_B; particle[i].mass_inv = 1.0/PARTICLE_MASS_B; } break; } case (IC_CIRCLE): { if (module2(particle[i].xc,particle[i].yc) < LAMBDA) particle[i].type = 1; else { particle[i].type = 2; particle[i].radius = MU_B; particle[i].mass_inv = 1.0/PARTICLE_MASS_B; } break; } case (IC_CATALYSIS): { if ((particle[i].xc > 0.0)||((double)rand()/(double)RAND_MAX < TYPE_PROPORTION)) { particle[i].type = 1; particle[i].radius = MU; particle[i].mass_inv = 1.0/PARTICLE_MASS; } else { particle[i].type = 2; particle[i].radius = MU_B; particle[i].mass_inv = 1.0/PARTICLE_MASS_B; } break; } case (IC_LAYERS): { if (particle[i].yc > 0.0) { particle[i].type = 1; particle[i].radius = MU; particle[i].eq_dist = EQUILIBRIUM_DIST; particle[i].mass_inv = 1.0/PARTICLE_MASS; } else if (particle[i].yc < -LAMBDA) { particle[i].type = 2; particle[i].radius = MU_B; particle[i].eq_dist = EQUILIBRIUM_DIST_B; particle[i].mass_inv = 1.0/PARTICLE_MASS_B; } else particle[i].active = 0; break; } case (IC_BZ): { rnd = (double)rand()/(double)RAND_MAX; if (rnd < 60.0/180.0) { particle[i].type = 4; particle[i].radius = MU; particle[i].eq_dist = EQUILIBRIUM_DIST; particle[i].mass_inv = 1.0/(PARTICLE_MASS + 3.0*PARTICLE_MASS_B); } else if (rnd < 100.0/180.0) { particle[i].type = 5; particle[i].radius = MU_B; particle[i].eq_dist = EQUILIBRIUM_DIST; particle[i].mass_inv = 1.0/(3.5*PARTICLE_MASS); } else { particle[i].type = 6; particle[i].radius = 1.5*MU_B; particle[i].eq_dist = EQUILIBRIUM_DIST; particle[i].mass_inv = 0.2/(PARTICLE_MASS); } break; } case (IC_SIGNX): { if (particle[i].xc < 0.0) particle[i].type = 1; else { particle[i].type = 2; particle[i].radius = MU_B; particle[i].mass_inv = 1.0/PARTICLE_MASS_B; } break; } case (IC_TWOROCKETS): { if (vabs(particle[i].xc) < SEGMENTS_X0) particle[i].type = 1; else { particle[i].type = 2; particle[i].radius = MU_B; particle[i].mass_inv = 1.0/PARTICLE_MASS_B; } break; } case (IC_TWOROCKETS_TWOFUELS): { if (vabs(particle[i].xc) < SEGMENTS_X0) particle[i].type = 1; else { if (particle[i].xc < 0) particle[i].type = 2; else particle[i].type = 3; particle[i].radius = MU_B; particle[i].mass_inv = 1.0/PARTICLE_MASS_B; } break; } case (IC_DNA_POLYMERASE): { /* do nothing? */ break; } case (IC_DNA_POLYMERASE_REC): { /* do nothing? */ break; } default: { /* do nothing */ } } } /* keep only active particles */ // ncircles = reorder_particles(particle, py, pangle); /* count number of active particles */ for (i=0; i< ncircles; i++) nactive += particle[i].active; printf("%i active particles\n", nactive); // for (i=0; i 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; if ((TRACER_PARTICLE)&&(n_tracers < N_TRACER_PARTICLES)) { tracer_n[n_tracers] = ncircles - 1; n_tracers++; printf("%i tracer particles\n", n_tracers); } // init_molecule_data(particle, molecule); return (ncircles); // add_particle(XMIN + 0.1, 0.0, 50.0, 0.0, 3.0, 0, particle); // px[ncircles - 1] = particle[ncircles - 1].vx; // py[ncircles - 1] = particle[ncircles - 1].vy; // particle[ncircles - 1].radius = 1.5*MU; // j = 0; // while (module2(particle[j].xc,particle[j].yc) > 0.7) j = rand()%ncircles; // x = particle[j].xc + 2.5*MU; // y = particle[j].yc; // x = XMIN + (XMAX - XMIN)*rand()/RAND_MAX; // y = YMAX + 0.01*rand()/RAND_MAX; // add_particle(x, y, 0.0, 0.0, 1.0, 0, particle); // x = XMIN + 0.25*(XMAX - XMIN); // y = YMAX + 0.01; // prop = 1.0 - (double)nadd_particle/5.0; // vx = 100.0*prop; // add_particle(x, y, vx, -10.0, 5.0*prop, 0, particle); // particle[ncircles - 1].radius = 10.0*MU*prop; // particle[ncircles - 1].eq_dist = 2.0; // particle[ncircles - 1].thermostat = 0; // px[ncircles - 1] = particle[ncircles - 1].vx; // py[ncircles - 1] = particle[ncircles - 1].vy; // add_particle(MU*(2.0*rand()/RAND_MAX - 1.0), YMAX + 2.0*MU, 0.0, 0.0, PARTICLE_MASS, 0, particle); // add_particle(XMIN - 0.5*MU, 0.0, 50.0 + 5.0*(double)i, 0.0, 2.0*PARTICLE_MASS, 0, particle); // x = INITXMIN + (INITXMAX - INITXMIN)*(double)rand()/(double)RAND_MAX; // y = INITYMIN + (INITYMAX - INITYMIN)*(double)rand()/(double)RAND_MAX; // x = BCXMIN + (BCXMAX - BCXMIN)*(double)rand()/(double)RAND_MAX; // y = YMAX + 0.5*(BCYMAX - YMAX)*(double)rand()/(double)RAND_MAX; } void center_momentum(double p[NMAXCIRCLES]) { int i; double ptot = 0.0, pmean; for (i=0; i PMAX) { p[i] = PMAX; floor = 1; } else if (p[i] < -PMAX) { p[i] = -PMAX; floor = 1; } } if (floor) printf("Flooring momentum\n"); return (floor); } int partial_thermostat_coupling(t_particle particle[NMAXCIRCLES], double xmin, t_segment segment[NMAXSEGMENTS], t_lj_parameters params) /* only couple particles satisfying condition PARTIAL_THERMO_REGION to thermostat */ { int condition, i, nthermo = 0; double x, y, height, a, b; static double maxheight; static int first = 1; if (first) { maxheight = YMIN + PARTIAL_THERMO_HEIGHT*(YMAX - YMIN); first = 0; } if (PARTIAL_THERMO_REGION == TH_LAYER_TYPE2) { height = YMIN; for (i=0; i height)&&(y <= maxheight)) height = y; } if (height > maxheight) height = maxheight; printf("Thermostat region y > %.3lg, max height = %.3lg\n", height, maxheight); } for (i=0; i xmin); break; } case (TH_INSEGMENT): { condition = (in_segment_region(particle[i].xc, particle[i].yc, segment)); // condition = (in_segment_region(particle[i].xc - xsegments[0], particle[i].yc - ysegments[0])); // // condition = (in_segment_region(particle[i].xc - xsegments[0], particle[i].yc - ysegments[0])); // if (TWO_OBSTACLES) // condition = condition||(in_segment_region(particle[i].xc - xsegments[1], particle[i].yc - ysegments[1])); break; } case (TH_INBOX): { x = particle[i].xc; y = particle[i].yc; condition = ((y < YMIN + PARTIAL_THERMO_HEIGHT*(YMAX - YMIN))&&(vabs(x) < PARTIAL_THERMO_WIDTH*XMAX)); break; } case (TH_LAYER): { y = particle[i].yc; condition = (y > PARTIAL_THERMO_HEIGHT); break; } case (TH_LAYER_TYPE2): { y = particle[i].yc; condition = (y > height); break; } case (TH_RING): { x = particle[i].xc; y = particle[i].yc; condition = x*x + y*y > PARTIAL_THERMO_WIDTH*PARTIAL_THERMO_WIDTH; break; } case (TH_RING_EXPAND): { x = particle[i].xc; y = particle[i].yc; condition = x*x + y*y > params.thermo_radius*params.thermo_radius; break; } case (TH_INIT): { x = particle[i].xc; y = particle[i].yc; condition = ((x > INITXMIN)&&(x < INITXMAX)&&(y > INITYMIN)&&(y < INITYMAX)); break; } case (TH_THERMO): { x = particle[i].xc; y = particle[i].yc; condition = ((x > THERMOXMIN)&&(x < THERMOXMAX)&&(y > THERMOYMIN)&&(y < THERMOYMAX)); break; } case (TH_CONE): { x = particle[i].xc; y = particle[i].yc; a = (LAMBDA + 0.2)/(LAMBDA - 0.05); b = LAMBDA*(1.0 - a); condition = ((y < LAMBDA)&&(y > 0.0)&&(y > a*vabs(x) + b)); break; } default: condition = 1; } if (condition) { particle[i].thermostat = 1; nthermo++; } else particle[i].thermostat = 0; } return(nthermo); } int partial_efield(double x, double y) /* */ { switch (EFIELD_REGION) { case (EF_CONST): return(1); case (EF_SQUARE): return(vabs(x) < YMAX); case (EF_LEFT): return(x < YMIN); } } int partial_bfield(double x, double y) /* */ { switch (BFIELD_REGION) { case (BF_CONST): return(1); case (BF_SQUARE): return(vabs(x) < YMAX); } } double compute_mean_energy(t_particle particle[NMAXCIRCLES]) { int i, nactive = 0; double total_energy = 0.0; for (i=0; ivx*(particle->vx) + particle->vy*(particle->vy))*(particle->mass_inv); if (old_energy + delta_ekin > 0.0) e_ratio = sqrt((old_energy + delta_ekin)/old_energy); else e_ratio = 0.0; particle->vx *= e_ratio; particle->vy *= e_ratio; } int chem_merge_AAB(int i, int newtype, t_particle particle[NMAXCIRCLES], t_collision *collisions, int ncollisions, double inv_masses[RD_TYPES+1], double radii[RD_TYPES+1]) /* merging particle i with a particle of same type into a particle of type newtype */ /* particular case of chem_merge */ { int j, k, type1; short int search = 1; double distance, rx, ry, mr, newmass_inv; type1 = particle[i].type; newmass_inv = inv_masses[newtype]; mr = inv_masses[newtype]/inv_masses[type1]; for (j=0; j 10) { printf("Error: need to increase size of k1 table in chem_multi_merge()\n"); exit(1); } type1 = particle[i].type; newmass_inv = inv_masses[newtype]; mr1 = newmass_inv/inv_masses[type1]; mr2 = newmass_inv/inv_masses[type2]; n1 = 0; n2 = 0; for (j=0; j= 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, 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; for (p=0; p 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; } } } } } } 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= cluster[j].nparticles) { il = i; is = j; pl = pi; ps = pj; } else { il = j; is = i; pl = pj; ps = pi; } ns = cluster[is].nparticles; nl = cluster[il].nparticles; if (ns + nl >= NMAXPARTINCLUSTER) { printf("Warning: NMAXPARTINCLUSTER is too small\n"); printf("Try increasing to %i\n", cluster[i].nparticles + cluster[j].nparticles); return(0); } printf("Merging clusters %i and %i, having %i and %i particles\n", il, is, nl, ns); fprintf(lj_log, "Merging clusters %i and %i, having %i and %i particles\n", il, is, nl, ns); fprintf(lj_log, "Large cluster %i:\n", il); for (p=0; p= alpha2) particle[q].angle -= alpha2; while (particle[q].angle < 0.0) particle[q].angle += alpha2; } q = cluster[il].particle[0]; newangle = particle[q].angle; cluster[il].angle = newangle; // if (NPOLY%2==1) newangle += alpha; for (p=nl; p maxpartners) return(ncollisions); m1 = 1.0/particle[i].mass_inv; for (p=0; p maxpartners) return(0); if ((particle[k].active)&&(particle[k].type == type2)&&(nq <= maxpartners)) { distance = module2(particle[i].deltax[p], particle[i].deltay[p]); rel_angle = argument(particle[i].deltax[p], particle[i].deltay[p]); poly_condition = (distance < REACTION_DIST*(ri+rk)); rorient = particle[k].angle - particle[i].angle; if (NPOLY%2 == 1) rorient -= alpha; norient = (int)(rorient/alpha2 + 0.5); delta = rorient - norient*alpha2; poly_condition *= (vabs(delta) < DELTAMAX); if ((poly_condition)&&((double)rand()/RAND_MAX < reaction_prob)) { /* find nearest orientation facing particle k */ rorient_new = norient*alpha2; nangle = (int)((particle[i].angle - rel_angle)/alpha2); anglei_new = rel_angle + (double)(2*nangle+1)*alpha; m2 = 1.0/particle[k].mass_inv; for (q=0; q (%.3lg, %.3lg)\n", particle[q].xc, particle[q].yc, newx, newy); dx = newx - particle[q].xc; dy = newy - particle[q].yc; r = module2(dx, dy); // printf("Particle %i distance to target = %.3lg\n", k, r); // if ((t == nsteps-1)&&(q == 470)&&(k == 472)) fprintf(lj_log, "Particle %i relative distance to target %i = %.5lg\n", q, k, r/eqdist); fx += r*dx; fy += r*dy; } else { eqdist2 = 1.0*(particle[k].radius + particle[q].radius)*sa; if (dist < eqdist2) { fx -= 0.1*(eqdist2 - dist)*dx/dist; fy -= 0.1*(eqdist2 - dist)*dy/dist; } } } if ((kill_overlays)&&(dist < REPAIR_MIN_DIST*eqdist)) { nq = particle[q].npartners; np = particle[k].npartners; if (nq == 1) particle[q].active = 0; if (np == 1) particle[k].active = 0; // if (np > nq) particle[q].active = 0; // else particle[k].active = 0; } } } particle[q].xc += kspring*fx*DT_PARTICLE; particle[q].yc += kspring*fy*DT_PARTICLE; } } } int chem_multi_glue_polygon2(int i, int maxpartners, t_particle particle[NMAXCIRCLES], t_molecule molecule[NMAXCIRCLES], t_cluster cluster[NMAXCIRCLES], t_collision *collisions, int ncollisions, double reaction_prob, int mergeclusters, int singlecluster) /* simplified version of chem_multi_glue_molecule without triangle and similar */ /* conditions, but taking polygon geometry into account */ /* version 2, for CHEM_POLYGON_CLUSTER interaction */ { int k, p, q, p1, q1, p2, np, nq, n, m, closeby = 0, reaction = 0, jp[NMAXPARTNERS], jq[NMAXPARTNERS], jtot[NMAXPARTNERS], r, kk, different, np2, moli, molk, mp, newpartner, pp, type1, poly_condition, nangle, norient, cl1, cl2, cp, cq, newcluster, nmin, clmin, clmax, newreaction; double distance, angle, move_factor, xnew, ynew, vxnew, vynew, m1, m2, mr1, mr2, deltav, x, y, ri, rk, alpha, alpha2, ca, rel_angle, rorient, delta, rorient_new, anglei_new, deltax, deltay, x1, y1, aratio, delta2, eqdistance, dx, dy, dalpha1, dalpha2, ssize1, ssize2; static short int first = 1; static int scluster; if ((singlecluster)&&(first)) { scluster = 0; while (cluster[scluster].active == 0) scluster = rand()%ncircles; cluster[scluster].selected = 1; particle[scluster].collision = 1; printf("Selected cluster %i as single cluster\n", scluster); first = 0; } type1 = particle[i].type; np = particle[i].npartners; if (np > maxpartners) return(ncollisions); cl1 = particle[i].cluster; cp = cluster[cl1].nparticles; m1 = cluster[cl1].mass; alpha = PI/(double)NPOLY; alpha2 = 2.0*alpha; ca = cos(alpha); ri = particle[i].radius*ca; rk = particle[k].radius*ca; p = 0; while ((p maxpartners) return(ncollisions); if ((!reaction)&&(particle[k].active)&&(np+nq+1 < maxpartners)&&(cl2!=cl1)&&(nmin <= SMALL_NP_MAXSIZE)&&(clmin <= SMALL_CLUSTER_MAXSIZE)) { cq = cluster[cl2].nparticles; // poly_condition = (cl1!=cl2); /* test distance */ distance = module2(particle[i].deltax[p], particle[i].deltay[p]); rel_angle = argument(particle[i].deltax[p], particle[i].deltay[p]); poly_condition = (distance < REACTION_DIST*(ri+rk)); /* test relative angle */ rorient = particle[k].angle - particle[i].angle; if (NPOLY%2 == 1) rorient -= alpha; while (rorient > alpha) rorient -= alpha2; while (rorient < -alpha) rorient += alpha2; poly_condition *= (vabs(rorient) < DELTAMAX); // norient = (int)(rorient/alpha2 + 0.5 + 2.0*(double)NPOLY); // delta = rorient - (double)(norient-2*NPOLY)*alpha2; // poly_condition *= (vabs(delta) < DELTAMAX); /* test orientation */ aratio = 0.5*vabs(rel_angle)/alpha; delta2 = aratio - (double)((int)aratio); poly_condition *= (vabs(delta2 - 0.5) < 0.5*DELTAMAX*alpha); /* option singlecluster: only allow merger if one of the clusters is selected cluster */ if (singlecluster) poly_condition *= ((cluster[cl1].selected)||(cluster[cl2].selected)||(clmax <= NOTSELECTED_CLUSTER_MAXSIZE)); newreaction = 0; if ((poly_condition)&&((double)rand()/RAND_MAX < reaction_prob)) { /* find nearest orientation facing particle k */ rorient_new = norient*alpha2; nangle = (int)((particle[i].angle - rel_angle)/alpha2); anglei_new = rel_angle + (double)(2*nangle+1)*alpha; m2 = cluster[cl2].mass; mr1 = m1/(m1 + m2); mr2 = 1.0 - mr1; if ((np == 0)&&(nq == 0)) { reaction = 1; newreaction = 1; printf("Merging molecule %i with particle %i\n", i, k); eqdistance = ri + rk; particle[i].angle = anglei_new; particle[k].angle = particle[i].angle + rorient_new; if (NPOLY%2 == 1) particle[k].angle += alpha; dx = 0.5*(distance - eqdistance)*cos(rel_angle); dy = 0.5*(distance - eqdistance)*sin(rel_angle); translate_cluster(cl1, cluster, particle, dx, dy); translate_cluster(cl2, cluster, particle, -dx, -dy); particle[i].npartners++; particle[i].partner[0] = k; particle[k].npartners++; particle[k].partner[0] = i; /* equalize speeds */ vxnew = mr1*particle[i].vx + mr2*particle[k].vx; vynew = mr1*particle[i].vy + mr2*particle[k].vy; particle[i].vx = vxnew; particle[i].vy = vynew; particle[k].vx = vxnew; particle[k].vy = vynew; if (mergeclusters) newcluster = merge_clusters(particle[i].cluster, particle[k].cluster, i, k, particle, cluster, 0, (particle[i].flip == particle[k].flip)); } // else if ((np+nq+2 < NMAXPARTNERS)&&(vabs(delta) < 0.5*DELTAMAX)) else if ((np+nq+2 < NMAXPARTNERS)&&(cp+cq <= CLUSTER_MAXSIZE)) { deltav = module2(particle[i].vx - particle[k].vx, particle[i].vy - particle[k].vy); // printf("Delta v is %.3lg\n", deltav); if (deltav < DELTAVMAX) { reaction = 1; newreaction = 1; printf("Pairing clusters containing particles %i and %i\n", i, k); printf("np = %i, nq = %i\n", np, nq); fprintf(lj_log, "Pairing clusters containing particles %i and %i\n", i, k); fprintf(lj_log, "np = %i, nq = %i\n", np, nq); eqdistance = ri + rk; dx = 0.5*(distance - eqdistance)*cos(rel_angle); dy = 0.5*(distance - eqdistance)*sin(rel_angle); dalpha1 = anglei_new - particle[i].angle; dalpha2 = anglei_new + rorient_new - particle[k].angle; if (NPOLY%2 == 1) dalpha2 += alpha; /* round to closest multiple of alpha2 */ while (dalpha1 > 0.5*DELTAMAX) dalpha1 -= alpha2; while (dalpha1 < -0.5*DELTAMAX) dalpha1 += alpha2; while (dalpha2 > 0.5*DELTAMAX) dalpha2 -= alpha2; while (dalpha2 < -0.5*DELTAMAX) dalpha2 += alpha2; translate_cluster(cl1, cluster, particle, dx, dy); translate_cluster(cl2, cluster, particle, -dx, -dy); ssize1 = 0.02*sqrt((double)cluster[cl1].nparticles); ssize2 = 0.02*sqrt((double)cluster[cl2].nparticles); if (vabs(dalpha1) + ssize1 < 0.5*DELTAMAX) rotate_cluster_around_particle(cl1, i, cluster, particle, dalpha1); else printf("Did not rotate cluster of size %.3lg\n", ssize1); if (vabs(dalpha2) + ssize2 < 0.5*DELTAMAX) rotate_cluster_around_particle(cl2, k, cluster, particle, dalpha2); else printf("Did not rotate cluster of size %.3lg\n", ssize2); particle[i].npartners++; particle[i].partner[np] = k; particle[k].npartners++; particle[k].partner[nq] = i; /* equalize speeds */ vxnew = mr1*particle[i].vx + mr2*particle[k].vx; vynew = mr1*particle[i].vy + mr2*particle[k].vy; for (p1=0; p1 SMALL_NP_MAXSIZE)||(clmin > SMALL_CLUSTER_MAXSIZE)) // { // printf("No merger - Parameters: nmin = %i, clmin = %i\n", nmin, clmin); // } p++; } return(ncollisions); } int chem_split_molecule(int i, t_particle particle[NMAXCIRCLES], t_molecule molecule[NMAXCIRCLES], t_collision *collisions, int ncollisions) /* split molecule containing particle i from all other molecules */ { int np, nmolp, mol, n, m, p, q, r, mm; int m_table[NMAXPARTNERMOLECULES]; short int reaction = 0; np = particle[i].npartners; if (np==0) return(ncollisions); mol = particle[i].molecule; nmolp = molecule[mol].npartners; printf("Molecule %i has %i partners: ", mol, nmolp); fprintf(lj_log, "[chem_split_molecule] Molecule %i has %i partners: ", mol, nmolp); for (mm=0; mm= 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 */ { int i, n0=0, n1=0, n2=0, ntot, nc=0, nmod, counter=0, cmax=1000; for (i=0; i 0) while ((nc prop*(double)ntot) { nmod = (int)((double)n0 - (double)ntot*prop); if (nmod > 0) while ((ncactive) return(-1); pstack[position]->cluster = cluster; for (k=0; knpartners; k++) { n = pstack[position]->partner[k]; if ((particle[n].active)&&(particle[n].cluster != cluster)) { particle[n].tested = 1; particle[n].cluster = cluster; particle[n].cactive = 1; nopen++; if (*stacksize < ncircles) { (*stacksize)++; pstack[*stacksize-1] = &particle[n]; } } } if (nopen == 0) { pstack[position]->cactive = 0; return(-1); } else return(nopen); } int update_cluster_color(t_particle particle[NMAXCIRCLES]) /* update the colors of clusters */ { int i, k, p, position, nactive, stacksize, nclusters = 0; t_particle **cstack; cstack = (t_particle* *)malloc(ncircles*sizeof(struct t_particle *)); #pragma omp parallel for private(i) for (i=0; i 0) { while (!cstack[position]->cactive) position++; if (position == stacksize) position = 0; /* find an active partner in stack */ nactive += find_partners(particle, position, cstack, &stacksize, particle[i].cluster); } nclusters++; } free(cstack); return(nclusters); } 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_cluster cluster[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, atype, btype, oldncollisions, delta_n; short int reacted; double distance, rnd, p1, p2, reac_dist; static double inv_masses[RD_TYPES+1], radii[RD_TYPES+1]; static int first = 1, repair_counter = 0; if (first) /* compute total mass and mass ratios */ { compute_inverse_masses(inv_masses); compute_radii(radii); first = 0; } if (EXOTHERMIC) oldncollisions = ncollisions; switch (RD_REACTION) { case (CHEM_RPS): { for (i=0; i 2)&&((double)rand()/RAND_MAX < DISSOCIATION_PROB)) ncollisions = chem_split(i, 1, particle[i].type - 1, particle, collisions, ncollisions, inv_masses, radii); } printf("%i collisions\n", ncollisions); return(ncollisions); } case (CHEM_POLYMER_STEP): { for (i=0; i 1)&&((double)rand()/RAND_MAX < DISSOCIATION_PROB)) { k = rand()%(type-1) + 1; // printf("Splitting type %i into type %i and type %i\n", type, k, type - k); ncollisions = chem_split(i, k, type - k, particle, collisions, ncollisions, inv_masses, radii); } } printf("%i collisions\n", ncollisions); return(ncollisions); } case (CHEM_CATALYTIC_A2D): { for (i=0; i X */ { ncollisions = chem_convert(i, 3, particle, collisions, ncollisions, inv_masses, radii, DISSOCIATION_PROB); break; } case (2): /* B + X -> Y + D */ { ncollisions = chem_transfer(i, 3, 4, 5, particle, collisions, ncollisions, inv_masses, radii, REACTION_DIST); break; } case (3): { /* X -> B, modified from Brusselator */ ncollisions = chem_convert(i, 2, particle, collisions, ncollisions, inv_masses, radii, 0.5*DISSOCIATION_PROB); /* 2X + Y -> 3X */ ncollisions = chem_catalytic_convert(i, 4, 3, particle, collisions, ncollisions, inv_masses, radii, 3.0*REACTION_DIST, 1.0); break; } case (5): /* D -> A or B if concentration of X or Y is small */ { n = time*(RD_TYPES+1); n3 = particle_numbers[n+3]; n4 = particle_numbers[n+4]; p1 = 1.0/((double)(n3*n3/10+1)); p2 = 5.0*DISSOCIATION_PROB + 1.0/((double)(n4*n4/200+1)); rnd = (double)rand()/RAND_MAX; if (rnd < p1/(p1+p2)) ncollisions = chem_convert(i, 1, particle, collisions, ncollisions, inv_masses, radii, p1); else /*if (n4 < 1000)*/ ncollisions = chem_convert(i, 2, particle, collisions, ncollisions, inv_masses, radii, p2); } } } return(ncollisions); } case (CHEM_ABDACBE): { 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_AGGREGATION_CHARGE): { 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_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= 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_POLYGON_CLUSTER): { for (i=0; i= AGREG_DECOUPLE) particle[i].thermostat = 0; } /* repair clusters */ if (!REPAIR_CLUSTERS) for (i=0; i= 2)) repair_cluster(i, particle, cluster, 1000, 0); printf("%i collisions\n", ncollisions); delta_n = ncollisions - oldncollisions; printf("delta_n = %i\n", delta_n); return(ncollisions); }case (CHEM_POLYGON_ONECLUSTER): { i = 0; reacted = 0; while ((i oldncollisions) reacted = 1; } /* decouple particles with several partners from thermostat */ if (particle[i].npartners >= AGREG_DECOUPLE) particle[i].thermostat = 0; } i++; } /* repair clusters */ if (!REPAIR_CLUSTERS) for (i=0; i= 2)) repair_cluster(i, particle, cluster, 1000, 1); printf("%i collisions\n", ncollisions); delta_n = ncollisions - oldncollisions; printf("delta_n = %i\n", delta_n); return(ncollisions); } } } double plot_coord(double x, double xmin, double xmax) { return(xmin + x*(xmax - xmin)); } void draw_speed_plot(t_group_data *group_speeds, int i) /* draw plot of obstacle speeds as a function of time */ { int j, group; char message[100]; static double xmin, xmax, ymin, ymax, xmid, ymid, dx, dy, plotxmin, plotxmax, plotymin, plotymax; double pos[2], x1, y1, x2, y2, rgb[3]; static int first = 1, gshift = INITIAL_TIME + NSTEPS; if (first) { // xmin = XMAX - 1.35; // xmax = XMAX - 0.05; // ymin = YMAX - 1.8; // ymax = YMAX - 0.5; xmin = XMAX - 1.8; xmax = XMAX - 0.066; ymin = YMAX - 2.5; ymax = YMAX - 0.76; xmid = 0.5*(xmin + xmax); ymid = 0.5*(ymin + ymax); dx = 0.5*(xmax - xmin); dy = 0.5*(ymax - ymin); plotxmin = xmin + 0.05; plotxmax = xmax - 0.1; plotymin = ymin + 0.07; plotymax = ymax - 0.15; first = 0; } // rgb[0] = 1.0; rgb[1] = 1.0; rgb[2] = 1.0; glLineWidth(2); /* plot angular speed */ for (group=1; group scaley) scaley = y0; x2 = plot_coord(0.5 + x0/scalex, plotxmin, plotxmax); y2 = plot_coord(y0/scaley + 0.05, plotymin, plotymax); // printf("yinitial = %.3lg, (x0, y0) = (%.3lg, %.3lg), (x2, y2) = (%.3lg, %.3lg)\n", yinitial, x0, y0, x2, y2); if ((j>0)&&(module2(x2-x1, y2-y1) < 0.1)) draw_line(x1, y1, x2, y2); x1 = x2; y1 = y2; } if (i>0) draw_colored_circle_precomp(x1, y1, 0.015, rgb); } glColor3f(1.0, 1.0, 1.0); /* axes and labels */ draw_line(plotxmin, plotymin, plotxmax + 0.05, plotymin); draw_line(plotxmin, plotymin, plotxmin, plotymax + 0.1); for (j=1; j<(int)scaley; j++) { y1 = plot_coord((double)j/scaley, plotymin, plotymax); draw_line(plotxmin - 0.02, y1, plotxmin + 0.02, y1); } sprintf(message, "%i", (int)scaley - 1); write_text_fixedwidth(plotxmin - 0.15, y1 - 0.025, message); // sprintf(message, "%.1f", VMAX_PLOT_SPEEDS); sprintf(message, "y"); write_text_fixedwidth(plotxmin - 0.1, plotymax - 0.005, message); sprintf(message, "x"); write_text_fixedwidth(plotxmax - 0.1, plotymin - 0.1, message); } void write_size_text(double x, double y, char *message, int largewin) { if (largewin) write_text(x, y, message); else write_text_fixedwidth(x, y, message); } void draw_particle_nb_plot(int *particle_numbers, int i) /* draw plot of number of particles as a function of time */ { int j, type, power; char message[100]; static double xmin, xmax, ymin, ymax, xmid, ymid, dx, dy, plotxmin, plotxmax, plotymin, plotymax; double pos[2], x1, y1, x2, y2, rgb[3]; static int second = 2, gshift = INITIAL_TIME + NSTEPS, nmax, ygrad, largewin; if (second > 0) { xmin = XMAX - 1.05; xmax = XMAX - 0.05; ymin = YMAX - 1.0; ymax = YMAX - 0.05; xmid = 0.5*(xmin + xmax); ymid = 0.5*(ymin + ymax); dx = 0.5*(xmax - xmin); dy = 0.5*(ymax - ymin); plotxmin = xmin + 0.18; plotxmax = xmax - 0.1; plotymin = ymin + 0.07; plotymax = ymax - 0.15; nmax = 0; for (j=1; j nmax) nmax = particle_numbers[RD_TYPES+1+j]; nmax *= PARTICLE_NB_PLOT_FACTOR; power = (int)(log((double)nmax)/log(10.0)) - 1.0; ygrad = (int)ipow(10.0, power); largewin = (WINWIDTH > 1280); second--; } erase_area_hsl(xmid, ymid, dx, dy, 0.0, 0.9, 0.0); glLineWidth(2); /* plot particle number */ for (type=1; type= 1000) write_size_text(plotxmin - 0.15, y1 - 0.015, message, largewin); // else write_size_text(plotxmin - 0.12, y1 - 0.015, message, largewin); if (j*ygrad >= 1000) write_text_fixedwidth(plotxmin - 0.15, y1 - 0.015, message); else write_text_fixedwidth(plotxmin - 0.12, y1 - 0.015, message); } j++; } sprintf(message, "time"); // write_size_text(plotxmax - 0.1, plotymin - 0.05, message, largewin); write_text_fixedwidth(plotxmax - 0.1, plotymin - 0.05, message); } void init_segment_group(t_segment segment[NMAXSEGMENTS], int group, t_group_segments segment_group[NMAXGROUPS]) /* initialize center of mass and similar data of grouped segments */ { int i, nseg_group = 0; double xc = 0.0, yc = 0.0; for (i=0; i emax) { printf("Particle %i at (%.3lg, %.3lg) has energy %.5lg, resetting to 0\n", i, particle[i].xc, particle[i].yc, particle[i].energy); particle[i].vx = 0.0; particle[i].vy = 0.0; px[i] = 0.0; py[i] = 0.0; particle[i].energy = 0.0; } printf("\n\n"); } int init_cluster_config(t_particle particle[NMAXCIRCLES], t_cluster cluster[NMAXCIRCLES]) /* initialize the clusters for option CLUSTER_PARTICLES */ /* returns number of active clusters */ { int i, j, tmp, nclusters = 0; for (i=0; i cluster %i -> particle %i\n", i, particle[i].cluster, cluster[particle[i].cluster].particle[0]); // } return(nclusters); } void compute_cluster_force(t_cluster cluster[NMAXCIRCLES], t_particle particle[NMAXCIRCLES]) /* compute force and torque on clusters from force and torque on their particles */ { int i, p, j; double fx, fy, torque, energy; for (i=0; i beltlength) newpos -= beltlength; belt[b].shovel_pos[shovel] = newpos; } } } void draw_frame(int i, int plot, int bg_color, int ncollisions, int traj_position, int traj_length, int wall, double pressure[N_PRESSURES], double pleft, double pright, int *particle_numbers, short int refresh, t_lj_parameters params, t_particle particle[NMAXCIRCLES], t_cluster cluster[NMAXCIRCLES], t_collision *collisions, t_hashgrid hashgrid[HASHX*HASHY], t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTICLES], t_obstacle obstacle[NMAXOBSTACLES], t_segment segment[NMAXSEGMENTS], t_group_data *group_speeds, t_group_segments *segment_group, t_belt *conveyor_belt, int *tracer_n, t_otriangle otriangle[NMAX_TRIANGLES_PER_OBSTACLE*NMAXOBSTACLES], t_ofacet ofacet[NMAXOBSTACLES]) /* draw a movie frame */ { printf("Drawing frame\n"); if ((COLOR_BACKGROUND)&&(bg_color > 0)) color_background(particle, obstacle, bg_color, hashgrid); // else if (!TRACER_PARTICLE) blank(); if (DRAW_OBSTACLE_LINKS) draw_obstacle_triangles(obstacle, otriangle, ofacet); // if (DRAW_OBSTACLE_LINKS) old_draw_obstacle_triangles(obstacle); if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length, particle, cluster, tracer_n, plot); draw_particles(particle, cluster, plot, params.beta, collisions, ncollisions, bg_color, hashgrid, params); draw_container(params.xmincontainer, params.xmaxcontainer, obstacle, segment, conveyor_belt, wall); if (PRINT_PARAMETERS) print_parameters(params, PRINT_LEFT, pressure, refresh); if (PLOT_SPEEDS) draw_speed_plot(group_speeds, i); if (PLOT_TRAJECTORIES) draw_trajectory_plot(group_speeds, i); if ((i > INITIAL_TIME)&&(PLOT_PARTICLE_NUMBER)) draw_particle_nb_plot(particle_numbers, i - INITIAL_TIME); // if (PLOT_PARTICLE_NUMBER) draw_particle_nb_plot(particle_numbers, i - INITIAL_TIME); if (BOUNDARY_COND == BC_EHRENFEST) print_ehrenfest_parameters(particle, pleft, pright); else if (PRINT_PARTICLE_NUMBER) { if (REACTION_DIFFUSION) print_particle_types_number(particle, RD_TYPES); else print_particle_number(ncircles); } else if (PRINT_PARTICLE_SPEEDS) print_particles_speeds(particle); else if (PRINT_SEGMENTS_SPEEDS) print_segment_group_speeds(segment_group); // printf("5\n"); }