From http://www.jwz.org/xscreensaver/xscreensaver-5.40.tar.gz
[xscreensaver] / hacks / discrete.c
1 /* -*- Mode: C; tab-width: 4 -*- */
2 /* discrete --- chaotic mappings */
3
4 #if 0
5 static const char sccsid[] = "@(#)discrete.c    5.00 2000/11/01 xlockmore";
6 #endif
7
8 /*-
9  * Copyright (c) 1996 by Tim Auckland <tda10.geo@yahoo.com>
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  * "discrete" shows a number of fractals based on the "discrete map"
24  * type of dynamical systems.  They include a different way of looking
25  * at the HOPALONG system, an inverse julia-set iteration, the "Standard
26  * Map" and the "Bird in a Thornbush" fractal.
27  *
28  * Revision History:
29  * 01-Nov-2000: Allocation checks
30  * 31-Jul-1997: Ported to xlockmore-4
31  * 08-Aug-1996: Adapted from hop.c Copyright (c) 1991 by Patrick J. Naughton.
32  */
33
34 #ifdef STANDALONE
35 # define MODE_discrete
36 #define DEFAULTS        "*delay: 20000 \n" \
37                                         "*count:  4096 \n" \
38                                         "*cycles: 2500 \n" \
39                                         "*ncolors: 100 \n" \
40                                         "*fpsSolid: true \n" \
41                                     "*ignoreRotation: True \n" \
42                                     "*lowrez: True \n" \
43
44 # define SMOOTH_COLORS
45 # define release_discrete 0
46 # define discrete_handle_event 0
47 # include "xlockmore.h"         /* in xscreensaver distribution */
48 #else /* STANDALONE */
49 # include "xlock.h"             /* in xlockmore distribution */
50 #endif /* STANDALONE */
51
52 #ifdef MODE_discrete
53
54 ENTRYPOINT ModeSpecOpt discrete_opts =
55 {0, (XrmOptionDescRec *) NULL, 0, (argtype *) NULL, (OptionStruct *) NULL};
56
57 #ifdef USE_MODULES
58 ModStruct   discrete_description =
59 {"discrete", "init_discrete", "draw_discrete", (char *) NULL,
60  "refresh_discrete", "init_discrete", "free_discrete", &discrete_opts,
61  1000, 4096, 2500, 1, 64, 1.0, "",
62  "Shows various discrete maps", 0, NULL};
63
64 #endif
65
66 enum ftypes {
67         SQRT, BIRDIE, STANDARD, TRIG, CUBIC, HENON, AILUJ, HSHOE, DELOG
68 };
69
70 /*#define TEST STANDARD */
71
72 #define BIASES 18
73 static enum ftypes bias[BIASES] =
74 {
75         STANDARD, STANDARD, STANDARD, STANDARD,
76         SQRT, SQRT, SQRT, SQRT,
77         BIRDIE, BIRDIE, BIRDIE,
78         AILUJ, AILUJ, AILUJ,
79         TRIG, TRIG,
80         CUBIC,
81         HENON,
82 };
83
84 typedef struct {
85         int         maxx;
86         int         maxy;       /* max of the screen */
87         double      a;
88         double      b;
89         double      c;
90         double      d;
91         double      e;
92         double      i;
93         double      j;          /* discrete parameters */
94         double      ic;
95         double      jc;
96         double      is;
97         double      js;
98         int         inc;
99         int         pix;
100         enum ftypes op;
101         int         count;
102         XPoint     *pointBuffer;        /* pointer for XDrawPoints */
103
104     int sqrt_sign, std_sign;
105
106 } discretestruct;
107
108 static discretestruct *discretes = (discretestruct *) NULL;
109
110 ENTRYPOINT void
111 init_discrete (ModeInfo * mi)
112 {
113         double      range;
114         discretestruct *hp;
115
116         MI_INIT (mi, discretes);
117         hp = &discretes[MI_SCREEN(mi)];
118
119         hp->maxx = MI_WIDTH(mi);
120         hp->maxy = MI_HEIGHT(mi);
121 #ifdef TEST
122         hp->op = TEST;
123 #else
124         hp->op = bias[LRAND() % BIASES];
125 #endif
126         switch (hp->op) {
127                 case HSHOE:
128                         hp->ic = 0;
129                         hp->jc = 0;
130                         hp->is = hp->maxx / (4);
131                         hp->js = hp->maxy / (4);
132                         hp->a = 0.5;
133                         hp->b = 0.5;
134                         hp->c = 0.2;
135                         hp->d = -1.25;
136                         hp->e = 1;
137                         hp->i = hp->j = 0.0;
138                         break;
139                 case DELOG:
140                         hp->ic = 0.5;
141                         hp->jc = 0.3;
142                         hp->is = hp->maxx / 1.5;
143                         hp->js = hp->maxy / 1.5;
144                         hp->a = 2.176399;
145                         hp->i = hp->j = 0.01;
146                         break;
147                 case HENON:
148                         hp->jc = ((LRAND() / MAXRAND) * 2.0 - 1.0) * 0.4;
149                         hp->ic = 1.3 * (1 - (hp->jc * hp->jc) / (0.4 * 0.4));
150                         hp->is = hp->maxx;
151                         hp->js = hp->maxy * 1.5;
152                         hp->a = 1;
153                         hp->b = 1.4;
154                         hp->c = 0.3;
155                         hp->i = hp->j = 0;
156                         break;
157                 case SQRT:
158                         hp->ic = 0;
159                         hp->jc = 0;
160                         hp->is = 1;
161                         hp->js = 1;
162                         range = sqrt((double) hp->maxx * 2 * hp->maxx * 2 +
163                                      (double) hp->maxy * 2 * hp->maxy * 2) /
164                                 (10.0 + LRAND() % 10);
165
166                         hp->a = (LRAND() / MAXRAND) * range - range / 2.0;
167                         hp->b = (LRAND() / MAXRAND) * range - range / 2.0;
168                         hp->c = (LRAND() / MAXRAND) * range - range / 2.0;
169                         if (!(LRAND() % 2))
170                                 hp->c = 0.0;
171                         hp->i = hp->j = 0.0;
172                         break;
173                 case STANDARD:
174                         hp->ic = M_PI;
175                         hp->jc = M_PI;
176                         hp->is = hp->maxx / (M_PI * 2);
177                         hp->js = hp->maxy / (M_PI * 2);
178                         hp->a = 0;      /* decay */
179                         hp->b = (LRAND() / MAXRAND) * 2.0;
180                         hp->c = 0;
181                         hp->i = M_PI;
182                         hp->j = M_PI;
183                         break;
184                 case BIRDIE:
185                         hp->ic = 0;
186                         hp->jc = 0;
187                         hp->is = hp->maxx / 2;
188                         hp->js = hp->maxy / 2;
189                         hp->a = 1.99 + ((LRAND() / MAXRAND) * 2.0 - 1.0) * 0.2;
190                         hp->b = 0;
191                         hp->c = 0.8 + ((LRAND() / MAXRAND) * 2.0 - 1.0) * 0.1;
192                         hp->i = hp->j = 0;
193                         break;
194                 case TRIG:
195                         hp->a = 5;
196                         hp->b = 0.5 + ((LRAND() / MAXRAND) * 2.0 - 1.0) * 0.3;
197                         hp->ic = hp->a;
198                         hp->jc = 0;
199                         hp->is = hp->maxx / (hp->b * 20);
200                         hp->js = hp->maxy / (hp->b * 20);
201                         hp->i = hp->j = 0;
202                         break;
203                 case CUBIC:
204                         hp->a = 2.77;
205                         hp->b = 0.1 + ((LRAND() / MAXRAND) * 2.0 - 1.0) * 0.1;
206                         hp->ic = 0;
207                         hp->jc = 0;
208                         hp->is = hp->maxx / 4;
209                         hp->js = hp->maxy / 4;
210                         hp->i = hp->j = 0.1;
211                         break;
212                 case AILUJ:
213                         {
214                                 int         i;
215                                 double      x, y, xn, yn;
216
217                                 hp->ic = 0;
218                                 hp->jc = 0;
219                                 hp->is = hp->maxx / 4;
220                                 hp->js = hp->maxx / 4;
221                                 do {
222                                         hp->a = ((LRAND() / MAXRAND) * 2.0 - 1.0) * 1.5 - 0.5;
223                                         hp->b = ((LRAND() / MAXRAND) * 2.0 - 1.0) * 1.5;
224                                         x = y = 0;
225 #define MAXITER 10
226                                         for (i = 0; i < MAXITER && x * x + y * y < 13; i++) {   /* 'Brot calc */
227                                                 xn = x * x - y * y + hp->a;
228                                                 yn = 2 * x * y + hp->b;
229                                                 x = xn;
230                                                 y = yn;
231                                         }
232                                 } while (i < MAXITER);  /* wait for a connected set */
233                                 hp->i = hp->j = 0.1;
234                                 break;
235                         }
236         }
237         hp->pix = 0;
238         hp->inc = 0;
239
240         if (hp->pointBuffer == NULL) {
241                 hp->pointBuffer = (XPoint *) malloc(sizeof (XPoint) * MI_COUNT(mi));
242                 /* if fails will check later */
243         }
244
245         /* Clear the background. */
246         MI_CLEARWINDOW(mi);
247
248         XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_WHITE_PIXEL(mi));
249         hp->count = 0;
250     hp->sqrt_sign = 1;
251     hp->std_sign = 1;
252 }
253
254
255 static void
256 draw_discrete_1 (ModeInfo * mi)
257 {
258         Display    *dsp = MI_DISPLAY(mi);
259         Window      win = MI_WINDOW(mi);
260         double      oldj, oldi;
261         int         count = MI_COUNT(mi);
262         int         cycles = MI_CYCLES(mi);
263         int         k;
264         XPoint     *xp;
265         GC          gc = MI_GC(mi);
266         discretestruct *hp;
267
268         if (discretes == NULL)
269                 return;
270         hp = &discretes[MI_SCREEN(mi)];
271         if (hp->pointBuffer == NULL)
272                 return;
273
274         k = count;
275         xp = hp->pointBuffer;
276
277         hp->inc++;
278
279         MI_IS_DRAWN(mi) = True;
280
281         if (MI_NPIXELS(mi) > 2) {
282                 XSetForeground(dsp, gc, MI_PIXEL(mi, hp->pix));
283                 if (++hp->pix >= MI_NPIXELS(mi))
284                         hp->pix = 0;
285         }
286         while (k--) {
287                 oldj = hp->j;
288                 oldi = hp->i;
289                 switch (hp->op) {
290                         case HSHOE:
291                                 {
292                                         int         i;
293
294 #if 0
295                                         if (!k) {
296                                                 XSetForeground(dsp, gc, MI_BLACK_PIXEL(mi));
297                                                 XFillRectangle(dsp, win, gc, 0, 0, hp->maxx, hp->maxy);
298                                                 XSetForeground(dsp, gc, MI_PIXEL(mi, hp->pix));
299                                         } else
300 #endif
301 #define HD
302 #ifdef HD
303                                         if (k < count / 4) {
304                                                 hp->i = ((double) k / count) * 8 - 1;
305                                                 hp->j = 1;
306                                         } else if (k < count / 2) {
307                                                 hp->i = 1;
308                                                 hp->j = 3 - ((double) k / count) * 8;
309                                         } else if (k < 3 * count / 4) {
310                                                 hp->i = 5 - ((double) k / count) * 8;
311                                                 hp->j = -1;
312                                         } else {
313                                                 hp->i = -1;
314                                                 hp->j = ((double) k / count) * 8 - 7;
315                                         }
316                                         for (i = 1; i < (hp->inc % 15); i++) {
317                                                 oldj = hp->j;
318                                                 oldi = hp->i;
319 #endif
320                                                 hp->i = (hp->a * oldi + hp->b) * oldj;
321                                                 hp->j = (hp->e - hp->d + hp->c * oldi) * oldj * oldj - hp->c * oldi + hp->d;
322 #ifdef HD
323                                         }
324 #endif
325                                         break;
326                                 }
327                         case DELOG:
328                                 hp->j = oldi;
329                                 hp->i = hp->a * oldi * (1 - oldj);
330                                 break;
331                         case HENON:
332                                 hp->i = oldj + hp->a - hp->b * oldi * oldi;
333                                 hp->j = hp->c * oldi;
334                                 break;
335                         case SQRT:
336                                 if (k) {
337                                         hp->j = hp->a + hp->i;
338                                         hp->i = -oldj + (hp->i < 0
339                                         ? sqrt(fabs(hp->b * (hp->i - hp->c)))
340                                                          : -sqrt(fabs(hp->b * (hp->i - hp->c))));
341                                 } else {
342                                         hp->i = (hp->sqrt_sign ? 1 : -1) * hp->inc * hp->maxx / cycles / 2;
343                                         hp->j = hp->a + hp->i;
344                                         hp->sqrt_sign = !hp->sqrt_sign;
345                                 }
346                                 break;
347                         case STANDARD:
348                                 if (k) {
349                                         hp->j = (1 - hp->a) * oldj + hp->b * sin(oldi) + hp->a * hp->c;
350                                         hp->j = fmod(hp->j + 2 * M_PI, 2 * M_PI);
351                                         hp->i = oldi + hp->j;
352                                         hp->i = fmod(hp->i + 2 * M_PI, 2 * M_PI);
353                                 } else {
354                                         hp->j = M_PI + fmod((hp->std_sign ? 1 : -1) * hp->inc * 2 * M_PI / (cycles - 0.5), M_PI);
355                                         hp->i = M_PI;
356                                         hp->std_sign = !hp->std_sign;
357                                 }
358                                 break;
359                         case BIRDIE:
360                                 hp->j = oldi;
361                                 hp->i = (1 - hp->c) * cos(M_PI * hp->a * oldj) + hp->c * hp->b;
362                                 hp->b = oldj;
363                                 break;
364                         case TRIG:
365                                 {
366                                         double      r2 = oldi * oldi + oldj * oldj;
367
368                                         hp->i = hp->a + hp->b * (oldi * cos(r2) - oldj * sin(r2));
369                                         hp->j = hp->b * (oldj * cos(r2) + oldi * sin(r2));
370                                 }
371                                 break;
372                         case CUBIC:
373                                 hp->i = oldj;
374                                 hp->j = hp->a * oldj - oldj * oldj * oldj - hp->b * oldi;
375                                 break;
376                         case AILUJ:
377                                 hp->i = ((LRAND() < MAXRAND / 2) ? -1 : 1) *
378                                         sqrt(((oldi - hp->a) +
379                                               sqrt((oldi - hp->a) * (oldi - hp->a) + (oldj - hp->b) * (oldj - hp->b))) / 2);
380                                 if (hp->i < 0.00000001 && hp->i > -0.00000001)
381                                         hp->i = (hp->i > 0.0) ? 0.00000001 : -0.00000001;
382                                 hp->j = (oldj - hp->b) / (2 * hp->i);
383                                 break;
384                 }
385                 xp->x = hp->maxx / 2 + (int) ((hp->i - hp->ic) * hp->is);
386                 xp->y = hp->maxy / 2 - (int) ((hp->j - hp->jc) * hp->js);
387                 xp++;
388         }
389         XDrawPoints(dsp, win, gc, hp->pointBuffer, count, CoordModeOrigin);
390 }
391
392 ENTRYPOINT void
393 draw_discrete (ModeInfo * mi)
394 {
395   discretestruct *hp = &discretes[MI_SCREEN(mi)];
396   int cycles = MI_CYCLES(mi);
397   int i;
398
399   for (i = 0; i < 10; i++) {
400     draw_discrete_1 (mi);
401     hp->count++;
402   }
403
404   if (hp->count > cycles) {
405     init_discrete(mi);
406   }
407 }
408
409
410 ENTRYPOINT void
411 reshape_discrete(ModeInfo * mi, int width, int height)
412 {
413   discretestruct *hp = &discretes[MI_SCREEN(mi)];
414   hp->maxx = width;
415   hp->maxy = height;
416   XClearWindow (MI_DISPLAY (mi), MI_WINDOW(mi));
417 }
418
419 ENTRYPOINT void
420 free_discrete(ModeInfo * mi)
421 {
422         discretestruct *hp = &discretes[MI_SCREEN(mi)];
423
424         if (hp->pointBuffer != NULL) {
425                 (void) free((void *) hp->pointBuffer);
426                 /* hp->pointBuffer = NULL; */
427         }
428 }
429
430 #ifndef STANDALONE
431 ENTRYPOINT void
432 refresh_discrete(ModeInfo * mi)
433 {
434         MI_CLEARWINDOW(mi);
435 }
436 #endif
437
438 XSCREENSAVER_MODULE ("Discrete", discrete)
439
440 #endif /* MODE_discrete */