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