QBJS Complex Numbers

[vince]

some basic complex math routines to get started

Code: (Select All)

const sw = 1024
const sh = 768

dim shared pi
pi = 4*atn(1)

zoom = 140

screen _newimage(sw, sh, 32)

dim z(1), w(1), p(1), q(1), f(1), g(1)

for yy=0 to sh
for xx=0 to sw

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

i=0
select case i
case 0
    cnum p, x + 1, y - 1
    cnum q, x + 1, y + 1

    cnum z, 1, 0
    cmul w, z, p

    cnum z, w(0), w(1)
    cmul w, z, q

    cnum g, x - 1, y
    cnum z, w(0), w(1)
    cdiv w, z, g

    'pset (xx, yy), checker(w)
    pset (xx, yy), hue1(w)
case 1

    'cnum z, exp(1), 0
    cnum z, x, y
    'cnum w, x, y
    ccos w, z
    'cexp w, z, w
   
    pset (xx, yy), checker(w)
   
case 2

    n=10

    cnum g, 0, 0
    for j=0 to n-1
        'C: z(t)
        p(0) = 1.5*cos(j*2*pi/n)
        p(1) = 1.5*sin(j*2*pi/n)
       
        'f(z(t))
        csin w, p
        'cnum f, exp(1), 0
        'cexp w, f, w
       
        'f(z)/(z - z0)^(n + 1)
        cnum q, p(0) - x, p(1) - y
        cnum f, 2, 0
        cexp q, q, f
        cdiv w, w, q
       
        'dz/dt
        cnum q, -1.5*sin(j*2*pi/n), 1.5*cos(j*2*pi/n)
        cmul w, w, q
       
        if j=0 or j=n - 1 then
            g(0) = g(0) + 0.5*w(0)
            g(1) = g(1) + 0.5*w(1)
        else
            g(0) = g(0) + 0.5*w(0)
            g(1) = g(1) + 0.5*w(1)
        end if
    next
   
    'dt
    w(0) = g(0)*2*pi/n
    w(1) = g(1)*2*pi/n
   
    '1/(2 pi i)
    cnum q, 0, -1/(2*pi)
    cmul w, w, q
   
    'n!
    w(0) = 1*w(0)
    w(1) = 1*w(1)
   
    pset (xx, yy), checker(w)

case 3

    cnum z, x + 1, y
    cnum w, x - 1, y
   
    zz = checker(z)
    if zz = 0 then
    else
        pset (xx, yy), zz
    end if
   
    ww = checker(w)
    if ww = 0 then
    else
        pset (xx, yy), ww
    end if
   
    'pset (xx, yy), checker(z)' + checker(w)
    'pset (xx, yy), hue1(w)

end select
next
next

'n=100
'for a=0 to 2*pi step 2*pi/n
'    x = 1.5*cos(a)
'    y = 1.5*sin(a)
'    circle (x*zoom + sw/2, sh/2 - y*zoom), 3, _rgb(255,255,0)
'next
'sleep
system

function hue1( z() )

    m = sqr(z(0)*z(0) + z(1)*z(1))
    a = (pi + _atan2(z(1), z(0))) / (2*_pi)

    'dim rr, gg, bb
    'hue(v) ( .6 + .6 * cos( 2.*PI*(v) + vec3(0,-2.*PI/3.,2.*PI/3.)))
    rr =  0.5 - 0.5*sin(2*pi*a - pi/2)
    gg = (0.5 + 0.5*sin(2*pi*a*1.5 - pi/2)) * (a < 0.66)
    bb = (0.5 + 0.5*sin(2*pi*a*1.5 + pi/2)) * (a > 0.33)

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

    rr = rr - 0.0005*mm - 0.14*pp
    gg = gg - 0.0005*mm - 0.14*pp
    bb = bb - 0.0005*mm - 0.14*pp

    hue1 = _rgb(255*rr, 255*gg, 255*bb)
end function

function checker(z())
    if 0 then
        x = z(0)
        y = z(1)
    else
        x = _atan2(z(1), z(0))/(pi/8)
        y = sqr(z(0)*z(0) + z(1)*z(1))
    end if

    a = 2*(abs(x - int(x)))
    'b = 2*(abs(y - int(y)))

    $If Javascript Then
        c = a ^ b;
    $End If

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

sub cnum( w(), x, y )
    w(0) = x
    w(1) = y
end sub

sub cmul( w(), z1(), z2() )
    x1 = z1(0)
    y1 = z1(1)
    a1 = z2(0)
    b1 = z2(1)

    w(0) = x1*a1 - y1*b1
    w(1) = x1*b1 + y1*a1
end sub

sub cdiv( w(), z1(), z2() )
    x1 = z1(0)
    y1 = z1(1)
    a1 = z2(0)
    b1 = z2(1)

    d1 = a1*a1 + b1*b1
    w(0) = (x1*a1 + y1*b1)/d1
    w(1) = (y1*a1 - x1*b1)/d1
end sub

sub cexp( w(), z1(), z2() )
    x1 = z1(0)
    y1 = z1(1)
    a1 = z2(0)
    b1 = z2(1)
   
    lnz = x1*x1 + y1*y1
   
    if lnz = 0 then
        w(0) = 0
        w(1) = 0
    else
        lnz = 0.5*log(lnz)
        argz = _atan2(y1, x1)
        mz = exp(a1*lnz - b1*argz)
        az = a1*argz + b1*lnz
        w(0) = mz*cos(az)
        w(1) = mz*sin(az)
    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( w(), z1() )
    x1 = z1(0)
    y1 = z1(1)
    w(0) = sin(x1)*cosh(y1)
    w(1) = cos(x1)*sinh(y1)
end sub

sub ccos( w(), z1() )
    x1 = z1(0)
    y1 = z1(1)
    w(0) = cos(x1)*cosh(y1)
    w(1) =-sin(x1)*sinh(y1)
end sub