Flying Ignatz Release v1.12 -Using quaternions to fly a plane.

1aTabOrder

Flying Ignatz Tab Order Version: 1.12
------------------------------
This file should not be included in the Codea project.
#Main
#Flight
#Joystick
#Plane
#Quaternion
#Sky
#VecExt

Flight.lua

--Flight class
 
Flight=class()
 


function Flight:init(t) 
    t = t or {}
    self.mode= t.mode or vec3(0,1,0) --- AAA
    self.sens=t.sens or vec3(1,1,1)*0.25
    self.limit=t.limit or 2*vec3(45,90,22.5)*math.pi/180  --- AAA
    self.steady=t.steady or 5
    self:Reset()
end
 
function Flight:Reset() 
    self.dv=vec3(0,0,0)
    self.pr = vec4(1,0,0,0)
    self.y = vec4(1,0,0,0)
    --CCCC initialise
    self.cam = {
        at = vec3(0,0,1),
        frame = vec4(1,0,0,0),
        vel = vec3(0,0,0),
        avel = vec3(0,0,0)}
end

--[[
self.turning needs to use a ternary logic so we can detect not only its state but also when it flips
--]]
 
function Flight:Turn(r)
    self.dv = self.mode * vec3(-r.value.y,-r.value.x,-r.value.x) * self.sens * DeltaTime
    self.pr = self.pr * qTangent((1 - self.mode) * vec3(-r.delta.y,-r.delta.x,-r.delta.x) * self.limit)
    self.y = self.y * qTangent(self.dv)
    self.turning = 1
end
 
function Flight:Position(s) --d is vec3 distance in pixels
    return s^(self.y * self.pr)
end
 
function Flight:Rotate()
    self:Interpolate()     --RRRR update interpolation when rotating
    rotate(self.y * self.pr)
end
 
function Flight:rotation()
    return self.y * self.pr 
end
 
function Flight:Interpolate()  --RRRR this has the interpolation code from Main
    --don't update if we've done it already in this timeslice
    if self.stime==ElapsedTime then return end 
    --interpolate
    if self.turning == 2 then
        self.interpolating = true
        self.stime = ElapsedTime
        self.slerpQ = self:InterpolateSlerp()
        self.turning = 0
    elseif self.turning == 1 then
        self.turning = 2
    elseif self.interpolating then
        self.interpolating = self.slerpQ(ElapsedTime)
    end
end
 
function Flight:InterpolateSlerp() --RRRR this is now called from Interpolate above
    local rt, stime, l, rl, dy, y
    stime = ElapsedTime
    dy = self.dv
    rt = vec4(self.pr.x,0,self.pr.z,0):normalise()
    l = self.pr:sdist(rt) * self.steady
    rl = self.pr:make_slerp(rt)
    return function(t)
        t = t - stime
        if t > l then
            self.pr = vec4(1,0,0,0)
            self.y = self.y * rt
            return false
        end
        self.y = self.y * qTangent(edge(t,l,0) * dy)
        self.pr = rl(smootherstep(t,0,l))
        return true
    end
end
 
 -- CCCC new function to manage camera position
--p = vec3, plane position
--offset = vec3, location of camera relative to self, in object space
--lag = scalar, camera follow lag, from 0 (no lag) to 5 (very slow)
function Flight:SetCamera(p,offset,lag)
    if not lag or lag==0 then --slavishly follow plane
        self.cam.at = p + self:Position(offset)
        self.cam.frame= self:rotation()
    else  --lag plane rotation
        local ca,cq = self.cam.at,self.cam.frame
        self.cam.at = self.cam.at + self.cam.vel * DeltaTime/lag
        self.cam.frame = self.cam.avel:exp(DeltaTime/lag) * self.cam.frame
        self.cam.vel = (p + self:Position(offset) - ca)
        self.cam.avel = (self:rotation() * cq^""):log()
    end
    local up = vec3(0,1,0)^self.cam.frame
    --NOTE - addition to end of next line, so camera looks in front of plane (and not at it)
    -- local look = self.cam.at - offset^self.cam.frame + self:Position(vec3(0,0,-50))
    local look = p + self:Position(vec3(0,0,-20))
    camera(self.cam.at,look,up)
end
 

Joystick.lua

 
JoyStick = class()
 
function JoyStick:init(t)
    t = t or {}
    self.radius = t.radius or 100
    self.stick = t.stick or 30
    self.centre = t.centre or self.radius * vec2(1,1) + vec2(5,5)
    self.position = vec2(0,0)
    self.target = vec2(0,0)
    self.value = vec2(0,0)
    self.delta = vec2(0,0)
    self.mspeed = 60
    self.moving = 0
end
 
function JoyStick:draw()
    -- Codea does not automatically call this method
    ortho()
    viewMatrix(matrix())
    pushStyle()
    fill(160, 182, 191, 50)
    stroke(118, 154, 195, 100)
    strokeWidth(1) 
    ellipse(self.centre.x,self.centre.y,2*self.radius)
    fill(78, 131, 153, 50)
    ellipse(self.centre.x+self.position.x, self.centre.y+self.position.y, self.stick*2)
    popStyle()
end
 
function JoyStick:touched(t)
    if t.state == BEGAN then
        local v = vec2(t.x,t.y)
        if v:dist(self.centre)<self.radius-self.stick then
            self.mytouch = t.id
        end
    end
    if t.id == self.mytouch then
        if t.state~=ENDED then
            local v = vec2(t.x,t.y)
            if v:dist(self.centre)>self.radius-self.stick then
                v = (v - self.centre):normalize()*(self.radius - self.stick) + self.centre
            end  --set x,y values for joy based on touch
            self.target=v - self.centre
        else --reset joystick to centre when touch ends
            self.target=vec2(0,0)
            self.mytouch = false
        end
    end
end
 
function JoyStick:update()
    local p = self.target - self.position
    if p:len() < DeltaTime * self.mspeed then
        self.position = self.target
        if not self.mytouch then
            if self.moving ~= 0 then
                self.moving = self.moving - 1
            end
        else
            self.moving = 2
        end
    else
        self.position = self.position + p:normalize() * DeltaTime * self.mspeed
        self.moving = 2
    end
    local v = self.position/(self.radius - self.stick)
    self.delta = v - self.value
    self.value = v
end

function JoyStick:isMoving()
    return self.moving
end
 

Main.lua

VERSION = "1.12"
 
function setup()
    -- [[
    if AutoGist then
        autogist = AutoGist("Flying Ignatz","Using quaternions to fly a plane.",VERSION)
        autogist:backup(true)
    end
    -- ]]

    --**Flight controls **
    speed=vec3(0,0,-60)  --pixels/sec 

    --** Plane settings **
    plane=Plane() --plane mesh
    planeQ=Flight()--plane quaternion
    planePos=vec3(0,0,0) --starting position
    camOffset={vec3(0,3,-9),vec3(0,0,50),vec3(0,40,100)} --camera position relative to plane
    parameter.integer("state",-1,5,0)
    parameter.integer("View",1,#camOffset,2)
    --** Joystick settings **
    rotjoy = JoyStick({centre = vec2(105,105)})
    mvjoy = JoyStick({centre = vec2(WIDTH - 105,105)})

    --set up sky globe
    sky=SkyGlobe("Documents:SkyDome","Documents:SkyDome")  --SSSS
end
 
function draw()
    rotjoy:update()
    mvjoy:update()
    background(179, 205, 220, 255)
    fill(255)
    perspective()
    camera(0,0,100,0,0,0)   
    RotateAndDraw()    
    rotjoy:draw()
    mvjoy:draw()
end
 
--this function adjusts the rotation based on joystick position, and draws the plane
function RotateAndDraw()
    if rotjoy.mytouch then
        planeQ:Turn(rotjoy)
    end
    if mvjoy:isMoving() then
        planePos=planePos+planeQ:Position(mvjoy.value.y*speed*DeltaTime)
    end
 
    -- CCCC -see changes in this set of if statements, and SetCamera function at bottom of Flight tab
    if state == -1 then -- stationary camera
        --user can set it any way they want, this is not part of the flight library
        camera(0,0,50,0,0,-1000) --this just looks straight ahead
    else -- slavishly follows plane
        planeQ:SetCamera(planePos,camOffset[View],state)
    end
    pushMatrix() 
    translate(planePos)
    sky:draw()             --SSSS
    planeQ:Rotate()
    plane:draw()
    popMatrix()
    
end
 
--** Utility functions - joystick and touch ***
 
--manage joystick
function touched(t)
    rotjoy:touched(t)
    mvjoy:touched(t)
end
 
 
 

Plane.lua

--plane
 
Plane=class()
 
function Plane:init()
    a={}
    c1=color(107, 139, 184, 255)
    c2=color(72, 148, 195, 255)
    c3=color(176, 159, 101, 255)
    c4=color(255,0,0)
    a[1]=Block(3,3,25,c1,vec3(0,0,0))   
    a[2]=Block(15,1,5,c2,vec3(-9,1,-6))
    a[3]=Block(15,1,5,c2,vec3(9,1,-6))
    a[4]=Block(1,4,3,c2,vec3(0,3.5,11))
    a[5]=Block(4,1,2,c2,vec3(-3.5,1,11))
    a[6]=Block(4,1,2,c2,vec3(3.5,1,11))
    a[7]=Block(3,2,5,c3,vec3(0,2.5,-6))
    a[8]=Block(3,3,1,c4,vec3(0,0,-13))
    self.blocks=a 
end
 
function Plane:draw()
    for i=1,#self.blocks do
        self.blocks[i]:draw()
    end 
end
 
Block = class() --taken from 3D lab project
 
function Block:init(w,h,d,c,p,t,r) --width,height,depth,colour,position,texture, (optional) texture range (see above)
    self.width=w
    self.height=h
    self.depth=d  
    self.color=c  
    self.pos=p
    self.tex=t 
    --if no limits specified on which part of image to draw, set to 0,1
    if r~=nil then self.texR=r else self.texR={0,0,1,1} end
    self.blk=self:createBlock()
end
 
function Block:createBlock()
    -- all the unique vertices that make up a block
    --There are only 8 corners in a cube - we define them as vertices
    --all measurements are taken from the centre of the block
    --so bottom left front has x of -1/2 width, y of -1/2 height, and z of 1/2 depth
    local w,h,d=self.width,self.height,self.depth
    local x,y,z=self.pos.x,self.pos.y,self.pos.z
    local v = {
      vec3(x-0.5*w, y-0.5*h, z+0.5*d), -- Left  bottom front
      vec3(x+0.5*w, y-0.5*h, z+0.5*d), -- Right bottom front
      vec3(x+0.5*w, y+0.5*h, z+0.5*d), -- Right top    front
      vec3(x-0.5*w, y+0.5*h, z+0.5*d), -- Left  top    front
      vec3(x-0.5*w, y-0.5*h, z-0.5*d), -- Left  bottom back
      vec3(x+0.5*w, y-0.5*h, z-0.5*d), -- Right bottom back
      vec3(x+0.5*w, y+0.5*h, z-0.5*d), -- Right top    back
      vec3(x-0.5*w, y+0.5*h, z-0.5*d), -- Left  top    back
    }
 
    local cubeverts = {
      -- Front, Right, Back, Left, Top, Bottom
      v[1], v[2], v[3], v[1], v[3], v[4],
      v[2], v[6], v[7], v[2], v[7], v[3],
      v[6], v[5], v[8], v[6], v[8], v[7],
      v[5], v[1], v[4], v[5], v[4], v[8],
      v[4], v[3], v[7], v[4], v[7], v[8],
      v[5], v[6], v[2], v[5], v[2], v[1],
    }
    local cc={}
     --assign colors
        local c=self.color
        local x=0.5
        local BL=color(c.r,c.g,c.b) --bottom left
        local BR=color(c.r,c.g,c.b) --bottom right  
        local TR=color(c.r*x,c.g*x,c.b*x) --top right
        local TL=color(c.r*x,c.g*x,c.b*x) --top left
        for i=1,6 do
            cc[#cc+1]=BL
            cc[#cc+1]=BR
            cc[#cc+1]=TR
            cc[#cc+1]=BL
            cc[#cc+1]=TR
            cc[#cc+1]=TL
        end
  
    --put it all together
    local ms = mesh()
    ms.vertices = cubeverts
    if self.tex then
        ms.texture = self.tex
        ms.texCoords = cubetexCoords
    else
        ms.colors=cc
    end
    --ms:setColors(self.color)
    return ms
end  
    
function Block:draw()
    self.blk:draw() 
end
 

Quaternion.lua

--[==[
-- Quaternion library wherein vec4s are promoted to quaternions
-- Author: Andrew Stacey
-- Website: http://www.math.ntnu.no/~stacey/HowDidIDoThat/iPad/Codea.html
-- Licence: CC0 http://wiki.creativecommons.org/CC0
 
--[[
This provides code for handling rotations.  Internally, the rotations
are implemented as quaternions which are in turn implemented by
extending the methods available for vec4 userdata.
 
Due to how vec4s are implemented, we cannot have methods beginning with
letters a, b, g, r, w, x, y, z
--]]
 
-- Simplistic error handling
local error = error or print
 
-- Localise the maths functions/constants that we use for faster lookup.
local abs = math.abs
local pow = math.pow
local sqrt = math.sqrt
local sin = math.sin
local cos = math.cos
local acos = math.acos
local asin = math.asin
local pi = math.pi
local floor = math.floor
local min = math.min
local max = math.max
 
--[[
The function "is_a" extends the capabilities of the method "is_a" which is automatically defined by Codea for classes.
 
Parameters:
a: object to be tested
b: test
 
The tests work as follows.
 
1. If the type of b is a string, it is taken as the name of a type to test a against.
2. If the type of b is a table, it is assumed to be a class to test if a is an instance thereof.
3. If the type of b is a userdata, the test is to see if a is the same type of object.
4. If b is a function, then it is replaced by the value of that function.
--]]
 
local function is_a(a,b)
    if type(b) == "function" then
        b = b()
    end
    if type(b) == "string" then
        return type(a) == b
    end
    if type(b) == "table" then
        if type(a) == "table"
            and a.is_a
            and a:is_a(b) then
                return true
        else
            return false
        end
    end
    if type(b) == "userdata" then
        if type(a) == "userdata" then
            a = getmetatable(a)
            b = getmetatable(b)
            return a == b
        else
            return false
        end
    end
    return false
end

function edge(t,a,b)
    a = a or 0
    b = b or 1
    t = (t-a)/(b-a)
    return min(1,max(0,t))
end

function smoothstep(t,a,b)
    a = a or 0
    b = b or 1
    t = (t-a)/(b-a)
    t = min(1,max(0,t))
    return t * t * (3 - 2 * t)
end

function smootherstep(t,a,b)
    a = a or 0
    b = b or 1
    t = (t-a)/(b-a)
    t = min(1,max(0,t))
    return t * t * t * (t * (t * 6 - 15) + 10)
end
 
-- Import the metatable of vec4s for extending
local mtq = getmetatable(vec4())
 
-- We run into difficulties if our vec4s contain NaNs or infs.  This
-- tests for those.
mtq["is_finite"] = function(q)
    if q.x < math.huge and q.x > -math.huge 
        and q.y < math.huge and q.y > -math.huge 
        and q.z < math.huge and q.z > -math.huge 
        and q.w < math.huge and q.w > -math.huge 
    then
        return true
    end
    return false
end
 
--[[
Test if we are real.
--]]
 
mtq["is_real"] = function (q)
    if q.y ~= 0 or q.z ~= 0 or q.w ~= 0 then
        return false
    end
    return true
end
 
--[[
Test if the real part is zero.
--]]
 
mtq["is_imaginary"] = function (q)
    if q.x ~= 0 then
        return false
    end
    return true
end
 
--[[
Normalise a quaternion to have length 1, safely.  Here "safely" means
that we ensure that we do not have NaNs or infs.  If we do, the
returned quaternion is the unit.
--]]
 
mtq["normalise"] = function (q)
    q = q:normalize()
    if q:is_finite() then
        return q
    else
        return vec4(1,0,0,0)
    end
end
 
--[[
Spherical length and distance.
--]]
 
mtq["slen"] = function(q)
    q = q:normalise()
    q.x = q.x - 1
    return 2*asin(q:len()/2)
end
 
mtq["sdist"] = function(q,qq)
    q = q:normalise()
    qq = qq:normalise()
    return 2*asin(q:dist(qq)/2)
end
 
--[[
Add two quaternions inline, including promotion of a number.
 
q + p
--]]
 
local __add = mtq["__add"]
 
mtq["__add"] = function (a,b)
    if is_a(a,"number") then
        a = vec4(a,0,0,0)
    end
    if is_a(b,"number") then
        b = vec4(b,0,0,0)
    end
    return __add(a,b)
end
 
--[[
Same for inline subtraction.
 
q - p
--]]
 
local __sub = mtq["__sub"]
 
mtq["__sub"] = function (a,b)
    if is_a(a,"number") then
        a = vec4(a,0,0,0)
    end
    if is_a(b,"number") then
        b = vec4(b,0,0,0)
    end
    return __sub(a,b)
end
 
--[[
For inline multiplication, we also allow multiplication by a matrix.
In this case, we use the fact that we're viewing quaternions as
rotations and matrices as affine transformations, so the best result
of multiplication is as the matrix representing the appropriate
composition.
 
q * p
--]]
 
local __mul = mtq["__mul"]
 
mtq["__mul"] = function (a,b)
    if is_a(a,"number") then
        return __mul(a,b)
    end
    if is_a(b,"number") then
        return __mul(a,b)
    end
    if is_a(a,matrix) then
        return a:__mul(b:tomatrixleft())
    end
    if is_a(b,matrix) then
        a = a:tomatrixleft()
        return a:__mul(b)
    end
    local x,y,z,w
    x = a.x * b.x - a.y * b.y - a.z * b.z - a.w * b.w
    y = a.x * b.y + a.y * b.x + a.z * b.w - a.w * b.z
    z = a.x * b.z - a.y * b.w + a.z * b.x + a.w * b.y
    w = a.x * b.w + a.y * b.z - a.z * b.y + a.w * b.x
    return vec4(x,y,z,w)
end
 
--[[
Conjugation (corresponds to inverting a rotation).
--]]
 
mtq["conjugate"] = function (q)
    return vec4(q.x, - q.y, - q.z, - q.w)
end
 
mtq["co"] = mtq["conjugate"]
 
--[[
Inline division, including division of and by a number.
--]]
 
local __div = mtq["__div"]
 
mtq["__div"] = function (a,b)
    if is_a(b,"number") then
        return __div(a,b)
    end
    local l = b:lenSqr()
    b = vec4(b.x/l,-b.y/l,-b.z/l,-b.w/l)
    if is_a(a,"number") then
        return vec4(a*b.x,a*b.y,a*b.z,a*b.w)
    end
    local x,y,z,w
    x = a.x * b.x - a.y * b.y - a.z * b.z - a.w * b.w
    y = a.x * b.y + a.y * b.x + a.z * b.w - a.w * b.z
    z = a.x * b.z - a.y * b.w + a.z * b.x + a.w * b.y
    w = a.x * b.w + a.y * b.z - a.z * b.y + a.w * b.x
    return vec4(x,y,z,w)
end
 
--[[
Powers.  Integer powers are defined iteratively.  Non-integer powers
use the slerp function.
 
TODO: test whether the slerp is faster than the iterative method for
integer powers as well.
--]]
 
local function intpower(q,n)
    if n == 0 then
        return vec4(1,0,0,0)
    elseif n > 0 then
       return q:__mul(intpower(q,n-1))
    elseif n < 0 then
       local l = q:lenSqr()
       q = vec4(q.x/l,-q.y/l,-q.z/l,-q.w/l)
       return q:intpower(-n)
    end
end
 
local function power(q,n)
   if n == floor(n) then
      return intpower(q,n)
   end
   local l = q:len()
   q = q:normalise()
   return l^n * q:slerp(n)
end
 
--[[
Inline exponentiation.
 
q^n
 
The behaviour depends on the type of the exponent:
 
* number: compute the power.
* vec4: conjugate by the exponent.
* other: return the conjugate.
--]]
 
mtq["__pow"] = function (q,n)
    if is_a(n,"number") then
        return power(q,n)
    elseif is_a(n,vec4) then
        return n:__mul(q):__div(n)
    else
        return q:conjugate()
    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.
 
Parameters:
q initial rotation
qq final rotation
t parameter (only for the direct functions)
 
If the input does not provide enough parameters, it is assumed that
the missing one is the initial rotation and that this should be taken
to be the identity rotation (represented by vec4(1,0,0,0)).
--]]
 
 
--[[
Linear interpolation, renormalised.
--]]
 
mtq["lerp"] = function (q,qq,t)
    if not t then
        return vec4(1,0,0,0):lerp(q,qq)
    end
    local v
    if (q + qq):len() == 0 then
        -- antipodal points, need a midpoint
        v = vec4(q.y,-q.x,q.w,-q.z)
        v = (1 - 2*t)*q + (1-abs(2*t-1))*v
    else
        v = (1-t)*q + t*qq
    end
    return v:normalise()
end
 
--[[
Spherical interpolation.
 
We have to be quite careful here not to emit anything with NaNs or
infs as there are several opportunities to create them due to the
limits of finite precision mathematics.
--]]
 
mtq["slerp"] = function (q,qq,t)
    if not t then
        return vec4(1,0,0,0):slerp(q,qq)
    end
    local v
    if (q + qq):len() == 0 then
        -- antipodal points, need a midpoint
        v = vec4(q.y,-q.x,q.w,-q.z)
        t = 2*t
    elseif (q - qq):len() == 0 then
        return q
    else
        v = qq
    end
    local ca = q:dot(v)
    local sa = sqrt(1 - pow(ca,2))
    if sa == 0 or sa ~= sa then
        return q
    end
    local a = acos(ca)
    sa = sin(a*t)/sa
    v = (cos(a*t) - ca*sa)*q+ sa*v
    return v
end
 
--[[
Constructor for normalised linear interpolation.
--]]
 
mtq["make_lerp"] = function (q,qq)
    if not qq then
        return vec4(1,0,0,0):make_lerp(q)
    end
    local v,w
    w = q:normalise()
    if (q + qq):len() == 0 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-abs(2*t-1))*v
            return u:normalise()
        end
    else
        v = qq:normalise()
        return function(t)
            local u = (1-t)*w + t*v
            return u:normalise()
        end
    end
end
 
--[[
Spherical interpolation
--]]
 
mtq["make_slerp"] = function (q,qq)
    if not qq then
        q,qq = vec4(1,0,0,0),q
    end
    local v,f,u
    if (q + qq):len() == 0 then
        -- antipodal points, need a midpoint
        v = vec4(q.y,-q.x,q.w,-q.z)
        f = 2
    elseif (q - qq):len() == 0 then
        return function(t)
            return q
        end
    else
        v = qq
        f = 1
    end
    v = v:normalise()
    u = q:normalise()
    local ca = u:dot(v)
    local sa = sqrt(1 - pow(ca,2))
    if sa == 0 or sa ~= sa then
        return function(t)
            return q
        end
    end
    local a = acos(ca)
    v = (v - ca*q)/sa
    return function(t)
        return cos(a*f*t)*q + sin(a*f*t)*v
    end
end
 
--[[
Returns the real part.
--]]
 
mtq["toreal"] = function (q)
    return q.x
end
 
--[[
Returns the vector (imaginary) part as a vec3 object.
--]]
 
mtq["vector"] = function (q)
    return vec3(q.y, q.z, q.w)
end
 
mtq["tovector"] = mtq["vector"]
 
mtq["log"] = function (q)
    local l = q:slen()
    local v = q:tovector()
    v = v:normalize()
    if not v:is_finite() then
        return vec3(0,0,0)
    else
        return v * l
    end
end
                                    
--[[
Represents a quaternion as a string.
--]]
 
mtq["tostring"] = function (q)
    local s
    local im = {{q.y,"i"},{q.z,"j"},{q.w,"k"}}
    if q.x ~= 0 then
        s = string.format("%.3f",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
 
 
mtq["__concat"] = function (q,s)
    if is_a(s,"string") then
        return q:tostring() .. s
    else
        return q .. s:tostring()
    end
end
 
--[[
Converts the quaternion to a matrix.
 
We distinguish between left and right actions.  A rotation matrix
acting on the right is the transpose of the matrix on the left.
--]]
 
mtq["tomatrixleft"] = function (q)
    q = q:normalise()
    local a,b,c,d = q.x,q.y,q.z,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
 
mtq["tomatrixright"] = function (q)
    q = q:normalise()
    local a,b,c,d = q.x,-q.y,-q.z,-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
 
mtq["tomatrix"] = mtq["tomatrixright"]
 
--[[
Converts the quaternion to an "angle-axis" representation.
--]]
 
mtq["toangleaxis"] = function (q)
    q = q:normalise()
    local a = 2*acos(q.x)
    local v = vec3(q.y,q.z,q.w)
    if v == vec3(0,0,0) then
        return 0,vec3(0,0,1)
    end
    return a,v:normalise()
end
 
--[[
The following code modifies various of Codea's functions to enable
them to take a quaternion as input.  We have to be careful as for some
passing nil is distinct to passing an empty input.
 
The other issue is as to the distinction between left and right
actions of matrices.  In OpenGL, matrices act on the right.  However,
in standard linear algebra, the convention is for matrices to act on
the left.
--]]
 
local __modelMatrix = modelMatrix
 
function modelMatrix(m)
    if m then
        if is_a(m,vec4) then
            m = m:tomatrixright()
        end
        return __modelMatrix(m)
    else
        return __modelMatrix()
    end
end
 
local __applyMatrix = applyMatrix
 
function applyMatrix(m)
    if m then
        if is_a(m,vec4) then
            m = m:tomatrixright()
        end
        return __applyMatrix(m)
    else
        return __applyMatrix()
    end
end
 
local __viewMatrix = viewMatrix
 
function viewMatrix(m)
    if m then
        if is_a(m,vec4) then
            m = m:tomatrixright()
        end
        return __viewMatrix(m)
    else
        return __viewMatrix()
    end
end
 
local __projectionMatrix = projectionMatrix
 
function projectionMatrix(m)
    if m then
        if is_a(m,vec4) then
            m = m:tomatrixright()
        end
        return __projectionMatrix(m)
    else
        return __projectionMatrix()
    end
end
 
local __rotate = rotate
 
function rotate(a,x,y,z)
    if is_a(a,vec4) then
        local v
        a,v = a:toangleaxis()
        x,y,z = v.x,v.y,v.z
        a = a*180/pi
    end
    return __rotate(a,x,y,z)
end
 
local __translate = translate
 
function translate(x,y,z)
    if not y then
        x,y,z = x.x,x.y,x.z
    end
    return __translate(x,y,z)
end
                                    
 
local __scale = scale
 
function scale(a,b,c)
    if is_a(a,vec3) then
        a,b,c = a.x,a.y,a.z
    end
    if c then
        return __scale(a,b,c)
    end
    if b then
        return __scale(a,b)
    end
    if a then
        return __scale(a)
    end
end
                                    
local __camera = camera
 
function camera(a,b,c,d,e,f,g,h,i)
    if is_a(a,vec3) then
        a,b,c,d,e,f,g,h,i = a.x,a.y,a.z,b,c,d,e,f,g                              
    end
    if is_a(d,vec3) then
        d,e,f,g,h,i = d.x,d.y,d.z,e,f,g
    end
    if is_a(g,vec3) then
        g,h,i = g.x,g.y,g.z
    end
    if g then
        return __camera(a,b,c,d,e,f,g,h,i)
    elseif d then
        return __camera(a,b,c,d,e,f)
    elseif a then
        return __camera(a,b,c)
    else
        return __camera()
    end
end
 
--[[
We also define some extensions to the vec3 type.
--]]
 
local mt = getmetatable(vec3())
 
-- Same as for vec4.
mt["is_finite"] = function(v)
    if v.x < math.huge and v.x > -math.huge 
        and v.y < math.huge and v.y > -math.huge 
        and v.z < math.huge and v.z > -math.huge 
    then
        return true
    end
    return false
end
 
--[[
Promote to a quaternion with 0 real part.
--]]
 
mt["toQuaternion"] = function (v)
    return vec4(0,v.x,v.y,v.z)
end
 
--[[
Apply a quaternion as a rotation (assumes unit quaternion for speed)
using conjugation.
--]]
 
mt["applyQuaternion"] = function (v,q)
   v = v:toQuaternion()
    v = q:__mul(v)
    v = v:__mul(q:conjugate())
    return v:vector()
end
 
-- There is no native rotate method for a vec3.
mt["rotate"] = function(v,q,x,y,z)
    if is_a(q,"number") then
        q = qRotation(q,x,y,z)
    end
    return v:applyQuaternion(q)
end
 
-- We use the exponential notation for writing the action of a
-- quaternion on a vector (this is consistent with group theory
-- notation).
mt["__pow"] = function (v,q)
    if is_a(q,vec4) then
        return v:applyQuaternion(q)
    end
    return false
end
 
mt["__concat"] = function (u,s)
    if is_a(s,"string") then
        return u:__tostring() .. s
    else
        return u .. s:__tostring()
    end
end
 
--[[
The rotateTo method produces a rotation that rotates the first vector
to the second about an axis that is perpendicular to both.  So long as
the two are not collinear, the axis is unique.  The angle is taken to
be the smallest angle.  If the vectors point in opposite directions,
an orthogonal axis is chosen in such a way as to minimise precision
error.
 
TODO: Allow for a third parameter to specify the axis in case of
ambiguity.
--]]
mt["rotateTo"] = function (u,v)
    if v:lenSqr() == 0 or u:lenSqr() == 0 then
        return vec4(1,0,0,0)
    end
    u = u:normalise()
    v = u + v:normalise()
    if v:lenSqr() == 0 then
        -- Opposite vectors, no canonical direction
        local a,b,c = abs(u.x), abs(u.y), 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:normalise()
    local d = u:dot(v)
    u = u:cross(v)
    return vec4(d,u.x,u.y,u.z)
end
 
--[[
Safe renormalisation, as for quaternions.  Except that there is no
"cannonical" unit vec3, so we return the unit vector in the z
direction.  The rationale for that is that if this is an axis, this
vector represents the axis out of the screen.
--]]
mt["normalise"] = function (v)
    v = v:normalize()
    if v:is_finite() then
        return v
    else
        return vec3(0,0,1)
    end
end
 
--[[
Inline multiplication extended to allow for multiplication by a matrix.
 
We interpret the vec3 as vertical or horizontal depending on whether
the matrix is on the right or left.  Multiplication by a matrix is
viewed as applying an affine transformation.
--]]
local __mulv = mt["__mul"]
 
mt["__mul"] = function(m,v)
    if is_a(m,vec3) and is_a(v,"number") then
        return __mulv(m,v)
    end
    if is_a(m,"number") and is_a(v,vec3) then
        return __mulv(m,v)
    end
    if is_a(m,vec3) and is_a(v,vec3) then
        return vec3(m.x*v.x,m.y*v.y,m.z*v.z)
    end
    if is_a(m,matrix) and is_a(v,vec3) then
        local l = m[13]*v.x + m[14]*v.y + m[15]*v.z + m[16]
        return vec3(
            (m[1]*v.x + m[2]*v.y + m[3]*v.z + m[4])/l,
            (m[5]*v.x + m[6]*v.y + m[7]*v.z + m[8])/l,
            (m[9]*v.x + m[10]*v.y + m[11]*v.z + m[12])/l
        )
    end
    if is_a(m,vec3) and is_a(v,matrix) then
       m,v = v,m
        local l = m[4]*v.x + m[8]*v.y + m[12]*v.z + m[16]
        return vec3(
            (m[1]*v.x + m[5]*v.y + m[9]*v.z + m[13])/l,
            (m[2]*v.x + m[6]*v.y + m[10]*v.z + m[14])/l,
            (m[3]*v.x + m[7]*v.y + m[11]*v.z + m[15])/l
        )
    end
end
 
local __addv = mt["__add"]
 
mt["__add"] = function(a,b)
    if is_a(a,"number") then
        a = vec3(a,a,a)
    end
    if is_a(b,"number") then
        b = vec3(b,b,b)
    end
    return __addv(a,b)
end
 
local __subv = mt["__sub"]
 
mt["__sub"] = function(a,b)
    if is_a(a,"number") then
        a = vec3(a,a,a)
    end
    if is_a(b,"number") then
        b = vec3(b,b,b)
    end
    return __subv(a,b)
end
                    
 
--[[
Extensions to the matrix class.
--]]
local mtm = getmetatable(matrix())
 
--[[
Inline multiplication by either quaternion or vector.
--]]
local __mulm = mtm["__mul"]
 
mtm["__mul"] = function (m,mm)
    if is_a(m,matrix) and is_a(mm,matrix) then
        return __mulm(m,mm)
    end
    if is_a(m,matrix) and is_a(mm,vec4) then
        return __mulm(m,mm:tomatrix())
    end
    if is_a(m,vec4) and is_a(mm,matrix) then
        return __mulm(m:tomatrix(),mm)
    end
    if is_a(m,matrix) and is_a(mm,vec3) then
        local l = m[13]*mm.x + m[14]*mm.y + m[15]*mm.z + m[16]
        return vec3(
            (m[1]*mm.x + m[2]*mm.y + m[3]*mm.z + m[4])/l,
            (m[5]*mm.x + m[6]*mm.y + m[7]*mm.z + m[8])/l,
            (m[9]*mm.x + m[10]*mm.y + m[11]*mm.z + m[12])/l
        )
    end
    if is_a(m,vec3) and is_a(mm,matrix) then
        local l = mm[4]*m.x + mm[8]*m.y + mm[12]*m.z + mm[16]
        return vec3(
            (mm[1]*m.x + mm[5]*m.y + mm[9]*m.z + mm[13])/l,
            (mm[2]*m.x + mm[6]*m.y + mm[10]*m.z + mm[14])/l,
            (mm[3]*m.x + mm[7]*m.y + mm[11]*m.z + mm[15])/l
        )
    end
end
 
-- Extending the rotate method to take a quaternion.
local __mrotate = mtm["rotate"]
 
mtm["rotate"] = function(m,a,x,y,z)
    if is_a(a,vec4) then
        a,x = a:toangleaxis()
        x,y,z = x.x,x.y,x.z
    end
    return __mrotate(m,a,x,y,z)
end
 
 
--[[
The following functions are intended to be more "user friendly" and
provide solutions to common problems.
 
The primary functions construct a rotation dependent on some initial
data.  The simplest are to define a rotation from some other method of
specifying rotations: angle-axis or Euler angles.  More complicated
are methods to define a rotation by giving two (orthogonal) frames and
computing the rotation to get from the first to the last.  In
particular, there are certain common frames that might be expected to
be used frequently:
 
1. The initial orientation of the iPad.
 
2. A gravitational frame wherein the Gravity vector is "straight
down".  This is not completely defined as there is currently no way to
choose a corresponding horizontal direction (access to the compass
would provide this).  Two possibilities are: to choose the x-axis so
that it is always in the plane of the iPad, and to use RotationRate to
try to keep track of the initial x-axis.
 
3. A "DeltaOrientation" in which the change in orientation from one
frame to the next is used (this uses RotationRate).
 
There are obvious limits in accuracy in using these.
--]]
 
--[[
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.
--]]
 
local function qGravity()
    local gxy, gy, gygxy, a, b, c, d
    if Gravity.x == 0 and Gravity.y == 0 then
        return vec4(1,0,0,0)
    else
        gy = - Gravity.y
        gxy = sqrt(pow(Gravity.x,2) + pow(Gravity.y,2))
        gygxy = gy/gxy
        a = sqrt(1 + gxy - gygxy - gy)/2
        b = sqrt(1 - gxy - gygxy + gy)/2
        c = sqrt(1 - gxy + gygxy - gy)/2
        d = 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 vec4(a,b,c,d)
    end
end
 
function qrGravity(q)
    local qg,qx
    if not q then
        qx = ReferenceFrame()
    else
        qx = q
    end
    local y = vec3(0,-1,0)^qx
    qg = y:rotateTo(Gravity)
    return qg*qx
end
 
mtq["gravity"] = qrGravity
 
local frame = {
        vec3(1,0,0),
        vec3(0,1,0),
        vec3(0,0,1)
    }
local qzyx = vec4(1,0,0,0)
 
--[[
Needs to be run once every frame!
--]]
 
function updateReferenceFrame(q)
    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
    --]]
    x,y,z = unpack(frame)
    local qz = qRotation(RotationRate.z*DeltaTime,z.x,z.y,z.z)
    local qy = qRotation(RotationRate.y*DeltaTime,y.x,y.y,y.z)
    local qx = qRotation(RotationRate.x*DeltaTime,x.x,x.y,x.z)
    if q then
        local qq = qz * qy * qx * q
        q.x,q.y,q.z,q.w = qq.x,qq.y,qq.z,qq.w
    else
        qzyx = qz * qy * qx * qzyx
        return qzyx
    end
end
 
function ReferenceFrame(q)
    return q or qzyx
end
 
mtq["updateReferenceFrame"] = updateReferenceFrame
mtq["ReferenceFrame"] = ReferenceFrame
 
--[[
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 qRotation(a,x,y,z)
    local q,c,s
    if not y then
        x,y,z = x.x,x.y,x.z
    end
    q = vec4(0,x,y,z)
    q = q:normalise()
    if q == vec4(1,0,0,0) then
        return q
    end
    c = cos(a/2)
    s = sin(a/2)
    q = q:__mul(s)
    q = q:__add(c)
    return q
end
 
--[[
The qEuler function handles conversion from Euler angles.  As there
are many ways to specify Euler angles, we have to allow for quite a
variety.  The table __euler contains the code necessary for allowing
all of this flexibility.
 
The last parameter is the specification, which is a table taken from
the alphabet "xXyYzZ".  This ought to correspond to the axes of
rotation, using lowercase to mean the absolute frame and uppercase to
mean the transformed frame.
 
TODO: Test this.
 
If only two parameters are passed, the first is taken to be a vec3 or
a table containing the Euler angles.
--]]
 
local __euler = {}
__euler.x = function(q)
    return vec3(1,0,0)^q
end
__euler.X = function(q)
    return vec3(1,0,0)
end
__euler.y = function(q)
    return vec3(0,1,0)^q
end
__euler.Y = function(q)
    return vec3(0,1,0)
end
__euler.z = function(q)
    return vec3(0,0,1)^q
end
__euler.Z = function(q)
    return vec3(0,0,1)
end
 
function qEuler(a,b,c,v)
    if c then
        a = {a,b,c}
    else
        if is_a(a,vec3) then
            a = {a.x,a.y,a.z}
        end
        v = b
    end
    v = v or {"x","y","z"}
    local q = vec4(1,0,0,0)
    local w
    for k,u in ipairs(v) do
        w = __euler[u](q)
        q = q * qRotation(a[k],w)
    end
    return q
end
 
-- v is a tangent vector at the identity
function qTangent(x,y,z,t)
    local q
    if is_a(x,"number") then
        q = vec4(0,x,y,z)
        t = t or 1
    else
        q = vec4(0,x.x,x.y,x.z)
        t = y or 1
    end
    local qn = q:normalise()
    if qn == vec4(1,0,0,0) then
        return qn
    end
    local a = t * q:len() --/2
    return cos(a)*vec4(1,0,0,0) + sin(a)*qn
end
 
mt["exp"] = qTangent
 
function qRotationRate()
    return qTangent(DeltaTime * RotationRate)
end
 
local iGravity = vec3(0,-1,0)
Frame = {
    Roll = vec4(1,0,0,0),
    Pitch = vec4(1,0,0,0),
    Yaw = vec4(1,0,0,0),
    Gravity = vec4(1,0,0,0),
    Rotation = vec4(1,0,0,0),
    AdjustedRotation = vec4(1,0,0,0),
    RotationRate = vec4(1,0,0,0)
}
 
local update_fn 
update_fn = function()
    local q = vec3(iGravity.x,iGravity.y,0):rotateTo(
            vec3(Gravity.x,Gravity.y,0))
    Frame.Roll.x = q.x
    Frame.Roll.y = q.y
    Frame.Roll.z = q.z
    Frame.Roll.w = q.w
    q = vec3(0,iGravity.y,iGravity.z):rotateTo(
            vec3(0,Gravity.y,Gravity.z))
    Frame.Pitch.x = q.x
    Frame.Pitch.y = q.y
    Frame.Pitch.z = q.z
    Frame.Pitch.w = q.w
    q = vec3(iGravity.x,0,iGravity.z):rotateTo(
            vec3(Gravity.x,0,Gravity.z))
    Frame.Yaw.x = q.x
    Frame.Yaw.y = q.y
    Frame.Yaw.z = q.z
    Frame.Yaw.w = q.w
    q = iGravity:rotateTo(Gravity)
    Frame.Gravity.x = q.x
    Frame.Gravity.y = q.y
    Frame.Gravity.z = q.z
    Frame.Gravity.w = q.w
    q = vec4(0,RotationRate.x,RotationRate.y,RotationRate.z)
    local qn = q:normalise()
    if qn ~= vec4(1,0,0,0) then
        local a = q:len()*DeltaTime/2
        qn = cos(a)*vec4(1,0,0,0) + sin(a)*qn
    end
    Frame.RotationRate.x = qn.x
    Frame.RotationRate.y = qn.y
    Frame.RotationRate.z = qn.z
    Frame.RotationRate.w = qn.w
    q = qn*Frame.Rotation
    Frame.Rotation.x = q.x
    Frame.Rotation.y = q.y
    Frame.Rotation.z = q.z
    Frame.Rotation.w = q.w
    q = qn*Frame.AdjustedRotation
    local g = iGravity^(q)
    qn = g:rotateTo(iGravity)
    q = q*qn
    Frame.AdjustedRotation.x = q.x
    Frame.AdjustedRotation.y = q.y
    Frame.AdjustedRotation.z = q.z
    Frame.AdjustedRotation.w = q.w
    tween.delay(0,update_fn)
end
    
 
tween.delay(.5,function()
     iGravity = vec3(Gravity.x,Gravity.y,Gravity.z)
     tween.delay(0,update_fn)
     end)
 
--[[
A suite of test functions.
--]]
function testRotations()
    print("Rotation tests:")
    local tolerance = 10^(-6)
    local vars = [[local q1,q2,i,j,k
    q1 = vec4(1,2,3,4)
    q2 = vec4(.5,.5,-.5,-.5)
    i = vec4(0,1,0,0)
    j = vec4(0,0,1,0)
    k = vec4(0,0,0,1)
    q3 = vec4(1/math.sqrt(2),0,0,1/math.sqrt(2))
    v1 = vec3(1,0,0)
    v2 = vec3(0,1,0)
    v3 = vec3(1,-1,2)
    m = matrix(
        10,10,5,1,
        4,6,4,1,
        1,3,3,1,
        0,1,2,1
    )
    return ]]
    local tests = {
        {"q1 * q2", vec4(3,2,4,-1)},
        {"q1 + q2", vec4(1.5,2.5,2.5,3.5)},
        {"q1 - q2", vec4(.5,1.5,3.5,4.5)},
        {"q1 / q2", vec4(-2,0,-1,5)},
        {"q1^q2", vec4(1,-4,-2,3)},
        {"q1^\"\"", vec4(1,-2,-3,-4)},
        {"i * j", vec4(0,0,0,1)},
        {"j * k", vec4(0,1,0,0)},
        {"k * i", vec4(0,0,1,0)},
        {"v1^q3", vec3(0,1,0)},
        {"q3:tomatrixleft()*v1", vec3(0,1,0)},
        {"v1:rotate(q3)", vec3(0,1,0)},
        {"v2^q3", vec3(-1,0,0)},
        {"q3:tomatrixleft()*v2", vec3(-1,0,0)},
        {"v2:rotate(q3)", vec3(-1,0,0)},
        {"qRotation(math.pi/2,0,0,1)", vec4(1/math.sqrt(2),0,0,1/math.sqrt(2))},
        {"v1:rotateTo(v2)", vec4(1/math.sqrt(2),0,0,1/math.sqrt(2))},
        {"v3 * m", vec3(8,11,9)/3},
        {"m * v3", vec3(11,7,5)/4}
    }
    local t,a,m,n,tn
    tn,n = 0,0
    for k,u in ipairs(tests) do
        tn = tn + 1
        t = loadstring(vars .. u[1])
        a = t()
        if a == u[2] then
            m = "OK"
            n = n + 1
        elseif a.dist and a:dist(u[2]) < tolerance then
            m = "OK (" .. a:dist(u[2]) .. ")"
            n = n + 1
        else
            m = "not OK (expected " .. u[2] .. ")"
        end
        print(u[1] .. " = " .. a .. " : " .. m)
    end
    print(n .. " tests passed out of " .. tn)
end
 
 
 
 
 
--]==]

Sky.lua

-- Sky globe
-- wraps two images around the top and bottom halves of a sphere to create sky and horizons
-- if you only want the sky, leave out the second image
-- the images need to be twice as wide as they are high, and specially stretched to fit on a sphere
-- so search for them on the internet
 
SkyGlobe = class()
 
function SkyGlobe:init(i1,i2)
    --create image to wrap on sphere
    local img1=readImage(i1)  --get first image
    local sky=image(img1.width,img1.width/2)  --set up image to wrap on sphere
    setContext(sky) --start drawing on it
    spriteMode(CORNER)
    sprite(img1,0,sky.height-img1.height)  --draw one or two images, as provided
    if i2 then sprite(i2,0,sky.height-img1.height*2) end
    setContext()
    
    --now the sphere
    --the code below comes from Jmv38, who explains it here
    --http://jmv38.comze.com/CODEAbis/server.php (look for 3D tutorial)
    self.r=500 --arbitrary radius, it doesn't matter how big you make it, it looks the same
    -- create mesh and colors
    local vertices,tc = {},{}
    vertices,tc = self:optimMesh({ nx=40, ny=20 })
    vertices = self:warpVertices({verts=vertices, xangle=180, yangle=180 }) 
 
    -- create the mesh itself
    self.ms = mesh()
    self.ms.vertices = vertices
    self.ms.texture = sky
    self.ms.texCoords = tc
end
 
function SkyGlobe:optimMesh(input)
    -- create the mesh tables
    local vertices = {}
    local texCoords = {}
    local k = 0
    local s = 1
    -- create a set of triangles with approx constant surface on a sphere
    local x,y
    local x1,x2 = {},{}
    local i1,i2 = 0,0
    local nx,ny = input.nx,input.ny
    local sx, sy = nx/ny, 1/ny
    local center = vec3(1,0.5,0)
    local m1,m2,c
    local flip = -1
    for y=0,ny-1 do -- for each horizontal band
        -- number of points on each side of the band
        local nx1 = math.floor( nx * math.abs(math.cos(( y*sy-0.5)*2 * math.pi/2)))
        if nx1<6 then nx1=6 end
        local nx2 = math.floor( nx * math.abs(math.cos(((y+1)*sy-0.5)*2 * math.pi/2))) 
        if nx2<6 then nx2=6 end
        -- points on each side of the band
        x1,x2 = {},{}
        for i1 = 1,nx1 do x1[i1] = (i1-1)/(nx1-1)*sx end
        for i2 = 1,nx2 do x2[i2] = (i2-1)/(nx2-1)*sx end
        x1[nx1+1] = x1[nx1] -- just a trick to manage last triangle without thinking
        x2[nx2+1] = x2[nx2]
        -- start on the left
        local i1,i2 = 1,1
        local continue = true
        local n,nMax = 0,0
        nMax = nx*2+1
        while continue do
            -- center of the 2 current segments
            m1 = (x1[i1]+x1[i1+1])/2
            m2 = (x2[i2]+x2[i2+1])/2
            if m1<=m2 then -- the less advanced base makes the triangle
                vertices[k+1] = vec3(   x1[i1], sy*y    , 1)  - center
                vertices[k+2] = vec3( x1[i1+1], sy*y    , 1)  - center
                vertices[k+3] = vec3(   x2[i2], sy*(y+1), 1)  - center
                texCoords[k+1] = vec2(   x1[i1]/2*flip, sy*y    ) 
                texCoords[k+2] = vec2( x1[i1+1]/2*flip, sy*y    )
                texCoords[k+3] = vec2(   x2[i2]/2*flip, sy*(y+1))
                if i1<nx1 then i1 = i1 +1 end
            else
                vertices[k+1] = vec3(   x1[i1], sy*y    , 1) - center
                vertices[k+2] = vec3(   x2[i2], sy*(y+1), 1) - center
                vertices[k+3] = vec3( x2[i2+1], sy*(y+1), 1) - center
                texCoords[k+1] = vec2(   x1[i1]/2*flip, sy*y    )
                texCoords[k+2] = vec2(   x2[i2]/2*flip, sy*(y+1))
                texCoords[k+3] = vec2( x2[i2+1]/2*flip, sy*(y+1))
                if i2<nx2 then i2 = i2 +1 end        
            end
            if c==1 then c=2 else c=1 end
            if i1==nx1 and i2==nx2 then continue=false end
            -- increment index for next triangle
            k = k + 3
            n = n + 1
            if n>nMax then continue=false  end -- just in case of infinite loop
        end
    end   
    return vertices,texCoords
end
 
 
function SkyGlobe:warpVertices(input)
    -- move each vector to its position on sphere
    local verts = input.verts
    local xangle = -input.xangle/180*math.pi
    local yangle = -input.yangle/180*math.pi   
    for i,v in ipairs(verts) do
        verts[i] = vec3(0,0,1):rotate(xangle*v[2],1,0,0):rotate(yangle*v[1],0,1,0)
    end
    return verts
end
 
function SkyGlobe:draw()
    pushMatrix()
    scale(self.r)
    self.ms:draw()  
    popMatrix()
end
 

VecExt.lua

-- Extensions to the native Codea vector and matrix types
-- Author: Andrew Stacey
-- Website: http://loopspace.mathforge.com
-- Licence: CC0 http://wiki.creativecommons.org/CC0

local ab,ac,ad,ae,af,ag,am,as,ay,ca,ch,co,cp,er,eu,ex,fl,ip,is,lo,m,ma,mi,mm,mr,mu,mv,mw,mx,pi,pm,po,pr,qp,ra,ro,rp,sc,sh,si,sq,su,sv,sw,sy,ta,tc,th,tp,tr,ts,vm
er,m = error or print,math
ab,po,sq,si,co,ac,as,pi,fl,mi,ma,ex,lo,ta,sh,ch,th = m.abs,m.pow,m.sqrt,m.sin,m.cos,m.acos,m.asin,m.pi,m.floor,m.min,m.max,m.exp,m.log,m.tan,m.sinh,m.cosh,m.tanh
mm,am,vm,pm,ro,tr,sc,ca = modelMatrix,applyMatrix,viewMatrix,projectionMatrix,rotate,translate,scale,camera
function is(a,b)
    if type(b) == "function" then b = b() end
    if type(b) == "string" then return type(a) == b end
    if type(b) == "table" and type(a) == "table" and a.is_a then return a:is_a(b) end
    if type(b) == "userdata" and type(a) == "userdata" then return  getmetatable(a) == getmetatable(b) end
    return false
end

function edge(t,a,b)
    a,b = a or 0,b or 1
    return mi(1,ma(0,(t-a)/(b-a)))
end
function smoothstep(t,a,b)
    a,b = a or 0,b or 1
    t = mi(1,ma(0,(t-a)/(b-a)))
    return t * t * (3 - 2 * t)
end
function smootherstep(t,a,b)
    a,b = a or 0,b or 1
    t = mi(1,ma(0,(t-a)/(b-a)))
    return t * t * t * (t * (t * 6 - 15) + 10)
end

sy = readLocalData("Complex Symbol","i")
ra = readLocalData("Complex Angle","rad")
if ra == "rad" then ag,ay = 1,"π" else ag,ay = 180,"°" end
pr = readLocalData("Complex Precision",2)

function setComplex(t)
    ag,ay,pr,sy,ts = t.angle or ag,t.angsym or ay,t.precision or pr,t.symbol or sy,t.tostring or ts
end
m = getmetatable(vec2())
m["clone"] = function (c) return vec2(c.x,c.y) end
m["is_finite"] = function (c) if c.x < math.huge and c.x > -math.huge and c.y < math.huge and c.y > -math.huge then return true end return false end
m["is_real"] = function (c) return c.y == 0 end
m["is_imaginary"] = function (c) return c.x == 0 end
m["normalise"] = function (c) c=c:normalize() if c:is_finite() then return c else return vec2(1,0) end end
ad,su,mu = m["__add"],m["__sub"],m["__mul"]
m["__add"] = function (a,b)
    if is(a,"number") then a = vec2(a,0) end
    if is(b,"number") then b = vec2(b,0) end
    return ad(a,b)
end
m["__sub"] = function (a,b)
    if is(a,"number") then a = vec2(a,0) end
    if is(b,"number") then b = vec2(b,0) end
    return su(a,b)
end
m["__mul"] = function (a,b)
    if is(a,"number") then a = vec2(a,0) end
    if is(b,"number") then b = vec2(b,0) end
    return vec2(a.x*b.x - a.y*b.y,a.x*b.y+a.y*b.x)
end
m["conjugate"] = function (c) return vec2(c.x, - c.y) end
m["co"] = m["conjugate"]
function rp(c,n,k)
    k = k or 0
    local r,t = pow(c:len(),n), (k*2*pi-c:angleBetween(vec2(1,0)))*n
    return vec2(r*cos(t),r*sin(t))
end
function cp(c,w,k)
    if is(w,"number") then return rp(c,w,k) end
    if c == vec2(0,0) then
        er("Taking powers of 0 is somewhat dubious")
        return false
    end
    local r,t = c:len(),-c:angleBetween(vec2(1,0))
    k = k or 0
    local nr,nt = pow(r,w.x) * exp(-w.y*t), (t + k * 2 * pi) * w.x + log(r) * w.y
    return vec2(nr*cos(nt),nr*sin(nt))
end
m["__pow"] = function (c,n)
    if is(n,"number") then return rp(c,n)
    elseif is(n,vec2) then return cp(c,n)
    else return c:co() end
end
m["__div"] = function (c,q)
    if is(q,"number") then return vec2(c.x/q,c.y/q)
    elseif is(c,"number") then return c/q:lenSqr()*vec2(q.x,-q.y)
    else return vec2(c.x*q.x+c.y*q.y,c.y*q.x-c.x*q.y)/q:lenSqr() end
end
m["real"] = function (c) return c.x end
m["imaginary"] = function (c) return c.y end
m["__concat"] = function (c,v) if is(v,vec2) then return c .. v:tostring() else return c:tostring() .. v end end
m["tostring"] = function (c) return ts(c) end
function tc(c)
    local s
    local x,y = fl(c.x * 10^pr +.5)/10^pr,fl(c.y * 10^pr +.5)/10^pr
    if x ~= 0 then s = x end
    if y ~= 0 then if s then if y > 0 then
        if y == 1 then s = s .. " + " .. sy
        else s = s .. " + " .. y .. sy end
    else
        if y == -1 then s = s .. " - " .. sy
        else s = s .. " - " .. (-y) .. sy end
    end
    else
        if y == 1 then s = sy
        elseif y == - 1 then s = "-" .. sy
        else s = y .. sy end
    end
    end
    s = s or "0"
    return s
end
function tp (c)
    local t,r = fl(ag *nc:arg() * 10^pr/pi +.5)/10^pr,fl(c:len() * 10^pr +.5)/10^pr
    return "(" .. r .. "," .. t .. angsym .. ")"
end
ts = tc
m["topolarstring"] = tp
m["tocartesianstring"] = tc
m["arg"] = function (c) return -c:angleBetween(vec2(1,0)) end
function Complex_unit() return vec2(1,0) end
function Complex_zero() return vec2(0,0) end
function Complex_i() return vec2(0,-1) end
function math.abs(n)
    if is(n,"number") then return ab(n) end
    if is(n,vec2) then return n:len() end
    er("Cannot take the length of " .. n)
end
function math.pow(n,e)
    if is(n,"number") then
        if is(e,"number") then return po(n,e)
        elseif is(e,vec2) then return rp(vec2(n,0),e,0) end
    end
    if is(n,vec2) then return cp(n,e,0) end
    er("Cannot take the power of " .. n .. " by " .. e)
end
function math.sqrt(n)
    if is(n,"number") then return sq(n) end
    if is(n,vec2) then return rp(n,.5,0) end
    er("Cannot take the square root of " .. n)
end
function math.exp(n)
    if is(n,"number") then return ex(n) end
    if is(n,vec2) then
        local r = ex(n.x)
        return vec2(r*cos(n.y),r*sin(n.y))
    end
    er("Cannot exponentiate " .. n)
end
function math.cos(n)
    if is(n,"number") then return co(n) end
    if is(n,vec2) then return vec2(co(n.x)*ch(n.y),-si(n.x)*sh(n.y)) end
    er("Cannot take the cosine of " .. n)
end
function math.sin(n)
    if is(n,"number") then return si(n) end
    if is(n,vec2) then return vec2(si(n.x)*ch(n.y),co(n.x)*sh(n.y)) end
    error("Cannot take the sine of " .. n)
end
function math.cosh(n)
    if is(n,"number") then return ch(n) end
    if is(n,vec2) then return vec2(ch(n.x)*co(n.y), sh(n.x)*si(n.y)) end
    er("Cannot take the hyperbolic cosine of " .. n)
end
function math.sinh(n)
    if is(n,"number") then return sh(n) end
    if is(n,vec2) then return vec2(sh(x)*co(y), ch(x)*si(y)) end
    er("Cannot take the hyperbolic sine of " .. n)
end
function math.tan(n)
    if is(n,"number") then return ta(n) end
    if is(n,vec2) then
            local cx,sx,chy,shy = co(n.x),si(n.x),ch(n.y),sh(n.y)
            local d = cx^2 * chy^2 + sx^2 * shy^2
            if d == 0 then return false end
            return vec2(sx*cx/d,shy*chy/d)
    end
    er("Cannot take the tangent of " .. n)
end
function math.tanh(n)
    if is(n,"number") then return th(n) end
    if is(n,vec2) then
            local cx,sx,chy,shy = co(n.x),si(n.x),ch(n.y),sh(n.y)
            local d = cx^2 * chy^2 + sx^2 * shy^2
            if d == 0 then return false end
            return vec2(shx*chx/d,sy*cy/d)
    end
    er("Cannot take the hyperbolic tangent of " .. n)
end
function math.log(n,k)
    if is(n,"number") then
        if k then return vec2(log(n), 2*k*pi)
        else return log(n) end
    end
    k = k or 0
    if is(n,vec2) then return vec2(lo(n:len()),n:arg() + 2*k*pi) end
    er("Cannot take the logarithm of " .. n)
end
m = getmetatable(vec4())
m["is_finite"] = function(q)
    if q.x < math.huge and q.x > -math.huge and q.y < math.huge and q.y > -math.huge and q.z < math.huge and q.z > -math.huge and q.w < math.huge and q.w > -math.huge then return true end
    return false
end
m["is_real"] = function (q)
    if q.y ~= 0 or q.z ~= 0 or q.w ~= 0 then return false end
    return true
end
m["is_imaginary"] = function (q) return q.x == 0 end
m["normalise"] = function (q) q = q:normalize() if q:is_finite() then return q else return vec4(1,0,0,0) end
end
m["slen"] = function(q) 
    q = q:normalise()
    q.x = q.x - 1
    return 2*as(q:len()/2)
end
m["sdist"] = function(q,qq)
    q = q:normalise()
    qq = qq:normalise()
    return 2*as(q:dist(qq)/2)
end
ae,sv,mv,dv = m["__add"],m["__sub"],m["__mul"],m["__div"]
m["__add"] = function (a,b)
    if is(a,"number") then a = vec4(a,0,0,0) end
    if is(b,"number") then b = vec4(b,0,0,0) end
    return ae(a,b)
end
m["__sub"] = function (a,b)
    if is(a,"number") then a = vec4(a,0,0,0) end
    if is(b,"number") then b = vec4(b,0,0,0) end
    return sv(a,b)
end
m["__mul"] = function (a,b)
    if is(a,"number") then return mv(a,b) end
    if is(b,"number") then return mv(a,b) end
    if is(a,matrix) then return a:__mul(b:tomatrixleft()) end
    if is(b,matrix) then return a:tomatrixleft():__mul(b) end
    return vec4(a.x*b.x-a.y*b.y-a.z*b.z-a.w*b.w,a.x*b.y+a.y*b.x+a.z*b.w-a.w*b.z,a.x*b.z-a.y*b.w+a.z*b.x+a.w*b.y,a.x*b.w+a.y*b.z-a.z*b.y+a.w*b.x)
end
m["conjugate"] = function (q) return vec4(q.x, - q.y, - q.z, - q.w) end
m["co"] = m["conjugate"]
m["__div"] = function (a,b)
    if is(b,"number") then return dv(a,b) end
    local l = b:lenSqr()
    b = vec4(b.x/l,-b.y/l,-b.z/l,-b.w/l)
    if is(a,"number") then return vec4(a*b.x,a*b.y,a*b.z,a*b.w) end
    return vec4(a.x*b.x-a.y*b.y-a.z*b.z-a.w*b.w,a.x*b.y+a.y*b.x+a.z*b.w-a.w*b.z,a.x*b.z-a.y*b.w+a.z*b.x+a.w*b.y,a.x*b.w+a.y*b.z-a.z*b.y+a.w*b.x)
end
function ip(q,n)
    if n == 0 then return vec4(1,0,0,0)
    elseif n > 0 then return q:__mul(ip(q,n-1))
    elseif n < 0 then
       local l = q:lenSqr()
       q = vec4(q.x/l,-q.y/l,-q.z/l,-q.w/l)
       return q:ip(-n)
    end
end
function qp(q,n)
   if n == fl(n) then return ip(q,n) end
   local l = q:len()
   q = q:normalise()
   return l^n * q:slerp(n)
end
m["__pow"] = function (q,n)
    if is(n,"number") then return qp(q,n)
    elseif is(n,vec4) then return n:__mul(q):__div(n)
    else return q:conjugate() end
end
m["lerp"] = function (q,qq,t)
    if not t then q,qq,t = vec4(1,0,0,0),q,qq end
    if (q + qq):len() == 0 then q = (1-2*t)*q+(1-ab(2*t-1))*vec4(q.y,-q.x,q.w,-q.z)
    else q = (1-t)*q + t*qq end
    return q:normalise()
end
m["slerp"] = function (q,qq,t)
    if not t then q,qq,t = vec4(1,0,0,0),q,qq end
    if (q + qq):len() == 0 then qq,t = vec4(q.y,-q.x,q.w,-q.z),2*t
    elseif (q - qq):len() == 0 then return q end
    local ca = q:dot(qq)
    local sa = sq(1 - po(ca,2))
    if sa == 0 or sa ~= sa then return q end
    local a = ac(ca)
    sa = si(a*t)/sa
    return (co(a*t)-ca*sa)*q+sa*qq
end
m["make_lerp"] = function (q,qq)
    if not qq then q,qq = vec4(1,0,0,0),q end
    q,qq = q:normalise(),qq:normalise()
    if (q + qq):len() == 0 then
        qq = vec4(q.y,-q.x,q.w,-q.z)
        return function(t) return ((1-2*t)*q+(1-ab(2*t-1))*qq):normalise() end
    else return function(t) return ((1-t)*q+t*qq):normalise() end
    end
end
m["make_slerp"] = function (q,qq)
    if not qq then q,qq = vec4(1,0,0,0),q end
    q,qq = q:normalise(),qq:normalise()
    local f
    if (q + qq):len() == 0 then qq,f = vec4(q.y,-q.x,q.w,-q.z),2
    elseif (q - qq):len() == 0 then return function(t) return q end
    else f = 1 end
    local ca = q:dot(qq)
    local sa = sq(1 - po(ca,2))
    if sa == 0 or sa ~= sa then return function(t) return q end end
    local a = ac(ca)
    qq = (qq - ca*q)/sa
    return function(t) return co(a*f*t)*q + si(a*f*t)*qq end
end
m["toreal"] = function (q) return q.x end
m["vector"] = function (q) return vec3(q.y, q.z, q.w) end
m["tovector"] = m["vector"]
m["log"] = function (q)
    local l = q:slen()
    q = q:tovector():normalize()
    if not q:is_finite() then return vec3(0,0,0)
    else return q * l end
end
m["tostring"] = function (q)
    local s
    local im = {{q.y,"i"},{q.z,"j"},{q.w,"k"}}
    if q.x ~= 0 then s = string.format("%.3f",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
    s = s or "0"
    return s
end
m["__concat"] = function (q,s)
    if is(s,"string") then return q:tostring() .. s else return q .. s:tostring() end
end
m["tomatrixleft"] = function (q)
    q = q:normalise()
    local a,b,c,d = q.x,q.y,q.z,q.w
    local ab,ac,ad,bb,bc,bd,cc,cd,dd = 2*a*b,2*a*c,2*a*d,2*b*b,2*b*c,2*b*d,2*c*c,2*c*d,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
m["tomatrixright"] = function (q)
    q = q:normalise()
    local a,b,c,d = q.x,-q.y,-q.z,-q.w
    local ab,ac,ad,bb,bc,bd,cc,cd,dd = 2*a*b,2*a*c,2*a*d,2*b*b,2*b*c,2*b*d,2*c*c,2*c*d,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
m["tomatrix"] = m["tomatrixright"]
m["toangleaxis"] = function (q)
    q = q:normalise()
    local a = q.x
    q = vec3(q.y,q.z,q.w)
    if q == vec3(0,0,0) then return 0,vec3(0,0,1) end
    return 2*ac(a),q:normalise()
end
function qGravity()
    if Gravity.x == 0 and Gravity.y == 0 then return vec4(1,0,0,0)
    else
        local gxy, gy, gygxy, a, b, c, d
        gy,gxy = - Gravity.y,sq(po(Gravity.x,2) + po(Gravity.y,2))
        gygxy = gy/gxy
        a,b,c,d = sq(1 + gxy - gygxy - gy)/2,sq(1 - gxy - gygxy + gy)/2,sq(1 - gxy + gygxy - gy)/2,sq(1 + gxy + gygxy + gy)/2
        if Gravity.z < 0 then b,c = - b,-c end
        if Gravity.x > 0 then c,d = - c,-d end
        return vec4(a,b,c,d)
    end
end
function qRotation(a,x,y,z)
    local q,c,s
    if not y then x,y,z = x.x,x.y,x.z end
    q = vec4(0,x,y,z):normalise()
    if q == vec4(1,0,0,0) then return q end
    return q:__mul(si(a/2)):__add(co(a/2))
end
eu = {}
eu.x = function(q) return vec3(1,0,0)^q end
eu.X = function(q) return vec3(1,0,0) end
eu.y = function(q) return vec3(0,1,0)^q end
eu.Y = function(q) return vec3(0,1,0) end
eu.z = function(q) return vec3(0,0,1)^q end
eu.Z = function(q) return vec3(0,0,1) end
function qEuler(a,b,c,v)
    if c then a = {a,b,c}
    else if is(a,vec3) then a = {a.x,a.y,a.z} end v = b end
    v = v or {"x","y","z"}
    local q = vec4(1,0,0,0)
    for k,u in ipairs(v) do
        q = q * qRotation(a[k],eu[u](q))
    end
    return q
end
function qTangent(x,y,z,t)
    local q
    if is(x,"number") then q,t = vec4(0,x,y,z),t or 1
    else q,t = vec4(0,x.x,x.y,x.z),y or 1 end
    local qn = q:normalise()
    if qn == vec4(1,0,0,0) then return qn end
    t = t * q:len()
    return co(t)*vec4(1,0,0,0) + si(t)*qn
end
function qRotationRate() return qTangent(DeltaTime * RotationRate) end
function modelMatrix(m)
    if m then if is(m,vec4) then m = m:tomatrixright() end return mm(m) else return mm() end
end
function applyMatrix(m)
    if m then if is(m,vec4) then m = m:tomatrixright() end return am(m) else return am() end
end
function viewMatrix(m)
    if m then if is(m,vec4) then m = m:tomatrixright() end return vm(m) else return vm() end
end
function projectionMatrix(m)
    if m then if is(m,vec4) then m = m:tomatrixright() end return pm(m) else return pm() end
end
function rotate(a,x,y,z)
    if is(a,vec4) then
        local v
        a,v = a:toangleaxis()
        x,y,z,a = v.x,v.y,v.z,a*180/pi
    end
    if x then return ro(a,x,y,z) end
    return ro(a)
end
function translate(x,y,z)
    if not y then x,y,z = x.x,x.y,x.z end
    if z then return tr(x,y,z) end
    return tr(x,y)
end
function scale(a,b,c)
    if is(a,vec3) then a,b,c = a.x,a.y,a.z end
    if c then return sc(a,b,c) end
    if b then return sc(a,b) end
    if a then return sc(a) end
    return sc()
end
function camera(a,b,c,d,e,f,g,h,i)
    if is(a,vec3) then a,b,c,d,e,f,g,h,i = a.x,a.y,a.z,b,c,d,e,f,g                              end
    if is(d,vec3) then d,e,f,g,h,i = d.x,d.y,d.z,e,f,g end
    if is(g,vec3) then g,h,i = g.x,g.y,g.z end
    if g then return ca(a,b,c,d,e,f,g,h,i)
    elseif d then return ca(a,b,c,d,e,f)
    elseif a then return ca(a,b,c)
    else return ca end
end
m = getmetatable(vec3())
m["is_finite"] = function(v)
    if v.x < math.huge and v.x > -math.huge and v.y < math.huge and v.y > -math.huge and v.z < math.huge and v.z > -math.huge then return true end
    return false
end
m["toQuaternion"] = function (v) return vec4(0,v.x,v.y,v.z) end
m["applyQuaternion"] = function (v,q) return q:__mul(v:toQuaternion()):__mul(q:conjugate()):vector() end
m["rotate"] = function(v,q,x,y,z)
    if is(q,"number") then q = qRotation(q,x,y,z) end
    return v:applyQuaternion(q)
end
m["__pow"] = function (v,q)
    if is(q,vec4) then return v:applyQuaternion(q) end
    return false
end
m["__concat"] = function (u,s)
    if is(s,"string") then return u:__tostring() .. s
    else return u .. s:__tostring() end
end
m["rotateTo"] = function (u,v)
    if v:lenSqr() == 0 or u:lenSqr() == 0 then return vec4(1,0,0,0) end
    u = u:normalise()
    v = u + v:normalise()
    if v:lenSqr() == 0 then
        local a,b,c = abs(u.x), abs(u.y), 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:normalise()
    local d = u:dot(v)
    u = u:cross(v)
    return vec4(d,u.x,u.y,u.z)
end
m["normalise"] = function (v)
    v = v:normalize()
    if v:is_finite() then return v
    else return vec3(0,0,1) end
end
mw,af,sw = m["__mul"],m["__add"],m["__sub"]
m["__mul"] = function(m,v)
    if is(m,vec3) and is(v,"number") then return mw(m,v) end
    if is(m,"number") and is(v,vec3) then return mw(m,v) end
    if is(m,vec3) and is(v,vec3) then return vec3(m.x*v.x,m.y*v.y,m.z*v.z) end
    if is(m,matrix) and is(v,vec3) then
        local l = m[13]*v.x+m[14]*v.y+m[15]*v.z+m[16]
        return vec3((m[1]*v.x + m[2]*v.y + m[3]*v.z + m[4])/l,(m[5]*v.x + m[6]*v.y + m[7]*v.z + m[8])/l,(m[9]*v.x + m[10]*v.y + m[11]*v.z + m[12])/l)
    end
    if is(m,vec3) and is(v,matrix) then
        local l = v[4]*m.x+v[8]*m.y+v[12]*m.z+v[16]
        return vec3((v[1]*m.x + v[5]*m.y + v[9]*m.z + v[13])/l,(v[2]*m.x + v[6]*m.y + v[10]*m.z + v[14])/l,(v[3]*m.x + v[7]*m.y + v[11]*m.z + v[15])/l)
    end
end
m["__add"] = function(a,b)
    if is(a,"number") then a = vec3(a,a,a) end
    if is(b,"number") then b = vec3(b,b,b) end
    return af(a,b)
end
m["__sub"] = function(a,b)
    if is(a,"number") then a = vec3(a,a,a) end
    if is(b,"number") then b = vec3(b,b,b) end
    return sw(a,b)
end
m["exp"] = qTangent
m = getmetatable(matrix())
mx,mr = m["__mul"],m["rotate"]
m["__mul"] = function (m,mm)
    if is(m,matrix) and is(mm,matrix) then return mx(m,mm) end
    if is(m,matrix) and is(mm,vec4) then return mx(m,mm:tomatrix()) end
    if is(m,vec4) and is(mm,matrix) then return mx(m:tomatrix(),mm) end
    if is(m,matrix) and is(mm,vec3) then
        local l = m[13]*mm.x + m[14]*mm.y + m[15]*mm.z + m[16]
        return vec3((m[1]*mm.x + m[2]*mm.y + m[3]*mm.z + m[4])/l,(m[5]*mm.x + m[6]*mm.y + m[7]*mm.z + m[8])/l,(m[9]*mm.x + m[10]*mm.y + m[11]*mm.z + m[12])/l)
    end
    if is(m,vec3) and is(mm,matrix) then
        local l = mm[4]*m.x + mm[8]*m.y + mm[12]*m.z + mm[16]
        return vec3((mm[1]*m.x + mm[5]*m.y + mm[9]*m.z + mm[13])/l,(mm[2]*m.x + mm[6]*m.y + mm[10]*m.z + mm[14])/l,(mm[3]*m.x + mm[7]*m.y + mm[11]*m.z + mm[15])/l)
    end
end
m["rotate"] = function(m,a,x,y,z)
    if is(a,vec4) then
        a,x = a:toangleaxis()
        x,y,z = x.x,x.y,x.z
    end
    return mr(m,a,x,y,z)
end