Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
Vince's Corner Takeout
#33
Complex Number Library

This is a work in progress collection of operators and functions for complex numbers.  Mostly in the form of SUBs and FUNCTIONs at the end of the code.  The demo shows the library in use in several complex number topics with domain coloring plots as well as a simple example of how to plot the Mandelbrot set at the end.

Code: (Select All)
defdbl a-z

const sw = 800
const sh = 600

dim shared pi
pi = 4*atn(1)

zoom = 140

dim as long i, j, k, xx, yy

screen _newimage(sw, sh, 32)

for i=0 to 9
        '''plots
        for yy=0 to sh
        for xx=0 to sw

                x = (xx - sw/2)/zoom
                y = (sh/2 - yy)/zoom

                select case i
                case 0
                        u = x
                        v = y
                        pset (xx, yy), hrgb(u, v)
                case 1
                        u = x
                        v = y
                        pset (xx, yy), checker(u, v)
                case 2
                        'w = 1/z
                        cdiv u, v, 1, 0, x, y

                        'pset (xx, yy), hrgb(u, v)
                        pset (xx, yy), checker(u, v)

                case 3
                        'w = sin(1/z)

                        'extra zoom
                        d = 0.35*(x*x + y*y)
                        if d<>0 then
                                u = sin(x/d)*cosh(-y/d)
                                v = cos(x/d)*sinh(-y/d)
                        else
                                u = 0
                                v = 0
                        end if

                        pset (xx, yy), hrgb(u, v)
                        'pset (xx, yy), checker(u, v)

                case 4
                        'extra zoom
                        u = 0.56*x
                        v = 0.56*y
                        for j=0 to 14
                                uu = u*u - v*v + 0.35
                                vv = 2*u*v + 0.0

                                u = uu
                                v = vv
                        next

                        pset (xx, yy), hrgb(u, v)
                        'pset (xx, yy), checker(u, v)

                case 5
                        cmul u, v, 1, 0, x - cos(2*pi/3), y + sin(2*pi/3)
                        cmul u, v, u, v, x - cos(2*pi/3), y - sin(2*pi/3)
                        cmul u, v, u, v, x - 1, y
                        'cdiv u, v, u, v, x - 1, y

                        pset (xx, yy), hrgb(u, v)
                        'pset (xx, yy), checker(u, v)

                case 6
                        'CIF numerical integration
                        'f(z_0) = (2 pi i)^-1 int_C f(z)/(z - z0) dz

                        n = 35

                        uu = 0
                        vv = 0

                        for j=0 to n - 1
                                'C: z(t)
                                p = 1.5*cos(j*2*pi/n)
                                q = 1.5*sin(j*2*pi/n)

                                'f(z(t)):
                                cmul u, v, 1, 0, p - cos(2*pi/3), q + sin(2*pi/3)
                                cmul u, v, u, v, p - cos(2*pi/3), q - sin(2*pi/3)
                                cmul u, v, u, v, p - 1, q

                                'f(z)/(z - z0)
                                cdiv u, v, u, v, p - x, q - y

                                'dz/dt
                                cmul u, v, u, v, -1.5*sin(j*2*pi/n), 1.5*cos(j*2*pi/n)

                                if j = 0 or j = n - 1 then
                                        uu = uu + 0.5*u
                                        vv = vv + 0.5*v
                                else
                                        uu = uu + u
                                        vv = vv + v
                                end if
                        next
                        'dt
                        u = uu*2*pi/n
                        v = vv*2*pi/n

                        '1/(2 pi i)
                        cmul u, v, u, v, 0, -1/(2*pi)

                        pset (xx, yy), hrgb(u, v)
                        'pset (xx, yy), checker(u, v)

                case 7
                        'extra zoom
                        x = x*0.5
                        y = y*0.5

                        p = 1
                        q = 0

                        for j=0 to 5
                                cmul uu, vv, 1, 0, -0.4, -0.18*(j - 2.1)
                                cmul p, q, p, q, x - uu - 0.2, y - vv
                                cdiv p, q, p, q, x - uu + 0.2, y - vv + 0.1
                        next
                        for j=0 to 2
                                cmul uu, vv, 1, 0, 0.4, -0.18*(j - 2.1) - 0.18*2.1/2
                                cdiv p, q, p, q, x - uu - 0.2, y - vv
                                cmul p, q, p, q, x - uu + 0.2, y - vv + 0.1
                        next

                        u = p
                        v = q

                        pset (xx, yy), grey(u, v)

                case 8
                        'extra zoom
                        u = 0.66*x - 0.5
                        v = 0.66*y
                        x0 = u
                        y0 = v
                        for j=0 to 3
                                uu = u*u - v*v + x0
                                vv = 2*u*v + y0

                                u = uu
                                v = vv
                        next

                        'pset (xx, yy), hrgb(u, v)
                        pset (xx, yy), checker(u, v)
                case 9
                        'extra zoom
                        u = 0.66*x - 0.5
                        v = 0.66*y
                        x0 = u
                        y0 = v
                        for j=0 to 70
                                uu = u*u - v*v + x0
                                vv = 2*u*v + y0

                                u = uu
                                v = vv
                        next

                        'pset (xx, yy), hrgb(u, v)
                        pset (xx, yy), checker(u, v)
                end select
        next
        next

        '''diagrams
        select case i
        case 0
                _title "w = z polar contouring"
        case 1
                _title "w = z checkerboard"
        case 2
                _title "w = 1/z singularity"
        case 3
                _title "w = sin(1/z) essential singularity"
        case 4
                _title "Julia fractal"
        case 5
                _title "tri mass"
        case 6
                _title "Cauchy integral formula"
                a = 0
                x = 1.5*cos(a)
                y = 1.5*sin(a)
                circle (x*zoom + sw/2, sh/2 - y*zoom), 3, _rgb(255,255,0)

                for a=0 to 2*pi step 2*pi/n
                        x = 1.5*cos(a)
                        y = 1.5*sin(a)

                        line -(x*zoom + sw/2, sh/2 - y*zoom), _rgb(255,255,0)
                        circle step(0,0), 3, _rgb(255,255,0)
                next

        case 7
                _title "mutual inductance"

                sleep

                'extra zoom
                m = zoom/0.5

                'this diagram is a JB original
                'left
                a = -pi/2
                x = sw/2/m + 0.2*cos(a) - 0.4
                y = sh/2/m + 0.18*sin(a) + (a*0.18*0.5/pi) - 2.1*0.18
                for t=x to 0 step -0.001
                        circlef (x - t)*m, y*m, 1, _rgb(0,0,0)
                next
                for a = -pi/2 to 5*2*pi + 2*pi + pi/2 step 0.01
                        x = sw/2/m + 0.2*cos(a) - 0.4
                        y = sh/2/m + 0.18*sin(a) + (a*0.18*0.5/pi) - 2.1*0.18
                        circlef x*m, y*m, 1, _rgb(0,0,0)
                next
                a = 5*2*pi + 2*pi + pi/2
                x = sw/2/m + 0.2*cos(a) - 0.4
                y = sh/2/m + 0.18*sin(a) + (a*0.18*0.5/pi) - 2.1*0.18
                for t=x to 0 step -0.001
                        circlef (x - t)*m, y*m, 1, _rgb(0,0,0)
                next
                'right
                a = -pi/2
                x = sw/2/m - 0.2*cos(a) + 0.4
                y = sh/2/m + 0.18*sin(a) + (a*0.18*0.5)/pi - 2.1*0.18 + 0.18*2.1/4
                for t=0 to 1.5 step 0.001
                        circlef (x + t)*m, y*m, 1, _rgb(0,0,0)
                next
                for a = -pi/2 to 2*2*pi + 2*pi + pi/2 step 0.01
                        x = sw/2/m - 0.2*cos(a) + 0.4
                        y = sh/2/m + 0.18*sin(a) + (a*0.18*0.5)/pi - 2.1*0.18 + 0.18*2.1/4
                        circlef x*m, y*m, 1, _rgb(0,0,0)
                next
                a = 2*2*pi + 2*pi + pi/2
                x = sw/2/m - 0.2*cos(a) + 0.4
                y = sh/2/m + 0.18*sin(a) + (a*0.18*0.5)/pi - 2.1*0.18 + 0.18*2.1/4
                for t=0 to 1.5 step 0.001
                        circlef (x + t)*m, y*m, 1, _rgb(0,0,0)
                next

        case 8
                _title "checkerboard Mandelbrot"
        end select

        sleep
next

_title "escape time example with nth order Mandelbrot"

zoom = 210

for n=2 to 6
        'line (0, 0)-(sw, sh), _rgb(255,255,255), bf

        for yy=0 to sh
        for xx=0 to sw
                x = (xx - sw/2)/zoom - 0.5
                y = (sh/2 - yy)/zoom

                u = 0
                v = 0
                for i=0 to 140
                        'f(z) = z^n + c
                        cexp u, v, u, v, n, 0
                        u = u + x
                        v = v + y
                        if u*u + v*v > 4 then exit for
                next
                if i>=140 then
                        pset (xx, yy), _rgb(0,0,0)
                else
                        pset (xx, yy), _rgb(255,255,255)
                end if
        next
        next

        sleep
next

system


sub circlef(x as long, y as long, r as long, c as long)
        dim as long x0, y0, e
        x0 = r
        y0 = 0
        e = -r
        do while y0 < x0
                if e <=0 then
                        y0 = y0 + 1
                        line (x - x0, y + y0)-(x + x0, y + y0), c, bf
                        line (x - x0, y - y0)-(x + x0, y - y0), c, bf
                        e = e + 2*y0
                else
                        line (x - y0, y - x0)-(x + y0, y - x0), c, bf
                        line (x - y0, y + x0)-(x + y0, y + x0), c, bf
                        x0 = x0 - 1
                        e = e - 2*x0
                end if
        loop
        line (x - r, y)-(x + r, y), c, bf
end sub

function grey~&(x, y)
        m = sqr(x*x + y*y)
        a = (pi + _atan2(y, x))/(2*pi)

        m = log(1 + 100*m)

        'polar contouring
        n = 16
        mm = m*5000 mod 500
        p = abs(a*n - int(a*n))

        g = 1 - 0.0007*mm - 0.21*p

        grey = _rgb(255*g, 255*g, 255*g)
end function

function checker~&(xx, yy)
        if 1 then
                x = xx
                y = yy
        else 'polar checkerboard
                x = _atan2(yy, xx)/(pi/4)
                y = sqr(xx*xx + yy*yy)

                y = log(1 + 1000*y)
        end if

        z = abs(x - int(x)) xor abs(y - int(y))

        if z then checker = _rgb(0,0,0) else checker = _rgb(255,255,255)
end function

function hrgb~&(x, y)
        m = sqr(x*x + y*y)
        a = (pi + _atan2(y, x))/(2*pi)

        'm = log(1 + 1000*m)

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

        'polar contouring
        n = 16
        mm = m*500 mod 500
        p = abs(a*n - int(a*n))

        r = r - 0.0005*mm - 0.14*p
        g = g - 0.0005*mm - 0.14*p
        b = b - 0.0005*mm - 0.14*p

        'cartesian shading
        if 0 then
                t = 0.03 'thickness
                xx = abs(x - int(x)) < t or abs(-x - int(-x)) < t
                yy = abs(y - int(y)) < t or abs(-y - int(-y)) < t
                if xx or yy then
                'if m > 1 then 'dont shade origin
                        r = r - 0.5
                        g = g - 0.5
                        b = b - 0.5
                'end if
                end if
        end if

        hrgb = _rgb(255*r, 255*g, 255*b)
end function

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

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

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)
                m = exp(a*lnz - b*argz)
                a = a*argz + b*lnz
                u = m*cos(a)
                v = m*sin(a)
        end if
end sub

sub clog(u, v, xx, yy)
        x = xx
        y = yy
        lnz = x*x + y*y
        if lnz=0 then
                u = 0
                v = 0
        else
                u = 0.5*log(lnz)
                v = _atan2(y, x)
        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

sub csin(u, v, xx, yy)
        x = xx
        y = yy
        u = sin(x)*cosh(y)
        v = cos(x)*sinh(y)
end sub

sub ccos(u, v, xx, yy)
        x = xx
        y = yy
        u = cos(x)*cosh(y)
        v =-sin(x)*sinh(y)
end sub

function factorial~&(n)
        if n = 0 then
                factorial = 1
        else
                factorial = n*factorial(n - 1)
        end if
end function
Reply


Messages In This Thread
Vince's Corner Takeout - by bplus - 04-29-2022, 02:12 PM
RE: Vince's Corner Takeout - by vince - 04-29-2022, 09:34 PM
RE: Vince's Corner Takeout - by vince - 05-02-2022, 03:10 AM
RE: Vince's Corner Takeout - by bplus - 05-02-2022, 04:25 AM
RE: Vince's Corner Takeout - by vince - 05-02-2022, 11:16 PM
RE: Vince's Corner Takeout - by vince - 05-03-2022, 01:10 AM
RE: Vince's Corner Takeout - by bplus - 05-03-2022, 01:15 AM
RE: Vince's Corner Takeout - by vince - 05-03-2022, 04:26 AM
RE: Vince's Corner Takeout - by bplus - 05-03-2022, 03:32 PM
RE: Vince's Corner Takeout - by vince - 05-10-2022, 03:41 AM
RE: Vince's Corner Takeout - by vince - 05-10-2022, 03:57 AM
RE: Vince's Corner Takeout - by dcromley - 05-10-2022, 02:57 PM
RE: Vince's Corner Takeout - by vince - 05-10-2022, 08:14 PM
RE: Vince's Corner Takeout - by SMcNeill - 05-10-2022, 02:59 PM
RE: Vince's Corner Takeout - by vince - 05-11-2022, 01:13 AM
RE: Vince's Corner Takeout - by dcromley - 05-11-2022, 01:58 AM
RE: Vince's Corner Takeout - by vince - 06-01-2022, 09:05 AM
RE: Vince's Corner Takeout - by vince - 08-11-2022, 02:51 AM
RE: Vince's Corner Takeout - by bplus - 06-03-2022, 02:47 PM
RE: Vince's Corner Takeout - by triggered - 06-04-2022, 02:00 AM
RE: Vince's Corner Takeout - by vince - 06-07-2022, 02:02 AM
RE: Vince's Corner Takeout - by bplus - 06-07-2022, 02:15 AM
RE: Vince's Corner Takeout - by vince - 07-13-2022, 05:23 AM
RE: Vince's Corner Takeout - by BSpinoza - 07-14-2022, 04:54 AM
RE: Vince's Corner Takeout - by bplus - 07-14-2022, 04:35 PM
RE: Vince's Corner Takeout - by aurel - 08-11-2022, 01:02 PM
RE: Vince's Corner Takeout - by bplus - 08-11-2022, 04:22 PM
RE: Vince's Corner Takeout - by aurel - 08-11-2022, 05:33 PM
RE: Vince's Corner Takeout - by BSpinoza - 08-12-2022, 03:44 AM
RE: Vince's Corner Takeout - by vince - 08-11-2022, 08:42 PM
RE: Vince's Corner Takeout - by vince - 08-19-2022, 05:00 AM
RE: Vince's Corner Takeout - by bplus - 08-19-2022, 06:33 PM
RE: Vince's Corner Takeout - by vince - 08-23-2022, 10:04 PM
RE: Vince's Corner Takeout - by vince - 11-04-2022, 01:48 AM
RE: Vince's Corner Takeout - by vince - 03-31-2023, 11:07 PM
RE: Vince's Corner Takeout - by vince - 09-18-2023, 11:45 PM
RE: Vince's Corner Takeout - by Dav - 09-19-2023, 12:54 AM
RE: Vince's Corner Takeout - by bplus - 09-19-2023, 01:37 AM
RE: Vince's Corner Takeout - by GareBear - 09-19-2023, 03:56 PM
RE: Vince's Corner Takeout - by bplus - 09-19-2023, 04:47 PM
RE: Vince's Corner Takeout - by vince - 09-19-2023, 06:54 PM
RE: Vince's Corner Takeout - by bplus - 09-19-2023, 09:02 PM
RE: Vince's Corner Takeout - by vince - 01-13-2024, 07:15 PM
RE: Vince's Corner Takeout - by bplus - 01-13-2024, 07:59 PM
RE: Vince's Corner Takeout - by GareBear - 01-13-2024, 10:54 PM
RE: Vince's Corner Takeout - by vince - 02-16-2024, 04:01 AM
RE: Vince's Corner Takeout - by bplus - 02-16-2024, 02:27 PM
RE: Vince's Corner Takeout - by vince - 02-16-2024, 07:16 PM
RE: Vince's Corner Takeout - by Sprezzo - 02-17-2024, 03:04 AM
RE: Vince's Corner Takeout - by bplus - 02-17-2024, 02:44 PM
RE: Vince's Corner Takeout - by vince - 10-07-2024, 08:11 AM
RE: Vince's Corner Takeout - by bplus - 10-07-2024, 01:32 PM
RE: Vince's Corner Takeout - by vince - 10-11-2024, 12:05 AM
RE: Vince's Corner Takeout - by vince - 10-11-2024, 12:16 AM
RE: Vince's Corner Takeout - by PhilOfPerth - 10-11-2024, 03:51 AM
RE: Vince's Corner Takeout - by vince - 10-11-2024, 07:55 PM
RE: Vince's Corner Takeout - by PhilOfPerth - 10-11-2024, 11:30 PM
RE: Vince's Corner Takeout - by vince - 10-15-2024, 10:54 AM
RE: Vince's Corner Takeout - by bplus - 10-15-2024, 01:40 PM



Users browsing this thread: 13 Guest(s)