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