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