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