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" \
51 # define reshape_euler2d 0
52 # define euler2d_handle_event 0
53 # define SMOOTH_COLORS
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", "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)];
532 jwxyz_XSetAntiAliasing (MI_DISPLAY(mi), MI_GC(mi), False);
535 sp->boundary_color = NRAND(MI_NPIXELS(mi));
536 sp->hide_vortex = NRAND(4) != 0;
540 sp->width = MI_WIDTH(mi);
541 sp->height = MI_HEIGHT(mi);
543 sp->N = MI_COUNT(mi)+number_of_vortex_points;
544 sp->Nvortex = number_of_vortex_points;
546 if (tail_len < 1) { /* minimum tail */
549 if (tail_len > MI_CYCLES(mi)) { /* maximum tail */
550 tail_len = MI_CYCLES(mi);
553 /* Clear the background. */
558 /* Allocate memory. */
560 if (sp->csegs == NULL) {
561 allocate(sp->csegs, XSegment, sp->N);
562 allocate(sp->old_segs, XSegment, sp->N * tail_len);
563 allocate(sp->nold_segs, int, tail_len);
564 allocate(sp->lastx, short, sp->N * 2);
565 allocate(sp->x, double, sp->N * 2);
566 allocate(sp->diffx, double, sp->N * 2);
567 allocate(sp->w, double, sp->Nvortex);
568 allocate(sp->olddiffx, double, sp->N * 2);
569 allocate(sp->tempdiffx, double, sp->N * 2);
570 allocate(sp->tempx, double, sp->N * 2);
571 allocate(sp->dead, short, sp->N);
572 allocate(sp->boundary, XSegment, n_bound_p);
573 allocate(sp->xs, double, sp->Nvortex * 2);
574 allocate(sp->x_is_zero, short, sp->Nvortex);
575 allocate(sp->p, double, sp->N * 2);
576 allocate(sp->mod_dp2, double, sp->N);
578 for (i=0;i<tail_len;i++) {
579 sp->nold_segs[i] = 0;
582 (void) memset(sp->dead,0,sp->N*sizeof(short));
584 if (variable_boundary)
586 /* Initialize polynomial p */
588 The polynomial p(z) = z + c_2 z^2 + ... c_n z^n needs to be
589 a bijection of the unit disk onto its image. This is achieved
590 by insisting that sum_{k=2}^n k |c_k| <= 1. Actually we set
591 the inequality to be equality (to get more interesting shapes).
594 for(k=2;k<=deg_p;k++)
596 r = positive_rand(1.0/k);
597 theta = balance_rand(2*M_PI);
598 sp->p_coef[2*(k-2)+0]=r*cos(theta);
599 sp->p_coef[2*(k-2)+1]=r*sin(theta);
602 if (mag > 0.0001) for(k=2;k<=deg_p;k++)
604 sp->p_coef[2*(k-2)+0] /= mag;
605 sp->p_coef[2*(k-2)+1] /= mag;
608 #if DEBUG_POINTED_REGION
609 for(k=2;k<=deg_p;k++){
610 sp->p_coef[2*(k-2)+0]=0;
611 sp->p_coef[2*(k-2)+1]=0;
613 sp->p_coef[2*(6-2)+0] = 1.0/6.0;
617 /* Here we figure out the best rotation of the domain so that it fills as
618 much of the screen as possible. The number of angles we look at is determined
619 by nr_rotates (we look every 180/nr_rotates degrees).
620 While we figure out the best angle to rotate, we also figure out the correct scaling factors.
623 for(k=0;k<nr_rotates;k++) {
628 for(k=0;k<n_bound_p;k++) {
629 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);
630 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);
631 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);
632 angle1 = nr_rotates/M_PI*atan2(p2-pp2,p1-pp1)-nr_rotates/2;
633 angle2 = nr_rotates/M_PI*atan2(pn2-p2,pn1-p1)-nr_rotates/2;
634 while (angle1<0) angle1+=nr_rotates*2;
635 while (angle2<0) angle2+=nr_rotates*2;
636 if (angle1>nr_rotates*1.75 && angle2<nr_rotates*0.25) angle2+=nr_rotates*2;
637 if (angle1<nr_rotates*0.25 && angle2>nr_rotates*1.75) angle1+=nr_rotates*2;
643 for(i=(int)floor(angle1);i<(int)ceil(angle2);i++) {
644 dist = cos((double)i*M_PI/nr_rotates)*p1 + sin((double)i*M_PI/nr_rotates)*p2;
645 if (i%(nr_rotates*2)<nr_rotates) {
646 if (dist>high[i%nr_rotates]) high[i%nr_rotates] = dist;
647 if (dist<low[i%nr_rotates]) low[i%nr_rotates] = dist;
650 if (-dist>high[i%nr_rotates]) high[i%nr_rotates] = -dist;
651 if (-dist<low[i%nr_rotates]) low[i%nr_rotates] = -dist;
656 for (i=0;i<nr_rotates;i++) {
657 xscale = (sp->width-5.0)/(high[i]-low[i]);
658 yscale = (sp->height-5.0)/(high[(i+nr_rotates/2)%nr_rotates]-low[(i+nr_rotates/2)%nr_rotates]);
659 scale = (xscale>yscale) ? yscale : xscale;
660 if (scale>bestscale) {
665 /* Here we do the rotation. The way we do this is to replace the
666 polynomial p(z) by a^{-1} p(a z) where a = exp(i best_angle).
670 for(k=2;k<=deg_p;k++)
672 mult(p1,p2,cos((double)besti*M_PI/nr_rotates),sin((double)besti*M_PI/nr_rotates));
673 mult(sp->p_coef[2*(k-2)+0],sp->p_coef[2*(k-2)+1],p1,p2);
676 sp->scale = bestscale;
677 sp->xshift = -(low[besti]+high[besti])/2.0*sp->scale+sp->width/2;
678 if (besti<nr_rotates/2)
679 sp->yshift = -(low[besti+nr_rotates/2]+high[besti+nr_rotates/2])/2.0*sp->scale+sp->height/2;
681 sp->yshift = (low[besti-nr_rotates/2]+high[besti-nr_rotates/2])/2.0*sp->scale+sp->height/2;
684 /* Initialize boundary */
686 for(k=0;k<n_bound_p;k++)
689 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);
690 sp->boundary[k].x1 = (short)(p1*sp->scale+sp->xshift);
691 sp->boundary[k].y1 = (short)(p2*sp->scale+sp->yshift);
693 for(k=1;k<n_bound_p;k++)
695 sp->boundary[k].x2 = sp->boundary[k-1].x1;
696 sp->boundary[k].y2 = sp->boundary[k-1].y1;
698 sp->boundary[0].x2 = sp->boundary[n_bound_p-1].x1;
699 sp->boundary[0].y2 = sp->boundary[n_bound_p-1].y1;
703 if (sp->width>sp->height)
704 sp->radius = sp->height/2.0-5.0;
706 sp->radius = sp->width/2.0-5.0;
709 /* Initialize point positions */
711 for (i=sp->Nvortex;i<sp->N;i++) {
713 r = sqrt(positive_rand(1.0));
714 theta = balance_rand(2*M_PI);
715 sp->x[2*i+0]=r*cos(theta);
716 sp->x[2*i+1]=r*sin(theta);
717 /* This is to make sure the initial distribution of points is uniform */
718 } while (variable_boundary &&
719 calc_mod_dp2(sp->x[2*i+0],sp->x[2*i+1],sp->p_coef)
724 /* number of vortex points with negative vorticity */
729 /* if n is even make sure that np==n/2 is twice as likely
730 as the other possibilities. */
736 r = sqrt(positive_rand(0.77));
737 theta = balance_rand(2*M_PI);
740 r = 0.02+positive_rand(0.1);
741 w = (2*(k<np)-1)*2.0/sp->Nvortex;
742 for (i=sp->Nvortex*k/n;i<sp->Nvortex*(k+1)/n;i++) {
743 theta = balance_rand(2*M_PI);
744 sp->x[2*i+0]=x + r*cos(theta);
745 sp->x[2*i+1]=y + r*sin(theta);
752 draw_euler2d (ModeInfo * mi)
754 Display *display = MI_DISPLAY(mi);
755 Window window = MI_WINDOW(mi);
757 int b, col, n_non_vortex_segs;
760 MI_IS_DRAWN(mi) = True;
762 if (euler2ds == NULL)
764 sp = &euler2ds[MI_SCREEN(mi)];
765 if (sp->csegs == NULL)
769 if (variable_boundary)
773 for(b=sp->Nvortex;b<sp->N;b++) if(!sp->dead[b])
775 sp->csegs[sp->cnsegs].x1 = sp->lastx[2*b+0];
776 sp->csegs[sp->cnsegs].y1 = sp->lastx[2*b+1];
777 if (variable_boundary)
779 sp->csegs[sp->cnsegs].x2 = (short)(sp->p[2*b+0]*sp->scale+sp->xshift);
780 sp->csegs[sp->cnsegs].y2 = (short)(sp->p[2*b+1]*sp->scale+sp->yshift);
784 sp->csegs[sp->cnsegs].x2 = (short)(sp->x[2*b+0]*sp->radius+sp->width/2);
785 sp->csegs[sp->cnsegs].y2 = (short)(sp->x[2*b+1]*sp->radius+sp->height/2);
787 sp->lastx[2*b+0] = sp->csegs[sp->cnsegs].x2;
788 sp->lastx[2*b+1] = sp->csegs[sp->cnsegs].y2;
791 n_non_vortex_segs = sp->cnsegs;
793 if (!sp->hide_vortex) for(b=0;b<sp->Nvortex;b++) if(!sp->dead[b])
795 sp->csegs[sp->cnsegs].x1 = sp->lastx[2*b+0];
796 sp->csegs[sp->cnsegs].y1 = sp->lastx[2*b+1];
797 if (variable_boundary)
799 sp->csegs[sp->cnsegs].x2 = (short)(sp->p[2*b+0]*sp->scale+sp->xshift);
800 sp->csegs[sp->cnsegs].y2 = (short)(sp->p[2*b+1]*sp->scale+sp->yshift);
804 sp->csegs[sp->cnsegs].x2 = (short)(sp->x[2*b+0]*sp->radius+sp->width/2);
805 sp->csegs[sp->cnsegs].y2 = (short)(sp->x[2*b+1]*sp->radius+sp->height/2);
807 sp->lastx[2*b+0] = sp->csegs[sp->cnsegs].x2;
808 sp->lastx[2*b+1] = sp->csegs[sp->cnsegs].y2;
813 XSetForeground(display, gc, MI_BLACK_PIXEL(mi));
815 XDrawSegments(display, window, gc, sp->old_segs+sp->c_old_seg*sp->N, sp->nold_segs[sp->c_old_seg]);
817 if (MI_NPIXELS(mi) > 2){ /* render colour */
818 for (col = 0; col < MI_NPIXELS(mi); col++) {
819 int start = col*n_non_vortex_segs/MI_NPIXELS(mi);
820 int finish = (col+1)*n_non_vortex_segs/MI_NPIXELS(mi);
821 XSetForeground(display, gc, MI_PIXEL(mi, col));
822 XDrawSegments(display, window, gc,sp->csegs+start, finish-start);
824 if (!sp->hide_vortex) {
825 XSetForeground(display, gc, MI_WHITE_PIXEL(mi));
826 XDrawSegments(display, window, gc,sp->csegs+n_non_vortex_segs, sp->cnsegs-n_non_vortex_segs);
829 } else { /* render mono */
830 XSetForeground(display, gc, MI_WHITE_PIXEL(mi));
831 XDrawSegments(display, window, gc,
832 sp->csegs, sp->cnsegs);
835 if (MI_NPIXELS(mi) > 2) /* render colour */
836 XSetForeground(display, gc, MI_PIXEL(mi, sp->boundary_color));
838 XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_WHITE_PIXEL(mi));
839 if (variable_boundary)
840 XDrawSegments(display, window, gc,
841 sp->boundary, n_bound_p);
843 XDrawArc(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
844 sp->width/2 - (int) sp->radius - 1, sp->height/2 - (int) sp->radius -1,
845 (int) (2*sp->radius) + 2, (int) (2* sp->radius) + 2, 0, 64*360);
847 /* Copy to erase-list */
848 (void) memcpy(sp->old_segs+sp->c_old_seg*sp->N, sp->csegs, sp->cnsegs*sizeof(XSegment));
849 sp->nold_segs[sp->c_old_seg] = sp->cnsegs;
851 if (sp->c_old_seg >= tail_len)
855 if (++sp->count > MI_CYCLES(mi)) /* pick a new flow */
861 release_euler2d (ModeInfo * mi)
863 if (euler2ds != NULL) {
866 for (screen = 0; screen < MI_NUM_SCREENS(mi); screen++)
867 free_euler2d(&euler2ds[screen]);
868 (void) free((void *) euler2ds);
869 euler2ds = (euler2dstruct *) NULL;
874 refresh_euler2d (ModeInfo * mi)
879 XSCREENSAVER_MODULE ("Euler2D", euler2d)
881 #endif /* MODE_euler2d */