Rotation demos v4

gistfile1.txt

--# Notes
--[[

    ROTATION WITH QUATERNIONS AND FLYING

This project demonstrates the use of quaternions to rotate objects and fly them

The Main tab only contains code to help you choose from a number of demo projects. You will not need it in
your own projects.

Each numbered tab is a separate demo project, becoming more advanced from left to right. The first line of code is used by Main to manage your choice of tab, and can be deleted if you copy the code elsewhere.

The other tabs (to the right of the numbered tabs) contain the rotation library code. They are as follows:

Quaternion - is the engine, with all the quaternion functions. This is the only essential code you need, but because the use of most of the functions is not obvious, there are other "helper" classes like Rotation and Flight, which are used in the demos. 

Rotation - has a simple class to help with rotations

Flight - has specific functionality for flying. If you use this, you don't need Rotations.

Plane - builds a simple box plane, also contains a box class that is used for some of the demos in this project

Model - imports and draws full 3D models in OBJ format

Sky - wraps an image on a sphere so you can look in all directions. As used in the flight projects, the sphere is centred on the aeroplane position but not rotated, so the plane is inside the sphere (not that you can tell) and can never fly outside it.

Joystick - provides controls for both an on screen joystick and iPad tilting

You need to download the sky image from here 
https://www.dropbox.com/s/48eykfqdzxlfcjn/SkyDome1.jpg
and save to Documents:SkyDome (or save it somewhere else and change the name throughout the project)

--]]
--# Main
--Main

--This code manages which Code tab is run
--it remembers your last choice, and if you select a different one, it runs that instead

local tabs = {}
local fnames = {"setup","draw","touched","collide","orientationChanged","close","restart","keyboard","cleanup"}
local fns = {}   
local tabDesc={}

function setup()
    for k,v in ipairs(fnames) do --store addresses of key event functions
        fns[v] = _G[v] 
    end 
    LastCode=readProjectData("RotationDemo") or 1  --load stored tab number
    Choose_a_Demo=LastCode
    RunCode() 
end

function ShowList()
    output.clear()
    print("Select a demo (drag this section down out of the way, to see the selection slider below the other parameters)")  
    for i=1,#tabs do
        print(i,tabDesc[i])
    end
end

--these two functions do all the tab switching magic (thanks to Andrew_Stacey)
function localise(n,d)
    if d then tabDesc[n]=d end
    local t= {}
    setmetatable(t,{__index = _G})
    setfenv(2,t)
    tabs[n] = t
end

--change tabs
function RunCode()
    output.clear()
    saveProjectData("RotationDemo",Choose_a_Demo)
    cleanup()
    local t = tabs[Choose_a_Demo]
    for k,v in ipairs(fnames) do
        if t[v] then _G[v] = t[v] else _G[v] = fns[v] end -- overwrite with the new code
    end
    setup()
    LastCode=readProjectData("RotationDemo") or 1  --load stored tab number
    parameter.integer("Choose_a_Demo",1,#tabs,LastCode,ShowList) --tab selector
    parameter.action("Run selected demo", RunCode)
end

--default empty function to avoid errors for tabs that don't have it
function cleanup()
    parameter.clear()
end

--# 1_Basic
--1. Basic rotation
if localise then localise(1,"Basic rotation") end

function setup()
    --create a block to play with
    --this uses the Block class in the Plane tab 
    b=Block(10,15,30,color(255),vec3(0,0,0),"Platformer Art:Block Brick",{0.1,0.1,0.9,0.9})
    r=Rotation() --initialise rotations
    --these parameters let you play with rotations
    --these angles will be added to the x,y,z rotation at each draw
    parameter.number("X",0,1,.1) 
    parameter.number("Y",0,1,.1) 
    parameter.number("Z",0,1,.1) 
    parameter.action("Reset_Rotation",Reset) --restart from a level position
end

function Reset()
    r=Rotation()
end

function draw()
    background(220)
    fill(0,0,0,150)
    DrawText()
    perspective()
    camera(0,0,0,0,0,-100)
    r:turn(X,Y,Z) --apply and store additional rotation
    pushMatrix()
    translate(0,0,-50)
    r:rotate() --rotate (special version of the function, requires the Quaternion tab
    b:draw()
    popMatrix()
end

function DrawText()
    textSize(12)
    textMode(CORNER)
    text("Play with the XYZ settings at left to rotate the block.",50,HEIGHT-50)
    text("The XYZ values are in degrees, and are added to the existing rotation when redrawing",50,HEIGHT-75)
    text("Drag the Output area down to see all the parameters",50,HEIGHT-100)
end

--# 2_Moving
--2. Moving
if localise then localise(2,"Moving and rotating") end

-- All changes are marked with ** 2

function setup()
    b=Block(10,15,30,color(255),vec3(0,0,0),"Platformer Art:Block Brick",{0.1,0.1,0.9,0.9})
    r=Rotation() --initialise rotations
    pos=vec3(0,-50,-200) --initial position of block --** 2
    parameter.number("Speed",0,100,10) --speed of movement (pixels/sec)  --** 2
    parameter.number("X",-1,1,.1) -- now you can move positive or negative ** 2
    parameter.number("Y",-1,1,0.05) 
    parameter.number("Z",-1,1,0) 
    parameter.action("Reset",Reset)
end

function Reset()
    r=Rotation()
end

function draw()
    background(220)
    DrawText()
    perspective()
    camera(0,0,0,0,0,-100)
    pushMatrix()
    r:turn(X,Y,Z) --apply and store additional rotation
    --update position of block. You can provide just the z distance as shown below  --** 2
    --or if you prefer, a vec3 giving a position in x,y and z
    pos=pos+r:Move(-Speed*DeltaTime) --use negative as z axis is reversed  --** 2
    translate(pos.x,pos.y,pos.z)
    r:rotate() --rotate (special version of the function, requires the Quaternion tab)
    b:draw()
    popMatrix()
end

function DrawText()
    textSize(12)
    textMode(CORNER)
    text("Now the block moves in the direction it is facing.",50,HEIGHT-50)
    text("You can change the speed as well as XYZ direction",50,HEIGHT-75)
    text("Drag the Output area down to see all the parameters",50,HEIGHT-100)
end

--# 3_Following
--3. Moving
if localise then localise(3,"Following") end

-- All changes are marked with ** 3

function setup()
    b=Block(10,15,30,color(255),vec3(0,0,0),"Platformer Art:Block Brick",{0.1,0.1,0.9,0.9})
    r=Rotation() --initialise rotations
    pos=vec3(0,-50,-200) --initial position of block
    --we'll make another block follow the first one around, sitting on top of it and turning with it ** 3
    --it doesn't need its own Rotation because it will follow the first block
    follow=Block(5,5,5,color(255),vec3(0,0,0),"Platformer Art:Block Special") -- ** 3
    followPos=vec3(0,10,0) --follower is 20 pixels above block         ** 3
    
    parameter.number("Speed",0,100,10) --speed of movement (pixels/sec) 
    parameter.number("X",-1,1,.1) 
    parameter.number("Y",-1,1,0.05) 
    parameter.number("Z",-1,1,0) 
    parameter.action("Reset",Reset)    
end

function Reset()
    r=Rotation()
end

function draw()
    background(220)
    DrawText()
    perspective()
    camera(0,0,0,0,0,-100)
    r:turn(X,Y,Z) --apply and store additional rotation
    --update position of block. You can provide just the z distance as shown below  
    --or if you prefer, a vec3 giving a position in x,y and z
    pos=pos+r:Move(-Speed*DeltaTime) --use negative as z axis is reversed
    local f=pos + r:Move(followPos) --calculate position of follower  ** 3
    pushMatrix()
    translate(pos.x,pos.y,pos.z)
    r:rotate() --rotate (special version of the function, requires the Quaternion tab)
    b:draw()
    popMatrix()
    -- now rotate and position the following block   ** 3
    pushMatrix()
    translate(f.x,f.y,f.z)
    r:rotate() --rotate the same way as the first block
    follow:draw()
    popMatrix()
end

function DrawText()
    textSize(12)
    textMode(CORNER)
    text("A second block is positioned on top of the first one",50,HEIGHT-50)
    text("This shows how to get different objects to move together",50,HEIGHT-75)
    text("Drag the Output area down to see all the parameters",50,HEIGHT-100)
end

--# 4_CameraTracking
--4. Camera tracking
if localise then localise(4,"Camera tracking") end

-- All changes are marked with ** 4

function setup()
    b=Block(10,15,30,color(255),vec3(0,0,0),"Platformer Art:Block Brick",{0.1,0.1,0.9,0.9})
    r=Rotation() --initialise rotations
    pos=vec3(0,-50,-200) --initial position of block
    Speed=20 --fix speed for this demo       ** 4
    parameter.boolean("CameraTrackingOn",false)
    parameter.integer("CameraSetting",1,3,1) --** 4
    --three cam settings below, for each of them, x,y,z give camera position relative to the centre of the block
    camSettings={vec3(0,0,100),vec3(0,0,100),vec3(0,10,0)} -- ** 4
    parameter.number("X",-1,1,.3) 
    parameter.number("Y",-1,1,0.05) 
    parameter.number("Z",-1,1,0.05) 
    parameter.action("Reset",Reset)
end

function Reset()
    r=Rotation()
end

function draw()
    background(220)
    DrawText()
    perspective()
    r:turn(X,Y,Z) --apply and store additional rotation
    --update position of block. 
    pos=pos+r:Move(-Speed*DeltaTime) 
    
    -- now rotate and position the camera    ** 4
    --this needs to be done before we draw anything
    local camPos,camLook,camUp
    if CameraTrackingOn==false then
        camPos=vec3(0,0,0)
        camLook=vec3(0,0,-100)
        camUp=vec3(0,1,0)
    else
        camPos=pos + r:Move(camSettings[CameraSetting]) --calculate position of camera  ** 4
        --calculate where camera is looking
        --if the camera is behind the plane, we can just look at the plane position
        --but if the camera is on top of the plane, we won't see much if we look at it!
        --so this code looks ahead. This means projecting forward in the current direction
        camLook =pos+r:Move(-100) --look forward 100 pixels ** 4
        if CameraSetting>1 then --rotate camera for settings 2 & 3
            camUp=r:Move(vec3(0,1,0))  --rotate up vector
        else  --don't rotate for setting 1
            camUp=vec3(0,1,0) --default upright position
        end
    end
    camera(camPos.x,camPos.y,camPos.z,camLook.x,camLook.y,camLook.z,camUp.x,camUp.y,camUp.z)
    pushMatrix()
    translate(pos.x,pos.y,pos.z)
    r:rotate() --rotate (special version of the function, requires the Quaternion tab)
    b:draw()
    popMatrix()

end

function DrawText()
    textSize(12)
    textMode(CORNER)
    text("Now the camera tracks the object (if you turn it on, alongside)",50,HEIGHT-50)
    text("Setting 1 = camera follows behind, upright (block rotates left/right)",50,HEIGHT-70)
    text("Setting 2 = camera follows behind, rotated with the block (block stays upright)",50,HEIGHT-90)
    text("Setting 3 = camera is on top of the block (POV)",50,HEIGHT-110)
    text("It's hard to see what's happening without a skyline - so that's the next demo",50,HEIGHT-130)
    text("Drag the Output area down to see all the parameters",50,HEIGHT-150)
end

--# 5_Sky
--5. Sky sphere
if localise then localise(5,"Sky") end

-- All changes are marked with ** 5

--This only requires two additional lines of code (and a suitable image)
--see the notes in the Sky tab for how this works

function setup()
    b=Plane() -- use a simple plane  -- ** 5
    r=Rotation() --initialise rotations
    pos=vec3(0,-50,-200) --initial position of block
    Speed=20 --fix speed for this demo      
    parameter.boolean("CameraTrackingOn",false)
    parameter.integer("CameraSetting",1,3,1) 
    --three cam settings below, for each of them, x,y,z give camera position relative to the centre of the block
    camSettings={vec3(0,0,100),vec3(0,40,200),vec3(0,3,-9)} 
    parameter.number("X",-1,1,0) 
    parameter.number("Y",-1,1,0.05) 
    parameter.number("Z",-1,1,0) 
    parameter.action("Reset",Reset)
    --set up sky sphere
    sky=SkyGlobe("Documents:SkyDome","Documents:SkyDome")  -- ** 5
end

function Reset()
    r=Rotation()
end

function draw()
    background(220)
    DrawText()
    perspective()
    r:turn(X,Y,Z) --apply and store additional rotation
    --update position of block. 
    pos=pos+r:Move(-Speed*DeltaTime) 
    
    -- now rotate and position the camera   
    --this needs to be done before we draw anything
    local camPos,camLook,camUp
    if CameraTrackingOn==false then
        camPos=vec3(0,0,0)
        camLook=vec3(0,0,-100)
        camUp=vec3(0,1,0)
    else
        camPos=pos + r:Move(camSettings[CameraSetting]) --calculate position of camera  
        --calculate where camera is looking
        --if the camera is behind the plane, we can just look at the plane position
        --but if the camera is on top of the plane, we won't see much if we look at it!
        --so this code looks ahead. This means projecting forward in the current direction
        camLook =pos+r:Move(-100) --look forward 100 pixels 
        if CameraSetting>1 then --rotate camera for settings 2 & 3
            camUp=r:Move(vec3(0,1,0))  --rotate up vector
        else  --don't rotate for setting 1
            camUp=vec3(0,1,0) --default upright position
        end
    end
    camera(camPos.x,camPos.y,camPos.z,camLook.x,camLook.y,camLook.z,camUp.x,camUp.y,camUp.z)
    pushMatrix()
    translate(camPos.x,camPos.y,camPos.z)
    sky:draw()  -- ** 5
    popMatrix()
    pushMatrix()
    translate(pos.x,pos.y,pos.z)
    r:rotate() --rotate (special version of the function, requires the Quaternion tab)
    b:draw()
    popMatrix()

end

function DrawText()
    textSize(12)
    textMode(CORNER)
    text("The sky makes it easier to see what is happening (If you don't see sky, re-start the program)",50,HEIGHT-50)
    text("And now we are using a simple boxy plane object",50,HEIGHT-70)
    text("You'll find it pretty hard to make the plane behave, though...",50,HEIGHT-90)
    text("Drag the Output area down to see all the parameters",50,HEIGHT-110)
end

--# 6_Flight
--6. Introducing flight
if localise then localise(6,"Flying") end

--There are quite a few changes because the Flight functions are organised differently
-- All changes are marked with ** 6

function setup()
    b=Plane()
    f=Flight() --initialise flight instance                -- ** 6
    pos=vec3(0,-50,-200) --initial position of block
    Speed=20 --fix speed for this demo      
    parameter.number("Pitch_Limit",30,90,90,AdjustSettings)
    parameter.number("Yaw_Sensitivity",0.1,1,0.5,AdjustSettings)
    parameter.number("Roll_Limit",30,120,90,AdjustSettings)
    parameter.integer("CameraSetting",1,3,1) 
    --three cam settings below, for each of them, x,y,z give camera position relative to the centre of the plane
    camSettings={vec3(0,0,50),vec3(-10,20,100),vec3(0,3,-9)} 
    parameter.action("Reset",Reset)
    --set up sky sphere
    sky=SkyGlobe("Documents:SkyDome","Documents:SkyDome") 
    --set up joystick (see Joystick tab for options)      
    j=JoyStick()                                           -- ** 6
end

function Reset()
    r=Flight()
end

function AdjustSettings()  --this is only here so you can play with settings interactively --** 6 
    f.limit=vec3(Pitch_Limit,0,Roll_Limit)*math.pi/180  --these assume default settings apply!
    f.sens=vec3(0,Yaw_Sensitivity,0)
end

function draw()
    background(255)
    DrawText()
    perspective()
    local v=j:update()  --get joystick movement   ** 6
    if j.touch then f:turn(v) end   -- store additional rotation      ** 6
    --update position of plane 
    pos=pos+f:position(-Speed*DeltaTime)     -- ** 6
    
    -- now rotate and position the camera   
    --Flight class does it for us, send it the plane position and camera offset vector
    f:setCamera(pos,camSettings[CameraSetting])    -- ** 6
    
    pushMatrix()
    translate(pos.x,pos.y,pos.z)
    sky:draw() --we centred the sky sphere on the camera before, do it on the plane for convenience ** 6
    f:rotate() 
    b:draw()
    popMatrix() 
    j:draw()           -- draw joystick    ** 6
end

 
--manage joystick          -- pass touch values to joystick class ** 6
function touched(t)
    j:touched(t)
end

function DrawText()
    textSize(12)
    textMode(CORNER)
    text("Now we'll use the Flight class",50,HEIGHT-50)
    text("Use the on screen joystick to fly the plane",50,HEIGHT-70)
    text("Try changing the limits and sensitivity settings",50,HEIGHT-90)
end

--# 7_3DModels
--7. Importing 3D models
if localise then localise(7,"3D models") end

--all changes marked ** 7
--you can use different vehicles, see the list in the Model tab

--there can be quite a lot of work involved in using 3D models. You should get rid of any vertices you
--don't need, such as cockpit interiors. In this demo, this is done in SetupPlane.

--For the spitfire, I had to delete the static propellor and draw a spinning propellor myself.
--I left this out of this demo to save you downloading it, but if you want to include it
--it is here: https://www.dropbox.com/s/4ehrkhzeyi5fgb5/Prop2.png
--save to Dropbox:Prop2, then uncomment the line that loads the image, in setup

--note also that when we draw the plane, we can choose which of its meshes to draw
--so when we are in the cockpit, we only want the mesh named "Front", we don't need wings, tail etc

function setup()
    SetupPlane()
    --imgProp=readImage("Dropbox:Prop2") --extra image needed for propellor
    f=Flight() --initialise flight instance               
    pos=vec3(0,-50,-200) --initial position of block
    Speed=20 --fix speed for this demo      
    parameter.integer("CameraSetting",1,3,1) 
    --three cam settings below, for each of them, x,y,z give camera position relative to the centre of the plane
    camSettings={vec3(0,0,25),vec3(0,8,25),vec3(0,0.6,1.5)} 
    parameter.action("Reset",Reset)
    --set up sky sphere
    sky=SkyGlobe("Documents:SkyDome","Documents:SkyDome") 
    --set up joystick (see Joystick tab for options)      
    j=JoyStick()                                       
end

function SetupPlane()   -- ** 7 all new
    plane=OBJ("plane",1,"http://goo.gl/yw0i8m")    
    local m=plane.m --set of meshes
    --remove the meshes we don't need
    for i=#m,1,-1 do
        if m[i].name=="Propellor" then table.remove(m,i) 
        elseif m[i].name=="Cockpit" then table.remove(m,i) 
        elseif m[i].name=="Pilot" then table.remove(m,i)
        elseif m[i].name=="propeller_rotation" then table.remove(m,i) end
    end
end

function Reset()
    r=Flight()
end

function draw()
    background(255)
    DrawText()
    perspective()
    local v=j:update()  --get joystick movement   
    if j.touch then f:turn(v) end   -- store additional rotation    
    --update position of plane 
    pos=pos+f:position(-Speed*DeltaTime)     
    
    -- now rotate and position the camera   
    --Flight class does it for us, send it the plane position and camera offset vector
    f:setCamera(pos,camSettings[CameraSetting])   
    
    pushMatrix()
    translate(pos.x,pos.y,pos.z)
    sky:draw() --we centred the sky sphere on the camera before, do it on the plane for convenience 
    f:rotate() 
    --if view is inside cockpit, only draw the front of the plane
    if CameraSetting==3 then plane:draw("Body_Front") else plane:draw("All") end -- ** 7
    plane:draw()
    pushStyle() 
    --draw prop separately if we have loaded it --just for the spitfire
    if imgProp then
        translate(0,0,-2.12)
        spriteMode(CENTER)
        sprite(imgProp,0,0,3)
    end 
    popStyle() 
    popMatrix()
    j:draw()           -- draw joystick   
end

 
--manage joystick          -- pass touch values to joystick class 
function touched(t)
    j:touched(t)
end

function DrawText()
    textSize(12)
    textMode(CORNER)
    text("Now we'll use an imported 3D model",50,HEIGHT-50)
end

--# 8_Tilting
--8. Controlling the plane using tilting
if localise then localise(8,"Using tilt") end

--all changes marked ** 8

--This demo uses iPad tilt to steer the plane
--it's extremely simple
--the most important change is setting all three directions to be incremental (see post at top of Fight for more)
--to give us finer control over how the pane reacts to tilting

function setup()
    SetupPlane()
    --imgProp=readImage("Dropbox:Prop2") --extra image needed for propellor
    f=Flight({mode=vec3(1,1,1),sens=vec3(0.5,0.5,0.5)}) --initialise flight instance  ** 8             
    pos=vec3(0,-50,-200) --initial position of block
    Speed=20 --fix speed for this demo      
    parameter.integer("CameraSetting",1,3,1) 
    --three cam settings below, for each of them, x,y,z give camera position relative to the centre of the plane
    camSettings={vec3(0,0,25),vec3(0,8,25),vec3(0,0.6,1.5)} 
    parameter.action("Reset",Reset)
    --set up sky sphere
    sky=SkyGlobe("Documents:SkyDome","Documents:SkyDome")                                         
end

function SetupPlane()  
    plane=OBJ("plane",1,"http://goo.gl/yw0i8m")      
    local m=plane.m --set of meshes
    --remove the meshes we don't need
    for i=#m,1,-1 do
        if m[i].name=="Propellor" then table.remove(m,i) 
        elseif m[i].name=="Cockpit" then table.remove(m,i) 
        elseif m[i].name=="Pilot" then table.remove(m,i)
        elseif m[i].name=="propeller_rotation" then table.remove(m,i) end
    end
end

function Reset()
    r=Flight()
end

function draw()
    background(255)
    DrawText()
    perspective()
    local v=vec2(Gravity.x,Gravity.y)  --get tilt movement   ** 8
    f:turn(v)   -- store additional rotation    
    --update position of plane 
    pos=pos+f:position(-Speed*DeltaTime)     
    
    -- now rotate and position the camera   
    --Flight class does it for us, send it the plane position and camera offset vector
    f:setCamera(pos,camSettings[CameraSetting])   
    
    pushMatrix()
    translate(pos.x,pos.y,pos.z)
    sky:draw() --we centred the sky sphere on the camera before, do it on the plane for convenience 
    f:rotate() 
    --if view is inside cockpit, only draw the front of the plane
    if CameraSetting==3 then plane:draw("Body_Front") else plane:draw("All") end -- ** 7
    plane:draw()
    pushStyle() 
    --draw prop separately if we have loaded it --just for the spitfire
    if imgProp then
        translate(0,0,-2.12)
        spriteMode(CENTER)
        sprite(imgProp,0,0,3)
    end 
    popStyle() 
    popMatrix() 
end

 
--manage joystick          -- pass touch values to joystick class 
function touched(t)
    j:touched(t)
end

function DrawText()
    textSize(12)
    textMode(CORNER)
    text("Now we'll use use iPad tilt to fy",50,HEIGHT-50)
end

--# 9_Dogfight
--9. Dogfight
if localise then localise(9,"Dogfight") end

--all changes marked ** 9

--we add another plane
--and a radar screen

function setup()
    flight=Flight() --initialise flight instance   
    SetupPlane()
    SetupEnemy()                                                           --** 9         
    parameter.integer("CameraSetting",1,3,1) 
    --three cam settings below, for each of them, x,y,z give camera position relative to the centre of the plane
    camSettings={vec3(0,0,25),vec3(0,8,25),vec3(0,0.6,1.5)} 
    parameter.action("Reset",Reset)
    --set up sky sphere
    sky=SkyGlobe("Documents:SkyDome","Documents:SkyDome") 
    --set up joystick (see Joystick tab for options)      
    j=JoyStick()  
    FPS=60                                        
end

function SetupPlane()  
    flight.pos=vec3(0,0,0) --initial position of plane
    flight.speed=30 --fix speed for this demo      
    plane=OBJ("plane",1,"http://goo.gl/yw0i8m")    
    --remove the meshes we don't need
    local m=plane.m --set of meshes
    for i=#m,1,-1 do
        if m[i].name=="Propellor" then table.remove(m,i) 
        elseif m[i].name=="Cockpit" then table.remove(m,i) 
        elseif m[i].name=="Pilot" then table.remove(m,i)
        elseif m[i].name=="propeller_rotation" then table.remove(m,i) end
    end 
    radar=Radar(flight) --set up radar screen centred on this plane  ** 9
end

function SetupEnemy()  -- all new - note adjustments to default settings ** 9
    ePlane=plane --use the same plane
    eFlight=Flight() --set up flight instance to control rotation
    eFlight.mode=vec3(1,1,1) --make all adjustments incremental
    eFlight.sens=vec3(1,1,1) --so any rotations we do are not adjusted (default is 0.25)
    eFlight.level=false --plane doesn't try to level itself out
    eFlight.pos=vec3(-100,10,-500)   --start out in front of us
    eFlight.speed=15 --pixels per sec
    eFlight:turn(vec3(0,-45,0),1) --turn sideways
    radar:AddTarget(1,eFlight) --add to radar table, pass id and ref to position
end

function Reset()
    r=Flight()
end

function draw()
    FPS=FPS*.9+.1/DeltaTime
    background(255)
    DrawText()
    perspective()
    local v=j:update()  --get joystick movement   
    if j.touch then flight:turn(v) end   -- store additional rotation    
    --update position of plane 
    flight.pos=flight.pos+flight:position(-flight.speed*DeltaTime)     
    
    -- now rotate and position the camera   
    --Flight class does it for us, send it the plane position and camera offset vector
    flight:setCamera(flight.pos,camSettings[CameraSetting])   
    
    --draw enemy plane first because it's usually behind our plane            ** 9
    pushMatrix()
    --update the enemy plane rotation here if you want to
    --eFlight:turn(vec3(pitch,yaw,roll))  --in degrees                   ** 9
    eFlight.pos=eFlight.pos+eFlight:position(-eFlight.speed*DeltaTime)  --** 9
    translate(eFlight.pos.x,eFlight.pos.y,eFlight.pos.z)               --** 9
    eFlight:rotate()                                                   --** 9
    ePlane:draw("All")                                                 --** 9
    popMatrix()    
    pushMatrix()
    translate(flight.pos.x,flight.pos.y,flight.pos.z)
    sky:draw() --we centred the sky sphere on the camera before, do it on the plane for convenience 
    flight:rotate() 
    --if view is inside cockpit, only draw the front of the plane
    plane:draw("All")
    popMatrix()
    j:draw()   -- joystick  
    radar:draw()       -- draw radar ** 9
end

 
--manage joystick          -- pass touch values to joystick class 
function touched(t)
    j:touched(t)
end

function DrawText()
    textSize(12)
    textMode(CORNER)
    text("We add another plane - try to get close to it",50,HEIGHT-50)
    text(FPS,50,HEIGHT-70)
end

--# Quaternion
-- Quaternion
-- 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") or is(n,"string") 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
--# Rotation
--rotation

Rotation=class()

function Rotation:init() --create one for each object you want to rotate
   self.q=vec4(1,0,0,0) --initialise quaternion
end

--adds to rotation (degrees), either pass through x,y,z or vec3(x,y,z)
function Rotation:turn(v,y,z) 
    if y then v=vec3(v,y,z) end
    self.q=self.q*qTangent(v*math.pi/180)
end

--actually rotates the current scene
--can be used for other objects than the object you attached this to, if you want them to rotate the same way
function Rotation:rotate()
    rotate(self.q)  --special rotate function, requires Quaternion tab
end

--rotates the position of a movement vector
--in plain English, if you want to move 5 pixels forward, this will work out the actual x,y,z change
--taking into account any rotation. You then need to add it to the current position, to get the new position.
--the parameter v can be a vec3 or just the z value
function Rotation:Move(v)
    if type(v)=="number" then v=vec3(0,0,v) end
    return v^self.q  --return rotated value of v
end


--# Flight
--Flight class

--See this post for an explanation of this class
--http://coolcodea.wordpress.com/2014/02/05/149-how-the-3d-flight-library-works/
 
Flight=class()
 
--Note all the options you can set below, by passing through a table. 
--Normally the defaults are fine, so you can omit the tabe
--However, for the explanation here, we'll assume you created a table T

--You can set the initial values of yaw, pitch and roll (although usually they will start out = 0)
--If you set yaw, it needs to be T.yaw = vec3(0,Y,0) where Y is the angle you want
--Similarly if you set pitch and/or roll, it is T.pitchroll = vec3(P,0,R) where P=pitch angel and R=roll angle

--mode tells the class whether to use incremental or fixed adjustments (see post above for details)
-- T.mode = vec3(P,Y,R) with values for P, Y and R of 0 for fixed, or 1 for incremental 

--You can set the sensitivity of the controls for incremental adjustments, with T.sens=vec3(P,Y,R) and the 
--values can be anything you want. The default is 1

--You can set the rotation limits for fixed adjustments, with T.limit=vec3(P,Y,R), and the
--values are in degrees. The defaults are 90, 180, 45

--When the plane returns to level flight automatically, you can set how long this takes with T.steady=S (seconds

function Flight:init(t) 
    t = t or {}
    --yaw is kept separate from pitch and roll
    --don't mess with this unless you REALLY know what you are doing
    --self.yaw = t.yaw or vec3(0,1,0)
    --self.pitchroll = t.pitchroll or vec3(1,0,1)
    self.mode= t.mode or vec3(0,1,0)
    self.sens=t.sens or vec3(1,1,1)*0.25 --sensitivity to incremental changes
    self.limit=t.limit or vec3(90,180,90)*math.pi/180 --rotation limit
    self.level=t.level or true --true if you want the plane to rotate back to level flight when controls released
    self.steady=t.steady or 5
    self.pos=t.pos or vec3(0,0,0)  --starting position
    self.z=t.z or -1
    self.speed=t.speed or 0
    self:reset()
end
 
--resets plane to level
function Flight:reset() 
    self.y=vec4(1,0,0,0)
    self.pr=vec4(1,0,0,0)
    self.dv=vec3(0,0,0)
    self.prevRot=vec3(0,0,0)
    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
--]]
 
--rotates the plane
-- r = vec2(x,y) if we are passing x,y control values in range -1 to +1 or 
-- r = vec3(P,Y,R) in degrees (fixed or incremental value)
--t is fraction to apply, default is DeltaTime
function Flight:turn(r,t)
    if not r.z then r=vec3(-r.y,-r.x,-r.x) else r=r*math.pi/180 end
    t=t or DeltaTime
    --calculate change in fixed values
    local delta = vec3(r.x - self.prevRot.x,r.y - self.prevRot.y,r.z - self.prevRot.z)
    self.prevRot=r
    --calculate both fixed and incremental results
    --use mode to apply the correct value to pitch, yaw and roll
    --incremental adjustments are also modified by the sensitivity vector and time
    self.dv = self.mode * r * self.sens * t
    self.pr = self.pr * qTangent((1 - self.mode) * delta * self.limit)
    self.y = self.y * qTangent(self.dv)
    self.turning = 1
end

function Flight:position(s) --s is distance in pixels, either a z value or vec3
    if type(s)=="number" then s=vec3(0,0,s) end
    return s^(self.y*self.pr)
end

function Flight:inversePosition(v)
    local q=self.y*self.pr
    return v^((self.y*self.pr):conjugate())
end
 
function Flight:rotate()
    if self.level~=1 then self:interpolate() end    --update interpolation when rotating
    rotate(self.y*self.pr)
end
 
function Flight:rotation()
    return self.y*self.pr
end
 
function Flight:interpolate() 
    --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
        self.prevRot=self.prevRot*self.mode
    elseif self.turning == 1 then
        self.turning = 2
    elseif self.interpolating then
        self.interpolating = self.slerpQ(ElapsedTime)
    end
end

function Flight:interpolateSlerp()
    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

--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
    --camera looks in front of plane (and not at it)
    local look = p + self:position(vec3(0,0,-500))
    camera(self.cam.at,look,up)
end



--# Controls
--Controls

--contains Joystick controls and Tilting controls 

--1. Tilting the iPad

--This is as easy as just returning the Gravity values, which range from -1 to +1 for both x and y
--You may want to make all three directions incremental [self.mode in Flight tab = (1,1,1)] so you 
--can control how much the plane reacts to tilting, and so you can do 360 degree rolls without
--dislocating your shoulders
function Tilt()
    return vec2(Gravity.x,Gravity.y)
end

--2. Joystick

JoyStick = class()
 
--Note all the options you can set below. Pass them through in a named table
function JoyStick:init(t)
    t = t or {}
    self.radius = t.radius or 100  --size of joystick on screen
    self.stick = t.stick or 30 --size of inner circle
    self.centre = t.centre or self.radius * vec2(1,1) + vec2(5,5)
    self.position = vec2(0,0) --initial position of inner circle
    self.target = vec2(0,0) --current position of inner circle (used when we interpolate movement)
    self.value = vec2(0,0)
    self.delta = vec2(0,0)
    self.mspeed = 60
    self.moving = 0
end
 
function JoyStick:draw()
    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.touch = t.id
        end
    end
    if t.id == self.touch 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.touch = 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.touch 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
    return self.position/(self.radius - self.stick)
end
 
function JoyStick:isMoving()
    return self.moving
end
 
--3. Radar

--call this in setup
--it is positioned at bottom right of screen, change if you like

Radar=class()

function Radar:init(f)  
    self.flight=f --instance of Flight for this plane
    self.range=1000 --number of pixels covered by scope radius
    self.radius=105
    self.screenPos=vec2(WIDTH-self.radius-5,self.radius+5)  --position
    self.image=self:CreateRadarImage()
    self.targets={}
end

function Radar:AddTarget(id,f)  --p is vec reference to flight instance
    self.targets[id]=f
end
function RemoveTarget(id)
    table.remove(self.targets,id)
end
 
--create radar screen
function Radar:CreateRadarImage()
    local img=image(self.radius*2,self.radius*2)
    local r=self.radius
    setContext(img)
    pushStyle()
    fill(160, 182, 191, 50)
    stroke(118, 154, 195, 100)
    strokeWidth(1) 
    ellipse(r,r,2*r)
    stroke(0,0,0,100)
    line(r,r-10,r,r+5)
    line(r-10,r,r+10,r)
    line(r-4,r-10,r+4,r-10)
    popStyle()
    setContext()
    return img
end
 
function Radar:draw()
    ortho()
    viewMatrix(matrix())
    pushStyle()
    spriteMode(CENTER)
    sprite(self.image,self.screenPos.x,self.screenPos.y)
    fill(0,0,0,100) 
    local pa=self.flight:position(vec3(0,0,-1))
    local pv=vec2(pa.x,pa.z):normalize()
    for i,o in pairs(self.targets) do
        local d=vec2(self.flight.pos.x,self.flight.pos.z):dist(vec2(o.pos.x,o.pos.z))  --distance
        if d<=self.range then  --only display if on screen
            d=d*self.radius/self.range --calculate distance scaled to radar screen
            --direction from plane to target (without taking plane orientation into account
            local ov=vec2(o.pos.x-self.flight.pos.x,o.pos.z-self.flight.pos.z):normalize()
            --angle between plane orientation and direction from plane to target
            local a=pv:angleBetween(ov)
            ellipse(self.screenPos.x+math.sin(a)*d,self.screenPos.y+math.cos(a)*d,6)
        end
    end
    popStyle()
end

--# Sky
-- 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 (ie top half), 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
--the code will uses the width of the first image to create an image with half that height

--for more info on skyglobes see here
-- http://coolcodea.wordpress.com/2013/12/05/139-3d-sky/

SkyGlobe = class()
 
function SkyGlobe:init(i1,i2) --i1 and i2 are the names of two images
    --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=1000 --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:setColors(color(255))
    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
 

--# Model
-- This class imports 3D models in the OBJ format.

--See this post for more info
--http://coolcodea.wordpress.com/2013/12/13/139-importing-full-3d-models-into-codea/

-- You should be able to import most OBJ files without alteration, except that
--1. you need to include the accompanying .MTL file at the top of the OBJ file and insert a tag [obj] on a
--blank line between the two sections
--2. if you are planning to use image textures, you need to alter the map entries in the .MTL file so they
--show the name you want to use for the image, and its internet URL so it can be downloaded
--3. the best way to understand this is to look at the spitfire file below

--Object class
O={}

function O.LoadModel(p)
    local mod={}
    for i=1,#p do
        mod[i]=OBJ(p[i][1],p[i][2],p[i][3])
    end
    return mod
end

--OBJ library

OBJ=class()
OBJ.DataPrefix="cfg_"
OBJ.imgPrefix="Documents:z3D"

function OBJ:init(name,scale,url)
    self.name=name
    self.scale=scale or 1
    self.centre=vec3(0,0,0)
    self.data=readGlobalData(OBJ.DataPrefix..name)
    if self.data then self:ProcessData()
    else http.request(url,function(d) self:DownloadData(d) end) end
end

function OBJ:DownloadData(data)    
    if data~=nil and string.find(data,"[obj]") then
        saveGlobalData(OBJ.DataPrefix..self.name,data)
        self.data=data
        self:ProcessData()
    else print("Error loading data for "..self.name) return end
end

function OBJ:ProcessData()
    self.mtl={}
    self.m={}
    local p, v, tx, t, np, n={},{},{},{},{},{}
    local s=self.data
    local mname
    local section="mtl"
    for line in s:gmatch("[^\r\n]+") do
        line=OBJ.trim(line)
        if string.find(line,"%[obj%]")~=nil then section="obj" mname=nil end                  
        --material definition section
        if section=="mtl" then
            if string.find(line,"newmtl") then 
              mname=OBJ.GetValue(line)
              --print(mname)
              self.mtl[mname]={}
            else
              local code=string.sub(line,1,2)
              if code=="Ka" then --ambient
                  self.mtl[mname].Ka=OBJ.GetColor(line)
              elseif code=="Kd" then --diffuse
                  self.mtl[mname].Kd=OBJ.GetColor(line)
              elseif code=="Ks" then --specular
                  self.mtl[mname].Ks=OBJ.GetColor(line)
              elseif code=="Ns" then --specular exponent
                  self.mtl[mname].Ns=OBJ.GetValue(line)
              elseif code=="ill" then --illumination code
                  self.mtl[mname].illum=OBJ.GetValue(line)
              elseif code=="ma" then --texture map name
                  local u=OBJ.split(OBJ.GetValue(line)," ")
                  if string.find(u[1],"%.") then
                      self.mtl[mname].map=string.sub(u[1],1,string.find(u[1],"%.")-1)
                  else
                      self.mtl[mname].map=u[1]
                  end
                  self.mtl[mname].path=u[2]
              end
            end
         
        --data section
        elseif section=="obj" then
          --read in groups of data into separate meshes
            local code=string.sub(line,1,2)
            --look for material settings, a separate mesh is used for each
            if string.find(line,"usemtl") then
                if mname then
                    local m=mesh()
                    if self.scale~=1 then for i=1,#v do v[i]=v[i]*self.scale end end
                    m.vertices=v
                    if #t>0 then m.texCoords=t end
                    if #n>0 then m.normals=n end
                    if self.mtl[mname] then m.settings=self.mtl[mname] end
                    if m.settings and m.settings.map then m.texture=OBJ.imgPrefix..m.settings.map end
                    m.name=mname
                    self.m[#self.m+1]=m
                end
                mname=OBJ.GetValue(line)
                v,t,n={},{},{}
                --print(mname)
            end
            if code=="v " then --point position
              p[#p+1]=OBJ.GetVec3(line)
            elseif code=="vn" then --point normal
              np[#np+1]=OBJ.GetVec3(line)    
            elseif code=="vt" then --texture co-ord
              tx[#tx+1]=OBJ.GetVec2(line)
            elseif code=="f " then --vertex
              local pts,ptex,pnorm=OBJ.GetList(line)
              if #pts==3 then
                for i=1,3 do v[#v+1]=p[tonumber(pts[i])] end
                if ptex then for i=1,3 do t[#t+1]=tx[tonumber(ptex[i])] end end
                if pnorm then for i=1,3 do n[#n+1]=np[tonumber(pnorm[i])] end end
              elseif #pts==4 then
                for i=1,3 do v[#v+1]=p[tonumber(pts[i])] end
                if ptex then for i=1,3 do t[#t+1]=tx[tonumber(ptex[i])] end end
                if pnorm then for i=1,3 do n[#n+1]=np[tonumber(pnorm[i])] end end
                v[#v+1]=p[tonumber(pts[3])]
                if ptex then t[#t+1]=tx[tonumber(ptex[3])] end
                if pnorm then n[#n+1]=np[tonumber(pnorm[3])] end
                v[#v+1]=p[tonumber(pts[4])]
                if ptex then t[#t+1]=tx[tonumber(ptex[4])] end
                if pnorm then n[#n+1]=np[tonumber(pnorm[4])] end
                v[#v+1]=p[tonumber(pts[1])]
                if ptex then t[#t+1]=tx[tonumber(ptex[1])] end
                if pnorm then n[#n+1]=np[tonumber(pnorm[1])] end
              elseif #pts>4 then
                local cx,cy,cz=0,0,0
                local ttx,tty=0,0
                local nx,ny,nz=0,0,0
                for i=1,#pts do
                    local u=p[tonumber(pts[i])] cx,cy,cz=cx+u.x,cy+u.y,cz+u.z
                    if ptex then local u=tx[tonumber(ptex[i])] ttx,tty=ttx+u.x,tty+u.y end
                end
                local cp=vec3(cx/#pts,cy/#pts,cz/#pts)
                if ptex then ct=vec2(ttx/#pts,tty/#pts) end
                local j
                for i=1,#pts do
                    if i<#pts then j=i+1 else j=1 end
                    v[#v+1]=p[tonumber(pts[i])]
                    if ptex then t[#t+1]=tx[tonumber(ptex[i])] end
                    v[#v+1]=p[tonumber(pts[j])]
                    if ptex then t[#t+1]=tx[tonumber(ptex[j])] end
                    v[#v+1]=cp
                    if ptex then t[#t+1]=ct end
                end
              end
            end
        end  
    end
    local m=mesh()
    if self.scale~=1 then for i=1,#v do v[i]=v[i]*self.scale end end
    m.vertices=v
    if #t>0 then m.texCoords=t end
    if #n>0 then m.normals=n end
    m.settings=self.mtl[mname]
    if m.settings and m.settings.map then m.texture=OBJ.imgPrefix..m.settings.map end
    m.name=mname
    self.m[#self.m+1]=m 
    self:GetStats()
    --download images if not stored locally
    self.MissingImages={}
    for i,O in pairs(self.mtl) do
        if O.map then
            local y=readImage(OBJ.imgPrefix..O.map)
            if not y then self.MissingImages[#self.MissingImages+1]={O.map,O.path} end
        end
    end
    if #self.MissingImages>0 then self:LoadImages() end
    propAngle=0 --propellor
end

function OBJ:LoadImages()
    --print("downloading"..self.MissingImages[1][1])
    http.request(self.MissingImages[1][2],function(d) self:StoreImage(d) end)
end

function OBJ:StoreImage(d)
    --print("saving"..self.MissingImages[1][1])
    saveImage(OBJ.imgPrefix..self.MissingImages[1][1],d)
    table.remove(self.MissingImages,1)
    if #self.MissingImages~=0 then self:LoadImages() end
end

function OBJ:DeleteData()
    saveGlobalData(OBJ.DataPrefix..self.name,nil)
    for i,O in pairs(self.mtl) do
        if O.map then
            ---print("deleting "..OBJ.imgPrefix..O.map)
            local y=saveImage(OBJ.imgPrefix..O.map,nil)
        end
    end
end

function OBJ.GetColor(s)
  local s1=string.find(s," ")
  local s2=string.find(s," ",s1+1)
  local s3=string.find(s," ",s2+1)
  return color(string.sub(s,s1+1,s2-1)*255,string.sub(s,s2+1,s3-1)*255,string.sub(s,s3+1,string.len(s))*255)
end

function OBJ.GetVec3(s)
  local s1=string.find(s," ")
  local s2=string.find(s," ",s1+1)
  local s3=string.find(s," ",s2+1)
  return vec3(math.floor(string.sub(s,s1+1,s2-1)*100)/100,
    math.floor(string.sub(s,s2+1,s3-1)*100)/100,
    math.floor(string.sub(s,s3+1,string.len(s))*100)/100)
end

function OBJ.GetVec2(s)
  local s1=string.find(s," ")
  local s2=string.find(s," ",s1+1)
  local s3=string.find(s," ",s2+1)
  if s3 then
      return vec3(math.floor(string.sub(s,s1+1,s2-1)*100)/100,
            math.floor(string.sub(s,s2+1,s3-1)*100)/100)
  else
      return vec2(math.floor(string.sub(s,s1+1,s2-1)*100)/100,
            math.floor(string.sub(s,s2+1,string.len(s))*100)/100)
  end
end

function OBJ:draw(v)
    for i=1,#self.m do
        if v=="All" or (v and string.find(v,self.m[i].name)) then
            self.m[i]:draw()
        end
    end
end

function OBJ.GetValue(s)
  return string.sub(s,string.find(s," ")+1,string.len(s))
 end
 
function OBJ.trim(s)
  while string.find(s,"  ") do s = string.gsub(s,"  "," ") end
  return s:match'^()%s*$' and '' or s:match'^%s*(.*%S)'
end

function OBJ.split(s,sep)
  sep=sep or "/"
  local p={}
  local pattern = string.format("([^%s]+)", sep)
  string.gsub(s,pattern, function(c) p[#p+1] = c end)
  return p
end

function OBJ.GetList(s)
   local p,t,n={},{},{}
   p=OBJ.split(s," ")
   table.remove(p,1)
   for i=1,#p do
      local a=OBJ.split(p[i])
      if #a==1 then
        p[i]=math.abs(a[1])
      elseif #a==2 then
        p[i]=math.abs(a[1])
        t[i]=math.abs(a[2])
      elseif #a==3 then
        p[i]=math.abs(a[1])
        t[i]=math.abs(a[2])
        n[i]=math.abs(a[3])
      end
   end
   return p,t,n
end

function OBJ:GetStats()
    local x,y,z,n=0,0,0,0
    local minx,maxx,miny,maxy,minz,maxz=999,-999,999,-999,999,-999
    for i=1,#self.m do
        local vv=self.m[i].vertices
        local nn=0
        for j=1,#vv do
            local v=vv[i]
            if v~=nil then
                x,y,z=x+v.x,y+v.y,z+v.z 
                nn=nn+1
                if v.x<minx then minx=v.x end
                if v.x>maxx then maxx=v.x end
                if v.y<miny then miny=v.y end
                if v.y>maxy then maxy=v.y end
                if v.z<minz then minz=v.z end
                if v.z>maxz then maxz=v.z end
            end
        end
        n=n+nn
    end
    self.centre=vec3(x/n,y/n,z/n)
    self.minvert=vec3(minx,miny,minz)
    self.maxvert=vec3(maxx,maxy,maxz)
    self.size=vec3(maxx-minx,maxy-miny,maxz-minz)
end

function CalculateNormals(vertices)
    --this assumes flat surfaces, and hard edges between triangles
    local norm = {}
    for i=1, #vertices,3 do --calculate normal for each set of 3 vertices
        local n = ((vertices[i+1] - vertices[i]):cross(vertices[i+2] - vertices[i])):normalize()
        norm[i] = n --then apply it to all 3
        norm[i+1] = n
        norm[i+2] = n
    end
    return norm
end   

function CalculateAverageNormals(vertices,f)
    --average normals at each vertex
    --first get a list of unique vertices, concatenate the x,y,z values as a key
    local norm,unique= {},{}
    for i=1, #vertices do
        unique[vertices[i].x ..vertices[i].y..vertices[i].z]=vec3(0,0,0)
    end
    --calculate normals, add them up for each vertex and keep count
    for i=1, #vertices,3 do --calculate normal for each set of 3 vertices
        local n = (vertices[i+1] - vertices[i]):cross(vertices[i+2] - vertices[i]) 
        for j=0,2 do
            local v=vertices[i+j].x ..vertices[i+j].y..vertices[i+j].z
            unique[v]=unique[v]+n  
        end
    end
    --calculate average for each unique vertex
    for i=1,#unique do
        unique[i] = unique[i]:normalize()
    end
    --now apply averages to list of vertices
    for i=1, #vertices,3 do --calculate average
        local n = (vertices[i+1] - vertices[i]):cross(vertices[i+2] - vertices[i]) 
        for j=0,2 do
            norm[i+j] = unique[vertices[i+j].x ..vertices[i+j].y..vertices[i+j].z]
        end
    end
    return norm 
end   

--# Plane
--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))  --fuselage 
    a[2]=Block(15,1,5,c2,vec3(-9,1,-6))  --left wing
    a[3]=Block(15,1,5,c2,vec3(9,1,-6))  --right wing
    a[4]=Block(1,4,3,c2,vec3(0,3.5,11)) --vertical tail 
    a[5]=Block(4,1,2,c2,vec3(-3.5,1,11))  --left tail
    a[6]=Block(4,1,2,c2,vec3(3.5,1,11))  --right tail
    a[7]=Block(3,2,5,c3,vec3(0,2.5,-6)) --cockpit
    a[8]=Block(3,3,1,c4,vec3(0,0,-13)) --nose
    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[1],r[2],r[3],r[4]} 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
  
    local cubetexCoords
    if self.tex then    
        local texvertices = { vec2(self.texR[1],self.texR[2]),
                              vec2(self.texR[3],self.texR[2]),
                              vec2(self.texR[1],self.texR[4]),
                              vec2(self.texR[3],self.texR[4]) }
                
        cubetexCoords = {
        texvertices[1], texvertices[2], texvertices[4],
        texvertices[1], texvertices[4], texvertices[3],
        -- Right
        texvertices[1], texvertices[2], texvertices[4],
        texvertices[1], texvertices[4], texvertices[3],
        -- Back
        texvertices[1], texvertices[2], texvertices[4],
        texvertices[1], texvertices[4], texvertices[3],
        -- Left
        texvertices[1], texvertices[2], texvertices[4],
        texvertices[1], texvertices[4], texvertices[3],
        -- Top
        texvertices[1], texvertices[2], texvertices[4],
        texvertices[1], texvertices[4], texvertices[3],
        -- Bottom
        texvertices[1], texvertices[2], texvertices[4],
        texvertices[1], texvertices[4], texvertices[3],
        }
    end
    
    --put it all together
    local ms = mesh()
    ms.vertices = cubeverts
    if self.tex then
        ms.texture = self.tex
        ms.texCoords = cubetexCoords
        ms:setColors(self.color)
    else
        ms.colors=cc
    end
    return ms
end  
    
function Block:draw()
    self.blk:draw() 
end