1 /*
2  * Copyright 2013, Blender Foundation.
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  */
19 /*
20  * Lyapunov Shader - in memory of great mathematician Aleksandr Mikhailovich Lyapunov
21  * Original code: Sylvio Sell - maitag.de - Rostock Germany 2013
22  * OSL port by Thomas Dinges
23  * More information: https://developer.blender.org/T32305
24  */
26 /* Fac_Type
27  * 0, SPREAD,   Spread indices for fac output
28  * 1, ABS,              Absolute values from indices
29  * 2, COLOR,    Get fac output from used colors
30  * 3, REAL,             Real indices
31  */
33 /* Render_Type
34  * 0, NEG,      Negative Lyapunov indices only
35  * 1, POS,      Positive Lyapunov indices only
36  * 2, ALL,      Positive and negative indices
37  */
39 float lyapunov(point p, float iteration_pre, float iteration_main, float p1, float p2)
40 {
41         /* Coordinates */
42         float a = p;
43         float b = p;
44         float c = p;
46         int iter_pre =  (int)floor(iteration_pre);
47         int iter_main = (int)floor(iteration_main);
48         float nabla_pre = iteration_pre - (float)iter_pre;
49         float nabla_main = iteration_main - (float)iter_main;
51         float x = 0.0;
52         float index = 0.0;
53         float derivation = 0.0;
54         int iter = 0;
56         /* Pre-iteration */
57         for (int i = 0; i < iter_pre; i++) {
58                 x = p1 * sin(x + a) * sin(x + a) + p2;
59                 x = p1 * sin(x + b) * sin(x + b) + p2;
60                 x = p1 * sin(x + c) * sin(x + c) + p2;
61         }
63         if (nabla_pre != 0.0) {
64                 float x_pre = x;
66                 x = p1 * sin(x + a) * sin(x + a) + p2;
67                 x = p1 * sin(x + b) * sin(x + b) + p2;
68                 x = p1 * sin(x + c) * sin(x + c) + p2;
69                 x = x * nabla_pre + x_pre * (1.0 - nabla_pre);
70         }
72         /* Main-iteration */
73         for (int i = 0; i < iter_main; i++) {
74                 x = p1 * sin(x + a) * sin(x + a) + p2;
75                 derivation = 2.0 *p1 *sin(x + a) * cos(x + a);
76                 if (derivation != 0.0) { index += log(fabs(derivation)); iter++; }
78                 x = p1 * sin(x + b) * sin(x + b) + p2;
79                 derivation = 2.0 *p1 *sin(x + b) * cos(x + b);
80                 if (derivation != 0.0) { index += log(fabs(derivation)); iter++; }
82                 x = p1 * sin(x + c) * sin(x + c) + p2;
83                 derivation = 2.0 *p1 *sin(x + c) * cos(x + c);
84                 if (derivation != 0.0) { index += log(fabs(derivation)); iter++; }
85         }
87         if (nabla_main == 0.0) {
88                 index = (iter != 0) ? index / (float)(iter) : 0.0;
89         }
90         else {
91                 float index_pre = (iter != 0) ? index / (float)(iter) : 0.0;
93                 x = p1 * sin(x + a) * sin(x + a) + p2;
94                 derivation = 2.0 *p1 *sin(x + a) * cos(x + a);
95                 if (derivation != 0.0) { index += log(fabs(derivation)); iter++; }
97                 x = p1 * sin(x + b) * sin(x + b) + p2;
98                 derivation = 2.0 *p1 *sin(x + b) * cos(x + b);
99                 if (derivation != 0.0) { index += log(fabs(derivation)); iter++; }
101                 x = p1 * sin(x + c) * sin(x + c) + p2;
102                 derivation = 2.0 *p1 *sin(x + c) * cos(x + c);
103                 if (derivation != 0.0) { index += log(fabs(derivation)); iter++; }
105                 index = (iter != 0) ? index / (float)(iter) : 0.0;
106                 index = index * nabla_main + index_pre * (1.0 - nabla_main);
107         }
109         return index;
110 }
113         color Pos_Color = color(1.0, 0.0, 0.0),
114         color Mid_Color = color(0.0, 0.0, 0.0),
115         color Neg_Color = color(0.0, 0.0, 1.0),
116         float Pre_Iteration = 0.0,
117         float Main_Iteration = 1.0,
118         float Pos_Scale = 0.5,
119         float Neg_Scale = 0.5,
120         float Param1 = 2.0,
121         float Param2 = 2.0,
122         int Fac_Type = 0,
123         int Render_Type = 2,
124         float Scale = 0.25,
125         point Pos = P,
126         output float Fac = 0.0,
127         output color Color = 0.0)
128 {
129         /* Calculate Texture */
130         float index = lyapunov(Pos * Scale, Pre_Iteration, Main_Iteration, Param1, Param2);
132         /* Calculate Color */
133         if (index > 0.0 && (Render_Type != 0)) {
134                 index *= Pos_Scale;
135                 if (index > 1.0) { index = 1.0; }
136                 Color = (Pos_Color - Mid_Color) * index + Mid_Color;
137         }
138         else if (index < 0.0 && (Render_Type != 1)) {
139                 index *= Neg_Scale;
140                 if (index < -1.0) { index = -1.0; }
141                 Color = (Mid_Color - Neg_Color) * index + Mid_Color;
142         }
143         else {
144                 Color = Mid_Color;
145         }
147         /* Adjust Index */
148         if (Fac_Type == 0) {
149                 index = 0.5 + index * 0.5;
150         }
151         else if (Fac_Type == 1) {
152                 index = fabs(index);
153         }
154         else if (Fac_Type == 2) {
155                 index = (Color + Color + Color) * (1.0 / 3.0);
156         }
158         Fac = index;
159 }