From http://www.jwz.org/xscreensaver/xscreensaver-5.30.tar.gz
[xscreensaver] / hacks / euler2d.c
1 /* -*- Mode: C; tab-width: 4 -*- */
2 /* euler2d --- 2 Dimensional Incompressible Inviscid Fluid Flow */
3
4 #if 0
5 static const char sccsid[] = "@(#)euler2d.c     5.00 2000/11/01 xlockmore";
6 #endif
7
8 /*
9  * Copyright (c) 2000 by Stephen Montgomery-Smith <stephen@math.missouri.edu>
10  *
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.
16  *
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.
22  *
23  * Revision History:
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)
36  */
37
38 /*
39  * The mathematical aspects of this program are discussed in the file
40  * euler2d.tex.
41  */
42
43 #ifdef STANDALONE
44 # define MODE_euler2d
45 # define DEFAULTS       "*delay:   10000 \n" \
46                                         "*count:   1024  \n" \
47                                         "*cycles:  3000  \n" \
48                                         "*ncolors: 64    \n" \
49                                         "*fpsSolid: true    \n" \
50                                         "*ignoreRotation: True \n" \
51
52 # define SMOOTH_COLORS
53 # include "xlockmore.h"         /* in xscreensaver distribution */
54 #else /* STANDALONE */
55 # include "xlock.h"             /* in xlockmore distribution */
56 #endif /* STANDALONE */
57
58 #ifdef MODE_euler2d
59
60 #define DEF_EULERTAIL "10"
61
62 #define DEBUG_POINTED_REGION    0
63
64 static int tail_len;
65 static int variable_boundary = 1;
66 static float power = 1;
67
68 static XrmOptionDescRec opts[] =
69 {
70   {"-eulertail", ".euler2d.eulertail",   XrmoptionSepArg, NULL},
71   {"-eulerpower", ".euler2d.eulerpower", XrmoptionSepArg, NULL},
72 };
73 static argtype vars[] =
74 {
75   {&tail_len, "eulertail",
76    "EulerTail", (char *) DEF_EULERTAIL, t_Int},
77   {&power, "eulerpower",
78    "EulerPower", "1", t_Float},
79 };
80 static OptionStruct desc[] =
81 {
82   {"-eulertail len", "Length of Euler2d tails"},
83   {"-eulerpower power", "power of interaction law for points for Euler2d"},
84 };
85
86 ENTRYPOINT ModeSpecOpt euler2d_opts =
87 {sizeof opts / sizeof opts[0], opts,
88  sizeof vars / sizeof vars[0], vars, desc};
89
90 #ifdef USE_MODULES
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
96 };
97
98 #endif
99
100 #define balance_rand(v) ((LRAND()/MAXRAND*(v))-((v)/2)) /* random around 0 */
101 #define positive_rand(v)        (LRAND()/MAXRAND*(v))   /* positive random */
102
103 #define number_of_vortex_points 20
104
105 #define n_bound_p 500
106 #define deg_p 6
107
108 static double delta_t;
109
110 typedef struct {
111         int        width;
112         int        height;
113         int        count;
114         double     xshift,yshift,scale;
115         double     radius;
116
117         int        N;
118         int        Nvortex;
119
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
123 */
124         double     *x;
125         double     *w;
126
127         double     *diffx;
128         double     *olddiffx;
129         double     *tempx;
130         double     *tempdiffx;
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
134     where
135     nx = x[2i+0]*x[2i+0] + x[2i+1]*x[2i+1]
136
137     x_is_zero[i] = (nx < 1e-10)
138 */
139         double     *xs;
140         short      *x_is_zero;
141
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]).
144 */
145         double     *p;
146         double     *mod_dp2;
147
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].
150 */
151         short      *dead;
152
153         XSegment   *csegs;
154         int         cnsegs;
155         XSegment   *old_segs;
156         int        *nold_segs;
157         int         c_old_seg;
158         int         boundary_color;
159         int         hide_vortex;
160         short      *lastx;
161
162         double      p_coef[2*(deg_p-1)];
163         XSegment   *boundary;
164
165 } euler2dstruct;
166
167 static euler2dstruct *euler2ds = (euler2dstruct *) NULL;
168
169 /*
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 
172   polynomial p, where 
173   p(z) = z + c_2 z^2 + ... + c_n z^n
174   where n = deg_p.  sp->p_coef contains the complex numbers
175   c_2, c_3, ... c_n.
176 */
177
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
181
182 static void
183 calc_p(double *p1, double *p2, double z1, double z2, double p_coef[])
184 {
185   int i;
186   double temp;
187
188   *p1=0;
189   *p2=0;
190   for(i=deg_p;i>=2;i--)
191   {
192     add(*p1,*p2,p_coef[(i-2)*2],p_coef[(i-2)*2+1]);
193     mult(*p1,*p2,z1,z2);
194   }
195   add(*p1,*p2,1,0);
196   mult(*p1,*p2,z1,z2);
197 }
198
199 /* Calculate |p'(z)|^2 */
200 static double
201 calc_mod_dp2(double z1, double z2, double p_coef[])
202 {
203   int i;
204   double temp,mp1,mp2;
205
206   mp1=0;
207   mp2=0;
208   for(i=deg_p;i>=2;i--)
209   {
210     add(mp1,mp2,i*p_coef[(i-2)*2],i*p_coef[(i-2)*2+1]);
211     mult(mp1,mp2,z1,z2);
212   }
213   add(mp1,mp2,1,0);
214   return mp1*mp1+mp2*mp2;
215 }
216
217 static void
218 calc_all_p(euler2dstruct *sp)
219 {
220   int i,j;
221   double temp,p1,p2,z1,z2;
222   for(j=(sp->hide_vortex?sp->Nvortex:0);j<sp->N;j++) if(!sp->dead[j])
223   {
224     p1=0;
225     p2=0;
226     z1=sp->x[2*j+0];
227     z2=sp->x[2*j+1];
228     for(i=deg_p;i>=2;i--)
229     {
230       add(p1,p2,sp->p_coef[(i-2)*2],sp->p_coef[(i-2)*2+1]);
231       mult(p1,p2,z1,z2);
232     }
233     add(p1,p2,1,0);
234     mult(p1,p2,z1,z2);
235     sp->p[2*j+0] = p1;
236     sp->p[2*j+1] = p2;
237   }
238 }
239
240 static void
241 calc_all_mod_dp2(double *x, euler2dstruct *sp)
242 {
243   int i,j;
244   double temp,mp1,mp2,z1,z2;
245   for(j=0;j<sp->N;j++) if(!sp->dead[j])
246   {
247     mp1=0;
248     mp2=0;
249     z1=x[2*j+0];
250     z2=x[2*j+1];
251     for(i=deg_p;i>=2;i--)
252     {
253       add(mp1,mp2,i*sp->p_coef[(i-2)*2],i*sp->p_coef[(i-2)*2+1]);
254       mult(mp1,mp2,z1,z2);
255     }
256     add(mp1,mp2,1,0);
257     sp->mod_dp2[j] = mp1*mp1+mp2*mp2;
258   }
259 }
260
261 static void
262 derivs(double *x, euler2dstruct *sp)
263 {
264   int i,j;
265   double u1,u2,x1,x2,xij1,xij2,nxij;
266   double nx;
267
268   if (variable_boundary)
269     calc_all_mod_dp2(sp->x,sp);
270
271   for (j=0;j<sp->Nvortex;j++) if (!sp->dead[j])
272   {
273     nx = x[2*j+0]*x[2*j+0] + x[2*j+1]*x[2*j+1];
274     if (nx < 1e-10)
275       sp->x_is_zero[j] = 1;
276     else {
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;
280     }
281   }
282
283   (void) memset(sp->diffx,0,sizeof(double)*2*sp->N);
284
285   for (i=0;i<sp->N;i++) if (!sp->dead[i])
286   {
287     x1 = x[2*i+0];
288     x2 = x[2*i+1];
289     for (j=0;j<sp->Nvortex;j++) if (!sp->dead[j])
290     {
291 /*
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).
295
296   In the plane, this is given by the formula
297
298   u = (x-a)/|x-a|^2  or zero if x=a.
299
300   However, in the unit disk we have to subtract from the
301   above:
302
303   (x-as)/|x-as|^2
304
305   where as = a/|a|^2 is the reflection of a about the unit circle.
306
307   If however power != 1, then
308
309   u = (x-a)/|x-a|^(power+1)  -  |a|^(1-power) (x-as)/|x-as|^(power+1)
310
311 */
312
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);
316
317       if(nxij >= 1e-4)  {
318         u1 =  xij2/nxij;
319         u2 = -xij1/nxij;
320       }
321       else
322         u1 = u2 = 0.0;
323
324       if (!sp->x_is_zero[j])
325       {
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);
329
330         if (nxij < 1e-5)
331         {
332           sp->dead[i] = 1;
333           u1 = u2 = 0.0;
334         }
335         else
336         {
337           u1 -= xij2/nxij;
338           u2 += xij1/nxij;
339         }
340       }
341
342       if (!sp->dead[i])
343       {
344         sp->diffx[2*i+0] += u1*sp->w[j];
345         sp->diffx[2*i+1] += u2*sp->w[j];
346       }
347     }
348
349     if (!sp->dead[i] && variable_boundary)
350     {
351       if (sp->mod_dp2[i] < 1e-5)
352         sp->dead[i] = 1;
353       else
354       {
355         sp->diffx[2*i+0] /= sp->mod_dp2[i];
356         sp->diffx[2*i+1] /= sp->mod_dp2[i];
357       }
358     }
359   }
360 }
361
362 /*
363   What perturb does is effectively
364     ret = x + k,
365   where k should be of order delta_t.  
366
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.
370
371   This might reduce (but does not remove) problems where particles near 
372   the edge of the boundary bounce around.
373
374   But it seems to be not that effective, so for now switch it off.
375 */
376
377 #define SUBTLE_PERTURB 0
378
379 static void
380 perturb(double ret[], double x[], double k[], euler2dstruct *sp)
381 {
382   int i;
383   double x1,x2,k1,k2;
384
385 #if SUBTLE_PERTURB
386   double d1,d2,t1,t2,mag,mag2,mlog1mmag,memmagdmag,xdotk;
387   for (i=0;i<sp->N;i++) if (!sp->dead[i])
388   {
389     x1 = x[2*i+0];
390     x2 = x[2*i+1];
391     k1 = k[2*i+0];
392     k2 = k[2*i+1];
393     mag2 = x1*x1 + x2*x2;
394     if (mag2 < 1e-10)
395     {
396       ret[2*i+0] = x1+k1;
397       ret[2*i+1] = x2+k2;
398     }
399     else if (mag2 > 1-1e-5)
400       sp->dead[i] = 1;
401     else
402     {
403       mag = sqrt(mag2);
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) */)
410         sp->dead[i] = 1;
411       else
412       {
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;
416       }
417     }
418     if (!sp->dead[i])
419     {
420       d1 = ret[2*i+0]-x1;
421       d2 = ret[2*i+1]-x2;
422       if (d1*d1+d2*d2 > 0.1)
423         sp->dead[i] = 1;
424     }
425   }
426
427 #else
428
429   for (i=0;i<sp->N;i++) if (!sp->dead[i])
430   {
431     x1 = x[2*i+0];
432     x2 = x[2*i+1];
433     k1 = k[2*i+0];
434     k2 = k[2*i+1];
435     if (k1*k1+k2*k2 > 0.1 || x1*x1+x2*x2 > 1-1e-5)
436       sp->dead[i] = 1;
437     else
438     {
439       ret[2*i+0] = x1+k1;
440       ret[2*i+1] = x2+k2;
441     }
442   }
443 #endif
444 }
445
446 static void
447 ode_solve(euler2dstruct *sp)
448 {
449   int i;
450   double *temp;
451
452   if (sp->count < 1) {
453     /* midpoint method */
454     derivs(sp->x,sp);
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];
459     }
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];
465     }
466     perturb(sp->x,sp->x,sp->tempdiffx,sp);
467   } else {
468     /* Adams Basforth */
469     derivs(sp->x,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]);
473     }
474     perturb(sp->x,sp->x,sp->tempdiffx,sp);
475     temp = sp->olddiffx;
476     sp->olddiffx = sp->diffx;
477     sp->diffx = temp;
478   }
479 }
480
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;}
484
485 static void
486 free_euler2d(euler2dstruct *sp)
487 {
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);
504 }
505
506 ENTRYPOINT void
507 init_euler2d (ModeInfo * mi)
508 {
509 #define nr_rotates 18 /* how many rotations to try to fill as much of screen as possible - must be even number */
510         euler2dstruct *sp;
511         int         i,k,n,np;
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;
515         int         besti = 0;
516
517         if (power<0.5) power = 0.5;
518         if (power>3.0) power = 3.0;
519         variable_boundary &= power == 1.0;
520         delta_t = 0.001;
521         if (power>1.0) delta_t *= pow(0.1,power-1);
522
523         if (euler2ds == NULL) {
524                 if ((euler2ds = (euler2dstruct *) calloc(MI_NUM_SCREENS(mi),
525                                                sizeof (euler2dstruct))) == NULL)
526                         return;
527         }
528         sp = &euler2ds[MI_SCREEN(mi)];
529
530 #ifdef HAVE_COCOA
531   jwxyz_XSetAntiAliasing (MI_DISPLAY(mi), MI_GC(mi),  False);
532 #endif
533
534         sp->boundary_color = NRAND(MI_NPIXELS(mi));
535         sp->hide_vortex = NRAND(4) != 0;
536
537         sp->count = 0;
538
539         sp->width = MI_WIDTH(mi);
540         sp->height = MI_HEIGHT(mi);
541
542         sp->N = MI_COUNT(mi)+number_of_vortex_points;
543         sp->Nvortex = number_of_vortex_points;
544
545         if (tail_len < 1) { /* minimum tail */
546           tail_len = 1;
547         }
548         if (tail_len > MI_CYCLES(mi)) { /* maximum tail */
549           tail_len = MI_CYCLES(mi);
550         }
551
552         /* Clear the background. */
553         MI_CLEARWINDOW(mi);
554          
555         free_euler2d(sp);
556
557         /* Allocate memory. */
558
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);
576         }
577         for (i=0;i<tail_len;i++) {
578                 sp->nold_segs[i] = 0;
579         }
580         sp->c_old_seg = 0;
581         (void) memset(sp->dead,0,sp->N*sizeof(short));
582
583         if (variable_boundary)
584         {
585         /* Initialize polynomial p */
586 /*
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).
591 */
592                 mag = 0;
593                 for(k=2;k<=deg_p;k++)
594                 {
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);
599                         mag += k*r;
600                 }
601                 if (mag > 0.0001) for(k=2;k<=deg_p;k++)
602                 {
603                         sp->p_coef[2*(k-2)+0] /= mag;
604                         sp->p_coef[2*(k-2)+1] /= mag;
605                 }
606
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;
611                 }
612                 sp->p_coef[2*(6-2)+0] = 1.0/6.0;
613 #endif
614
615
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.
620 */
621
622                 for(k=0;k<nr_rotates;k++) {
623                         low[k] = 1e5;
624                         high[k] = -1e5;
625                 }
626
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;
637                         if (angle2<angle1) {
638                                 tempangle=angle1;
639                                 angle1=angle2;
640                                 angle2=tempangle;
641                         }
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;
647                                 }
648                                 else {
649                                         if (-dist>high[i%nr_rotates]) high[i%nr_rotates] = -dist;
650                                         if (-dist<low[i%nr_rotates]) low[i%nr_rotates] = -dist;
651                                 }
652                         }
653                 }
654                 bestscale = 0;
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) {
660                                 bestscale = scale;
661                                 besti = i;
662                         }
663                 }
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).
666 */
667                 p1 = 1;
668                 p2 = 0;
669                 for(k=2;k<=deg_p;k++)
670                 {
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);
673                 }
674
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;
679                 else
680                         sp->yshift = (low[besti-nr_rotates/2]+high[besti-nr_rotates/2])/2.0*sp->scale+sp->height/2;
681
682
683 /* Initialize boundary */
684
685                 for(k=0;k<n_bound_p;k++)
686                 {
687
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);
691                 }
692                 for(k=1;k<n_bound_p;k++)
693                 {
694                         sp->boundary[k].x2 = sp->boundary[k-1].x1;
695                         sp->boundary[k].y2 = sp->boundary[k-1].y1;
696                 }
697                 sp->boundary[0].x2 = sp->boundary[n_bound_p-1].x1;
698                 sp->boundary[0].y2 = sp->boundary[n_bound_p-1].y1;
699         }
700         else
701         {
702                 if (sp->width>sp->height)
703                         sp->radius = sp->height/2.0-5.0;
704                 else
705                         sp->radius = sp->width/2.0-5.0;
706         }
707
708         /* Initialize point positions */
709
710         for (i=sp->Nvortex;i<sp->N;i++) {
711                 do {
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)
719                           < positive_rand(4));
720         }
721
722         n = NRAND(4)+2;
723         /* number of vortex points with negative vorticity */
724         if (n%2) {
725                 np = NRAND(n+1); 
726         }
727         else {
728                 /* if n is even make sure that np==n/2 is twice as likely
729                    as the other possibilities. */
730                 np = NRAND(n+2);
731                 if (np==n+1) np=n/2;
732         }
733         for(k=0;k<n;k++)
734         {
735                 r = sqrt(positive_rand(0.77));
736                 theta = balance_rand(2*M_PI);
737                 x=r*cos(theta);
738                 y=r*sin(theta);
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);
745                         sp->w[i]=w;
746                 }
747         }
748 }
749
750 ENTRYPOINT void
751 draw_euler2d (ModeInfo * mi)
752 {
753         Display    *display = MI_DISPLAY(mi);
754         Window      window = MI_WINDOW(mi);
755         GC          gc = MI_GC(mi);
756         int         b, col, n_non_vortex_segs;
757         euler2dstruct *sp;
758
759         MI_IS_DRAWN(mi) = True;
760
761         if (euler2ds == NULL)
762                 return;
763         sp = &euler2ds[MI_SCREEN(mi)];
764         if (sp->csegs == NULL)
765                 return;
766
767         ode_solve(sp);
768         if (variable_boundary)
769                 calc_all_p(sp);
770
771         sp->cnsegs = 0;
772         for(b=sp->Nvortex;b<sp->N;b++) if(!sp->dead[b])
773         {
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)
777                 {
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);
780                 }
781                 else
782                 {
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);
785                 }
786                 sp->lastx[2*b+0] = sp->csegs[sp->cnsegs].x2;
787                 sp->lastx[2*b+1] = sp->csegs[sp->cnsegs].y2;
788                 sp->cnsegs++;
789         }
790         n_non_vortex_segs = sp->cnsegs;
791
792         if (!sp->hide_vortex) for(b=0;b<sp->Nvortex;b++) if(!sp->dead[b])
793         {
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)
797                 {
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);
800                 }
801                 else
802                 {
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);
805                 }
806                 sp->lastx[2*b+0] = sp->csegs[sp->cnsegs].x2;
807                 sp->lastx[2*b+1] = sp->csegs[sp->cnsegs].y2;
808                 sp->cnsegs++;
809         }
810
811         if (sp->count) {
812                 XSetForeground(display, gc, MI_BLACK_PIXEL(mi));
813
814                 XDrawSegments(display, window, gc, sp->old_segs+sp->c_old_seg*sp->N, sp->nold_segs[sp->c_old_seg]);
815
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);
822                         }
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);
826                         }
827
828                 } else {                /* render mono */
829                   XSetForeground(display, gc, MI_WHITE_PIXEL(mi));
830                   XDrawSegments(display, window, gc,
831                                                 sp->csegs, sp->cnsegs);
832                 }
833
834                 if (MI_NPIXELS(mi) > 2) /* render colour */
835                         XSetForeground(display, gc, MI_PIXEL(mi, sp->boundary_color));
836                 else
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);
841                 else
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);
845
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;
849                 sp->c_old_seg++;
850                 if (sp->c_old_seg >= tail_len)
851                   sp->c_old_seg = 0;
852         }
853
854         if (++sp->count > MI_CYCLES(mi)) /* pick a new flow */
855                 init_euler2d(mi);
856
857 }
858
859 ENTRYPOINT void
860 reshape_euler2d(ModeInfo * mi, int width, int height)
861 {
862   XClearWindow (MI_DISPLAY (mi), MI_WINDOW(mi));
863   init_euler2d (mi);
864 }
865
866 ENTRYPOINT void
867 release_euler2d (ModeInfo * mi)
868 {
869         if (euler2ds != NULL) {
870                 int         screen;
871
872                 for (screen = 0; screen < MI_NUM_SCREENS(mi); screen++)
873                         free_euler2d(&euler2ds[screen]);
874                 (void) free((void *) euler2ds);
875                 euler2ds = (euler2dstruct *) NULL;
876         }
877 }
878
879 ENTRYPOINT void
880 refresh_euler2d (ModeInfo * mi)
881 {
882         MI_CLEARWINDOW(mi);
883 }
884
885 ENTRYPOINT Bool
886 euler2d_handle_event (ModeInfo *mi, XEvent *event)
887 {
888   if (screenhack_event_helper (MI_DISPLAY(mi), MI_WINDOW(mi), event))
889     {
890       init_euler2d (mi);
891       return True;
892     }
893   return False;
894 }
895
896
897
898 XSCREENSAVER_MODULE ("Euler2D", euler2d)
899
900 #endif /* MODE_euler2d */