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