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 PROGCLASS "Euler2d"
  46 #define HACK_INIT init_euler2d
  47 #define HACK_DRAW draw_euler2d
  48 #define euler2d_opts xlockmore_opts
  49 #define DEFAULTS "*delay: 10000 \n" \
  50 "*count: 1024 \n" \
  51 "*cycles: 3000 \n" \
  52 "*ncolors: 64 \n"
  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 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 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         sp->boundary_color = NRAND(MI_NPIXELS(mi));
 532         sp->hide_vortex = NRAND(4) != 0;
 533
 534         sp->count = 0;
 535
 536         sp->width = MI_WIDTH(mi);
 537         sp->height = MI_HEIGHT(mi);
 538
 539         sp->N = MI_COUNT(mi)+number_of_vortex_points;
 540         sp->Nvortex = number_of_vortex_points;
 541
 542         if (tail_len < 1) { /* minimum tail */
 543           tail_len = 1;
 544         }
 545         if (tail_len > MI_CYCLES(mi)) { /* maximum tail */
 546           tail_len = MI_CYCLES(mi);
 547         }
 548
 549         /* Clear the background. */
 550         MI_CLEARWINDOW(mi);
 551
 552         free_euler2d(sp);
 553
 554         /* Allocate memory. */
 555
 556         if (sp->csegs == NULL) {
 557                 allocate(sp->csegs, XSegment, sp->N);
 558                 allocate(sp->old_segs, XSegment, sp->N * tail_len);
 559                 allocate(sp->nold_segs, int, tail_len);
 560                 allocate(sp->lastx, short, sp->N * 2);
 561                 allocate(sp->x, double, sp->N * 2);
 562                 allocate(sp->diffx, double, sp->N * 2);
 563                 allocate(sp->w, double, sp->Nvortex);
 564                 allocate(sp->olddiffx, double, sp->N * 2);
 565                 allocate(sp->tempdiffx, double, sp->N * 2);
 566                 allocate(sp->tempx, double, sp->N * 2);
 567                 allocate(sp->dead, short, sp->N);
 568                 allocate(sp->boundary, XSegment, n_bound_p);
 569                 allocate(sp->xs, double, sp->Nvortex * 2);
 570                 allocate(sp->x_is_zero, short, sp->Nvortex);
 571                 allocate(sp->p, double, sp->N * 2);
 572                 allocate(sp->mod_dp2, double, sp->N);
 573         }
 574         for (i=0;i<tail_len;i++) {
 575                 sp->nold_segs[i] = 0;
 576         }
 577         sp->c_old_seg = 0;
 578         (void) memset(sp->dead,0,sp->N*sizeof(short));
 579
 580         if (variable_boundary)
 581         {
 582         /* Initialize polynomial p */
 583 /*
 584   The polynomial p(z) = z + c_2 z^2 + ... c_n z^n needs to be
 585   a bijection of the unit disk onto its image.  This is achieved
 586   by insisting that sum_{k=2}^n k |c_k| <= 1.  Actually we set
 587   the inequality to be equality (to get more interesting shapes).
 588 */
 589                 mag = 0;
 590                 for(k=2;k<=deg_p;k++)
 591                 {
 592                         r = positive_rand(1.0/k);
 593                         theta = balance_rand(2*M_PI);
 594                         sp->p_coef[2*(k-2)+0]=r*cos(theta);
 595                         sp->p_coef[2*(k-2)+1]=r*sin(theta);
 596                         mag += k*r;
 597                 }
 598                 if (mag > 0.0001) for(k=2;k<=deg_p;k++)
 599                 {
 600                         sp->p_coef[2*(k-2)+0] /= mag;
 601                         sp->p_coef[2*(k-2)+1] /= mag;
 602                 }
 603
 604 #if DEBUG_POINTED_REGION
 605                 for(k=2;k<=deg_p;k++){
 606                         sp->p_coef[2*(k-2)+0]=0;
 607                         sp->p_coef[2*(k-2)+1]=0;
 608                 }
 609                 sp->p_coef[2*(6-2)+0] = 1.0/6.0;
 610 #endif
 611
 612
 613 /* Here we figure out the best rotation of the domain so that it fills as
 614    much of the screen as possible.  The number of angles we look at is determined
 615    by nr_rotates (we look every 180/nr_rotates degrees).
 616    While we figure out the best angle to rotate, we also figure out the correct scaling factors.
 617 */
 618
 619                 for(k=0;k<nr_rotates;k++) {
 620                         low[k] = 1e5;
 621                         high[k] = -1e5;
 622                 }
 623
 624                 for(k=0;k<n_bound_p;k++) {
 625                         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);
 626                         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);
 627                         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);
 628                         angle1 = nr_rotates/M_PI*atan2(p2-pp2,p1-pp1)-nr_rotates/2;
 629                         angle2 = nr_rotates/M_PI*atan2(pn2-p2,pn1-p1)-nr_rotates/2;
 630                         while (angle1<0) angle1+=nr_rotates*2;
 631                         while (angle2<0) angle2+=nr_rotates*2;
 632                         if (angle1>nr_rotates*1.75 && angle2<nr_rotates*0.25) angle2+=nr_rotates*2;
 633                         if (angle1<nr_rotates*0.25 && angle2>nr_rotates*1.75) angle1+=nr_rotates*2;
 634                         if (angle2<angle1) {
 635                                 tempangle=angle1;
 636                                 angle1=angle2;
 637                                 angle2=tempangle;
 638                         }
 639                         for(i=(int)floor(angle1);i<(int)ceil(angle2);i++) {
 640                                 dist = cos((double)i*M_PI/nr_rotates)*p1 + sin((double)i*M_PI/nr_rotates)*p2;
 641                                 if (i%(nr_rotates*2)<nr_rotates) {
 642                                         if (dist>high[i%nr_rotates]) high[i%nr_rotates] = dist;
 643                                         if (dist<low[i%nr_rotates]) low[i%nr_rotates] = dist;
 644                                 }
 645                                 else {
 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                         }
 650                 }
 651                 bestscale = 0;
 652                 for (i=0;i<nr_rotates;i++) {
 653                         xscale = (sp->width-5.0)/(high[i]-low[i]);
 654                         yscale = (sp->height-5.0)/(high[(i+nr_rotates/2)%nr_rotates]-low[(i+nr_rotates/2)%nr_rotates]);
 655                         scale = (xscale>yscale) ? yscale : xscale;
 656                         if (scale>bestscale) {
 657                                 bestscale = scale;
 658                                 besti = i;
 659                         }
 660                 }
 661 /* Here we do the rotation.  The way we do this is to replace the
 662    polynomial p(z) by a^{-1} p(a z) where a = exp(i best_angle).
 663 */
 664                 p1 = 1;
 665                 p2 = 0;
 666                 for(k=2;k<=deg_p;k++)
 667                 {
 668                         mult(p1,p2,cos((double)besti*M_PI/nr_rotates),sin((double)besti*M_PI/nr_rotates));
 669                         mult(sp->p_coef[2*(k-2)+0],sp->p_coef[2*(k-2)+1],p1,p2);
 670                 }
 671
 672                 sp->scale = bestscale;
 673                 sp->xshift = -(low[besti]+high[besti])/2.0*sp->scale+sp->width/2;
 674                 if (besti<nr_rotates/2)
 675                         sp->yshift = -(low[besti+nr_rotates/2]+high[besti+nr_rotates/2])/2.0*sp->scale+sp->height/2;
 676                 else
 677                         sp->yshift = (low[besti-nr_rotates/2]+high[besti-nr_rotates/2])/2.0*sp->scale+sp->height/2;
 678
 679
 680 /* Initialize boundary */
 681
 682                 for(k=0;k<n_bound_p;k++)
 683                 {
 684
 685                         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);
 686                         sp->boundary[k].x1 = (short)(p1*sp->scale+sp->xshift);
 687                         sp->boundary[k].y1 = (short)(p2*sp->scale+sp->yshift);
 688                 }
 689                 for(k=1;k<n_bound_p;k++)
 690                 {
 691                         sp->boundary[k].x2 = sp->boundary[k-1].x1;
 692                         sp->boundary[k].y2 = sp->boundary[k-1].y1;
 693                 }
 694                 sp->boundary[0].x2 = sp->boundary[n_bound_p-1].x1;
 695                 sp->boundary[0].y2 = sp->boundary[n_bound_p-1].y1;
 696         }
 697         else
 698         {
 699                 if (sp->width>sp->height)
 700                         sp->radius = sp->height/2.0-5.0;
 701                 else
 702                         sp->radius = sp->width/2.0-5.0;
 703         }
 704
 705         /* Initialize point positions */
 706
 707         for (i=sp->Nvortex;i<sp->N;i++) {
 708                 do {
 709                         r = sqrt(positive_rand(1.0));
 710                         theta = balance_rand(2*M_PI);
 711                         sp->x[2*i+0]=r*cos(theta);
 712                         sp->x[2*i+1]=r*sin(theta);
 713                 /* This is to make sure the initial distribution of points is uniform */
 714                 } while (variable_boundary &&
 715                           calc_mod_dp2(sp->x[2*i+0],sp->x[2*i+1],sp->p_coef)
 716                           < positive_rand(4));
 717         }
 718
 719         n = NRAND(4)+2;
 720         /* number of vortex points with negative vorticity */
 721         if (n%2) {
 722                 np = NRAND(n+1);
 723         }
 724         else {
 725                 /* if n is even make sure that np==n/2 is twice as likely
 726                    as the other possibilities. */
 727                 np = NRAND(n+2);
 728                 if (np==n+1) np=n/2;
 729         }
 730         for(k=0;k<n;k++)
 731         {
 732                 r = sqrt(positive_rand(0.77));
 733                 theta = balance_rand(2*M_PI);
 734                 x=r*cos(theta);
 735                 y=r*sin(theta);
 736                 r = 0.02+positive_rand(0.1);
 737                 w = (2*(k<np)-1)*2.0/sp->Nvortex;
 738                 for (i=sp->Nvortex*k/n;i<sp->Nvortex*(k+1)/n;i++) {
 739                         theta = balance_rand(2*M_PI);
 740                         sp->x[2*i+0]=x + r*cos(theta);
 741                         sp->x[2*i+1]=y + r*sin(theta);
 742                         sp->w[i]=w;
 743                 }
 744         }
 745 }
 746
 747 void
 748 draw_euler2d(ModeInfo * mi)
 749 {
 750         Display    *display = MI_DISPLAY(mi);
 751         Window      window = MI_WINDOW(mi);
 752         GC          gc = MI_GC(mi);
 753         int         b, col, n_non_vortex_segs;
 754         euler2dstruct *sp;
 755
 756         MI_IS_DRAWN(mi) = True;
 757
 758         if (euler2ds == NULL)
 759                 return;
 760         sp = &euler2ds[MI_SCREEN(mi)];
 761         if (sp->csegs == NULL)
 762                 return;
 763
 764         ode_solve(sp);
 765         if (variable_boundary)
 766                 calc_all_p(sp);
 767
 768         sp->cnsegs = 0;
 769         for(b=sp->Nvortex;b<sp->N;b++) if(!sp->dead[b])
 770         {
 771                 sp->csegs[sp->cnsegs].x1 = sp->lastx[2*b+0];
 772                 sp->csegs[sp->cnsegs].y1 = sp->lastx[2*b+1];
 773                 if (variable_boundary)
 774                 {
 775                         sp->csegs[sp->cnsegs].x2 = (short)(sp->p[2*b+0]*sp->scale+sp->xshift);
 776                         sp->csegs[sp->cnsegs].y2 = (short)(sp->p[2*b+1]*sp->scale+sp->yshift);
 777                 }
 778                 else
 779                 {
 780                         sp->csegs[sp->cnsegs].x2 = (short)(sp->x[2*b+0]*sp->radius+sp->width/2);
 781                         sp->csegs[sp->cnsegs].y2 = (short)(sp->x[2*b+1]*sp->radius+sp->height/2);
 782                 }
 783                 sp->lastx[2*b+0] = sp->csegs[sp->cnsegs].x2;
 784                 sp->lastx[2*b+1] = sp->csegs[sp->cnsegs].y2;
 785                 sp->cnsegs++;
 786         }
 787         n_non_vortex_segs = sp->cnsegs;
 788
 789         if (!sp->hide_vortex) for(b=0;b<sp->Nvortex;b++) if(!sp->dead[b])
 790         {
 791                 sp->csegs[sp->cnsegs].x1 = sp->lastx[2*b+0];
 792                 sp->csegs[sp->cnsegs].y1 = sp->lastx[2*b+1];
 793                 if (variable_boundary)
 794                 {
 795                         sp->csegs[sp->cnsegs].x2 = (short)(sp->p[2*b+0]*sp->scale+sp->xshift);
 796                         sp->csegs[sp->cnsegs].y2 = (short)(sp->p[2*b+1]*sp->scale+sp->yshift);
 797                 }
 798                 else
 799                 {
 800                         sp->csegs[sp->cnsegs].x2 = (short)(sp->x[2*b+0]*sp->radius+sp->width/2);
 801                         sp->csegs[sp->cnsegs].y2 = (short)(sp->x[2*b+1]*sp->radius+sp->height/2);
 802                 }
 803                 sp->lastx[2*b+0] = sp->csegs[sp->cnsegs].x2;
 804                 sp->lastx[2*b+1] = sp->csegs[sp->cnsegs].y2;
 805                 sp->cnsegs++;
 806         }
 807
 808         if (sp->count) {
 809                 XSetForeground(display, gc, MI_BLACK_PIXEL(mi));
 810
 811                 XDrawSegments(display, window, gc, sp->old_segs+sp->c_old_seg*sp->N, sp->nold_segs[sp->c_old_seg]);
 812
 813                 if (MI_NPIXELS(mi) > 2){ /* render colour */
 814                         for (col = 0; col < MI_NPIXELS(mi); col++) {
 815                           int start = col*n_non_vortex_segs/MI_NPIXELS(mi);
 816                           int finish = (col+1)*n_non_vortex_segs/MI_NPIXELS(mi);
 817                           XSetForeground(display, gc, MI_PIXEL(mi, col));
 818                           XDrawSegments(display, window, gc,sp->csegs+start, finish-start);
 819                         }
 820                         if (!sp->hide_vortex) {
 821                           XSetForeground(display, gc, MI_WHITE_PIXEL(mi));
 822                           XDrawSegments(display, window, gc,sp->csegs+n_non_vortex_segs, sp->cnsegs-n_non_vortex_segs);
 823                         }
 824
 825                 } else {                /* render mono */
 826                   XSetForeground(display, gc, MI_WHITE_PIXEL(mi));
 827                   XDrawSegments(display, window, gc,
 828                                                 sp->csegs, sp->cnsegs);
 829                 }
 830
 831                 if (MI_NPIXELS(mi) > 2) /* render colour */
 832                         XSetForeground(display, gc, MI_PIXEL(mi, sp->boundary_color));
 833                 else
 834                         XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_WHITE_PIXEL(mi));
 835                 if (variable_boundary)
 836                         XDrawSegments(display, window, gc,
 837                                                 sp->boundary, n_bound_p);
 838                 else
 839                         XDrawArc(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
 840                                  sp->width/2 - (int) sp->radius - 1, sp->height/2 - (int) sp->radius -1,
 841                                  (int) (2*sp->radius) + 2, (int) (2* sp->radius) + 2, 0, 64*360);
 842
 843                 /* Copy to erase-list */
 844                 (void) memcpy(sp->old_segs+sp->c_old_seg*sp->N, sp->csegs, sp->cnsegs*sizeof(XSegment));
 845                 sp->nold_segs[sp->c_old_seg] = sp->cnsegs;
 846                 sp->c_old_seg++;
 847                 if (sp->c_old_seg >= tail_len)
 848                   sp->c_old_seg = 0;
 849         }
 850
 851         if (++sp->count > MI_CYCLES(mi)) /* pick a new flow */
 852                 init_euler2d(mi);
 853
 854 }
 855
 856 void
 857 release_euler2d(ModeInfo * mi)
 858 {
 859         if (euler2ds != NULL) {
 860                 int         screen;
 861
 862                 for (screen = 0; screen < MI_NUM_SCREENS(mi); screen++)
 863                         free_euler2d(&euler2ds[screen]);
 864                 (void) free((void *) euler2ds);
 865                 euler2ds = (euler2dstruct *) NULL;
 866         }
 867 }
 868
 869 void
 870 refresh_euler2d(ModeInfo * mi)
 871 {
 872         MI_CLEARWINDOW(mi);
 873 }
 874
 875 #endif /* MODE_euler2d */