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 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 reshape_euler2d(ModeInfo * mi, int width, int height)
 862 {
 863   XClearWindow (MI_DISPLAY (mi), MI_WINDOW(mi));
 864   init_euler2d (mi);
 865 }
 866
 867 ENTRYPOINT void
 868 release_euler2d (ModeInfo * mi)
 869 {
 870         if (euler2ds != NULL) {
 871                 int         screen;
 872
 873                 for (screen = 0; screen < MI_NUM_SCREENS(mi); screen++)
 874                         free_euler2d(&euler2ds[screen]);
 875                 (void) free((void *) euler2ds);
 876                 euler2ds = (euler2dstruct *) NULL;
 877         }
 878 }
 879
 880 ENTRYPOINT void
 881 refresh_euler2d (ModeInfo * mi)
 882 {
 883         MI_CLEARWINDOW(mi);
 884 }
 885
 886 XSCREENSAVER_MODULE ("Euler2D", euler2d)
 887
 888 #endif /* MODE_euler2d */