Skip to content

Instantly share code, notes, and snippets.

@loopspace loopspace/1aTabOrder
Created May 20, 2013

Embed
What would you like to do?
Library Maths Release v2.1 -A library of classes and functions for mathematical objects.
Library Maths Tab Order Version: 2.1
------------------------------
This file should not be included in the Codea project.
#ChangeLog
#Main
#Matrix
#Complex
#Quaternion
#Vec3
#Vector
#ccConfig
--[[
###########################################
##Codea Community Project Config Settings##
###########################################
##You can use # to comment out a line
##Use 1 for true and 0 for false
###########################################
# Add project info below #
#==========================================
ProjectName: Place Project Name Here
Version: Alpha 1.0
Comments:
Author:
##License Info: http://choosealicense.com
##Supported Licneses: MIT, GPL, Apache, NoLicense
License: MIT
#==========================================
###########################################
# Settings #
[Settings]=================================
##Codea Community Configuration settings
##Format: Setting state
Button 1
NotifyCCUpdate 1
ResetUserOption 0
AddHeaderInfo 1
[/Settings]================================
###########################################
# Screenshots #
[Screenshots]==============================
##Screenshots from your project.
##Format: url
##Example: http://www.dropbox.com/screenshot.jpg
[/Screenshots]=============================
###########################################
# Video #
[Video]====================================
##Link to a YouTube.com video.
##Format: url
##Example: http://www.youtube.com/videolink
[/Video]===================================
###########################################
# Dependencies #
[Dependencies]=============================
##Include the names of any dependencies here
##Format: Dependency
##Example: Codea Community
[/Dependencies]============================
############################################
# Tabs #
[Tabs]======================================
##Select which tabs are to be uploaded.
##Keyword 'not' excludes a tab or tabs. Keyword 'add' includes a tab or tabs.
##not * will exclude all tabs to allow for individual selection
##not *tab1 will look for any tab with tab1 on the end.
##not tab1* will look for any tab with tab1 at the beginning.
##Format: (add/not)(*)tab(*)
##Example: not Main --this will remove main.
##Example: not *tab1 --this will remove any tab with "tab1" on the end.
##Example: add Main --This will add Main back in.
[/Tabs]=====================================
#############################################
# Assets #
[Assets]=====================================
##Directory, path and url info for any assets besides the standard Codea assets.
##Format: Folder:sprite URL
##Example: Documents:sprite1 http://www.somewebsite.com/img.jpg
[/Assets]====================================
--]]
--[[
ChangeLog
=========
v2.1 Fixed bug in Quaternion "tomatrix" function.
v2.0 Split into separate libraries: "Library Maths" contains the
mathematical modules.
v1.0 Converted to use toadkick's cmodule for importing
and Briarfox's AutoGist for versionning.
It needs toadkick's code from
https://gist.github.com/apendley/5411561
tested with version 0.0.8
To use without AutoGist, comment out the lines in setup.
--]]
--[==[
-- Complex numbers class
-- Author: Andrew Stacey
-- Website: http://www.math.ntnu.no/~stacey/HowDidIDoThat/iPad/Codea.html
-- Licence: CC0 http://wiki.creativecommons.org/CC0
--[[
This is a class for handling complex numbers.
--]]
local Complex = class()
Complex.symbol = readLocalData("Complex Symbol","i")
local radeg = readLocalData("Complex Angle","rad")
if radeg == "rad" then
Complex.angle = 1
Complex.angsym = "π"
else
Complex.angle = 180
Complex.angsym = "°"
end
Complex.precision = readLocalData("Complex Precision",2)
--[[
A complex number can either be specified by giving the two
coordinates as real numbers or by giving a vector.
--]]
function Complex:init(...)
-- you can accept and set parameters here
if arg.n == 2 then
-- two numbers
self.z = vec2(arg[1],arg[2])
elseif arg.n == 1 then
-- vector
self.z = arg[1]
else
print("Incorrect number of arguments to Complex")
end
end
function Complex:clone(z)
return Complex(z:real(),z:imaginary())
end
--[[
Test if we are zero.
--]]
function Complex:is_zero()
-- are we the zero vector
if self.z ~= vec2(0,0) then
return false
end
return true
end
--[[
Test if we are real.
--]]
function Complex:is_real()
-- are we the zero vector
if self.z.y ~= 0 then
return false
end
return true
end
--[[
Test if the real part is zero.
--]]
function Complex:is_imaginary()
-- are we the zero vector
if self.z.x ~= 0 then
return false
end
return true
end
--[[
Test for equality.
--]]
function Complex:is_eq(q)
if self.z ~= q.z then
return false
end
return true
end
--[[
Defines the "==" shortcut.
--]]
function Complex:__eq(q)
return self:is_eq(q)
end
--[[
The inner product of two complex numbers.
Why did I program this?
--]]
function Complex:dot(q)
return self.z:dot(q.z)
end
--[[
Length of a complex number.
--]]
function Complex:len()
return self.z:len()
end
--[[
Often enough to know the length squared, which is quicker.
--]]
function Complex:lensq()
return self.z:lenSqr()
end
--[[
Distance between two complex numbers.
--]]
function Complex:dist(w)
return self.z:dist(w.z)
end
--[[
Often enough to know the distance squared, which is quicker.
--]]
function Complex:distSqr(w)
return self.z:distSqr(w.z)
end
--[[
Normalise a complex number to have length 1, if possible.
--]]
function Complex:normalise()
local l
if self:is_zero() then
print("Unable to normalise a zero-length complex number")
return false
end
l = 1/self:len()
return self:scale(l)
end
--[[
Scale the complex number.
--]]
function Complex:scale(l)
return Complex(l * self.z)
end
--[[
Add two complex numbers. Or add a real number to a complex one.
--]]
function Complex:add(q)
if type(q) == "number" then
return Complex(self.z + vec2(q,0))
else
return Complex(self.z + q.z)
end
end
--[[
q + p
--]]
function Complex:__add(q)
if type(self) == "number" then
return q:add(self)
else
return self:add(q)
end
end
--[[
Subtraction
--]]
function Complex:subtract(q)
if type(q) == "number" then
return Complex(self.z - vec2(q,0))
else
return Complex(self.z - q.z)
end
end
--[[
q - p
--]]
function Complex:__sub(q)
if type(self) == "number" then
return q:subtract(self):scale(-1)
else
return self:subtract(q)
end
end
--[[
Negation (-q)
--]]
function Complex:__unm()
return self:scale(-1)
end
--[[
Length (#q)
--]]
function Complex:__len()
return self:len()
end
--[[
Multiplication.
--]]
function Complex:multiply(q)
local a,b,c,d
a = self.z.x
b = self.z.y
c = q.z.x
d = q.z.y
return Complex(a*c - b*d,a*d + b*c)
end
--[[
q * p
--]]
function Complex:__mul(q)
if type(q) == "number" then
return self:scale(q)
elseif type(self) == "number" then
return q:scale(self)
elseif type(q) == "table" then
if q:is_a(Complex) then
return self:multiply(q)
end
end
end
--[[
Conjugation.
--]]
function Complex:conjugate()
return Complex(self.z.x, - self.z.y)
end
function Complex:co()
return self:conjugate()
end
--[[
Reciprocal: 1/q
--]]
function Complex:reciprocal()
if self:is_zero() then
print("Cannot reciprocate a zero complex number")
return false
end
local q = self:conjugate()
local l = self:lensq()
q = q:scale(1/l)
return q
end
--[[
Real powers.
--]]
function Complex:repower(n,k)
local r = self.z:len()
local t = -self.z:angleBetween(vec2(1,0))
k = k or 0
r = math.pow(r,n)
t = (t + k * 2 * math.pi) *n
return Complex:polar(r,t)
end
--[[
Complex powers.
--]]
function Complex:power(w,k)
if type(w) == "number" then
return self:repower(w,k)
end
if self:is_zero() then
print("Taking powers of 0 is somewhat dubious")
return false
end
local r = self.z:len()
local t = -self.z:angleBetween(vec2(1,0))
local u = w.z.x
local v = w.z.y
k = k or 0
local nr = math.pow(r,u) * math.exp(-v*t)
local nt = (t + k * 2 * math.pi) * u + math.log(r) * v
return Complex:polar(nr,nt)
end
--[[
q^n
This is overloaded so that a non-(complex) number exponent returns
the conjugate. This means that one can write things like q^"" to
get the conjugate of a complex number.
--]]
function Complex:__pow(n)
if type(n) == "number" then
return self:repower(n)
elseif type(n) == "table" and n:is_a(Complex) then
return self:power(n)
else
return self:conjugate()
end
end
--[[
Division: q/p
--]]
function Complex:__div(q)
if type(q) == "number" then
return self:scale(1/q)
elseif type(self) == "number" then
return q:scale(1/self):reciprocal()
elseif type(q) == "table" then
if q:is_a(Complex) then
return self:multiply(q:reciprocal())
end
end
end
--[[
Returns the real part.
--]]
function Complex:real()
return self.z.x
end
--[[
Returns the imaginary part.
--]]
function Complex:imaginary()
return self.z.y
end
--[[
Represents a complex number as a string.
--]]
Complex.precision = 2
function Complex:__concat(v)
if type(v) == "table"
and v:is_a(Complex) then
return self .. v:tostring()
else
return self:tostring() .. v
end
end
function Complex:__tostring()
return self:tostring()
end
function Complex:tostring()
local s
local x = math.floor(
self.z.x * 10^Complex.precision +.5
)/10^Complex.precision
local y = math.floor(
self.z.y * 10^Complex.precision +.5
)/10^Complex.precision
if x ~= 0 then
s = x
end
if y ~= 0 then
if s then
if y > 0 then
if y == 1 then
s = s .. " + " .. Complex.symbol
else
s = s .. " + " .. y .. Complex.symbol
end
else
if y == -1 then
s = s .. " - " .. Complex.symbol
else
s = s .. " - " .. (-y) .. Complex.symbol
end
end
else
if y == 1 then
s = Complex.symbol
elseif y == - 1 then
s = "-" .. Complex.symbol
else
s = y .. Complex.symbol
end
end
end
if s then
return s
else
return "0"
end
end
function Complex:topolarstring()
local t = math.floor(Complex.angle *
self:arg() * 10^Complex.precision/math.pi +.5
)/10^Complex.precision
local r = math.floor(
self:len() * 10^Complex.precision +.5
)/10^Complex.precision
local s = ""
--[[
-- this is exponential notation
if t == 0 then
s = r
else
if r ~= 1 then
s = r
end
s = s .. "e^(" .. t .."i)"
end
--]]
s = "(" .. r .. "," .. t .. Complex.angsym .. ")"
return s
end
function Complex:arg()
return -self.z:angleBetween(vec2(1,0))
end
function Complex:polar(r,t)
return Complex(r * math.cos(t), r * math.sin(t))
end
--[[
The unit complex number.
--]]
function Complex.unit()
return Complex(1,0)
end
function Complex.zero()
return Complex(0,0)
end
function Complex.i()
return Complex(0,-1)
end
--[[
Overload the maths functions
--]]
local __math = {}
for _,v in ipairs({
"abs",
"pow",
"sqrt",
"exp",
"log",
"sin",
"cos",
"sinh",
"cosh",
"tan",
"tanh"
}) do
__math[v] = math[v]
end
function math.abs(n)
if type(n) == "number" then return __math.abs(n) end
if type(n) == "table"
and n.is_a
and n:is_a(Complex)
then
return n:len()
end
print("Cannot take the length of " .. n)
end
function math.pow(n,e)
if type(n) == "number" then
if type(e) == "number" then
return __math.pow(n,e)
elseif type(e) == "table"
and e.is_a
and e:is_a(Complex)
then
local w = Complex(n,0)
return w:repower(e,0)
end
end
if type(n) == "table"
and n.is_a
and n:is_a(Complex)
then
return n:power(e,0)
end
print("Cannot take the power of " .. n .. " by " .. e)
end
function math.sqrt(n)
if type(n) == "number" then return __math.sqrt(n) end
if type(n) == "table"
and n.is_a
and n:is_a(Complex)
then
return n:repower(.5,0)
end
print("Cannot take the square root of " .. n)
end
function math.exp(n)
if type(n) == "number" then return __math.exp(n) end
if type(n) == "table"
and n.is_a
and n:is_a(Complex)
then
return Complex:polar(math.exp(n:real()),n:imaginary())
end
print("Cannot exponentiate " .. n)
end
--[[
cos(x+iy) = cos(x) cos(iy) - sin(x) sin(iy)
= cos(x) cosh(y) - i sin(x) sinh(y)
--]]
function math.cos(n)
if type(n) == "number" then return __math.cos(n) end
if type(n) == "table"
and n.is_a
and n:is_a(Complex)
then
local x = n:real()
local y = n:imaginary()
return Complex(__math.cos(x)*__math.cosh(y),-__math.sin(x)*__math.sinh(y))
end
print("Cannot take the cosine of " .. n)
end
--[[
sin(x+iy) = sin(x) cos(iy) + cos(x) sin(iy)
= sin(x) cosh(y) + i cos(x) sinh(y)
--]]
function math.sin(n)
if type(n) == "number" then return __math.sin(n) end
if type(n) == "table"
and n.is_a
and n:is_a(Complex)
then
local x = n:real()
local y = n:imaginary()
return Complex(__math.sin(x)*__math.cosh(y), __math.cos(x)*__math.sinh(y))
end
print("Cannot take the sine of " .. n)
end
--[[
cosh(x+iy) = cosh(x) cosh(iy) + sinh(x) sinh(iy)
= cosh(x) cos(y) + i sinh(x) sin(y)
--]]
function math.cosh(n)
if type(n) == "number" then return __math.cosh(n) end
if type(n) == "table"
and n.is_a
and n:is_a(Complex)
then
local x = n:real()
local y = n:imaginary()
return Complex(__math.cosh(x)*__math.cos(y), __math.sinh(x)*__math.sin(y))
end
print("Cannot take the hyperbolic cosine of " .. n)
end
--[[
sinh(x+iy) = sinh(x) cosh(iy) + cosh(x) sinh(iy)
= sinh(x) cos(y) + i cosh(x) sin(y)
--]]
function math.sinh(n)
if type(n) == "number" then return __math.sinh(n) end
if type(n) == "table"
and n.is_a
and n:is_a(Complex)
then
local x = n:real()
local y = n:imaginary()
return Complex(__math.sinh(x)*__math.cos(y), __math.cosh(x)*__math.sin(y))
end
print("Cannot take the hyperbolic sine of " .. n)
end
--[[
tan(x+iy) = (sin(x) cos(x) + i sinh(y) cosh(y))
/(cos^2(x) cosh^2(y) + sin^2(x) sinh^2(y))
--]]
function math.tan(n)
if type(n) == "number" then return __math.tan(n) end
if type(n) == "table"
and n.is_a
and n:is_a(Complex)
then
local x = n:real()
local y = n:imaginary()
local cx = __math.cos(x)
local sx = __math.sin(x)
local chy = __math.cosh(y)
local shy = __math.sinh(y)
local d = cx^2 * chy^2 + sx^2 * shy^2
if d == 0 then
return false
end
return Complex(sx*cx/d,shy*chy/d)
end
print("Cannot take the tangent of " .. n)
end
--[[
tanh(x+iy) = i tan(y - ix)
= (sin(x) cos(x) + i sinh(y) cosh(y))
/(cos^2(x) cosh^2(y) + sin^2(x) sinh^2(y))
= (sinh(x) cosh(x) + i sin(y) cos(y))
/(cos^2(y) cosh^2(x) + sin^2(y) sinh^2(x))
--]]
function math.tanh(n)
if type(n) == "number" then return __math.tanh(n) end
if type(n) == "table"
and n.is_a
and n:is_a(Complex)
then
local x = n:real()
local y = n:imaginary()
local cy = __math.cos(y)
local sy = __math.sin(y)
local chx = __math.cosh(x)
local shx = __math.sinh(x)
local d = cy^2 * chx^2 + sy^2 * shx^2
if d == 0 then
return false
end
return Complex(shx*chx/d,sy*cy/d)
end
print("Cannot take the hyperbolic tangent of " .. n)
end
--[[
log(r e^(i a)) = log(r) + i (a + 2 k pi)
--]]
function math.log(n,k)
if type(n) == "number" then
if k then
return Complex(__math.log(n), 2*k*math.pi)
else
return __math.log(n)
end
end
k = k or 0
if type(n) == "table"
and n.is_a
and n:is_a(Complex)
then
return Complex(__math.log(n:len()),n:arg() + 2*k*math.pi)
end
print("Cannot take the logarithm of " .. n)
end
if cmodule.loaded "TestSuite" then
testsuite.addTest({
name = "Complex",
setup = function()
local z = Complex(2,4)
local w = Complex(-1,1)
for k,v in ipairs({
{"z", z},
{"w", w},
{"Muliplication" , z*w},
{"ScalingLeft" , 2*z},
{"ScalingRight" , z*2},
{"Addition" , z+w},
{"Addition" , 3+w},
{"Addition" , w+3},
{"Subtraction" , z-w},
{"Subtraction" , 3-w},
{"Subtraction" , w-3},
{"Division" , z/w},
{"DivisionLeft" , 2/z},
{"DivisionRight" , z/2},
{"Reciprocal", 1/z},
{"Negation", -z},
{"Powers", z^2},
{"Roots", z^.5},
{"Complex Powers", z^w},
{"Conjugation", w^""},
{"Polar Form", z:topolarstring()},
{"Length", z:len()},
{"Overloading",""},
{"Absolute value", math.abs(z)},
{"Powers", math.pow(z,w)},
{"Square Root", math.sqrt(z)},
{"Cosine", math.cos(z)},
{"Sine", math.sin(z)},
{"Tangent", math.tan(z)},
{"Hyperbolic Cosine", math.cosh(z)},
{"Hyperbolic Sine", math.sinh(z)},
{"Hyperbolic Tangent", math.tanh(z)},
{"Logarithm", math.log(z,1)},
}) do
print(v[1] .. ": " .. v[2])
end
end,
draw = function()
end
})
end
return Complex
--]==]
--Project: Library Maths
--Version: 2.1
--Dependencies:
--Comments:
VERSION = 2.1
clearProjectData()
-- DEBUG = true
-- Use this function to perform your initial setup
function setup()
if AutoGist then
autogist = AutoGist("Library Maths","A library of classes and functions for mathematical objects.",VERSION)
autogist:backup(true)
end
if not cmodule then
openURL("http://loopspace.mathforge.org/discussion/36/my-codea-libraries")
print("You need to enable the module loading mechanism.")
print("See http://loopspace.mathforge.org/discussion/36/my-codea-libraries")
print("Touch the screen to exit the program.")
draw = function()
end
touched = function()
close()
end
return
end
--displayMode(FULLSCREEN_NO_BUTTONS)
cmodule "Library Maths"
cmodule.path("Library Base", "Library UI", "Library Utilities")
cimport "TestSuite"
Complex = cimport "Complex"
Quaternion = cimport "Quaternion"
--displayMode(FULLSCREEN)
local Touches = cimport "Touch"
local UI = cimport "UI"
cimport "Menu"
touches = Touches()
ui = UI(touches)
ui:systemmenu()
fps = {}
for k=1,20 do
table.insert(fps,1/60)
end
afps = 60
parameter.watch("math.floor(20/afps)")
sl = Quaternion.unit():make_slerp(Quaternion(1.00001,0,0,0))
print(sl(1))
qc = vec3(1,0,0):rotateTo(vec3(0,0,1))
print(qc)
print(qc:tomatrix())
print(vec3(1,0,0)^qc)
print(vec3(0,1,0)^qc)
print(vec3(0,0,1)^qc)
print(qc*Quaternion(0,1,0,0)*qc^"")
print(qc^""*Quaternion(0,1,0,0)*qc)
qq = Quaternion(qc)
end
function draw()
touches:draw()
table.remove(fps,1)
table.insert(fps,DeltaTime)
afps = 0
for k,v in ipairs(fps) do
afps = afps + v
end
background(75, 104, 90, 255)
ui:draw()
end
function touched(touch)
touches:addTouch(touch)
end
function orientationChanged(o)
end
function fullscreen()
end
function reset()
end
--[==[
-- Matrix
local Vector = cimport "Vector"
local _,Sentence = unpack(cimport "Font")
local Matrix = class()
function Matrix:init(t)
local r = 0
local c = 0
for k,v in ipairs(t) do
table.insert(self,Vector(v))
r = r + 1
end
c = self[1].size
for k,v in ipairs(self) do
if v.size ~= c then
return false
end
end
self.rows = r
self.cols = c
end
function Matrix:add(m)
if self.rows ~= m.rows then
return false
end
if self.cols ~= m.cols then
return false
end
local u = {}
for k,v in ipairs(self) do
table.insert(u,v + m[k])
end
return Matrix(u)
end
function Matrix:subtract(m)
if self.rows ~= m.rows then
return false
end
if self.cols ~= m.cols then
return false
end
local u = {}
for k,v in ipairs(self) do
table.insert(u,v - m[k])
end
return Matrix(u)
end
function Matrix:scale(l)
local u = {}
for k,v in ipairs(self) do
table.insert(u,l*v)
end
return Matrix(u)
end
function Matrix:multiplyRight(m)
if self.cols ~= m.rows then
return false
end
local u = {}
for k,v in ipairs(self) do
table.insert(u,v*m)
end
return Matrix(u)
end
function Matrix:multiplyLeft(m)
return m:multiplyRight(self)
end
function Matrix:transpose()
local nrows = {}
for i = 1,self.cols do
nrows[i] = {}
end
for k,v in ipairs(self) do
for i = 1,self.cols do
table.insert(nrows[i],v[i])
end
end
return Matrix(nrows)
end
--[[
Inline operators:
u + v
u - v
-u
u * v : scaling
u / v : scaling
u == v : equality
--]]
function Matrix:__add(v)
return self:add(v)
end
function Matrix:__sub(v)
return self:subtract(v)
end
function Matrix:__unm()
return self:scale(-1)
end
function Matrix:__mul(v)
if type(self) == "number" then
return v:scale(self)
elseif type(v) == "number" then
return self:scale(v)
elseif type(self) == "userdata" then
return Vector(self):applyMatrixRight(v):tovec()
elseif type(v) == "userdata" then
return Vector(v):applyMatrixLeft(self):tovec()
else
if self.is_a and self:is_a(Vector) then
return self:applyMatrixRight(v)
elseif v.is_a and v:is_a(Vector) then
return v:applyMatrixLeft(self)
else
return self:multiplyRight(v)
end
end
return false
end
function Matrix:__div(l)
if type(l) == "number" then
return self:scale(1/l)
else
return false
end
end
function Matrix:__eq(v)
return self:is_eq(v)
end
function Matrix:__concat(v)
if type(v) == "table"
and v:is_a(Matrix) then
return self .. v:tostring()
else
return self:tostring() .. v
end
end
function Matrix:tostring()
local u = {}
for k,v in ipairs(self) do
table.insert(u,v:tostring())
end
return "[" .. table.concat(u,";") .. "]"
end
function Matrix:__tostring()
return self:tostring()
end
function Matrix:toModelMatrix()
local m = {}
for j=1,4 do
for i=1,4 do
if self[i] and self[i][j] then
table.insert(m,self[i][j])
elseif i == j then
table.insert(m,1)
else
table.insert(m,0)
end
end
end
return matrix(unpack(m))
end
function Matrix:setParameters(t)
t = t or {}
for k,v in pairs({
"font",
"colour",
"ui",
"pos",
"rowSep",
"colSep",
"anchor"
}) do
if t[v] then
self[v] = t[v]
end
end
if t.touchHandler then
t.touchHandler:pushHandler(self)
end
end
function Matrix:prepare()
if not self.prepared then
local f = self.font
local t = {}
local r
local w = {}
local s
local col = self.colour or Colour.svg.White
for k,v in ipairs(self) do
r = {}
for l,u in ipairs(v) do
s = Sentence(f,math.floor(100*u+.5)/100)
s:prepare()
if not w[l] then
w[l] = s.width
else
w[l] = math.max(w[l],s.width)
end
s:setColour(col)
table.insert(r,s)
end
table.insert(t,r)
end
self.entries = t
self.widths = w
self.prepared = true
end
end
function Matrix:draw()
pushStyle()
resetStyle()
self:prepare()
local rs = self.rowSep or 0
local cs = self.colSep or 0
local x,y = WIDTH/2,HEIGHT/2
if self.pos then
x,y = self.pos()
end
local a = self.anchor or "centre"
local widths = self.widths
local w = cs * self.cols
for k,v in ipairs(widths) do
w = w + v
end
local lh = self.font:lineheight() + rs
local h = self.rows * lh
x,y = RectAnchorAt(x,y,w,h,a)
local yy = y + h + rs/2
local xx = x + cs/2
y = yy + rs/2 + self.font.descent
for k,v in ipairs(self.entries) do
x = xx
y = y - lh
for l,u in ipairs(v) do
u:draw(x,y)
x = x + widths[l] + cs
end
end
xx = xx - cs/2
yy = yy - rs/2
smooth()
strokeWidth(4)
local col = self.colour or Colour.svg.White
stroke(col)
lineCapMode(PROJECT)
line(xx-5,yy,xx,yy)
line(xx-5,yy,xx-5,yy-h)
line(xx-5,yy-h,xx,yy-h)
line(xx+w+5,yy,xx+w,yy)
line(xx+w+5,yy,xx+w+5,yy-h)
line(xx+w+5,yy-h,xx+w,yy-h)
popStyle()
end
function Matrix:isTouchedBy(touch)
local rs = self.rowSep or 0
local cs = self.colSep or 0
local x,y = WIDTH/2,HEIGHT/2
if self.pos then
x,y = self.pos()
end
local a = self.anchor or "centre"
local widths = self.widths
local w = cs * self.cols
for k,v in ipairs(widths) do
w = w + v
end
local lh = self.font:lineheight() + rs
local h = self.rows * lh
x,y = RectAnchorAt(x,y,w,h,a)
if touch.x < x then
return false
end
if touch.x > x + w then
return false
end
if touch.y < y then
return false
end
if touch.y > y + h then
return false
end
return true
end
function Matrix:processTouches(g)
if g.type.ended and g.type.tap then
local t = g.touchesArr[1]
local rs = self.rowSep or 0
local cs = self.colSep or 0
local x,y = WIDTH/2,HEIGHT/2
if self.pos then
x,y = self.pos()
end
local a = self.anchor or "centre"
local lh = self.font:lineheight() + rs
local h = self.rows * lh
local widths = self.widths
local w = cs * self.cols
for k,v in ipairs(widths) do
w = w + v
end
x,y = RectAnchorAt(x,y,w,h,a)
local r = math.floor((y + h - t.touch.y)/lh) + 1
w = x
local c
for k,v in ipairs(widths) do
w = w + v + cs
if t.touch.x < w then
c = k
break
end
end
self.ui:getNumberSpinner({
action = function(v)
self[r][c] = v
self.prepared = false
return true
end,
value = self[r][c]
})
end
if g.updated then
g:noted()
end
if g.type.ended then
g:reset()
end
end
return Matrix
--]==]
--[==[
-- Path object
local Path = class()
local PATH_STEP = 1
local PATH_RES = .05*.05
local PATH_RENDER = 0
local PATH_CREATE = 1
local PATH_POINTS = 2
function Path:init(t)
t = t or {}
self.style = t.style or {}
self.lastPoint = vec2(0,0)
self.elements = t.elements or {}
if t.touchHandler then
t.touchHandler:pushHandler(self)
end
self.touchpt = {}
self.m = mesh()
end
function Path:deletePoints()
self.elements = {}
self.lastPoint = vec2(0,0)
end
function Path:draw()
pushStyle()
resetStyle()
self:applyStyle(PATH_RENDER)
self.m:draw()
if self.edit then
self:applyStyle(PATH_POINTS)
local r = self.ptRadius
local lpt
pushStyle()
noStroke()
for k,v in ipairs(self.elements) do
ellipse(v[2].x,v[2].y,r)
if v[1] == "curve" then
popStyle()
line(lpt.x,lpt.y,v[2].x,v[2].y)
line(v[3].x,v[3].y,v[4].x,v[4].y)
pushStyle()
noStroke()
ellipse(v[3].x,v[3].y,r)
ellipse(v[4].x,v[4].y,r)
lpt = v[4]
else
lpt = v[2]
end
end
popStyle()
end
popStyle()
end
function Path:generate()
self:applyStyle(PATH_CREATE)
local ver = {}
local s,h
for k,v in ipairs(self.elements) do
ver,s,h = self:getPoints(v,ver,s,h)
end
if h then
ver,s = HobbyPoints(ver,s,h)
end
self.lastPoint = s
--debug:log({name = "path", message = "got called"})
print("got called")
self:makeMesh(ver)
end
function Path:applyStyle(t)
local s = self.style or {}
if t == PATH_CREATE then
self.linewidth = s.linewidth or 5
self.blur = s.blur or 1
self.smooth = s.smooth or false
self.drawColour = s.drawColour or Colour.svg.Black
elseif t == PATH_RENDER then
elseif t == PATH_POINTS then
self.ptRadius = s.pointRadius or 15
local l = s.pointLineWidth or 3
strokeWidth(l)
local l = s.pointLineCap or SQUARE
lineCapMode(l)
end
end
function Path:addElement(e)
table.insert(self.elements,e)
end
function Path:moveto(v)
self:addElement({"move",v})
self.lastPoint = v
end
function Path:lineto(v)
self:addElement({"line",v})
self.lastPoint = v
end
function Path:curveto(b,c,d)
self:addElement({"curve",b,c,d})
self.lastPoint = d
end
function Path:curvethrough(v)
self:addElement({"hobby",v})
self.lastPoint = v
end
function Path:getPoints(e,ver,s,h)
local t = e[1]
if h and t ~= "hobby" then
ver,s = HobbyPoints(ver,s,h)
end
if t == "move" then
table.insert(ver,{e[2]})
return ver, e[2]
elseif t == "line" then
local n = e[2] - s
n = n:rotate90()
n = n/n:len()
table.insert(ver,{s,n})
table.insert(ver,{e[2],n})
return ver,e[2]
elseif e[1] == "curve" then
return BezierPoints(ver,s,e[2],e[3],e[4])
elseif e[1] == "hobby" then
h = h or {}
table.insert(h,e[2])
return ver,s,h
else
return ver,s
end
end
function Path:makeMesh(pts)
local ver = {}
local col = {}
local l = self.linewidth/2
local b = self.blur + l
local s = false -- self.smooth
local c = self.drawColour
local ct = Colour.opacity(c,0)
local p,n
local m = 0
for k,v in ipairs(pts) do
if v[2] and n then
table.insert(ver,p + l*n)
table.insert(ver,p - l*n)
table.insert(ver,v[1] + l*v[2])
table.insert(ver,v[1] + l*v[2])
table.insert(ver,v[1] - l*v[2])
table.insert(ver,p - l*n)
for i=1,6 do
table.insert(col,c)
end
m = m + 6
if s then
table.insert(ver,p + l*n)
table.insert(ver,p + b*n)
table.insert(ver,v[1] + l*v[2])
table.insert(ver,v[1] + l*v[2])
table.insert(ver,v[1] + b*v[2])
table.insert(ver,p + b*n)
table.insert(ver,p - l*n)
table.insert(ver,p - b*n)
table.insert(ver,v[1] - l*v[2])
table.insert(ver,v[1] - l*v[2])
table.insert(ver,v[1] - b*v[2])
table.insert(ver,p - b*n)
for i=1,12 do
table.insert(col,ct)
end
end
end
p,n = unpack(v)
end
self.m.vertices = ver
self.m.colors = col
end
function Path:isTouchedBy(touch)
local p = vec2(touch.x,touch.y)
for k,v in ipairs(self.elements) do
if p:distSqr(v[2]) < 625 then
self.touchpt[touch.id] = v[2]
self.edit = true
return true
end
if v[1] == "curve" then
if p:distSqr(v[3]) < 625 then
self.touchpt[touch.id] = v[3]
self.edit = true
return true
end
if p:distSqr(v[4]) < 625 then
self.touchpt[touch.id] = v[4]
self.edit = true
return true
end
end
end
self.edit = false
return false
end
function Path:processTouches(g)
local regenerate = false
if g.updated then
for k,t in ipairs(g.touchesArr) do
if t.updated then
regenerate = true
self.touchpt[t.touch.id].x =
self.touchpt[t.touch.id].x
+ t.touch.deltaX
self.touchpt[t.touch.id].y =
self.touchpt[t.touch.id].y
+ t.touch.deltaY
end
if t.touch.state == ENDED then
self.touchpt[t.touch.id] = nil
end
end
g:noted()
if regenerate then
self:generate()
end
end
if g.type.ended then
g:reset()
end
end
function Path:setLineWidth(l)
self.style.linewidth = l
end
function Path:setBlur(l)
self.style.blur = l
end
function Path:setSmooth(l)
self.style.smooth = l
end
function Path:setColour(c)
self.style.drawColour = c
end
function BezierPoints(pts,a,b,c,d)
local t = 0
local r = PATH_RES
local s = PATH_STEP
local tpt = a
table.insert(pts,{tpt,Bezierunml(0,a,b,c,d)})
local dis
local p
while t < 1 do
dis = 0
while dis < r do
t = t + s
p = Bezierpt(t,a,b,c,d)
dis = p:distSqr(tpt)
end
if t > 1 then
t = 1
p = d
end
table.insert(pts,{p,Bezierunml(t,a,b,c,d)})
tpt = p
end
return pts,d
end
function Bezierpt(t,a,b,c,d)
local s = 1 - t
return s^3 * a
+ 3*s*s*t * b
+ 3*s*t*t * c
+ t^3 * d
end
function Beziertgt(t,a,b,c,d)
local s = 1 - t
return 3*s^2 * (b - a)
+ 6*s*t * (c - b)
+ 3*t^2 * (d - c)
end
function Beziernml(t,a,b,c,d)
return Beziertgt(t,a,b,c,d):rotate90()
end
function Bezierunml(t,a,b,c,d)
local pt = Beziernml(t,a,b,c,d)
local l = pt:len()
if l == 0 then
return vec2(0,0)
else
return pt/l
end
end
function Bezierutgt(t,a,b,c,d)
local pt = Beziertgt(t,a,b,c,d)
local l = pt:len()
if l == 0 then
return vec2(0,0)
else
return pt/l
end
end
function QuickHobby(a,b,c,tha)
if type(a) == "table" then
tha = b
a,b,c = unpack(a)
end
local da = a:dist(b)
local db = b:dist(c)
local wa = vec2(1,0):angleBetween(b-a)
local wb = vec2(1,0):angleBetween(c-b)
local psi = wb - wa
if psi > math.pi then
psi = psi - 2*math.pi
end
if psi <= -math.pi then
psi = psi + 2*math.pi
end
local thb,phb,phc
if tha then
thb = -(2*psi + tha) * db / (2*db + da)
phb = - psi - thb
phc = thb
else
thb = - psi * db / (da + db)
tha = - psi - thb
phb = tha
phc = thb
end
local alpha = math.sqrt(2) * (math.sin(tha) - 1/16 * math.sin(phb)) * (math.sin(phb) - 1/16 * math.sin(tha)) * (math.cos(tha) - math.cos(phb))
local rho = (2 + alpha)/(1 + (1 - (3 - math.sqrt(5))/2) * math.cos(tha) + (3 - math.sqrt(5))/2 * math.cos(phb))
local sigma = (2 - alpha)/(1 + (1 - (3 - math.sqrt(5))/2) * math.cos(phb) + (3 - math.sqrt(5))/2 * math.cos(tha))
local ctrlaa = a + da * rho * vec2(math.cos(tha + wa), math.sin(tha + wa))/3
local ctrlab = b - da * sigma * vec2(math.cos(wa - phb), math.sin(wa - phb))/3
alpha = math.sqrt(2) * (math.sin(thb) - 1/16 * math.sin(phc)) * (math.sin(phc) - 1/16 * math.sin(thb)) * (math.cos(thb) - math.cos(phc))
rho = (2 + alpha)/(1 + (1 - (3 - math.sqrt(5))/2) * math.cos(thb) + (3 - math.sqrt(5))/2 * math.cos(phc))
sigma = (2 - alpha)/(1 + (1 - (3 - math.sqrt(5))/2) * math.cos(phc) + (3 - math.sqrt(5))/2 * math.cos(thb))
local ctrlba = b + db * rho * vec2(math.cos(thb + wb), math.sin(thb + wb))/3
local ctrlbb = c - db * sigma * vec2(math.cos(wb - phc), math.sin(wb - phc))/3
return {a,ctrlaa,ctrlab,b},{b,ctrlba,ctrlbb,c},thb
end
function QHobbyGenerator(a,b)
local th
return function(c)
local p,q
p,q,th = QuickHobby(a,b,c,th)
a = b
b = c
return p,q
end
end
function HobbyPoints(ver,s,pts,extra)
if #pts == 1 then
local v = pts[1]
if type(v) == "table" then
v = v[1]
end
local n = v - s
n = n:rotate90()
n = n/n:len()
table.insert(ver,{s,n})
table.insert(ver,{v,n})
return ver,v
end
table.insert(pts,1,s)
local bcurves = Hobby(pts,extra)
for _,v in ipairs(bcurves) do
ver,s = BezierPoints(ver,unpack(v))
end
return ver,s
end
function Hobby(pts,extra)
local z = {}
local d = {}
local omega = {}
local psi = {}
local it = {}
local ot = {}
local n = 0
local A = {}
local B = {}
local C = {}
local D = {}
local icurl = 1
local ocurl = 1
local theta = {}
local phi = {}
local rho ={}
local sigma = {}
local a = math.sqrt(2)
local b = 1/16
local c = (3 - math.sqrt(5))/2
local cpta = {}
local cptb = {}
local dten = 1
local curves = {}
ot[0] = 1
if extra then
dten = extra.tension or 1
icurl = extra.inCurl or 1
ocurl = extra.outCurl or 1
ot[0] = extra.initialTension or 1
end
local ang,i,o,v
z[0] = table.remove(pts,1)
for _,t in ipairs(pts) do
if type(t) == "table" then
v = t[1]
i = t[2] or dten
o = t[3] or dten
else
v = t
i = dten
o = dten
end
d[n] = z[n]:dist(v)
omega[n] = vec2(1,0):angleBetween(v-z[n])
if n > 0 then
ang = omega[n] - omega[n-1]
if ang > math.pi then
ang = ang - 2*math.pi
end
if ang < -math.pi then
ang = ang + 2*math.pi
end
psi[n] = ang
end
n = n + 1
z[n] = v
it[n] = i
ot[n] = o
end
if extra then
theta[0] = extra.inAngle
phi[n] = extra.outAngle
end
for k=1,n-1 do
A[k] = d[k] * it[k+1] * it[k]^2
end
if theta[0] then
B[0] = 1
else
B[0] = ot[0]^3 * (3 * it[1] - 1) + icurl * it[1]^3
end
for k = 1,n-2 do
B[k] = d[k] * it[k+1] * it[k]^2 * (3 * ot[k-1] - 1) + d[k-1] * ot[k-1] * ot[k]^2 * (3 * it[k+1] - 1)
end
B[n-1] = d[n-1] * it[n] * it[n-1]^2 * (3 * ot[n-2] - 1) + d[n-2] * ot[n-2] * ot[n-1]^2 * (3 * it[n] - 1)
if not phi[n] then
B[n-1] = B[n-1] - d[n-2] * ot[n-2] * ot[n-1]^2 * (it[n]^3 + ocurl * ot[n-1]^3 * (3 * it[n] - 1)) / (it[n]^3 * (3 * ot[n-1] - 1) + ocurl * ot[n-1]^3)
end
if theta[0] then
C[0] = 0
else
C[0] = ot[0]^3 + icurl * it[1]^3 * (3 * ot[0] - 1)
end
for i=1,n do
C[i] = d[i-1] * ot[i-1] * ot[i]^2
end
if theta[0] then
D[0] = theta[0]
else
D[0] = - (ot[0]^3 + icurl * it[1]^3 * (3 * ot[0] - 1)) * psi[1]
end
for i=1,n-2 do
D[i] = - d[i] * it[i+1] * it[i]^2 * (3 * ot[i-1] - 1) * psi[i] - d[i-1] * ot[i-1] * ot[i]^2 * psi[i+1]
end
D[n-1] = - d[n-1] * it[n] * it[n-1]^2 * (3 * ot[n-2] - 1) * psi[n-1]
if phi[n] then
D[n-1] = D[n-1] - d[n-2] * ot[n-2] * ot[n-1]^2 * phi[n]
end
for i=1,n-1 do
B[i] = B[i-1] * B[i] - A[i] * C[i-1]
C[i] = B[i-1] * C[i]
D[i] = B[i-1] * D[i] - A[i] * D[i-1]
end
theta[n-1] = D[n-1]/B[n-1]
for i=n-2,1,-1 do
theta[i] = (D[i] - C[i] * theta[i+1])/B[i]
end
for i=1,n-1 do
phi[i] = -psi[i] - theta[i]
end
if not theta[0] then
theta[0] = (ot[0]^3 + icurl * it[1]^3 * (3 * ot[0] - 1)) / (ot[0]^3 * (3 * it[1] - 1) + icurl * it[1]^3) * phi[1]
end
if not phi[n] then
phi[n] = (it[n]^3 + ocurl * it[n-1]^3 * (3 * it[n] - 1)) / (it[n]^3 * (3 * ot[n-1] - 1) + ocurl * ot[n-1]^3) * theta[n-1]
end
local alpha
for i = 0,n-1 do
alpha = a * (math.sin(theta[i]) - b * math.sin(phi[i+1])) * (math.sin(phi[i+1]) - b * math.sin(theta[i])) * (math.cos(theta[i]) - math.cos(phi[i+1]))
rho[i] = (2 + alpha) / (1 + (1 - c) * math.cos(theta[i]) + c * math.cos(phi[i+1]))
sigma[i+1] = (2 - alpha) / (1 + (1 - c) * math.cos(phi[i+1]) + c * math.cos(theta[i]))
end
for i = 0,n-1 do
table.insert(curves,{
z[i],
z[i] + d[i]*rho[i] * vec2(math.cos(theta[i] + omega[i]), math.sin(theta[i] + omega[i]))/3,
z[i+1] - d[i] * sigma[i+1] * vec2(math.cos(omega[i] - phi[i+1]), math.sin(omega[i] - phi[i+1]))/3,
z[i+1]
})
end
return curves
end
cmodule.gexport {
BezierPoints = BezierPoints,
Bezierpt = Bezierpt,
Beziertgt = Beziertgt,
Beziernml = Beziernml,
Bezierunml = Bezierunml,
Bezierutgt = Bezierutgt,
QuickHobby = QuickHobby,
HobbyPoints = HobbyPoints,
QHobbyGenerator = QHobbyGenerator,
Hobby = Hobby
}
return Path
--]==]
--[==[
-- Quaternions
-- Author: Andrew Stacey
-- Website: http://www.math.ntnu.no/~stacey/HowDidIDoThat/iPad/Codea.html
-- Licence: CC0 http://wiki.creativecommons.org/CC0
--[[
This is a class for handling quaternion numbers. It was originally
designed as a way of encoding rotations of 3 dimensional space.
--]]
local Quaternion = class()
--[[
A quaternion can either be specified by giving the four coordinates as
real numbers or by giving the scalar part and the vector part.
--]]
local error = error or print
function Quaternion:init(...)
-- you can accept and set parameters here
if arg.n == 4 then
-- four numbers
self.q = vec4(arg[1],arg[2],arg[3],arg[4])
elseif arg.n == 2 then
-- real number plus vector
self.q = vec4(arg[1],arg[2].x,arg[2].y,arg[2].z)
elseif arg.n == 1 then
if type(arg[1]) == "userdata" then
self.q = vec4(arg[1].x,arg[1].y,arg[1].z,arg[1].w)
elseif type(arg[1]) == "table" and arg[1].is_a and arg[1]:is_a(Quaternion) then
self.q = vec4(arg[1].q.x,arg[1].q.y,arg[1].q.z,arg[1].q.w)
else
error("Incorrect type of argument to Quaternion")
end
else
error("Incorrect number of arguments to Quaternion")
end
end
--[[
Test if we are zero.
--]]
function Quaternion:is_zero()
return self.q == vec4(0,0,0,0)
end
--[[
Test if we are real.
--]]
function Quaternion:is_real()
-- are we the zero vector
if self.q.y ~= 0 or self.q.z ~= 0 or self.q.w ~= 0 then
return false
end
return true
end
--[[
Test if the real part is zero.
--]]
function Quaternion:is_imaginary()
-- are we the zero vector
if self.q.x ~= 0 then
return false
end
return true
end
--[[
Test for equality.
--]]
function Quaternion:is_eq(q)
return self.q == q.q
end
--[[
Defines the "==" shortcut.
--]]
function Quaternion:__eq(q)
return self:is_eq(q)
end
--[[
The inner product of two quaternions.
--]]
function Quaternion:dot(q)
return self.q:dot(q.q)
end
--[[
Makes "q .. p" return the inner product.
Probably a bad choice and likely to be removed in future versions.
function Quaternion:__concat(q)
return self:dot(q)
end
--]]
--[[
Length of a quaternion.
--]]
function Quaternion:len()
return self.q:len()
end
--[[
Often enough to know the length squared, which is quicker.
--]]
function Quaternion:lenSqr()
return self.q:lenSqr()
end
--[[
Distance between two quaternions.
--]]
function Quaternion:dist(q)
return self.q:dist(q.q)
end
--[[
Often enough to know the distance squared, which is quicker.
--]]
function Quaternion:distSqr(q)
return self.q:distSqr(q.q)
end
--[[
Normalise a quaternion to have length 1, if possible.
--]]
function Quaternion:normalise()
return Quaternion(self.q:normalize())
end
function Quaternion:normalize()
return Quaternion(self.q:normalize())
end
--[[
Scale the quaternion.
--]]
function Quaternion:scale(l)
return Quaternion(self.q * l)
end
--[[
Add two quaternions.
--]]
function Quaternion:add(q)
if type(q) == "number" then
return Quaternion(self.q.x + q, self.q.y, self.q.z, self.q.w)
else
return Quaternion(self.q + q.q)
end
end
--[[
q + p
--]]
function Quaternion:__add(q)
if type(q) == "number" then
return Quaternion(self.q.x + q, self.q.y, self.q.z, self.q.w)
elseif type(self) == "number" then
return Quaternion(self + q.q.x, q.q.y, q.q.z, q.q.w)
else
return Quaternion(self.q + q.q)
end
end
--[[
Subtraction
--]]
function Quaternion:subtract(q)
return Quaternion(self.q - q.q)
end
--[[
q - p
--]]
function Quaternion:__sub(q)
if type(q) == "number" then
return Quaternion(self.q.x - q, self.q.y, self.q.z, self.q.w)
elseif type(self) == "number" then
return Quaternion(self - q.q.x, - q.q.y, - q.q.z, - q.q.w)
else
return Quaternion(self.q - q.q)
end
end
--[[
Negation (-q)
--]]
function Quaternion:__unm()
return Quaternion(-self.q)
end
--[[
Length (#q)
--]]
function Quaternion:__len()
return self:len()
end
--[[
Multiply the current quaternion on the right.
Corresponds to composition of rotations.
--]]
function Quaternion:multiplyRight(q)
local a,b,c,d
a = self.q.x * q.q.x - self.q.y * q.q.y - self.q.z * q.q.z - self.q.w * q.q.w
b = self.q.x * q.q.y + self.q.y * q.q.x + self.q.z * q.q.w - self.q.w * q.q.z
c = self.q.x * q.q.z - self.q.y * q.q.w + self.q.z * q.q.x + self.q.w * q.q.y
d = self.q.x * q.q.w + self.q.y * q.q.z - self.q.z * q.q.y + self.q.w * q.q.x
return Quaternion(a,b,c,d)
end
--[[
q * p
--]]
function Quaternion:__mul(q)
if type(q) == "number" then
return self:scale(q)
elseif type(self) == "number" then
return q:scale(self)
elseif type(q) == "table" and q.is_a and q:is_a(Quaternion) then
return self:multiplyRight(q)
end
end
--[[
Multiply the current quaternion on the left.
Corresponds to composition of rotations.
--]]
function Quaternion:multiplyLeft(q)
return q:multiplyRight(self)
end
--[[
Conjugation (corresponds to inverting a rotation).
--]]
function Quaternion:conjugate()
return Quaternion(self.q.x, - self.q.y, - self.q.z, - self.q.w)
end
function Quaternion:co()
return self:conjugate()
end
--[[
Reciprocal: 1/q
--]]
function Quaternion:reciprocal()
if self:is_zero() then
error("Cannot reciprocate a zero quaternion")
return false
end
local q = self:conjugate()
local l = self:lenSqr()
q = q:scale(1/l)
return q
end
--[[
Integral powers.
--]]
function Quaternion:power(n)
if n ~= math.floor(n) then
error("Only able to do integer powers")
return false
end
if n == 0 then
return Quaternion(1,0,0,0)
elseif n > 0 then
return self:multiplyRight(self:power(n-1))
elseif n < 0 then
return self:reciprocal():power(-n)
end
end
--[[
q^n
This is overloaded so that a non-number exponent returns the
conjugate. This means that one can write things like q^* or q^"" to
get the conjugate of a quaternion.
--]]
function Quaternion:__pow(n)
if type(n) == "number" then
return self:power(n)
elseif type(n) == "table" and n.is_a and n:is_a(Quaternion) then
return self:multiplyLeft(n):multiplyRight(n:reciprocal())
else
return self:conjugate()
end
end
--[[
Division: q/p
--]]
function Quaternion:__div(q)
if type(q) == "number" then
return self:scale(1/q)
elseif type(self) == "number" then
return q:reciprocal():scale(self)
elseif type(q) == "table" and q.is_a and q:is_a(Quaternion) then
return self:multiplyRight(q:reciprocal())
end
end
--[[
Interpolation functions, we assume the input to be already normalised
for speed. If you cannot guarantee this, renormalise the input first.
The constructor functions do do the renormalisation.
--]]
--[[
Linear interpolation, renormalised
--]]
function Quaternion:lerp(q,t)
if not t then
return Quaternion.unit():lerp(self,q)
end
local v
if self.q == -q.q then
-- antipodal points, need a midpoint
v = vec4(self.q.y,-self.q.x,self.q.w,-self.q.z)
v = (1 - 2*t)*self.q + (1-math.abs(2*t-1))*v
else
v = (1-t)*self.q + t*q.q
end
return Quaternion(v:normalize())
end
--[[
Spherical interpolation
--]]
function Quaternion:slerp(q,t)
if not t then
return Quaternion.unit():slerp(self,q)
end
local v
if self.q == -q.q then
-- antipodal points, need a midpoint
v = vec4(self.q.y,-self.q.x,self.q.w,-self.q.z)
t = 2*t
elseif self.q == q.q then
return Quaternion(self)
else
v = q.q
end
local ca = self.q:dot(v)
local sa = math.sqrt(1 - math.pow(ca,2))
if sa == 0 then
return Quaternion(self)
end
local a = math.acos(ca)
sa = math.sin(a*t)/sa
v = (math.cos(a*t) - ca*sa)*self.q+ sa*v
return Quaternion(v)
end
--[[
Constructor for normalised linear interpolation.
--]]
function Quaternion:make_lerp(q)
if not q then
return Quaternion.unit():make_lerp(self)
end
local v,w
w = self.q:normalize()
if self.q == -q.q then
-- antipodal points, need a midpoint
v = vec4(w.y,-w.x,w.w,-w.z)
return function(t)
local u = (1 - 2*t)*w + (1-math.abs(2*t-1))*v
return Quaternion(u:normalize())
end
else
v = q.q:normalize()
return function(t)
local u = (1-t)*w + t*v
return Quaternion(u:normalize())
end
end
end
--[[
Spherical interpolation
--]]
function Quaternion:make_slerp(q)
if not q then
return Quaternion.unit():make_slerp(self)
end
local v,f,u
if self.q == -q.q then
-- antipodal points, need a midpoint
v = vec4(self.q.y,-self.q.x,self.q.w,-self.q.z)
f = 2
elseif self.q == q.q then
return function(t)
return Quaternion(self)
end
else
v = q.q
f = 1
end
v = v:normalize()
u = self.q:normalize()
local ca = u:dot(v)
local sa = math.sqrt(1 - math.pow(ca,2))
if sa == 0 then
return function(t)
return Quaternion(self)
end
end
local a = math.acos(ca)
v = (v - ca*self.q)/sa
return function(t)
local u = math.cos(a*f*t)*self.q + math.sin(a*f*t)*v
return Quaternion(u)
end
end
--[[
Returns the real part.
--]]
function Quaternion:real()
return self.q.x
end
--[[
Returns the vector (imaginary) part as a Vec3 object.
--]]
function Quaternion:vector()
return vec3(self.q.y, self.q.z, self.q.w)
end
--[[
Represents a quaternion as a string.
--]]
function Quaternion:tostring()
local s
local im ={{self.q.y,"i"},{self.q.z,"j"},{self.q.w,"k"}}
if self.q.x ~= 0 then
s = string.format("%.3f",self.q.x)
end
for k,v in pairs(im) do
if v[1] ~= 0 then
if s then
if v[1] > 0 then
if v[1] == 1 then
s = s .. " + " .. v[2]
else
s = s .. " + " ..
string.format("%.3f",v[1]) .. v[2]
end
else
if v[1] == -1 then
s = s .. " - " .. v[2]
else
s = s .. " - " ..
string.format("%.3f",-v[1]) .. v[2]
end
end
else
if v[1] == 1 then
s = v[2]
elseif v[1] == - 1 then
s = "-" .. v[2]
else
s = string.format("%.3f",v[1]) .. v[2]
end
end
end
end
if s then
return s
else
return "0"
end
end
function Quaternion:__tostring()
return self:tostring()
end
function Quaternion:__concat(s)
if type(s) == "string" then
return self:tostring() .. s
elseif type(s) == "table" and s.is_a and s:is_a(Quaternion) then
return self .. s:tostring()
end
end
function Quaternion:tomatrix()
local a,b,c,d = self.q.x,-self.q.y,-self.q.z,-self.q.w
local ab = 2*a*b
local ac = 2*a*c
local ad = 2*a*d
local bb = 2*b*b
local bc = 2*b*c
local bd = 2*b*d
local cc = 2*c*c
local cd = 2*c*d
local dd = 2*d*d
return matrix(
1 - cc - dd,
bc - ad,
ac + bd,
0,
bc + ad,
1 - bb - dd,
cd - ab,
0,
bd - ac,
cd + ab,
1 - bb - cc,
0,
0,0,0,1
)
end
--[[
(Not a class function)
Returns a quaternion corresponding to the current gravitational vector
so that after applying the corresponding rotation, the y-axis points
in the gravitational direction and the x-axis is in the plane of the
iPad screen.
When we have access to the compass, the x-axis behaviour might change.
--]]
--[[
function Quaternion.Gravity()
local gxy, gy, gygxy, a, b, c, d
if Gravity.x == 0 and Gravity.y == 0 then
return Quaternion(1,0,0,0)
else
gy = - Gravity.y
gxy = math.sqrt(math.pow(Gravity.x,2) + math.pow(Gravity.y,2))
gygxy = gy/gxy
a = math.sqrt(1 + gxy - gygxy - gy)/2
b = math.sqrt(1 - gxy - gygxy + gy)/2
c = math.sqrt(1 - gxy + gygxy - gy)/2
d = math.sqrt(1 + gxy + gygxy + gy)/2
if Gravity.y > 0 then
a = a
b = b
end
if Gravity.z < 0 then
b = - b
c = - c
end
if Gravity.x > 0 then
c = - c
d = - d
end
return Quaternion(a,b,c,d)
end
end
--]]
function Quaternion:Gravity()
local qg,qx
if not self then
qx = Quaternion.ReferenceFrame()
else
qx = self
end
local y = vec3(0,-1,0)^qx
qg = y:rotateTo(Gravity)
return qg*qx
end
local frame = {
vec3(1,0,0),
vec3(0,1,0),
vec3(0,0,1)
}
local qzyx = Quaternion(1,0,0,0)
--[[
Needs to be run once every frame!
--]]
function Quaternion:updateReferenceFrame()
local x,y,z
if CurrentOrientation == PORTRAIT then
x,y,z = -frame[1],-frame[2],-frame[3]
elseif CurrentOrientation == PORTRAIT_UPSIDE_DOWN then
x,y,z = frame[1],frame[2],-frame[3]
elseif CurrentOrientation == LANDSCAPE_LEFT then
x,y,z = frame[2],-frame[1],-frame[3]
elseif CurrentOrientation == LANDSCAPE_RIGHT then
x,y,z = -frame[2],frame[1],-frame[3]
end
local qz = Quaternion.Rotation(RotationRate.z*DeltaTime,z.x,z.y,z.z)
local qy = Quaternion.Rotation(RotationRate.y*DeltaTime,y.x,y.y,y.z)
local qx = Quaternion.Rotation(RotationRate.x*DeltaTime,x.x,x.y,x.z)
if self then
local q = qz * qy * qx * self
self.q = q.q
else
qzyx = qz * qy * qx * qzyx
return qzyx
end
end
function Quaternion:ReferenceFrame()
return self or qzyx
end
--[[
Converts a rotation to a quaternion. The first argument is the angle
to rotate, the rest must specify an axis, either as a Vec3 object or
as three numbers.
--]]
function Quaternion.Rotation(a,...)
local q,c,s
q = Quaternion(0,...)
q = q:normalise()
c = math.cos(a/2)
s = math.sin(a/2)
q = q:scale(s)
q = q:add(c)
return q
end
--[[
The unit quaternion.
--]]
function Quaternion.unit()
return Quaternion(1,0,0,0)
end
--[[
Extensions to vec3 type.
--]]
do
local mt = getmetatable(vec3())
--[[
Promote to a quaternion with 0 real part.
--]]
mt["toQuaternion"] = function (self)
return Quaternion(0,self.x,self.y,self.z)
end
--[[
Apply a quaternion as a rotation (assumes unit quaternion for speed).
--]]
mt["applyQuaternion"] = function (self,q)
local x = self:toQuaternion()
x = q:multiplyRight(x)
x = x:multiplyRight(q:conjugate())
return x:vector()
end
mt["__pow"] = function (self,q)
if type(q) == "table" then
if q:is_a(Quaternion) then
return self:applyQuaternion(q)
end
end
return false
end
mt["rotateTo"] = function (self,v)
if v:lenSqr() == 0 or self:lenSqr() == 0 then
return Quaternion(1,0,0,0)
end
local u = self:normalize()
v = u + v:normalize()
if v:lenSqr() == 0 then
-- Opposite vectors, no canonical direction
local a,b,c = math.abs(u.x), math.abs(u.y), math.abs(u.z)
if a < b and a < c then
v = vec3(0,-u.z,u.y)
elseif b < c then
v = vec3(u.z,0,-u.x)
else
v = vec3(u.y,-u.x,0)
end
end
v = v:normalize()
return Quaternion(u:dot(v),u:cross(v))
end
end
if cmodule.loaded "TestSuite" then
testsuite.addTest({
name = "Quaternion",
setup = function()
local q = Quaternion(.3,.4,.5,.6)
local qq = Quaternion(.2,vec3(.4,-.5,.1))
for k,v in ipairs({
{"Quaternion", q},
{"Sum", q + qq},
{"Sum (number)", q + 5},
{"Sum (number)", 5 + q},
{"Subtract", q - qq},
{"Subtract (number)", q - 5},
{"Subtract (number)", 5 - q},
{"Product", q * qq},
{"Scale (right)", q * 5},
{"Scale (left)", 5 * q},
{"Divison", q / qq},
{"by number", q / 5},
{"of number", 5 / q},
{"Length", q:len()},
{"Length Squared", q:lenSqr()},
}) do
print(v[1] .. ": " .. v[2])
end
end,
draw = function()
end
})
end
return Quaternion
--]==]
--[==[
-- Vec3 class
-- Author: Andrew Stacey
-- Website: http://www.math.ntnu.no/~stacey/HowDidIDoThat/iPad/Codea.html
-- Licence: CC0 http://wiki.creativecommons.org/CC0
--[[
The "Vec3" class is for handling 3 dimensional vectors and defining a
variety of methods on them.
This is largely superceded by the vec3 userdata. There are a few extra
methods, but these should really be added to vec3.
--]]
local Vec3 = class()
--[[
A 3-vector is three numbers.
--]]
function Vec3:init(x,y,z)
self.x = x
self.y = y
self.z = z
end
--[[
Test for zero vector.
--]]
function Vec3:is_zero()
if self.x ~= 0 or self.y ~= 0 or self.z ~= 0 then
return false
end
return true
end
--[[
Test for equality.
--]]
function Vec3:is_eq(v)
if self.x ~= v.x or self.y ~= v.y or self.z ~= v.z then
return false
end
return true
end
--[[
Inner product.
--]]
function Vec3:dot(v)
return self.x * v.x + self.y * v.y + self.z * v.z
end
--[[
Cross product.
--]]
function Vec3:cross(v)
local x,y,z
x = self.y * v.z - self.z * v.y
y = self.z * v.x - self.x * v.z
z = self.x * v.y - self.y * v.x
return Vec3(x,y,z)
end
--[[
Apply a given matrix (which is specified as a triple of vectors).
--]]
function Vec3:applyMatrix(a,b,c)
local u,v,w
u = a:scale(self.x)
v = b:scale(self.y)
w = c:scale(self.z)
u = u:add(v)
u = u:add(w)
return u
end
--[[
Length of the vector
--]]
function Vec3:len()
return math.sqrt(math.pow(self.x,2) + math.pow(self.y,2) + math.pow(self.z,2))
end
--[[
Squared length of the vector.
--]]
function Vec3:lenSqr()
return math.pow(self.x,2) + math.pow(self.y,2) + math.pow(self.z,2)
end
--[[
Normalise the vector (if possible) to length 1.
--]]
function Vec3:normalise()
local l
if self:is_zero() then
print("Unable to normalise a zero-length vector")
return false
end
l = 1/self:len()
return self:scale(l)
end
--[[
Scale the vector.
--]]
function Vec3:scale(l)
return Vec3(self.x * l,self.y * l,self.z * l)
end
--[[
Add vectors.
--]]
function Vec3:add(v)
return Vec3(self.x + v.x, self.y + v.y, self.z + v.z)
end
--[[
Subtract vectors.
--]]
function Vec3:subtract(v)
return Vec3(self.x - v.x, self.y - v.y, self.z - v.z)
end
--[[
Apply a transformation between "absolute" coordinates and "relative"
coordinates. In the "relative" system, xy are in the iPad screen and
z points straight out. In the "absolute" system, y is in the
direction of the gravity vector, x is in the plane of the iPad screen,
and z is orthogonal to those two.
This function interprets the vector as being with respect to the
absolute coordinates and returns the corresponding vector in the
relative system.
--]]
function Vec3:absCoords()
local gxy,l,ga,gb,gc
gxy = Vec3(-Gravity.y,Gravity.x,0)
l = gxy:len()
if l == 0 then
print("Unable to compute coordinate system, gravity vector is (" .. Gravity.x .. "," .. Gravity.y .. "," .. Gravity.z .. ")")
return false
end
ga = gxy:scale(1/l)
gb = Vec3(-Gravity.x,-Gravity.y,-Gravity.z)
gc = Vec3(-Gravity.x * Gravity.z /l, -Gravity.y * Gravity.z /l, l)
return self:applyMatrix(ga,gb,gc)
end
--[[
Determine whether or not the vector is in front of the "eye" (crude
test for visibility).
--]]
function Vec3:isInFront(e)
if not e then
e = Vec3.eye
end
if self:dot(e) < e:dot(e) then
return true
else
return false
end
end
--[[
Project the vector onto the screen using stereographic projection from
the "eye".
--]]
function Vec3:stereoProject(e)
local t,v
if not e then
e = Vec3.eye
end
if self.z == e.z then
-- can't project
return false
end
t = 1 / (1 - self.z / e.z)
v = self:subtract(e)
v = v:scale(t)
v = v:add(e)
-- hopefully v.z is now 0!
return vec2(v.x,v.y)
end
--[[
Partial inverse to stereographic projection: given a point on the
screen and a height, we find the point in space at that height which
would project to the point on the screen.
--]]
function Vec3.stereoInvProject(v,e,h)
local t,u
if not e then
e = Vec3.eye
end
u = Vec3(v.x,v.y,0)
t = h / e.z
u = (1 - t) * u + t * e
-- hopefully u.z is now h!
return u
end
--[[
Returns the distance from the eye; useful for sorting objects.
--]]
function Vec3:stereoLevel(e)
local v
if not e then
e = Vec3.eye
end
v = self:subtract(e)
return e:len() - v:len()
end
--[[
Applies a rotation as specified by another 3-vector, with direction
being the axis and magnitude the angle.
--]]
function Vec3:rotate(w)
local theta, u, v, a, b, c, x
if w:is_zero() then
return self
else
theta = w:len()
w = w:normalise()
if w.x ~= 0 then
u = Vec3(-w.y/w.x,1,0)
u = u:normalise()
else
u = Vec3(1,0,0)
end
v = w:cross(u)
a = self:dot(u)
b = self:dot(v)
c = self:dot(w)
x = w:scale(c)
u = u:scale(a * math.cos(theta) + b * math.sin(theta))
v = v:scale(-a * math.sin(theta) + b * math.cos(theta))
x = x:add(u)
x = x:add(v)
return x
end
end
--[[
Promote to a quaternion with 0 real part.
--]]
function Vec3:toQuaternion()
return Quaternion(0,self.x,self.y,self.z)
end
--[[
Apply a quaternion as a rotation.
--]]
function Vec3:applyQuaternion(q)
local x = self:toQuaternion()
x = q:multiplyRight(x)
x = x:multiplyRight(q:conjugate())
return x:vector()
end
--[[
Inline operators:
u + v
u - v
-u
u * v : cross product (possibly bad choice) or scaling
u / v : scaling
u ^ q : apply quaternion as rotation
u == v : equality
u .. v : dot product (possibly bad choice)
The notation for cross product and dot product may be removed in a
later version.
--]]
function Vec3:__add(v)
return self:add(v)
end
function Vec3:__sub(v)
return self:subtract(v)
end
function Vec3:__unm()
return self:scale(-1)
end
function Vec3:__mul(v)
if type(self) == "number" then
return v:scale(self)
elseif type(v) == "number" then
return self:scale(v)
elseif type(v) == "table" then
if v:is_a(Vec3) then
return self:cross(v)
end
end
return false
end
function Vec3:__div(l)
if type(l) == "number" then
return self:scale(1/l)
else
return false
end
end
function Vec3:__pow(q)
if type(q) == "table" then
if q:is_a(Quaternion) then
return self:applyQuaternion(q)
end
end
return false
end
function Vec3:__eq(v)
return self:is_eq(v)
end
function Vec3:__concat(v)
if type(v) == "table"
and v:is_a(Vec3)
and type(self) == "table"
and self:is_a(Vec3)
then
return self:dot(v)
else
if type(v) == "table"
and v:is_a(Vec3) then
return self .. v:tostring()
else
return self:tostring() .. v
end
end
end
function Vec3:tostring()
return "(" .. self.x .. "," .. self.y .. "," .. self.z .. ")"
end
function Vec3:__tostring()
return self:tostring()
end
function Vec3:tovec3()
return vec3(self.x,self.y,self.z)
end
--[[
The following functions are not class methods but are still related to
vectors.
--]]
--[[
Sets the "eye" for stereographic projection. Input can either be a
Vec3 object or the information required to specify one.
--]]
function Vec3.SetEye(...)
if arg.n == 1 then
Vec3.eye = arg[1]
elseif arg.n == 3 then
Vec3.eye = Vec3(unpack(arg))
else
print("Wrong number of arguments to Vec3.SetEye (1 or 3 expected, got " .. arg.n .. ")")
return false
end
return true
end
--[[
Some useful Vec3 objects.
--]]
Vec3.eye = Vec3(0,0,1)
Vec3.origin = Vec3(0,0,0)
Vec3.e1 = Vec3(1,0,0)
Vec3.e2 = Vec3(0,1,0)
Vec3.e3 = Vec3(0,0,1)
--[[
Is the line segment a-b over c-d when seen from e?
--]]
function Vec3.isOverLine(a,b,c,d,e)
if not e then
e = Vec3.eye
end
-- rebase at a
b = b - a
c = c - a
d = d - a
e = e - a
-- test signs of various determinants
a = c:cross(d)
if a:dot(b) * a:dot(e) < 0 then
return false
end
a = d:cross(e)
if a:dot(b) * a:dot(c) < 0 then
return false
end
a = e:cross(c)
if a:dot(b) * a:dot(d) < 0 then
return false
end
-- right direction, is it far enough?
c = c:subtract(e)
d = d:subtract(e)
a = c:cross(d)
local l = a:dot(b)
local m = a:dot(e)
if l * m > 0 and math.abs(l) > math.abs(m) then
return true
else
return false
end
end
--[[
Is the line segment a-b over the point c when seen from e?
r is the "significant distance"
--]]
function Vec3.isOverPoint(a,b,c,r,e)
if not e then
e = Vec3.eye
end
local aa,bb,ab,ac,bc,d,l
-- rebase at e
a = a - e
b = b - e
c = c - e
d = a:cross(b)
if math.abs(d:dot(c)) > r * d:len() then
return false
end
aa = a:lenSqr()
bb = b:lenSqr()
ab = a:dot(b)
ac = a:dot(c)
bc = b:dot(c)
if aa * bc < ab * ac then
return false
end
if bb * ac < ab * bc then
return false
end
l = math.sqrt((aa * bb - ab * ab) * (aa + bb - 2 * ab))
if (bb - ab) * ac + (aa - ab) * bc < aa * bb - ab * ab + r * l then
return false
end
return true
end
return Vec3
--]==]
--[==[
-- Vector class
-- Author: Andrew Stacey
-- Website: http://www.math.ntnu.no/~stacey/HowDidIDoThat/iPad/Codea.html
-- Licence: CC0 http://wiki.creativecommons.org/CC0
--[[
The "Vector" class is for handling arbitrary dimensional vectors
and defining a variety of methods on them.
--]]
local Vector = class()
--[[
A vector is a size and an array of numbers.
--]]
function Vector:init(...)
local v
if arg.n == 1 then
if type(arg[1]) == "table" then
if arg[1].is_a and arg[1]:is_a(Complex) then
v = {arg[1]}
else
v = arg[1]
end
elseif type(arg[1]) == "userdata" then
local mt = getmetatable(arg[1])
if mt == getmetatable(vec2()) then
v = {arg[1].x,arg[1].y}
elseif mt == getmetatable(vec3()) then
v = {arg[1].x,arg[1].y,arg[1].z}
elseif mt == getmetatable(vec4()) then
v = {arg[1].x,arg[1].y,arg[1].z,arg[1].w}
end
else
v = {tonumber(arg[1])}
end
else
v = arg
end
--local u = {}
local n = 0
for k,l in ipairs(v) do
table.insert(self,l)
--table.insert(u,l)
n = n + 1
end
--self.vec = u
self.size = n
-- Shortcuts for low dimension
self.x = self[1]
self.y = self[2]
self.z = self[3]
self.w = self[4]
end
--[[
Test for zero vector.
--]]
function Vector:is_zero()
for k,v in ipairs(self) do
if v ~= 0 then
return false
end
end
return true
end
--[[
Test for equality.
--]]
function Vector:is_eq(u)
if self.size ~= u.size then
return false
end
for k,v in ipairs(self) do
if v ~= u[k] then
return false
end
end
return true
end
--[[
Inner product.
--]]
function Vector:dot(u)
if self.size ~= u.size then
return false
end
local d = 0
for k,v in ipairs(self) do
d = d + v * u[k]
end
return d
end
--[[
Apply a given matrix (which is specified as a list of row vectors).
--]]
function Vector:applyMatrixLeft(m)
if m.cols ~= self.size then
return false
end
local u = {}
for k,v in ipairs(m) do
table.insert(u,self:dot(v))
end
return Vector(u)
end
function Vector:applyMatrixRight(m)
if m.rows ~= self.size then
return false
end
local u = {}
local a = Vector.zero(m.cols)
for k,v in ipairs(m) do
a = a + self[k]*v
end
return a
end
--[[
Length of the vector
--]]
function Vector:len()
local d = 0
for k,v in ipairs(self) do
d = d + math.pow(v,2)
end
return math.sqrt(d)
end
--[[
Squared length of the vector.
--]]
function Vector:lenSqr()
local d = 0
for k,v in ipairs(self) do
d = d + math.pow(v,2)
end
return d
end
--[[
Norm infinity
--]]
function Vector:linfty()
local m = 0
for k,v in ipairs(self) do
m = math.max(m,math.abs(v))
end
return m
end
--[[
Norm one
--]]
function Vector:lone()
local m = 0
for k,v in ipairs(self) do
m = m + math.abs(v)
end
return m
end
--[[
Normalise the vector (if possible) to length 1.
--]]
function Vector:normalise()
local l
if self:is_zero() then
print("Unable to normalise a zero-length vector")
return false
end
l = 1/self:len()
return self:scale(l)
end
--[[
Scale the vector.
--]]
function Vector:scale(l)
local u = {}
for k,v in ipairs(self) do
table.insert(u,l*v)
end
return Vector(u)
end
--[[
Add vectors.
--]]
function Vector:add(u)
if self.size ~= u.size then
return false
end
local w = {}
for k,v in ipairs(self) do
table.insert(w, v + u[k])
end
return Vector(w)
end
--[[
Subtract vectors.
--]]
function Vector:subtract(u)
if self.size ~= u.size then
return false
end
local w = {}
for k,v in ipairs(self) do
table.insert(w, v - u[k])
end
return Vector(w)
end
--[[
Inline operators:
u + v
u - v
-u
u * v : scaling
u / v : scaling
u == v : equality
--]]
function Vector:__add(v)
return self:add(v)
end
function Vector:__sub(v)
return self:subtract(v)
end
function Vector:__unm()
return self:scale(-1)
end
function Vector:__mul(v)
if type(self) == "number" then
return v:scale(self)
elseif type(v) == "number" then
return self:scale(v)
else
if self.is_a and self:is_a(Matrix) then
return v:applyMatrixLeft(self)
elseif v.is_a and v:is_a(Matrix) then
return self:applyMatrixRight(v)
end
end
return false
end
function Vector:__div(l)
if type(l) == "number" then
return self:scale(1/l)
else
return false
end
end
function Vector:__eq(v)
return self:is_eq(v)
end
function Vector:__concat(v)
if type(v) == "table"
and v:is_a(Vector) then
return self .. v:tostring()
else
return self:tostring() .. v
end
end
function Vector:tostring()
local t = {}
for k,v in ipairs(self) do
if type(v) == "table" and v.tostring then
table.insert(t,v:tostring())
else
table.insert(t,v)
end
end
return "(" .. table.concat(t,",") .. ")"
end
function Vector:__tostring()
return self:tostring()
end
function Vector:tovec()
if self.size == 2 then
return vec2(self[1],self[2])
elseif self.size == 3 then
return vec3(self[1],self[2],self[3])
else
return vec4(self[1],self[2],self[3],self[4])
end
end
function Vector.zero(n)
if type(n) == "table" then
n = n.size
end
u = {}
for i = 1,n do
table.insert(u,0)
end
return Vector(u)
end
function Vector:bitwiseReorder(b)
local v
if b then
v = self
else
v = Vector({})
end
local h = BinLength(self.size)
local m = math.pow(2,h)
local l
local w = {}
for k = 1,m do
if not w[k] then
l = BinReverse(k-1,h)+1
v[k],v[l] = self[l] or 0,self[k] or 0
w[k],w[l] = 1,1
end
end
v.size = m
return v
end
function Vector:FFT(t)
t = t or {}
local v
v = self:bitwiseReorder(t.inplace)
local pi = math.pi
if not t.inverse then
pi = - pi
end
local r = v.size
local n = BinLength(r) -1
local i,j,d,s,m,fi,p,pr
for k = 0,n do
i = math.pow(2,k)
j = math.pow(2,k+1)
d = pi/i
s = math.sin(d/2)
m = Complex(-2 * s * s, math.sin(d))
fi = Complex(1,0)
for l = 1,i do
for ll = l,r-1,j do
p = ll + i
pr = fi * v[p]
v[p] = v[ll] - pr
v[ll] = v[ll] + pr
end
fi = m*fi + fi
end
end
return v
end
if cmodule.loaded "TestSuite" then
testsuite.addTest({
name = "Vector",
setup = function()
local w = Vector({1,2,3,4,5,6,7,8})
local z = Vector({Complex(1,1), Complex(2,3), Complex(4,5)})
for k,v in ipairs({
{"Vector", w},
{"Complex Vector", z},
{"Bitwise Reverse", w:bitwiseReorder()},
{"Bitwise Reverse", z:bitwiseReorder()},
{"FFT", w:FFT()},
{"FFT", z:FFT()}
}) do
print(v[1] .. ": " .. v[2])
end
end,
draw = function()
end
})
end
return Vector
--]==]
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.