git.hungrycats.org Git - xscreensaver/blob - 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 */