Domain Coloring

[_vince]

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 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

[bplus]

Very nice screen shots with the shading added.

[bplus]

@_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?

[_vince]

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


[bplus]

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.

[bplus]

Bugs found in mods by b+

[_vince]

do the images take forever to load?

[bplus]

do the images take forever to load?

About the same amount of time, you really have to download to see what I mean.

[_vince]

img

[bplus]

Spoiler ;-))

Looks like font didn't load?
 

[Richard Frost]

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.

[bplus]

And Richard's reminds me of moving Plasma points around in my plasma experiments.

specially fig. d

[_vince]

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

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


[_vince]

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


[bplus]

Is that an edit of the post above?