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 {"-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 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", "release_euler2d",
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(sp);return;}
487 free_euler2d(euler2dstruct *sp)
489 deallocate(sp->csegs, XSegment);
490 deallocate(sp->old_segs, XSegment);
491 deallocate(sp->nold_segs, int);
492 deallocate(sp->lastx, short);
493 deallocate(sp->x, double);
494 deallocate(sp->diffx, double);
495 deallocate(sp->w, double);
496 deallocate(sp->olddiffx, double);
497 deallocate(sp->tempdiffx, double);
498 deallocate(sp->tempx, double);
499 deallocate(sp->dead, short);
500 deallocate(sp->boundary, XSegment);
501 deallocate(sp->xs, double);
502 deallocate(sp->x_is_zero, short);
503 deallocate(sp->p, double);
504 deallocate(sp->mod_dp2, double);
508 init_euler2d(ModeInfo * mi)
510 #define nr_rotates 18 /* how many rotations to try to fill as much of screen as possible - must be even number */
513 double r,theta,x,y,w;
514 double mag,xscale,yscale,p1,p2;
515 double low[nr_rotates],high[nr_rotates],pp1,pp2,pn1,pn2,angle1,angle2,tempangle,dist,scale,bestscale,temp;
518 if (power<0.5) power = 0.5;
519 if (power>3.0) power = 3.0;
520 variable_boundary &= power == 1.0;
522 if (power>1.0) delta_t *= pow(0.1,power-1);
524 if (euler2ds == NULL) {
525 if ((euler2ds = (euler2dstruct *) calloc(MI_NUM_SCREENS(mi),
526 sizeof (euler2dstruct))) == NULL)
529 sp = &euler2ds[MI_SCREEN(mi)];
531 sp->boundary_color = NRAND(MI_NPIXELS(mi));
532 sp->hide_vortex = NRAND(4) != 0;
536 sp->width = MI_WIDTH(mi);
537 sp->height = MI_HEIGHT(mi);
539 sp->N = MI_COUNT(mi)+number_of_vortex_points;
540 sp->Nvortex = number_of_vortex_points;
542 if (tail_len < 1) { /* minimum tail */
545 if (tail_len > MI_CYCLES(mi)) { /* maximum tail */
546 tail_len = MI_CYCLES(mi);
549 /* Clear the background. */
554 /* Allocate memory. */
556 if (sp->csegs == NULL) {
557 allocate(sp->csegs, XSegment, sp->N);
558 allocate(sp->old_segs, XSegment, sp->N * tail_len);
559 allocate(sp->nold_segs, int, tail_len);
560 allocate(sp->lastx, short, sp->N * 2);
561 allocate(sp->x, double, sp->N * 2);
562 allocate(sp->diffx, double, sp->N * 2);
563 allocate(sp->w, double, sp->Nvortex);
564 allocate(sp->olddiffx, double, sp->N * 2);
565 allocate(sp->tempdiffx, double, sp->N * 2);
566 allocate(sp->tempx, double, sp->N * 2);
567 allocate(sp->dead, short, sp->N);
568 allocate(sp->boundary, XSegment, n_bound_p);
569 allocate(sp->xs, double, sp->Nvortex * 2);
570 allocate(sp->x_is_zero, short, sp->Nvortex);
571 allocate(sp->p, double, sp->N * 2);
572 allocate(sp->mod_dp2, double, sp->N);
574 for (i=0;i<tail_len;i++) {
575 sp->nold_segs[i] = 0;
578 (void) memset(sp->dead,0,sp->N*sizeof(short));
580 if (variable_boundary)
582 /* Initialize polynomial p */
584 The polynomial p(z) = z + c_2 z^2 + ... c_n z^n needs to be
585 a bijection of the unit disk onto its image. This is achieved
586 by insisting that sum_{k=2}^n k |c_k| <= 1. Actually we set
587 the inequality to be equality (to get more interesting shapes).
590 for(k=2;k<=deg_p;k++)
592 r = positive_rand(1.0/k);
593 theta = balance_rand(2*M_PI);
594 sp->p_coef[2*(k-2)+0]=r*cos(theta);
595 sp->p_coef[2*(k-2)+1]=r*sin(theta);
598 if (mag > 0.0001) for(k=2;k<=deg_p;k++)
600 sp->p_coef[2*(k-2)+0] /= mag;
601 sp->p_coef[2*(k-2)+1] /= mag;
604 #if DEBUG_POINTED_REGION
605 for(k=2;k<=deg_p;k++){
606 sp->p_coef[2*(k-2)+0]=0;
607 sp->p_coef[2*(k-2)+1]=0;
609 sp->p_coef[2*(6-2)+0] = 1.0/6.0;
613 /* Here we figure out the best rotation of the domain so that it fills as
614 much of the screen as possible. The number of angles we look at is determined
615 by nr_rotates (we look every 180/nr_rotates degrees).
616 While we figure out the best angle to rotate, we also figure out the correct scaling factors.
619 for(k=0;k<nr_rotates;k++) {
624 for(k=0;k<n_bound_p;k++) {
625 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);
626 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);
627 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);
628 angle1 = nr_rotates/M_PI*atan2(p2-pp2,p1-pp1)-nr_rotates/2;
629 angle2 = nr_rotates/M_PI*atan2(pn2-p2,pn1-p1)-nr_rotates/2;
630 while (angle1<0) angle1+=nr_rotates*2;
631 while (angle2<0) angle2+=nr_rotates*2;
632 if (angle1>nr_rotates*1.75 && angle2<nr_rotates*0.25) angle2+=nr_rotates*2;
633 if (angle1<nr_rotates*0.25 && angle2>nr_rotates*1.75) angle1+=nr_rotates*2;
639 for(i=(int)floor(angle1);i<(int)ceil(angle2);i++) {
640 dist = cos((double)i*M_PI/nr_rotates)*p1 + sin((double)i*M_PI/nr_rotates)*p2;
641 if (i%(nr_rotates*2)<nr_rotates) {
642 if (dist>high[i%nr_rotates]) high[i%nr_rotates] = dist;
643 if (dist<low[i%nr_rotates]) low[i%nr_rotates] = dist;
646 if (-dist>high[i%nr_rotates]) high[i%nr_rotates] = -dist;
647 if (-dist<low[i%nr_rotates]) low[i%nr_rotates] = -dist;
652 for (i=0;i<nr_rotates;i++) {
653 xscale = (sp->width-5.0)/(high[i]-low[i]);
654 yscale = (sp->height-5.0)/(high[(i+nr_rotates/2)%nr_rotates]-low[(i+nr_rotates/2)%nr_rotates]);
655 scale = (xscale>yscale) ? yscale : xscale;
656 if (scale>bestscale) {
661 /* Here we do the rotation. The way we do this is to replace the
662 polynomial p(z) by a^{-1} p(a z) where a = exp(i best_angle).
666 for(k=2;k<=deg_p;k++)
668 mult(p1,p2,cos((double)besti*M_PI/nr_rotates),sin((double)besti*M_PI/nr_rotates));
669 mult(sp->p_coef[2*(k-2)+0],sp->p_coef[2*(k-2)+1],p1,p2);
672 sp->scale = bestscale;
673 sp->xshift = -(low[besti]+high[besti])/2.0*sp->scale+sp->width/2;
674 if (besti<nr_rotates/2)
675 sp->yshift = -(low[besti+nr_rotates/2]+high[besti+nr_rotates/2])/2.0*sp->scale+sp->height/2;
677 sp->yshift = (low[besti-nr_rotates/2]+high[besti-nr_rotates/2])/2.0*sp->scale+sp->height/2;
680 /* Initialize boundary */
682 for(k=0;k<n_bound_p;k++)
685 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);
686 sp->boundary[k].x1 = (short)(p1*sp->scale+sp->xshift);
687 sp->boundary[k].y1 = (short)(p2*sp->scale+sp->yshift);
689 for(k=1;k<n_bound_p;k++)
691 sp->boundary[k].x2 = sp->boundary[k-1].x1;
692 sp->boundary[k].y2 = sp->boundary[k-1].y1;
694 sp->boundary[0].x2 = sp->boundary[n_bound_p-1].x1;
695 sp->boundary[0].y2 = sp->boundary[n_bound_p-1].y1;
699 if (sp->width>sp->height)
700 sp->radius = sp->height/2.0-5.0;
702 sp->radius = sp->width/2.0-5.0;
705 /* Initialize point positions */
707 for (i=sp->Nvortex;i<sp->N;i++) {
709 r = sqrt(positive_rand(1.0));
710 theta = balance_rand(2*M_PI);
711 sp->x[2*i+0]=r*cos(theta);
712 sp->x[2*i+1]=r*sin(theta);
713 /* This is to make sure the initial distribution of points is uniform */
714 } while (variable_boundary &&
715 calc_mod_dp2(sp->x[2*i+0],sp->x[2*i+1],sp->p_coef)
720 /* number of vortex points with negative vorticity */
725 /* if n is even make sure that np==n/2 is twice as likely
726 as the other possibilities. */
732 r = sqrt(positive_rand(0.77));
733 theta = balance_rand(2*M_PI);
736 r = 0.02+positive_rand(0.1);
737 w = (2*(k<np)-1)*2.0/sp->Nvortex;
738 for (i=sp->Nvortex*k/n;i<sp->Nvortex*(k+1)/n;i++) {
739 theta = balance_rand(2*M_PI);
740 sp->x[2*i+0]=x + r*cos(theta);
741 sp->x[2*i+1]=y + r*sin(theta);
748 draw_euler2d(ModeInfo * mi)
750 Display *display = MI_DISPLAY(mi);
751 Window window = MI_WINDOW(mi);
753 int b, col, n_non_vortex_segs;
756 MI_IS_DRAWN(mi) = True;
758 if (euler2ds == NULL)
760 sp = &euler2ds[MI_SCREEN(mi)];
761 if (sp->csegs == NULL)
765 if (variable_boundary)
769 for(b=sp->Nvortex;b<sp->N;b++) if(!sp->dead[b])
771 sp->csegs[sp->cnsegs].x1 = sp->lastx[2*b+0];
772 sp->csegs[sp->cnsegs].y1 = sp->lastx[2*b+1];
773 if (variable_boundary)
775 sp->csegs[sp->cnsegs].x2 = (short)(sp->p[2*b+0]*sp->scale+sp->xshift);
776 sp->csegs[sp->cnsegs].y2 = (short)(sp->p[2*b+1]*sp->scale+sp->yshift);
780 sp->csegs[sp->cnsegs].x2 = (short)(sp->x[2*b+0]*sp->radius+sp->width/2);
781 sp->csegs[sp->cnsegs].y2 = (short)(sp->x[2*b+1]*sp->radius+sp->height/2);
783 sp->lastx[2*b+0] = sp->csegs[sp->cnsegs].x2;
784 sp->lastx[2*b+1] = sp->csegs[sp->cnsegs].y2;
787 n_non_vortex_segs = sp->cnsegs;
789 if (!sp->hide_vortex) for(b=0;b<sp->Nvortex;b++) if(!sp->dead[b])
791 sp->csegs[sp->cnsegs].x1 = sp->lastx[2*b+0];
792 sp->csegs[sp->cnsegs].y1 = sp->lastx[2*b+1];
793 if (variable_boundary)
795 sp->csegs[sp->cnsegs].x2 = (short)(sp->p[2*b+0]*sp->scale+sp->xshift);
796 sp->csegs[sp->cnsegs].y2 = (short)(sp->p[2*b+1]*sp->scale+sp->yshift);
800 sp->csegs[sp->cnsegs].x2 = (short)(sp->x[2*b+0]*sp->radius+sp->width/2);
801 sp->csegs[sp->cnsegs].y2 = (short)(sp->x[2*b+1]*sp->radius+sp->height/2);
803 sp->lastx[2*b+0] = sp->csegs[sp->cnsegs].x2;
804 sp->lastx[2*b+1] = sp->csegs[sp->cnsegs].y2;
809 XSetForeground(display, gc, MI_BLACK_PIXEL(mi));
811 XDrawSegments(display, window, gc, sp->old_segs+sp->c_old_seg*sp->N, sp->nold_segs[sp->c_old_seg]);
813 if (MI_NPIXELS(mi) > 2){ /* render colour */
814 for (col = 0; col < MI_NPIXELS(mi); col++) {
815 int start = col*n_non_vortex_segs/MI_NPIXELS(mi);
816 int finish = (col+1)*n_non_vortex_segs/MI_NPIXELS(mi);
817 XSetForeground(display, gc, MI_PIXEL(mi, col));
818 XDrawSegments(display, window, gc,sp->csegs+start, finish-start);
820 if (!sp->hide_vortex) {
821 XSetForeground(display, gc, MI_WHITE_PIXEL(mi));
822 XDrawSegments(display, window, gc,sp->csegs+n_non_vortex_segs, sp->cnsegs-n_non_vortex_segs);
825 } else { /* render mono */
826 XSetForeground(display, gc, MI_WHITE_PIXEL(mi));
827 XDrawSegments(display, window, gc,
828 sp->csegs, sp->cnsegs);
831 if (MI_NPIXELS(mi) > 2) /* render colour */
832 XSetForeground(display, gc, MI_PIXEL(mi, sp->boundary_color));
834 XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_WHITE_PIXEL(mi));
835 if (variable_boundary)
836 XDrawSegments(display, window, gc,
837 sp->boundary, n_bound_p);
839 XDrawArc(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
840 sp->width/2 - (int) sp->radius - 1, sp->height/2 - (int) sp->radius -1,
841 (int) (2*sp->radius) + 2, (int) (2* sp->radius) + 2, 0, 64*360);
843 /* Copy to erase-list */
844 (void) memcpy(sp->old_segs+sp->c_old_seg*sp->N, sp->csegs, sp->cnsegs*sizeof(XSegment));
845 sp->nold_segs[sp->c_old_seg] = sp->cnsegs;
847 if (sp->c_old_seg >= tail_len)
851 if (++sp->count > MI_CYCLES(mi)) /* pick a new flow */
857 release_euler2d(ModeInfo * mi)
859 if (euler2ds != NULL) {
862 for (screen = 0; screen < MI_NUM_SCREENS(mi); screen++)
863 free_euler2d(&euler2ds[screen]);
864 (void) free((void *) euler2ds);
865 euler2ds = (euler2dstruct *) NULL;
870 refresh_euler2d(ModeInfo * mi)
875 #endif /* MODE_euler2d */