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 DEFAULTS "*delay: 10000 \n" \
49 "*fpsSolid: true \n" \
50 "*ignoreRotation: True \n" \
52 # define SMOOTH_COLORS
53 # define release_euler2d 0
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 {"-eulertail", ".euler2d.eulertail", XrmoptionSepArg, NULL},
72 {"-eulerpower", ".euler2d.eulerpower", XrmoptionSepArg, NULL},
74 static argtype vars[] =
76 {&tail_len, "eulertail",
77 "EulerTail", (char *) DEF_EULERTAIL, t_Int},
78 {&power, "eulerpower",
79 "EulerPower", "1", t_Float},
81 static OptionStruct desc[] =
83 {"-eulertail len", "Length of Euler2d tails"},
84 {"-eulerpower power", "power of interaction law for points for Euler2d"},
87 ENTRYPOINT ModeSpecOpt euler2d_opts =
88 {sizeof opts / sizeof opts[0], opts,
89 sizeof vars / sizeof vars[0], vars, desc};
92 ModStruct euler2d_description = {
93 "euler2d", "init_euler2d", "draw_euler2d", (char *) NULL,
94 "refresh_euler2d", "init_euler2d", (char *) NULL, &euler2d_opts,
95 1000, 1024, 3000, 1, 64, 1.0, "",
96 "Simulates 2D incompressible invisid fluid.", 0, NULL
101 #define balance_rand(v) ((LRAND()/MAXRAND*(v))-((v)/2)) /* random around 0 */
102 #define positive_rand(v) (LRAND()/MAXRAND*(v)) /* positive random */
104 #define number_of_vortex_points 20
106 #define n_bound_p 500
109 static double delta_t;
115 double xshift,yshift,scale;
121 /* x[2i+0] = x coord for nth point
122 x[2i+1] = y coord for nth point
123 w[i] = vorticity at nth point
132 /* (xs[2i+0],xs[2i+1]) is reflection of (x[2i+0],x[2i+1]) about unit circle
133 xs[2i+0] = x[2i+0]/nx
134 xs[2i+1] = x[2i+1]/nx
136 nx = x[2i+0]*x[2i+0] + x[2i+1]*x[2i+1]
138 x_is_zero[i] = (nx < 1e-10)
143 /* (p[2i+0],p[2i+1]) is image of (x[2i+0],x[2i+1]) under polynomial p.
144 mod_dp2 is |p'(z)|^2 when z = (x[2i+0],x[2i+1]).
149 /* Sometimes in our calculations we get overflow or numbers that are too big.
150 If that happens with the point x[2*i+0], x[2*i+1], we set dead[i].
163 double p_coef[2*(deg_p-1)];
168 static euler2dstruct *euler2ds = (euler2dstruct *) NULL;
171 If variable_boundary == 1, then we make a variable boundary.
172 The way this is done is to map the unit disk under a
174 p(z) = z + c_2 z^2 + ... + c_n z^n
175 where n = deg_p. sp->p_coef contains the complex numbers
179 #define add(a1,a2,b1,b2) (a1)+=(b1);(a2)+=(b2)
180 #define mult(a1,a2,b1,b2) temp=(a1)*(b1)-(a2)*(b2); \
181 (a2)=(a1)*(b2)+(a2)*(b1);(a1)=temp
184 calc_p(double *p1, double *p2, double z1, double z2, double p_coef[])
191 for(i=deg_p;i>=2;i--)
193 add(*p1,*p2,p_coef[(i-2)*2],p_coef[(i-2)*2+1]);
200 /* Calculate |p'(z)|^2 */
202 calc_mod_dp2(double z1, double z2, double p_coef[])
209 for(i=deg_p;i>=2;i--)
211 add(mp1,mp2,i*p_coef[(i-2)*2],i*p_coef[(i-2)*2+1]);
215 return mp1*mp1+mp2*mp2;
219 calc_all_p(euler2dstruct *sp)
222 double temp,p1,p2,z1,z2;
223 for(j=(sp->hide_vortex?sp->Nvortex:0);j<sp->N;j++) if(!sp->dead[j])
229 for(i=deg_p;i>=2;i--)
231 add(p1,p2,sp->p_coef[(i-2)*2],sp->p_coef[(i-2)*2+1]);
242 calc_all_mod_dp2(double *x, euler2dstruct *sp)
245 double temp,mp1,mp2,z1,z2;
246 for(j=0;j<sp->N;j++) if(!sp->dead[j])
252 for(i=deg_p;i>=2;i--)
254 add(mp1,mp2,i*sp->p_coef[(i-2)*2],i*sp->p_coef[(i-2)*2+1]);
258 sp->mod_dp2[j] = mp1*mp1+mp2*mp2;
263 derivs(double *x, euler2dstruct *sp)
266 double u1,u2,x1,x2,xij1,xij2,nxij;
269 if (variable_boundary)
270 calc_all_mod_dp2(sp->x,sp);
272 for (j=0;j<sp->Nvortex;j++) if (!sp->dead[j])
274 nx = x[2*j+0]*x[2*j+0] + x[2*j+1]*x[2*j+1];
276 sp->x_is_zero[j] = 1;
278 sp->x_is_zero[j] = 0;
279 sp->xs[2*j+0] = x[2*j+0]/nx;
280 sp->xs[2*j+1] = x[2*j+1]/nx;
284 (void) memset(sp->diffx,0,sizeof(double)*2*sp->N);
286 for (i=0;i<sp->N;i++) if (!sp->dead[i])
290 for (j=0;j<sp->Nvortex;j++) if (!sp->dead[j])
293 Calculate the Biot-Savart kernel, that is, effect of a
294 vortex point at a = (x[2*j+0],x[2*j+1]) at the point
295 x = (x1,x2), returning the vector field in (u1,u2).
297 In the plane, this is given by the formula
299 u = (x-a)/|x-a|^2 or zero if x=a.
301 However, in the unit disk we have to subtract from the
306 where as = a/|a|^2 is the reflection of a about the unit circle.
308 If however power != 1, then
310 u = (x-a)/|x-a|^(power+1) - |a|^(1-power) (x-as)/|x-as|^(power+1)
314 xij1 = x1 - x[2*j+0];
315 xij2 = x2 - x[2*j+1];
316 nxij = (power==1.0) ? xij1*xij1+xij2*xij2 : pow(xij1*xij1+xij2*xij2,(power+1)/2.0);
325 if (!sp->x_is_zero[j])
327 xij1 = x1 - sp->xs[2*j+0];
328 xij2 = x2 - sp->xs[2*j+1];
329 nxij = (power==1.0) ? xij1*xij1+xij2*xij2 : pow(xij1*xij1+xij2*xij2,(power+1)/2.0);
345 sp->diffx[2*i+0] += u1*sp->w[j];
346 sp->diffx[2*i+1] += u2*sp->w[j];
350 if (!sp->dead[i] && variable_boundary)
352 if (sp->mod_dp2[i] < 1e-5)
356 sp->diffx[2*i+0] /= sp->mod_dp2[i];
357 sp->diffx[2*i+1] /= sp->mod_dp2[i];
364 What perturb does is effectively
366 where k should be of order delta_t.
368 We have the option to do this more subtly by mapping points x
369 in the unit disk to points y in the plane, where y = f(|x|) x,
370 with f(t) = -log(1-t)/t.
372 This might reduce (but does not remove) problems where particles near
373 the edge of the boundary bounce around.
375 But it seems to be not that effective, so for now switch it off.
378 #define SUBTLE_PERTURB 0
381 perturb(double ret[], double x[], double k[], euler2dstruct *sp)
387 double d1,d2,t1,t2,mag,mag2,mlog1mmag,memmagdmag,xdotk;
388 for (i=0;i<sp->N;i++) if (!sp->dead[i])
394 mag2 = x1*x1 + x2*x2;
400 else if (mag2 > 1-1e-5)
405 mlog1mmag = -log(1-mag);
406 xdotk = x1*k1 + x2*k2;
407 t1 = (x1 + k1)*mlog1mmag/mag + x1*xdotk*(1.0/(1-mag)-mlog1mmag/mag)/mag/mag;
408 t2 = (x2 + k2)*mlog1mmag/mag + x2*xdotk*(1.0/(1-mag)-mlog1mmag/mag)/mag/mag;
409 mag = sqrt(t1*t1+t2*t2);
410 if (mag > 11.5 /* log(1e5) */)
414 memmagdmag = (mag>1e-5) ? ((1.0-exp(-mag))/mag) : (1-mag/2.0);
415 ret[2*i+0] = t1*memmagdmag;
416 ret[2*i+1] = t2*memmagdmag;
423 if (d1*d1+d2*d2 > 0.1)
430 for (i=0;i<sp->N;i++) if (!sp->dead[i])
436 if (k1*k1+k2*k2 > 0.1 || x1*x1+x2*x2 > 1-1e-5)
448 ode_solve(euler2dstruct *sp)
454 /* midpoint method */
456 (void) memcpy(sp->olddiffx,sp->diffx,sizeof(double)*2*sp->N);
457 for (i=0;i<sp->N;i++) if (!sp->dead[i]) {
458 sp->tempdiffx[2*i+0] = 0.5*delta_t*sp->diffx[2*i+0];
459 sp->tempdiffx[2*i+1] = 0.5*delta_t*sp->diffx[2*i+1];
461 perturb(sp->tempx,sp->x,sp->tempdiffx,sp);
462 derivs(sp->tempx,sp);
463 for (i=0;i<sp->N;i++) if (!sp->dead[i]) {
464 sp->tempdiffx[2*i+0] = delta_t*sp->diffx[2*i+0];
465 sp->tempdiffx[2*i+1] = delta_t*sp->diffx[2*i+1];
467 perturb(sp->x,sp->x,sp->tempdiffx,sp);
471 for (i=0;i<sp->N;i++) if (!sp->dead[i]) {
472 sp->tempdiffx[2*i+0] = delta_t*(1.5*sp->diffx[2*i+0] - 0.5*sp->olddiffx[2*i+0]);
473 sp->tempdiffx[2*i+1] = delta_t*(1.5*sp->diffx[2*i+1] - 0.5*sp->olddiffx[2*i+1]);
475 perturb(sp->x,sp->x,sp->tempdiffx,sp);
477 sp->olddiffx = sp->diffx;
482 #define deallocate(p,t) if (p!=NULL) {(void) free((void *) p); p=(t*)NULL; }
483 #define allocate(p,t,s) if ((p=(t*)malloc(sizeof(t)*s))==NULL)\
484 {free_euler2d(mi);return;}
487 free_euler2d(ModeInfo * mi)
489 euler2dstruct *sp = &euler2ds[MI_SCREEN(mi)];
490 deallocate(sp->csegs, XSegment);
491 deallocate(sp->old_segs, XSegment);
492 deallocate(sp->nold_segs, int);
493 deallocate(sp->lastx, short);
494 deallocate(sp->x, double);
495 deallocate(sp->diffx, double);
496 deallocate(sp->w, double);
497 deallocate(sp->olddiffx, double);
498 deallocate(sp->tempdiffx, double);
499 deallocate(sp->tempx, double);
500 deallocate(sp->dead, short);
501 deallocate(sp->boundary, XSegment);
502 deallocate(sp->xs, double);
503 deallocate(sp->x_is_zero, short);
504 deallocate(sp->p, double);
505 deallocate(sp->mod_dp2, double);
509 init_euler2d (ModeInfo * mi)
511 #define nr_rotates 18 /* how many rotations to try to fill as much of screen as possible - must be even number */
514 double r,theta,x,y,w;
515 double mag,xscale,yscale,p1,p2;
516 double low[nr_rotates],high[nr_rotates],pp1,pp2,pn1,pn2,angle1,angle2,tempangle,dist,scale,bestscale,temp;
519 if (power<0.5) power = 0.5;
520 if (power>3.0) power = 3.0;
521 variable_boundary &= power == 1.0;
523 if (power>1.0) delta_t *= pow(0.1,power-1);
525 MI_INIT (mi, euler2ds, free_euler2d);
526 sp = &euler2ds[MI_SCREEN(mi)];
529 jwxyz_XSetAntiAliasing (MI_DISPLAY(mi), MI_GC(mi), False);
532 sp->boundary_color = NRAND(MI_NPIXELS(mi));
533 sp->hide_vortex = NRAND(4) != 0;
537 sp->width = MI_WIDTH(mi);
538 sp->height = MI_HEIGHT(mi);
540 sp->N = MI_COUNT(mi)+number_of_vortex_points;
541 sp->Nvortex = number_of_vortex_points;
543 if (tail_len < 1) { /* minimum tail */
546 if (tail_len > MI_CYCLES(mi)) { /* maximum tail */
547 tail_len = MI_CYCLES(mi);
550 /* Clear the background. */
553 /* Allocate memory. */
555 if (sp->csegs == NULL) {
556 allocate(sp->csegs, XSegment, sp->N);
557 allocate(sp->old_segs, XSegment, sp->N * tail_len);
558 allocate(sp->nold_segs, int, tail_len);
559 allocate(sp->lastx, short, sp->N * 2);
560 allocate(sp->x, double, sp->N * 2);
561 allocate(sp->diffx, double, sp->N * 2);
562 allocate(sp->w, double, sp->Nvortex);
563 allocate(sp->olddiffx, double, sp->N * 2);
564 allocate(sp->tempdiffx, double, sp->N * 2);
565 allocate(sp->tempx, double, sp->N * 2);
566 allocate(sp->dead, short, sp->N);
567 allocate(sp->boundary, XSegment, n_bound_p);
568 allocate(sp->xs, double, sp->Nvortex * 2);
569 allocate(sp->x_is_zero, short, sp->Nvortex);
570 allocate(sp->p, double, sp->N * 2);
571 allocate(sp->mod_dp2, double, sp->N);
573 for (i=0;i<tail_len;i++) {
574 sp->nold_segs[i] = 0;
577 (void) memset(sp->dead,0,sp->N*sizeof(short));
579 if (variable_boundary)
581 /* Initialize polynomial p */
583 The polynomial p(z) = z + c_2 z^2 + ... c_n z^n needs to be
584 a bijection of the unit disk onto its image. This is achieved
585 by insisting that sum_{k=2}^n k |c_k| <= 1. Actually we set
586 the inequality to be equality (to get more interesting shapes).
589 for(k=2;k<=deg_p;k++)
591 r = positive_rand(1.0/k);
592 theta = balance_rand(2*M_PI);
593 sp->p_coef[2*(k-2)+0]=r*cos(theta);
594 sp->p_coef[2*(k-2)+1]=r*sin(theta);
597 if (mag > 0.0001) for(k=2;k<=deg_p;k++)
599 sp->p_coef[2*(k-2)+0] /= mag;
600 sp->p_coef[2*(k-2)+1] /= mag;
603 #if DEBUG_POINTED_REGION
604 for(k=2;k<=deg_p;k++){
605 sp->p_coef[2*(k-2)+0]=0;
606 sp->p_coef[2*(k-2)+1]=0;
608 sp->p_coef[2*(6-2)+0] = 1.0/6.0;
612 /* Here we figure out the best rotation of the domain so that it fills as
613 much of the screen as possible. The number of angles we look at is determined
614 by nr_rotates (we look every 180/nr_rotates degrees).
615 While we figure out the best angle to rotate, we also figure out the correct scaling factors.
618 for(k=0;k<nr_rotates;k++) {
623 for(k=0;k<n_bound_p;k++) {
624 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);
625 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);
626 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);
627 angle1 = nr_rotates/M_PI*atan2(p2-pp2,p1-pp1)-nr_rotates/2;
628 angle2 = nr_rotates/M_PI*atan2(pn2-p2,pn1-p1)-nr_rotates/2;
629 while (angle1<0) angle1+=nr_rotates*2;
630 while (angle2<0) angle2+=nr_rotates*2;
631 if (angle1>nr_rotates*1.75 && angle2<nr_rotates*0.25) angle2+=nr_rotates*2;
632 if (angle1<nr_rotates*0.25 && angle2>nr_rotates*1.75) angle1+=nr_rotates*2;
638 for(i=(int)floor(angle1);i<(int)ceil(angle2);i++) {
639 dist = cos((double)i*M_PI/nr_rotates)*p1 + sin((double)i*M_PI/nr_rotates)*p2;
640 if (i%(nr_rotates*2)<nr_rotates) {
641 if (dist>high[i%nr_rotates]) high[i%nr_rotates] = dist;
642 if (dist<low[i%nr_rotates]) low[i%nr_rotates] = dist;
645 if (-dist>high[i%nr_rotates]) high[i%nr_rotates] = -dist;
646 if (-dist<low[i%nr_rotates]) low[i%nr_rotates] = -dist;
651 for (i=0;i<nr_rotates;i++) {
652 xscale = (sp->width-5.0)/(high[i]-low[i]);
653 yscale = (sp->height-5.0)/(high[(i+nr_rotates/2)%nr_rotates]-low[(i+nr_rotates/2)%nr_rotates]);
654 scale = (xscale>yscale) ? yscale : xscale;
655 if (scale>bestscale) {
660 /* Here we do the rotation. The way we do this is to replace the
661 polynomial p(z) by a^{-1} p(a z) where a = exp(i best_angle).
665 for(k=2;k<=deg_p;k++)
667 mult(p1,p2,cos((double)besti*M_PI/nr_rotates),sin((double)besti*M_PI/nr_rotates));
668 mult(sp->p_coef[2*(k-2)+0],sp->p_coef[2*(k-2)+1],p1,p2);
671 sp->scale = bestscale;
672 sp->xshift = -(low[besti]+high[besti])/2.0*sp->scale+sp->width/2;
673 if (besti<nr_rotates/2)
674 sp->yshift = -(low[besti+nr_rotates/2]+high[besti+nr_rotates/2])/2.0*sp->scale+sp->height/2;
676 sp->yshift = (low[besti-nr_rotates/2]+high[besti-nr_rotates/2])/2.0*sp->scale+sp->height/2;
679 /* Initialize boundary */
681 for(k=0;k<n_bound_p;k++)
684 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);
685 sp->boundary[k].x1 = (short)(p1*sp->scale+sp->xshift);
686 sp->boundary[k].y1 = (short)(p2*sp->scale+sp->yshift);
688 for(k=1;k<n_bound_p;k++)
690 sp->boundary[k].x2 = sp->boundary[k-1].x1;
691 sp->boundary[k].y2 = sp->boundary[k-1].y1;
693 sp->boundary[0].x2 = sp->boundary[n_bound_p-1].x1;
694 sp->boundary[0].y2 = sp->boundary[n_bound_p-1].y1;
698 if (sp->width>sp->height)
699 sp->radius = sp->height/2.0-5.0;
701 sp->radius = sp->width/2.0-5.0;
704 /* Initialize point positions */
706 for (i=sp->Nvortex;i<sp->N;i++) {
708 r = sqrt(positive_rand(1.0));
709 theta = balance_rand(2*M_PI);
710 sp->x[2*i+0]=r*cos(theta);
711 sp->x[2*i+1]=r*sin(theta);
712 /* This is to make sure the initial distribution of points is uniform */
713 } while (variable_boundary &&
714 calc_mod_dp2(sp->x[2*i+0],sp->x[2*i+1],sp->p_coef)
719 /* number of vortex points with negative vorticity */
724 /* if n is even make sure that np==n/2 is twice as likely
725 as the other possibilities. */
731 r = sqrt(positive_rand(0.77));
732 theta = balance_rand(2*M_PI);
735 r = 0.02+positive_rand(0.1);
736 w = (2*(k<np)-1)*2.0/sp->Nvortex;
737 for (i=sp->Nvortex*k/n;i<sp->Nvortex*(k+1)/n;i++) {
738 theta = balance_rand(2*M_PI);
739 sp->x[2*i+0]=x + r*cos(theta);
740 sp->x[2*i+1]=y + r*sin(theta);
747 draw_euler2d (ModeInfo * mi)
749 Display *display = MI_DISPLAY(mi);
750 Window window = MI_WINDOW(mi);
752 int b, col, n_non_vortex_segs;
755 MI_IS_DRAWN(mi) = True;
757 if (euler2ds == NULL)
759 sp = &euler2ds[MI_SCREEN(mi)];
760 if (sp->csegs == NULL)
764 if (variable_boundary)
768 for(b=sp->Nvortex;b<sp->N;b++) if(!sp->dead[b])
770 sp->csegs[sp->cnsegs].x1 = sp->lastx[2*b+0];
771 sp->csegs[sp->cnsegs].y1 = sp->lastx[2*b+1];
772 if (variable_boundary)
774 sp->csegs[sp->cnsegs].x2 = (short)(sp->p[2*b+0]*sp->scale+sp->xshift);
775 sp->csegs[sp->cnsegs].y2 = (short)(sp->p[2*b+1]*sp->scale+sp->yshift);
779 sp->csegs[sp->cnsegs].x2 = (short)(sp->x[2*b+0]*sp->radius+sp->width/2);
780 sp->csegs[sp->cnsegs].y2 = (short)(sp->x[2*b+1]*sp->radius+sp->height/2);
782 sp->lastx[2*b+0] = sp->csegs[sp->cnsegs].x2;
783 sp->lastx[2*b+1] = sp->csegs[sp->cnsegs].y2;
786 n_non_vortex_segs = sp->cnsegs;
788 if (!sp->hide_vortex) for(b=0;b<sp->Nvortex;b++) if(!sp->dead[b])
790 sp->csegs[sp->cnsegs].x1 = sp->lastx[2*b+0];
791 sp->csegs[sp->cnsegs].y1 = sp->lastx[2*b+1];
792 if (variable_boundary)
794 sp->csegs[sp->cnsegs].x2 = (short)(sp->p[2*b+0]*sp->scale+sp->xshift);
795 sp->csegs[sp->cnsegs].y2 = (short)(sp->p[2*b+1]*sp->scale+sp->yshift);
799 sp->csegs[sp->cnsegs].x2 = (short)(sp->x[2*b+0]*sp->radius+sp->width/2);
800 sp->csegs[sp->cnsegs].y2 = (short)(sp->x[2*b+1]*sp->radius+sp->height/2);
802 sp->lastx[2*b+0] = sp->csegs[sp->cnsegs].x2;
803 sp->lastx[2*b+1] = sp->csegs[sp->cnsegs].y2;
808 XSetForeground(display, gc, MI_BLACK_PIXEL(mi));
810 XDrawSegments(display, window, gc, sp->old_segs+sp->c_old_seg*sp->N, sp->nold_segs[sp->c_old_seg]);
812 if (MI_NPIXELS(mi) > 2){ /* render colour */
813 for (col = 0; col < MI_NPIXELS(mi); col++) {
814 int start = col*n_non_vortex_segs/MI_NPIXELS(mi);
815 int finish = (col+1)*n_non_vortex_segs/MI_NPIXELS(mi);
816 XSetForeground(display, gc, MI_PIXEL(mi, col));
817 XDrawSegments(display, window, gc,sp->csegs+start, finish-start);
819 if (!sp->hide_vortex) {
820 XSetForeground(display, gc, MI_WHITE_PIXEL(mi));
821 XDrawSegments(display, window, gc,sp->csegs+n_non_vortex_segs, sp->cnsegs-n_non_vortex_segs);
824 } else { /* render mono */
825 XSetForeground(display, gc, MI_WHITE_PIXEL(mi));
826 XDrawSegments(display, window, gc,
827 sp->csegs, sp->cnsegs);
830 if (MI_NPIXELS(mi) > 2) /* render colour */
831 XSetForeground(display, gc, MI_PIXEL(mi, sp->boundary_color));
833 XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_WHITE_PIXEL(mi));
834 if (variable_boundary)
835 XDrawSegments(display, window, gc,
836 sp->boundary, n_bound_p);
838 XDrawArc(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
839 sp->width/2 - (int) sp->radius - 1, sp->height/2 - (int) sp->radius -1,
840 (int) (2*sp->radius) + 2, (int) (2* sp->radius) + 2, 0, 64*360);
842 /* Copy to erase-list */
843 (void) memcpy(sp->old_segs+sp->c_old_seg*sp->N, sp->csegs, sp->cnsegs*sizeof(XSegment));
844 sp->nold_segs[sp->c_old_seg] = sp->cnsegs;
846 if (sp->c_old_seg >= tail_len)
850 if (++sp->count > MI_CYCLES(mi)) /* pick a new flow */
856 reshape_euler2d(ModeInfo * mi, int width, int height)
858 XClearWindow (MI_DISPLAY (mi), MI_WINDOW(mi));
863 refresh_euler2d (ModeInfo * mi)
869 euler2d_handle_event (ModeInfo *mi, XEvent *event)
871 if (screenhack_event_helper (MI_DISPLAY(mi), MI_WINDOW(mi), event))
881 XSCREENSAVER_MODULE ("Euler2D", euler2d)
883 #endif /* MODE_euler2d */