--- /dev/null
+/* -*- 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: 10000 \n" \
+"*count: 1024 \n" \
+"*cycles: 3000 \n" \
+"*ncolors: 64 \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 DEBUG_POINTED_REGION 0
+
+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[] =
+{
+ {(caddr_t *) &tail_len, (char *) "eulertail",
+ (char *) "EulerTail", (char *) DEF_EULERTAIL, t_Int},
+ {(caddr_t *) &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) {(void) free((void *) 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,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;
+ }
+
+#if DEBUG_POINTED_REGION
+ 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 finish = (col+1)*n_non_vortex_segs/MI_NPIXELS(mi);
+ XSetForeground(display, gc, MI_PIXEL(mi, col));
+ XDrawSegments(display, window, gc,sp->csegs+start, finish-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]);
+ (void) free((void *) euler2ds);
+ euler2ds = (euler2dstruct *) NULL;
+ }
+}
+
+void
+refresh_euler2d(ModeInfo * mi)
+{
+ MI_CLEARWINDOW(mi);
+}
+
+#endif /* MODE_euler2d */
projects / xscreensaver / blobdiff