1 /* -*- Mode: C; tab-width: 4 -*- */
2 /* euler2d --- 2 Dimensional Incompressible Inviscid Fluid Flow */
4 #if !defined( lint ) && !defined( SABER )
5 static const char sccsid[] = "@(#)euler2d.c 5.00 2000/11/01 xlockmore";
10 * Copyright (c) 2000 by Stephen Montgomery-Smith <stephen@math.missouri.edu>
12 * Permission to use, copy, modify, and distribute this software and its
13 * documentation for any purpose and without fee is hereby granted,
14 * provided that the above copyright notice appear in all copies and that
15 * both that copyright notice and this permission notice appear in
16 * supporting documentation.
18 * This file is provided AS IS with no warranties of any kind. The author
19 * shall have no liability with respect to the infringement of copyrights,
20 * trade secrets or any patents by this file or any part thereof. In no
21 * event will the author be liable for any lost revenue or profits or
22 * other special, indirect and consequential damages.
25 * 04-Nov-2000: Added an option eulerpower. This allows for example the
26 * quasi-geostrophic equation by setting eulerpower to 2.
27 * 01-Nov-2000: Allocation checks.
28 * 10-Sep-2000: Added optimizations, and removed subtle_perturb, by stephen.
29 * 03-Sep-2000: Changed method of solving ode to Adams-Bashforth of order 2.
30 * Previously used a rather compilcated method of order 4.
31 * This doubles the speed of the program. Also it seems
32 * to have improved numerical stability. Done by stephen.
33 * 27-Aug-2000: Added rotation of region to maximize screen fill by stephen.
34 * 05-Jun-2000: Adapted from flow.c Copyright (c) 1996 by Tim Auckland
35 * 18-Jul-1996: Adapted from swarm.c Copyright (c) 1991 by Patrick J. Naughton.
36 * 31-Aug-1990: Adapted from xswarm by Jeff Butterworth. (butterwo@ncsc.org)
40 * The mathematical aspects of this program are discussed in the file
46 #define PROGCLASS "Euler2d"
47 #define HACK_INIT init_euler2d
48 #define HACK_DRAW draw_euler2d
49 #define euler2d_opts xlockmore_opts
50 #define DEFAULTS "*delay: 10000 \n" \
55 #include "xlockmore.h" /* in xscreensaver distribution */
56 #else /* STANDALONE */
57 #include "xlock.h" /* in xlockmore distribution */
58 #endif /* STANDALONE */
62 #define DEF_EULERTAIL "10"
64 #define DEBUG_POINTED_REGION 0
67 static int variable_boundary = 1;
68 static float power = 1;
70 static XrmOptionDescRec opts[] =
72 {(char* ) "-eulertail", (char *) ".euler2d.eulertail",
73 XrmoptionSepArg, (caddr_t) NULL},
74 {(char* ) "-eulerpower", (char *) ".euler2d.eulerpower",
75 XrmoptionSepArg, (caddr_t) NULL},
77 static argtype vars[] =
79 {(caddr_t *) &tail_len, (char *) "eulertail",
80 (char *) "EulerTail", (char *) DEF_EULERTAIL, t_Int},
81 {(caddr_t *) &power, (char *) "eulerpower",
82 (char *) "EulerPower", (char *) "1", t_Float},
84 static OptionStruct desc[] =
86 {(char *) "-eulertail len", (char *) "Length of Euler2d tails"},
87 {(char *) "-eulerpower power", (char *) "power of interaction law for points for Euler2d"},
90 ModeSpecOpt euler2d_opts =
91 {sizeof opts / sizeof opts[0], opts,
92 sizeof vars / sizeof vars[0], vars, desc};
95 ModStruct euler2d_description = {
96 "euler2d", "init_euler2d", "draw_euler2d", "release_euler2d",
97 "refresh_euler2d", "init_euler2d", (char *) NULL, &euler2d_opts,
98 1000, 1024, 3000, 1, 64, 1.0, "",
99 "Simulates 2D incompressible invisid fluid.", 0, NULL
104 #define balance_rand(v) ((LRAND()/MAXRAND*(v))-((v)/2)) /* random around 0 */
105 #define positive_rand(v) (LRAND()/MAXRAND*(v)) /* positive random */
107 #define number_of_vortex_points 20
109 #define n_bound_p 500
112 static double delta_t;
118 double xshift,yshift,scale;
124 /* x[2i+0] = x coord for nth point
125 x[2i+1] = y coord for nth point
126 w[i] = vorticity at nth point
135 /* (xs[2i+0],xs[2i+1]) is reflection of (x[2i+0],x[2i+1]) about unit circle
136 xs[2i+0] = x[2i+0]/nx
137 xs[2i+1] = x[2i+1]/nx
139 nx = x[2i+0]*x[2i+0] + x[2i+1]*x[2i+1]
141 x_is_zero[i] = (nx < 1e-10)
146 /* (p[2i+0],p[2i+1]) is image of (x[2i+0],x[2i+1]) under polynomial p.
147 mod_dp2 is |p'(z)|^2 when z = (x[2i+0],x[2i+1]).
152 /* Sometimes in our calculations we get overflow or numbers that are too big.
153 If that happens with the point x[2*i+0], x[2*i+1], we set dead[i].
166 double p_coef[2*(deg_p-1)];
171 static euler2dstruct *euler2ds = (euler2dstruct *) NULL;
174 If variable_boundary == 1, then we make a variable boundary.
175 The way this is done is to map the unit disk under a
177 p(z) = z + c_2 z^2 + ... + c_n z^n
178 where n = deg_p. sp->p_coef contains the complex numbers
182 #define add(a1,a2,b1,b2) (a1)+=(b1);(a2)+=(b2)
183 #define mult(a1,a2,b1,b2) temp=(a1)*(b1)-(a2)*(b2); \
184 (a2)=(a1)*(b2)+(a2)*(b1);(a1)=temp
187 calc_p(double *p1, double *p2, double z1, double z2, double p_coef[])
194 for(i=deg_p;i>=2;i--)
196 add(*p1,*p2,p_coef[(i-2)*2],p_coef[(i-2)*2+1]);
203 /* Calculate |p'(z)|^2 */
205 calc_mod_dp2(double z1, double z2, double p_coef[])
212 for(i=deg_p;i>=2;i--)
214 add(mp1,mp2,i*p_coef[(i-2)*2],i*p_coef[(i-2)*2+1]);
218 return mp1*mp1+mp2*mp2;
222 calc_all_p(euler2dstruct *sp)
225 double temp,p1,p2,z1,z2;
226 for(j=(sp->hide_vortex?sp->Nvortex:0);j<sp->N;j++) if(!sp->dead[j])
232 for(i=deg_p;i>=2;i--)
234 add(p1,p2,sp->p_coef[(i-2)*2],sp->p_coef[(i-2)*2+1]);
245 calc_all_mod_dp2(double *x, euler2dstruct *sp)
248 double temp,mp1,mp2,z1,z2;
249 for(j=0;j<sp->N;j++) if(!sp->dead[j])
255 for(i=deg_p;i>=2;i--)
257 add(mp1,mp2,i*sp->p_coef[(i-2)*2],i*sp->p_coef[(i-2)*2+1]);
261 sp->mod_dp2[j] = mp1*mp1+mp2*mp2;
266 derivs(double *x, euler2dstruct *sp)
269 double u1,u2,x1,x2,xij1,xij2,nxij;
272 if (variable_boundary)
273 calc_all_mod_dp2(sp->x,sp);
275 for (j=0;j<sp->Nvortex;j++) if (!sp->dead[j])
277 nx = x[2*j+0]*x[2*j+0] + x[2*j+1]*x[2*j+1];
279 sp->x_is_zero[j] = 1;
281 sp->x_is_zero[j] = 0;
282 sp->xs[2*j+0] = x[2*j+0]/nx;
283 sp->xs[2*j+1] = x[2*j+1]/nx;
287 (void) memset(sp->diffx,0,sizeof(double)*2*sp->N);
289 for (i=0;i<sp->N;i++) if (!sp->dead[i])
293 for (j=0;j<sp->Nvortex;j++) if (!sp->dead[j])
296 Calculate the Biot-Savart kernel, that is, effect of a
297 vortex point at a = (x[2*j+0],x[2*j+1]) at the point
298 x = (x1,x2), returning the vector field in (u1,u2).
300 In the plane, this is given by the formula
302 u = (x-a)/|x-a|^2 or zero if x=a.
304 However, in the unit disk we have to subtract from the
309 where as = a/|a|^2 is the reflection of a about the unit circle.
311 If however power != 1, then
313 u = (x-a)/|x-a|^(power+1) - |a|^(1-power) (x-as)/|x-as|^(power+1)
317 xij1 = x1 - x[2*j+0];
318 xij2 = x2 - x[2*j+1];
319 nxij = (power==1.0) ? xij1*xij1+xij2*xij2 : pow(xij1*xij1+xij2*xij2,(power+1)/2.0);
328 if (!sp->x_is_zero[j])
330 xij1 = x1 - sp->xs[2*j+0];
331 xij2 = x2 - sp->xs[2*j+1];
332 nxij = (power==1.0) ? xij1*xij1+xij2*xij2 : pow(xij1*xij1+xij2*xij2,(power+1)/2.0);
348 sp->diffx[2*i+0] += u1*sp->w[j];
349 sp->diffx[2*i+1] += u2*sp->w[j];
353 if (!sp->dead[i] && variable_boundary)
355 if (sp->mod_dp2[i] < 1e-5)
359 sp->diffx[2*i+0] /= sp->mod_dp2[i];
360 sp->diffx[2*i+1] /= sp->mod_dp2[i];
367 What perturb does is effectively
369 where k should be of order delta_t.
371 We have the option to do this more subtly by mapping points x
372 in the unit disk to points y in the plane, where y = f(|x|) x,
373 with f(t) = -log(1-t)/t.
375 This might reduce (but does not remove) problems where particles near
376 the edge of the boundary bounce around.
378 But it seems to be not that effective, so for now switch it off.
381 #define SUBTLE_PERTURB 0
384 perturb(double ret[], double x[], double k[], euler2dstruct *sp)
390 double d1,d2,t1,t2,mag,mag2,mlog1mmag,memmagdmag,xdotk;
391 for (i=0;i<sp->N;i++) if (!sp->dead[i])
397 mag2 = x1*x1 + x2*x2;
403 else if (mag2 > 1-1e-5)
408 mlog1mmag = -log(1-mag);
409 xdotk = x1*k1 + x2*k2;
410 t1 = (x1 + k1)*mlog1mmag/mag + x1*xdotk*(1.0/(1-mag)-mlog1mmag/mag)/mag/mag;
411 t2 = (x2 + k2)*mlog1mmag/mag + x2*xdotk*(1.0/(1-mag)-mlog1mmag/mag)/mag/mag;
412 mag = sqrt(t1*t1+t2*t2);
413 if (mag > 11.5 /* log(1e5) */)
417 memmagdmag = (mag>1e-5) ? ((1.0-exp(-mag))/mag) : (1-mag/2.0);
418 ret[2*i+0] = t1*memmagdmag;
419 ret[2*i+1] = t2*memmagdmag;
426 if (d1*d1+d2*d2 > 0.1)
433 for (i=0;i<sp->N;i++) if (!sp->dead[i])
439 if (k1*k1+k2*k2 > 0.1 || x1*x1+x2*x2 > 1-1e-5)
451 ode_solve(euler2dstruct *sp)
457 /* midpoint method */
459 (void) memcpy(sp->olddiffx,sp->diffx,sizeof(double)*2*sp->N);
460 for (i=0;i<sp->N;i++) if (!sp->dead[i]) {
461 sp->tempdiffx[2*i+0] = 0.5*delta_t*sp->diffx[2*i+0];
462 sp->tempdiffx[2*i+1] = 0.5*delta_t*sp->diffx[2*i+1];
464 perturb(sp->tempx,sp->x,sp->tempdiffx,sp);
465 derivs(sp->tempx,sp);
466 for (i=0;i<sp->N;i++) if (!sp->dead[i]) {
467 sp->tempdiffx[2*i+0] = delta_t*sp->diffx[2*i+0];
468 sp->tempdiffx[2*i+1] = delta_t*sp->diffx[2*i+1];
470 perturb(sp->x,sp->x,sp->tempdiffx,sp);
474 for (i=0;i<sp->N;i++) if (!sp->dead[i]) {
475 sp->tempdiffx[2*i+0] = delta_t*(1.5*sp->diffx[2*i+0] - 0.5*sp->olddiffx[2*i+0]);
476 sp->tempdiffx[2*i+1] = delta_t*(1.5*sp->diffx[2*i+1] - 0.5*sp->olddiffx[2*i+1]);
478 perturb(sp->x,sp->x,sp->tempdiffx,sp);
480 sp->olddiffx = sp->diffx;
485 #define deallocate(p,t) if (p!=NULL) {(void) free((void *) p); p=(t*)NULL; }
486 #define allocate(p,t,s) if ((p=(t*)malloc(sizeof(t)*s))==NULL)\
487 {free_euler2d(sp);return;}
490 free_euler2d(euler2dstruct *sp)
492 deallocate(sp->csegs, XSegment);
493 deallocate(sp->old_segs, XSegment);
494 deallocate(sp->nold_segs, int);
495 deallocate(sp->lastx, short);
496 deallocate(sp->x, double);
497 deallocate(sp->diffx, double);
498 deallocate(sp->w, double);
499 deallocate(sp->olddiffx, double);
500 deallocate(sp->tempdiffx, double);
501 deallocate(sp->tempx, double);
502 deallocate(sp->dead, short);
503 deallocate(sp->boundary, XSegment);
504 deallocate(sp->xs, double);
505 deallocate(sp->x_is_zero, short);
506 deallocate(sp->p, double);
507 deallocate(sp->mod_dp2, double);
511 init_euler2d(ModeInfo * mi)
513 #define nr_rotates 18 /* how many rotations to try to fill as much of screen as possible - must be even number */
516 double r,theta,x,y,w;
517 double mag,xscale,yscale,p1,p2;
518 double low[nr_rotates],high[nr_rotates],pp1,pp2,pn1,pn2,angle1,angle2,tempangle,dist,scale,bestscale,temp;
521 if (power<0.5) power = 0.5;
522 if (power>3.0) power = 3.0;
523 variable_boundary &= power == 1.0;
525 if (power>1.0) delta_t *= pow(0.1,power-1);
527 if (euler2ds == NULL) {
528 if ((euler2ds = (euler2dstruct *) calloc(MI_NUM_SCREENS(mi),
529 sizeof (euler2dstruct))) == NULL)
532 sp = &euler2ds[MI_SCREEN(mi)];
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 release_euler2d(ModeInfo * mi)
862 if (euler2ds != NULL) {
865 for (screen = 0; screen < MI_NUM_SCREENS(mi); screen++)
866 free_euler2d(&euler2ds[screen]);
867 (void) free((void *) euler2ds);
868 euler2ds = (euler2dstruct *) NULL;
873 refresh_euler2d(ModeInfo * mi)
878 #endif /* MODE_euler2d */