git.hungrycats.org Git - xscreensaver/blob - hacks/euler2d.c

   1 /* -*- Mode: C; tab-width: 4 -*- */
   2 /* euler2d --- 2 Dimensional Incompressible Inviscid Fluid Flow */
   3
   4 #if 0
   5 static const char sccsid[] = "@(#)euler2d.c     5.00 2000/11/01 xlockmore";
   6 #endif
   7
   8 /*
   9  * Copyright (c) 2000 by Stephen Montgomery-Smith <stephen@math.missouri.edu>
  10  *
  11  * Permission to use, copy, modify, and distribute this software and its
  12  * documentation for any purpose and without fee is hereby granted,
  13  * provided that the above copyright notice appear in all copies and that
  14  * both that copyright notice and this permission notice appear in
  15  * supporting documentation.
  16  *
  17  * This file is provided AS IS with no warranties of any kind.  The author
  18  * shall have no liability with respect to the infringement of copyrights,
  19  * trade secrets or any patents by this file or any part thereof.  In no
  20  * event will the author be liable for any lost revenue or profits or
  21  * other special, indirect and consequential damages.
  22  *
  23  * Revision History:
  24  * 04-Nov-2000: Added an option eulerpower.  This allows for example the
  25  *              quasi-geostrophic equation by setting eulerpower to 2.
  26  * 01-Nov-2000: Allocation checks.
  27  * 10-Sep-2000: Added optimizations, and removed subtle_perturb, by stephen.
  28  * 03-Sep-2000: Changed method of solving ode to Adams-Bashforth of order 2.
  29  *              Previously used a rather compilcated method of order 4.
  30  *              This doubles the speed of the program.  Also it seems
  31  *              to have improved numerical stability.  Done by stephen.
  32  * 27-Aug-2000: Added rotation of region to maximize screen fill by stephen.
  33  * 05-Jun-2000: Adapted from flow.c Copyright (c) 1996 by Tim Auckland
  34  * 18-Jul-1996: Adapted from swarm.c Copyright (c) 1991 by Patrick J. Naughton.
  35  * 31-Aug-1990: Adapted from xswarm by Jeff Butterworth. (butterwo@ncsc.org)
  36  */
  37
  38 /*
  39  * The mathematical aspects of this program are discussed in the file
  40  * euler2d.tex.
  41  */
  42
  43 #ifdef STANDALONE
  44 # define MODE_euler2d
  45 # define DEFAULTS       "*delay:   10000 \n" \
  46                                         "*count:   1024  \n" \
  47                                         "*cycles:  3000  \n" \
  48                                         "*ncolors: 64    \n" \
  49                                         "*fpsSolid: true    \n" \
  50                                         "*ignoreRotation: True \n" \
  51
  52 # define SMOOTH_COLORS
  53 # define release_euler2d 0
  54 # define reshape_euler2d 0
  55 # define euler2d_handle_event 0
  56 # include "xlockmore.h"         /* in xscreensaver distribution */
  57 #else /* STANDALONE */
  58 # include "xlock.h"             /* in xlockmore distribution */
  59 #endif /* STANDALONE */
  60
  61 #ifdef MODE_euler2d
  62
  63 #define DEF_EULERTAIL "10"
  64
  65 #define DEBUG_POINTED_REGION    0
  66
  67 static int tail_len;
  68 static int variable_boundary = 1;
  69 static float power = 1;
  70
  71 static XrmOptionDescRec opts[] =
  72 {
  73   {"-eulertail", ".euler2d.eulertail",   XrmoptionSepArg, NULL},
  74   {"-eulerpower", ".euler2d.eulerpower", XrmoptionSepArg, NULL},
  75 };
  76 static argtype vars[] =
  77 {
  78   {&tail_len, "eulertail",
  79    "EulerTail", (char *) DEF_EULERTAIL, t_Int},
  80   {&power, "eulerpower",
  81    "EulerPower", "1", t_Float},
  82 };
  83 static OptionStruct desc[] =
  84 {
  85   {"-eulertail len", "Length of Euler2d tails"},
  86   {"-eulerpower power", "power of interaction law for points for Euler2d"},
  87 };
  88
  89 ENTRYPOINT 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", (char *) NULL,
  96         "refresh_euler2d", "init_euler2d", "free_euler2d", &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     xshift2,yshift2;
 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(mi);return;}
 488
 489 ENTRYPOINT void
 490 free_euler2d(ModeInfo * mi)
 491 {
 492         euler2dstruct *sp = &euler2ds[MI_SCREEN(mi)];
 493         deallocate(sp->csegs, XSegment);
 494         deallocate(sp->old_segs, XSegment);
 495         deallocate(sp->nold_segs, int);
 496         deallocate(sp->lastx, short);
 497         deallocate(sp->x, double);
 498         deallocate(sp->diffx, double);
 499         deallocate(sp->w, double);
 500         deallocate(sp->olddiffx, double);
 501         deallocate(sp->tempdiffx, double);
 502         deallocate(sp->tempx, double);
 503         deallocate(sp->dead, short);
 504         deallocate(sp->boundary, XSegment);
 505         deallocate(sp->xs, double);
 506         deallocate(sp->x_is_zero, short);
 507         deallocate(sp->p, double);
 508         deallocate(sp->mod_dp2, double);
 509 }
 510
 511 ENTRYPOINT void
 512 init_euler2d (ModeInfo * mi)
 513 {
 514 #define nr_rotates 18 /* how many rotations to try to fill as much of screen as possible - must be even number */
 515         euler2dstruct *sp;
 516         int         i,k,n,np;
 517         double      r,theta,x,y,w;
 518         double      mag,xscale,yscale,p1,p2;
 519         double      low[nr_rotates],high[nr_rotates],pp1,pp2,pn1,pn2,angle1,angle2,tempangle,dist,scale,bestscale,temp;
 520         int         besti = 0;
 521
 522         if (power<0.5) power = 0.5;
 523         if (power>3.0) power = 3.0;
 524         variable_boundary &= power == 1.0;
 525         delta_t = 0.001;
 526         if (power>1.0) delta_t *= pow(0.1,power-1);
 527
 528         MI_INIT (mi, euler2ds);
 529         sp = &euler2ds[MI_SCREEN(mi)];
 530
 531 #ifdef HAVE_JWXYZ
 532   jwxyz_XSetAntiAliasing (MI_DISPLAY(mi), MI_GC(mi),  False);
 533 #endif
 534
 535         sp->boundary_color = NRAND(MI_NPIXELS(mi));
 536         sp->hide_vortex = NRAND(4) != 0;
 537
 538         sp->count = 0;
 539
 540     sp->xshift2 = sp->yshift2 = 0;
 541
 542         sp->width = MI_WIDTH(mi);
 543         sp->height = MI_HEIGHT(mi);
 544
 545     if (sp->width > sp->height * 5 ||  /* window has weird aspect */
 546         sp->height > sp->width * 5)
 547     {
 548       if (sp->width > sp->height)
 549         {
 550           sp->height = sp->width * 0.8;
 551           sp->yshift2 = -sp->height/2;
 552         }
 553       else
 554         {
 555           sp->width = sp->height * 0.8;
 556           sp->xshift2 = -sp->width/2;
 557         }
 558     }
 559
 560     if (sp->width > 2560 || sp->height > 2560)  /* Retina displays */
 561       XSetLineAttributes (MI_DISPLAY(mi), MI_GC(mi),
 562                           3, LineSolid, CapRound, JoinRound);
 563
 564
 565         sp->N = MI_COUNT(mi)+number_of_vortex_points;
 566         sp->Nvortex = number_of_vortex_points;
 567
 568         if (tail_len < 1) { /* minimum tail */
 569           tail_len = 1;
 570         }
 571         if (tail_len > MI_CYCLES(mi)) { /* maximum tail */
 572           tail_len = MI_CYCLES(mi);
 573         }
 574
 575         /* Clear the background. */
 576         MI_CLEARWINDOW(mi);
 577
 578         /* Allocate memory. */
 579
 580         if (sp->csegs == NULL) {
 581                 allocate(sp->csegs, XSegment, sp->N);
 582                 allocate(sp->old_segs, XSegment, sp->N * tail_len);
 583                 allocate(sp->nold_segs, int, tail_len);
 584                 allocate(sp->lastx, short, sp->N * 2);
 585                 allocate(sp->x, double, sp->N * 2);
 586                 allocate(sp->diffx, double, sp->N * 2);
 587                 allocate(sp->w, double, sp->Nvortex);
 588                 allocate(sp->olddiffx, double, sp->N * 2);
 589                 allocate(sp->tempdiffx, double, sp->N * 2);
 590                 allocate(sp->tempx, double, sp->N * 2);
 591                 allocate(sp->dead, short, sp->N);
 592                 allocate(sp->boundary, XSegment, n_bound_p);
 593                 allocate(sp->xs, double, sp->Nvortex * 2);
 594                 allocate(sp->x_is_zero, short, sp->Nvortex);
 595                 allocate(sp->p, double, sp->N * 2);
 596                 allocate(sp->mod_dp2, double, sp->N);
 597         }
 598         for (i=0;i<tail_len;i++) {
 599                 sp->nold_segs[i] = 0;
 600         }
 601         sp->c_old_seg = 0;
 602         (void) memset(sp->dead,0,sp->N*sizeof(short));
 603
 604         if (variable_boundary)
 605         {
 606         /* Initialize polynomial p */
 607 /*
 608   The polynomial p(z) = z + c_2 z^2 + ... c_n z^n needs to be
 609   a bijection of the unit disk onto its image.  This is achieved
 610   by insisting that sum_{k=2}^n k |c_k| <= 1.  Actually we set
 611   the inequality to be equality (to get more interesting shapes).
 612 */
 613                 mag = 0;
 614                 for(k=2;k<=deg_p;k++)
 615                 {
 616                         r = positive_rand(1.0/k);
 617                         theta = balance_rand(2*M_PI);
 618                         sp->p_coef[2*(k-2)+0]=r*cos(theta);
 619                         sp->p_coef[2*(k-2)+1]=r*sin(theta);
 620                         mag += k*r;
 621                 }
 622                 if (mag > 0.0001) for(k=2;k<=deg_p;k++)
 623                 {
 624                         sp->p_coef[2*(k-2)+0] /= mag;
 625                         sp->p_coef[2*(k-2)+1] /= mag;
 626                 }
 627
 628 #if DEBUG_POINTED_REGION
 629                 for(k=2;k<=deg_p;k++){
 630                         sp->p_coef[2*(k-2)+0]=0;
 631                         sp->p_coef[2*(k-2)+1]=0;
 632                 }
 633                 sp->p_coef[2*(6-2)+0] = 1.0/6.0;
 634 #endif
 635
 636
 637 /* Here we figure out the best rotation of the domain so that it fills as
 638    much of the screen as possible.  The number of angles we look at is determined
 639    by nr_rotates (we look every 180/nr_rotates degrees).
 640    While we figure out the best angle to rotate, we also figure out the correct scaling factors.
 641 */
 642
 643                 for(k=0;k<nr_rotates;k++) {
 644                         low[k] = 1e5;
 645                         high[k] = -1e5;
 646                 }
 647
 648                 for(k=0;k<n_bound_p;k++) {
 649                         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);
 650                         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);
 651                         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);
 652                         angle1 = nr_rotates/M_PI*atan2(p2-pp2,p1-pp1)-nr_rotates/2;
 653                         angle2 = nr_rotates/M_PI*atan2(pn2-p2,pn1-p1)-nr_rotates/2;
 654                         while (angle1<0) angle1+=nr_rotates*2;
 655                         while (angle2<0) angle2+=nr_rotates*2;
 656                         if (angle1>nr_rotates*1.75 && angle2<nr_rotates*0.25) angle2+=nr_rotates*2;
 657                         if (angle1<nr_rotates*0.25 && angle2>nr_rotates*1.75) angle1+=nr_rotates*2;
 658                         if (angle2<angle1) {
 659                                 tempangle=angle1;
 660                                 angle1=angle2;
 661                                 angle2=tempangle;
 662                         }
 663                         for(i=(int)floor(angle1);i<(int)ceil(angle2);i++) {
 664                                 dist = cos((double)i*M_PI/nr_rotates)*p1 + sin((double)i*M_PI/nr_rotates)*p2;
 665                                 if (i%(nr_rotates*2)<nr_rotates) {
 666                                         if (dist>high[i%nr_rotates]) high[i%nr_rotates] = dist;
 667                                         if (dist<low[i%nr_rotates]) low[i%nr_rotates] = dist;
 668                                 }
 669                                 else {
 670                                         if (-dist>high[i%nr_rotates]) high[i%nr_rotates] = -dist;
 671                                         if (-dist<low[i%nr_rotates]) low[i%nr_rotates] = -dist;
 672                                 }
 673                         }
 674                 }
 675                 bestscale = 0;
 676                 for (i=0;i<nr_rotates;i++) {
 677                         xscale = (sp->width-5.0)/(high[i]-low[i]);
 678                         yscale = (sp->height-5.0)/(high[(i+nr_rotates/2)%nr_rotates]-low[(i+nr_rotates/2)%nr_rotates]);
 679                         scale = (xscale>yscale) ? yscale : xscale;
 680                         if (scale>bestscale) {
 681                                 bestscale = scale;
 682                                 besti = i;
 683                         }
 684                 }
 685 /* Here we do the rotation.  The way we do this is to replace the
 686    polynomial p(z) by a^{-1} p(a z) where a = exp(i best_angle).
 687 */
 688                 p1 = 1;
 689                 p2 = 0;
 690                 for(k=2;k<=deg_p;k++)
 691                 {
 692                         mult(p1,p2,cos((double)besti*M_PI/nr_rotates),sin((double)besti*M_PI/nr_rotates));
 693                         mult(sp->p_coef[2*(k-2)+0],sp->p_coef[2*(k-2)+1],p1,p2);
 694                 }
 695
 696                 sp->scale = bestscale;
 697                 sp->xshift = -(low[besti]+high[besti])/2.0*sp->scale+sp->width/2;
 698                 if (besti<nr_rotates/2)
 699                         sp->yshift = -(low[besti+nr_rotates/2]+high[besti+nr_rotates/2])/2.0*sp->scale+sp->height/2;
 700                 else
 701                         sp->yshift = (low[besti-nr_rotates/2]+high[besti-nr_rotates/2])/2.0*sp->scale+sp->height/2;
 702
 703         sp->xshift += sp->xshift2;
 704         sp->yshift += sp->yshift2;
 705
 706 /* Initialize boundary */
 707
 708                 for(k=0;k<n_bound_p;k++)
 709                 {
 710
 711                         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);
 712                         sp->boundary[k].x1 = (short)(p1*sp->scale+sp->xshift);
 713                         sp->boundary[k].y1 = (short)(p2*sp->scale+sp->yshift);
 714                 }
 715                 for(k=1;k<n_bound_p;k++)
 716                 {
 717                         sp->boundary[k].x2 = sp->boundary[k-1].x1;
 718                         sp->boundary[k].y2 = sp->boundary[k-1].y1;
 719                 }
 720                 sp->boundary[0].x2 = sp->boundary[n_bound_p-1].x1;
 721                 sp->boundary[0].y2 = sp->boundary[n_bound_p-1].y1;
 722         }
 723         else
 724         {
 725                 if (sp->width>sp->height)
 726                         sp->radius = sp->height/2.0-5.0;
 727                 else
 728                         sp->radius = sp->width/2.0-5.0;
 729         }
 730
 731         /* Initialize point positions */
 732
 733         for (i=sp->Nvortex;i<sp->N;i++) {
 734                 do {
 735                         r = sqrt(positive_rand(1.0));
 736                         theta = balance_rand(2*M_PI);
 737                         sp->x[2*i+0]=r*cos(theta);
 738                         sp->x[2*i+1]=r*sin(theta);
 739                 /* This is to make sure the initial distribution of points is uniform */
 740                 } while (variable_boundary &&
 741                           calc_mod_dp2(sp->x[2*i+0],sp->x[2*i+1],sp->p_coef)
 742                           < positive_rand(4));
 743         }
 744
 745         n = NRAND(4)+2;
 746         /* number of vortex points with negative vorticity */
 747         if (n%2) {
 748                 np = NRAND(n+1);
 749         }
 750         else {
 751                 /* if n is even make sure that np==n/2 is twice as likely
 752                    as the other possibilities. */
 753                 np = NRAND(n+2);
 754                 if (np==n+1) np=n/2;
 755         }
 756         for(k=0;k<n;k++)
 757         {
 758                 r = sqrt(positive_rand(0.77));
 759                 theta = balance_rand(2*M_PI);
 760                 x=r*cos(theta);
 761                 y=r*sin(theta);
 762                 r = 0.02+positive_rand(0.1);
 763                 w = (2*(k<np)-1)*2.0/sp->Nvortex;
 764                 for (i=sp->Nvortex*k/n;i<sp->Nvortex*(k+1)/n;i++) {
 765                         theta = balance_rand(2*M_PI);
 766                         sp->x[2*i+0]=x + r*cos(theta);
 767                         sp->x[2*i+1]=y + r*sin(theta);
 768                         sp->w[i]=w;
 769                 }
 770         }
 771 }
 772
 773 ENTRYPOINT void
 774 draw_euler2d (ModeInfo * mi)
 775 {
 776         Display    *display = MI_DISPLAY(mi);
 777         Window      window = MI_WINDOW(mi);
 778         GC          gc = MI_GC(mi);
 779         int         b, col, n_non_vortex_segs;
 780         euler2dstruct *sp;
 781
 782         MI_IS_DRAWN(mi) = True;
 783
 784         if (euler2ds == NULL)
 785                 return;
 786         sp = &euler2ds[MI_SCREEN(mi)];
 787         if (sp->csegs == NULL)
 788                 return;
 789
 790         ode_solve(sp);
 791         if (variable_boundary)
 792                 calc_all_p(sp);
 793
 794         sp->cnsegs = 0;
 795         for(b=sp->Nvortex;b<sp->N;b++) if(!sp->dead[b])
 796         {
 797                 sp->csegs[sp->cnsegs].x1 = sp->lastx[2*b+0];
 798                 sp->csegs[sp->cnsegs].y1 = sp->lastx[2*b+1];
 799                 if (variable_boundary)
 800                 {
 801                         sp->csegs[sp->cnsegs].x2 = (short)(sp->p[2*b+0]*sp->scale+sp->xshift);
 802                         sp->csegs[sp->cnsegs].y2 = (short)(sp->p[2*b+1]*sp->scale+sp->yshift);
 803                 }
 804                 else
 805                 {
 806                         sp->csegs[sp->cnsegs].x2 = (short)(sp->x[2*b+0]*sp->radius+sp->width/2);
 807                         sp->csegs[sp->cnsegs].y2 = (short)(sp->x[2*b+1]*sp->radius+sp->height/2);
 808                 }
 809                 sp->lastx[2*b+0] = sp->csegs[sp->cnsegs].x2;
 810                 sp->lastx[2*b+1] = sp->csegs[sp->cnsegs].y2;
 811                 sp->cnsegs++;
 812         }
 813         n_non_vortex_segs = sp->cnsegs;
 814
 815         if (!sp->hide_vortex) for(b=0;b<sp->Nvortex;b++) if(!sp->dead[b])
 816         {
 817                 sp->csegs[sp->cnsegs].x1 = sp->lastx[2*b+0];
 818                 sp->csegs[sp->cnsegs].y1 = sp->lastx[2*b+1];
 819                 if (variable_boundary)
 820                 {
 821                         sp->csegs[sp->cnsegs].x2 = (short)(sp->p[2*b+0]*sp->scale+sp->xshift);
 822                         sp->csegs[sp->cnsegs].y2 = (short)(sp->p[2*b+1]*sp->scale+sp->yshift);
 823                 }
 824                 else
 825                 {
 826                         sp->csegs[sp->cnsegs].x2 = (short)(sp->x[2*b+0]*sp->radius+sp->width/2);
 827                         sp->csegs[sp->cnsegs].y2 = (short)(sp->x[2*b+1]*sp->radius+sp->height/2);
 828                 }
 829                 sp->lastx[2*b+0] = sp->csegs[sp->cnsegs].x2;
 830                 sp->lastx[2*b+1] = sp->csegs[sp->cnsegs].y2;
 831                 sp->cnsegs++;
 832         }
 833
 834         if (sp->count) {
 835                 XSetForeground(display, gc, MI_BLACK_PIXEL(mi));
 836
 837                 XDrawSegments(display, window, gc, sp->old_segs+sp->c_old_seg*sp->N, sp->nold_segs[sp->c_old_seg]);
 838
 839                 if (MI_NPIXELS(mi) > 2){ /* render colour */
 840                         for (col = 0; col < MI_NPIXELS(mi); col++) {
 841                           int start = col*n_non_vortex_segs/MI_NPIXELS(mi);
 842                           int finish = (col+1)*n_non_vortex_segs/MI_NPIXELS(mi);
 843                           XSetForeground(display, gc, MI_PIXEL(mi, col));
 844                           XDrawSegments(display, window, gc,sp->csegs+start, finish-start);
 845                         }
 846                         if (!sp->hide_vortex) {
 847                           XSetForeground(display, gc, MI_WHITE_PIXEL(mi));
 848                           XDrawSegments(display, window, gc,sp->csegs+n_non_vortex_segs, sp->cnsegs-n_non_vortex_segs);
 849                         }
 850
 851                 } else {                /* render mono */
 852                   XSetForeground(display, gc, MI_WHITE_PIXEL(mi));
 853                   XDrawSegments(display, window, gc,
 854                                                 sp->csegs, sp->cnsegs);
 855                 }
 856
 857                 if (MI_NPIXELS(mi) > 2) /* render colour */
 858                         XSetForeground(display, gc, MI_PIXEL(mi, sp->boundary_color));
 859                 else
 860                         XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_WHITE_PIXEL(mi));
 861                 if (variable_boundary)
 862                         XDrawSegments(display, window, gc,
 863                                                 sp->boundary, n_bound_p);
 864                 else
 865                         XDrawArc(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
 866                                  sp->width/2 - (int) sp->radius - 1, sp->height/2 - (int) sp->radius -1,
 867                                  (int) (2*sp->radius) + 2, (int) (2* sp->radius) + 2, 0, 64*360);
 868
 869                 /* Copy to erase-list */
 870                 (void) memcpy(sp->old_segs+sp->c_old_seg*sp->N, sp->csegs, sp->cnsegs*sizeof(XSegment));
 871                 sp->nold_segs[sp->c_old_seg] = sp->cnsegs;
 872                 sp->c_old_seg++;
 873                 if (sp->c_old_seg >= tail_len)
 874                   sp->c_old_seg = 0;
 875         }
 876
 877         if (++sp->count > MI_CYCLES(mi)) /* pick a new flow */
 878                 init_euler2d(mi);
 879
 880 }
 881
 882 #ifndef STANDALONE
 883 ENTRYPOINT void
 884 refresh_euler2d (ModeInfo * mi)
 885 {
 886         MI_CLEARWINDOW(mi);
 887 }
 888 #endif
 889
 890 XSCREENSAVER_MODULE ("Euler2D", euler2d)
 891
 892 #endif /* MODE_euler2d */