QB64.org Forum

Active Forums => Programs => Topic started by: _vince on August 15, 2020, 02:13:58 pm

Title: Domain Coloring
Post by: _vince on August 15, 2020, 02:13:58 pm
A kind of method for plotting complex functions:
https://en.wikipedia.org/wiki/Domain_coloring (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
Title: Re: Domain Coloring
Post by: _vince on August 15, 2020, 03:23:42 pm
Pretty enough to take pictures
Title: Re: Domain Coloring
Post by: bplus on August 15, 2020, 03:33:36 pm
@_vince  Nice! How did you come across this?
Title: Re: Domain Coloring
Post by: _vince on August 15, 2020, 03:42:55 pm
I found the greatest book ever https://www.springer.com/gp/book/9783034801799 (https://www.springer.com/gp/book/9783034801799)

Lots of pretty pictures, I can email you a pdf copy if you want
Title: Re: Domain Coloring
Post by: euklides on August 16, 2020, 04:18:04 am
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...
Title: Re: Domain Coloring
Post by: Ashish on August 16, 2020, 10:24:26 am
Beautilful! I like it.
Title: Re: Domain Coloring
Post by: bplus on August 16, 2020, 11:59:04 am
:) Almost has Mandelbrot looking like exotic butterfly.

and sounds just like "f(f(f(f(f(f(f(f(f(f(z))))))))))"
Title: Re: Domain Coloring
Post by: _vince on August 17, 2020, 03:01:42 pm
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
Title: Re: Domain Coloring
Post by: _vince on August 17, 2020, 03:12:50 pm
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
Title: Re: Domain Coloring
Post by: bplus on August 17, 2020, 03:50:11 pm
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.  
Title: Re: Domain Coloring
Post by: _vince on August 17, 2020, 03:54:50 pm
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 (https://en.wikipedia.org/wiki/File:Domain_coloring_x2-1_x-2-i_x-2-i_d_x2%2B2%2B2i.xcf)
Title: Re: Domain Coloring
Post by: bplus on August 17, 2020, 04:10:29 pm
I discovered an effect while messing with the colors, it's a moving experience ;-))
Code: QB64: [Select]
  1. DEFDBL A-Z
  2. pi = 4 * ATN(1)
  3.  
  4. CONST sw = 800
  5. CONST sh = 600
  6. _TITLE "b+ fiddles... press spacebar "
  7.  
  8. SCREEN _NEWIMAGE(sw, sh, 32)
  9.  
  10.  
  11. 'dim as double x, y, u, v, m, a, uu, vv, x0, y0
  12.  
  13.  
  14. FOR j = 6 TO 6
  15.     FOR yy = 0 TO sh
  16.         FOR xx = 0 TO sw
  17.             SELECT CASE j
  18.                 CASE 0
  19.                     u = (xx - sw / 2) * 0.015
  20.                     v = (sh / 2 - yy) * 0.015
  21.  
  22.                     x = u
  23.                     y = v
  24.  
  25.                     p = "z"
  26.                 CASE 1
  27.                     u = (xx - sw / 2) * 0.01 + 2
  28.                     v = (sh / 2 - yy) * 0.01
  29.  
  30.                     x = EXP(u) * COS(v)
  31.                     y = EXP(u) * SIN(v)
  32.  
  33.                     p = "exp(z)"
  34.                 CASE 2
  35.                     u = (xx - sw / 2) * 0.015
  36.                     v = (sh / 2 - yy) * 0.015
  37.  
  38.                     x = SIN(u) * _COSH(v)
  39.                     y = COS(u) * _SINH(v)
  40.  
  41.                     p = "sin(z)"
  42.                 CASE 3
  43.                     u = (xx - sw / 2) * 0.005 - 0.5
  44.                     v = (sh / 2 - yy) * 0.005
  45.                     x0 = u
  46.                     y0 = v
  47.                     FOR i = 0 TO 0
  48.                         uu = u * u - v * v + x0
  49.                         vv = 2 * u * v + y0
  50.  
  51.                         u = uu
  52.                         v = vv
  53.                     NEXT
  54.                     x = u
  55.                     y = v
  56.  
  57.                     p = "f(z) = z^2 + c"
  58.                 CASE 4
  59.                     u = (xx - sw / 2) * 0.005 - 0.5
  60.                     v = (sh / 2 - yy) * 0.005
  61.                     x0 = u
  62.                     y0 = v
  63.                     FOR i = 0 TO 1
  64.                         uu = u * u - v * v + x0
  65.                         vv = 2 * u * v + y0
  66.  
  67.                         u = uu
  68.                         v = vv
  69.                     NEXT
  70.                     x = u
  71.                     y = v
  72.  
  73.                     p = "f(f(z))"
  74.                 CASE 5
  75.                     u = (xx - sw / 2) * 0.005 - 0.5
  76.                     v = (sh / 2 - yy) * 0.005
  77.                     x0 = u
  78.                     y0 = v
  79.                     FOR i = 0 TO 3
  80.                         uu = u * u - v * v + x0
  81.                         vv = 2 * u * v + y0
  82.  
  83.                         u = uu
  84.                         v = vv
  85.                     NEXT
  86.                     x = u
  87.                     y = v
  88.  
  89.                     p = "f(f(f(f(z))))"
  90.                 CASE 6
  91.                     u = (xx - sw / 2) * 0.005 - 0.5
  92.                     v = (sh / 2 - yy) * 0.005
  93.                     x0 = u
  94.                     y0 = v
  95.                     FOR i = 0 TO 9
  96.                         uu = u * u - v * v + x0
  97.                         vv = 2 * u * v + y0
  98.  
  99.                         u = uu
  100.                         v = vv
  101.                     NEXT
  102.                     x = u
  103.                     y = v
  104.  
  105.                     p = "f(f(f(f(f(f(f(f(f(f(z))))))))))"
  106.                 CASE 7
  107.                     u = (xx - sw / 2) * 0.0015
  108.                     v = (sh / 2 - yy) * 0.0015
  109.  
  110.                     x = SIN(u / (u * u + v * v)) * _COSH(-v / (u * u + v * v))
  111.                     y = COS(u / (u * u + v * v)) * _SINH(-v / (u * u + v * v))
  112.  
  113.                     p = "sin(1/z)"
  114.             END SELECT
  115.  
  116.             m = SQR(x * x + y * y)
  117.             a = (pi + _ATAN2(y, x)) / (2 * pi)
  118.  
  119.             PSET (xx, yy), hrgb(a, m)
  120.  
  121.             LOCATE 1, 1: PRINT p
  122.         NEXT
  123.     NEXT
  124.     DO
  125.         _SCREENMOVE 200 + RND * 4 - 2, 80 + RND * 3 - 1.5
  126.         _LIMIT 30
  127.     LOOP UNTIL _KEYDOWN(32)
  128.     SLEEP
  129.  
  130.  
  131. FUNCTION hrgb (h AS DOUBLE, m AS DOUBLE)
  132.     'dim as double r, g, b, mm
  133.  
  134.     r = 0.5 - 0.5 * SIN(2 * pi * h - pi / 2)
  135.     g = (0.5 + 0.5 * SIN(2 * pi * h * 1.5 - pi / 2)) * -(h < 0.66)
  136.     b = (0.5 + 0.5 * SIN(2 * pi * h * 1.5 + pi / 2)) * -(h > 0.33)
  137.  
  138.     mm = (m * 100 MOD 100) * 0.01
  139.  
  140.     IF ((m * 2 * pi * 100 MOD 2 * pi * 100) < 50) OR ((h * 2 * pi * 100 MOD 4 * 2 * pi) < 6) THEN
  141.         'hrgb = _rgb(0,0,0)
  142.         hrgb = _RGB(115 * r, 115 * g, 115 * b)
  143.     ELSE
  144.         hrgb = _RGB(255 * r, 255 * g, 255 * b)
  145.     END IF
  146.  
  147.  
Title: Re: Domain Coloring
Post by: _vince on August 17, 2020, 07:24:11 pm
edit:
Title: Re: Domain Coloring
Post by: _vince on August 19, 2020, 12:54:08 pm
https://en.wikipedia.org/wiki/Cauchy%27s_integral_formula (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 (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 (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. pi = 4 * ATN(1)
  8.  
  9. 'declare function cosh(x as double) as double
  10. 'declare function sinh(x as double) as double
  11.  
  12. 'screenres sw, sh, 32
  13. SCREEN _NEWIMAGE(sw, sh, 32)
  14.  
  15.  
  16. n = 30
  17. DIM sum AS DOUBLE, delta AS DOUBLE, p AS DOUBLE, q AS DOUBLE, pp AS DOUBLE, qq AS DOUBLE
  18.  
  19.  
  20. FOR i = 0 TO 2
  21.  
  22.     FOR yy = 0 TO sh
  23.         FOR xx = 0 TO sw
  24.    
  25.             SELECT CASE i
  26.                 CASE 0
  27.                     u = (xx - sw / 2) * 0.015
  28.                     v = (sh / 2 - yy) * 0.015
  29.                     x = u
  30.                     y = v
  31.                 CASE 1
  32.                     u = (xx - sw / 2) * 0.01 + 3 * pi / 2
  33.                     v = (sh / 2 - yy) * 0.01
  34.  
  35.                     x = u
  36.                     y = v
  37.        
  38.                     cexp uu, vv, u, v, 3, 0
  39.                     x = x - uu / 6
  40.                     y = y - vv / 6
  41.        
  42.                     cexp uu, vv, u, v, 5, 0
  43.                     x = x + uu / 120
  44.                     y = y + vv / 120
  45.        
  46.                     cexp uu, vv, u, v, 7, 0
  47.                     x = x - uu / 5040
  48.                     y = y - vv / 5040
  49.                 CASE 2
  50.                     u = (xx - sw / 2) * 0.01 + 3 * pi / 2
  51.                     v = (sh / 2 - yy) * 0.01
  52.  
  53.                     x = u
  54.                     y = v
  55.        
  56.                     cexp uu, vv, u, v, 3, 0
  57.                     x = x - uu / 6
  58.                     y = y - vv / 6
  59.        
  60.                     cexp uu, vv, u, v, 5, 0
  61.                     x = x + uu / 120
  62.                     y = y + vv / 120
  63.  
  64.                     cexp uu, vv, u, v, 7, 0
  65.                     x = x - uu / 5040
  66.                     y = y - vv / 5040
  67.        
  68.                     sum = 0
  69.                     FOR j = 0 TO n - 1
  70.                         p = pi * COS(j * 2 * pi / n) + 3 * pi / 2
  71.                         q = pi * SIN(j * 2 * pi / n)
  72.            
  73.                         pp = p
  74.                         qq = q
  75.                         cexp uu, vv, p, q, 3, 0
  76.                         pp = pp - uu / 6
  77.                         qq = qq - vv / 6
  78.                         cexp uu, vv, p, q, 5, 0
  79.                         pp = pp + uu / 120
  80.                         qq = qq + vv / 120
  81.                         cexp uu, vv, p, q, 7, 0
  82.                         pp = pp - uu / 5040
  83.                         qq = qq - vv / 5040
  84.            
  85.                         cdiv uu, vv, pp, qq, p - u, q - v
  86.            
  87.                         IF j = 0 OR j = n - 1 THEN
  88.                             sum = sum + 0.5 * (-uu * pi * SIN(j * 2 * pi / n) - vv * pi * COS(j * 2 * pi / n))
  89.                         ELSE
  90.                             sum = sum + (-uu * pi * SIN(j * 2 * pi / n) - vv * pi * COS(j * 2 * pi / n))
  91.                         END IF
  92.                     NEXT
  93.                     x = sum * 2 * pi / n
  94.        
  95.                     sum = 0
  96.                     FOR j = 0 TO n - 1
  97.                         p = pi * COS(j * 2 * pi / n) + 3 * pi / 2
  98.                         q = pi * SIN(j * 2 * pi / n)
  99.            
  100.                         pp = p
  101.                         qq = q
  102.                         cexp uu, vv, p, q, 3, 0
  103.                         pp = pp - uu / 6
  104.                         qq = qq - vv / 6
  105.                         cexp uu, vv, p, q, 5, 0
  106.                         pp = pp + uu / 120
  107.                         qq = qq + vv / 120
  108.                         cexp uu, vv, p, q, 7, 0
  109.                         pp = pp - uu / 5040
  110.                         qq = qq - vv / 5040
  111.            
  112.                         cdiv uu, vv, pp, qq, p - u, q - v
  113.            
  114.                         IF j = 0 OR j = n - 1 THEN
  115.                             sum = sum + 0.5 * (-vv * pi * SIN(j * 2 * pi / n) + uu * pi * COS(j * 2 * pi / n))
  116.                         ELSE
  117.                             sum = sum + (-vv * pi * SIN(j * 2 * pi / n) + uu * pi * COS(j * 2 * pi / n))
  118.                         END IF
  119.                     NEXT
  120.                     y = sum * 2 * pi / n
  121.        
  122.                     cmul uu, vv, x, y, 0, -1 / (2 * pi)
  123.        
  124.                     x = uu
  125.                     y = vv
  126.  
  127.  
  128.             END SELECT
  129.  
  130.             m = SQR(x * x + y * y)
  131.             a = (pi + _ATAN2(y, x)) / (2 * pi)
  132.  
  133.             PSET (xx, yy), hrgb(a, m)
  134.         NEXT
  135.     NEXT
  136.  
  137.     CIRCLE (sw / 2, sh / 2), pi / 0.01 - 2, _RGB(0, 0, 0)
  138.     CIRCLE (sw / 2, sh / 2), pi / 0.01 - 1, _RGB(0, 0, 0)
  139.     CIRCLE (sw / 2, sh / 2), pi / 0.01, _RGB(255, 255, 0)
  140.     CIRCLE (sw / 2, sh / 2), pi / 0.01 + 1, _RGB(0, 0, 0)
  141.     CIRCLE (sw / 2, sh / 2), pi / 0.01 + 2, _RGB(0, 0, 0)
  142.  
  143.     'for a = 0 to 2*pi step 0.01
  144.     '    x = pi*cos(a) + 3*pi/2
  145.     '    y = pi*sin(a)
  146.  
  147.     '    pset ((x/0.015 + sw/2), (sh/2 - y/0.015)), rgb(255,255,255)
  148.     'next
  149.  
  150.     SLEEP
  151.  
  152.  
  153. FUNCTION hrgb (h AS DOUBLE, m AS DOUBLE)
  154.     DIM r AS DOUBLE, g AS DOUBLE, b AS DOUBLE, t AS DOUBLE
  155.  
  156.     r = 0.5 - 0.5 * SIN(2 * pi * h - pi / 2)
  157.     g = (0.5 + 0.5 * SIN(2 * pi * h * 1.5 - pi / 2)) * -(h < 0.66)
  158.     b = (0.5 + 0.5 * SIN(2 * pi * h * 1.5 + pi / 2)) * -(h > 0.33)
  159.  
  160.     n = 16
  161.     't = 0.3/m
  162.     t = 0.2
  163.     h = h + (0.5 * t / n)
  164.  
  165.     IF ((m * 500 MOD 500) < 70) OR (ABS((h * n) - INT(h * n)) < t) THEN
  166.         'hrgb = rgb(0,0,0)
  167.         hrgb = _RGB(125 * r, 125 * g, 125 * b)
  168.     ELSE
  169.         hrgb = _RGB(255 * r, 255 * g, 255 * b)
  170.     END IF
  171.  
  172. SUB cexp (u AS DOUBLE, v AS DOUBLE, x AS DOUBLE, y AS DOUBLE, a AS DOUBLE, b AS DOUBLE)
  173.     DIM lnz AS DOUBLE, argz AS DOUBLE, mag AS DOUBLE, ang AS DOUBLE
  174.  
  175.     lnz = x * x + y * y
  176.  
  177.     IF lnz = 0 THEN
  178.         u = 0
  179.         v = 0
  180.     ELSE
  181.         lnz = 0.5 * LOG(lnz)
  182.         argz = _ATAN2(y, x)
  183.  
  184.         mag = EXP(a * lnz - b * argz)
  185.         ang = a * argz + b * lnz
  186.  
  187.         u = mag * COS(ang)
  188.         v = mag * SIN(ang)
  189.     END IF
  190.  
  191. SUB cmul (u AS DOUBLE, v AS DOUBLE, x AS DOUBLE, y AS DOUBLE, a AS DOUBLE, b AS DOUBLE)
  192.     u = x * a - y * b
  193.     v = x * b + y * a
  194.  
  195. SUB cdiv (u AS DOUBLE, v AS DOUBLE, x AS DOUBLE, y AS DOUBLE, a AS DOUBLE, b AS DOUBLE)
  196.     DIM d AS DOUBLE
  197.     d = a * a + b * b
  198.     u = (x * a + y * b) / d
  199.     v = (y * a - x * b) / d
  200.  
Title: Re: Domain Coloring
Post by: bplus on August 19, 2020, 01:10:47 pm
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
Title: Re: Domain Coloring
Post by: _vince on September 05, 2021, 08:52:41 am
Going to add a couple of pictures here for archival purposes

zp11.png :  this is the simple rational function f(z) = (z + 1)/(z - 1).  So if you solve f(z) = 0 you get:  z + 1 = 0 -> z = -1 which means there is a 'zero' at z = -1.  Then you solve f(z) = infinity you get z-1 = 0 -> z = 1 which is the so-called 'pole' at z = 1.  In the plot, you see perpendicular contours, which look like field lines produced by a magnet or a magnetic dipole -- these are the phase contours.  Then you see the concentric circles around the pole getting closer and closer together as the function magnitude tends towards infinity -- these are the magnitude contours.

zp21.png :  here I add an extra zero for a total of two zeros and one pole (simply multiply the rational function f(z) by (z - a) to add a zero at z=a,,, or divide by (z - b) to add a pole at z=b). You see the interesting and pretty behavior of the field lines

zpspiral.png and zpspiral2.png :  here I create a rational function now placing the zeros and poles in a spiral to produce these beautiful floral patterns.  Something like,  for k=0 to 1000 (or whatever),  multiply or divide by (z - {k cos k + i k sin k}) with, of course, all the appropriate scaling constants etc

* * *

cgamman10000.png :  is an attempt to approximate the complex gamma function by numerically integrating the following with trapezoid rule:

\int_0^1 (log(1/t))^(z - 1) dt, for Re{z} > 0

then using the reflection formula:

Gamma(z) Gamma(-z) = -pi/(z sin(pi z)) to get the negative real half of the complex plane

and if you compare to https://upload.wikimedia.org/wikipedia/commons/3/34/Gamma1.png (https://upload.wikimedia.org/wikipedia/commons/3/34/Gamma1.png) you see that there is a lot of integration error around the y-axis (or the imaginary axis), and that is likely due logarithmic or rational nature of the integrand, but further down the x-axis (or the real axis) you see it better resemble the URL.  Just as the post above on CIF, these integration errors look interesting and pretty

Title: Re: Domain Coloring
Post by: bplus on September 05, 2021, 10:37:11 am
Very nice screen shots with the shading added.
Title: Re: Domain Coloring
Post by: bplus on September 05, 2021, 11:18:41 am
@_vince  let me guess you are coding in FB and showing screen shots until perfected then will translate to QB64 so we have code for QB64 programs board?
Title: Re: Domain Coloring
Post by: _vince on September 05, 2021, 02:47:24 pm
the code isn't that interesting (just produces the image) but if anyone's interested:

Code: [Select]
defdbl a-z

const sw = 1024
const sh = 768

dim shared pi
pi = 4*atn(1)

screen _newimage(sw, sh, 32)

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

                        x = u
                        y = v

                case 1
                        u = (xx - sw/2)*0.003
                        v = (sh/2 - yy)*0.003

                        cdiv uu, vv, u + 1, v, u - 1, v

                        x = uu
                        y = vv

                case 2
                        u = (xx - sw/2)*0.008
                        v = (sh/2 - yy)*0.008

                        cmul uu, vv, u + 0.707, v + 0.707, u + 0.707, v - 0.707
                        cdiv x, y, uu, vv, u - 1, v

                case 3
                        u = (xx - sw/2)*0.035
                        v = (sh/2 - yy)*0.035

                        uu = u
                        vv = v

                        a = 0
                        for t = 0 to 100 step 1
                                dx = 0.1*t*cos(100*t)
                                dy = 0.1*t*sin(100*t)

                                if a then
                                        cmul uu, vv, uu, vv, u + dx, v + dy
                                else
                                        cdiv uu, vv, uu, vv, u + dx, v + dy
                                end if

                                a = a xor 1
                        next

                        x = uu
                        y = vv

                end select

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

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

        sleep
next

system

function hrgb&(h, m)
        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)

        dim n as long
    n = 16
   
    mm = m*500 mod 500
        p = abs((h*n) - int(h*n))

        rr = 255*r - 0.15*mm - 35*p
        gg = 255*g - 0.15*mm - 35*p
        bb = 255*b - 0.15*mm - 35*p

        if rr < 0 then rr = 0
        if gg < 0 then gg = 0
        if bb < 0 then bb = 0

        hrgb& = _rgb(rr, gg, bb)
end function

function cosh(x)
    cosh = 0.5*(exp(x) + exp(-x))
end function

function sinh(x)
    sinh = 0.5*(exp(x) - exp(-x))
end function

'u + iv = (x + iy)*(a + ib)
sub cmul(u, v, xx, yy, aa, bb)
        x = xx
        y = yy
        a = aa
        b = bb
    u = x*a - y*b
    v = x*b + y*a
end sub

'u + iv = (x + iy)/(a + ib)
sub cdiv(u, v, xx, yy, aa, bb)
        x = xx
        y = yy
        a = aa
        b = bb
    d = a*a + b*b
    u = (x*a + y*b)/d
    v = (y*a - x*b)/d
end sub
Title: Re: Domain Coloring
Post by: bplus on September 05, 2021, 09:19:05 pm
Thanks @_vince

I am looking at the first image and thinking, that's it? Then I am reading over the code, see a Sleep and say "Oh there's more!" very nice, I was hoping to see those more complex images on computer screen.
Title: Re: Domain Coloring
Post by: bplus on September 05, 2021, 10:26:22 pm
Bugs found in mods by b+
Title: Re: Domain Coloring
Post by: _vince on September 06, 2021, 05:50:54 pm
do the images take forever to load?
Title: Re: Domain Coloring
Post by: bplus on September 06, 2021, 06:06:58 pm
do the images take forever to load?

About the same amount of time, you really have to download to see what I mean.
Title: Re: Domain Coloring
Post by: _vince on September 06, 2021, 06:42:31 pm
img
Title: Re: Domain Coloring
Post by: bplus on September 06, 2021, 07:40:01 pm
Spoiler ;-))

Looks like font didn't load?
  [ This attachment cannot be displayed inline in 'Print Page' view ]  
Title: Re: Domain Coloring
Post by: Richard Frost on September 08, 2021, 09:23:21 am
Some of those images remind me of "method of moments", which is all
about why marshmallows get bigger in a microwave oven.  Check the Wiki
article on the quoted phrase.
Title: Re: Domain Coloring
Post by: bplus on September 08, 2021, 12:10:12 pm
And Richard's reminds me of moving Plasma points around in my plasma experiments.

specially fig. d
Title: Re: Domain Coloring
Post by: _vince on September 12, 2021, 03:32:51 pm
More plots, this time using logarithmic magnitude contouring, or log magnitude.

notes:

Code: [Select]
defdbl a-z

const sw = 1024
const sh = 768

dim shared pi
pi = 4*atn(1)

screen _newimage(sw, sh, 32)

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

                        x = u
                        y = v
                case 1
                        u = (xx - sw/2)*0.015
                        v = (sh/2 - yy)*0.015

                        'x = exp(u)*cos(v)
                        'y = exp(u)*sin(v)
                        cexp x, y, exp(1), 0, u, v

                case 2
                        u = (xx - sw/2)*0.005
                        v = (sh/2 - yy)*0.005

                        cdiv uu, vv, u + 1, v, u - 1, v

                        x = uu
                        y = vv

                case 3
                        u = (xx - sw/2)*0.008
                        v = (sh/2 - yy)*0.008

                        cmul uu, vv, u + 1, v + 1, u + 1, v - 1
                        cdiv x, y, uu, vv, u - 1, v

                case 4
                        u = (xx - sw/2)*0.015
                        v = (sh/2 - yy)*0.015

                        x = sin(u)*cosh(v)
                        y = cos(u)*sinh(v)

                case 5
                        u = (xx - sw/2)*0.0035
                        v = (sh/2 - yy)*0.0035

                        'cdiv u, v, 1, 0, u, v
                        x = sin(u)*cosh(v)
                        y = cos(u)*sinh(v)

                case 6
                        u = (xx - sw/2)*0.01
                        v = (sh/2 - yy)*0.01

                        uu = exp(u)*cos(v)
                        vv = exp(u)*sin(v)

                        x = sin(uu)*cosh(vv)
                        y = cos(uu)*sinh(vv)
                case 7
                        u = (xx - sw/2)*0.01
                        v = (sh/2 - yy)*0.01

                        uu = sin(u)*cosh(v)
                        vv = cos(u)*sinh(v)

                        x = exp(uu)*cos(vv)
                        y = exp(uu)*sin(vv)

                case 8
                        u = (xx - sw/2)*0.0035
                        v = (sh/2 - yy)*0.0035

                        uu = exp(u)*cos(v)
                        vv = exp(u)*sin(v)

                        x = sin(uu)*cosh(vv)
                        y = cos(uu)*sinh(vv)

                        uu = exp(x)*cos(y)
                        vv = exp(x)*sin(y)

                        cdiv x, y, 1, 0, uu, vv
                end select

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

                'log magnitude
                pset (xx, yy), hrgb&(a, log(5000*m + 1))
        next
        next

        locate 1,1
        select case i
        case 0: ? "w = z"
        case 1: ? "w = exp(z)"
        case 2: ? "w = (z + 1)/(z - 1)"
        case 3: ? "w = (z + 1 + i)*(z + 1 - i)/(z - 1)"
        case 4: ? "w = sin(z)"
        case 5: ? "w = sin(z) zoomed in"
        case 6: ? "w = sin(exp(z))"
        case 7: ? "w = exp(sin(z))"
        case 8: ? "w = exp(sin(exp(z)))"
        end select

        sleep
        if inkey$ = chr$(27) then exit for
next

system

function hrgb&(h as double, m as double)
        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*500 mod 500

        dim n as long
    n = 16
        p = abs((h*n) - int(h*n))

        rr = 255*r - 0.15*mm - 35*p
        gg = 255*g - 0.15*mm - 35*p
        bb = 255*b - 0.15*mm - 35*p

        if rr < 0 then rr = 0
        if gg < 0 then gg = 0
        if bb < 0 then bb = 0

        hrgb& = _rgb(rr, gg, bb)
end function

'u + iv = (x + iy)^(a + ib)
sub cexp(u, v, xx, yy, aa, bb)
        x = xx
        y = yy
        a = aa
        b = bb
   
    lnz = x*x + y*y
   
    if lnz = 0 then
        u = 0
        v = 0
    else
        lnz = 0.5*log(lnz)
        argz = _atan2(y, x)
       
        mag = exp(a*lnz - b*argz)
        ang = a*argz + b*lnz
       
        u = mag*cos(ang)
        v = mag*sin(ang)
    end if
end sub

function cosh(x)
    cosh = 0.5*(exp(x) + exp(-x))
end function

function sinh(x)
    sinh = 0.5*(exp(x) - exp(-x))
end function

'u + iv = (x + iy)*(a + ib)
sub cmul(u, v, xx, yy, aa, bb)
        x = xx
        y = yy
        a = aa
        b = bb
    u = x*a - y*b
    v = x*b + y*a
end sub

'u + iv = (x + iy)/(a + ib)
sub cdiv(u, v, xx, yy, aa, bb)
        x = xx
        y = yy
        a = aa
        b = bb
    d = a*a + b*b
    u = (x*a + y*b)/d
    v = (y*a - x*b)/d
end sub
Title: Re: Domain Coloring
Post by: _vince on September 19, 2021, 06:34:27 pm
More plots, this time using logarithmic magnitude contouring, or log magnitude.

notes:
  • in the second plot w=exp(z) you see that the phase and mag contours are parallel and aligned with x and y axis, like a transformation to cartesian coordinates from polar
  • in the 'dipole' plot, w=(z + 1)/(z - 1) you see in the log mag that phase contours or concentric circles are roughly equally spaced apart around both the zero and pole, but the gradient is reversed (almost like a real magnetic dipole)
  • the zoomed in w=sin(z) plot resembles w=z showing the identity z=sin(z) for small z, better seen in log mag, but then farther away from the real axis you see it more resemble w=exp(z), interesting: expand (where z=x+iy, and w=u+iv):  w=sin(z) = {sin x cosh y} + i {cos x sinh y} = {(-i/2) * (exp(ix) - exp(-ix)) * (1/2) * (exp(y) + exp(-y))  +  (-1/2) * (exp(ix) + exp(-ix)) * (i/2) * (exp(y) - exp(-y))} = ...

Code: [Select]
defdbl a-z

const sw = 1024
const sh = 768

dim shared pi
pi = 4*atn(1)

screen _newimage(sw, sh, 32)

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

                        x = u
                        y = v
                case 1
                        u = (xx - sw/2)*0.015
                        v = (sh/2 - yy)*0.015

                        'x = exp(u)*cos(v)
                        'y = exp(u)*sin(v)
                        cexp x, y, exp(1), 0, u, v

                case 2
                        u = (xx - sw/2)*0.005
                        v = (sh/2 - yy)*0.005

                        cdiv uu, vv, u + 1, v, u - 1, v

                        x = uu
                        y = vv

                case 3
                        u = (xx - sw/2)*0.008
                        v = (sh/2 - yy)*0.008

                        cmul uu, vv, u + 1, v + 1, u + 1, v - 1
                        cdiv x, y, uu, vv, u - 1, v

                case 4
                        u = (xx - sw/2)*0.015
                        v = (sh/2 - yy)*0.015

                        x = sin(u)*cosh(v)
                        y = cos(u)*sinh(v)

                case 5
                        u = (xx - sw/2)*0.0035
                        v = (sh/2 - yy)*0.0035

                        'cdiv u, v, 1, 0, u, v
                        x = sin(u)*cosh(v)
                        y = cos(u)*sinh(v)

                case 6
                        u = (xx - sw/2)*0.01
                        v = (sh/2 - yy)*0.01

                        uu = exp(u)*cos(v)
                        vv = exp(u)*sin(v)

                        x = sin(uu)*cosh(vv)
                        y = cos(uu)*sinh(vv)
                case 7
                        u = (xx - sw/2)*0.01
                        v = (sh/2 - yy)*0.01

                        uu = sin(u)*cosh(v)
                        vv = cos(u)*sinh(v)

                        x = exp(uu)*cos(vv)
                        y = exp(uu)*sin(vv)

                case 8
                        u = (xx - sw/2)*0.0035
                        v = (sh/2 - yy)*0.0035

                        uu = exp(u)*cos(v)
                        vv = exp(u)*sin(v)

                        x = sin(uu)*cosh(vv)
                        y = cos(uu)*sinh(vv)

                        uu = exp(x)*cos(y)
                        vv = exp(x)*sin(y)

                        cdiv x, y, 1, 0, uu, vv
                end select

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

                'log magnitude
                pset (xx, yy), hrgb&(a, log(5000*m + 1))
        next
        next

        locate 1,1
        select case i
        case 0: ? "w = z"
        case 1: ? "w = exp(z)"
        case 2: ? "w = (z + 1)/(z - 1)"
        case 3: ? "w = (z + 1 + i)*(z + 1 - i)/(z - 1)"
        case 4: ? "w = sin(z)"
        case 5: ? "w = sin(z) zoomed in"
        case 6: ? "w = sin(exp(z))"
        case 7: ? "w = exp(sin(z))"
        case 8: ? "w = exp(sin(exp(z)))"
        end select

        sleep
        if inkey$ = chr$(27) then exit for
next

system

function hrgb&(h as double, m as double)
        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*500 mod 500

        dim n as long
    n = 16
        p = abs((h*n) - int(h*n))

        rr = 255*r - 0.15*mm - 35*p
        gg = 255*g - 0.15*mm - 35*p
        bb = 255*b - 0.15*mm - 35*p

        if rr < 0 then rr = 0
        if gg < 0 then gg = 0
        if bb < 0 then bb = 0

        hrgb& = _rgb(rr, gg, bb)
end function

'u + iv = (x + iy)^(a + ib)
sub cexp(u, v, xx, yy, aa, bb)
        x = xx
        y = yy
        a = aa
        b = bb
   
    lnz = x*x + y*y
   
    if lnz = 0 then
        u = 0
        v = 0
    else
        lnz = 0.5*log(lnz)
        argz = _atan2(y, x)
       
        mag = exp(a*lnz - b*argz)
        ang = a*argz + b*lnz
       
        u = mag*cos(ang)
        v = mag*sin(ang)
    end if
end sub

function cosh(x)
    cosh = 0.5*(exp(x) + exp(-x))
end function

function sinh(x)
    sinh = 0.5*(exp(x) - exp(-x))
end function

'u + iv = (x + iy)*(a + ib)
sub cmul(u, v, xx, yy, aa, bb)
        x = xx
        y = yy
        a = aa
        b = bb
    u = x*a - y*b
    v = x*b + y*a
end sub

'u + iv = (x + iy)/(a + ib)
sub cdiv(u, v, xx, yy, aa, bb)
        x = xx
        y = yy
        a = aa
        b = bb
    d = a*a + b*b
    u = (x*a + y*b)/d
    v = (y*a - x*b)/d
end sub
Title: Re: Domain Coloring
Post by: bplus on September 19, 2021, 07:59:38 pm
Is that an edit of the post above?
Title: Re: Domain Coloring
Post by: _vince on September 19, 2021, 08:04:07 pm
exactly, just a few personal notes for ~archival~ reasons, but new rule is you can only modify within 5 minutes of posting (not eternally)
Title: Re: Domain Coloring
Post by: SMcNeill on September 19, 2021, 09:51:59 pm
exactly, just a few personal notes for ~archival~ reasons, but new rule is you can only modify within 5 minutes of posting (not eternally)

You do whatcha have to do, to keep things updated.  😘
Title: Re: Domain Coloring
Post by: _vince on September 28, 2021, 07:27:00 am
A few more pictures I want backed up in this thread. Also an example of a black and white coloring which i like