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 # include "xlockmore.h" /* in xscreensaver distribution */
54 #else /* STANDALONE */
55 # include "xlock.h" /* in xlockmore distribution */
56 #endif /* STANDALONE */
60 #define DEF_EULERTAIL "10"
62 #define DEBUG_POINTED_REGION 0
65 static int variable_boundary = 1;
66 static float power = 1;
68 static XrmOptionDescRec opts[] =
70 {"-eulertail", ".euler2d.eulertail", XrmoptionSepArg, NULL},
71 {"-eulerpower", ".euler2d.eulerpower", XrmoptionSepArg, NULL},
73 static argtype vars[] =
75 {&tail_len, "eulertail",
76 "EulerTail", (char *) DEF_EULERTAIL, t_Int},
77 {&power, "eulerpower",
78 "EulerPower", "1", t_Float},
80 static OptionStruct desc[] =
82 {"-eulertail len", "Length of Euler2d tails"},
83 {"-eulerpower power", "power of interaction law for points for Euler2d"},
86 ENTRYPOINT ModeSpecOpt euler2d_opts =
87 {sizeof opts / sizeof opts[0], opts,
88 sizeof vars / sizeof vars[0], vars, desc};
91 ModStruct euler2d_description = {
92 "euler2d", "init_euler2d", "draw_euler2d", "release_euler2d",
93 "refresh_euler2d", "init_euler2d", (char *) NULL, &euler2d_opts,
94 1000, 1024, 3000, 1, 64, 1.0, "",
95 "Simulates 2D incompressible invisid fluid.", 0, NULL
100 #define balance_rand(v) ((LRAND()/MAXRAND*(v))-((v)/2)) /* random around 0 */
101 #define positive_rand(v) (LRAND()/MAXRAND*(v)) /* positive random */
103 #define number_of_vortex_points 20
105 #define n_bound_p 500
108 static double delta_t;
114 double xshift,yshift,scale;
120 /* x[2i+0] = x coord for nth point
121 x[2i+1] = y coord for nth point
122 w[i] = vorticity at nth point
131 /* (xs[2i+0],xs[2i+1]) is reflection of (x[2i+0],x[2i+1]) about unit circle
132 xs[2i+0] = x[2i+0]/nx
133 xs[2i+1] = x[2i+1]/nx
135 nx = x[2i+0]*x[2i+0] + x[2i+1]*x[2i+1]
137 x_is_zero[i] = (nx < 1e-10)
142 /* (p[2i+0],p[2i+1]) is image of (x[2i+0],x[2i+1]) under polynomial p.
143 mod_dp2 is |p'(z)|^2 when z = (x[2i+0],x[2i+1]).
148 /* Sometimes in our calculations we get overflow or numbers that are too big.
149 If that happens with the point x[2*i+0], x[2*i+1], we set dead[i].
162 double p_coef[2*(deg_p-1)];
167 static euler2dstruct *euler2ds = (euler2dstruct *) NULL;
170 If variable_boundary == 1, then we make a variable boundary.
171 The way this is done is to map the unit disk under a
173 p(z) = z + c_2 z^2 + ... + c_n z^n
174 where n = deg_p. sp->p_coef contains the complex numbers
178 #define add(a1,a2,b1,b2) (a1)+=(b1);(a2)+=(b2)
179 #define mult(a1,a2,b1,b2) temp=(a1)*(b1)-(a2)*(b2); \
180 (a2)=(a1)*(b2)+(a2)*(b1);(a1)=temp
183 calc_p(double *p1, double *p2, double z1, double z2, double p_coef[])
190 for(i=deg_p;i>=2;i--)
192 add(*p1,*p2,p_coef[(i-2)*2],p_coef[(i-2)*2+1]);
199 /* Calculate |p'(z)|^2 */
201 calc_mod_dp2(double z1, double z2, double p_coef[])
208 for(i=deg_p;i>=2;i--)
210 add(mp1,mp2,i*p_coef[(i-2)*2],i*p_coef[(i-2)*2+1]);
214 return mp1*mp1+mp2*mp2;
218 calc_all_p(euler2dstruct *sp)
221 double temp,p1,p2,z1,z2;
222 for(j=(sp->hide_vortex?sp->Nvortex:0);j<sp->N;j++) if(!sp->dead[j])
228 for(i=deg_p;i>=2;i--)
230 add(p1,p2,sp->p_coef[(i-2)*2],sp->p_coef[(i-2)*2+1]);
241 calc_all_mod_dp2(double *x, euler2dstruct *sp)
244 double temp,mp1,mp2,z1,z2;
245 for(j=0;j<sp->N;j++) if(!sp->dead[j])
251 for(i=deg_p;i>=2;i--)
253 add(mp1,mp2,i*sp->p_coef[(i-2)*2],i*sp->p_coef[(i-2)*2+1]);
257 sp->mod_dp2[j] = mp1*mp1+mp2*mp2;
262 derivs(double *x, euler2dstruct *sp)
265 double u1,u2,x1,x2,xij1,xij2,nxij;
268 if (variable_boundary)
269 calc_all_mod_dp2(sp->x,sp);
271 for (j=0;j<sp->Nvortex;j++) if (!sp->dead[j])
273 nx = x[2*j+0]*x[2*j+0] + x[2*j+1]*x[2*j+1];
275 sp->x_is_zero[j] = 1;
277 sp->x_is_zero[j] = 0;
278 sp->xs[2*j+0] = x[2*j+0]/nx;
279 sp->xs[2*j+1] = x[2*j+1]/nx;
283 (void) memset(sp->diffx,0,sizeof(double)*2*sp->N);
285 for (i=0;i<sp->N;i++) if (!sp->dead[i])
289 for (j=0;j<sp->Nvortex;j++) if (!sp->dead[j])
292 Calculate the Biot-Savart kernel, that is, effect of a
293 vortex point at a = (x[2*j+0],x[2*j+1]) at the point
294 x = (x1,x2), returning the vector field in (u1,u2).
296 In the plane, this is given by the formula
298 u = (x-a)/|x-a|^2 or zero if x=a.
300 However, in the unit disk we have to subtract from the
305 where as = a/|a|^2 is the reflection of a about the unit circle.
307 If however power != 1, then
309 u = (x-a)/|x-a|^(power+1) - |a|^(1-power) (x-as)/|x-as|^(power+1)
313 xij1 = x1 - x[2*j+0];
314 xij2 = x2 - x[2*j+1];
315 nxij = (power==1.0) ? xij1*xij1+xij2*xij2 : pow(xij1*xij1+xij2*xij2,(power+1)/2.0);
324 if (!sp->x_is_zero[j])
326 xij1 = x1 - sp->xs[2*j+0];
327 xij2 = x2 - sp->xs[2*j+1];
328 nxij = (power==1.0) ? xij1*xij1+xij2*xij2 : pow(xij1*xij1+xij2*xij2,(power+1)/2.0);
344 sp->diffx[2*i+0] += u1*sp->w[j];
345 sp->diffx[2*i+1] += u2*sp->w[j];
349 if (!sp->dead[i] && variable_boundary)
351 if (sp->mod_dp2[i] < 1e-5)
355 sp->diffx[2*i+0] /= sp->mod_dp2[i];
356 sp->diffx[2*i+1] /= sp->mod_dp2[i];
363 What perturb does is effectively
365 where k should be of order delta_t.
367 We have the option to do this more subtly by mapping points x
368 in the unit disk to points y in the plane, where y = f(|x|) x,
369 with f(t) = -log(1-t)/t.
371 This might reduce (but does not remove) problems where particles near
372 the edge of the boundary bounce around.
374 But it seems to be not that effective, so for now switch it off.
377 #define SUBTLE_PERTURB 0
380 perturb(double ret[], double x[], double k[], euler2dstruct *sp)
386 double d1,d2,t1,t2,mag,mag2,mlog1mmag,memmagdmag,xdotk;
387 for (i=0;i<sp->N;i++) if (!sp->dead[i])
393 mag2 = x1*x1 + x2*x2;
399 else if (mag2 > 1-1e-5)
404 mlog1mmag = -log(1-mag);
405 xdotk = x1*k1 + x2*k2;
406 t1 = (x1 + k1)*mlog1mmag/mag + x1*xdotk*(1.0/(1-mag)-mlog1mmag/mag)/mag/mag;
407 t2 = (x2 + k2)*mlog1mmag/mag + x2*xdotk*(1.0/(1-mag)-mlog1mmag/mag)/mag/mag;
408 mag = sqrt(t1*t1+t2*t2);
409 if (mag > 11.5 /* log(1e5) */)
413 memmagdmag = (mag>1e-5) ? ((1.0-exp(-mag))/mag) : (1-mag/2.0);
414 ret[2*i+0] = t1*memmagdmag;
415 ret[2*i+1] = t2*memmagdmag;
422 if (d1*d1+d2*d2 > 0.1)
429 for (i=0;i<sp->N;i++) if (!sp->dead[i])
435 if (k1*k1+k2*k2 > 0.1 || x1*x1+x2*x2 > 1-1e-5)
447 ode_solve(euler2dstruct *sp)
453 /* midpoint method */
455 (void) memcpy(sp->olddiffx,sp->diffx,sizeof(double)*2*sp->N);
456 for (i=0;i<sp->N;i++) if (!sp->dead[i]) {
457 sp->tempdiffx[2*i+0] = 0.5*delta_t*sp->diffx[2*i+0];
458 sp->tempdiffx[2*i+1] = 0.5*delta_t*sp->diffx[2*i+1];
460 perturb(sp->tempx,sp->x,sp->tempdiffx,sp);
461 derivs(sp->tempx,sp);
462 for (i=0;i<sp->N;i++) if (!sp->dead[i]) {
463 sp->tempdiffx[2*i+0] = delta_t*sp->diffx[2*i+0];
464 sp->tempdiffx[2*i+1] = delta_t*sp->diffx[2*i+1];
466 perturb(sp->x,sp->x,sp->tempdiffx,sp);
470 for (i=0;i<sp->N;i++) if (!sp->dead[i]) {
471 sp->tempdiffx[2*i+0] = delta_t*(1.5*sp->diffx[2*i+0] - 0.5*sp->olddiffx[2*i+0]);
472 sp->tempdiffx[2*i+1] = delta_t*(1.5*sp->diffx[2*i+1] - 0.5*sp->olddiffx[2*i+1]);
474 perturb(sp->x,sp->x,sp->tempdiffx,sp);
476 sp->olddiffx = sp->diffx;
481 #define deallocate(p,t) if (p!=NULL) {(void) free((void *) p); p=(t*)NULL; }
482 #define allocate(p,t,s) if ((p=(t*)malloc(sizeof(t)*s))==NULL)\
483 {free_euler2d(sp);return;}
486 free_euler2d(euler2dstruct *sp)
488 deallocate(sp->csegs, XSegment);
489 deallocate(sp->old_segs, XSegment);
490 deallocate(sp->nold_segs, int);
491 deallocate(sp->lastx, short);
492 deallocate(sp->x, double);
493 deallocate(sp->diffx, double);
494 deallocate(sp->w, double);
495 deallocate(sp->olddiffx, double);
496 deallocate(sp->tempdiffx, double);
497 deallocate(sp->tempx, double);
498 deallocate(sp->dead, short);
499 deallocate(sp->boundary, XSegment);
500 deallocate(sp->xs, double);
501 deallocate(sp->x_is_zero, short);
502 deallocate(sp->p, double);
503 deallocate(sp->mod_dp2, double);
507 init_euler2d (ModeInfo * mi)
509 #define nr_rotates 18 /* how many rotations to try to fill as much of screen as possible - must be even number */
512 double r,theta,x,y,w;
513 double mag,xscale,yscale,p1,p2;
514 double low[nr_rotates],high[nr_rotates],pp1,pp2,pn1,pn2,angle1,angle2,tempangle,dist,scale,bestscale,temp;
517 if (power<0.5) power = 0.5;
518 if (power>3.0) power = 3.0;
519 variable_boundary &= power == 1.0;
521 if (power>1.0) delta_t *= pow(0.1,power-1);
523 if (euler2ds == NULL) {
524 if ((euler2ds = (euler2dstruct *) calloc(MI_NUM_SCREENS(mi),
525 sizeof (euler2dstruct))) == NULL)
528 sp = &euler2ds[MI_SCREEN(mi)];
531 jwxyz_XSetAntiAliasing (MI_DISPLAY(mi), MI_GC(mi), False);
534 sp->boundary_color = NRAND(MI_NPIXELS(mi));
535 sp->hide_vortex = NRAND(4) != 0;
539 sp->width = MI_WIDTH(mi);
540 sp->height = MI_HEIGHT(mi);
542 sp->N = MI_COUNT(mi)+number_of_vortex_points;
543 sp->Nvortex = number_of_vortex_points;
545 if (tail_len < 1) { /* minimum tail */
548 if (tail_len > MI_CYCLES(mi)) { /* maximum tail */
549 tail_len = MI_CYCLES(mi);
552 /* Clear the background. */
557 /* Allocate memory. */
559 if (sp->csegs == NULL) {
560 allocate(sp->csegs, XSegment, sp->N);
561 allocate(sp->old_segs, XSegment, sp->N * tail_len);
562 allocate(sp->nold_segs, int, tail_len);
563 allocate(sp->lastx, short, sp->N * 2);
564 allocate(sp->x, double, sp->N * 2);
565 allocate(sp->diffx, double, sp->N * 2);
566 allocate(sp->w, double, sp->Nvortex);
567 allocate(sp->olddiffx, double, sp->N * 2);
568 allocate(sp->tempdiffx, double, sp->N * 2);
569 allocate(sp->tempx, double, sp->N * 2);
570 allocate(sp->dead, short, sp->N);
571 allocate(sp->boundary, XSegment, n_bound_p);
572 allocate(sp->xs, double, sp->Nvortex * 2);
573 allocate(sp->x_is_zero, short, sp->Nvortex);
574 allocate(sp->p, double, sp->N * 2);
575 allocate(sp->mod_dp2, double, sp->N);
577 for (i=0;i<tail_len;i++) {
578 sp->nold_segs[i] = 0;
581 (void) memset(sp->dead,0,sp->N*sizeof(short));
583 if (variable_boundary)
585 /* Initialize polynomial p */
587 The polynomial p(z) = z + c_2 z^2 + ... c_n z^n needs to be
588 a bijection of the unit disk onto its image. This is achieved
589 by insisting that sum_{k=2}^n k |c_k| <= 1. Actually we set
590 the inequality to be equality (to get more interesting shapes).
593 for(k=2;k<=deg_p;k++)
595 r = positive_rand(1.0/k);
596 theta = balance_rand(2*M_PI);
597 sp->p_coef[2*(k-2)+0]=r*cos(theta);
598 sp->p_coef[2*(k-2)+1]=r*sin(theta);
601 if (mag > 0.0001) for(k=2;k<=deg_p;k++)
603 sp->p_coef[2*(k-2)+0] /= mag;
604 sp->p_coef[2*(k-2)+1] /= mag;
607 #if DEBUG_POINTED_REGION
608 for(k=2;k<=deg_p;k++){
609 sp->p_coef[2*(k-2)+0]=0;
610 sp->p_coef[2*(k-2)+1]=0;
612 sp->p_coef[2*(6-2)+0] = 1.0/6.0;
616 /* Here we figure out the best rotation of the domain so that it fills as
617 much of the screen as possible. The number of angles we look at is determined
618 by nr_rotates (we look every 180/nr_rotates degrees).
619 While we figure out the best angle to rotate, we also figure out the correct scaling factors.
622 for(k=0;k<nr_rotates;k++) {
627 for(k=0;k<n_bound_p;k++) {
628 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);
629 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);
630 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);
631 angle1 = nr_rotates/M_PI*atan2(p2-pp2,p1-pp1)-nr_rotates/2;
632 angle2 = nr_rotates/M_PI*atan2(pn2-p2,pn1-p1)-nr_rotates/2;
633 while (angle1<0) angle1+=nr_rotates*2;
634 while (angle2<0) angle2+=nr_rotates*2;
635 if (angle1>nr_rotates*1.75 && angle2<nr_rotates*0.25) angle2+=nr_rotates*2;
636 if (angle1<nr_rotates*0.25 && angle2>nr_rotates*1.75) angle1+=nr_rotates*2;
642 for(i=(int)floor(angle1);i<(int)ceil(angle2);i++) {
643 dist = cos((double)i*M_PI/nr_rotates)*p1 + sin((double)i*M_PI/nr_rotates)*p2;
644 if (i%(nr_rotates*2)<nr_rotates) {
645 if (dist>high[i%nr_rotates]) high[i%nr_rotates] = dist;
646 if (dist<low[i%nr_rotates]) low[i%nr_rotates] = dist;
649 if (-dist>high[i%nr_rotates]) high[i%nr_rotates] = -dist;
650 if (-dist<low[i%nr_rotates]) low[i%nr_rotates] = -dist;
655 for (i=0;i<nr_rotates;i++) {
656 xscale = (sp->width-5.0)/(high[i]-low[i]);
657 yscale = (sp->height-5.0)/(high[(i+nr_rotates/2)%nr_rotates]-low[(i+nr_rotates/2)%nr_rotates]);
658 scale = (xscale>yscale) ? yscale : xscale;
659 if (scale>bestscale) {
664 /* Here we do the rotation. The way we do this is to replace the
665 polynomial p(z) by a^{-1} p(a z) where a = exp(i best_angle).
669 for(k=2;k<=deg_p;k++)
671 mult(p1,p2,cos((double)besti*M_PI/nr_rotates),sin((double)besti*M_PI/nr_rotates));
672 mult(sp->p_coef[2*(k-2)+0],sp->p_coef[2*(k-2)+1],p1,p2);
675 sp->scale = bestscale;
676 sp->xshift = -(low[besti]+high[besti])/2.0*sp->scale+sp->width/2;
677 if (besti<nr_rotates/2)
678 sp->yshift = -(low[besti+nr_rotates/2]+high[besti+nr_rotates/2])/2.0*sp->scale+sp->height/2;
680 sp->yshift = (low[besti-nr_rotates/2]+high[besti-nr_rotates/2])/2.0*sp->scale+sp->height/2;
683 /* Initialize boundary */
685 for(k=0;k<n_bound_p;k++)
688 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);
689 sp->boundary[k].x1 = (short)(p1*sp->scale+sp->xshift);
690 sp->boundary[k].y1 = (short)(p2*sp->scale+sp->yshift);
692 for(k=1;k<n_bound_p;k++)
694 sp->boundary[k].x2 = sp->boundary[k-1].x1;
695 sp->boundary[k].y2 = sp->boundary[k-1].y1;
697 sp->boundary[0].x2 = sp->boundary[n_bound_p-1].x1;
698 sp->boundary[0].y2 = sp->boundary[n_bound_p-1].y1;
702 if (sp->width>sp->height)
703 sp->radius = sp->height/2.0-5.0;
705 sp->radius = sp->width/2.0-5.0;
708 /* Initialize point positions */
710 for (i=sp->Nvortex;i<sp->N;i++) {
712 r = sqrt(positive_rand(1.0));
713 theta = balance_rand(2*M_PI);
714 sp->x[2*i+0]=r*cos(theta);
715 sp->x[2*i+1]=r*sin(theta);
716 /* This is to make sure the initial distribution of points is uniform */
717 } while (variable_boundary &&
718 calc_mod_dp2(sp->x[2*i+0],sp->x[2*i+1],sp->p_coef)
723 /* number of vortex points with negative vorticity */
728 /* if n is even make sure that np==n/2 is twice as likely
729 as the other possibilities. */
735 r = sqrt(positive_rand(0.77));
736 theta = balance_rand(2*M_PI);
739 r = 0.02+positive_rand(0.1);
740 w = (2*(k<np)-1)*2.0/sp->Nvortex;
741 for (i=sp->Nvortex*k/n;i<sp->Nvortex*(k+1)/n;i++) {
742 theta = balance_rand(2*M_PI);
743 sp->x[2*i+0]=x + r*cos(theta);
744 sp->x[2*i+1]=y + r*sin(theta);
751 draw_euler2d (ModeInfo * mi)
753 Display *display = MI_DISPLAY(mi);
754 Window window = MI_WINDOW(mi);
756 int b, col, n_non_vortex_segs;
759 MI_IS_DRAWN(mi) = True;
761 if (euler2ds == NULL)
763 sp = &euler2ds[MI_SCREEN(mi)];
764 if (sp->csegs == NULL)
768 if (variable_boundary)
772 for(b=sp->Nvortex;b<sp->N;b++) if(!sp->dead[b])
774 sp->csegs[sp->cnsegs].x1 = sp->lastx[2*b+0];
775 sp->csegs[sp->cnsegs].y1 = sp->lastx[2*b+1];
776 if (variable_boundary)
778 sp->csegs[sp->cnsegs].x2 = (short)(sp->p[2*b+0]*sp->scale+sp->xshift);
779 sp->csegs[sp->cnsegs].y2 = (short)(sp->p[2*b+1]*sp->scale+sp->yshift);
783 sp->csegs[sp->cnsegs].x2 = (short)(sp->x[2*b+0]*sp->radius+sp->width/2);
784 sp->csegs[sp->cnsegs].y2 = (short)(sp->x[2*b+1]*sp->radius+sp->height/2);
786 sp->lastx[2*b+0] = sp->csegs[sp->cnsegs].x2;
787 sp->lastx[2*b+1] = sp->csegs[sp->cnsegs].y2;
790 n_non_vortex_segs = sp->cnsegs;
792 if (!sp->hide_vortex) for(b=0;b<sp->Nvortex;b++) if(!sp->dead[b])
794 sp->csegs[sp->cnsegs].x1 = sp->lastx[2*b+0];
795 sp->csegs[sp->cnsegs].y1 = sp->lastx[2*b+1];
796 if (variable_boundary)
798 sp->csegs[sp->cnsegs].x2 = (short)(sp->p[2*b+0]*sp->scale+sp->xshift);
799 sp->csegs[sp->cnsegs].y2 = (short)(sp->p[2*b+1]*sp->scale+sp->yshift);
803 sp->csegs[sp->cnsegs].x2 = (short)(sp->x[2*b+0]*sp->radius+sp->width/2);
804 sp->csegs[sp->cnsegs].y2 = (short)(sp->x[2*b+1]*sp->radius+sp->height/2);
806 sp->lastx[2*b+0] = sp->csegs[sp->cnsegs].x2;
807 sp->lastx[2*b+1] = sp->csegs[sp->cnsegs].y2;
812 XSetForeground(display, gc, MI_BLACK_PIXEL(mi));
814 XDrawSegments(display, window, gc, sp->old_segs+sp->c_old_seg*sp->N, sp->nold_segs[sp->c_old_seg]);
816 if (MI_NPIXELS(mi) > 2){ /* render colour */
817 for (col = 0; col < MI_NPIXELS(mi); col++) {
818 int start = col*n_non_vortex_segs/MI_NPIXELS(mi);
819 int finish = (col+1)*n_non_vortex_segs/MI_NPIXELS(mi);
820 XSetForeground(display, gc, MI_PIXEL(mi, col));
821 XDrawSegments(display, window, gc,sp->csegs+start, finish-start);
823 if (!sp->hide_vortex) {
824 XSetForeground(display, gc, MI_WHITE_PIXEL(mi));
825 XDrawSegments(display, window, gc,sp->csegs+n_non_vortex_segs, sp->cnsegs-n_non_vortex_segs);
828 } else { /* render mono */
829 XSetForeground(display, gc, MI_WHITE_PIXEL(mi));
830 XDrawSegments(display, window, gc,
831 sp->csegs, sp->cnsegs);
834 if (MI_NPIXELS(mi) > 2) /* render colour */
835 XSetForeground(display, gc, MI_PIXEL(mi, sp->boundary_color));
837 XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_WHITE_PIXEL(mi));
838 if (variable_boundary)
839 XDrawSegments(display, window, gc,
840 sp->boundary, n_bound_p);
842 XDrawArc(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
843 sp->width/2 - (int) sp->radius - 1, sp->height/2 - (int) sp->radius -1,
844 (int) (2*sp->radius) + 2, (int) (2* sp->radius) + 2, 0, 64*360);
846 /* Copy to erase-list */
847 (void) memcpy(sp->old_segs+sp->c_old_seg*sp->N, sp->csegs, sp->cnsegs*sizeof(XSegment));
848 sp->nold_segs[sp->c_old_seg] = sp->cnsegs;
850 if (sp->c_old_seg >= tail_len)
854 if (++sp->count > MI_CYCLES(mi)) /* pick a new flow */
860 reshape_euler2d(ModeInfo * mi, int width, int height)
862 XClearWindow (MI_DISPLAY (mi), MI_WINDOW(mi));
867 release_euler2d (ModeInfo * mi)
869 if (euler2ds != NULL) {
872 for (screen = 0; screen < MI_NUM_SCREENS(mi); screen++)
873 free_euler2d(&euler2ds[screen]);
874 (void) free((void *) euler2ds);
875 euler2ds = (euler2dstruct *) NULL;
880 refresh_euler2d (ModeInfo * mi)
886 euler2d_handle_event (ModeInfo *mi, XEvent *event)
888 if (screenhack_event_helper (MI_DISPLAY(mi), MI_WINDOW(mi), event))
898 XSCREENSAVER_MODULE ("Euler2D", euler2d)
900 #endif /* MODE_euler2d */