Importing xlockmore 5.22

1 files changed, 879 insertions, 0 deletions

diff --git a/app/xlockmore/modes/euler2d.c b/app/xlockmore/modes/euler2d.c
new file mode 100644
index 000000000..f744387f3
--- /dev/null
+++ b/app/xlockmore/modes/euler2d.c
@@ -0,0 +1,879 @@
+/* -*- Mode: C; tab-width: 4 -*- */
+/* euler2d --- 2 Dimensional Incompressible Inviscid Fluid Flow */
+
+#if !defined( lint ) && !defined( SABER )
+static const char sccsid[] = "@(#)euler2d.c	5.00 2000/11/01 xlockmore";
+
+#endif
+
+/*
+ * Copyright (c) 2000 by Stephen Montgomery-Smith <stephen@math.missouri.edu>
+ *
+ * Permission to use, copy, modify, and distribute this software and its
+ * documentation for any purpose and without fee is hereby granted,
+ * provided that the above copyright notice appear in all copies and that
+ * both that copyright notice and this permission notice appear in
+ * supporting documentation.
+ *
+ * This file is provided AS IS with no warranties of any kind.  The author
+ * shall have no liability with respect to the infringement of copyrights,
+ * trade secrets or any patents by this file or any part thereof.  In no
+ * event will the author be liable for any lost revenue or profits or
+ * other special, indirect and consequential damages.
+ *
+ * Revision History:
+ * 04-Nov-2000: Added an option eulerpower.  This allows for example the
+ *              quasi-geostrophic equation by setting eulerpower to 2.
+ * 01-Nov-2000: Allocation checks.
+ * 10-Sep-2000: Added optimizations, and removed subtle_perturb, by stephen.
+ * 03-Sep-2000: Changed method of solving ode to Adams-Bashforth of order 2.
+ *              Previously used a rather compilcated method of order 4.
+ *              This doubles the speed of the program.  Also it seems
+ *              to have improved numerical stability.  Done by stephen.
+ * 27-Aug-2000: Added rotation of region to maximize screen fill by stephen.
+ * 05-Jun-2000: Adapted from flow.c Copyright (c) 1996 by Tim Auckland
+ * 18-Jul-1996: Adapted from swarm.c Copyright (c) 1991 by Patrick J. Naughton.
+ * 31-Aug-1990: Adapted from xswarm by Jeff Butterworth. (butterwo@ncsc.org)
+ */
+
+/*
+ * The mathematical aspects of this program are discussed in the file
+ * euler2d.tex.
+ */
+
+#ifdef STANDALONE
+#define MODE_euler2d
+#define PROGCLASS "Euler2d"
+#define HACK_INIT init_euler2d
+#define HACK_DRAW draw_euler2d
+#define euler2d_opts xlockmore_opts
+#define DEFAULTS "*delay: 1000 \n" \
+"*count: 1024 \n" \
+"*cycles: 3000 \n" \
+"*ncolors: 200 \n"
+#define SMOOTH_COLORS
+#include "xlockmore.h"		/* in xscreensaver distribution */
+#else /* STANDALONE */
+#include "xlock.h"		/* in xlockmore distribution */
+#endif /* STANDALONE */
+
+#ifdef MODE_euler2d
+
+#define DEF_EULERTAIL "10"
+
+/* #define USE_POINTED_REGION    1 */
+
+static int tail_len;
+static int variable_boundary = 1;
+static float power = 1;
+
+static XrmOptionDescRec opts[] =
+{
+  {(char* ) "-eulertail", (char *) ".euler2d.eulertail",
+   XrmoptionSepArg, (caddr_t) NULL},
+  {(char* ) "-eulerpower", (char *) ".euler2d.eulerpower",
+   XrmoptionSepArg, (caddr_t) NULL},
+};
+static argtype vars[] =
+{
+  {(void *) &tail_len, (char *) "eulertail",
+   (char *) "EulerTail", (char *) DEF_EULERTAIL, t_Int},
+  {(void *) &power, (char *) "eulerpower",
+   (char *) "EulerPower", (char *) "1", t_Float},
+};
+static OptionStruct desc[] =
+{
+  {(char *) "-eulertail len", (char *) "Length of Euler2d tails"},
+  {(char *) "-eulerpower power", (char *) "power of interaction law for points for Euler2d"},
+};
+
+ModeSpecOpt euler2d_opts =
+{sizeof opts / sizeof opts[0], opts,
+ sizeof vars / sizeof vars[0], vars, desc};
+
+#ifdef USE_MODULES
+ModStruct   euler2d_description = {
+	"euler2d", "init_euler2d", "draw_euler2d", "release_euler2d",
+	"refresh_euler2d", "init_euler2d", (char *) NULL, &euler2d_opts,
+	1000, 1024, 3000, 1, 64, 1.0, "",
+	"Simulates 2D incompressible invisid fluid.", 0, NULL
+};
+
+#endif
+
+#define balance_rand(v)	((LRAND()/MAXRAND*(v))-((v)/2))	/* random around 0 */
+#define positive_rand(v)	(LRAND()/MAXRAND*(v))	/* positive random */
+
+#define number_of_vortex_points 20
+
+#define n_bound_p 500
+#define deg_p 6
+
+static double delta_t;
+
+typedef struct {
+	int        width;
+	int        height;
+	int        count;
+	double     xshift,yshift,scale;
+	double     radius;
+
+        int        N;
+        int        Nvortex;
+
+/*  x[2i+0] = x coord for nth point
+    x[2i+1] = y coord for nth point
+    w[i] = vorticity at nth point
+*/
+	double     *x;
+	double     *w;
+
+	double     *diffx;
+	double     *olddiffx;
+	double     *tempx;
+	double     *tempdiffx;
+/*  (xs[2i+0],xs[2i+1]) is reflection of (x[2i+0],x[2i+1]) about unit circle
+    xs[2i+0] = x[2i+0]/nx
+    xs[2i+1] = x[2i+1]/nx
+    where
+    nx = x[2i+0]*x[2i+0] + x[2i+1]*x[2i+1]
+
+    x_is_zero[i] = (nx < 1e-10)
+*/
+	double     *xs;
+	short      *x_is_zero;
+
+/*  (p[2i+0],p[2i+1]) is image of (x[2i+0],x[2i+1]) under polynomial p.
+    mod_dp2 is |p'(z)|^2 when z = (x[2i+0],x[2i+1]).
+*/
+	double	   *p;
+	double	   *mod_dp2;
+
+/* Sometimes in our calculations we get overflow or numbers that are too big.
+   If that happens with the point x[2*i+0], x[2*i+1], we set dead[i].
+*/
+	short      *dead;
+
+	XSegment   *csegs;
+	int         cnsegs;
+	XSegment   *old_segs;
+	int        *nold_segs;
+	int         c_old_seg;
+	int	    boundary_color;
+	int         hide_vortex;
+	short      *lastx;
+
+	double      p_coef[2*(deg_p-1)];
+	XSegment   *boundary;
+
+} euler2dstruct;
+
+static euler2dstruct *euler2ds = (euler2dstruct *) NULL;
+
+/*
+  If variable_boundary == 1, then we make a variable boundary.
+  The way this is done is to map the unit disk under a
+  polynomial p, where
+  p(z) = z + c_2 z^2 + ... + c_n z^n
+  where n = deg_p.  sp->p_coef contains the complex numbers
+  c_2, c_3, ... c_n.
+*/
+
+#define add(a1,a2,b1,b2) (a1)+=(b1);(a2)+=(b2)
+#define mult(a1,a2,b1,b2) temp=(a1)*(b1)-(a2)*(b2); \
+                          (a2)=(a1)*(b2)+(a2)*(b1);(a1)=temp
+
+static void
+calc_p(double *p1, double *p2, double z1, double z2, double p_coef[])
+{
+  int i;
+  double temp;
+
+  *p1=0;
+  *p2=0;
+  for(i=deg_p;i>=2;i--)
+  {
+    add(*p1,*p2,p_coef[(i-2)*2],p_coef[(i-2)*2+1]);
+    mult(*p1,*p2,z1,z2);
+  }
+  add(*p1,*p2,1,0);
+  mult(*p1,*p2,z1,z2);
+}
+
+/* Calculate |p'(z)|^2 */
+static double
+calc_mod_dp2(double z1, double z2, double p_coef[])
+{
+  int i;
+  double temp,mp1,mp2;
+
+  mp1=0;
+  mp2=0;
+  for(i=deg_p;i>=2;i--)
+  {
+    add(mp1,mp2,i*p_coef[(i-2)*2],i*p_coef[(i-2)*2+1]);
+    mult(mp1,mp2,z1,z2);
+  }
+  add(mp1,mp2,1,0);
+  return mp1*mp1+mp2*mp2;
+}
+
+static void
+calc_all_p(euler2dstruct *sp)
+{
+  int i,j;
+  double temp,p1,p2,z1,z2;
+  for(j=(sp->hide_vortex?sp->Nvortex:0);j<sp->N;j++) if(!sp->dead[j])
+  {
+    p1=0;
+    p2=0;
+    z1=sp->x[2*j+0];
+    z2=sp->x[2*j+1];
+    for(i=deg_p;i>=2;i--)
+    {
+      add(p1,p2,sp->p_coef[(i-2)*2],sp->p_coef[(i-2)*2+1]);
+      mult(p1,p2,z1,z2);
+    }
+    add(p1,p2,1,0);
+    mult(p1,p2,z1,z2);
+    sp->p[2*j+0] = p1;
+    sp->p[2*j+1] = p2;
+  }
+}
+
+static void
+calc_all_mod_dp2(double *x, euler2dstruct *sp)
+{
+  int i,j;
+  double temp,mp1,mp2,z1,z2;
+  for(j=0;j<sp->N;j++) if(!sp->dead[j])
+  {
+    mp1=0;
+    mp2=0;
+    z1=x[2*j+0];
+    z2=x[2*j+1];
+    for(i=deg_p;i>=2;i--)
+    {
+      add(mp1,mp2,i*sp->p_coef[(i-2)*2],i*sp->p_coef[(i-2)*2+1]);
+      mult(mp1,mp2,z1,z2);
+    }
+    add(mp1,mp2,1,0);
+    sp->mod_dp2[j] = mp1*mp1+mp2*mp2;
+  }
+}
+
+static void
+derivs(double *x, euler2dstruct *sp)
+{
+  int i,j;
+  double u1,u2,x1,x2,xij1,xij2,nxij;
+  double nx;
+
+  if (variable_boundary)
+    calc_all_mod_dp2(sp->x,sp);
+
+  for (j=0;j<sp->Nvortex;j++) if (!sp->dead[j])
+  {
+    nx = x[2*j+0]*x[2*j+0] + x[2*j+1]*x[2*j+1];
+    if (nx < 1e-10)
+      sp->x_is_zero[j] = 1;
+    else {
+      sp->x_is_zero[j] = 0;
+      sp->xs[2*j+0] = x[2*j+0]/nx;
+      sp->xs[2*j+1] = x[2*j+1]/nx;
+    }
+  }
+
+  (void) memset(sp->diffx,0,sizeof(double)*2*sp->N);
+
+  for (i=0;i<sp->N;i++) if (!sp->dead[i])
+  {
+    x1 = x[2*i+0];
+    x2 = x[2*i+1];
+    for (j=0;j<sp->Nvortex;j++) if (!sp->dead[j])
+    {
+/*
+  Calculate the Biot-Savart kernel, that is, effect of a
+  vortex point at a = (x[2*j+0],x[2*j+1]) at the point
+  x = (x1,x2), returning the vector field in (u1,u2).
+
+  In the plane, this is given by the formula
+
+  u = (x-a)/|x-a|^2  or zero if x=a.
+
+  However, in the unit disk we have to subtract from the
+  above:
+
+  (x-as)/|x-as|^2
+
+  where as = a/|a|^2 is the reflection of a about the unit circle.
+
+  If however power != 1, then
+
+  u = (x-a)/|x-a|^(power+1)  -  |a|^(1-power) (x-as)/|x-as|^(power+1)
+
+*/
+
+      xij1 = x1 - x[2*j+0];
+      xij2 = x2 - x[2*j+1];
+      nxij = (power==1.0) ? xij1*xij1+xij2*xij2 : pow(xij1*xij1+xij2*xij2,(power+1)/2.0);
+
+      if(nxij >= 1e-4)  {
+        u1 =  xij2/nxij;
+        u2 = -xij1/nxij;
+      }
+      else
+        u1 = u2 = 0.0;
+
+      if (!sp->x_is_zero[j])
+      {
+        xij1 = x1 - sp->xs[2*j+0];
+        xij2 = x2 - sp->xs[2*j+1];
+        nxij = (power==1.0) ? xij1*xij1+xij2*xij2 : pow(xij1*xij1+xij2*xij2,(power+1)/2.0);
+
+        if (nxij < 1e-5)
+        {
+          sp->dead[i] = 1;
+          u1 = u2 = 0.0;
+        }
+        else
+        {
+          u1 -= xij2/nxij;
+          u2 += xij1/nxij;
+        }
+      }
+
+      if (!sp->dead[i])
+      {
+        sp->diffx[2*i+0] += u1*sp->w[j];
+        sp->diffx[2*i+1] += u2*sp->w[j];
+      }
+    }
+
+    if (!sp->dead[i] && variable_boundary)
+    {
+      if (sp->mod_dp2[i] < 1e-5)
+        sp->dead[i] = 1;
+      else
+      {
+        sp->diffx[2*i+0] /= sp->mod_dp2[i];
+        sp->diffx[2*i+1] /= sp->mod_dp2[i];
+      }
+    }
+  }
+}
+
+/*
+  What perturb does is effectively
+    ret = x + k,
+  where k should be of order delta_t.
+
+  We have the option to do this more subtly by mapping points x
+  in the unit disk to points y in the plane, where y = f(|x|) x,
+  with f(t) = -log(1-t)/t.
+
+  This might reduce (but does not remove) problems where particles near
+  the edge of the boundary bounce around.
+
+  But it seems to be not that effective, so for now switch it off.
+*/
+
+#define SUBTLE_PERTURB 0
+
+static void
+perturb(double ret[], double x[], double k[], euler2dstruct *sp)
+{
+  int i;
+  double x1,x2,k1,k2;
+
+#if SUBTLE_PERTURB
+  double d1,d2,t1,t2,mag,mag2,mlog1mmag,memmagdmag,xdotk;
+  for (i=0;i<sp->N;i++) if (!sp->dead[i])
+  {
+    x1 = x[2*i+0];
+    x2 = x[2*i+1];
+    k1 = k[2*i+0];
+    k2 = k[2*i+1];
+    mag2 = x1*x1 + x2*x2;
+    if (mag2 < 1e-10)
+    {
+      ret[2*i+0] = x1+k1;
+      ret[2*i+1] = x2+k2;
+    }
+    else if (mag2 > 1-1e-5)
+      sp->dead[i] = 1;
+    else
+    {
+      mag = sqrt(mag2);
+      mlog1mmag = -log(1-mag);
+      xdotk = x1*k1 + x2*k2;
+      t1 = (x1 + k1)*mlog1mmag/mag + x1*xdotk*(1.0/(1-mag)-mlog1mmag/mag)/mag/mag;
+      t2 = (x2 + k2)*mlog1mmag/mag + x2*xdotk*(1.0/(1-mag)-mlog1mmag/mag)/mag/mag;
+      mag = sqrt(t1*t1+t2*t2);
+      if (mag > 11.5 /* log(1e5) */)
+        sp->dead[i] = 1;
+      else
+      {
+        memmagdmag = (mag>1e-5) ? ((1.0-exp(-mag))/mag) : (1-mag/2.0);
+        ret[2*i+0] = t1*memmagdmag;
+        ret[2*i+1] = t2*memmagdmag;
+      }
+    }
+    if (!sp->dead[i])
+    {
+      d1 = ret[2*i+0]-x1;
+      d2 = ret[2*i+1]-x2;
+      if (d1*d1+d2*d2 > 0.1)
+        sp->dead[i] = 1;
+    }
+  }
+
+#else
+
+  for (i=0;i<sp->N;i++) if (!sp->dead[i])
+  {
+    x1 = x[2*i+0];
+    x2 = x[2*i+1];
+    k1 = k[2*i+0];
+    k2 = k[2*i+1];
+    if (k1*k1+k2*k2 > 0.1 || x1*x1+x2*x2 > 1-1e-5)
+      sp->dead[i] = 1;
+    else
+    {
+      ret[2*i+0] = x1+k1;
+      ret[2*i+1] = x2+k2;
+    }
+  }
+#endif
+}
+
+static void
+ode_solve(euler2dstruct *sp)
+{
+  int i;
+  double *temp;
+
+  if (sp->count < 1) {
+    /* midpoint method */
+    derivs(sp->x,sp);
+    (void) memcpy(sp->olddiffx,sp->diffx,sizeof(double)*2*sp->N);
+    for (i=0;i<sp->N;i++) if (!sp->dead[i]) {
+      sp->tempdiffx[2*i+0] = 0.5*delta_t*sp->diffx[2*i+0];
+      sp->tempdiffx[2*i+1] = 0.5*delta_t*sp->diffx[2*i+1];
+    }
+    perturb(sp->tempx,sp->x,sp->tempdiffx,sp);
+    derivs(sp->tempx,sp);
+    for (i=0;i<sp->N;i++) if (!sp->dead[i]) {
+      sp->tempdiffx[2*i+0] = delta_t*sp->diffx[2*i+0];
+      sp->tempdiffx[2*i+1] = delta_t*sp->diffx[2*i+1];
+    }
+    perturb(sp->x,sp->x,sp->tempdiffx,sp);
+  } else {
+    /* Adams Basforth */
+    derivs(sp->x,sp);
+    for (i=0;i<sp->N;i++) if (!sp->dead[i]) {
+      sp->tempdiffx[2*i+0] = delta_t*(1.5*sp->diffx[2*i+0] - 0.5*sp->olddiffx[2*i+0]);
+      sp->tempdiffx[2*i+1] = delta_t*(1.5*sp->diffx[2*i+1] - 0.5*sp->olddiffx[2*i+1]);
+    }
+    perturb(sp->x,sp->x,sp->tempdiffx,sp);
+    temp = sp->olddiffx;
+    sp->olddiffx = sp->diffx;
+    sp->diffx = temp;
+  }
+}
+
+#define deallocate(p,t) if (p!=NULL) {free(p); p=(t*)NULL; }
+#define allocate(p,t,s) if ((p=(t*)malloc(sizeof(t)*s))==NULL)\
+{free_euler2d(sp);return;}
+
+static void
+free_euler2d(euler2dstruct *sp)
+{
+	deallocate(sp->csegs, XSegment);
+	deallocate(sp->old_segs, XSegment);
+	deallocate(sp->nold_segs, int);
+	deallocate(sp->lastx, short);
+	deallocate(sp->x, double);
+	deallocate(sp->diffx, double);
+	deallocate(sp->w, double);
+	deallocate(sp->olddiffx, double);
+	deallocate(sp->tempdiffx, double);
+	deallocate(sp->tempx, double);
+	deallocate(sp->dead, short);
+	deallocate(sp->boundary, XSegment);
+	deallocate(sp->xs, double);
+	deallocate(sp->x_is_zero, short);
+	deallocate(sp->p, double);
+	deallocate(sp->mod_dp2, double);
+}
+
+void
+init_euler2d(ModeInfo * mi)
+{
+#define nr_rotates 18 /* how many rotations to try to fill as much of screen as possible - must be even number */
+	euler2dstruct *sp;
+	int         i,k,n,np;
+	double      r,theta,x,y,w;
+	double      mag,xscale,yscale,p1,p2;
+	double      low[nr_rotates],high[nr_rotates],pp1,pp2,pn1,pn2,angle1,angle2,tempangle,dist,scale,bestscale,temp;
+	int	    besti = 0;
+
+        if (power<0.5) power = 0.5;
+        if (power>3.0) power = 3.0;
+	variable_boundary &= power == 1.0;
+	delta_t = 0.001;
+        if (power>1.0) delta_t *= pow(0.1,(double) power-1);
+
+	if (euler2ds == NULL) {
+		if ((euler2ds = (euler2dstruct *) calloc(MI_NUM_SCREENS(mi),
+					       sizeof (euler2dstruct))) == NULL)
+			return;
+	}
+	sp = &euler2ds[MI_SCREEN(mi)];
+
+	sp->boundary_color = NRAND(MI_NPIXELS(mi));
+	sp->hide_vortex = NRAND(4) != 0;
+
+	sp->count = 0;
+
+	sp->width = MI_WIDTH(mi);
+	sp->height = MI_HEIGHT(mi);
+
+	sp->N = MI_COUNT(mi)+number_of_vortex_points;
+	sp->Nvortex = number_of_vortex_points;
+
+	if (tail_len < 1) { /* minimum tail */
+	  tail_len = 1;
+	}
+	if (tail_len > MI_CYCLES(mi)) { /* maximum tail */
+	  tail_len = MI_CYCLES(mi);
+	}
+
+	/* Clear the background. */
+	MI_CLEARWINDOW(mi);
+
+	free_euler2d(sp);
+
+	/* Allocate memory. */
+
+	if (sp->csegs == NULL) {
+		allocate(sp->csegs, XSegment, sp->N);
+		allocate(sp->old_segs, XSegment, sp->N * tail_len);
+		allocate(sp->nold_segs, int, tail_len);
+		allocate(sp->lastx, short, sp->N * 2);
+		allocate(sp->x, double, sp->N * 2);
+		allocate(sp->diffx, double, sp->N * 2);
+		allocate(sp->w, double, sp->Nvortex);
+		allocate(sp->olddiffx, double, sp->N * 2);
+		allocate(sp->tempdiffx, double, sp->N * 2);
+		allocate(sp->tempx, double, sp->N * 2);
+		allocate(sp->dead, short, sp->N);
+		allocate(sp->boundary, XSegment, n_bound_p);
+		allocate(sp->xs, double, sp->Nvortex * 2);
+		allocate(sp->x_is_zero, short, sp->Nvortex);
+		allocate(sp->p, double, sp->N * 2);
+		allocate(sp->mod_dp2, double, sp->N);
+	}
+	for (i=0;i<tail_len;i++) {
+		sp->nold_segs[i] = 0;
+	}
+	sp->c_old_seg = 0;
+	(void) memset(sp->dead,0,sp->N*sizeof(short));
+
+	if (variable_boundary)
+	{
+	/* Initialize polynomial p */
+/*
+  The polynomial p(z) = z + c_2 z^2 + ... c_n z^n needs to be
+  a bijection of the unit disk onto its image.  This is achieved
+  by insisting that sum_{k=2}^n k |c_k| <= 1.  Actually we set
+  the inequality to be equality (to get more interesting shapes).
+*/
+		mag = 0;
+		for(k=2;k<=deg_p;k++)
+		{
+			r = positive_rand(1.0/k);
+			theta = balance_rand(2*M_PI);
+			sp->p_coef[2*(k-2)+0]=r*cos(theta);
+			sp->p_coef[2*(k-2)+1]=r*sin(theta);
+			mag += k*r;
+		}
+		if (mag > 0.0001) for(k=2;k<=deg_p;k++)
+		{
+			sp->p_coef[2*(k-2)+0] /= mag;
+			sp->p_coef[2*(k-2)+1] /= mag;
+		}
+
+#ifdef USE_POINTED_REGION
+		/* Five symmetric buldges */
+		for(k=2;k<=deg_p;k++){
+			sp->p_coef[2*(k-2)+0]=0;
+			sp->p_coef[2*(k-2)+1]=0;
+		}
+		sp->p_coef[2*(6-2)+0] = 1.0/6.0;
+#endif
+
+
+/* Here we figure out the best rotation of the domain so that it fills as
+   much of the screen as possible.  The number of angles we look at is determined
+   by nr_rotates (we look every 180/nr_rotates degrees).
+   While we figure out the best angle to rotate, we also figure out the correct scaling factors.
+*/
+
+		for(k=0;k<nr_rotates;k++) {
+			low[k] = 1e5;
+			high[k] = -1e5;
+		}
+
+		for(k=0;k<n_bound_p;k++) {
+			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);
+			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);
+			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);
+			angle1 = nr_rotates/M_PI*atan2(p2-pp2,p1-pp1)-nr_rotates/2;
+			angle2 = nr_rotates/M_PI*atan2(pn2-p2,pn1-p1)-nr_rotates/2;
+			while (angle1<0) angle1+=nr_rotates*2;
+			while (angle2<0) angle2+=nr_rotates*2;
+			if (angle1>nr_rotates*1.75 && angle2<nr_rotates*0.25) angle2+=nr_rotates*2;
+			if (angle1<nr_rotates*0.25 && angle2>nr_rotates*1.75) angle1+=nr_rotates*2;
+                        if (angle2<angle1) {
+				tempangle=angle1;
+				angle1=angle2;
+				angle2=tempangle;
+			}
+			for(i=(int)floor(angle1);i<(int)ceil(angle2);i++) {
+				dist = cos((double)i*M_PI/nr_rotates)*p1 + sin((double)i*M_PI/nr_rotates)*p2;
+				if (i%(nr_rotates*2)<nr_rotates) {
+					if (dist>high[i%nr_rotates]) high[i%nr_rotates] = dist;
+					if (dist<low[i%nr_rotates]) low[i%nr_rotates] = dist;
+				}
+				else {
+					if (-dist>high[i%nr_rotates]) high[i%nr_rotates] = -dist;
+					if (-dist<low[i%nr_rotates]) low[i%nr_rotates] = -dist;
+				}
+			}
+		}
+		bestscale = 0;
+		for (i=0;i<nr_rotates;i++) {
+			xscale = (sp->width-5.0)/(high[i]-low[i]);
+			yscale = (sp->height-5.0)/(high[(i+nr_rotates/2)%nr_rotates]-low[(i+nr_rotates/2)%nr_rotates]);
+			scale = (xscale>yscale) ? yscale : xscale;
+			if (scale>bestscale) {
+				bestscale = scale;
+				besti = i;
+			}
+		}
+/* Here we do the rotation.  The way we do this is to replace the
+   polynomial p(z) by a^{-1} p(a z) where a = exp(i best_angle).
+*/
+		p1 = 1;
+		p2 = 0;
+		for(k=2;k<=deg_p;k++)
+		{
+			mult(p1,p2,cos((double)besti*M_PI/nr_rotates),sin((double)besti*M_PI/nr_rotates));
+			mult(sp->p_coef[2*(k-2)+0],sp->p_coef[2*(k-2)+1],p1,p2);
+		}
+
+		sp->scale = bestscale;
+		sp->xshift = -(low[besti]+high[besti])/2.0*sp->scale+sp->width/2;
+                if (besti<nr_rotates/2)
+			sp->yshift = -(low[besti+nr_rotates/2]+high[besti+nr_rotates/2])/2.0*sp->scale+sp->height/2;
+		else
+			sp->yshift = (low[besti-nr_rotates/2]+high[besti-nr_rotates/2])/2.0*sp->scale+sp->height/2;
+
+
+/* Initialize boundary */
+
+		for(k=0;k<n_bound_p;k++)
+		{
+
+			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);
+			sp->boundary[k].x1 = (short)(p1*sp->scale+sp->xshift);
+			sp->boundary[k].y1 = (short)(p2*sp->scale+sp->yshift);
+		}
+		for(k=1;k<n_bound_p;k++)
+		{
+			sp->boundary[k].x2 = sp->boundary[k-1].x1;
+			sp->boundary[k].y2 = sp->boundary[k-1].y1;
+		}
+		sp->boundary[0].x2 = sp->boundary[n_bound_p-1].x1;
+		sp->boundary[0].y2 = sp->boundary[n_bound_p-1].y1;
+	}
+	else
+	{
+		if (sp->width>sp->height)
+			sp->radius = sp->height/2.0-5.0;
+		else
+			sp->radius = sp->width/2.0-5.0;
+	}
+
+	/* Initialize point positions */
+
+	for (i=sp->Nvortex;i<sp->N;i++) {
+		do {
+			r = sqrt(positive_rand(1.0));
+			theta = balance_rand(2*M_PI);
+			sp->x[2*i+0]=r*cos(theta);
+			sp->x[2*i+1]=r*sin(theta);
+		/* This is to make sure the initial distribution of points is uniform */
+		} while (variable_boundary &&
+                          calc_mod_dp2(sp->x[2*i+0],sp->x[2*i+1],sp->p_coef)
+                          < positive_rand(4));
+	}
+
+	n = NRAND(4)+2;
+	/* number of vortex points with negative vorticity */
+	if (n%2) {
+		np = NRAND(n+1);
+	}
+	else {
+		/* if n is even make sure that np==n/2 is twice as likely
+		   as the other possibilities. */
+		np = NRAND(n+2);
+		if (np==n+1) np=n/2;
+	}
+	for(k=0;k<n;k++)
+	{
+		r = sqrt(positive_rand(0.77));
+		theta = balance_rand(2*M_PI);
+		x=r*cos(theta);
+		y=r*sin(theta);
+		r = 0.02+positive_rand(0.1);
+		w = (2*(k<np)-1)*2.0/sp->Nvortex;
+		for (i=sp->Nvortex*k/n;i<sp->Nvortex*(k+1)/n;i++) {
+			theta = balance_rand(2*M_PI);
+			sp->x[2*i+0]=x + r*cos(theta);
+			sp->x[2*i+1]=y + r*sin(theta);
+			sp->w[i]=w;
+		}
+	}
+}
+
+void
+draw_euler2d(ModeInfo * mi)
+{
+	Display    *display = MI_DISPLAY(mi);
+	Window      window = MI_WINDOW(mi);
+	GC          gc = MI_GC(mi);
+	int         b, col, n_non_vortex_segs;
+	euler2dstruct *sp;
+
+	MI_IS_DRAWN(mi) = True;
+
+	if (euler2ds == NULL)
+		return;
+	sp = &euler2ds[MI_SCREEN(mi)];
+	if (sp->csegs == NULL)
+		return;
+
+	ode_solve(sp);
+	if (variable_boundary)
+		calc_all_p(sp);
+
+	sp->cnsegs = 0;
+	for(b=sp->Nvortex;b<sp->N;b++) if(!sp->dead[b])
+	{
+		sp->csegs[sp->cnsegs].x1 = sp->lastx[2*b+0];
+		sp->csegs[sp->cnsegs].y1 = sp->lastx[2*b+1];
+		if (variable_boundary)
+		{
+			sp->csegs[sp->cnsegs].x2 = (short)(sp->p[2*b+0]*sp->scale+sp->xshift);
+			sp->csegs[sp->cnsegs].y2 = (short)(sp->p[2*b+1]*sp->scale+sp->yshift);
+		}
+		else
+		{
+			sp->csegs[sp->cnsegs].x2 = (short)(sp->x[2*b+0]*sp->radius+sp->width/2);
+			sp->csegs[sp->cnsegs].y2 = (short)(sp->x[2*b+1]*sp->radius+sp->height/2);
+		}
+		sp->lastx[2*b+0] = sp->csegs[sp->cnsegs].x2;
+		sp->lastx[2*b+1] = sp->csegs[sp->cnsegs].y2;
+		sp->cnsegs++;
+	}
+	n_non_vortex_segs = sp->cnsegs;
+
+	if (!sp->hide_vortex) for(b=0;b<sp->Nvortex;b++) if(!sp->dead[b])
+	{
+		sp->csegs[sp->cnsegs].x1 = sp->lastx[2*b+0];
+		sp->csegs[sp->cnsegs].y1 = sp->lastx[2*b+1];
+		if (variable_boundary)
+		{
+			sp->csegs[sp->cnsegs].x2 = (short)(sp->p[2*b+0]*sp->scale+sp->xshift);
+			sp->csegs[sp->cnsegs].y2 = (short)(sp->p[2*b+1]*sp->scale+sp->yshift);
+		}
+		else
+		{
+			sp->csegs[sp->cnsegs].x2 = (short)(sp->x[2*b+0]*sp->radius+sp->width/2);
+			sp->csegs[sp->cnsegs].y2 = (short)(sp->x[2*b+1]*sp->radius+sp->height/2);
+		}
+		sp->lastx[2*b+0] = sp->csegs[sp->cnsegs].x2;
+		sp->lastx[2*b+1] = sp->csegs[sp->cnsegs].y2;
+		sp->cnsegs++;
+	}
+
+	if (sp->count) {
+		XSetForeground(display, gc, MI_BLACK_PIXEL(mi));
+
+		XDrawSegments(display, window, gc, sp->old_segs+sp->c_old_seg*sp->N, sp->nold_segs[sp->c_old_seg]);
+
+		if (MI_NPIXELS(mi) > 2){ /* render colour */
+			for (col = 0; col < MI_NPIXELS(mi); col++) {
+			  int start = col*n_non_vortex_segs/MI_NPIXELS(mi);
+			  int end = (col+1)*n_non_vortex_segs/MI_NPIXELS(mi);
+			  XSetForeground(display, gc, MI_PIXEL(mi, col));
+			  XDrawSegments(display, window, gc,sp->csegs+start, end-start);
+			}
+			if (!sp->hide_vortex) {
+			  XSetForeground(display, gc, MI_WHITE_PIXEL(mi));
+			  XDrawSegments(display, window, gc,sp->csegs+n_non_vortex_segs, sp->cnsegs-n_non_vortex_segs);
+			}
+
+		} else { 		/* render mono */
+		  XSetForeground(display, gc, MI_WHITE_PIXEL(mi));
+		  XDrawSegments(display, window, gc,
+						sp->csegs, sp->cnsegs);
+		}
+
+		if (MI_NPIXELS(mi) > 2) /* render colour */
+			XSetForeground(display, gc, MI_PIXEL(mi, sp->boundary_color));
+		else
+			XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_WHITE_PIXEL(mi));
+		if (variable_boundary)
+			XDrawSegments(display, window, gc,
+						sp->boundary, n_bound_p);
+		else
+			XDrawArc(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
+				 sp->width/2 - (int) sp->radius - 1, sp->height/2 - (int) sp->radius -1,
+				 (int) (2*sp->radius) + 2, (int) (2* sp->radius) + 2, 0, 64*360);
+
+		/* Copy to erase-list */
+		(void) memcpy(sp->old_segs+sp->c_old_seg*sp->N, sp->csegs, sp->cnsegs*sizeof(XSegment));
+		sp->nold_segs[sp->c_old_seg] = sp->cnsegs;
+		sp->c_old_seg++;
+		if (sp->c_old_seg >= tail_len)
+		  sp->c_old_seg = 0;
+	}
+
+	if (++sp->count > MI_CYCLES(mi)) /* pick a new flow */
+		init_euler2d(mi);
+
+}
+
+void
+release_euler2d(ModeInfo * mi)
+{
+	if (euler2ds != NULL) {
+		int         screen;
+
+		for (screen = 0; screen < MI_NUM_SCREENS(mi); screen++)
+			free_euler2d(&euler2ds[screen]);
+		free(euler2ds);
+		euler2ds = (euler2dstruct *) NULL;
+	}
+}
+
+void
+refresh_euler2d(ModeInfo * mi)
+{
+	MI_CLEARWINDOW(mi);
+}
+
+#endif /* MODE_euler2d */