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   {(char* ) "-eulertail", (char *) ".euler2d.eulertail",
  72    XrmoptionSepArg, (caddr_t) NULL},
  73   {(char* ) "-eulerpower", (char *) ".euler2d.eulerpower",
  74    XrmoptionSepArg, (caddr_t) NULL},
  75 };
  76 static argtype vars[] =
  77 {
  78   {&tail_len, "eulertail",
  79    (char *) "EulerTail", (char *) DEF_EULERTAIL, t_Int},
  80   {&power, "eulerpower",
  81    (char *) "EulerPower", (char *) "1", t_Float},
  82 };
  83 static OptionStruct desc[] =
  84 {
  85   {(char *) "-eulertail len", (char *) "Length of Euler2d tails"},
  86   {(char *) "-eulerpower power", (char *) "power of interaction law for points for Euler2d"},
  87 };
  88
  89 ModeSpecOpt euler2d_opts =
  90 {sizeof opts / sizeof opts[0], opts,
  91  sizeof vars / sizeof vars[0], vars, desc};
  92
  93 #ifdef USE_MODULES
  94 ModStruct   euler2d_description = {
  95         "euler2d", "init_euler2d", "draw_euler2d", "release_euler2d",
  96         "refresh_euler2d", "init_euler2d", (char *) NULL, &euler2d_opts,
  97         1000, 1024, 3000, 1, 64, 1.0, "",
  98         "Simulates 2D incompressible invisid fluid.", 0, NULL
  99 };
 100
 101 #endif
 102
 103 #define balance_rand(v) ((LRAND()/MAXRAND*(v))-((v)/2)) /* random around 0 */
 104 #define positive_rand(v)        (LRAND()/MAXRAND*(v))   /* positive random */
 105
 106 #define number_of_vortex_points 20
 107
 108 #define n_bound_p 500
 109 #define deg_p 6
 110
 111 static double delta_t;
 112
 113 typedef struct {
 114         int        width;
 115         int        height;
 116         int        count;
 117         double     xshift,yshift,scale;
 118         double     radius;
 119
 120         int        N;
 121         int        Nvortex;
 122
 123 /*  x[2i+0] = x coord for nth point
 124     x[2i+1] = y coord for nth point
 125     w[i] = vorticity at nth point
 126 */
 127         double     *x;
 128         double     *w;
 129
 130         double     *diffx;
 131         double     *olddiffx;
 132         double     *tempx;
 133         double     *tempdiffx;
 134 /*  (xs[2i+0],xs[2i+1]) is reflection of (x[2i+0],x[2i+1]) about unit circle
 135     xs[2i+0] = x[2i+0]/nx
 136     xs[2i+1] = x[2i+1]/nx
 137     where
 138     nx = x[2i+0]*x[2i+0] + x[2i+1]*x[2i+1]
 139
 140     x_is_zero[i] = (nx < 1e-10)
 141 */
 142         double     *xs;
 143         short      *x_is_zero;
 144
 145 /*  (p[2i+0],p[2i+1]) is image of (x[2i+0],x[2i+1]) under polynomial p.
 146     mod_dp2 is |p'(z)|^2 when z = (x[2i+0],x[2i+1]).
 147 */
 148         double     *p;
 149         double     *mod_dp2;
 150
 151 /* Sometimes in our calculations we get overflow or numbers that are too big.
 152    If that happens with the point x[2*i+0], x[2*i+1], we set dead[i].
 153 */
 154         short      *dead;
 155
 156         XSegment   *csegs;
 157         int         cnsegs;
 158         XSegment   *old_segs;
 159         int        *nold_segs;
 160         int         c_old_seg;
 161         int         boundary_color;
 162         int         hide_vortex;
 163         short      *lastx;
 164
 165         double      p_coef[2*(deg_p-1)];
 166         XSegment   *boundary;
 167
 168 } euler2dstruct;
 169
 170 static euler2dstruct *euler2ds = (euler2dstruct *) NULL;
 171
 172 /*
 173   If variable_boundary == 1, then we make a variable boundary.
 174   The way this is done is to map the unit disk under a
 175   polynomial p, where
 176   p(z) = z + c_2 z^2 + ... + c_n z^n
 177   where n = deg_p.  sp->p_coef contains the complex numbers
 178   c_2, c_3, ... c_n.
 179 */
 180
 181 #define add(a1,a2,b1,b2) (a1)+=(b1);(a2)+=(b2)
 182 #define mult(a1,a2,b1,b2) temp=(a1)*(b1)-(a2)*(b2); \
 183                           (a2)=(a1)*(b2)+(a2)*(b1);(a1)=temp
 184
 185 static void
 186 calc_p(double *p1, double *p2, double z1, double z2, double p_coef[])
 187 {
 188   int i;
 189   double temp;
 190
 191   *p1=0;
 192   *p2=0;
 193   for(i=deg_p;i>=2;i--)
 194   {
 195     add(*p1,*p2,p_coef[(i-2)*2],p_coef[(i-2)*2+1]);
 196     mult(*p1,*p2,z1,z2);
 197   }
 198   add(*p1,*p2,1,0);
 199   mult(*p1,*p2,z1,z2);
 200 }
 201
 202 /* Calculate |p'(z)|^2 */
 203 static double
 204 calc_mod_dp2(double z1, double z2, double p_coef[])
 205 {
 206   int i;
 207   double temp,mp1,mp2;
 208
 209   mp1=0;
 210   mp2=0;
 211   for(i=deg_p;i>=2;i--)
 212   {
 213     add(mp1,mp2,i*p_coef[(i-2)*2],i*p_coef[(i-2)*2+1]);
 214     mult(mp1,mp2,z1,z2);
 215   }
 216   add(mp1,mp2,1,0);
 217   return mp1*mp1+mp2*mp2;
 218 }
 219
 220 static void
 221 calc_all_p(euler2dstruct *sp)
 222 {
 223   int i,j;
 224   double temp,p1,p2,z1,z2;
 225   for(j=(sp->hide_vortex?sp->Nvortex:0);j<sp->N;j++) if(!sp->dead[j])
 226   {
 227     p1=0;
 228     p2=0;
 229     z1=sp->x[2*j+0];
 230     z2=sp->x[2*j+1];
 231     for(i=deg_p;i>=2;i--)
 232     {
 233       add(p1,p2,sp->p_coef[(i-2)*2],sp->p_coef[(i-2)*2+1]);
 234       mult(p1,p2,z1,z2);
 235     }
 236     add(p1,p2,1,0);
 237     mult(p1,p2,z1,z2);
 238     sp->p[2*j+0] = p1;
 239     sp->p[2*j+1] = p2;
 240   }
 241 }
 242
 243 static void
 244 calc_all_mod_dp2(double *x, euler2dstruct *sp)
 245 {
 246   int i,j;
 247   double temp,mp1,mp2,z1,z2;
 248   for(j=0;j<sp->N;j++) if(!sp->dead[j])
 249   {
 250     mp1=0;
 251     mp2=0;
 252     z1=x[2*j+0];
 253     z2=x[2*j+1];
 254     for(i=deg_p;i>=2;i--)
 255     {
 256       add(mp1,mp2,i*sp->p_coef[(i-2)*2],i*sp->p_coef[(i-2)*2+1]);
 257       mult(mp1,mp2,z1,z2);
 258     }
 259     add(mp1,mp2,1,0);
 260     sp->mod_dp2[j] = mp1*mp1+mp2*mp2;
 261   }
 262 }
 263
 264 static void
 265 derivs(double *x, euler2dstruct *sp)
 266 {
 267   int i,j;
 268   double u1,u2,x1,x2,xij1,xij2,nxij;
 269   double nx;
 270
 271   if (variable_boundary)
 272     calc_all_mod_dp2(sp->x,sp);
 273
 274   for (j=0;j<sp->Nvortex;j++) if (!sp->dead[j])
 275   {
 276     nx = x[2*j+0]*x[2*j+0] + x[2*j+1]*x[2*j+1];
 277     if (nx < 1e-10)
 278       sp->x_is_zero[j] = 1;
 279     else {
 280       sp->x_is_zero[j] = 0;
 281       sp->xs[2*j+0] = x[2*j+0]/nx;
 282       sp->xs[2*j+1] = x[2*j+1]/nx;
 283     }
 284   }
 285
 286   (void) memset(sp->diffx,0,sizeof(double)*2*sp->N);
 287
 288   for (i=0;i<sp->N;i++) if (!sp->dead[i])
 289   {
 290     x1 = x[2*i+0];
 291     x2 = x[2*i+1];
 292     for (j=0;j<sp->Nvortex;j++) if (!sp->dead[j])
 293     {
 294 /*
 295   Calculate the Biot-Savart kernel, that is, effect of a
 296   vortex point at a = (x[2*j+0],x[2*j+1]) at the point
 297   x = (x1,x2), returning the vector field in (u1,u2).
 298
 299   In the plane, this is given by the formula
 300
 301   u = (x-a)/|x-a|^2  or zero if x=a.
 302
 303   However, in the unit disk we have to subtract from the
 304   above:
 305
 306   (x-as)/|x-as|^2
 307
 308   where as = a/|a|^2 is the reflection of a about the unit circle.
 309
 310   If however power != 1, then
 311
 312   u = (x-a)/|x-a|^(power+1)  -  |a|^(1-power) (x-as)/|x-as|^(power+1)
 313
 314 */
 315
 316       xij1 = x1 - x[2*j+0];
 317       xij2 = x2 - x[2*j+1];
 318       nxij = (power==1.0) ? xij1*xij1+xij2*xij2 : pow(xij1*xij1+xij2*xij2,(power+1)/2.0);
 319
 320       if(nxij >= 1e-4)  {
 321         u1 =  xij2/nxij;
 322         u2 = -xij1/nxij;
 323       }
 324       else
 325         u1 = u2 = 0.0;
 326
 327       if (!sp->x_is_zero[j])
 328       {
 329         xij1 = x1 - sp->xs[2*j+0];
 330         xij2 = x2 - sp->xs[2*j+1];
 331         nxij = (power==1.0) ? xij1*xij1+xij2*xij2 : pow(xij1*xij1+xij2*xij2,(power+1)/2.0);
 332
 333         if (nxij < 1e-5)
 334         {
 335           sp->dead[i] = 1;
 336           u1 = u2 = 0.0;
 337         }
 338         else
 339         {
 340           u1 -= xij2/nxij;
 341           u2 += xij1/nxij;
 342         }
 343       }
 344
 345       if (!sp->dead[i])
 346       {
 347         sp->diffx[2*i+0] += u1*sp->w[j];
 348         sp->diffx[2*i+1] += u2*sp->w[j];
 349       }
 350     }
 351
 352     if (!sp->dead[i] && variable_boundary)
 353     {
 354       if (sp->mod_dp2[i] < 1e-5)
 355         sp->dead[i] = 1;
 356       else
 357       {
 358         sp->diffx[2*i+0] /= sp->mod_dp2[i];
 359         sp->diffx[2*i+1] /= sp->mod_dp2[i];
 360       }
 361     }
 362   }
 363 }
 364
 365 /*
 366   What perturb does is effectively
 367     ret = x + k,
 368   where k should be of order delta_t.
 369
 370   We have the option to do this more subtly by mapping points x
 371   in the unit disk to points y in the plane, where y = f(|x|) x,
 372   with f(t) = -log(1-t)/t.
 373
 374   This might reduce (but does not remove) problems where particles near
 375   the edge of the boundary bounce around.
 376
 377   But it seems to be not that effective, so for now switch it off.
 378 */
 379
 380 #define SUBTLE_PERTURB 0
 381
 382 static void
 383 perturb(double ret[], double x[], double k[], euler2dstruct *sp)
 384 {
 385   int i;
 386   double x1,x2,k1,k2;
 387
 388 #if SUBTLE_PERTURB
 389   double d1,d2,t1,t2,mag,mag2,mlog1mmag,memmagdmag,xdotk;
 390   for (i=0;i<sp->N;i++) if (!sp->dead[i])
 391   {
 392     x1 = x[2*i+0];
 393     x2 = x[2*i+1];
 394     k1 = k[2*i+0];
 395     k2 = k[2*i+1];
 396     mag2 = x1*x1 + x2*x2;
 397     if (mag2 < 1e-10)
 398     {
 399       ret[2*i+0] = x1+k1;
 400       ret[2*i+1] = x2+k2;
 401     }
 402     else if (mag2 > 1-1e-5)
 403       sp->dead[i] = 1;
 404     else
 405     {
 406       mag = sqrt(mag2);
 407       mlog1mmag = -log(1-mag);
 408       xdotk = x1*k1 + x2*k2;
 409       t1 = (x1 + k1)*mlog1mmag/mag + x1*xdotk*(1.0/(1-mag)-mlog1mmag/mag)/mag/mag;
 410       t2 = (x2 + k2)*mlog1mmag/mag + x2*xdotk*(1.0/(1-mag)-mlog1mmag/mag)/mag/mag;
 411       mag = sqrt(t1*t1+t2*t2);
 412       if (mag > 11.5 /* log(1e5) */)
 413         sp->dead[i] = 1;
 414       else
 415       {
 416         memmagdmag = (mag>1e-5) ? ((1.0-exp(-mag))/mag) : (1-mag/2.0);
 417         ret[2*i+0] = t1*memmagdmag;
 418         ret[2*i+1] = t2*memmagdmag;
 419       }
 420     }
 421     if (!sp->dead[i])
 422     {
 423       d1 = ret[2*i+0]-x1;
 424       d2 = ret[2*i+1]-x2;
 425       if (d1*d1+d2*d2 > 0.1)
 426         sp->dead[i] = 1;
 427     }
 428   }
 429
 430 #else
 431
 432   for (i=0;i<sp->N;i++) if (!sp->dead[i])
 433   {
 434     x1 = x[2*i+0];
 435     x2 = x[2*i+1];
 436     k1 = k[2*i+0];
 437     k2 = k[2*i+1];
 438     if (k1*k1+k2*k2 > 0.1 || x1*x1+x2*x2 > 1-1e-5)
 439       sp->dead[i] = 1;
 440     else
 441     {
 442       ret[2*i+0] = x1+k1;
 443       ret[2*i+1] = x2+k2;
 444     }
 445   }
 446 #endif
 447 }
 448
 449 static void
 450 ode_solve(euler2dstruct *sp)
 451 {
 452   int i;
 453   double *temp;
 454
 455   if (sp->count < 1) {
 456     /* midpoint method */
 457     derivs(sp->x,sp);
 458     (void) memcpy(sp->olddiffx,sp->diffx,sizeof(double)*2*sp->N);
 459     for (i=0;i<sp->N;i++) if (!sp->dead[i]) {
 460       sp->tempdiffx[2*i+0] = 0.5*delta_t*sp->diffx[2*i+0];
 461       sp->tempdiffx[2*i+1] = 0.5*delta_t*sp->diffx[2*i+1];
 462     }
 463     perturb(sp->tempx,sp->x,sp->tempdiffx,sp);
 464     derivs(sp->tempx,sp);
 465     for (i=0;i<sp->N;i++) if (!sp->dead[i]) {
 466       sp->tempdiffx[2*i+0] = delta_t*sp->diffx[2*i+0];
 467       sp->tempdiffx[2*i+1] = delta_t*sp->diffx[2*i+1];
 468     }
 469     perturb(sp->x,sp->x,sp->tempdiffx,sp);
 470   } else {
 471     /* Adams Basforth */
 472     derivs(sp->x,sp);
 473     for (i=0;i<sp->N;i++) if (!sp->dead[i]) {
 474       sp->tempdiffx[2*i+0] = delta_t*(1.5*sp->diffx[2*i+0] - 0.5*sp->olddiffx[2*i+0]);
 475       sp->tempdiffx[2*i+1] = delta_t*(1.5*sp->diffx[2*i+1] - 0.5*sp->olddiffx[2*i+1]);
 476     }
 477     perturb(sp->x,sp->x,sp->tempdiffx,sp);
 478     temp = sp->olddiffx;
 479     sp->olddiffx = sp->diffx;
 480     sp->diffx = temp;
 481   }
 482 }
 483
 484 #define deallocate(p,t) if (p!=NULL) {(void) free((void *) p); p=(t*)NULL; }
 485 #define allocate(p,t,s) if ((p=(t*)malloc(sizeof(t)*s))==NULL)\
 486 {free_euler2d(sp);return;}
 487
 488 static void
 489 free_euler2d(euler2dstruct *sp)
 490 {
 491         deallocate(sp->csegs, XSegment);
 492         deallocate(sp->old_segs, XSegment);
 493         deallocate(sp->nold_segs, int);
 494         deallocate(sp->lastx, short);
 495         deallocate(sp->x, double);
 496         deallocate(sp->diffx, double);
 497         deallocate(sp->w, double);
 498         deallocate(sp->olddiffx, double);
 499         deallocate(sp->tempdiffx, double);
 500         deallocate(sp->tempx, double);
 501         deallocate(sp->dead, short);
 502         deallocate(sp->boundary, XSegment);
 503         deallocate(sp->xs, double);
 504         deallocate(sp->x_is_zero, short);
 505         deallocate(sp->p, double);
 506         deallocate(sp->mod_dp2, double);
 507 }
 508
 509 void
 510 init_euler2d(ModeInfo * mi)
 511 {
 512 #define nr_rotates 18 /* how many rotations to try to fill as much of screen as possible - must be even number */
 513         euler2dstruct *sp;
 514         int         i,k,n,np;
 515         double      r,theta,x,y,w;
 516         double      mag,xscale,yscale,p1,p2;
 517         double      low[nr_rotates],high[nr_rotates],pp1,pp2,pn1,pn2,angle1,angle2,tempangle,dist,scale,bestscale,temp;
 518         int         besti = 0;
 519
 520         if (power<0.5) power = 0.5;
 521         if (power>3.0) power = 3.0;
 522         variable_boundary &= power == 1.0;
 523         delta_t = 0.001;
 524         if (power>1.0) delta_t *= pow(0.1,power-1);
 525
 526         if (euler2ds == NULL) {
 527                 if ((euler2ds = (euler2dstruct *) calloc(MI_NUM_SCREENS(mi),
 528                                                sizeof (euler2dstruct))) == NULL)
 529                         return;
 530         }
 531         sp = &euler2ds[MI_SCREEN(mi)];
 532
 533         sp->boundary_color = NRAND(MI_NPIXELS(mi));
 534         sp->hide_vortex = NRAND(4) != 0;
 535
 536         sp->count = 0;
 537
 538         sp->width = MI_WIDTH(mi);
 539         sp->height = MI_HEIGHT(mi);
 540
 541         sp->N = MI_COUNT(mi)+number_of_vortex_points;
 542         sp->Nvortex = number_of_vortex_points;
 543
 544         if (tail_len < 1) { /* minimum tail */
 545           tail_len = 1;
 546         }
 547         if (tail_len > MI_CYCLES(mi)) { /* maximum tail */
 548           tail_len = MI_CYCLES(mi);
 549         }
 550
 551         /* Clear the background. */
 552         MI_CLEARWINDOW(mi);
 553
 554         free_euler2d(sp);
 555
 556         /* Allocate memory. */
 557
 558         if (sp->csegs == NULL) {
 559                 allocate(sp->csegs, XSegment, sp->N);
 560                 allocate(sp->old_segs, XSegment, sp->N * tail_len);
 561                 allocate(sp->nold_segs, int, tail_len);
 562                 allocate(sp->lastx, short, sp->N * 2);
 563                 allocate(sp->x, double, sp->N * 2);
 564                 allocate(sp->diffx, double, sp->N * 2);
 565                 allocate(sp->w, double, sp->Nvortex);
 566                 allocate(sp->olddiffx, double, sp->N * 2);
 567                 allocate(sp->tempdiffx, double, sp->N * 2);
 568                 allocate(sp->tempx, double, sp->N * 2);
 569                 allocate(sp->dead, short, sp->N);
 570                 allocate(sp->boundary, XSegment, n_bound_p);
 571                 allocate(sp->xs, double, sp->Nvortex * 2);
 572                 allocate(sp->x_is_zero, short, sp->Nvortex);
 573                 allocate(sp->p, double, sp->N * 2);
 574                 allocate(sp->mod_dp2, double, sp->N);
 575         }
 576         for (i=0;i<tail_len;i++) {
 577                 sp->nold_segs[i] = 0;
 578         }
 579         sp->c_old_seg = 0;
 580         (void) memset(sp->dead,0,sp->N*sizeof(short));
 581
 582         if (variable_boundary)
 583         {
 584         /* Initialize polynomial p */
 585 /*
 586   The polynomial p(z) = z + c_2 z^2 + ... c_n z^n needs to be
 587   a bijection of the unit disk onto its image.  This is achieved
 588   by insisting that sum_{k=2}^n k |c_k| <= 1.  Actually we set
 589   the inequality to be equality (to get more interesting shapes).
 590 */
 591                 mag = 0;
 592                 for(k=2;k<=deg_p;k++)
 593                 {
 594                         r = positive_rand(1.0/k);
 595                         theta = balance_rand(2*M_PI);
 596                         sp->p_coef[2*(k-2)+0]=r*cos(theta);
 597                         sp->p_coef[2*(k-2)+1]=r*sin(theta);
 598                         mag += k*r;
 599                 }
 600                 if (mag > 0.0001) for(k=2;k<=deg_p;k++)
 601                 {
 602                         sp->p_coef[2*(k-2)+0] /= mag;
 603                         sp->p_coef[2*(k-2)+1] /= mag;
 604                 }
 605
 606 #if DEBUG_POINTED_REGION
 607                 for(k=2;k<=deg_p;k++){
 608                         sp->p_coef[2*(k-2)+0]=0;
 609                         sp->p_coef[2*(k-2)+1]=0;
 610                 }
 611                 sp->p_coef[2*(6-2)+0] = 1.0/6.0;
 612 #endif
 613
 614
 615 /* Here we figure out the best rotation of the domain so that it fills as
 616    much of the screen as possible.  The number of angles we look at is determined
 617    by nr_rotates (we look every 180/nr_rotates degrees).
 618    While we figure out the best angle to rotate, we also figure out the correct scaling factors.
 619 */
 620
 621                 for(k=0;k<nr_rotates;k++) {
 622                         low[k] = 1e5;
 623                         high[k] = -1e5;
 624                 }
 625
 626                 for(k=0;k<n_bound_p;k++) {
 627                         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);
 628                         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);
 629                         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);
 630                         angle1 = nr_rotates/M_PI*atan2(p2-pp2,p1-pp1)-nr_rotates/2;
 631                         angle2 = nr_rotates/M_PI*atan2(pn2-p2,pn1-p1)-nr_rotates/2;
 632                         while (angle1<0) angle1+=nr_rotates*2;
 633                         while (angle2<0) angle2+=nr_rotates*2;
 634                         if (angle1>nr_rotates*1.75 && angle2<nr_rotates*0.25) angle2+=nr_rotates*2;
 635                         if (angle1<nr_rotates*0.25 && angle2>nr_rotates*1.75) angle1+=nr_rotates*2;
 636                         if (angle2<angle1) {
 637                                 tempangle=angle1;
 638                                 angle1=angle2;
 639                                 angle2=tempangle;
 640                         }
 641                         for(i=(int)floor(angle1);i<(int)ceil(angle2);i++) {
 642                                 dist = cos((double)i*M_PI/nr_rotates)*p1 + sin((double)i*M_PI/nr_rotates)*p2;
 643                                 if (i%(nr_rotates*2)<nr_rotates) {
 644                                         if (dist>high[i%nr_rotates]) high[i%nr_rotates] = dist;
 645                                         if (dist<low[i%nr_rotates]) low[i%nr_rotates] = dist;
 646                                 }
 647                                 else {
 648                                         if (-dist>high[i%nr_rotates]) high[i%nr_rotates] = -dist;
 649                                         if (-dist<low[i%nr_rotates]) low[i%nr_rotates] = -dist;
 650                                 }
 651                         }
 652                 }
 653                 bestscale = 0;
 654                 for (i=0;i<nr_rotates;i++) {
 655                         xscale = (sp->width-5.0)/(high[i]-low[i]);
 656                         yscale = (sp->height-5.0)/(high[(i+nr_rotates/2)%nr_rotates]-low[(i+nr_rotates/2)%nr_rotates]);
 657                         scale = (xscale>yscale) ? yscale : xscale;
 658                         if (scale>bestscale) {
 659                                 bestscale = scale;
 660                                 besti = i;
 661                         }
 662                 }
 663 /* Here we do the rotation.  The way we do this is to replace the
 664    polynomial p(z) by a^{-1} p(a z) where a = exp(i best_angle).
 665 */
 666                 p1 = 1;
 667                 p2 = 0;
 668                 for(k=2;k<=deg_p;k++)
 669                 {
 670                         mult(p1,p2,cos((double)besti*M_PI/nr_rotates),sin((double)besti*M_PI/nr_rotates));
 671                         mult(sp->p_coef[2*(k-2)+0],sp->p_coef[2*(k-2)+1],p1,p2);
 672                 }
 673
 674                 sp->scale = bestscale;
 675                 sp->xshift = -(low[besti]+high[besti])/2.0*sp->scale+sp->width/2;
 676                 if (besti<nr_rotates/2)
 677                         sp->yshift = -(low[besti+nr_rotates/2]+high[besti+nr_rotates/2])/2.0*sp->scale+sp->height/2;
 678                 else
 679                         sp->yshift = (low[besti-nr_rotates/2]+high[besti-nr_rotates/2])/2.0*sp->scale+sp->height/2;
 680
 681
 682 /* Initialize boundary */
 683
 684                 for(k=0;k<n_bound_p;k++)
 685                 {
 686
 687                         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);
 688                         sp->boundary[k].x1 = (short)(p1*sp->scale+sp->xshift);
 689                         sp->boundary[k].y1 = (short)(p2*sp->scale+sp->yshift);
 690                 }
 691                 for(k=1;k<n_bound_p;k++)
 692                 {
 693                         sp->boundary[k].x2 = sp->boundary[k-1].x1;
 694                         sp->boundary[k].y2 = sp->boundary[k-1].y1;
 695                 }
 696                 sp->boundary[0].x2 = sp->boundary[n_bound_p-1].x1;
 697                 sp->boundary[0].y2 = sp->boundary[n_bound_p-1].y1;
 698         }
 699         else
 700         {
 701                 if (sp->width>sp->height)
 702                         sp->radius = sp->height/2.0-5.0;
 703                 else
 704                         sp->radius = sp->width/2.0-5.0;
 705         }
 706
 707         /* Initialize point positions */
 708
 709         for (i=sp->Nvortex;i<sp->N;i++) {
 710                 do {
 711                         r = sqrt(positive_rand(1.0));
 712                         theta = balance_rand(2*M_PI);
 713                         sp->x[2*i+0]=r*cos(theta);
 714                         sp->x[2*i+1]=r*sin(theta);
 715                 /* This is to make sure the initial distribution of points is uniform */
 716                 } while (variable_boundary &&
 717                           calc_mod_dp2(sp->x[2*i+0],sp->x[2*i+1],sp->p_coef)
 718                           < positive_rand(4));
 719         }
 720
 721         n = NRAND(4)+2;
 722         /* number of vortex points with negative vorticity */
 723         if (n%2) {
 724                 np = NRAND(n+1);
 725         }
 726         else {
 727                 /* if n is even make sure that np==n/2 is twice as likely
 728                    as the other possibilities. */
 729                 np = NRAND(n+2);
 730                 if (np==n+1) np=n/2;
 731         }
 732         for(k=0;k<n;k++)
 733         {
 734                 r = sqrt(positive_rand(0.77));
 735                 theta = balance_rand(2*M_PI);
 736                 x=r*cos(theta);
 737                 y=r*sin(theta);
 738                 r = 0.02+positive_rand(0.1);
 739                 w = (2*(k<np)-1)*2.0/sp->Nvortex;
 740                 for (i=sp->Nvortex*k/n;i<sp->Nvortex*(k+1)/n;i++) {
 741                         theta = balance_rand(2*M_PI);
 742                         sp->x[2*i+0]=x + r*cos(theta);
 743                         sp->x[2*i+1]=y + r*sin(theta);
 744                         sp->w[i]=w;
 745                 }
 746         }
 747 }
 748
 749 void
 750 draw_euler2d(ModeInfo * mi)
 751 {
 752         Display    *display = MI_DISPLAY(mi);
 753         Window      window = MI_WINDOW(mi);
 754         GC          gc = MI_GC(mi);
 755         int         b, col, n_non_vortex_segs;
 756         euler2dstruct *sp;
 757
 758         MI_IS_DRAWN(mi) = True;
 759
 760         if (euler2ds == NULL)
 761                 return;
 762         sp = &euler2ds[MI_SCREEN(mi)];
 763         if (sp->csegs == NULL)
 764                 return;
 765
 766         ode_solve(sp);
 767         if (variable_boundary)
 768                 calc_all_p(sp);
 769
 770         sp->cnsegs = 0;
 771         for(b=sp->Nvortex;b<sp->N;b++) if(!sp->dead[b])
 772         {
 773                 sp->csegs[sp->cnsegs].x1 = sp->lastx[2*b+0];
 774                 sp->csegs[sp->cnsegs].y1 = sp->lastx[2*b+1];
 775                 if (variable_boundary)
 776                 {
 777                         sp->csegs[sp->cnsegs].x2 = (short)(sp->p[2*b+0]*sp->scale+sp->xshift);
 778                         sp->csegs[sp->cnsegs].y2 = (short)(sp->p[2*b+1]*sp->scale+sp->yshift);
 779                 }
 780                 else
 781                 {
 782                         sp->csegs[sp->cnsegs].x2 = (short)(sp->x[2*b+0]*sp->radius+sp->width/2);
 783                         sp->csegs[sp->cnsegs].y2 = (short)(sp->x[2*b+1]*sp->radius+sp->height/2);
 784                 }
 785                 sp->lastx[2*b+0] = sp->csegs[sp->cnsegs].x2;
 786                 sp->lastx[2*b+1] = sp->csegs[sp->cnsegs].y2;
 787                 sp->cnsegs++;
 788         }
 789         n_non_vortex_segs = sp->cnsegs;
 790
 791         if (!sp->hide_vortex) for(b=0;b<sp->Nvortex;b++) if(!sp->dead[b])
 792         {
 793                 sp->csegs[sp->cnsegs].x1 = sp->lastx[2*b+0];
 794                 sp->csegs[sp->cnsegs].y1 = sp->lastx[2*b+1];
 795                 if (variable_boundary)
 796                 {
 797                         sp->csegs[sp->cnsegs].x2 = (short)(sp->p[2*b+0]*sp->scale+sp->xshift);
 798                         sp->csegs[sp->cnsegs].y2 = (short)(sp->p[2*b+1]*sp->scale+sp->yshift);
 799                 }
 800                 else
 801                 {
 802                         sp->csegs[sp->cnsegs].x2 = (short)(sp->x[2*b+0]*sp->radius+sp->width/2);
 803                         sp->csegs[sp->cnsegs].y2 = (short)(sp->x[2*b+1]*sp->radius+sp->height/2);
 804                 }
 805                 sp->lastx[2*b+0] = sp->csegs[sp->cnsegs].x2;
 806                 sp->lastx[2*b+1] = sp->csegs[sp->cnsegs].y2;
 807                 sp->cnsegs++;
 808         }
 809
 810         if (sp->count) {
 811                 XSetForeground(display, gc, MI_BLACK_PIXEL(mi));
 812
 813                 XDrawSegments(display, window, gc, sp->old_segs+sp->c_old_seg*sp->N, sp->nold_segs[sp->c_old_seg]);
 814
 815                 if (MI_NPIXELS(mi) > 2){ /* render colour */
 816                         for (col = 0; col < MI_NPIXELS(mi); col++) {
 817                           int start = col*n_non_vortex_segs/MI_NPIXELS(mi);
 818                           int finish = (col+1)*n_non_vortex_segs/MI_NPIXELS(mi);
 819                           XSetForeground(display, gc, MI_PIXEL(mi, col));
 820                           XDrawSegments(display, window, gc,sp->csegs+start, finish-start);
 821                         }
 822                         if (!sp->hide_vortex) {
 823                           XSetForeground(display, gc, MI_WHITE_PIXEL(mi));
 824                           XDrawSegments(display, window, gc,sp->csegs+n_non_vortex_segs, sp->cnsegs-n_non_vortex_segs);
 825                         }
 826
 827                 } else {                /* render mono */
 828                   XSetForeground(display, gc, MI_WHITE_PIXEL(mi));
 829                   XDrawSegments(display, window, gc,
 830                                                 sp->csegs, sp->cnsegs);
 831                 }
 832
 833                 if (MI_NPIXELS(mi) > 2) /* render colour */
 834                         XSetForeground(display, gc, MI_PIXEL(mi, sp->boundary_color));
 835                 else
 836                         XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_WHITE_PIXEL(mi));
 837                 if (variable_boundary)
 838                         XDrawSegments(display, window, gc,
 839                                                 sp->boundary, n_bound_p);
 840                 else
 841                         XDrawArc(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
 842                                  sp->width/2 - (int) sp->radius - 1, sp->height/2 - (int) sp->radius -1,
 843                                  (int) (2*sp->radius) + 2, (int) (2* sp->radius) + 2, 0, 64*360);
 844
 845                 /* Copy to erase-list */
 846                 (void) memcpy(sp->old_segs+sp->c_old_seg*sp->N, sp->csegs, sp->cnsegs*sizeof(XSegment));
 847                 sp->nold_segs[sp->c_old_seg] = sp->cnsegs;
 848                 sp->c_old_seg++;
 849                 if (sp->c_old_seg >= tail_len)
 850                   sp->c_old_seg = 0;
 851         }
 852
 853         if (++sp->count > MI_CYCLES(mi)) /* pick a new flow */
 854                 init_euler2d(mi);
 855
 856 }
 857
 858 void
 859 release_euler2d(ModeInfo * mi)
 860 {
 861         if (euler2ds != NULL) {
 862                 int         screen;
 863
 864                 for (screen = 0; screen < MI_NUM_SCREENS(mi); screen++)
 865                         free_euler2d(&euler2ds[screen]);
 866                 (void) free((void *) euler2ds);
 867                 euler2ds = (euler2dstruct *) NULL;
 868         }
 869 }
 870
 871 void
 872 refresh_euler2d(ModeInfo * mi)
 873 {
 874         MI_CLEARWINDOW(mi);
 875 }
 876
 877 #endif /* MODE_euler2d */