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 # define reshape_euler2d 0
50 # define euler2d_handle_event 0
51 # define SMOOTH_COLORS
52 # include "xlockmore.h" /* in xscreensaver distribution */
53 #else /* STANDALONE */
54 # include "xlock.h" /* in xlockmore distribution */
55 #endif /* STANDALONE */
59 #define DEF_EULERTAIL "10"
61 #define DEBUG_POINTED_REGION 0
64 static int variable_boundary = 1;
65 static float power = 1;
67 static XrmOptionDescRec opts[] =
69 {"-eulertail", ".euler2d.eulertail", XrmoptionSepArg, NULL},
70 {"-eulerpower", ".euler2d.eulerpower", XrmoptionSepArg, NULL},
72 static argtype vars[] =
74 {&tail_len, "eulertail",
75 "EulerTail", (char *) DEF_EULERTAIL, t_Int},
76 {&power, "eulerpower",
77 "EulerPower", "1", t_Float},
79 static OptionStruct desc[] =
81 {"-eulertail len", "Length of Euler2d tails"},
82 {"-eulerpower power", "power of interaction law for points for Euler2d"},
85 ENTRYPOINT ModeSpecOpt euler2d_opts =
86 {sizeof opts / sizeof opts[0], opts,
87 sizeof vars / sizeof vars[0], vars, desc};
90 ModStruct euler2d_description = {
91 "euler2d", "init_euler2d", "draw_euler2d", "release_euler2d",
92 "refresh_euler2d", "init_euler2d", (char *) NULL, &euler2d_opts,
93 1000, 1024, 3000, 1, 64, 1.0, "",
94 "Simulates 2D incompressible invisid fluid.", 0, NULL
99 #define balance_rand(v) ((LRAND()/MAXRAND*(v))-((v)/2)) /* random around 0 */
100 #define positive_rand(v) (LRAND()/MAXRAND*(v)) /* positive random */
102 #define number_of_vortex_points 20
104 #define n_bound_p 500
107 static double delta_t;
113 double xshift,yshift,scale;
119 /* x[2i+0] = x coord for nth point
120 x[2i+1] = y coord for nth point
121 w[i] = vorticity at nth point
130 /* (xs[2i+0],xs[2i+1]) is reflection of (x[2i+0],x[2i+1]) about unit circle
131 xs[2i+0] = x[2i+0]/nx
132 xs[2i+1] = x[2i+1]/nx
134 nx = x[2i+0]*x[2i+0] + x[2i+1]*x[2i+1]
136 x_is_zero[i] = (nx < 1e-10)
141 /* (p[2i+0],p[2i+1]) is image of (x[2i+0],x[2i+1]) under polynomial p.
142 mod_dp2 is |p'(z)|^2 when z = (x[2i+0],x[2i+1]).
147 /* Sometimes in our calculations we get overflow or numbers that are too big.
148 If that happens with the point x[2*i+0], x[2*i+1], we set dead[i].
161 double p_coef[2*(deg_p-1)];
166 static euler2dstruct *euler2ds = (euler2dstruct *) NULL;
169 If variable_boundary == 1, then we make a variable boundary.
170 The way this is done is to map the unit disk under a
172 p(z) = z + c_2 z^2 + ... + c_n z^n
173 where n = deg_p. sp->p_coef contains the complex numbers
177 #define add(a1,a2,b1,b2) (a1)+=(b1);(a2)+=(b2)
178 #define mult(a1,a2,b1,b2) temp=(a1)*(b1)-(a2)*(b2); \
179 (a2)=(a1)*(b2)+(a2)*(b1);(a1)=temp
182 calc_p(double *p1, double *p2, double z1, double z2, double p_coef[])
189 for(i=deg_p;i>=2;i--)
191 add(*p1,*p2,p_coef[(i-2)*2],p_coef[(i-2)*2+1]);
198 /* Calculate |p'(z)|^2 */
200 calc_mod_dp2(double z1, double z2, double p_coef[])
207 for(i=deg_p;i>=2;i--)
209 add(mp1,mp2,i*p_coef[(i-2)*2],i*p_coef[(i-2)*2+1]);
213 return mp1*mp1+mp2*mp2;
217 calc_all_p(euler2dstruct *sp)
220 double temp,p1,p2,z1,z2;
221 for(j=(sp->hide_vortex?sp->Nvortex:0);j<sp->N;j++) if(!sp->dead[j])
227 for(i=deg_p;i>=2;i--)
229 add(p1,p2,sp->p_coef[(i-2)*2],sp->p_coef[(i-2)*2+1]);
240 calc_all_mod_dp2(double *x, euler2dstruct *sp)
243 double temp,mp1,mp2,z1,z2;
244 for(j=0;j<sp->N;j++) if(!sp->dead[j])
250 for(i=deg_p;i>=2;i--)
252 add(mp1,mp2,i*sp->p_coef[(i-2)*2],i*sp->p_coef[(i-2)*2+1]);
256 sp->mod_dp2[j] = mp1*mp1+mp2*mp2;
261 derivs(double *x, euler2dstruct *sp)
264 double u1,u2,x1,x2,xij1,xij2,nxij;
267 if (variable_boundary)
268 calc_all_mod_dp2(sp->x,sp);
270 for (j=0;j<sp->Nvortex;j++) if (!sp->dead[j])
272 nx = x[2*j+0]*x[2*j+0] + x[2*j+1]*x[2*j+1];
274 sp->x_is_zero[j] = 1;
276 sp->x_is_zero[j] = 0;
277 sp->xs[2*j+0] = x[2*j+0]/nx;
278 sp->xs[2*j+1] = x[2*j+1]/nx;
282 (void) memset(sp->diffx,0,sizeof(double)*2*sp->N);
284 for (i=0;i<sp->N;i++) if (!sp->dead[i])
288 for (j=0;j<sp->Nvortex;j++) if (!sp->dead[j])
291 Calculate the Biot-Savart kernel, that is, effect of a
292 vortex point at a = (x[2*j+0],x[2*j+1]) at the point
293 x = (x1,x2), returning the vector field in (u1,u2).
295 In the plane, this is given by the formula
297 u = (x-a)/|x-a|^2 or zero if x=a.
299 However, in the unit disk we have to subtract from the
304 where as = a/|a|^2 is the reflection of a about the unit circle.
306 If however power != 1, then
308 u = (x-a)/|x-a|^(power+1) - |a|^(1-power) (x-as)/|x-as|^(power+1)
312 xij1 = x1 - x[2*j+0];
313 xij2 = x2 - x[2*j+1];
314 nxij = (power==1.0) ? xij1*xij1+xij2*xij2 : pow(xij1*xij1+xij2*xij2,(power+1)/2.0);
323 if (!sp->x_is_zero[j])
325 xij1 = x1 - sp->xs[2*j+0];
326 xij2 = x2 - sp->xs[2*j+1];
327 nxij = (power==1.0) ? xij1*xij1+xij2*xij2 : pow(xij1*xij1+xij2*xij2,(power+1)/2.0);
343 sp->diffx[2*i+0] += u1*sp->w[j];
344 sp->diffx[2*i+1] += u2*sp->w[j];
348 if (!sp->dead[i] && variable_boundary)
350 if (sp->mod_dp2[i] < 1e-5)
354 sp->diffx[2*i+0] /= sp->mod_dp2[i];
355 sp->diffx[2*i+1] /= sp->mod_dp2[i];
362 What perturb does is effectively
364 where k should be of order delta_t.
366 We have the option to do this more subtly by mapping points x
367 in the unit disk to points y in the plane, where y = f(|x|) x,
368 with f(t) = -log(1-t)/t.
370 This might reduce (but does not remove) problems where particles near
371 the edge of the boundary bounce around.
373 But it seems to be not that effective, so for now switch it off.
376 #define SUBTLE_PERTURB 0
379 perturb(double ret[], double x[], double k[], euler2dstruct *sp)
385 double d1,d2,t1,t2,mag,mag2,mlog1mmag,memmagdmag,xdotk;
386 for (i=0;i<sp->N;i++) if (!sp->dead[i])
392 mag2 = x1*x1 + x2*x2;
398 else if (mag2 > 1-1e-5)
403 mlog1mmag = -log(1-mag);
404 xdotk = x1*k1 + x2*k2;
405 t1 = (x1 + k1)*mlog1mmag/mag + x1*xdotk*(1.0/(1-mag)-mlog1mmag/mag)/mag/mag;
406 t2 = (x2 + k2)*mlog1mmag/mag + x2*xdotk*(1.0/(1-mag)-mlog1mmag/mag)/mag/mag;
407 mag = sqrt(t1*t1+t2*t2);
408 if (mag > 11.5 /* log(1e5) */)
412 memmagdmag = (mag>1e-5) ? ((1.0-exp(-mag))/mag) : (1-mag/2.0);
413 ret[2*i+0] = t1*memmagdmag;
414 ret[2*i+1] = t2*memmagdmag;
421 if (d1*d1+d2*d2 > 0.1)
428 for (i=0;i<sp->N;i++) if (!sp->dead[i])
434 if (k1*k1+k2*k2 > 0.1 || x1*x1+x2*x2 > 1-1e-5)
446 ode_solve(euler2dstruct *sp)
452 /* midpoint method */
454 (void) memcpy(sp->olddiffx,sp->diffx,sizeof(double)*2*sp->N);
455 for (i=0;i<sp->N;i++) if (!sp->dead[i]) {
456 sp->tempdiffx[2*i+0] = 0.5*delta_t*sp->diffx[2*i+0];
457 sp->tempdiffx[2*i+1] = 0.5*delta_t*sp->diffx[2*i+1];
459 perturb(sp->tempx,sp->x,sp->tempdiffx,sp);
460 derivs(sp->tempx,sp);
461 for (i=0;i<sp->N;i++) if (!sp->dead[i]) {
462 sp->tempdiffx[2*i+0] = delta_t*sp->diffx[2*i+0];
463 sp->tempdiffx[2*i+1] = delta_t*sp->diffx[2*i+1];
465 perturb(sp->x,sp->x,sp->tempdiffx,sp);
469 for (i=0;i<sp->N;i++) if (!sp->dead[i]) {
470 sp->tempdiffx[2*i+0] = delta_t*(1.5*sp->diffx[2*i+0] - 0.5*sp->olddiffx[2*i+0]);
471 sp->tempdiffx[2*i+1] = delta_t*(1.5*sp->diffx[2*i+1] - 0.5*sp->olddiffx[2*i+1]);
473 perturb(sp->x,sp->x,sp->tempdiffx,sp);
475 sp->olddiffx = sp->diffx;
480 #define deallocate(p,t) if (p!=NULL) {(void) free((void *) p); p=(t*)NULL; }
481 #define allocate(p,t,s) if ((p=(t*)malloc(sizeof(t)*s))==NULL)\
482 {free_euler2d(sp);return;}
485 free_euler2d(euler2dstruct *sp)
487 deallocate(sp->csegs, XSegment);
488 deallocate(sp->old_segs, XSegment);
489 deallocate(sp->nold_segs, int);
490 deallocate(sp->lastx, short);
491 deallocate(sp->x, double);
492 deallocate(sp->diffx, double);
493 deallocate(sp->w, double);
494 deallocate(sp->olddiffx, double);
495 deallocate(sp->tempdiffx, double);
496 deallocate(sp->tempx, double);
497 deallocate(sp->dead, short);
498 deallocate(sp->boundary, XSegment);
499 deallocate(sp->xs, double);
500 deallocate(sp->x_is_zero, short);
501 deallocate(sp->p, double);
502 deallocate(sp->mod_dp2, double);
506 init_euler2d (ModeInfo * mi)
508 #define nr_rotates 18 /* how many rotations to try to fill as much of screen as possible - must be even number */
511 double r,theta,x,y,w;
512 double mag,xscale,yscale,p1,p2;
513 double low[nr_rotates],high[nr_rotates],pp1,pp2,pn1,pn2,angle1,angle2,tempangle,dist,scale,bestscale,temp;
516 if (power<0.5) power = 0.5;
517 if (power>3.0) power = 3.0;
518 variable_boundary &= power == 1.0;
520 if (power>1.0) delta_t *= pow(0.1,power-1);
522 if (euler2ds == NULL) {
523 if ((euler2ds = (euler2dstruct *) calloc(MI_NUM_SCREENS(mi),
524 sizeof (euler2dstruct))) == NULL)
527 sp = &euler2ds[MI_SCREEN(mi)];
530 jwxyz_XSetAntiAliasing (MI_DISPLAY(mi), MI_GC(mi), False);
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 XSCREENSAVER_MODULE ("Euler2d", euler2d)
879 #endif /* MODE_euler2d */