+- QB64 Phoenix Edition (https://staging.qb64phoenix.com)
+-- Forum: QB64 Rising (https://staging.qb64phoenix.com/forumdisplay.php?fid=1)
+--- Forum: QBJS, BAM, and Other BASICs (https://staging.qb64phoenix.com/forumdisplay.php?fid=50)
+--- Thread: QBJS Complex Numbers (/showthread.php?tid=1068)
QBJS Complex Numbers - vince - 11-07-2022
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