Domain Coloring

[_vince]

A kind of method for plotting complex functions:
https://en.wikipedia.org/wiki/Domain_coloring

So pretty I had to share this

Code: [Select]

defdbl a-z
dim shared pi as double
pi = 4*atn(1)

const sw = 800
const sh = 600

declare function hrgb(h as double, m as double)

screen _newimage(sw, sh, 32)


'dim as double x, y, u, v, m, a, uu, vv, x0, y0

dim p as string

for j=0 to 7
for yy=0 to sh
for xx=0 to sw
        select case j
        case 0
                u = (xx - sw/2)*0.015
                v = (sh/2 - yy)*0.015

                x = u
                y = v

                p = "z"
        case 1
                u = (xx - sw/2)*0.01 + 2
                v = (sh/2 - yy)*0.01

                x = exp(u)*cos(v)
                y = exp(u)*sin(v)

                p = "exp(z)"
        case 2
                u = (xx - sw/2)*0.015
                v = (sh/2 - yy)*0.015

                x = sin(u)*_cosh(v)
                y = cos(u)*_sinh(v)

                p = "sin(z)"
        case 3
                u = (xx - sw/2)*0.005 - 0.5
                v = (sh/2 - yy)*0.005
                x0 = u
                y0 = v
                for i=0 to 0
                        uu = u*u - v*v + x0
                        vv = 2*u*v + y0

                        u = uu
                        v = vv
                next
                x = u
                y = v

                p = "f(z) = z^2 + c"
        case 4
                u = (xx - sw/2)*0.005 - 0.5
                v = (sh/2 - yy)*0.005
                x0 = u
                y0 = v
                for i=0 to 1
                        uu = u*u - v*v + x0
                        vv = 2*u*v + y0

                        u = uu
                        v = vv
                next
                x = u
                y = v

                p = "f(f(z))"
        case 5
                u = (xx - sw/2)*0.005 - 0.5
                v = (sh/2 - yy)*0.005
                x0 = u
                y0 = v
                for i=0 to 3
                        uu = u*u - v*v + x0
                        vv = 2*u*v + y0

                        u = uu
                        v = vv
                next
                x = u
                y = v

                p = "f(f(f(f(z))))"
        case 6
                u = (xx - sw/2)*0.005 - 0.5
                v = (sh/2 - yy)*0.005
                x0 = u
                y0 = v
                for i=0 to 9
                        uu = u*u - v*v + x0
                        vv = 2*u*v + y0

                        u = uu
                        v = vv
                next
                x = u
                y = v

                p = "f(f(f(f(f(f(f(f(f(f(z))))))))))"
case 7
u = (xx - sw/2)*0.0015
v = (sh/2 - yy)*0.0015

x = sin(u/(u*u + v*v))*_cosh(-v/(u*u + v*v))
y = cos(u/(u*u + v*v))*_sinh(-v/(u*u + v*v))

p = "sin(1/z)"
        end select

        m = sqr(x*x + y*y)
        a = (pi + _atan2(y, x))/(2*pi)

        pset (xx, yy), hrgb(a, m)

        locate 1,1: ? p
next
next

sleep
next

system

function hrgb(h as double, m as double)
'dim as double r, g, b, mm

r = 0.5 - 0.5*sin(2*pi*h - pi/2)
    g = (0.5 + 0.5*sin(2*pi*h*1.5 - pi/2)) * -(h < 0.66)
    b = (0.5 + 0.5*sin(2*pi*h*1.5 + pi/2)) * -(h > 0.33)

mm = (m*100 mod 100)*0.01

if ((m*2*pi*100 mod 2*pi*100) < 50) or ((h*2*pi*100 mod 4*2*pi) < 6) then
'hrgb = _rgb(0,0,0)
hrgb = _rgb(15*r, 155*g, 155*b)
else
hrgb = _rgb(255*r, 255*g, 255*b)
end if
end function


[_vince]

Pretty enough to take pictures

[bplus]

@_vince  Nice! How did you come across this?

[_vince]

I found the greatest book ever https://www.springer.com/gp/book/9783034801799

Lots of pretty pictures, I can email you a pdf copy if you want

[euklides]

Nice to see !
And you can change the cases, like...

= (xx - sw / 2) * 0.035 + SIN(RND * 6)
v = (sh / 2 - yy) * 0.015
x = SIN(u + 1) * _COSH(v)
y = COS(u - 1) * _SINH(v)

or

u = (xx - sw / 3) * 0.035 + SIN(RND * 6)
v = (sh / 2 - yy + u * 15) * 0.025
x = SIN(u + 1) * 1 / _COSH(v)
y= 1 / (COS(u - 1) * _SINH(v)) + 2

and so on...

[Ashish]

Beautilful! I like it.

[bplus]

:) Almost has Mandelbrot looking like exotic butterfly.

and sounds just like "f(f(f(f(f(f(f(f(f(f(z))))))))))"

[_vince]

Some more picture, Taylor expansion of sine

You can clearly see all the complex zeros and the number of them is the order of the polynomial, not something you'd see in a real plot

[_vince]

Moving average filter

Again, a fun image of how zeros and poles look, you see how magnitude contours get closer and closer together around poles as the magnitude quickly increases to infinity

[bplus]

Hi @_vince

Looking over code today. Why is red discriminated? Could this be a typo? There seems to be no lack of color?

From the hrgb sub in original post:

Code: QB64: [Select]

1.         if ((m*2*pi*100 mod 2*pi*100) < 50) or ((h*2*pi*100 mod 4*2*pi) < 6) then
2.                 'hrgb = _rgb(0,0,0)
3.                 hrgb = _rgb(15*r, 155*g, 155*b)  '  <<< what's with red, typo ?
4.         else
5.                 hrgb = _rgb(255*r, 255*g, 255*b)
6.         end if
7.  

[_vince]

Yes, just a typo.  That's why the contours don't have the red tint around red.

This colouring scheme isn't particularly good, At some point I'm intending to make a really nice one like in pictures like this, https://en.wikipedia.org/wiki/File:Domain_coloring_x2-1_x-2-i_x-2-i_d_x2%2B2%2B2i.xcf

[bplus]

I discovered an effect while messing with the colors, it's a moving experience ;-))

Code: QB64: [Select]

1. DEFDBL A-Z
2. DIM SHARED pi AS DOUBLE
3. pi = 4 * ATN(1)
4.  
5. CONST sw = 800
6. CONST sh = 600
7. _TITLE "b+ fiddles... press spacebar "
8. declare function hrgb(h as double, m as double)
9.  
10. SCREEN _NEWIMAGE(sw, sh, 32)
11.  
12.  
13. 'dim as double x, y, u, v, m, a, uu, vv, x0, y0
14.  
15. DIM p AS STRING
16.  
17. FOR j = 6 TO 6
18.     FOR yy = 0 TO sh
19.         FOR xx = 0 TO sw
20.             SELECT CASE j
21.                 CASE 0
22.                     u = (xx - sw / 2) * 0.015
23.                     v = (sh / 2 - yy) * 0.015
24.  
25.                     x = u
26.                     y = v
27.  
28.                     p = "z"
29.                 CASE 1
30.                     u = (xx - sw / 2) * 0.01 + 2
31.                     v = (sh / 2 - yy) * 0.01
32.  
33.                     x = EXP(u) * COS(v)
34.                     y = EXP(u) * SIN(v)
35.  
36.                     p = "exp(z)"
37.                 CASE 2
38.                     u = (xx - sw / 2) * 0.015
39.                     v = (sh / 2 - yy) * 0.015
40.  
41.                     x = SIN(u) * _COSH(v)
42.                     y = COS(u) * _SINH(v)
43.  
44.                     p = "sin(z)"
45.                 CASE 3
46.                     u = (xx - sw / 2) * 0.005 - 0.5
47.                     v = (sh / 2 - yy) * 0.005
48.                     x0 = u
49.                     y0 = v
50.                     FOR i = 0 TO 0
51.                         uu = u * u - v * v + x0
52.                         vv = 2 * u * v + y0
53.  
54.                         u = uu
55.                         v = vv
56.                     NEXT
57.                     x = u
58.                     y = v
59.  
60.                     p = "f(z) = z^2 + c"
61.                 CASE 4
62.                     u = (xx - sw / 2) * 0.005 - 0.5
63.                     v = (sh / 2 - yy) * 0.005
64.                     x0 = u
65.                     y0 = v
66.                     FOR i = 0 TO 1
67.                         uu = u * u - v * v + x0
68.                         vv = 2 * u * v + y0
69.  
70.                         u = uu
71.                         v = vv
72.                     NEXT
73.                     x = u
74.                     y = v
75.  
76.                     p = "f(f(z))"
77.                 CASE 5
78.                     u = (xx - sw / 2) * 0.005 - 0.5
79.                     v = (sh / 2 - yy) * 0.005
80.                     x0 = u
81.                     y0 = v
82.                     FOR i = 0 TO 3
83.                         uu = u * u - v * v + x0
84.                         vv = 2 * u * v + y0
85.  
86.                         u = uu
87.                         v = vv
88.                     NEXT
89.                     x = u
90.                     y = v
91.  
92.                     p = "f(f(f(f(z))))"
93.                 CASE 6
94.                     u = (xx - sw / 2) * 0.005 - 0.5
95.                     v = (sh / 2 - yy) * 0.005
96.                     x0 = u
97.                     y0 = v
98.                     FOR i = 0 TO 9
99.                         uu = u * u - v * v + x0
100.                         vv = 2 * u * v + y0
101.  
102.                         u = uu
103.                         v = vv
104.                     NEXT
105.                     x = u
106.                     y = v
107.  
108.                     p = "f(f(f(f(f(f(f(f(f(f(z))))))))))"
109.                 CASE 7
110.                     u = (xx - sw / 2) * 0.0015
111.                     v = (sh / 2 - yy) * 0.0015
112.  
113.                     x = SIN(u / (u * u + v * v)) * _COSH(-v / (u * u + v * v))
114.                     y = COS(u / (u * u + v * v)) * _SINH(-v / (u * u + v * v))
115.  
116.                     p = "sin(1/z)"
117.             END SELECT
118.  
119.             m = SQR(x * x + y * y)
120.             a = (pi + _ATAN2(y, x)) / (2 * pi)
121.  
122.             PSET (xx, yy), hrgb(a, m)
123.  
124.             LOCATE 1, 1: PRINT p
125.         NEXT
126.     NEXT
127.     DO
128.         _SCREENMOVE 200 + RND * 4 - 2, 80 + RND * 3 - 1.5
129.         _LIMIT 30
130.     LOOP UNTIL _KEYDOWN(32)
131.     SLEEP
132. NEXT
133.  
134. SYSTEM
135.  
136. FUNCTION hrgb (h AS DOUBLE, m AS DOUBLE)
137.     'dim as double r, g, b, mm
138.  
139.     r = 0.5 - 0.5 * SIN(2 * pi * h - pi / 2)
140.     g = (0.5 + 0.5 * SIN(2 * pi * h * 1.5 - pi / 2)) * -(h < 0.66)
141.     b = (0.5 + 0.5 * SIN(2 * pi * h * 1.5 + pi / 2)) * -(h > 0.33)
142.  
143.     mm = (m * 100 MOD 100) * 0.01
144.  
145.     IF ((m * 2 * pi * 100 MOD 2 * pi * 100) < 50) OR ((h * 2 * pi * 100 MOD 4 * 2 * pi) < 6) THEN
146.         'hrgb = _rgb(0,0,0)
147.         hrgb = _RGB(115 * r, 115 * g, 115 * b)
148.     ELSE
149.         hrgb = _RGB(255 * r, 255 * g, 255 * b)
150.     END IF
151. END FUNCTION
152.  
153.  

[_vince]

edit:

[_vince]

https://en.wikipedia.org/wiki/Cauchy%27s_integral_formula

Cauchy's integral formula is a powerful result in complex analysis that says you can evaluate a function f(z0) inside a contour just by taking the contour integral of f(z)/(z - z0).  It essentially says that the function inside the contour can be defined entirely by the contour.

To begin, in the first image I plot the function

f(z) = z - z^3/3! + z^5/5! - z^7/7!

and draw the contour:

x(t) = pi*cos(t) + 3*pi/2
y(t) = pi*sin(t)

Then I numerically evaluate the following line integral using trapezoidal rule: https://en.wikipedia.org/wiki/Trapezoidal_rule

f(z0) = (1/(2*pi*i)) \cint_C f(z)/(z - z0) dz

Firstly for N = 30 and you see all the distortion around the contour and then for N = 100 and N = 1000 which is dead on.  And you see that it converges to the original function as you'd expect. The pretty artifacts around the contour are the errors due to numerical integration, it is like a complex visual of https://upload.wikimedia.org/wikipedia/commons/d/d1/Integration_num_trapezes_notation.svg

Code: QB64: [Select]

1. '#lang "fblite"
2. DEFLNG A-Z
3.  
4. CONST sw = 800
5. CONST sh = 600
6.  
7. DIM SHARED pi AS DOUBLE
8. pi = 4 * ATN(1)
9.  
10. declare function hrgb(h as double, m as double)
11. 'declare function cosh(x as double) as double
12. 'declare function sinh(x as double) as double
13. declare sub cexp(u as double, v as double, x as double, y as double, a as double, b as double)
14. declare sub cmul(u as double, v as double, x as double, y as double, a as double, b as double)
15. declare sub cdiv(u as double, v as double, x as double, y as double, a as double, b as double)
16.  
17. 'screenres sw, sh, 32
18. SCREEN _NEWIMAGE(sw, sh, 32)
19.  
20. DIM m AS DOUBLE, a AS DOUBLE
21. DIM x AS DOUBLE, y AS DOUBLE, u AS DOUBLE, v AS DOUBLE, uu AS DOUBLE, vv AS DOUBLE
22.  
23. n = 30
24. DIM sum AS DOUBLE, delta AS DOUBLE, p AS DOUBLE, q AS DOUBLE, pp AS DOUBLE, qq AS DOUBLE
25.  
26.  
27. FOR i = 0 TO 2
28.  
29.     FOR yy = 0 TO sh
30.         FOR xx = 0 TO sw
31.    
32.             SELECT CASE i
33.                 CASE 0
34.                     u = (xx - sw / 2) * 0.015
35.                     v = (sh / 2 - yy) * 0.015
36.                     x = u
37.                     y = v
38.                 CASE 1
39.                     u = (xx - sw / 2) * 0.01 + 3 * pi / 2
40.                     v = (sh / 2 - yy) * 0.01
41.  
42.                     x = u
43.                     y = v
44.        
45.                     cexp uu, vv, u, v, 3, 0
46.                     x = x - uu / 6
47.                     y = y - vv / 6
48.        
49.                     cexp uu, vv, u, v, 5, 0
50.                     x = x + uu / 120
51.                     y = y + vv / 120
52.        
53.                     cexp uu, vv, u, v, 7, 0
54.                     x = x - uu / 5040
55.                     y = y - vv / 5040
56.                 CASE 2
57.                     u = (xx - sw / 2) * 0.01 + 3 * pi / 2
58.                     v = (sh / 2 - yy) * 0.01
59.  
60.                     x = u
61.                     y = v
62.        
63.                     cexp uu, vv, u, v, 3, 0
64.                     x = x - uu / 6
65.                     y = y - vv / 6
66.        
67.                     cexp uu, vv, u, v, 5, 0
68.                     x = x + uu / 120
69.                     y = y + vv / 120
70.  
71.                     cexp uu, vv, u, v, 7, 0
72.                     x = x - uu / 5040
73.                     y = y - vv / 5040
74.        
75.                     sum = 0
76.                     FOR j = 0 TO n - 1
77.                         p = pi * COS(j * 2 * pi / n) + 3 * pi / 2
78.                         q = pi * SIN(j * 2 * pi / n)
79.            
80.                         pp = p
81.                         qq = q
82.                         cexp uu, vv, p, q, 3, 0
83.                         pp = pp - uu / 6
84.                         qq = qq - vv / 6
85.                         cexp uu, vv, p, q, 5, 0
86.                         pp = pp + uu / 120
87.                         qq = qq + vv / 120
88.                         cexp uu, vv, p, q, 7, 0
89.                         pp = pp - uu / 5040
90.                         qq = qq - vv / 5040
91.            
92.                         cdiv uu, vv, pp, qq, p - u, q - v
93.            
94.                         IF j = 0 OR j = n - 1 THEN
95.                             sum = sum + 0.5 * (-uu * pi * SIN(j * 2 * pi / n) - vv * pi * COS(j * 2 * pi / n))
96.                         ELSE
97.                             sum = sum + (-uu * pi * SIN(j * 2 * pi / n) - vv * pi * COS(j * 2 * pi / n))
98.                         END IF
99.                     NEXT
100.                     x = sum * 2 * pi / n
101.        
102.                     sum = 0
103.                     FOR j = 0 TO n - 1
104.                         p = pi * COS(j * 2 * pi / n) + 3 * pi / 2
105.                         q = pi * SIN(j * 2 * pi / n)
106.            
107.                         pp = p
108.                         qq = q
109.                         cexp uu, vv, p, q, 3, 0
110.                         pp = pp - uu / 6
111.                         qq = qq - vv / 6
112.                         cexp uu, vv, p, q, 5, 0
113.                         pp = pp + uu / 120
114.                         qq = qq + vv / 120
115.                         cexp uu, vv, p, q, 7, 0
116.                         pp = pp - uu / 5040
117.                         qq = qq - vv / 5040
118.            
119.                         cdiv uu, vv, pp, qq, p - u, q - v
120.            
121.                         IF j = 0 OR j = n - 1 THEN
122.                             sum = sum + 0.5 * (-vv * pi * SIN(j * 2 * pi / n) + uu * pi * COS(j * 2 * pi / n))
123.                         ELSE
124.                             sum = sum + (-vv * pi * SIN(j * 2 * pi / n) + uu * pi * COS(j * 2 * pi / n))
125.                         END IF
126.                     NEXT
127.                     y = sum * 2 * pi / n
128.        
129.                     cmul uu, vv, x, y, 0, -1 / (2 * pi)
130.        
131.                     x = uu
132.                     y = vv
133.  
134.  
135.             END SELECT
136.  
137.             m = SQR(x * x + y * y)
138.             a = (pi + _ATAN2(y, x)) / (2 * pi)
139.  
140.             PSET (xx, yy), hrgb(a, m)
141.         NEXT
142.     NEXT
143.  
144.     CIRCLE (sw / 2, sh / 2), pi / 0.01 - 2, _RGB(0, 0, 0)
145.     CIRCLE (sw / 2, sh / 2), pi / 0.01 - 1, _RGB(0, 0, 0)
146.     CIRCLE (sw / 2, sh / 2), pi / 0.01, _RGB(255, 255, 0)
147.     CIRCLE (sw / 2, sh / 2), pi / 0.01 + 1, _RGB(0, 0, 0)
148.     CIRCLE (sw / 2, sh / 2), pi / 0.01 + 2, _RGB(0, 0, 0)
149.  
150.     'for a = 0 to 2*pi step 0.01
151.     '    x = pi*cos(a) + 3*pi/2
152.     '    y = pi*sin(a)
153.  
154.     '    pset ((x/0.015 + sw/2), (sh/2 - y/0.015)), rgb(255,255,255)
155.     'next
156.  
157.     SLEEP
158. NEXT
159.  
160. SYSTEM
161.  
162. FUNCTION hrgb (h AS DOUBLE, m AS DOUBLE)
163.     DIM r AS DOUBLE, g AS DOUBLE, b AS DOUBLE, t AS DOUBLE
164.  
165.     r = 0.5 - 0.5 * SIN(2 * pi * h - pi / 2)
166.     g = (0.5 + 0.5 * SIN(2 * pi * h * 1.5 - pi / 2)) * -(h < 0.66)
167.     b = (0.5 + 0.5 * SIN(2 * pi * h * 1.5 + pi / 2)) * -(h > 0.33)
168.  
169.     n = 16
170.     't = 0.3/m
171.     t = 0.2
172.     h = h + (0.5 * t / n)
173.  
174.     IF ((m * 500 MOD 500) < 70) OR (ABS((h * n) - INT(h * n)) < t) THEN
175.         'hrgb = rgb(0,0,0)
176.         hrgb = _RGB(125 * r, 125 * g, 125 * b)
177.     ELSE
178.         hrgb = _RGB(255 * r, 255 * g, 255 * b)
179.     END IF
180. END FUNCTION
181.  
182. SUB cexp (u AS DOUBLE, v AS DOUBLE, x AS DOUBLE, y AS DOUBLE, a AS DOUBLE, b AS DOUBLE)
183.     DIM lnz AS DOUBLE, argz AS DOUBLE, mag AS DOUBLE, ang AS DOUBLE
184.  
185.     lnz = x * x + y * y
186.  
187.     IF lnz = 0 THEN
188.         u = 0
189.         v = 0
190.     ELSE
191.         lnz = 0.5 * LOG(lnz)
192.         argz = _ATAN2(y, x)
193.  
194.         mag = EXP(a * lnz - b * argz)
195.         ang = a * argz + b * lnz
196.  
197.         u = mag * COS(ang)
198.         v = mag * SIN(ang)
199.     END IF
200. END SUB
201.  
202. SUB cmul (u AS DOUBLE, v AS DOUBLE, x AS DOUBLE, y AS DOUBLE, a AS DOUBLE, b AS DOUBLE)
203.     u = x * a - y * b
204.     v = x * b + y * a
205. END SUB
206.  
207. SUB cdiv (u AS DOUBLE, v AS DOUBLE, x AS DOUBLE, y AS DOUBLE, a AS DOUBLE, b AS DOUBLE)
208.     DIM d AS DOUBLE
209.     d = a * a + b * b
210.     u = (x * a + y * b) / d
211.     v = (y * a - x * b) / d
212. END SUB
213.  

[bplus]

Hey _vince your graphs remind me of my study of Inversion Points of 2D Circle:
https://www.qb64.org/forum/index.php?topic=2933.msg121929#msg121929