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