1 /* -*- Mode: C; tab-width: 4 -*- */
2 /* euler2d --- 2 Dimensional Incompressible Inviscid Fluid Flow */
5 static const char sccsid[] = "@(#)euler2d.c 5.00 2000/11/01 xlockmore";
9 * Copyright (c) 2000 by Stephen Montgomery-Smith <stephen@math.missouri.edu>
11 * Permission to use, copy, modify, and distribute this software and its
12 * documentation for any purpose and without fee is hereby granted,
13 * provided that the above copyright notice appear in all copies and that
14 * both that copyright notice and this permission notice appear in
15 * supporting documentation.
17 * This file is provided AS IS with no warranties of any kind. The author
18 * shall have no liability with respect to the infringement of copyrights,
19 * trade secrets or any patents by this file or any part thereof. In no
20 * event will the author be liable for any lost revenue or profits or
21 * other special, indirect and consequential damages.
24 * 04-Nov-2000: Added an option eulerpower. This allows for example the
25 * quasi-geostrophic equation by setting eulerpower to 2.
26 * 01-Nov-2000: Allocation checks.
27 * 10-Sep-2000: Added optimizations, and removed subtle_perturb, by stephen.
28 * 03-Sep-2000: Changed method of solving ode to Adams-Bashforth of order 2.
29 * Previously used a rather compilcated method of order 4.
30 * This doubles the speed of the program. Also it seems
31 * to have improved numerical stability. Done by stephen.
32 * 27-Aug-2000: Added rotation of region to maximize screen fill by stephen.
33 * 05-Jun-2000: Adapted from flow.c Copyright (c) 1996 by Tim Auckland
34 * 18-Jul-1996: Adapted from swarm.c Copyright (c) 1991 by Patrick J. Naughton.
35 * 31-Aug-1990: Adapted from xswarm by Jeff Butterworth. (butterwo@ncsc.org)
39 * The mathematical aspects of this program are discussed in the file
45 #define PROGCLASS "Euler2d"
46 #define HACK_INIT init_euler2d
47 #define HACK_DRAW draw_euler2d
48 #define euler2d_opts xlockmore_opts
49 #define DEFAULTS "*delay: 10000 \n" \
54 #include "xlockmore.h" /* in xscreensaver distribution */
55 #else /* STANDALONE */
56 #include "xlock.h" /* in xlockmore distribution */
57 #endif /* STANDALONE */
61 #define DEF_EULERTAIL "10"
63 #define DEBUG_POINTED_REGION 0
66 static int variable_boundary = 1;
67 static float power = 1;
69 static XrmOptionDescRec opts[] =
71 {(char* ) "-eulertail", (char *) ".euler2d.eulertail",
72 XrmoptionSepArg, (caddr_t) NULL},
73 {(char* ) "-eulerpower", (char *) ".euler2d.eulerpower",
74 XrmoptionSepArg, (caddr_t) NULL},
76 static argtype vars[] =
78 {(caddr_t *) &tail_len, (char *) "eulertail",
79 (char *) "EulerTail", (char *) DEF_EULERTAIL, t_Int},
80 {(caddr_t *) &power, (char *) "eulerpower",
81 (char *) "EulerPower", (char *) "1", t_Float},
83 static OptionStruct desc[] =
85 {(char *) "-eulertail len", (char *) "Length of Euler2d tails"},
86 {(char *) "-eulerpower power", (char *) "power of interaction law for points for Euler2d"},
89 ModeSpecOpt euler2d_opts =
90 {sizeof opts / sizeof opts[0], opts,
91 sizeof vars / sizeof vars[0], vars, desc};
94 ModStruct euler2d_description = {
95 "euler2d", "init_euler2d", "draw_euler2d", "release_euler2d",
96 "refresh_euler2d", "init_euler2d", (char *) NULL, &euler2d_opts,
97 1000, 1024, 3000, 1, 64, 1.0, "",
98 "Simulates 2D incompressible invisid fluid.", 0, NULL
103 #define balance_rand(v) ((LRAND()/MAXRAND*(v))-((v)/2)) /* random around 0 */
104 #define positive_rand(v) (LRAND()/MAXRAND*(v)) /* positive random */
106 #define number_of_vortex_points 20
108 #define n_bound_p 500
111 static double delta_t;
117 double xshift,yshift,scale;
123 /* x[2i+0] = x coord for nth point
124 x[2i+1] = y coord for nth point
125 w[i] = vorticity at nth point
134 /* (xs[2i+0],xs[2i+1]) is reflection of (x[2i+0],x[2i+1]) about unit circle
135 xs[2i+0] = x[2i+0]/nx
136 xs[2i+1] = x[2i+1]/nx
138 nx = x[2i+0]*x[2i+0] + x[2i+1]*x[2i+1]
140 x_is_zero[i] = (nx < 1e-10)
145 /* (p[2i+0],p[2i+1]) is image of (x[2i+0],x[2i+1]) under polynomial p.
146 mod_dp2 is |p'(z)|^2 when z = (x[2i+0],x[2i+1]).
151 /* Sometimes in our calculations we get overflow or numbers that are too big.
152 If that happens with the point x[2*i+0], x[2*i+1], we set dead[i].
165 double p_coef[2*(deg_p-1)];
170 static euler2dstruct *euler2ds = (euler2dstruct *) NULL;
173 If variable_boundary == 1, then we make a variable boundary.
174 The way this is done is to map the unit disk under a
176 p(z) = z + c_2 z^2 + ... + c_n z^n
177 where n = deg_p. sp->p_coef contains the complex numbers
181 #define add(a1,a2,b1,b2) (a1)+=(b1);(a2)+=(b2)
182 #define mult(a1,a2,b1,b2) temp=(a1)*(b1)-(a2)*(b2); \
183 (a2)=(a1)*(b2)+(a2)*(b1);(a1)=temp
186 calc_p(double *p1, double *p2, double z1, double z2, double p_coef[])
193 for(i=deg_p;i>=2;i--)
195 add(*p1,*p2,p_coef[(i-2)*2],p_coef[(i-2)*2+1]);
202 /* Calculate |p'(z)|^2 */
204 calc_mod_dp2(double z1, double z2, double p_coef[])
211 for(i=deg_p;i>=2;i--)
213 add(mp1,mp2,i*p_coef[(i-2)*2],i*p_coef[(i-2)*2+1]);
217 return mp1*mp1+mp2*mp2;
221 calc_all_p(euler2dstruct *sp)
224 double temp,p1,p2,z1,z2;
225 for(j=(sp->hide_vortex?sp->Nvortex:0);j<sp->N;j++) if(!sp->dead[j])
231 for(i=deg_p;i>=2;i--)
233 add(p1,p2,sp->p_coef[(i-2)*2],sp->p_coef[(i-2)*2+1]);
244 calc_all_mod_dp2(double *x, euler2dstruct *sp)
247 double temp,mp1,mp2,z1,z2;
248 for(j=0;j<sp->N;j++) if(!sp->dead[j])
254 for(i=deg_p;i>=2;i--)
256 add(mp1,mp2,i*sp->p_coef[(i-2)*2],i*sp->p_coef[(i-2)*2+1]);
260 sp->mod_dp2[j] = mp1*mp1+mp2*mp2;
265 derivs(double *x, euler2dstruct *sp)
268 double u1,u2,x1,x2,xij1,xij2,nxij;
271 if (variable_boundary)
272 calc_all_mod_dp2(sp->x,sp);
274 for (j=0;j<sp->Nvortex;j++) if (!sp->dead[j])
276 nx = x[2*j+0]*x[2*j+0] + x[2*j+1]*x[2*j+1];
278 sp->x_is_zero[j] = 1;
280 sp->x_is_zero[j] = 0;
281 sp->xs[2*j+0] = x[2*j+0]/nx;
282 sp->xs[2*j+1] = x[2*j+1]/nx;
286 (void) memset(sp->diffx,0,sizeof(double)*2*sp->N);
288 for (i=0;i<sp->N;i++) if (!sp->dead[i])
292 for (j=0;j<sp->Nvortex;j++) if (!sp->dead[j])
295 Calculate the Biot-Savart kernel, that is, effect of a
296 vortex point at a = (x[2*j+0],x[2*j+1]) at the point
297 x = (x1,x2), returning the vector field in (u1,u2).
299 In the plane, this is given by the formula
301 u = (x-a)/|x-a|^2 or zero if x=a.
303 However, in the unit disk we have to subtract from the
308 where as = a/|a|^2 is the reflection of a about the unit circle.
310 If however power != 1, then
312 u = (x-a)/|x-a|^(power+1) - |a|^(1-power) (x-as)/|x-as|^(power+1)
316 xij1 = x1 - x[2*j+0];
317 xij2 = x2 - x[2*j+1];
318 nxij = (power==1.0) ? xij1*xij1+xij2*xij2 : pow(xij1*xij1+xij2*xij2,(power+1)/2.0);
327 if (!sp->x_is_zero[j])
329 xij1 = x1 - sp->xs[2*j+0];
330 xij2 = x2 - sp->xs[2*j+1];
331 nxij = (power==1.0) ? xij1*xij1+xij2*xij2 : pow(xij1*xij1+xij2*xij2,(power+1)/2.0);
347 sp->diffx[2*i+0] += u1*sp->w[j];
348 sp->diffx[2*i+1] += u2*sp->w[j];
352 if (!sp->dead[i] && variable_boundary)
354 if (sp->mod_dp2[i] < 1e-5)
358 sp->diffx[2*i+0] /= sp->mod_dp2[i];
359 sp->diffx[2*i+1] /= sp->mod_dp2[i];
366 What perturb does is effectively
368 where k should be of order delta_t.
370 We have the option to do this more subtly by mapping points x
371 in the unit disk to points y in the plane, where y = f(|x|) x,
372 with f(t) = -log(1-t)/t.
374 This might reduce (but does not remove) problems where particles near
375 the edge of the boundary bounce around.
377 But it seems to be not that effective, so for now switch it off.
380 #define SUBTLE_PERTURB 0
383 perturb(double ret[], double x[], double k[], euler2dstruct *sp)
389 double d1,d2,t1,t2,mag,mag2,mlog1mmag,memmagdmag,xdotk;
390 for (i=0;i<sp->N;i++) if (!sp->dead[i])
396 mag2 = x1*x1 + x2*x2;
402 else if (mag2 > 1-1e-5)
407 mlog1mmag = -log(1-mag);
408 xdotk = x1*k1 + x2*k2;
409 t1 = (x1 + k1)*mlog1mmag/mag + x1*xdotk*(1.0/(1-mag)-mlog1mmag/mag)/mag/mag;
410 t2 = (x2 + k2)*mlog1mmag/mag + x2*xdotk*(1.0/(1-mag)-mlog1mmag/mag)/mag/mag;
411 mag = sqrt(t1*t1+t2*t2);
412 if (mag > 11.5 /* log(1e5) */)
416 memmagdmag = (mag>1e-5) ? ((1.0-exp(-mag))/mag) : (1-mag/2.0);
417 ret[2*i+0] = t1*memmagdmag;
418 ret[2*i+1] = t2*memmagdmag;
425 if (d1*d1+d2*d2 > 0.1)
432 for (i=0;i<sp->N;i++) if (!sp->dead[i])
438 if (k1*k1+k2*k2 > 0.1 || x1*x1+x2*x2 > 1-1e-5)
450 ode_solve(euler2dstruct *sp)
456 /* midpoint method */
458 (void) memcpy(sp->olddiffx,sp->diffx,sizeof(double)*2*sp->N);
459 for (i=0;i<sp->N;i++) if (!sp->dead[i]) {
460 sp->tempdiffx[2*i+0] = 0.5*delta_t*sp->diffx[2*i+0];
461 sp->tempdiffx[2*i+1] = 0.5*delta_t*sp->diffx[2*i+1];
463 perturb(sp->tempx,sp->x,sp->tempdiffx,sp);
464 derivs(sp->tempx,sp);
465 for (i=0;i<sp->N;i++) if (!sp->dead[i]) {
466 sp->tempdiffx[2*i+0] = delta_t*sp->diffx[2*i+0];
467 sp->tempdiffx[2*i+1] = delta_t*sp->diffx[2*i+1];
469 perturb(sp->x,sp->x,sp->tempdiffx,sp);
473 for (i=0;i<sp->N;i++) if (!sp->dead[i]) {
474 sp->tempdiffx[2*i+0] = delta_t*(1.5*sp->diffx[2*i+0] - 0.5*sp->olddiffx[2*i+0]);
475 sp->tempdiffx[2*i+1] = delta_t*(1.5*sp->diffx[2*i+1] - 0.5*sp->olddiffx[2*i+1]);
477 perturb(sp->x,sp->x,sp->tempdiffx,sp);
479 sp->olddiffx = sp->diffx;
484 #define deallocate(p,t) if (p!=NULL) {(void) free((void *) p); p=(t*)NULL; }
485 #define allocate(p,t,s) if ((p=(t*)malloc(sizeof(t)*s))==NULL)\
486 {free_euler2d(sp);return;}
489 free_euler2d(euler2dstruct *sp)
491 deallocate(sp->csegs, XSegment);
492 deallocate(sp->old_segs, XSegment);
493 deallocate(sp->nold_segs, int);
494 deallocate(sp->lastx, short);
495 deallocate(sp->x, double);
496 deallocate(sp->diffx, double);
497 deallocate(sp->w, double);
498 deallocate(sp->olddiffx, double);
499 deallocate(sp->tempdiffx, double);
500 deallocate(sp->tempx, double);
501 deallocate(sp->dead, short);
502 deallocate(sp->boundary, XSegment);
503 deallocate(sp->xs, double);
504 deallocate(sp->x_is_zero, short);
505 deallocate(sp->p, double);
506 deallocate(sp->mod_dp2, double);
510 init_euler2d(ModeInfo * mi)
512 #define nr_rotates 18 /* how many rotations to try to fill as much of screen as possible - must be even number */
515 double r,theta,x,y,w;
516 double mag,xscale,yscale,p1,p2;
517 double low[nr_rotates],high[nr_rotates],pp1,pp2,pn1,pn2,angle1,angle2,tempangle,dist,scale,bestscale,temp;
520 if (power<0.5) power = 0.5;
521 if (power>3.0) power = 3.0;
522 variable_boundary &= power == 1.0;
524 if (power>1.0) delta_t *= pow(0.1,power-1);
526 if (euler2ds == NULL) {
527 if ((euler2ds = (euler2dstruct *) calloc(MI_NUM_SCREENS(mi),
528 sizeof (euler2dstruct))) == NULL)
531 sp = &euler2ds[MI_SCREEN(mi)];
533 sp->boundary_color = NRAND(MI_NPIXELS(mi));
534 sp->hide_vortex = NRAND(4) != 0;
538 sp->width = MI_WIDTH(mi);
539 sp->height = MI_HEIGHT(mi);
541 sp->N = MI_COUNT(mi)+number_of_vortex_points;
542 sp->Nvortex = number_of_vortex_points;
544 if (tail_len < 1) { /* minimum tail */
547 if (tail_len > MI_CYCLES(mi)) { /* maximum tail */
548 tail_len = MI_CYCLES(mi);
551 /* Clear the background. */
556 /* Allocate memory. */
558 if (sp->csegs == NULL) {
559 allocate(sp->csegs, XSegment, sp->N);
560 allocate(sp->old_segs, XSegment, sp->N * tail_len);
561 allocate(sp->nold_segs, int, tail_len);
562 allocate(sp->lastx, short, sp->N * 2);
563 allocate(sp->x, double, sp->N * 2);
564 allocate(sp->diffx, double, sp->N * 2);
565 allocate(sp->w, double, sp->Nvortex);
566 allocate(sp->olddiffx, double, sp->N * 2);
567 allocate(sp->tempdiffx, double, sp->N * 2);
568 allocate(sp->tempx, double, sp->N * 2);
569 allocate(sp->dead, short, sp->N);
570 allocate(sp->boundary, XSegment, n_bound_p);
571 allocate(sp->xs, double, sp->Nvortex * 2);
572 allocate(sp->x_is_zero, short, sp->Nvortex);
573 allocate(sp->p, double, sp->N * 2);
574 allocate(sp->mod_dp2, double, sp->N);
576 for (i=0;i<tail_len;i++) {
577 sp->nold_segs[i] = 0;
580 (void) memset(sp->dead,0,sp->N*sizeof(short));
582 if (variable_boundary)
584 /* Initialize polynomial p */
586 The polynomial p(z) = z + c_2 z^2 + ... c_n z^n needs to be
587 a bijection of the unit disk onto its image. This is achieved
588 by insisting that sum_{k=2}^n k |c_k| <= 1. Actually we set
589 the inequality to be equality (to get more interesting shapes).
592 for(k=2;k<=deg_p;k++)
594 r = positive_rand(1.0/k);
595 theta = balance_rand(2*M_PI);
596 sp->p_coef[2*(k-2)+0]=r*cos(theta);
597 sp->p_coef[2*(k-2)+1]=r*sin(theta);
600 if (mag > 0.0001) for(k=2;k<=deg_p;k++)
602 sp->p_coef[2*(k-2)+0] /= mag;
603 sp->p_coef[2*(k-2)+1] /= mag;
606 #if DEBUG_POINTED_REGION
607 for(k=2;k<=deg_p;k++){
608 sp->p_coef[2*(k-2)+0]=0;
609 sp->p_coef[2*(k-2)+1]=0;
611 sp->p_coef[2*(6-2)+0] = 1.0/6.0;
615 /* Here we figure out the best rotation of the domain so that it fills as
616 much of the screen as possible. The number of angles we look at is determined
617 by nr_rotates (we look every 180/nr_rotates degrees).
618 While we figure out the best angle to rotate, we also figure out the correct scaling factors.
621 for(k=0;k<nr_rotates;k++) {
626 for(k=0;k<n_bound_p;k++) {
627 calc_p(&p1,&p2,cos((double)k/(n_bound_p)*2*M_PI),sin((double)k/(n_bound_p)*2*M_PI),sp->p_coef);
628 calc_p(&pp1,&pp2,cos((double)(k-1)/(n_bound_p)*2*M_PI),sin((double)(k-1)/(n_bound_p)*2*M_PI),sp->p_coef);
629 calc_p(&pn1,&pn2,cos((double)(k+1)/(n_bound_p)*2*M_PI),sin((double)(k+1)/(n_bound_p)*2*M_PI),sp->p_coef);
630 angle1 = nr_rotates/M_PI*atan2(p2-pp2,p1-pp1)-nr_rotates/2;
631 angle2 = nr_rotates/M_PI*atan2(pn2-p2,pn1-p1)-nr_rotates/2;
632 while (angle1<0) angle1+=nr_rotates*2;
633 while (angle2<0) angle2+=nr_rotates*2;
634 if (angle1>nr_rotates*1.75 && angle2<nr_rotates*0.25) angle2+=nr_rotates*2;
635 if (angle1<nr_rotates*0.25 && angle2>nr_rotates*1.75) angle1+=nr_rotates*2;
641 for(i=(int)floor(angle1);i<(int)ceil(angle2);i++) {
642 dist = cos((double)i*M_PI/nr_rotates)*p1 + sin((double)i*M_PI/nr_rotates)*p2;
643 if (i%(nr_rotates*2)<nr_rotates) {
644 if (dist>high[i%nr_rotates]) high[i%nr_rotates] = dist;
645 if (dist<low[i%nr_rotates]) low[i%nr_rotates] = dist;
648 if (-dist>high[i%nr_rotates]) high[i%nr_rotates] = -dist;
649 if (-dist<low[i%nr_rotates]) low[i%nr_rotates] = -dist;
654 for (i=0;i<nr_rotates;i++) {
655 xscale = (sp->width-5.0)/(high[i]-low[i]);
656 yscale = (sp->height-5.0)/(high[(i+nr_rotates/2)%nr_rotates]-low[(i+nr_rotates/2)%nr_rotates]);
657 scale = (xscale>yscale) ? yscale : xscale;
658 if (scale>bestscale) {
663 /* Here we do the rotation. The way we do this is to replace the
664 polynomial p(z) by a^{-1} p(a z) where a = exp(i best_angle).
668 for(k=2;k<=deg_p;k++)
670 mult(p1,p2,cos((double)besti*M_PI/nr_rotates),sin((double)besti*M_PI/nr_rotates));
671 mult(sp->p_coef[2*(k-2)+0],sp->p_coef[2*(k-2)+1],p1,p2);
674 sp->scale = bestscale;
675 sp->xshift = -(low[besti]+high[besti])/2.0*sp->scale+sp->width/2;
676 if (besti<nr_rotates/2)
677 sp->yshift = -(low[besti+nr_rotates/2]+high[besti+nr_rotates/2])/2.0*sp->scale+sp->height/2;
679 sp->yshift = (low[besti-nr_rotates/2]+high[besti-nr_rotates/2])/2.0*sp->scale+sp->height/2;
682 /* Initialize boundary */
684 for(k=0;k<n_bound_p;k++)
687 calc_p(&p1,&p2,cos((double)k/(n_bound_p)*2*M_PI),sin((double)k/(n_bound_p)*2*M_PI),sp->p_coef);
688 sp->boundary[k].x1 = (short)(p1*sp->scale+sp->xshift);
689 sp->boundary[k].y1 = (short)(p2*sp->scale+sp->yshift);
691 for(k=1;k<n_bound_p;k++)
693 sp->boundary[k].x2 = sp->boundary[k-1].x1;
694 sp->boundary[k].y2 = sp->boundary[k-1].y1;
696 sp->boundary[0].x2 = sp->boundary[n_bound_p-1].x1;
697 sp->boundary[0].y2 = sp->boundary[n_bound_p-1].y1;
701 if (sp->width>sp->height)
702 sp->radius = sp->height/2.0-5.0;
704 sp->radius = sp->width/2.0-5.0;
707 /* Initialize point positions */
709 for (i=sp->Nvortex;i<sp->N;i++) {
711 r = sqrt(positive_rand(1.0));
712 theta = balance_rand(2*M_PI);
713 sp->x[2*i+0]=r*cos(theta);
714 sp->x[2*i+1]=r*sin(theta);
715 /* This is to make sure the initial distribution of points is uniform */
716 } while (variable_boundary &&
717 calc_mod_dp2(sp->x[2*i+0],sp->x[2*i+1],sp->p_coef)
722 /* number of vortex points with negative vorticity */
727 /* if n is even make sure that np==n/2 is twice as likely
728 as the other possibilities. */
734 r = sqrt(positive_rand(0.77));
735 theta = balance_rand(2*M_PI);
738 r = 0.02+positive_rand(0.1);
739 w = (2*(k<np)-1)*2.0/sp->Nvortex;
740 for (i=sp->Nvortex*k/n;i<sp->Nvortex*(k+1)/n;i++) {
741 theta = balance_rand(2*M_PI);
742 sp->x[2*i+0]=x + r*cos(theta);
743 sp->x[2*i+1]=y + r*sin(theta);
750 draw_euler2d(ModeInfo * mi)
752 Display *display = MI_DISPLAY(mi);
753 Window window = MI_WINDOW(mi);
755 int b, col, n_non_vortex_segs;
758 MI_IS_DRAWN(mi) = True;
760 if (euler2ds == NULL)
762 sp = &euler2ds[MI_SCREEN(mi)];
763 if (sp->csegs == NULL)
767 if (variable_boundary)
771 for(b=sp->Nvortex;b<sp->N;b++) if(!sp->dead[b])
773 sp->csegs[sp->cnsegs].x1 = sp->lastx[2*b+0];
774 sp->csegs[sp->cnsegs].y1 = sp->lastx[2*b+1];
775 if (variable_boundary)
777 sp->csegs[sp->cnsegs].x2 = (short)(sp->p[2*b+0]*sp->scale+sp->xshift);
778 sp->csegs[sp->cnsegs].y2 = (short)(sp->p[2*b+1]*sp->scale+sp->yshift);
782 sp->csegs[sp->cnsegs].x2 = (short)(sp->x[2*b+0]*sp->radius+sp->width/2);
783 sp->csegs[sp->cnsegs].y2 = (short)(sp->x[2*b+1]*sp->radius+sp->height/2);
785 sp->lastx[2*b+0] = sp->csegs[sp->cnsegs].x2;
786 sp->lastx[2*b+1] = sp->csegs[sp->cnsegs].y2;
789 n_non_vortex_segs = sp->cnsegs;
791 if (!sp->hide_vortex) for(b=0;b<sp->Nvortex;b++) if(!sp->dead[b])
793 sp->csegs[sp->cnsegs].x1 = sp->lastx[2*b+0];
794 sp->csegs[sp->cnsegs].y1 = sp->lastx[2*b+1];
795 if (variable_boundary)
797 sp->csegs[sp->cnsegs].x2 = (short)(sp->p[2*b+0]*sp->scale+sp->xshift);
798 sp->csegs[sp->cnsegs].y2 = (short)(sp->p[2*b+1]*sp->scale+sp->yshift);
802 sp->csegs[sp->cnsegs].x2 = (short)(sp->x[2*b+0]*sp->radius+sp->width/2);
803 sp->csegs[sp->cnsegs].y2 = (short)(sp->x[2*b+1]*sp->radius+sp->height/2);
805 sp->lastx[2*b+0] = sp->csegs[sp->cnsegs].x2;
806 sp->lastx[2*b+1] = sp->csegs[sp->cnsegs].y2;
811 XSetForeground(display, gc, MI_BLACK_PIXEL(mi));
813 XDrawSegments(display, window, gc, sp->old_segs+sp->c_old_seg*sp->N, sp->nold_segs[sp->c_old_seg]);
815 if (MI_NPIXELS(mi) > 2){ /* render colour */
816 for (col = 0; col < MI_NPIXELS(mi); col++) {
817 int start = col*n_non_vortex_segs/MI_NPIXELS(mi);
818 int finish = (col+1)*n_non_vortex_segs/MI_NPIXELS(mi);
819 XSetForeground(display, gc, MI_PIXEL(mi, col));
820 XDrawSegments(display, window, gc,sp->csegs+start, finish-start);
822 if (!sp->hide_vortex) {
823 XSetForeground(display, gc, MI_WHITE_PIXEL(mi));
824 XDrawSegments(display, window, gc,sp->csegs+n_non_vortex_segs, sp->cnsegs-n_non_vortex_segs);
827 } else { /* render mono */
828 XSetForeground(display, gc, MI_WHITE_PIXEL(mi));
829 XDrawSegments(display, window, gc,
830 sp->csegs, sp->cnsegs);
833 if (MI_NPIXELS(mi) > 2) /* render colour */
834 XSetForeground(display, gc, MI_PIXEL(mi, sp->boundary_color));
836 XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_WHITE_PIXEL(mi));
837 if (variable_boundary)
838 XDrawSegments(display, window, gc,
839 sp->boundary, n_bound_p);
841 XDrawArc(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
842 sp->width/2 - (int) sp->radius - 1, sp->height/2 - (int) sp->radius -1,
843 (int) (2*sp->radius) + 2, (int) (2* sp->radius) + 2, 0, 64*360);
845 /* Copy to erase-list */
846 (void) memcpy(sp->old_segs+sp->c_old_seg*sp->N, sp->csegs, sp->cnsegs*sizeof(XSegment));
847 sp->nold_segs[sp->c_old_seg] = sp->cnsegs;
849 if (sp->c_old_seg >= tail_len)
853 if (++sp->count > MI_CYCLES(mi)) /* pick a new flow */
859 release_euler2d(ModeInfo * mi)
861 if (euler2ds != NULL) {
864 for (screen = 0; screen < MI_NUM_SCREENS(mi); screen++)
865 free_euler2d(&euler2ds[screen]);
866 (void) free((void *) euler2ds);
867 euler2ds = (euler2dstruct *) NULL;
872 refresh_euler2d(ModeInfo * mi)
877 #endif /* MODE_euler2d */