thermo_mechanical_subduction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## 2D thermo mechanical subduction models\n",
    "\n",
    "This notebook develops a simple, flexible model for 2d subduction using Underworld2. All functionality is parallel-compatible. The code has been run on 48 processors, at a vertical resolution of 256 elements (Q1).\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "This code was written by Dan Sandiford, Louis Moresi and the Underworld Team. It is licensed under the [Creative Commons Attribution 4.0 International License](https://creativecommons.org/licenses/by/4.0/) . We offer this licence to encourage you to modify and share the examples and use them to help you in your research.\n",
    "\n",
    "\n",
    "<hr>\n",
    "<a rel=\"license\" href=\"http://creativecommons.org/licenses/by/4.0/\"><img alt=\"Creative Commons License\" style=\"border-width:0\" src=\"https://i.creativecommons.org/l/by/4.0/88x31.png\" /></a><br />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Notes\n",
    "\n",
    "### General\n",
    "\n",
    "Many of the functions and scripts used here (various parts of the Underworld2 unsupported module) were written by Dan Sandiford during the bloodshed and mayhem of a PhD. This code has minimal safeguards. When things break (they will), post and issue, and we'll try to provide a fix!\n",
    "\n",
    "### Checkpointing\n",
    "\n",
    "All parts of the code that relate to checkpointing have the following identifier:\n",
    "\n",
    "`#*************CHECKPOINT-BLOCK**************#`\n",
    "\n",
    "### Materials\n",
    "\n",
    "Currently a weak material is added to the top of the lower plate to decouple the plates. For simplicity, in this notebook, there is no functionality for updating the distribution of materials (i.e. creating new crust). An example of one way you could add this functionality can be found here: https://github.com/dansand/materialTransformations\n",
    "\n",
    "\n",
    "### Solver\n",
    "\n",
    "We have found the penalty method works well for this problem. This is the default setting. (Unfortunately, it doesn't scale well to 3D.)\n",
    "\n",
    "### Resolution and particles\n",
    "\n",
    "Properly resolving the weak crust is critical to this model, which means that resolution (and particle numbers) need to be quite high (i.e. don't expect quick models on your laptop). The default crust thickness is 15 km. By switching the depth of the model to something like 660 km (upper mantle), you might get away with running at a resolution of ~ 96 elements (vertical). Note that the dyanamics are quite sensitive to the depth of the model. The default resolution and particle density are set low, essentially for model setup purposes. \n",
    "\n",
    "\n",
    "### nearest neighbour evaluation\n",
    "\n",
    "In a couple of places a function `nn_evaluation` (nearest neighbour evaluation) is called. This is used whenever we want to map a swarm variable to a mesh variable. This is an alternative to the `uw.utils.MeshVariable_Projection` functionality provided in underworld2.\n",
    "\n",
    "### Visualisation\n",
    "\n",
    "Examples of the both 'native' Underworld2 visualisation, through gLucifer, as well as xdmf file output are provided. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Package requirements\n",
    "\n",
    "Assuming you're running underworld2 through a docker image, you will need the following python packages:\n",
    "\n",
    "* easydict  https://pypi.python.org/pypi/easydict/\n",
    "\n",
    "It should be simple enough to do `pip install easydict` (though you'll need to do this for each container you start)\n",
    "\n",
    "You'll also need a copy of this module: https://github.com/dansand/unsupported.git\n",
    "\n",
    "This is occasionally merged with the official underworld2 module of the same name, but it's safest to grab a copy of the repo: \n",
    "\n",
    "`git clone https://github.com/dansand/unsupported.git`, \n",
    "\n",
    "(easiest to clone straight into in the same directory as this notebook.)\n",
    "\n",
    "To add the module to the pythonpath do:\n",
    "\n",
    "```python\n",
    "import os\n",
    "import sys\n",
    "\n",
    "if os.getcwd() == '/workspace/newSlab':\n",
    "    sys.path.append('./unsupported')\n",
    "    \n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "\n",
    "if '/workspace' in os.getcwd(): #a simple test for the standard Docker image\n",
    "    #sys.path.append('./unsupported') #unsupported module located in current working dir\n",
    "    sys.path.append('../unsupported') #unsupported module located in parent dir.\n",
    "\n",
    "else:\n",
    "    pass #put your non-docker path here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import underworld as uw\n",
    "import math\n",
    "from underworld import function as fn\n",
    "import glucifer\n",
    "from easydict import EasyDict as edict\n",
    "import operator\n",
    "import pickle\n",
    "\n",
    "\n",
    "#\n",
    "from unsupported_dan.utilities.interpolation import nn_evaluation\n",
    "from unsupported_dan.utilities.subduction import slab_top\n",
    "from unsupported_dan.interfaces.marker2D import markerLine2D\n",
    "from unsupported_dan.easymodels import checkpoint\n",
    "from unsupported_dan.easymodels import easy_args\n",
    "from unsupported_dan.utilities.misc import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Setup output directories"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "###########\n",
    "#Standard output directory setup\n",
    "###########\n",
    "\n",
    "\n",
    "#Model letter identifier demarker\n",
    "Model = \"T\"\n",
    "\n",
    "#Model number identifier demarker:\n",
    "ModNum = 1\n",
    "\n",
    "#Any isolated letter / integer command line args are interpreted as Model/ModelNum and read in here\n",
    "#We use this for doing parameter sweeps\n",
    "\n",
    "if len(sys.argv) == 1:\n",
    "    ModNum = ModNum \n",
    "elif sys.argv[1] == '-f': #\n",
    "    ModNum = ModNum \n",
    "else:\n",
    "    for farg in sys.argv[1:]:\n",
    "        if not '=' in farg: #then assume it's a not a parameter argument\n",
    "            try:\n",
    "                ModNum = int(farg) #try to convert everingthing to a float, else remains string\n",
    "            except ValueError:\n",
    "                Model  = farg\n",
    "                \n",
    "                \n",
    "###########\n",
    "#Standard output directory setup\n",
    "###########\n",
    "\n",
    "outputPath = \"results\" + \"/\" +  str(Model) + \"/\" + str(ModNum) + \"/\" \n",
    "imagePath = outputPath + 'images/'\n",
    "filePath = outputPath + 'files/'\n",
    "#checkpointPath = outputPath + 'checkpoint/'\n",
    "dbPath = outputPath + 'gldbs/'\n",
    "xdmfPath = outputPath + 'xdmf/'\n",
    "outputFile = 'results_model' + Model + '_' + str(ModNum) + '.dat'\n",
    "\n",
    "if uw.rank()==0:\n",
    "    # make directories if they don't exist\n",
    "    if not os.path.isdir(outputPath):\n",
    "        os.makedirs(outputPath)\n",
    "    if not os.path.isdir(imagePath):\n",
    "        os.makedirs(imagePath)\n",
    "    if not os.path.isdir(dbPath):\n",
    "        os.makedirs(dbPath)\n",
    "    if not os.path.isdir(filePath):\n",
    "        os.makedirs(filePath)\n",
    "    if not os.path.isdir(xdmfPath):\n",
    "        os.makedirs(xdmfPath)\n",
    "        \n",
    "uw.barrier() #Barrier here so no procs run the check in the next cell too early        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Checkpointing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#*************CHECKPOINT-BLOCK**************#\n",
    "\n",
    "#cp = checkpoint(outputPath + 'checkpoint/', loadpath='./results/A/1/checkpoint/10')\n",
    "cp = checkpoint(outputPath + 'checkpoint/')\n",
    "\n",
    "#*************CHECKPOINT-BLOCK**************#\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Model parameters and scaling\n",
    "\n",
    "Model versatility is controlled by 4 dictionaries, \n",
    "\n",
    "* 'dp': dimensional parameters\n",
    "* 'md': model dictionary (changing numerics, switching certain processes on / off)\n",
    "* 'sf': scaling factors\n",
    "* 'ndp': non dimensional parameters\n",
    "\n",
    "The dictionaries use the easyDict format, so that the following format can be applied to items:\n",
    "\n",
    "```\n",
    "\n",
    "key.value = blah\n",
    "\n",
    "```\n",
    "\n",
    "This approach was chosen so that:\n",
    "\n",
    "1. the non-dimensionalisation process is fully developed and explicit in each model\n",
    "2. parameters can be easy saved (dictionaries can be 'pickled', written to .csv etc)\n",
    "3. parameters can be easily altered through command line args. \n",
    "\n",
    "Command line arguments can only be provided to the `dp` and `md` dictionaries. The logic is that the scaling process (non-dimesionalisation) should basically be static, and simply provides a mapping between physical units and model units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "dp = edict({})\n",
    "#Main physical paramters\n",
    "dp.depth=1000e3                         #model Depth (width set by md.aspectRatio)\n",
    "dp.refDensity=3300.                     #reference density\n",
    "dp.refGravity=9.8                       #surface gravity\n",
    "dp.viscosityScale=1e20                  #reference upper mantle visc., \n",
    "dp.refDiffusivity=1e-6                  #thermal diffusivity\n",
    "dp.refExpansivity=3e-5                  #surface thermal expansivity\n",
    "dp.gasConstant=8.314                    #gas constant\n",
    "dp.specificHeat=1250.                   #Specific heat (Jkg-1K-1)\n",
    "dp.potentialTemp=1573.                  #mantle potential temp (K)\n",
    "dp.surfaceTemp=273.                     #surface temp (K)\n",
    "#Rheology - flow law paramters\n",
    "dp.cohesionMantle=20e6                  #mantle cohesion in Byerlee law\n",
    "dp.cohesionCrust=1e6                    #crust cohesion in Byerlee law\n",
    "dp.frictionMantle=0.2                   #mantle friction coefficient in Byerlee law (tan(phi))\n",
    "dp.frictionCrust=0.02                   #crust friction coefficient \n",
    "dp.diffusionPreExp=5.34e-10             #or 1./1.87e9, pre-exp factor for diffusion creep\n",
    "dp.diffusionEnergy=3e5 \n",
    "dp.diffusionVolume=5e-6\n",
    "dp.lowerMantlePreExp=4.23e-15           \n",
    "dp.lowerMantleEnergy=2.0e5\n",
    "dp.lowerMantleVolume=1.5e-6\n",
    "dp.lowerMantleViscFac = 30.\n",
    "#Rheology - cutoff values\n",
    "dp.viscosityMin=1e18\n",
    "dp.viscosityMax=1e25                    #viscosity max in the mantle material\n",
    "dp.viscosityMinCrust=1e20               #viscosity min in the weak-crust material\n",
    "dp.viscosityMaxCrust=1e20               #viscosity max in the weak-crust material\n",
    "dp.yieldStressMax=300*1e6              \n",
    "dp.crustViscCutoffDepth = 100e3\n",
    "dp.crustViscEndWidth = 20e3\n",
    "#Intrinsic Lengths\n",
    "dp.crustThickness = 15.*1e3             #weak layer thickness on top of slabs\n",
    "dp.crustMantleDepth=250.*1e3 \n",
    "dp.lowerMantleDepth=660.*1e3  \n",
    "dp.crustLimitDepth=650.*1e3             #Deeper than this, crust material rheology reverts to mantle rheology\n",
    "#Slab and plate init. parameters\n",
    "dp.subZoneLoc=-100e3                    #X position of subduction zone...km\n",
    "dp.leftRidge=-1.*(5000e3)               \n",
    "dp.rightRidge=(5000e3)\n",
    "dp.maxDepth=150e3\n",
    "dp.theta=40                             #Angle of slab\n",
    "dp.radiusOfCurv = 250e3                 #radius of curvature\n",
    "dp.slabMaxAge=70e6                      #age of subduction plate at trench\n",
    "dp.plateMaxAge=100e6                    #max age of slab (Plate model)\n",
    "dp.opMaxAge=35e6                        #age of op\n",
    "#Misc\n",
    "dp.stickyAirDepth=100e3                 #depth of sticky air layer\n",
    "dp.viscosityStickyAir=1e19              #stick air viscosity, normal\n",
    "#derived params\n",
    "dp.deltaTemp = dp.potentialTemp-dp.surfaceTemp\n",
    "\n",
    "\n",
    "\n",
    "#Modelling and Physics switches\n",
    "\n",
    "md = edict({})\n",
    "md.refineMeshStatic=True\n",
    "md.stickyAir=False\n",
    "md.aspectRatio=5.\n",
    "md.res=32\n",
    "md.ppc=25                               #particles per cell\n",
    "md.elementType=\"Q1/dQ0\"                 #\"Q2/DPC1\"\n",
    "md.secInvFac=math.sqrt(1.)\n",
    "md.courantFac=0.5                       #extra limitation on timestepping may need to be decreased if using \"Q2/DPC1\" )\n",
    "md.thermal = True                       #thermal system or compositional\n",
    "md.swarmInitialFac = 0.6                #initial swarm layout will be int(md.ppc*md.swarmInitialFac), popControl will densify later\n",
    "md.compBuoyancy = False\n",
    "md.nltol = 0.01\n",
    "md.maxSteps = 5000\n",
    "md.checkpointEvery = 100\n",
    "md.swarmUpdate = 10\n",
    "md.penaltyMethod = True\n",
    "md.opuniform = False\n",
    "md.spuniform = False\n",
    "md.opfixed = False\n",
    "md.spfixed = False\n",
    "md.buoyancyFac = 1.0\n",
    "#time-based actions\n",
    "md.filesFreqYears = 1.0e6 #dimensional time in years\n",
    "\n",
    "\n",
    "\n",
    "uw.barrier()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "##Parse any command-line args\n",
    "\n",
    "sysArgs = sys.argv\n",
    "\n",
    "#We want to run this on both the parameter dict, and the model dict\n",
    "easy_args(sysArgs, dp)\n",
    "easy_args(sysArgs, md)\n",
    "\n",
    "\n",
    "uw.barrier()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "sf = edict({})\n",
    "\n",
    "sf.lengthScale=2900e3\n",
    "sf.viscosityScale = dp.viscosityScale\n",
    "sf.stress = (dp.refDiffusivity*sf.viscosityScale)/sf.lengthScale**2\n",
    "#sf.lithGrad = dp.refDensity*dp.refGravity*(sf.lengthScale)**3/(sf.viscosityScale*dp.refDiffusivity) \n",
    "sf.lithGrad = (sf.viscosityScale*dp.refDiffusivity) /(dp.refDensity*dp.refGravity*(sf.lengthScale)**3)\n",
    "sf.velocity = dp.refDiffusivity/sf.lengthScale\n",
    "sf.strainRate = dp.refDiffusivity/(sf.lengthScale**2)\n",
    "sf.time = 1./sf.strainRate\n",
    "sf.actVolume = (dp.gasConstant*dp.deltaTemp)/(dp.refDensity*dp.refGravity*sf.lengthScale)\n",
    "sf.actEnergy = (dp.gasConstant*dp.deltaTemp)\n",
    "sf.diffusionPreExp = 1./sf.viscosityScale\n",
    "sf.deltaTemp  = dp.deltaTemp\n",
    "sf.pressureDepthGrad = (dp.refDensity*dp.refGravity*sf.lengthScale**3)/(dp.viscosityScale*dp.refDiffusivity)\n",
    "\n",
    "\n",
    "#dimesionless params\n",
    "ndp  = edict({})\n",
    "ndp.rayleigh = md.buoyancyFac*(dp.refExpansivity*dp.refDensity*dp.refGravity*dp.deltaTemp*sf.lengthScale**3)/(dp.viscosityScale*dp.refDiffusivity)\n",
    "ndp.dissipation = (dp.refExpansivity*sf.lengthScale*dp.refGravity)/dp.specificHeat\n",
    "ndp.surfaceTemp = dp.surfaceTemp/sf.deltaTemp  #Ts\n",
    "ndp.potentialTemp = dp.potentialTemp/sf.deltaTemp - ndp.surfaceTemp #Tp' = Tp - TS\n",
    "#lengths / distances\n",
    "ndp.depth = dp.depth/sf.lengthScale\n",
    "ndp.leftLim = -0.5*ndp.depth*md.aspectRatio\n",
    "ndp.rightLim = 0.5*ndp.depth*md.aspectRatio\n",
    "ndp.crustThickness = dp.crustThickness/sf.lengthScale\n",
    "ndp.leftRidge = max(ndp.leftLim,  dp.leftRidge/sf.lengthScale)\n",
    "ndp.rightRidge = min(ndp.rightLim, dp.rightRidge/sf.lengthScale)\n",
    "ndp.crustLimitDepth = dp.crustLimitDepth/sf.lengthScale\n",
    "ndp.lowerMantleDepth = dp.lowerMantleDepth/sf.lengthScale\n",
    "#times - for convenience and sanity the dimensional values are in years, conversion to seconds happens here\n",
    "ndp.slabMaxAge =  dp.slabMaxAge*(3600*24*365)/sf.time\n",
    "ndp.plateMaxAge =  dp.plateMaxAge*(3600*24*365)/sf.time\n",
    "ndp.opMaxAge = dp.opMaxAge*(3600*24*365)/sf.time\n",
    "#Rheology - flow law paramters\n",
    "ndp.cohesionMantle=dp.cohesionMantle/sf.stress                  #mantle cohesion in Byerlee law\n",
    "ndp.cohesionCrust=dp.cohesionCrust/sf.stress                  #crust cohesion in Byerlee law\n",
    "ndp.frictionMantle=dp.frictionMantle/sf.lithGrad                  #mantle friction coefficient in Byerlee law (tan(phi))\n",
    "ndp.frictionCrust=dp.frictionCrust/sf.lithGrad                  #crust friction coefficient \n",
    "ndp.diffusionPreExp=dp.diffusionPreExp/sf.diffusionPreExp                #pre-exp factor for diffusion creep\n",
    "ndp.diffusionEnergy=dp.diffusionEnergy/sf.actEnergy\n",
    "ndp.diffusionVolume=dp.diffusionVolume/sf.actVolume\n",
    "ndp.lowerMantlePreExp=dp.lowerMantlePreExp/sf.diffusionPreExp \n",
    "ndp.lowerMantleEnergy=dp.lowerMantleEnergy/sf.actEnergy\n",
    "ndp.lowerMantleVolume=dp.lowerMantleVolume/sf.actVolume\n",
    "ndp.yieldStressMax=dp.yieldStressMax/sf.stress \n",
    "#Rheology - cutoff values\n",
    "ndp.viscosityMin= dp.viscosityMin /sf.viscosityScale\n",
    "ndp.viscosityMax=dp.viscosityMax/sf.viscosityScale\n",
    "ndp.viscosityMinCrust= dp.viscosityMinCrust /sf.viscosityScale\n",
    "ndp.viscosityMaxCrust = dp.viscosityMaxCrust/sf.viscosityScale\n",
    "ndp.lowerMantleViscFac = dp.lowerMantleViscFac\n",
    "ndp.crustViscCutoffDepth = dp.crustViscCutoffDepth/sf.lengthScale\n",
    "ndp.crustViscEndWidth = dp.crustViscEndWidth/sf.lengthScale\n",
    "#Slab and plate init. parameters\n",
    "ndp.subZoneLoc = dp.subZoneLoc/sf.lengthScale\n",
    "ndp.maxDepth = dp.maxDepth/sf.lengthScale\n",
    "ndp.radiusOfCurv = dp.radiusOfCurv/sf.lengthScale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#*************CHECKPOINT-BLOCK**************#\n",
    "\n",
    "\n",
    "#if restart, attempt to read in saved dicts. \n",
    "if cp.restart:\n",
    "    try:\n",
    "        with open(os.path.join(cp.loadpath, 'dp.pkl'), 'rb') as fp:\n",
    "                            dp = pickle.load(fp)\n",
    "        with open(os.path.join(cp.loadpath, 'sf.pkl'), 'rb') as fp:\n",
    "                            sf = pickle.load(fp)\n",
    "        with open(os.path.join(cp.loadpath, 'md.pkl'), 'rb') as fp:\n",
    "                            md = pickle.load(fp)\n",
    "\n",
    "    except:\n",
    "        print(\"couldn't load paramter dictionaries on restart\")\n",
    "\n",
    "\n",
    "    \n",
    "#add dicts to the checkpointinng object\n",
    "cp.addDict(dp, 'dp')\n",
    "cp.addDict(sf, 'sf')\n",
    "cp.addDict(md, 'md')\n",
    "\n",
    "#*************CHECKPOINT-BLOCK**************#\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Build Mesh and FE variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#Domain and Mesh paramters\n",
    "yres = int(md.res)\n",
    "xres = int(md.res*8) \n",
    "\n",
    "\n",
    "\n",
    "mesh = uw.mesh.FeMesh_Cartesian( elementType = (md.elementType),\n",
    "                                 elementRes  = (xres, yres), \n",
    "                                 minCoord    = (ndp.leftLim, 1. - ndp.depth), \n",
    "                                 maxCoord    = (ndp.rightLim, 1.)) \n",
    "\n",
    "velocityField   = uw.mesh.MeshVariable( mesh=mesh,         nodeDofCount=2 )\n",
    "pressureField   = uw.mesh.MeshVariable( mesh=mesh.subMesh, nodeDofCount=1 )\n",
    "temperatureField    = uw.mesh.MeshVariable( mesh=mesh,         nodeDofCount=1 )\n",
    "\n",
    "if md.thermal:\n",
    "    temperatureDotField = uw.mesh.MeshVariable( mesh=mesh,         nodeDofCount=1 ) #create this only if Adv-diff\n",
    "    diffusivityFn = fn.misc.constant(1.)\n",
    "    \n",
    "    \n",
    "    \n",
    "# Any extra mesh vars. we want to define (mostly to facilite saving as xdmf)\n",
    "strainRateField    = uw.mesh.MeshVariable( mesh=mesh,         nodeDofCount=1 )\n",
    "viscosityField    = uw.mesh.MeshVariable( mesh=mesh,         nodeDofCount=1 )\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Here we provide a basic mesh refinement strategy: This refinement function is:\n",
    "\n",
    "$y_{new} = \\frac{ln((\\alpha y_n + e) -1)}{ln((\\alpha + e) - 1)}$\n",
    "\n",
    "where $y_n$ are the normalised original coordinates\n",
    "\n",
    "Note that the horizonal mesh resolution has been set to balance this vertical refinement."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "mesh.reset()\n",
    "alpha = 5.\n",
    "am = mesh.minCoord[1]\n",
    "bm = mesh.maxCoord[1]\n",
    "\n",
    "yorigscaled = (mesh.data[:,1].copy()  - am) /(bm - am)\n",
    "ynew = (np.log(alpha*yorigscaled + np.e ) -1.)/(np.log(alpha+ np.e ) -1.)\n",
    "ynew *= (bm - am)\n",
    "ynew += am\n",
    "\n",
    "#do the mesh deformation\n",
    "with mesh.deform_mesh():\n",
    "    mesh.data[:,1] = ynew\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#figMesh= glucifer.Figure(quality=3)\n",
    "#figMesh.append( glucifer.objects.Mesh(mesh))\n",
    "#figMesh.show()\n",
    "#figMesh.save_database('test.gldb')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#*************CHECKPOINT-BLOCK**************#\n",
    "cp.addObject(velocityField,'velocityField')\n",
    "cp.addObject(pressureField,'pressureField')\n",
    "if md.thermal:\n",
    "    cp.addObject(temperatureField,'temperatureField')\n",
    "    cp.addObject(temperatureDotField,'temperatureDotField')\n",
    "    \n",
    "\n",
    "#*************CHECKPOINT-BLOCK**************#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#print(cp.objDict.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#*************CHECKPOINT-BLOCK**************#\n",
    "\n",
    "if cp.restart:\n",
    "    velocityField.load(cp.loadpath + '/velocityField.h5')\n",
    "    pressureField.load(cp.loadpath + '/pressureField.h5')\n",
    "    if md.thermal:\n",
    "        temperatureField.load(cp.loadpath + '/temperatureField.h5')\n",
    "        temperatureDotField.load(cp.loadpath + '/temperatureDotField.h5')\n",
    "        \n",
    "#*************CHECKPOINT-BLOCK**************#\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#miscellaneous Uw functions functions\n",
    "\n",
    "coordinate = fn.input()\n",
    "depthFn = mesh.maxCoord[1] - coordinate[1] #a function providing the depth\n",
    "\n",
    "\n",
    "xFn = coordinate[0]  #a function providing the x-coordinate\n",
    "yFn = coordinate[1]\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "#Create a binary circle\n",
    "def inCircleFnGenerator(centre, radius):\n",
    "    coord = fn.input()\n",
    "    offsetFn = coord - centre\n",
    "    return fn.math.dot( offsetFn, offsetFn ) < radius**2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Boundary Conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#Stokes BCs\n",
    "\n",
    "iWalls = mesh.specialSets[\"MinI_VertexSet\"] + mesh.specialSets[\"MaxI_VertexSet\"]\n",
    "jWalls = mesh.specialSets[\"MinJ_VertexSet\"] + mesh.specialSets[\"MaxJ_VertexSet\"]\n",
    "tWalls = mesh.specialSets[\"MaxJ_VertexSet\"]\n",
    "bWalls =mesh.specialSets[\"MinJ_VertexSet\"]\n",
    "      \n",
    "        \n",
    "freeslipBC = uw.conditions.DirichletCondition( variable      = velocityField, \n",
    "                                               indexSetsPerDof = ( iWalls, jWalls) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#Energy BCs\n",
    "\n",
    "if md.thermal:\n",
    "    dirichTempBC = uw.conditions.DirichletCondition(     variable=temperatureField, \n",
    "                                              indexSetsPerDof=(tWalls,) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Swarm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#Materials\n",
    "mantleID = 0\n",
    "crustID = 1\n",
    "airID = 2      #in case we use sticky air\n",
    "\n",
    "#list of all material indexes\n",
    "material_list = [mantleID, crustID, airID]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#*************CHECKPOINT-BLOCK**************#\n",
    "\n",
    "swarm = uw.swarm.Swarm(mesh=mesh, particleEscape=True)\n",
    "cp.addObject(swarm,'swarm')\n",
    "materialVariable      = swarm.add_variable( dataType=\"int\", count=1 )\n",
    "cp.addObject(materialVariable,'materialVariable')\n",
    "proxyTempVariable = swarm.add_variable( dataType=\"double\", count=1 )\n",
    "if not md.thermal:\n",
    "    cp.addObject(proxyTempVariable,'proxyTempVariable')\n",
    "\n",
    "\n",
    "if cp.restart:\n",
    "    swarm.load(cp.loadpath + '/swarm.h5')\n",
    "    materialVariable.load(cp.loadpath + '/materialVariable.h5')\n",
    "    if not md.thermal:\n",
    "        proxyTempVariable.load(cp.loadpath + '/proxyTempVariable.h5')   \n",
    "\n",
    "\n",
    "else:\n",
    "    layout = uw.swarm.layouts.PerCellRandomLayout(swarm=swarm, particlesPerCell=int(md.ppc*md.swarmInitialFac))\n",
    "    swarm.populate_using_layout( layout=layout ) # Now use it to populate.\n",
    "    proxyTempVariable.data[:] = 0.0\n",
    "    materialVariable.data[:] = mantleID\n",
    "\n",
    "#*************CHECKPOINT-BLOCK**************#\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#These variables don;t need checkpointing. They can / should be rebuilt\n",
    "\n",
    "signedDistanceVariable = swarm.add_variable( dataType=\"double\", count=1 )\n",
    "#directorVector   = swarm.add_variable( dataType=\"double\", count=2)\n",
    "\n",
    "#directorVector.data[:,:] = 0.0\n",
    "signedDistanceVariable.data[:] = 0.0\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "\n",
    "#Pass this to Figures to see full extent\n",
    "bBox=((mesh.minCoord[0], mesh.minCoord[1]),(mesh.maxCoord[0], mesh.maxCoord[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Initial Conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#plate depth\n",
    "\n",
    "#Flags to control plate behaviour: md.opuniform, md.spuniform, md.opfixed, md.spfixed\n",
    "\n",
    "#T&S 4.126\n",
    "thicknessAtTrench = 2.32*math.sqrt(1.*ndp.slabMaxAge)\n",
    "\n",
    "sig = 150e3/sf.lengthScale\n",
    "ridgeFn = 1. -  \\\n",
    "                fn.math.exp(-1.*(xFn - ndp.leftRidge)**2/(2 * sig**2))- \\\n",
    "                fn.math.exp(-1.*(xFn - ndp.rightRidge)**2/(2 * sig**2))\n",
    "\n",
    "spAge = ndp.slabMaxAge*fn.math.abs((ndp.leftRidge - xFn)/(ndp.subZoneLoc - ndp.leftRidge))\n",
    "opAge = ndp.opMaxAge*fn.math.abs((ndp.rightRidge - xFn)/(ndp.subZoneLoc - ndp.rightRidge))\n",
    "\n",
    "if md.spuniform:\n",
    "    if not md.spfixed:\n",
    "        spAge = ridgeFn*ndp.slabMaxAge\n",
    "    else:\n",
    "        spAge = fn.misc.constant(ndp.slabMaxAge)\n",
    "        \n",
    "if md.opuniform:\n",
    "    if not md.opfixed:\n",
    "        opAge = ridgeFn*ndp.opMaxAge\n",
    "    else:\n",
    "        opAge = fn.misc.constant(ndp.opMaxAge)\n",
    "\n",
    "proxyageFn = fn.branching.conditional([(xFn <= ndp.subZoneLoc, spAge), #idea is to make this arbitrarily complex\n",
    "                                  (True, opAge)])\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Marker lines  for slab, fault, tracking"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### slab top"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function slab_top in module unsupported_dan.utilities.subduction:\n",
      "\n",
      "slab_top(trench, normal, gradientFn, ds, maxDepth, mesh)\n",
      "    Create points representing the top of a slab from trench to maxDepth\n",
      "    \n",
      "    Parameter\n",
      "    ---------\n",
      "    trench : list or list like \n",
      "            Points represnting trench location, \n",
      "    normal: list or list like\n",
      "            vector in the horizontal plane normal to trench\n",
      "    gradientFn: function\n",
      "             function that returns the dip or the slab dz/ds \n",
      "             where s is the distance along the normal vector\n",
      "    ds: float\n",
      "            distance between points, in model coordinates\n",
      "    \n",
      "    max depth: float, or list or list like\n",
      "            Maximum depth of slab\n",
      "    mesh: uw 2 mesh\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(slab_top)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#Create some slab gradient functions to use with slab_top()\n",
    "\n",
    "\n",
    "def linearGradientFn(S):\n",
    "    return np.tan(np.deg2rad(-45.))\n",
    "\n",
    "\n",
    "def circGradientFn(S):\n",
    "    if S == 0.:\n",
    "        return 0.\n",
    "    elif S < ndp.radiusOfCurv:\n",
    "        return -S/np.sqrt((ndp.radiusOfCurv**2 - S**2))\n",
    "    else:\n",
    "        return -1e5\n",
    "    \n",
    "    \n",
    "def polyGradientFn(S):\n",
    "    if S == 0.:\n",
    "        return 0.\n",
    "    else:\n",
    "        return -1*(S/ndp.radiusOfCurv)**2\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "ds = 5e3/sf.lengthScale\n",
    "normal = [1.,0.]\n",
    "\n",
    "\n",
    "\n",
    "#data = slab_top([ndp.subZoneLoc, 1.0], normal, linearGradientFn, ds, ndp.maxDepth, mesh)\n",
    "data = slab_top([ndp.subZoneLoc, 1.0], normal, polyGradientFn, ds, ndp.maxDepth, mesh)\n",
    "#data = slab_top([ndp.subZoneLoc, 1.0], normal, circGradientFn, ds, ndp.maxDepth, mesh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "slabxs = data[:,0]\n",
    "slabys = data[:,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "slabLine = markerLine2D(mesh, velocityField, slabxs, slabys, thicknessAtTrench, 1.)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n"
     ]
    }
   ],
   "source": [
    "#Assign the signed distance for the slab - in this case we only want the portion where the signed distance is positive\n",
    "sd, pts = slabLine.compute_signed_distance(swarm.particleCoordinates.data, distance=2.*thicknessAtTrench)\n",
    "signedDistanceVariable.data[np.logical_and(sd>0, sd<=slabLine.thickness)] = sd[np.logical_and(sd>0, sd<=slabLine.thickness)]\n",
    "\n",
    "\n",
    "#Note distance=2.*thicknessAtTrench: we actually want to allow distance greater than thicknessAtTrench in the kDTree query, \n",
    "#as some of these distances will not be orthogonal to the marker line, the dot product in the function will project these distances onto the normal vector\n",
    "#We'll cull distances greater than thicknessAtTrench with a numpy boolean slice - this helps things work parallel\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "slabXConds = operator.and_(xFn > slabxs.min(), xFn < slabxs.max())\n",
    "slabYConds = depthFn < 1. - slabys.min()\n",
    "\n",
    "#Two functions we'll use to limit the region of the initial thermal stancil\n",
    "slabRegion =  fn.branching.conditional([(operator.and_(slabXConds,slabYConds), True),\n",
    "                          (True, False)])\n",
    "\n",
    "#slabCirc = inCircleFnGenerator((ndp.subZoneLoc, 1.0 - ndp.maxDepth), ndp.maxDepth)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "bufferlength = 1e3/sf.lengthScale\n",
    "\n",
    "plateDepthFn = fn.branching.conditional([(depthFn < thicknessAtTrench, depthFn),\n",
    "                                        (True, 1.)])\n",
    "\n",
    "plateTempProxFn = ndp.potentialTemp*fn.math.erf((plateDepthFn)/(2.*fn.math.sqrt(1.*proxyageFn)))\n",
    "\n",
    "slabTempProx  = ndp.potentialTemp*fn.math.erf((signedDistanceVariable)/(2.*np.sqrt(1.*ndp.slabMaxAge)))\n",
    "\n",
    "\n",
    "proxytempConds = fn.branching.conditional([(signedDistanceVariable < bufferlength, plateTempProxFn),\n",
    "                          #(operator.and_(slabRegion, slabCirc), fn.misc.min(slabTempProx , plateTempProxFn)),\n",
    "                          (slabRegion,  fn.misc.min(slabTempProx , plateTempProxFn)),                 \n",
    "\n",
    "                          (True, plateTempProxFn)]) \n",
    "\n",
    "\n",
    "#*************CHECKPOINT-BLOCK**************#\n",
    "\n",
    "if not cp.restart:\n",
    "    proxyTempVariable.data[:] = proxytempConds.evaluate(swarm)\n",
    "\n",
    "#*************CHECKPOINT-BLOCK**************#\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### marker\n",
    "\n",
    "In this notebook, the markerLine helps us set up the initial distribution of weak material. We take the points that desribe the top of the slab, and we add extra particles along the top of the lower plate (so that there is weak material covering the entirety of the lower plate and slab)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "morexs = np.arange(mesh.minCoord[0] + 100e3/sf.lengthScale, ndp.subZoneLoc, ds)[:-1]\n",
    "moreys = mesh.maxCoord[1]*np.ones(morexs.shape)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#Build marker: copy the slab line, then move using the normal vector (director)\n",
    "\n",
    "markerxs = np.concatenate((morexs,slabxs[:-2]))\n",
    "markerys = np.concatenate((moreys,slabys[:-2]))\n",
    "marker = markerLine2D(mesh, velocityField,markerxs, markerys, ndp.crustThickness,  1)\n",
    "\n",
    "with marker.swarm.deform_swarm():\n",
    "    marker.swarm.particleCoordinates.data[:] += marker.director.data*ndp.crustThickness\n",
    "    \n",
    "\n",
    "marker.rebuild()\n",
    "marker.swarm.update_particle_owners()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n"
     ]
    }
   ],
   "source": [
    "#inform the mesh of the marker\n",
    "\n",
    "sd, pts0 = marker.compute_signed_distance(swarm.particleCoordinates.data, distance=thicknessAtTrench)\n",
    "sp, pts0 = marker.compute_marker_proximity(swarm.particleCoordinates.data, distance=ndp.crustThickness)\n",
    "\n",
    "#*************CHECKPOINT-BLOCK**************#\n",
    "if not cp.restart:\n",
    "    materialVariable.data[np.logical_and(sd<0,sp == marker.ID)] = sp[np.logical_and(sd<0,sp == marker.ID)]\n",
    "#*************CHECKPOINT-BLOCK**************#\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Interpolate to temperature field\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def swarmToTemp():\n",
    "\n",
    "    _ix, _weights, _dist = nn_evaluation(swarm.particleCoordinates.data, mesh.data, n=4, weighted=True)\n",
    "\n",
    "\n",
    "    #_dist.shape, mesh.data.shape\n",
    "    #if \n",
    "    tempMapTol = 0.2\n",
    "    tempMapMask = _dist.min(axis=1) < tempMapTol*(1. - mesh.minCoord[1])/mesh.elementRes[1] \n",
    "    \n",
    "    #temperatureField.data[:] = 0.\n",
    "    temperatureField.data[:] = ndp.potentialTemp #first set to dimensionless potential temp\n",
    "\n",
    "    #now used IDW to assign temp from particles to Field\n",
    "    #this is looking pretty ugly; nn_evaluation could use some grooming\n",
    "    temperatureField.data[:,0][tempMapMask] = np.average(proxyTempVariable.evaluate(swarm)[_ix][tempMapMask][:,:,0],weights=_weights[tempMapMask], axis=1)\n",
    "\n",
    "    #now cleanup any values that have fallen outside the Bcs\n",
    "\n",
    "    temperatureField.data[temperatureField.data > 1.] = ndp.potentialTemp\n",
    "    temperatureField.data[temperatureField.data < 0.] = 0.\n",
    "    \n",
    "    #and cleanup the BCs\n",
    "    \n",
    "    temperatureField.data[bWalls.data] = ndp.potentialTemp\n",
    "    temperatureField.data[tWalls.data] = 0.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#map proxy temp (swarm var) to mesh variable\n",
    "\n",
    "if not cp.restart:\n",
    "    swarmToTemp()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## choose temp field to use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "if md.thermal:\n",
    "    temperatureFn = temperatureField\n",
    "else:\n",
    "    temperatureFn = proxyTempVariable"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## adiabatic temp correction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'ndp' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m\u001b[0m",
      "\u001b[1;31mNameError\u001b[0mTraceback (most recent call last)",
      "\u001b[1;32m<ipython-input-1-0b69d7170544>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 5\u001b[1;33m \u001b[0madiabaticCorrectFn\u001b[0m \u001b[1;33m=\u001b[0m  \u001b[0mndp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpotentialTemp\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mfn\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmath\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mexp\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mndp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdissipation\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mdepthFn\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mndp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpotentialTemp\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      6\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'ndp' is not defined"
     ]
    }
   ],
   "source": [
    "#Adiabatic correction: this is added to the arrhenius laws to simulate the adiabatic component\n",
    "#note that this applies everywhere the heating that a parcel at the surface at the potential temp. would experience\n",
    "\n",
    "\n",
    "adiabaticCorrectFn =  ndp.potentialTemp*fn.math.exp(ndp.dissipation*depthFn) - ndp.potentialTemp\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#figTemp= glucifer.Figure(quality=3, boundingBox= bBox)\n",
    "#figTemp.append( glucifer.objects.Points(swarm, temperatureFn + adiabaticCorrectFn, pointSize=1))\n",
    "#figTempappend( glucifer.objects.Points(swarm, materialVariable, pointSize=1))\n",
    "#figTemp.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Rheology"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "\n",
    "symStrainrate = fn.tensor.symmetric( \n",
    "                            velocityField.fn_gradient )\n",
    "\n",
    "#Set up any functions required by the rheology\n",
    "strainRate_2ndInvariant = fn.tensor.second_invariant( \n",
    "                            fn.tensor.symmetric( \n",
    "                            velocityField.fn_gradient ))\n",
    "\n",
    "\n",
    "\n",
    "def safe_visc(func, viscmin=ndp.viscosityMin, viscmax=ndp.viscosityMax):\n",
    "    return fn.misc.max(viscmin, fn.misc.min(viscmax, func))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "##Diffusion Creep\n",
    "diffusionUM = (1./ndp.diffusionPreExp)*\\\n",
    "            fn.math.exp( ((ndp.diffusionEnergy + (depthFn*ndp.diffusionVolume))/((temperatureFn+ adiabaticCorrectFn + ndp.surfaceTemp))))\n",
    "\n",
    "diffusionLM = ndp.lowerMantleViscFac*(1./ndp.lowerMantlePreExp)*\\\n",
    "            fn.math.exp( ((ndp.lowerMantleEnergy + (depthFn*ndp.lowerMantleVolume))/((temperatureFn+ adiabaticCorrectFn + ndp.surfaceTemp))))\n",
    "\n",
    "\n",
    "    \n",
    "    \n",
    "diffusion = fn.branching.conditional( ((depthFn < ndp.lowerMantleDepth, diffusionUM ), \n",
    "                                           (True,                      diffusionLM )  ))\n",
    "\n",
    "\n",
    "\n",
    "    \n",
    "diffusion = safe_visc(diffusion, viscmax=1e5)\n",
    "\n",
    "\n",
    "    \n",
    "#mantle Plasticity\n",
    "ys =  ndp.cohesionMantle + (depthFn*ndp.frictionMantle)\n",
    "ysf = fn.misc.min(ys, ndp.yieldStressMax)\n",
    "yielding = ysf/(2.*(strainRate_2ndInvariant) + 1e-15) \n",
    "\n",
    "##Crust plasticity\n",
    "crustys =  ndp.cohesionCrust + (depthFn*ndp.frictionCrust)\n",
    "crustysf = fn.misc.min(crustys, ndp.yieldStressMax)\n",
    "crustyielding0 = crustysf/(2.*(strainRate_2ndInvariant) + 1e-15) \n",
    "\n",
    "\n",
    "#This bit phases out the weak crust (effective plastic viscosity and visc. max) over specified depths\n",
    "depthTaperFn = cosine_taper(depthFn, ndp.crustViscCutoffDepth, ndp.crustViscEndWidth)\n",
    "crustyielding = crustyielding0*(1. - depthTaperFn) + depthTaperFn*yielding\n",
    "viscmaxCrustFn = ndp.viscosityMaxCrust*(1. - depthTaperFn) + depthTaperFn*ndp.viscosityMax\n",
    "\n",
    "\n",
    "#combined rheologies for mantle and weak layer (crust / interface)\n",
    "mantleViscosityFn = safe_visc(fn.misc.min(diffusion, yielding), viscmin=ndp.viscosityMin, viscmax=ndp.viscosityMax)\n",
    "crustViscosityFn = safe_visc(crustyielding, viscmin=ndp.viscosityMinCrust, viscmax=viscmaxCrustFn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "viscosityMapFn = fn.branching.map( fn_key = materialVariable,\n",
    "                         mapping = {0:mantleViscosityFn,\n",
    "                                    1:crustViscosityFn} )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Buoyancy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#Thermal Buoyancy\n",
    "\n",
    "z_hat = ( 0.0, -1.0 )\n",
    "\n",
    "\n",
    "thermalBuoyancyFn = ndp.rayleigh*(1. - temperatureFn)\n",
    "thermalBuoyancyFn *=z_hat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Any other functions we'll need"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "###################\n",
    "#Create integral, max/min templates \n",
    "###################\n",
    "\n",
    "globRestFn = 1.\n",
    "\n",
    "def surfint(Fn = 1., rFn=globRestFn, surfaceIndexSet=mesh.specialSets[\"MaxJ_VertexSet\"]):\n",
    "    return uw.utils.Integral( Fn*rFn, mesh=mesh, integrationType='Surface', surfaceIndexSet=surfaceIndexSet)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Stokes system and solver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "got to Stokes\n"
     ]
    }
   ],
   "source": [
    "print('got to Stokes')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "stokesPIC = uw.systems.Stokes( velocityField  = velocityField, \n",
    "                                   pressureField  = pressureField,\n",
    "                                   conditions     = [freeslipBC,],\n",
    "                                   fn_viscosity   = viscosityMapFn, \n",
    "                                   fn_bodyforce   = thermalBuoyancyFn )\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[1;35m\n",
      " \n",
      "Pressure iterations:   3\n",
      "Velocity iterations:   1 (presolve)      \n",
      "Velocity iterations:  -1 (pressure solve)\n",
      "Velocity iterations:   1 (backsolve)     \n",
      "Velocity iterations:   1 (total solve)   \n",
      " \n",
      "SCR RHS  solve time: 2.1599e-01\n",
      "Pressure solve time: 2.2039e-02\n",
      "Velocity solve time: 2.1076e-01 (backsolve)\n",
      "Total solve time   : 4.9757e-01\n",
      " \n",
      "Velocity solution min/max: 0.0000e+00/0.0000e+00\n",
      "Pressure solution min/max: 0.0000e+00/0.0000e+00\n",
      " \n",
      "\u001b[00m\n"
     ]
    }
   ],
   "source": [
    "solver = uw.systems.Solver(stokesPIC)\n",
    "\n",
    "\n",
    "if md.penaltyMethod:\n",
    "    solver.set_inner_method(\"mumps\")\n",
    "    solver.options.scr.ksp_type=\"cg\"\n",
    "    solver.set_penalty(1.0e7)\n",
    "    solver.options.scr.ksp_rtol = 1.0e-4\n",
    "\n",
    "else:\n",
    "    solver.options.main.Q22_pc_type='gkgdiag'\n",
    "    solver.options.scr.ksp_rtol=5e-5\n",
    "    solver.set_inner_method('mg')\n",
    "    solver.options.mg.levels = 4\n",
    "    \n",
    "    \n",
    "    \n",
    "#avoid this initial solve if restarting\n",
    "if not cp.restart:\n",
    "    solver.solve(nonLinearIterate=True, nonLinearTolerance=md.nltol)\n",
    "    solver.print_stats()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#remove drift in pressure\n",
    "_pressure = surfint(pressureField)\n",
    "_surfLength = surfint()\n",
    "surfLength = _surfLength.evaluate()[0]\n",
    "\n",
    "pressureSurf = _pressure.evaluate()[0]   \n",
    "pressureField.data[:] -= pressureSurf/surfLength"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#figVel= glucifer.Figure(quality=3, boundingBox= bBox)\n",
    "#figVel.append( glucifer.objects.Points(swarm, fn.math.dot(velocityField,velocityField), pointSize=1))\n",
    "#figVel.append( glucifer.objects.VectorArrows(mesh, velocityField, arrowHead=5, scaling=0.0005))\n",
    "#figVel.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Setup advection-diffusion, swarm advection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "advector = uw.systems.SwarmAdvector( swarm=swarm, velocityField=velocityField, order=2 )\n",
    "\n",
    "population_control = uw.swarm.PopulationControl(swarm, deleteThreshold=0.006, splitThreshold=0.25, maxDeletions=1, maxSplits=3, aggressive=True,aggressiveThreshold=0.9, particlesPerCell=int(md.ppc))\n",
    "\n",
    "\n",
    "\n",
    "if md.thermal:\n",
    "    advDiff = uw.systems.AdvectionDiffusion( phiField       = temperatureField, \n",
    "                                         phiDotField    = temperatureDotField, \n",
    "                                         velocityField  = velocityField,\n",
    "                                         fn_sourceTerm    = 0.0,\n",
    "                                         fn_diffusivity = 1., \n",
    "                                         #conditions     = [neumannTempBC, dirichTempBC] )\n",
    "                                         conditions     = [ dirichTempBC] )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Viz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#Build a depth dependent mask for the vizualisation\n",
    "#This helps reduce the size of glucifer databases\n",
    "\n",
    "depthVariable      = swarm.add_variable( dataType=\"float\", count=1 )\n",
    "depthVariable.data[:] = depthFn.evaluate(swarm)\n",
    "\n",
    "vizVariable      = swarm.add_variable( dataType=\"int\", count=1 )\n",
    "vizVariable.data[:] = 0\n",
    "\n",
    "for index, value in enumerate(depthVariable.data[:]):\n",
    "    #print index, value\n",
    "    if np.random.rand(1)**30 > value/(mesh.maxCoord[1] - mesh.minCoord[1]):\n",
    "        vizVariable.data[index] = 1\n",
    "\n",
    "        \n",
    "        \n",
    "del index, value    #get rid of any variables that might be pointing at the .data handles (these are!)\n",
    "\n",
    "#Now randomly cull more particles \n",
    "removeRandom = True\n",
    "if removeRandom:\n",
    "    reducFac = 0.6  #0.9 > remove 90%\n",
    "    nonzs = np.where(vizVariable.data[:,0] == 1)[0].copy()\n",
    "    nstart =  nonzs.shape[0]\n",
    "    nend = int(np.ceil(nstart*reducFac))\n",
    "    np.random.shuffle(nonzs)\n",
    "    nstart, nend\n",
    "\n",
    "    vizVariable.data[nonzs[:nend]] = 0\n",
    "    del nonzs, nstart, nend"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#Set up the gLucifer stores\n",
    "\n",
    "fullpath = os.path.join(outputPath + \"gldbs/\")\n",
    "store1 = glucifer.Store(fullpath + 'subduction1.gldb')\n",
    "store2 = glucifer.Store(fullpath + 'subduction2.gldb')\n",
    "\n",
    "\n",
    "\n",
    "fig1 = glucifer.Figure(store1,figsize=(300*np.round(md.aspectRatio,2),300))\n",
    "if md.thermal:\n",
    "    fig1.append( glucifer.objects.Points(swarm, temperatureField, pointSize=2,  valueRange=[0.0, 1.0], fn_mask=vizVariable))\n",
    "else:\n",
    "    fig1.append( glucifer.objects.Points(swarm, proxyTempVariable, pointSize=2, valueRange=[0.0, 1.0],  fn_mask=vizVariable))\n",
    "\n",
    " \n",
    "\n",
    "\n",
    "fig2 = glucifer.Figure(store2,figsize=(300*np.round(md.aspectRatio,2),300))\n",
    "fig2.append( glucifer.objects.Points(swarm, viscosityMapFn, pointSize=2, fn_mask=vizVariable, logScale=True, valueRange=[100.*ndp.viscosityMin, ndp.viscosityMax]))\n",
    "fig2.append( glucifer.objects.VectorArrows(mesh, velocityField, arrowHead=1, scaling=0.00005))\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Update functions for main loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def main_update(next_image_step):\n",
    "    \n",
    "    \"\"\"\n",
    "    This includes some functionality for image / file writing at a specified frequency,\n",
    "    Assumes global variables:\n",
    "        time, step, files_freq, next_image_step\n",
    "    \n",
    "    if numerical dt exceeds next specified writing point\n",
    "    override dt make sure we hit that point\n",
    "    Set some flags so that image / file writing proceeds\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    \n",
    "    if md.thermal:\n",
    "        dt = advDiff.get_max_dt()*md.courantFac #additional md.courantFac helps stabilise advDiff\n",
    "        advDiff.integrate(dt)\n",
    "        \n",
    "    else:\n",
    "        dt = advector.get_max_dt()\n",
    "        \n",
    "    #This relates to file writing at set period:\n",
    "    #override dt make sure we hit certain time values\n",
    "    #Set some flags so that image / file writing proceeds\n",
    "    \n",
    "    if step == 0:\n",
    "        files_this_step = True\n",
    "    else:\n",
    "        files_this_step = False\n",
    "    \n",
    "    if time + dt >= next_image_step:\n",
    "        dt = next_image_step - time\n",
    "        files_this_step = True\n",
    "        next_image_step += files_freq #increment time for our next image / file dump\n",
    "        \n",
    "        \n",
    "    #Do advection    \n",
    "        \n",
    "    advector.integrate(dt)\n",
    "    marker.advection(dt)\n",
    "    \n",
    "    #remove drift in pressure\n",
    "    pressureSurf = _pressure.evaluate()[0]   \n",
    "    pressureField.data[:] -= pressureSurf/surfLength\n",
    "    \n",
    "    \n",
    "    return time+dt, step+1, files_this_step, next_image_step\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def viz_update():\n",
    "    \n",
    "    #Rebuild the viz. mask\n",
    "    vizVariable.data[:] = 0\n",
    "\n",
    "    for index, value in enumerate(depthVariable.data[:]):\n",
    "        #print index, value\n",
    "        if np.random.rand(1)**5 > value/(mesh.maxCoord[1] - mesh.minCoord[1]):\n",
    "            vizVariable.data[index] = 1\n",
    "\n",
    "    del index, value    #get rid of any variables that might be pointing at the .data handles\n",
    "    \n",
    "    if removeRandom:\n",
    "        reducFac = 0.6  #0.9 > remove 90%\n",
    "        nonzs = np.where(vizVariable.data[:,0] == 1)[0].copy()\n",
    "        nstart =  nonzs.shape[0]\n",
    "        nend = int(np.ceil(nstart*reducFac))\n",
    "        np.random.shuffle(nonzs)\n",
    "        nstart, nend\n",
    "\n",
    "        vizVariable.data[nonzs[:nend]] = 0\n",
    "        del nonzs, nstart, nend   \n",
    "    \n",
    "    #save gldbs\n",
    "    fullpath = os.path.join(outputPath + \"gldbs/\")\n",
    "    \n",
    "    store1.step = step\n",
    "    fig1.save( fullpath + \"Temp\" + str(step).zfill(5))\n",
    "    \n",
    "    store2.step = step\n",
    "    fig2.save( fullpath + \"visc\" + str(step).zfill(5))\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def swarm_update():\n",
    "    \n",
    "    #run swarm repopulation\n",
    "    population_control.repopulate()\n",
    "    \n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def markerLine_update():\n",
    "    \n",
    "    #cull particles once they reach a given depth\n",
    "    \n",
    "    cutoffDepth = 250e3/sf.lengthScale\n",
    "    mask = marker.swarm.particleCoordinates.data[:,1] < (mesh.maxCoord[1] - cutoffDepth)\n",
    "\n",
    "    with marker.swarm.deform_swarm():\n",
    "        marker.swarm.particleCoordinates.data[mask] = (999999.,999999.)\n",
    "        \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def xdmfs_update():\n",
    "    \n",
    "    #define any NN interps we'll need\n",
    "    ix1, weights1, d1 = nn_evaluation(swarm.particleCoordinates.data, mesh.data, n=5, weighted=True)\n",
    "    \n",
    "    \n",
    "    #rebuild any mesh vars that are not self-updating\n",
    "    viscosityField.data[:,0] =  np.average(viscosityMapFn.evaluate(swarm)[:,0][ix1], weights=weights1, axis=len((weights1.shape)) - 1)\n",
    "    strainRateField.data[:] = strainRate_2ndInvariant.evaluate(mesh)\n",
    "     \n",
    "    \n",
    "    fullpath = os.path.join(outputPath + \"xdmf/\")\n",
    "    #if not os.path.exists(fullpath+\"mesh.h5\"):\n",
    "    #    _mH = mesh.save(fullpath+\"mesh.h5\")\n",
    "    \n",
    "    try:\n",
    "        _mH\n",
    "    except:\n",
    "        _mH = mesh.save(fullpath+\"mesh.h5\")\n",
    "    \n",
    "    \n",
    "    #Part 1\n",
    "    mh = _mH\n",
    "    vH = velocityField.save(fullpath + \"velocity_\" + str(step) +\".h5\")\n",
    "    tH = temperatureField.save(fullpath + \"temp_\" + str(step) + \".h5\")\n",
    "    srH = strainRateField.save(fullpath + \"strainrate_\" + str(step) +\".h5\")\n",
    "    viscH = viscosityField.save(fullpath + \"visc_\" + str(step) + \".h5\")\n",
    "    \n",
    "    #part2\n",
    "    \n",
    "    velocityField.xdmf(fullpath + \"velocity_\" + str(step), vH, 'velocity', mh, 'mesh', modeltime=time)\n",
    "    temperatureField.xdmf(fullpath + \"temp_\" + str(step), tH, 'temperature', mh, 'mesh', modeltime=time)\n",
    "    strainRateField.xdmf(fullpath + \"strainrate_\" + str(step), srH, 'strainrate', mh, 'mesh', modeltime=time)\n",
    "    viscosityField.xdmf(fullpath + \"visc_\" + str(step), viscH, 'visc', mh, 'mesh', modeltime=time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Main loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "time = cp.time()  # Initial time\n",
    "step = cp.step()   # Initial timestep\n",
    "steps_output = 5   # output every 10 timesteps\n",
    "metrics_output = 5\n",
    "files_output = 10\n",
    "\n",
    "files_freq  = md.filesFreqYears*(3600.*365.*24.)/sf.time  #applies to files and gldbs\n",
    "files_this_step = False\n",
    "next_image_step = (np.floor(time/files_freq)+ 1.) *files_freq \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#checkpoint at time zero\n",
    "if not cp.restart:\n",
    "    cp.saveObjs(step, time)\n",
    "    cp.saveDicts(step, time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       ".lvctrl input { display: inline; }\n",
       ".lvctrl input[type=range] { width: 200px; }\n",
       ".lvctrl input[type=button] { width: 100px; }\n",
       ".lvctrl input[type=number] { width: 120px; height: 28px; border-radius: 0; border: 1px dotted #999;}\n",
       ".lvctrl pre { display: inline; }\n",
       ".lvctrl p  { color: gray; margin: 0px}\n",
       "\n",
       "/* Colour select */\n",
       ".colourbg {display: inline-block; margin:5px 0px; width: 100px; height: 20px; border: 1px dotted #999; }\n",
       ".colour {display: inline-block; width: 100%; height: 100%; margin: 0px; padding: 0px; border: none;}\n",
       "\n",
       ".checkerboard {background-image: url(\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAABAAAAAQCAYAAAAf8/9hAAAAIElEQVQ4jWP4TwAcOHAAL2YYNWBYGEBIASEwasCwMAAALvidroqDalkAAAAASUVORK5CYII=\");}\n",
       "\n",
       ".palette {margin: 3px; border:1px solid #000; border: 1px dotted #000;}\n",
       "\n",
       "\n",
       "</style>\n",
       "\n",
       "<script id=\"line-fs\" type=\"x-shader/x-fragment\">\n",
       "precision highp float;\n",
       "varying vec4 vColour;\n",
       "void main(void)\n",
       "{\n",
       "  gl_FragColor = vColour;\n",
       "}\n",
       "</script>\n",
       "\n",
       "<script id=\"line-vs\" type=\"x-shader/x-vertex\">\n",
       "precision highp float;\n",
       "//Line vertex shader\n",
       "attribute vec3 aVertexPosition;\n",
       "uniform mat4 uMVMatrix;\n",
       "uniform mat4 uPMatrix;\n",
       "uniform vec4 uColour;\n",
       "varying vec4 vColour;\n",
       "void main(void)\n",
       "{\n",
       "  vec4 mvPosition = uMVMatrix * vec4(aVertexPosition, 1.0);\n",
       "  gl_Position = uPMatrix * mvPosition;\n",
       "  vColour = uColour;\n",
       "}\n",
       "</script>\n",
       "<script>\n",
       "// gl-matrix 1.3.7 - https://github.com/toji/gl-matrix/blob/master/LICENSE.md\n",
       "(function(w,D){\"object\"===typeof exports?module.exports=D(global):\"function\"===false/*OK: removed require.js support.. typeof define&&define.amd*/?define([],function(){return D(w)}):D(w)})(this,function(w){function D(a){return o=a}function G(){return o=\"undefined\"!==typeof Float32Array?Float32Array:Array}var E={};(function(){if(\"undefined\"!=typeof Float32Array){var a=new Float32Array(1),b=new Int32Array(a.buffer);E.invsqrt=function(c){a[0]=c;b[0]=1597463007-(b[0]>>1);var d=a[0];return d*(1.5-0.5*c*d*d)}}else E.invsqrt=function(a){return 1/\n",
       "Math.sqrt(a)}})();var o=null;G();var r={create:function(a){var b=new o(3);a?(b[0]=a[0],b[1]=a[1],b[2]=a[2]):b[0]=b[1]=b[2]=0;return b},createFrom:function(a,b,c){var d=new o(3);d[0]=a;d[1]=b;d[2]=c;return d},set:function(a,b){b[0]=a[0];b[1]=a[1];b[2]=a[2];return b},equal:function(a,b){return a===b||1.0E-6>Math.abs(a[0]-b[0])&&1.0E-6>Math.abs(a[1]-b[1])&&1.0E-6>Math.abs(a[2]-b[2])},add:function(a,b,c){if(!c||a===c)return a[0]+=b[0],a[1]+=b[1],a[2]+=b[2],a;c[0]=a[0]+b[0];c[1]=a[1]+b[1];c[2]=a[2]+b[2];\n",
       "return c},subtract:function(a,b,c){if(!c||a===c)return a[0]-=b[0],a[1]-=b[1],a[2]-=b[2],a;c[0]=a[0]-b[0];c[1]=a[1]-b[1];c[2]=a[2]-b[2];return c},multiply:function(a,b,c){if(!c||a===c)return a[0]*=b[0],a[1]*=b[1],a[2]*=b[2],a;c[0]=a[0]*b[0];c[1]=a[1]*b[1];c[2]=a[2]*b[2];return c},negate:function(a,b){b||(b=a);b[0]=-a[0];b[1]=-a[1];b[2]=-a[2];return b},scale:function(a,b,c){if(!c||a===c)return a[0]*=b,a[1]*=b,a[2]*=b,a;c[0]=a[0]*b;c[1]=a[1]*b;c[2]=a[2]*b;return c},normalize:function(a,b){b||(b=a);var c=\n",
       "a[0],d=a[1],e=a[2],g=Math.sqrt(c*c+d*d+e*e);if(!g)return b[0]=0,b[1]=0,b[2]=0,b;if(1===g)return b[0]=c,b[1]=d,b[2]=e,b;g=1/g;b[0]=c*g;b[1]=d*g;b[2]=e*g;return b},cross:function(a,b,c){c||(c=a);var d=a[0],e=a[1],a=a[2],g=b[0],f=b[1],b=b[2];c[0]=e*b-a*f;c[1]=a*g-d*b;c[2]=d*f-e*g;return c},length:function(a){var b=a[0],c=a[1],a=a[2];return Math.sqrt(b*b+c*c+a*a)},squaredLength:function(a){var b=a[0],c=a[1],a=a[2];return b*b+c*c+a*a},dot:function(a,b){return a[0]*b[0]+a[1]*b[1]+a[2]*b[2]},direction:function(a,\n",
       "b,c){c||(c=a);var d=a[0]-b[0],e=a[1]-b[1],a=a[2]-b[2],b=Math.sqrt(d*d+e*e+a*a);if(!b)return c[0]=0,c[1]=0,c[2]=0,c;b=1/b;c[0]=d*b;c[1]=e*b;c[2]=a*b;return c},lerp:function(a,b,c,d){d||(d=a);d[0]=a[0]+c*(b[0]-a[0]);d[1]=a[1]+c*(b[1]-a[1]);d[2]=a[2]+c*(b[2]-a[2]);return d},dist:function(a,b){var c=b[0]-a[0],d=b[1]-a[1],e=b[2]-a[2];return Math.sqrt(c*c+d*d+e*e)}},H=null,y=new o(4);r.unproject=function(a,b,c,d,e){e||(e=a);H||(H=x.create());var g=H;y[0]=2*(a[0]-d[0])/d[2]-1;y[1]=2*(a[1]-d[1])/d[3]-1;y[2]=\n",
       "2*a[2]-1;y[3]=1;x.multiply(c,b,g);if(!x.inverse(g))return null;x.multiplyVec4(g,y);if(0===y[3])return null;e[0]=y[0]/y[3];e[1]=y[1]/y[3];e[2]=y[2]/y[3];return e};var L=r.createFrom(1,0,0),M=r.createFrom(0,1,0),N=r.createFrom(0,0,1),z=r.create();r.rotationTo=function(a,b,c){c||(c=k.create());var d=r.dot(a,b);if(1<=d)k.set(O,c);else if(-0.999999>d)r.cross(L,a,z),1.0E-6>r.length(z)&&r.cross(M,a,z),1.0E-6>r.length(z)&&r.cross(N,a,z),r.normalize(z),k.fromAngleAxis(Math.PI,z,c);else{var d=Math.sqrt(2*(1+\n",
       "d)),e=1/d;r.cross(a,b,z);c[0]=z[0]*e;c[1]=z[1]*e;c[2]=z[2]*e;c[3]=0.5*d;k.normalize(c)}1<c[3]?c[3]=1:-1>c[3]&&(c[3]=-1);return c};r.str=function(a){return\"[\"+a[0]+\", \"+a[1]+\", \"+a[2]+\"]\"};var A={create:function(a){var b=new o(9);a?(b[0]=a[0],b[1]=a[1],b[2]=a[2],b[3]=a[3],b[4]=a[4],b[5]=a[5],b[6]=a[6],b[7]=a[7],b[8]=a[8]):b[0]=b[1]=b[2]=b[3]=b[4]=b[5]=b[6]=b[7]=b[8]=0;return b},createFrom:function(a,b,c,d,e,g,f,h,j){var i=new o(9);i[0]=a;i[1]=b;i[2]=c;i[3]=d;i[4]=e;i[5]=g;i[6]=f;i[7]=h;i[8]=j;return i},\n",
       "determinant:function(a){var b=a[3],c=a[4],d=a[5],e=a[6],g=a[7],f=a[8];return a[0]*(f*c-d*g)+a[1]*(-f*b+d*e)+a[2]*(g*b-c*e)},inverse:function(a,b){var c=a[0],d=a[1],e=a[2],g=a[3],f=a[4],h=a[5],j=a[6],i=a[7],m=a[8],l=m*f-h*i,C=-m*g+h*j,q=i*g-f*j,n=c*l+d*C+e*q;if(!n)return null;n=1/n;b||(b=A.create());b[0]=l*n;b[1]=(-m*d+e*i)*n;b[2]=(h*d-e*f)*n;b[3]=C*n;b[4]=(m*c-e*j)*n;b[5]=(-h*c+e*g)*n;b[6]=q*n;b[7]=(-i*c+d*j)*n;b[8]=(f*c-d*g)*n;return b},multiply:function(a,b,c){c||(c=a);var d=a[0],e=a[1],g=a[2],\n",
       "f=a[3],h=a[4],j=a[5],i=a[6],m=a[7],a=a[8],l=b[0],C=b[1],q=b[2],n=b[3],k=b[4],p=b[5],o=b[6],s=b[7],b=b[8];c[0]=l*d+C*f+q*i;c[1]=l*e+C*h+q*m;c[2]=l*g+C*j+q*a;c[3]=n*d+k*f+p*i;c[4]=n*e+k*h+p*m;c[5]=n*g+k*j+p*a;c[6]=o*d+s*f+b*i;c[7]=o*e+s*h+b*m;c[8]=o*g+s*j+b*a;return c},multiplyVec2:function(a,b,c){c||(c=b);var d=b[0],b=b[1];c[0]=d*a[0]+b*a[3]+a[6];c[1]=d*a[1]+b*a[4]+a[7];return c},multiplyVec3:function(a,b,c){c||(c=b);var d=b[0],e=b[1],b=b[2];c[0]=d*a[0]+e*a[3]+b*a[6];c[1]=d*a[1]+e*a[4]+b*a[7];c[2]=\n",
       "d*a[2]+e*a[5]+b*a[8];return c},set:function(a,b){b[0]=a[0];b[1]=a[1];b[2]=a[2];b[3]=a[3];b[4]=a[4];b[5]=a[5];b[6]=a[6];b[7]=a[7];b[8]=a[8];return b},equal:function(a,b){return a===b||1.0E-6>Math.abs(a[0]-b[0])&&1.0E-6>Math.abs(a[1]-b[1])&&1.0E-6>Math.abs(a[2]-b[2])&&1.0E-6>Math.abs(a[3]-b[3])&&1.0E-6>Math.abs(a[4]-b[4])&&1.0E-6>Math.abs(a[5]-b[5])&&1.0E-6>Math.abs(a[6]-b[6])&&1.0E-6>Math.abs(a[7]-b[7])&&1.0E-6>Math.abs(a[8]-b[8])},identity:function(a){a||(a=A.create());a[0]=1;a[1]=0;a[2]=0;a[3]=0;\n",
       "a[4]=1;a[5]=0;a[6]=0;a[7]=0;a[8]=1;return a},transpose:function(a,b){if(!b||a===b){var c=a[1],d=a[2],e=a[5];a[1]=a[3];a[2]=a[6];a[3]=c;a[5]=a[7];a[6]=d;a[7]=e;return a}b[0]=a[0];b[1]=a[3];b[2]=a[6];b[3]=a[1];b[4]=a[4];b[5]=a[7];b[6]=a[2];b[7]=a[5];b[8]=a[8];return b},toMat4:function(a,b){b||(b=x.create());b[15]=1;b[14]=0;b[13]=0;b[12]=0;b[11]=0;b[10]=a[8];b[9]=a[7];b[8]=a[6];b[7]=0;b[6]=a[5];b[5]=a[4];b[4]=a[3];b[3]=0;b[2]=a[2];b[1]=a[1];b[0]=a[0];return b},str:function(a){return\"[\"+a[0]+\", \"+a[1]+\n",
       "\", \"+a[2]+\", \"+a[3]+\", \"+a[4]+\", \"+a[5]+\", \"+a[6]+\", \"+a[7]+\", \"+a[8]+\"]\"}},x={create:function(a){var b=new o(16);a&&(b[0]=a[0],b[1]=a[1],b[2]=a[2],b[3]=a[3],b[4]=a[4],b[5]=a[5],b[6]=a[6],b[7]=a[7],b[8]=a[8],b[9]=a[9],b[10]=a[10],b[11]=a[11],b[12]=a[12],b[13]=a[13],b[14]=a[14],b[15]=a[15]);return b},createFrom:function(a,b,c,d,e,g,f,h,j,i,m,l,C,q,n,k){var p=new o(16);p[0]=a;p[1]=b;p[2]=c;p[3]=d;p[4]=e;p[5]=g;p[6]=f;p[7]=h;p[8]=j;p[9]=i;p[10]=m;p[11]=l;p[12]=C;p[13]=q;p[14]=n;p[15]=k;return p},set:function(a,\n",
       "b){b[0]=a[0];b[1]=a[1];b[2]=a[2];b[3]=a[3];b[4]=a[4];b[5]=a[5];b[6]=a[6];b[7]=a[7];b[8]=a[8];b[9]=a[9];b[10]=a[10];b[11]=a[11];b[12]=a[12];b[13]=a[13];b[14]=a[14];b[15]=a[15];return b},equal:function(a,b){return a===b||1.0E-6>Math.abs(a[0]-b[0])&&1.0E-6>Math.abs(a[1]-b[1])&&1.0E-6>Math.abs(a[2]-b[2])&&1.0E-6>Math.abs(a[3]-b[3])&&1.0E-6>Math.abs(a[4]-b[4])&&1.0E-6>Math.abs(a[5]-b[5])&&1.0E-6>Math.abs(a[6]-b[6])&&1.0E-6>Math.abs(a[7]-b[7])&&1.0E-6>Math.abs(a[8]-b[8])&&1.0E-6>Math.abs(a[9]-b[9])&&1.0E-6>\n",
       "Math.abs(a[10]-b[10])&&1.0E-6>Math.abs(a[11]-b[11])&&1.0E-6>Math.abs(a[12]-b[12])&&1.0E-6>Math.abs(a[13]-b[13])&&1.0E-6>Math.abs(a[14]-b[14])&&1.0E-6>Math.abs(a[15]-b[15])},identity:function(a){a||(a=x.create());a[0]=1;a[1]=0;a[2]=0;a[3]=0;a[4]=0;a[5]=1;a[6]=0;a[7]=0;a[8]=0;a[9]=0;a[10]=1;a[11]=0;a[12]=0;a[13]=0;a[14]=0;a[15]=1;return a},transpose:function(a,b){if(!b||a===b){var c=a[1],d=a[2],e=a[3],g=a[6],f=a[7],h=a[11];a[1]=a[4];a[2]=a[8];a[3]=a[12];a[4]=c;a[6]=a[9];a[7]=a[13];a[8]=d;a[9]=g;a[11]=\n",
       "a[14];a[12]=e;a[13]=f;a[14]=h;return a}b[0]=a[0];b[1]=a[4];b[2]=a[8];b[3]=a[12];b[4]=a[1];b[5]=a[5];b[6]=a[9];b[7]=a[13];b[8]=a[2];b[9]=a[6];b[10]=a[10];b[11]=a[14];b[12]=a[3];b[13]=a[7];b[14]=a[11];b[15]=a[15];return b},determinant:function(a){var b=a[0],c=a[1],d=a[2],e=a[3],g=a[4],f=a[5],h=a[6],j=a[7],i=a[8],m=a[9],l=a[10],C=a[11],q=a[12],n=a[13],k=a[14],a=a[15];return q*m*h*e-i*n*h*e-q*f*l*e+g*n*l*e+i*f*k*e-g*m*k*e-q*m*d*j+i*n*d*j+q*c*l*j-b*n*l*j-i*c*k*j+b*m*k*j+q*f*d*C-g*n*d*C-q*c*h*C+b*n*h*C+\n",
       "g*c*k*C-b*f*k*C-i*f*d*a+g*m*d*a+i*c*h*a-b*m*h*a-g*c*l*a+b*f*l*a},inverse:function(a,b){b||(b=a);var c=a[0],d=a[1],e=a[2],g=a[3],f=a[4],h=a[5],j=a[6],i=a[7],m=a[8],l=a[9],k=a[10],q=a[11],n=a[12],o=a[13],p=a[14],r=a[15],s=c*h-d*f,v=c*j-e*f,t=c*i-g*f,u=d*j-e*h,w=d*i-g*h,x=e*i-g*j,y=m*o-l*n,z=m*p-k*n,F=m*r-q*n,A=l*p-k*o,D=l*r-q*o,E=k*r-q*p,B=s*E-v*D+t*A+u*F-w*z+x*y;if(!B)return null;B=1/B;b[0]=(h*E-j*D+i*A)*B;b[1]=(-d*E+e*D-g*A)*B;b[2]=(o*x-p*w+r*u)*B;b[3]=(-l*x+k*w-q*u)*B;b[4]=(-f*E+j*F-i*z)*B;b[5]=\n",
       "(c*E-e*F+g*z)*B;b[6]=(-n*x+p*t-r*v)*B;b[7]=(m*x-k*t+q*v)*B;b[8]=(f*D-h*F+i*y)*B;b[9]=(-c*D+d*F-g*y)*B;b[10]=(n*w-o*t+r*s)*B;b[11]=(-m*w+l*t-q*s)*B;b[12]=(-f*A+h*z-j*y)*B;b[13]=(c*A-d*z+e*y)*B;b[14]=(-n*u+o*v-p*s)*B;b[15]=(m*u-l*v+k*s)*B;return b},toRotationMat:function(a,b){b||(b=x.create());b[0]=a[0];b[1]=a[1];b[2]=a[2];b[3]=a[3];b[4]=a[4];b[5]=a[5];b[6]=a[6];b[7]=a[7];b[8]=a[8];b[9]=a[9];b[10]=a[10];b[11]=a[11];b[12]=0;b[13]=0;b[14]=0;b[15]=1;return b},toMat3:function(a,b){b||(b=A.create());b[0]=\n",
       "a[0];b[1]=a[1];b[2]=a[2];b[3]=a[4];b[4]=a[5];b[5]=a[6];b[6]=a[8];b[7]=a[9];b[8]=a[10];return b},toInverseMat3:function(a,b){var c=a[0],d=a[1],e=a[2],g=a[4],f=a[5],h=a[6],j=a[8],i=a[9],m=a[10],l=m*f-h*i,k=-m*g+h*j,q=i*g-f*j,n=c*l+d*k+e*q;if(!n)return null;n=1/n;b||(b=A.create());b[0]=l*n;b[1]=(-m*d+e*i)*n;b[2]=(h*d-e*f)*n;b[3]=k*n;b[4]=(m*c-e*j)*n;b[5]=(-h*c+e*g)*n;b[6]=q*n;b[7]=(-i*c+d*j)*n;b[8]=(f*c-d*g)*n;return b},multiply:function(a,b,c){c||(c=a);var d=a[0],e=a[1],g=a[2],f=a[3],h=a[4],j=a[5],\n",
       "i=a[6],m=a[7],l=a[8],k=a[9],q=a[10],n=a[11],o=a[12],p=a[13],r=a[14],a=a[15],s=b[0],v=b[1],t=b[2],u=b[3];c[0]=s*d+v*h+t*l+u*o;c[1]=s*e+v*j+t*k+u*p;c[2]=s*g+v*i+t*q+u*r;c[3]=s*f+v*m+t*n+u*a;s=b[4];v=b[5];t=b[6];u=b[7];c[4]=s*d+v*h+t*l+u*o;c[5]=s*e+v*j+t*k+u*p;c[6]=s*g+v*i+t*q+u*r;c[7]=s*f+v*m+t*n+u*a;s=b[8];v=b[9];t=b[10];u=b[11];c[8]=s*d+v*h+t*l+u*o;c[9]=s*e+v*j+t*k+u*p;c[10]=s*g+v*i+t*q+u*r;c[11]=s*f+v*m+t*n+u*a;s=b[12];v=b[13];t=b[14];u=b[15];c[12]=s*d+v*h+t*l+u*o;c[13]=s*e+v*j+t*k+u*p;c[14]=s*g+\n",
       "v*i+t*q+u*r;c[15]=s*f+v*m+t*n+u*a;return c},multiplyVec3:function(a,b,c){c||(c=b);var d=b[0],e=b[1],b=b[2];c[0]=a[0]*d+a[4]*e+a[8]*b+a[12];c[1]=a[1]*d+a[5]*e+a[9]*b+a[13];c[2]=a[2]*d+a[6]*e+a[10]*b+a[14];return c},multiplyVec4:function(a,b,c){c||(c=b);var d=b[0],e=b[1],g=b[2],b=b[3];c[0]=a[0]*d+a[4]*e+a[8]*g+a[12]*b;c[1]=a[1]*d+a[5]*e+a[9]*g+a[13]*b;c[2]=a[2]*d+a[6]*e+a[10]*g+a[14]*b;c[3]=a[3]*d+a[7]*e+a[11]*g+a[15]*b;return c},translate:function(a,b,c){var d=b[0],e=b[1],b=b[2],g,f,h,j,i,m,l,k,q,\n",
       "n,o,p;if(!c||a===c)return a[12]=a[0]*d+a[4]*e+a[8]*b+a[12],a[13]=a[1]*d+a[5]*e+a[9]*b+a[13],a[14]=a[2]*d+a[6]*e+a[10]*b+a[14],a[15]=a[3]*d+a[7]*e+a[11]*b+a[15],a;g=a[0];f=a[1];h=a[2];j=a[3];i=a[4];m=a[5];l=a[6];k=a[7];q=a[8];n=a[9];o=a[10];p=a[11];c[0]=g;c[1]=f;c[2]=h;c[3]=j;c[4]=i;c[5]=m;c[6]=l;c[7]=k;c[8]=q;c[9]=n;c[10]=o;c[11]=p;c[12]=g*d+i*e+q*b+a[12];c[13]=f*d+m*e+n*b+a[13];c[14]=h*d+l*e+o*b+a[14];c[15]=j*d+k*e+p*b+a[15];return c},scale:function(a,b,c){var d=b[0],e=b[1],b=b[2];if(!c||a===c)return a[0]*=\n",
       "d,a[1]*=d,a[2]*=d,a[3]*=d,a[4]*=e,a[5]*=e,a[6]*=e,a[7]*=e,a[8]*=b,a[9]*=b,a[10]*=b,a[11]*=b,a;c[0]=a[0]*d;c[1]=a[1]*d;c[2]=a[2]*d;c[3]=a[3]*d;c[4]=a[4]*e;c[5]=a[5]*e;c[6]=a[6]*e;c[7]=a[7]*e;c[8]=a[8]*b;c[9]=a[9]*b;c[10]=a[10]*b;c[11]=a[11]*b;c[12]=a[12];c[13]=a[13];c[14]=a[14];c[15]=a[15];return c},rotate:function(a,b,c,d){var e=c[0],g=c[1],c=c[2],f=Math.sqrt(e*e+g*g+c*c),h,j,i,m,l,k,q,n,o,p,r,s,v,t,u,w,x,y,z,A;if(!f)return null;1!==f&&(f=1/f,e*=f,g*=f,c*=f);h=Math.sin(b);j=Math.cos(b);i=1-j;b=a[0];\n",
       "f=a[1];m=a[2];l=a[3];k=a[4];q=a[5];n=a[6];o=a[7];p=a[8];r=a[9];s=a[10];v=a[11];t=e*e*i+j;u=g*e*i+c*h;w=c*e*i-g*h;x=e*g*i-c*h;y=g*g*i+j;z=c*g*i+e*h;A=e*c*i+g*h;e=g*c*i-e*h;g=c*c*i+j;d?a!==d&&(d[12]=a[12],d[13]=a[13],d[14]=a[14],d[15]=a[15]):d=a;d[0]=b*t+k*u+p*w;d[1]=f*t+q*u+r*w;d[2]=m*t+n*u+s*w;d[3]=l*t+o*u+v*w;d[4]=b*x+k*y+p*z;d[5]=f*x+q*y+r*z;d[6]=m*x+n*y+s*z;d[7]=l*x+o*y+v*z;d[8]=b*A+k*e+p*g;d[9]=f*A+q*e+r*g;d[10]=m*A+n*e+s*g;d[11]=l*A+o*e+v*g;return d},rotateX:function(a,b,c){var d=Math.sin(b),\n",
       "b=Math.cos(b),e=a[4],g=a[5],f=a[6],h=a[7],j=a[8],i=a[9],m=a[10],l=a[11];c?a!==c&&(c[0]=a[0],c[1]=a[1],c[2]=a[2],c[3]=a[3],c[12]=a[12],c[13]=a[13],c[14]=a[14],c[15]=a[15]):c=a;c[4]=e*b+j*d;c[5]=g*b+i*d;c[6]=f*b+m*d;c[7]=h*b+l*d;c[8]=e*-d+j*b;c[9]=g*-d+i*b;c[10]=f*-d+m*b;c[11]=h*-d+l*b;return c},rotateY:function(a,b,c){var d=Math.sin(b),b=Math.cos(b),e=a[0],g=a[1],f=a[2],h=a[3],j=a[8],i=a[9],m=a[10],l=a[11];c?a!==c&&(c[4]=a[4],c[5]=a[5],c[6]=a[6],c[7]=a[7],c[12]=a[12],c[13]=a[13],c[14]=a[14],c[15]=\n",
       "a[15]):c=a;c[0]=e*b+j*-d;c[1]=g*b+i*-d;c[2]=f*b+m*-d;c[3]=h*b+l*-d;c[8]=e*d+j*b;c[9]=g*d+i*b;c[10]=f*d+m*b;c[11]=h*d+l*b;return c},rotateZ:function(a,b,c){var d=Math.sin(b),b=Math.cos(b),e=a[0],g=a[1],f=a[2],h=a[3],j=a[4],i=a[5],m=a[6],l=a[7];c?a!==c&&(c[8]=a[8],c[9]=a[9],c[10]=a[10],c[11]=a[11],c[12]=a[12],c[13]=a[13],c[14]=a[14],c[15]=a[15]):c=a;c[0]=e*b+j*d;c[1]=g*b+i*d;c[2]=f*b+m*d;c[3]=h*b+l*d;c[4]=e*-d+j*b;c[5]=g*-d+i*b;c[6]=f*-d+m*b;c[7]=h*-d+l*b;return c},frustum:function(a,b,c,d,e,g,f){f||\n",
       "(f=x.create());var h=b-a,j=d-c,i=g-e;f[0]=2*e/h;f[1]=0;f[2]=0;f[3]=0;f[4]=0;f[5]=2*e/j;f[6]=0;f[7]=0;f[8]=(b+a)/h;f[9]=(d+c)/j;f[10]=-(g+e)/i;f[11]=-1;f[12]=0;f[13]=0;f[14]=-(2*g*e)/i;f[15]=0;return f},perspective:function(a,b,c,d,e){a=c*Math.tan(a*Math.PI/360);b*=a;return x.frustum(-b,b,-a,a,c,d,e)},ortho:function(a,b,c,d,e,g,f){f||(f=x.create());var h=b-a,j=d-c,i=g-e;f[0]=2/h;f[1]=0;f[2]=0;f[3]=0;f[4]=0;f[5]=2/j;f[6]=0;f[7]=0;f[8]=0;f[9]=0;f[10]=-2/i;f[11]=0;f[12]=-(a+b)/h;f[13]=-(d+c)/j;f[14]=\n",
       "-(g+e)/i;f[15]=1;return f},lookAt:function(a,b,c,d){d||(d=x.create());var e,g,f,h,j,i,m,l,k=a[0],o=a[1],a=a[2];f=c[0];h=c[1];g=c[2];m=b[0];c=b[1];e=b[2];if(k===m&&o===c&&a===e)return x.identity(d);b=k-m;c=o-c;m=a-e;l=1/Math.sqrt(b*b+c*c+m*m);b*=l;c*=l;m*=l;e=h*m-g*c;g=g*b-f*m;f=f*c-h*b;(l=Math.sqrt(e*e+g*g+f*f))?(l=1/l,e*=l,g*=l,f*=l):f=g=e=0;h=c*f-m*g;j=m*e-b*f;i=b*g-c*e;(l=Math.sqrt(h*h+j*j+i*i))?(l=1/l,h*=l,j*=l,i*=l):i=j=h=0;d[0]=e;d[1]=h;d[2]=b;d[3]=0;d[4]=g;d[5]=j;d[6]=c;d[7]=0;d[8]=f;d[9]=\n",
       "i;d[10]=m;d[11]=0;d[12]=-(e*k+g*o+f*a);d[13]=-(h*k+j*o+i*a);d[14]=-(b*k+c*o+m*a);d[15]=1;return d},fromRotationTranslation:function(a,b,c){c||(c=x.create());var d=a[0],e=a[1],g=a[2],f=a[3],h=d+d,j=e+e,i=g+g,a=d*h,m=d*j,d=d*i,k=e*j,e=e*i,g=g*i,h=f*h,j=f*j,f=f*i;c[0]=1-(k+g);c[1]=m+f;c[2]=d-j;c[3]=0;c[4]=m-f;c[5]=1-(a+g);c[6]=e+h;c[7]=0;c[8]=d+j;c[9]=e-h;c[10]=1-(a+k);c[11]=0;c[12]=b[0];c[13]=b[1];c[14]=b[2];c[15]=1;return c},str:function(a){return\"[\"+a[0]+\", \"+a[1]+\", \"+a[2]+\", \"+a[3]+\", \"+a[4]+\", \"+\n",
       "a[5]+\", \"+a[6]+\", \"+a[7]+\", \"+a[8]+\", \"+a[9]+\", \"+a[10]+\", \"+a[11]+\", \"+a[12]+\", \"+a[13]+\", \"+a[14]+\", \"+a[15]+\"]\"}},k={create:function(a){var b=new o(4);a?(b[0]=a[0],b[1]=a[1],b[2]=a[2],b[3]=a[3]):b[0]=b[1]=b[2]=b[3]=0;return b},createFrom:function(a,b,c,d){var e=new o(4);e[0]=a;e[1]=b;e[2]=c;e[3]=d;return e},set:function(a,b){b[0]=a[0];b[1]=a[1];b[2]=a[2];b[3]=a[3];return b},equal:function(a,b){return a===b||1.0E-6>Math.abs(a[0]-b[0])&&1.0E-6>Math.abs(a[1]-b[1])&&1.0E-6>Math.abs(a[2]-b[2])&&1.0E-6>\n",
       "Math.abs(a[3]-b[3])},identity:function(a){a||(a=k.create());a[0]=0;a[1]=0;a[2]=0;a[3]=1;return a}},O=k.identity();k.calculateW=function(a,b){var c=a[0],d=a[1],e=a[2];if(!b||a===b)return a[3]=-Math.sqrt(Math.abs(1-c*c-d*d-e*e)),a;b[0]=c;b[1]=d;b[2]=e;b[3]=-Math.sqrt(Math.abs(1-c*c-d*d-e*e));return b};k.dot=function(a,b){return a[0]*b[0]+a[1]*b[1]+a[2]*b[2]+a[3]*b[3]};k.inverse=function(a,b){var c=a[0],d=a[1],e=a[2],g=a[3],c=(c=c*c+d*d+e*e+g*g)?1/c:0;if(!b||a===b)return a[0]*=-c,a[1]*=-c,a[2]*=-c,a[3]*=\n",
       "c,a;b[0]=-a[0]*c;b[1]=-a[1]*c;b[2]=-a[2]*c;b[3]=a[3]*c;return b};k.conjugate=function(a,b){if(!b||a===b)return a[0]*=-1,a[1]*=-1,a[2]*=-1,a;b[0]=-a[0];b[1]=-a[1];b[2]=-a[2];b[3]=a[3];return b};k.length=function(a){var b=a[0],c=a[1],d=a[2],a=a[3];return Math.sqrt(b*b+c*c+d*d+a*a)};k.normalize=function(a,b){b||(b=a);var c=a[0],d=a[1],e=a[2],g=a[3],f=Math.sqrt(c*c+d*d+e*e+g*g);if(0===f)return b[0]=0,b[1]=0,b[2]=0,b[3]=0,b;f=1/f;b[0]=c*f;b[1]=d*f;b[2]=e*f;b[3]=g*f;return b};k.add=function(a,b,c){if(!c||\n",
       "a===c)return a[0]+=b[0],a[1]+=b[1],a[2]+=b[2],a[3]+=b[3],a;c[0]=a[0]+b[0];c[1]=a[1]+b[1];c[2]=a[2]+b[2];c[3]=a[3]+b[3];return c};k.multiply=function(a,b,c){c||(c=a);var d=a[0],e=a[1],g=a[2],a=a[3],f=b[0],h=b[1],j=b[2],b=b[3];c[0]=d*b+a*f+e*j-g*h;c[1]=e*b+a*h+g*f-d*j;c[2]=g*b+a*j+d*h-e*f;c[3]=a*b-d*f-e*h-g*j;return c};k.multiplyVec3=function(a,b,c){c||(c=b);var d=b[0],e=b[1],g=b[2],b=a[0],f=a[1],h=a[2],a=a[3],j=a*d+f*g-h*e,i=a*e+h*d-b*g,k=a*g+b*e-f*d,d=-b*d-f*e-h*g;c[0]=j*a+d*-b+i*-h-k*-f;c[1]=i*a+\n",
       "d*-f+k*-b-j*-h;c[2]=k*a+d*-h+j*-f-i*-b;return c};k.scale=function(a,b,c){if(!c||a===c)return a[0]*=b,a[1]*=b,a[2]*=b,a[3]*=b,a;c[0]=a[0]*b;c[1]=a[1]*b;c[2]=a[2]*b;c[3]=a[3]*b;return c};k.toMat3=function(a,b){b||(b=A.create());var c=a[0],d=a[1],e=a[2],g=a[3],f=c+c,h=d+d,j=e+e,i=c*f,k=c*h,c=c*j,l=d*h,d=d*j,e=e*j,f=g*f,h=g*h,g=g*j;b[0]=1-(l+e);b[1]=k+g;b[2]=c-h;b[3]=k-g;b[4]=1-(i+e);b[5]=d+f;b[6]=c+h;b[7]=d-f;b[8]=1-(i+l);return b};k.toMat4=function(a,b){b||(b=x.create());var c=a[0],d=a[1],e=a[2],g=\n",
       "a[3],f=c+c,h=d+d,j=e+e,i=c*f,k=c*h,c=c*j,l=d*h,d=d*j,e=e*j,f=g*f,h=g*h,g=g*j;b[0]=1-(l+e);b[1]=k+g;b[2]=c-h;b[3]=0;b[4]=k-g;b[5]=1-(i+e);b[6]=d+f;b[7]=0;b[8]=c+h;b[9]=d-f;b[10]=1-(i+l);b[11]=0;b[12]=0;b[13]=0;b[14]=0;b[15]=1;return b};k.slerp=function(a,b,c,d){d||(d=a);var e=a[0]*b[0]+a[1]*b[1]+a[2]*b[2]+a[3]*b[3],g,f;if(1<=Math.abs(e))return d!==a&&(d[0]=a[0],d[1]=a[1],d[2]=a[2],d[3]=a[3]),d;g=Math.acos(e);f=Math.sqrt(1-e*e);if(0.001>Math.abs(f))return d[0]=0.5*a[0]+0.5*b[0],d[1]=0.5*a[1]+0.5*b[1],\n",
       "d[2]=0.5*a[2]+0.5*b[2],d[3]=0.5*a[3]+0.5*b[3],d;e=Math.sin((1-c)*g)/f;c=Math.sin(c*g)/f;d[0]=a[0]*e+b[0]*c;d[1]=a[1]*e+b[1]*c;d[2]=a[2]*e+b[2]*c;d[3]=a[3]*e+b[3]*c;return d};k.fromRotationMatrix=function(a,b){b||(b=k.create());var c=a[0]+a[4]+a[8],d;if(0<c)d=Math.sqrt(c+1),b[3]=0.5*d,d=0.5/d,b[0]=(a[7]-a[5])*d,b[1]=(a[2]-a[6])*d,b[2]=(a[3]-a[1])*d;else{d=k.fromRotationMatrix.s_iNext=k.fromRotationMatrix.s_iNext||[1,2,0];c=0;a[4]>a[0]&&(c=1);a[8]>a[3*c+c]&&(c=2);var e=d[c],g=d[e];d=Math.sqrt(a[3*c+\n",
       "c]-a[3*e+e]-a[3*g+g]+1);b[c]=0.5*d;d=0.5/d;b[3]=(a[3*g+e]-a[3*e+g])*d;b[e]=(a[3*e+c]+a[3*c+e])*d;b[g]=(a[3*g+c]+a[3*c+g])*d}return b};A.toQuat4=k.fromRotationMatrix;(function(){var a=A.create();k.fromAxes=function(b,c,d,e){a[0]=c[0];a[3]=c[1];a[6]=c[2];a[1]=d[0];a[4]=d[1];a[7]=d[2];a[2]=b[0];a[5]=b[1];a[8]=b[2];return k.fromRotationMatrix(a,e)}})();k.identity=function(a){a||(a=k.create());a[0]=0;a[1]=0;a[2]=0;a[3]=1;return a};k.fromAngleAxis=function(a,b,c){c||(c=k.create());var a=0.5*a,d=Math.sin(a);\n",
       "c[3]=Math.cos(a);c[0]=d*b[0];c[1]=d*b[1];c[2]=d*b[2];return c};k.toAngleAxis=function(a,b){b||(b=a);var c=a[0]*a[0]+a[1]*a[1]+a[2]*a[2];0<c?(b[3]=2*Math.acos(a[3]),c=E.invsqrt(c),b[0]=a[0]*c,b[1]=a[1]*c,b[2]=a[2]*c):(b[3]=0,b[0]=1,b[1]=0,b[2]=0);return b};k.str=function(a){return\"[\"+a[0]+\", \"+a[1]+\", \"+a[2]+\", \"+a[3]+\"]\"};var J={create:function(a){var b=new o(2);a?(b[0]=a[0],b[1]=a[1]):(b[0]=0,b[1]=0);return b},createFrom:function(a,b){var c=new o(2);c[0]=a;c[1]=b;return c},add:function(a,b,c){c||\n",
       "(c=b);c[0]=a[0]+b[0];c[1]=a[1]+b[1];return c},subtract:function(a,b,c){c||(c=b);c[0]=a[0]-b[0];c[1]=a[1]-b[1];return c},multiply:function(a,b,c){c||(c=b);c[0]=a[0]*b[0];c[1]=a[1]*b[1];return c},divide:function(a,b,c){c||(c=b);c[0]=a[0]/b[0];c[1]=a[1]/b[1];return c},scale:function(a,b,c){c||(c=a);c[0]=a[0]*b;c[1]=a[1]*b;return c},dist:function(a,b){var c=b[0]-a[0],d=b[1]-a[1];return Math.sqrt(c*c+d*d)},set:function(a,b){b[0]=a[0];b[1]=a[1];return b},equal:function(a,b){return a===b||1.0E-6>Math.abs(a[0]-\n",
       "b[0])&&1.0E-6>Math.abs(a[1]-b[1])},negate:function(a,b){b||(b=a);b[0]=-a[0];b[1]=-a[1];return b},normalize:function(a,b){b||(b=a);var c=a[0]*a[0]+a[1]*a[1];0<c?(c=Math.sqrt(c),b[0]=a[0]/c,b[1]=a[1]/c):b[0]=b[1]=0;return b},cross:function(a,b,c){a=a[0]*b[1]-a[1]*b[0];if(!c)return a;c[0]=c[1]=0;c[2]=a;return c},length:function(a){var b=a[0],a=a[1];return Math.sqrt(b*b+a*a)},squaredLength:function(a){var b=a[0],a=a[1];return b*b+a*a},dot:function(a,b){return a[0]*b[0]+a[1]*b[1]},direction:function(a,\n",
       "b,c){c||(c=a);var d=a[0]-b[0],a=a[1]-b[1],b=d*d+a*a;if(!b)return c[0]=0,c[1]=0,c[2]=0,c;b=1/Math.sqrt(b);c[0]=d*b;c[1]=a*b;return c},lerp:function(a,b,c,d){d||(d=a);d[0]=a[0]+c*(b[0]-a[0]);d[1]=a[1]+c*(b[1]-a[1]);return d},str:function(a){return\"[\"+a[0]+\", \"+a[1]+\"]\"}},I={create:function(a){var b=new o(4);a?(b[0]=a[0],b[1]=a[1],b[2]=a[2],b[3]=a[3]):b[0]=b[1]=b[2]=b[3]=0;return b},createFrom:function(a,b,c,d){var e=new o(4);e[0]=a;e[1]=b;e[2]=c;e[3]=d;return e},set:function(a,b){b[0]=a[0];b[1]=a[1];\n",
       "b[2]=a[2];b[3]=a[3];return b},equal:function(a,b){return a===b||1.0E-6>Math.abs(a[0]-b[0])&&1.0E-6>Math.abs(a[1]-b[1])&&1.0E-6>Math.abs(a[2]-b[2])&&1.0E-6>Math.abs(a[3]-b[3])},identity:function(a){a||(a=I.create());a[0]=1;a[1]=0;a[2]=0;a[3]=1;return a},transpose:function(a,b){if(!b||a===b){var c=a[1];a[1]=a[2];a[2]=c;return a}b[0]=a[0];b[1]=a[2];b[2]=a[1];b[3]=a[3];return b},determinant:function(a){return a[0]*a[3]-a[2]*a[1]},inverse:function(a,b){b||(b=a);var c=a[0],d=a[1],e=a[2],g=a[3],f=c*g-e*\n",
       "d;if(!f)return null;f=1/f;b[0]=g*f;b[1]=-d*f;b[2]=-e*f;b[3]=c*f;return b},multiply:function(a,b,c){c||(c=a);var d=a[0],e=a[1],g=a[2],a=a[3];c[0]=d*b[0]+e*b[2];c[1]=d*b[1]+e*b[3];c[2]=g*b[0]+a*b[2];c[3]=g*b[1]+a*b[3];return c},rotate:function(a,b,c){c||(c=a);var d=a[0],e=a[1],g=a[2],a=a[3],f=Math.sin(b),b=Math.cos(b);c[0]=d*b+e*f;c[1]=d*-f+e*b;c[2]=g*b+a*f;c[3]=g*-f+a*b;return c},multiplyVec2:function(a,b,c){c||(c=b);var d=b[0],b=b[1];c[0]=d*a[0]+b*a[1];c[1]=d*a[2]+b*a[3];return c},scale:function(a,\n",
       "b,c){c||(c=a);var d=a[1],e=a[2],g=a[3],f=b[0],b=b[1];c[0]=a[0]*f;c[1]=d*b;c[2]=e*f;c[3]=g*b;return c},str:function(a){return\"[\"+a[0]+\", \"+a[1]+\", \"+a[2]+\", \"+a[3]+\"]\"}},K={create:function(a){var b=new o(4);a?(b[0]=a[0],b[1]=a[1],b[2]=a[2],b[3]=a[3]):(b[0]=0,b[1]=0,b[2]=0,b[3]=0);return b},createFrom:function(a,b,c,d){var e=new o(4);e[0]=a;e[1]=b;e[2]=c;e[3]=d;return e},add:function(a,b,c){c||(c=b);c[0]=a[0]+b[0];c[1]=a[1]+b[1];c[2]=a[2]+b[2];c[3]=a[3]+b[3];return c},subtract:function(a,b,c){c||(c=\n",
       "b);c[0]=a[0]-b[0];c[1]=a[1]-b[1];c[2]=a[2]-b[2];c[3]=a[3]-b[3];return c},multiply:function(a,b,c){c||(c=b);c[0]=a[0]*b[0];c[1]=a[1]*b[1];c[2]=a[2]*b[2];c[3]=a[3]*b[3];return c},divide:function(a,b,c){c||(c=b);c[0]=a[0]/b[0];c[1]=a[1]/b[1];c[2]=a[2]/b[2];c[3]=a[3]/b[3];return c},scale:function(a,b,c){c||(c=a);c[0]=a[0]*b;c[1]=a[1]*b;c[2]=a[2]*b;c[3]=a[3]*b;return c},set:function(a,b){b[0]=a[0];b[1]=a[1];b[2]=a[2];b[3]=a[3];return b},equal:function(a,b){return a===b||1.0E-6>Math.abs(a[0]-b[0])&&1.0E-6>\n",
       "Math.abs(a[1]-b[1])&&1.0E-6>Math.abs(a[2]-b[2])&&1.0E-6>Math.abs(a[3]-b[3])},negate:function(a,b){b||(b=a);b[0]=-a[0];b[1]=-a[1];b[2]=-a[2];b[3]=-a[3];return b},length:function(a){var b=a[0],c=a[1],d=a[2],a=a[3];return Math.sqrt(b*b+c*c+d*d+a*a)},squaredLength:function(a){var b=a[0],c=a[1],d=a[2],a=a[3];return b*b+c*c+d*d+a*a},lerp:function(a,b,c,d){d||(d=a);d[0]=a[0]+c*(b[0]-a[0]);d[1]=a[1]+c*(b[1]-a[1]);d[2]=a[2]+c*(b[2]-a[2]);d[3]=a[3]+c*(b[3]-a[3]);return d},str:function(a){return\"[\"+a[0]+\", \"+\n",
       "a[1]+\", \"+a[2]+\", \"+a[3]+\"]\"}};w&&(w.glMatrixArrayType=o,w.MatrixArray=o,w.setMatrixArrayType=D,w.determineMatrixArrayType=G,w.glMath=E,w.vec2=J,w.vec3=r,w.vec4=K,w.mat2=I,w.mat3=A,w.mat4=x,w.quat4=k);return{glMatrixArrayType:o,MatrixArray:o,setMatrixArrayType:D,determineMatrixArrayType:G,glMath:E,vec2:J,vec3:r,vec4:K,mat2:I,mat3:A,mat4:x,quat4:k}});/*\n",
       " Javascript graphics utility library\n",
       " Helper functions, WebGL classes, Mouse input, Colours and Gradients UI\n",
       " Copyright (c) 2014, Owen Kaluza\n",
       " Released into public domain:\n",
       " This program is free software. It comes without any warranty, to\n",
       " the extent permitted by applicable law. You can redistribute it\n",
       " and/or modify it as long as this header remains intact\n",
       "*/\n",
       "var OK=function(){var a={debug_on:!1,debug:function(b){if(a.debug_on){var c=document.getElementById(\"console\");c?c.innerHTML=\"<div style=\\\"font-family: 'monospace'; font-size: 8pt;\\\">\"+b+\"</div>\"+c.innerHTML:console.log(b)}},clear:function(){var a=document.getElementById(\"console\");a&&(a.innerHTML=\"\")}};return a}();function getSearchVariable(a,b){for(var c=window.location.search.substring(1).split(\"&\"),d=0;d<c.length;d++){var e=c[d].split(\"=\");if(unescape(e[0])==a)return unescape(e[1])}return b}\n",
       "function getImageDataURL(a){var b=document.createElement(\"canvas\");b.width=a.width;b.height=a.height;b.getContext(\"2d\").drawImage(a,0,0);return b.toDataURL(\"image/png\")}window.ELEMENT||(window.ELEMENT=function(a,b){return(\"object\"==typeof b?b:document).getElementById(a)});window.STYLE||(window.STYLE=function(a){if(a=document.getElementById(a))return a.style});\n",
       "window.toggle||(window.toggle=function(a){var b=document.getElementById(a).style.display;\"none\"!=b&&b?document.getElementById(a).style.display=\"none\":document.getElementById(a).style.display=\"block\"});function setAll(a,b){for(var c=document.getElementsByClassName(b),d=0;d<c.length;d++)c[d].style.display=a}function getSourceFromElement(a){var b=document.getElementById(a);if(!b)return null;a=\"\";for(b=b.firstChild;b;)3==b.nodeType&&(a+=b.textContent),b=b.nextSibling;return a}\n",
       "function removeChildren(a){if(a.hasChildNodes())for(;0<a.childNodes.length;)a.removeChild(a.firstChild)}window.requestAnimationFrame||(window.requestAnimationFrame=function(){return window.webkitRequestAnimationFrame||window.mozRequestAnimationFrame||window.oRequestAnimationFrame||window.msRequestAnimationFrame}());\n",
       "function requestFullScreen(a){a=document.getElementById(a);a.requestFullscreen?a.requestFullscreen():a.mozRequestFullScreen?a.mozRequestFullScreen():a.webkitRequestFullScreen&&a.webkitRequestFullScreen()}function typeOf(a){var b=typeof a;\"object\"===b&&(a?\"number\"!==typeof a.length||a.propertyIsEnumerable(\"length\")||\"function\"!==typeof a.splice||(b=\"array\"):b=\"null\");return b}\n",
       "function isEmpty(a){var b,c;if(\"object\"===typeOf(a))for(b in a)if(c=a[b],void 0!==c&&\"function\"!==typeOf(c))return!1;return!0}\n",
       "function ajaxReadFile(a,b,c,d,e){var f=new XMLHttpRequest,g=0;void 0!=d&&(\"number\"==typeof d?g=d:f.onprogress=d);f.onreadystatechange=function(){if(0<g&&2<f.readyState){var c=parseInt(f.responseText.length);d&&setProgress(c/g*100)}4==f.readyState&&(200==f.status?(d&&setProgress(100),OK.debug(\"RECEIVED: \"+a),b&&b(f.responseText,a)):b?b(\"Error: \"+f.status+\" : \"+a):OK.debug(\"Ajax Read File Error: returned status code \"+f.status+\" \"+f.statusText))};c?f.open(\"GET\",a+\"?d=\"+(new Date).getTime(),!0):f.open(\"GET\",\n",
       "a,!0);for(var h in e)f.setRequestHeader(h,e[h]);f.send(null)}function readURL(a,b,c){var d=new XMLHttpRequest,e=0;void 0!=c&&(\"number\"==typeof c?e=c:d.onprogress=c);d.onreadystatechange=function(){if(0<e&&2<d.readyState){var a=parseInt(d.responseText.length);c&&setProgress(a/e*100)}};b?d.open(\"GET\",a+\"?d=\"+(new Date).getTime(),!1):d.open(\"GET\",a,!1);d.overrideMimeType(\"text/plain; charset=x-user-defined\");d.send(null);if(200!=d.status)return\"\";c&&setProgress(100);return d.responseText}\n",
       "function updateProgress(a){a.lengthComputable&&(setProgress(a.loaded/a.total*100),OK.debug(a.loaded+\" / \"+a.total))}function setProgress(a){a=Math.round(a);document.getElementById(\"progressbar\").style.width=3*a+\"px\";document.getElementById(\"progressstatus\").innerHTML=a+\"%\"}\n",
       "function ajaxPost(a,b,c,d,e){var f=new XMLHttpRequest;void 0!=d&&(f.upload.onprogress=d);f.onreadystatechange=function(){4==f.readyState&&(200==f.status?(d&&setProgress(100),OK.debug(\"POST: \"+a),c&&c(f.responseText)):c?c(\"Error, status:\"+f.status):OK.debug(\"Ajax Post Error: returned status code \"+f.status+\" \"+f.statusText))};f.open(\"POST\",a,!0);\"string\"==typeof b&&f.setRequestHeader(\"Content-type\",\"application/x-www-form-urlencoded\");if(e)for(key in e)f.setRequestHeader(key,e[key]);f.send(b)}\n",
       "var defaultMouse,dragMouse;function MouseEventHandler(a,b,c,d,e,f,g){this.click=a;this.wheel=b;this.move=c;this.down=d;this.up=e;this.leave=f;this.pinch=g}\n",
       "function Mouse(a,b,c){this.element=a;this.handler=b;this.isdown=this.disabled=!1;this.button=null;this.dragged=!1;this.lastY=this.lastX=this.absoluteY=this.absoluteX=this.x=this.x=0;this.slider=null;this.spin=0;this.moveUpdate=!1;this.enableContext=c?!0:!1;a.addEventListener(\"onwheel\"in document?\"wheel\":\"mousewheel\",handleMouseWheel,!1);a.onmousedown=handleMouseDown;a.onmouseout=handleMouseLeave;document.onmouseup=handleMouseUp;document.onmousemove=handleMouseMove;a.addEventListener(\"touchstart\",\n",
       "touchHandler,!0);a.addEventListener(\"touchmove\",touchHandler,!0);a.addEventListener(\"touchend\",touchHandler,!0);a.oncontextmenu=function(){return this.mouse.enableContext}}Mouse.prototype.setDefault=function(){defaultMouse=document.mouse=this};Mouse.prototype.update=function(a){a||(a=window.event);var b=mousePageCoord(a);this.x=b[0];this.y=b[1];this.absoluteX=this.x;this.absoluteY=this.y;b=findElementPos(this.element);this.x-=b[0];this.y-=b[1];this.clientx=a.clientX-b[0];this.clienty=a.clientY-b[1]};\n",
       "function mousePageCoord(a){var b;a.pageX||a.pageY?(b=a.pageX,a=a.pageY):(b=a.clientX+document.body.scrollLeft+document.documentElement.scrollLeft,a=a.clientY+document.body.scrollTop+document.documentElement.scrollTop);return[b,a]}function elementRelativeCoord(a,b){var c=findElementPos(a);b[0]-=c[0];b[1]-=c[1]}function findElementPos(a){var b=curtop=0;do b+=a.offsetLeft,curtop+=a.offsetTop;while(a=a.offsetParent);return[b,curtop]}\n",
       "function getMouse(a){a||(a=window.event);var b=a.target.mouse;if(b)return b;for(a=a.target;a!=document;)if(a=a.parentNode,a.mouse)return a.mouse;return null}function handleMouseDown(a){var b=getMouse(a);if(!b||b.disabled)return!0;b.target=(a||window.event).target;b.dragged=!1;b.update(a);b.isdown||(b.lastX=b.absoluteX,b.lastY=b.absoluteY);b.isdown=!0;dragMouse=b;b.button=a.button;document.mouse=b;var c=!0;b.handler.down&&(c=b.handler.down(a,b));!c&&a.preventDefault&&a.preventDefault();return c}\n",
       "function handleMouseUp(a){var b=document.mouse;if(!b||b.disabled)return!0;var c=!0;b.isdown&&(b.update(a),b.handler.click&&(c=b.handler.click(a,b)),b.isdown=!1,dragMouse=null,b.button=null,b.dragged=!1);b.handler.up&&(c=c&&b.handler.up(a,b));document.mouse=defaultMouse;!c&&a.preventDefault&&a.preventDefault();return c}\n",
       "function handleMouseMove(a){var b=dragMouse?dragMouse:getMouse(a);if(!b||b.disabled)return!0;b.update(a);b.deltaX=b.absoluteX-b.lastX;b.deltaY=b.absoluteY-b.lastY;var c=!0;!b.dragged&&b.isdown&&3<Math.abs(b.deltaX)+Math.abs(b.deltaY)&&(b.dragged=!0);b.handler.move&&(c=b.handler.move(a,b));b.moveUpdate&&(b.lastX=b.absoluteX,b.lastY=b.absoluteY);!c&&a.preventDefault&&a.preventDefault();return c}\n",
       "function handleMouseWheel(a){var b=getMouse(a);if(!b||b.disabled)return!0;b.update(a);var c=!1;a.spin=0<(a.deltaY?-a.deltaY:a.wheelDelta)?1:-1;b.handler.wheel&&(c=b.handler.wheel(a,b));!c&&a.preventDefault&&a.preventDefault();return c}function handleMouseLeave(a){var b=getMouse(a);if(!b||b.disabled)return!0;var c=!0;b.handler.leave&&(c=b.handler.leave(a,b));!c&&a.preventDefault&&a.preventDefault();return a.returnValue=c}\n",
       "function touchHandler(a){var b=a.changedTouches[0],c=null,d=getMouse(a);switch(a.type){case \"touchstart\":2==a.touches.length?(d.isdown=!1,d.scaling=0):c=\"mousedown\";break;case \"touchmove\":if(null!=d.scaling&&2==a.touches.length){var e=Math.sqrt((a.touches[0].pageX-a.touches[1].pageX)*(a.touches[0].pageX-a.touches[1].pageX)+(a.touches[0].pageY-a.touches[1].pageY)*(a.touches[0].pageY-a.touches[1].pageY));0<d.scaling?(a.distance=e-d.scaling,d.handler.pinch&&d.handler.pinch(a,d),a.returnValue=!0):d.scaling=\n",
       "e}else c=\"mousemove\";break;case \"touchend\":null!=d.scaling?d.scaling=0==d.scaling?null:0:c=\"mouseup\";break;default:return}1<a.touches.length&&(c=null);c&&(d=document.createEvent(\"MouseEvent\"),d.initMouseEvent(c,!0,!0,window,1,b.screenX,b.screenY,b.clientX,b.clientY,a.ctrlKey,a.altKey,a.shiftKey,a.metaKey,0,null),b.target.dispatchEvent(d),a.preventDefault())}function Viewport(a,b,c,d){this.x=a;this.y=b;this.width=c;this.height=d}\n",
       "function WebGL(a,b){this.program=null;this.modelView=new ViewMatrix;this.perspective=new ViewMatrix;this.textures=[];this.timer=null;if(!window.WebGLRenderingContext)throw\"No browser WebGL support\";try{this.gl=a.getContext(\"webgl\",b)||a.getContext(\"experimental-webgl\",b)}catch(c){throw OK.debug(\"detectGL exception: \"+c),\"No context\";}this.viewport=new Viewport(0,0,a.width,a.height);if(!this.gl)throw\"Failed to get context\";}\n",
       "WebGL.prototype.setMatrices=function(){this.gl.uniformMatrix4fv(this.program.mvMatrixUniform,!1,this.modelView.matrix);this.gl.uniformMatrix4fv(this.program.pMatrixUniform,!1,this.perspective.matrix);if(this.program.nMatrixUniform){var a=mat4.create(this.modelView.matrix);mat4.inverse(a);mat4.transpose(a);this.gl.uniformMatrix4fv(this.program.nMatrixUniform,!1,a)}};\n",
       "WebGL.prototype.initDraw2d=function(){this.gl.viewport(this.viewport.x,this.viewport.y,this.viewport.width,this.viewport.height);this.gl.enableVertexAttribArray(this.program.attributes.aVertexPosition);this.gl.bindBuffer(this.gl.ARRAY_BUFFER,this.vertexPositionBuffer);this.gl.vertexAttribPointer(this.program.attributes.aVertexPosition,this.vertexPositionBuffer.itemSize,this.gl.FLOAT,!1,0,0);this.program.attributes.aTextureCoord&&(this.gl.enableVertexAttribArray(this.program.attributes.aTextureCoord),\n",
       "this.gl.bindBuffer(this.gl.ARRAY_BUFFER,this.textureCoordBuffer),this.gl.vertexAttribPointer(this.program.attributes.aTextureCoord,this.textureCoordBuffer.itemSize,this.gl.FLOAT,!1,0,0));this.setMatrices()};WebGL.prototype.updateTexture=function(a,b,c){void 0==c&&(c=this.gl.TEXTURE0);this.gl.activeTexture(c);this.gl.bindTexture(this.gl.TEXTURE_2D,a);this.gl.texImage2D(this.gl.TEXTURE_2D,0,this.gl.RGBA,this.gl.RGBA,this.gl.UNSIGNED_BYTE,b);this.gl.bindTexture(this.gl.TEXTURE_2D,null)};\n",
       "WebGL.prototype.init2dBuffers=function(a){void 0==a&&(a=this.gl.TEXTURE0);this.vertexPositionBuffer=this.gl.createBuffer();this.gl.bindBuffer(this.gl.ARRAY_BUFFER,this.vertexPositionBuffer);this.gl.bufferData(this.gl.ARRAY_BUFFER,new Float32Array([1,1,-1,1,1,-1,-1,-1]),this.gl.STATIC_DRAW);this.vertexPositionBuffer.itemSize=2;this.vertexPositionBuffer.numItems=4;this.gl.activeTexture(a);this.gradientTexture=this.gl.createTexture();this.gl.bindTexture(this.gl.TEXTURE_2D,this.gradientTexture);this.gl.texParameteri(this.gl.TEXTURE_2D,\n",
       "this.gl.TEXTURE_MAG_FILTER,this.gl.NEAREST);this.gl.texParameteri(this.gl.TEXTURE_2D,this.gl.TEXTURE_MIN_FILTER,this.gl.NEAREST);this.textureCoordBuffer=this.gl.createBuffer();this.gl.bindBuffer(this.gl.ARRAY_BUFFER,this.textureCoordBuffer);this.gl.bufferData(this.gl.ARRAY_BUFFER,new Float32Array([1,1,0,1,1,0,0,0]),this.gl.STATIC_DRAW);this.textureCoordBuffer.itemSize=2;this.textureCoordBuffer.numItems=4};\n",
       "WebGL.prototype.loadTexture=function(a,b,c,d){void 0==b&&(b=this.gl.NEAREST);void 0==c&&(c=this.gl.RGBA);this.texid=this.textures.length;this.textures.push(this.gl.createTexture());this.gl.bindTexture(this.gl.TEXTURE_2D,this.textures[this.texid]);d&&this.gl.pixelStorei(this.gl.UNPACK_FLIP_Y_WEBGL,!0);this.gl.texImage2D(this.gl.TEXTURE_2D,0,c,c,this.gl.UNSIGNED_BYTE,a);this.gl.texParameteri(this.gl.TEXTURE_2D,this.gl.TEXTURE_MAG_FILTER,b);this.gl.texParameteri(this.gl.TEXTURE_2D,this.gl.TEXTURE_MIN_FILTER,\n",
       "b);this.gl.texParameteri(this.gl.TEXTURE_2D,this.gl.TEXTURE_WRAP_S,this.gl.CLAMP_TO_EDGE);this.gl.texParameteri(this.gl.TEXTURE_2D,this.gl.TEXTURE_WRAP_T,this.gl.CLAMP_TO_EDGE);this.gl.bindTexture(this.gl.TEXTURE_2D,null);return this.textures[this.texid]};WebGL.prototype.setPerspective=function(a,b,c,d){this.perspective.matrix=mat4.perspective(a,b,c,d)};WebGL.prototype.use=function(a){this.program=a;this.program.program&&this.gl.useProgram(this.program.program)};\n",
       "WebGL.prototype.view=function(a){if(this.gl){this.gl.viewport(0,0,this.gl.viewportWidth,this.gl.viewportHeight);this.gl.clear(this.gl.COLOR_BUFFER_BIT|this.gl.DEPTH_BUFFER_BIT);this.setPerspective(a.fov,this.gl.viewportWidth/this.gl.viewportHeight,a.near_clip,a.far_clip);this.modelView.identity();this.modelView.translate([a.translate[0],a.translate[1],a.translate[2]]);var b=[-(a.focus[0]-a.centre[0]),-(a.focus[1]-a.centre[1]),-(a.focus[2]-a.centre[2])];this.modelView.translate(b);this.modelView.mult(quat4.toMat4(a.rotate));\n",
       "this.modelView.scale(a.scale);this.modelView.translate([-b[0],-b[1],-b[2]]);this.modelView.translate([-a.focus[0],-a.focus[1],-a.focus[2]*a.orientation]);1==a.orientation?this.gl.frontFace(a.gl.CCW):(this.gl.frontFace(a.gl.CW),this.modelView.scale([1,1,-1]))}};\n",
       "function WebGLProgram(a,b,c){this.program=null;0>b.indexOf(\"main\")&&(b=getSourceFromElement(b));0>c.indexOf(\"main\")&&(c=getSourceFromElement(c));this.gl=a;this.program&&this.gl.isProgram(this.program)&&(this.gl.isShader(this.vshader)&&(this.gl.detachShader(this.program,this.vshader),this.gl.deleteShader(this.vshader)),this.gl.isShader(this.fshader)&&(this.gl.detachShader(this.program,this.fshader),this.gl.deleteShader(this.fshader)),this.gl.deleteProgram(this.program));this.program=this.gl.createProgram();\n",
       "this.vshader=this.compileShader(b,this.gl.VERTEX_SHADER);this.fshader=this.compileShader(c,this.gl.FRAGMENT_SHADER);this.gl.attachShader(this.program,this.vshader);this.gl.attachShader(this.program,this.fshader);this.gl.linkProgram(this.program);if(!this.gl.getProgramParameter(this.program,this.gl.LINK_STATUS))throw\"Could not initialise shaders: \"+this.gl.getProgramInfoLog(this.program);}\n",
       "WebGLProgram.prototype.compileShader=function(a,b){var c=this.gl.createShader(b);this.gl.shaderSource(c,a);this.gl.compileShader(c);if(!this.gl.getShaderParameter(c,this.gl.COMPILE_STATUS))throw this.gl.getShaderInfoLog(c);return c};\n",
       "WebGLProgram.prototype.setup=function(a,b,c){if(this.program){void 0==a&&(a=[\"aVertexPosition\",\"aTextureCoord\"]);this.attributes={};for(var d in a)this.attributes[a[d]]=this.gl.getAttribLocation(this.program,a[d]),c||this.gl.enableVertexAttribArray(this.attributes[a[d]]);this.uniforms={};for(d in b)this.uniforms[b[d]]=this.gl.getUniformLocation(this.program,b[d]);this.mvMatrixUniform=this.gl.getUniformLocation(this.program,\"uMVMatrix\");this.pMatrixUniform=this.gl.getUniformLocation(this.program,\"uPMatrix\");\n",
       "this.nMatrixUniform=this.gl.getUniformLocation(this.program,\"uNMatrix\")}};function ViewMatrix(){this.matrix=mat4.create();mat4.identity(this.matrix);this.stack=[]}ViewMatrix.prototype.toString=function(){return JSON.stringify(this.toArray())};ViewMatrix.prototype.toArray=function(){return JSON.parse(mat4.str(this.matrix))};ViewMatrix.prototype.push=function(a){a?(this.stack.push(mat4.create(a)),this.matrix=mat4.create(a)):this.stack.push(mat4.create(this.matrix))};\n",
       "ViewMatrix.prototype.pop=function(){if(0==this.stack.length)throw\"Matrix stack underflow\";return this.matrix=this.stack.pop()};ViewMatrix.prototype.mult=function(a){mat4.multiply(this.matrix,a)};ViewMatrix.prototype.identity=function(){mat4.identity(this.matrix)};ViewMatrix.prototype.scale=function(a){mat4.scale(this.matrix,a)};ViewMatrix.prototype.translate=function(a){mat4.translate(this.matrix,a)};ViewMatrix.prototype.rotate=function(a,b){mat4.rotate(this.matrix,a*Math.PI/180,b)};\n",
       "function Palette(a,b){this.premultiply=b;this.background=new Colour(\"rgba(0,0,0,0)\");this.colours=[];this.slider=new Image;this.slider.src=\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAAkAAAAPCAYAAAA2yOUNAAAAj0lEQVQokWNIjHT8/+zZs//Pnj37/+TJk/9XLp/+f+bEwf9HDm79v2Prqv9aKrz/GUYVEaeoMDMQryJXayWIoi0bFmFV1NWS+z/E1/Q/AwMDA0NVcez/LRsWoSia2luOUAADVcWx/xfO6/1/5fLp/1N7y//HhlmhKoCBgoyA/w3Vyf8jgyyxK4CBUF8zDAUAAJRXY0G1eRgAAAAASUVORK5CYII=\";if(a){var c=!1;if(\"string\"==typeof a)for(var d=a.split(/[\\n;]/),e,f=0;f<\n",
       "d.length;f++){var g=d[f].trim();if(g){var h=g.split(\"=\");if(\"Background\"==h[0])this.background=new Colour(h[1]);else if(\"P\"==h[0][0])e=parseFloat(h[1]);else if(\"C\"==h[0][0]){if(this.colours.push(new ColourPos(h[1],e)),1==e)break}else 2==h.length?this.colours.push(new ColourPos(h[1],h[0])):(c=!0,this.colours.push(new ColourPos(g)))}}else{for(d=0;d<a.length;d++)void 0==a[d].position&&(c=!0),this.colours.push(new ColourPos(a[d].colour,a[d].position));a.background&&(this.background=new Colour(a.background))}if(c)for(d=\n",
       "0;d<this.colours.length;d++)this.colours[d].position=1/(this.colours.length-1)*d;this.sort();c=!1;for(d=0;d<this.colours.length;d++)0<this.colours[d].colour.alpha&&(c=!0),1<this.colours[d].colour.alpha&&(this.colours[d].colour.alpha=1);if(!c)for(d=0;d<this.colours.length;d++)this.colours[d].colour.alpha=1}else this.colours.push(new ColourPos(\"rgba(255,255,255,1)\",0)),this.colours.push(new ColourPos(\"rgba(0,0,0,1)\",1))}\n",
       "Palette.prototype.sort=function(){this.colours.sort(function(a,b){return a.position-b.position})};Palette.prototype.newColour=function(a,b){var c=new ColourPos(b,a);this.colours.push(c);this.sort();for(c=1;c<this.colours.length-1;c++)if(this.colours[c].position==a)return c;return-1};Palette.prototype.inRange=function(a,b,c){for(var d=0;d<this.colours.length;d++){var e=this.colours[d].position*c;if(a==e||1<b&&a>=e-b/2&&a<=e+b/2)return d}return-1};\n",
       "Palette.prototype.inDragRange=function(a,b,c){for(var d=1;d<this.colours.length-1;d++){var e=this.colours[d].position*c;if(a==e||1<b&&a>=e-b/2&&a<=e+b/2)return d}return 0};Palette.prototype.remove=function(a){this.colours.splice(a,1)};Palette.prototype.toString=function(){for(var a=\"Background=\"+this.background.html(),b=0;b<this.colours.length;b++)a+=\"\\n\"+this.colours[b].position.toFixed(6)+\"=\"+this.colours[b].colour.html();return a};\n",
       "Palette.prototype.get=function(){var a={};a.background=this.background.html();a.colours=[];for(var b=0;b<this.colours.length;b++)a.colours.push({position:this.colours[b].position,colour:this.colours[b].colour.html()});return a};Palette.prototype.toJSON=function(){return JSON.stringify(this.get())};\n",
       "Palette.prototype.draw=function(a,b){if(!this.slider.width&&b){var c=this;setTimeout(function(){c.draw(a,b)},150)}else if(a){var d=/webkit/.test(navigator.userAgent.toLowerCase());0==this.colours.length&&(this.background=new Colour(\"#ffffff\"),this.colours.push(new ColourPos(\"#000000\",0)),this.colours.push(new ColourPos(\"#ffffff\",1)));list=this.colours.slice(0);list.sort(function(a,b){return a.position-b.position});if(a.getContext){var e=a.width,f=a.height,g=a.getContext(\"2d\");g.clearRect(0,0,e,f);\n",
       "if(d)for(var h=0,d=1;d<list.length;d++){var k=Math.round(e*list[d].position);g.fillStyle=g.createLinearGradient(h,0,k,0);var l=list[d-1].colour,m=list[d].colour;this.premultiply&&!b&&(l=this.background.blend(l),m=this.background.blend(m));g.fillStyle.addColorStop(0,l.html());g.fillStyle.addColorStop(1,m.html());g.fillRect(h,0,k-h,f);h=k}else{g.fillStyle=g.createLinearGradient(0,0,e,0);for(d=0;d<list.length;d++)h=list[d].colour,this.premultiply&&!b&&(h=this.background.blend(h)),g.fillStyle.addColorStop(list[d].position,\n",
       "h.html());g.fillRect(0,0,e,f)}if(f=document.getElementById(\"backgroundCUR\"))f.style.background=this.background.html();if(b)for(d=1;d<list.length-1;d++)f=Math.floor(e*list[d].position)+.5,50<list[d].colour.HSV().V?g.strokeStyle=\"black\":g.strokeStyle=\"white\",g.beginPath(),g.moveTo(f,0),g.lineTo(f,a.height),g.closePath(),g.stroke(),f-=this.slider.width/2,g.drawImage(this.slider,f,0)}else alert(\"getContext failed!\")}else alert(\"Invalid canvas!\")};\n",
       "function ColourPos(a,b){this.position=void 0==b?0:parseFloat(b);if(0<=this.position&&1>=this.position)this.colour=a?\"object\"==typeof a?a:new Colour(a):new Colour(\"#000000\");else throw\"Invalid Colour Position: \"+b;}\n",
       "function Colour(a){\"undefined\"==typeof a?this.set(\"#ffffff\"):\"string\"==typeof a?this.set(a):\"object\"==typeof a?\"undefined\"!=typeof a.H?this.setHSV(a):\"undefined\"!=typeof a.red?(this.red=a.red,this.green=a.green,this.blue=a.blue,this.alpha=a.alpha):a.R?(this.red=a.R,this.green=a.G,this.blue=a.B,this.alpha=\"undefined\"==typeof a.A?1:a.A):(this.red=a[0],this.green=a[1],this.blue=a[2],1>=this.red&&1>=this.green&&1>=this.blue&&(this.red=Math.round(255*this.red),this.green=Math.round(255*this.green),this.blue=\n",
       "Math.round(255*this.blue)),this.alpha=\"undefined\"==typeof a[3]?1:a[3]):this.fromInt(a)}\n",
       "Colour.prototype.set=function(a){a||(a=\"#ffffff\");var b=/^rgba?\\((\\d{1,3})\\s*,\\s*(\\d{1,3})\\s*,\\s*(\\d{1,3})\\s*,?\\s*(\\d\\.?\\d*)?\\)$/.exec(a);b?(this.red=parseInt(b[1]),this.green=parseInt(b[2]),this.blue=parseInt(b[3]),this.alpha=\"undefined\"==typeof b[4]?1:parseFloat(b[4])):\"#\"==a.charAt(0)?(a=a.substring(1,7),this.alpha=1,this.red=parseInt(a.substring(0,2),16),this.green=parseInt(a.substring(2,4),16),this.blue=parseInt(a.substring(4,6),16)):this.fromInt(parseInt(a))};\n",
       "Colour.prototype.fromInt=function(a){this.red=a&255;this.green=(a&65280)>>>8;this.blue=(a&16711680)>>>16;this.alpha=((a&4278190080)>>>24)/255};Colour.prototype.toInt=function(){var a=this.red,a=a+(this.green<<8),a=a+(this.blue<<16);return a+=Math.round(255*this.alpha)<<24};Colour.prototype.toString=function(){return this.html()};Colour.prototype.html=function(){return\"rgba(\"+this.red+\",\"+this.green+\",\"+this.blue+\",\"+this.alpha.toFixed(2)+\")\"};\n",
       "Colour.prototype.rgbaGL=function(){return new Float32Array([this.red/255,this.green/255,this.blue/255,this.alpha])};Colour.prototype.rgbaGLSL=function(){var a=this.rgbaGL();return\"rgba(\"+a[0].toFixed(4)+\",\"+a[1].toFixed(4)+\",\"+a[2].toFixed(4)+\",\"+a[3].toFixed(4)+\")\"};Colour.prototype.rgba=function(){return[this.red/255,this.green/255,this.blue/255,this.alpha]};Colour.prototype.rgbaObj=function(){return{R:this.red,G:this.green,B:this.blue,A:this.alpha}};Colour.prototype.print=function(){OK.debug(this.printString(!0))};\n",
       "Colour.prototype.printString=function(a){return\"R:\"+this.red+\" G:\"+this.green+\" B:\"+this.blue+(a?\" A:\"+this.alpha:\"\")};Colour.prototype.HEX=function(a){a=Math.round(Math.min(Math.max(0,a),255));return\"0123456789ABCDEF\".charAt((a-a%16)/16)+\"0123456789ABCDEF\".charAt(a%16)};Colour.prototype.htmlHex=function(a){return\"#\"+this.HEX(this.red)+this.HEX(this.green)+this.HEX(this.blue)};Colour.prototype.hex=function(a){return this.HEX(this.red)+this.HEX(this.green)+this.HEX(this.blue)+this.HEX(255*this.alpha)};\n",
       "Colour.prototype.hexGL=function(a){return this.HEX(255*this.alpha)+this.HEX(this.blue)+this.HEX(this.green)+this.HEX(this.red)};\n",
       "Colour.prototype.setHSV=function(a){var b,c,d,e,f;f=a.S/100;var g=a.V/100,h=a.H/360;if(0<f){1<=h&&(h=0);h*=6;F=h-Math.floor(h);d=Math.round(255*g*(1-f));e=Math.round(255*g*(1-f*F));f=Math.round(255*g*(1-f*(1-F)));g=Math.round(255*g);switch(Math.floor(h)){case 0:b=g;c=f;e=d;break;case 1:b=e;c=g;e=d;break;case 2:b=d;c=g;e=f;break;case 3:b=d;c=e;e=g;break;case 4:b=f;c=d;e=g;break;case 5:b=g,c=d}this.red=b?b:0;this.green=c?c:0;this.blue=e?e:0}else this.blue=this.green=this.red=g=Math.round(255*g);this.alpha=\n",
       "\"undefined\"==typeof a.A?1:a.A};Colour.prototype.HSV=function(){var a=this.red/255,b=this.green/255,c=this.blue/255,d=Math.max(a,b,c);deltaMax=d-Math.min(a,b,c);var e,f,g,h,k;0==deltaMax?e=f=0:(e=deltaMax/d,g=((d-a)/6+deltaMax/2)/deltaMax,h=((d-b)/6+deltaMax/2)/deltaMax,k=((d-c)/6+deltaMax/2)/deltaMax,a==d?f=k-h:b==d?f=1/3+g-k:c==d&&(f=2/3+h-g),0>f&&(f+=1),1<f&&--f);return{H:360*f,S:100*e,V:100*d}};Colour.prototype.HSVA=function(){var a=this.HSV();a.A=this.alpha;return a};\n",
       "Colour.prototype.interpolate=function(a,b){this.red=Math.round(this.red+b*(a.red-this.red));this.green=Math.round(this.green+b*(a.green-this.green));this.blue=Math.round(this.blue+b*(a.blue-this.blue));this.alpha=Math.round(this.alpha+b*(a.alpha-this.alpha))};Colour.prototype.blend=function(a){return new Colour([Math.round((1-a.alpha)*this.red+a.alpha*a.red),Math.round((1-a.alpha)*this.green+a.alpha*a.green),Math.round((1-a.alpha)*this.blue+a.alpha*a.blue),(1-a.alpha)*this.alpha+a.alpha*a.alpha])};\n",
       "function MoveWindow(a){if(a){this.element=document.getElementById(a);if(!this.element)return alert(\"No such element: \"+a),null;this.mouse=new Mouse(this.element,this);this.mouse.moveUpdate=!0;this.element.mouse=this.mouse}}\n",
       "MoveWindow.prototype.open=function(a,b){var c=this.element.style;0>a&&(a=0);0>b&&(b=0);void 0!=a&&(c.left=a+\"px\");void 0!=b&&(c.top=b+\"px\");c.display=\"block\";var d=this.element.offsetWidth,e=this.element.offsetHeight;a+d>window.innerWidth-20&&(c.left=window.innerWidth-d-20+\"px\");b+e>window.innerHeight-20&&(c.top=window.innerHeight-e-20+\"px\")};MoveWindow.prototype.close=function(){this.element.style.display=\"none\"};\n",
       "MoveWindow.prototype.move=function(a,b){if(b.isdown&&!(0<b.button)){var c=b.element.style;c.left=parseInt(c.left)+b.deltaX+\"px\";c.top=parseInt(c.top)+b.deltaY+\"px\"}};MoveWindow.prototype.down=function(a,b){return!1};function scale(a,b,c,d){return clamp(d*a/b,c,d)}function clamp(a,b,c){return Math.max(b,Math.min(c,a))}\n",
       "function ColourPicker(a,b){function c(a,b,c){var d=document.createElement(\"div\");d.id=a;b&&(d.innerHTML=b);c&&(d.style.cssText=c);return d}var d=document.getElementById(\"picker\");if(d&&d.picker)return d.picker.savefn=a,d.picker.abortfn=b,d.picker;d=document.body;this.element=c(\"picker\",null,\"display:none; top: 58px; z-index: 20; background: #0d0d0d; color: #aaa; cursor: move; font-family: arial; font-size: 11px; padding: 7px 10px 11px 10px; position: fixed; width: 248px; border-radius: 5px; border: 1px solid #444;\");\n",
       "var e=c(\"pickCURBG\",null,'background-image: url(\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAABAAAAAQCAYAAAAf8/9hAAAAIElEQVQ4jWP4TwAcOHAAL2YYNWBYGEBIASEwasCwMAAALvidroqDalkAAAAASUVORK5CYII=\"); float: left; width: 12px; height: 12px; margin-right: 3px;');e.appendChild(c(\"pickCUR\",null,\"float: left; width: 12px; height: 12px; background: #fff; margin-right: 3px;\"));this.element.appendChild(e);e=c(\"pickRGB\",\"R: 255 G: 255 B: 255\",\"float: left; position: relative; top: -1px;\");e.onclick=\"colours.picker.updateString()\";\n",
       "this.element.appendChild(e);this.element.appendChild(c(\"pickCLOSE\",\"X\",\"float: right; cursor: pointer; margin: 0 8px 3px;\"));this.element.appendChild(c(\"pickOK\",\"OK\",\"float: right; cursor: pointer; margin: 0 8px 3px;\"));e=c(\"SV\",null,\"position: relative; cursor: crosshair; float: left; height: 170px; width: 170px; margin-right: 10px; background: url('data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAABAAAAAQCAYAAAAf8/9hAAAEG0lEQVQ4jQEQBO/7APz8/Pz7+/vx+/v75Pr6+tb6+vrF+Pj4tPf396H4+PiO9/f3e/X19Wfz8/NU8PDwQuvr6zLi4uIjzs7OFZmZmQoA8PDw/O/v7/Ht7e3l7Ozs2Ozs7Mjq6uq35ubmpeXl5ZLf39+A3NzcbtXV1VvMzMxLvr6+O6ioqCyEhIQfQEBAFADk5OT84eHh8uDg4Obe3t7Z3Nzcy9nZ2brV1dWq0NDQmcrKyofCwsJ2uLi4ZKqqqlSYmJhFfX19N1lZWSsnJychANPT0/zT09Pz0NDQ6c3NzdzKysrNx8fHv8DAwK+6urqfsrKyj6mpqX+cnJxvjIyMX3l5eVBeXl5EPz8/ORsbGy8Aw8PD/MHBwfS+vr7qurq63ra2ttKxsbHErKystaOjo6eampqXj4+PiYODg3lycnJrXl5eX0hISFIuLi5IEBAQPwCwsLD9r6+v9aysrOynp6fioqKi1p2dncmVlZW8jo6OroODg6F5eXmUa2trhl1dXXlLS0ttNzc3YiIiIlkNDQ1RAJ6env2bm5v2l5eX7pSUlOWPj4/aiIiIz4GBgcN5eXm3cHBwq2RkZJ5XV1eSSkpKhzk5OX0qKipzGBgYawgICGMAioqK/YeHh/eDg4PvgICA6Hp6et90dHTVbW1ty2VlZcBcXFy1UVFRqkZGRqA6OjqWLS0tjSEhIYQSEhJ9BgYGdwB2dnb+c3Nz+HFxcfJra2vrZmZm42JiYttaWlrRUlJSyUtLS79CQkK2Nzc3rS0tLaQiIiKdGBgYlQ4ODo8EBASKAGNjY/5gYGD5XV1d9FpaWu5VVVXnTk5O4UlJSdlCQkLRPDw8yTQ0NMEqKiq7IiIisxkZGa0RERGmCgoKoQMDA5wAUFBQ/k9PT/pKSkr3R0dH8kNDQ+w+Pj7mOTk54DMzM9otLS3TJycnzSAgIMgZGRnBExMTvA0NDbcHBweyAwMDrwA9PT3+PDw8+zo6Ovg2Njb0MzMz8DAwMOwqKirnJSUl4iEhId4cHBzYFxcX1BISEtAODg7KCQkJxwQEBMQBAQHBAC0tLf4rKyv9Kioq+iYmJvclJSX0ISEh8R4eHu4aGhrqFhYW5xMTE+MQEBDgDQ0N3AgICNkGBgbWBAQE0wAAANEAHh4e/h0dHf0bGxv7Ghoa+hgYGPcWFhb2FBQU8xEREfEPDw/uDAwM7AoKCuoICAjoBgYG5gMDA+MBAQHiAAAA4QARERH+EBAQ/g8PD/0NDQ38DQ0N+wsLC/kKCgr4CAgI9wcHB/YFBQX0BAQE8wICAvIBAQHwAQEB7wAAAO8AAADuAAUFBf4FBQX+BAQE/gQEBP4DAwP+AwMD/QMDA/0CAgL8AQEB/AEBAfsAAAD7AAAA+wAAAPoAAAD6AAAA+QAAAPmq2NbsCl2m4wAAAABJRU5ErkJggg==') no-repeat; background-size: 100%;\");\n",
       "e.appendChild(c(\"SVslide\",null,\"background: url('data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAAkAAAAJCAYAAADgkQYQAAAALUlEQVQYlWNgQAX/kTBW8B8ZYFMIk0ARQFaIoQCbQuopIspNRPsOrpABSzgBAFHzU61KjdKlAAAAAElFTkSuQmCC'); height: 9px; width: 9px; position: absolute; cursor: crosshair\"));this.element.appendChild(e);e=c(\"H\",null,'cursor: crosshair; float: left; height: 170px; position: relative; width: 19px; padding: 0;background-image: url(\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAABAAAAAQCAYAAAAf8/9hAAAAIElEQVQ4jWP4TwAcOHAAL2YYNWBYGEBIASEwasCwMAAALvidroqDalkAAAAASUVORK5CYII=\");');\n",
       "e.appendChild(c(\"Hmodel\",null,\"position: relative;\"));e.appendChild(c(\"Hslide\",null,'top: 0px; left: -5px; background: url(\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAB0AAAAFCAYAAAC5Fuf5AAAAKklEQVQokWP4////fwY6gv////9n+A8F9LIQxVJaW4xiz4D5lB4WIlsMAPjER7mTpG/OAAAAAElFTkSuQmCC\"); height: 5px; width: 29px; position: absolute; '));this.element.appendChild(e);e=c(\"O\",null,'cursor: crosshair; float: left; height: 170px; position: relative; width: 19px; padding: 0;background-image: url(\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAABAAAAAQCAYAAAAf8/9hAAAAIElEQVQ4jWP4TwAcOHAAL2YYNWBYGEBIASEwasCwMAAALvidroqDalkAAAAASUVORK5CYII=\");border: 1px solid #888; left: 9px;');\n",
       "e.appendChild(c(\"Omodel\",null,\"position: relative;\"));e.appendChild(c(\"Oslide\",null,'top: 0px; left: -5px; background: url(\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAB0AAAAFCAYAAAC5Fuf5AAAAKklEQVQokWP4////fwY6gv////9n+A8F9LIQxVJaW4xiz4D5lB4WIlsMAPjER7mTpG/OAAAAAElFTkSuQmCC\"); height: 5px; width: 29px; position: absolute; '));this.element.appendChild(e);d.appendChild(this.element);d=document.createElement(\"style\");d.styleSheet?d.styleSheet.cssText=\"#pickRGB:hover {color: #FFD000;} #pickCLOSE:hover {color: #FFD000;} #pickOK:hover {color: #FFD000;}\":\n",
       "d.appendChild(document.createTextNode(\"#pickRGB:hover {color: #FFD000;} #pickCLOSE:hover {color: #FFD000;} #pickOK:hover {color: #FFD000;}\"));document.getElementsByTagName(\"head\")[0].appendChild(d);MoveWindow.call(this,\"picker\");this.savefn=a;this.abortfn=b;this.size=170;this.sv=5;this.oh=2;this.picked={H:360,S:100,V:100,A:1};this.max={H:360,S:100,V:100,A:1};this.colour=new Colour;for(var e=\"\",f,d=0;d<=this.size;d++)f=new Colour({H:Math.round(360/this.size*d),S:100,V:100,A:1}),e+=\"<div class='hue' style='height: 1px; width: 19px; margin: 0; padding: 0; background: \"+\n",
       "f.htmlHex()+\";'> </div>\";document.getElementById(\"Hmodel\").innerHTML=e;e=\"\";for(d=0;d<=this.size;d++)f=1-d/this.size,e+=\"<div class='opacity' style='height: 1px; width: 19px; margin: 0; padding: 0; background: #000;opacity: \"+f.toFixed(2)+\";'> </div>\";document.getElementById(\"Omodel\").innerHTML=e;this.element.picker=this}ColourPicker.prototype=new MoveWindow;ColourPicker.prototype.constructor=MoveWindow;\n",
       "ColourPicker.prototype.pick=function(a,b,c){this.update(a.HSVA());\"block\"!=this.element.style.display&&MoveWindow.prototype.open.call(this,b,c)};ColourPicker.prototype.select=function(a,b,c){if(!b||!c){var d=findElementPos(a);b=b?b:d[0]+32;c=c?c:d[1]+32}d=new Colour(a.style.backgroundColor);this.update(d.HSVA());\"block\"!=this.element.style.display&&(MoveWindow.prototype.open.call(this,b,c),this.target=a)};\n",
       "ColourPicker.prototype.click=function(a,b){if(\"pickCLOSE\"==b.target.id)this.abortfn&&this.abortfn(),toggle(\"picker\");else if(\"pickOK\"==b.target.id){this.savefn&&this.savefn(this.picked);if(this.target){var c=new Colour(this.picked);this.target.style.backgroundColor=c.html()}toggle(\"picker\")}else\"SV\"==b.target.id?this.setSV(b):\"Hslide\"==b.target.id||\"hue\"==b.target.className?this.setHue(b):\"Oslide\"!=b.target.id&&\"opacity\"!=b.target.className||this.setOpacity(b)};\n",
       "ColourPicker.prototype.move=function(a,b){b.isdown&&0==b.button&&(\"picker\"==b.target.id||\"pickCUR\"==b.target.id||\"pickRGB\"==b.target.id?MoveWindow.prototype.move.call(this,a,b):b.target&&this.click(a,b))};ColourPicker.prototype.wheel=function(a,b){this.incHue(-a.spin)};\n",
       "ColourPicker.prototype.setSV=function(a){var b=a.clientx-parseInt(document.getElementById(\"SV\").offsetLeft);a=a.clienty-parseInt(document.getElementById(\"SV\").offsetTop);this.picked.S=scale(b,this.size,0,this.max.S);this.picked.V=this.max.V-scale(a,this.size,0,this.max.V);this.update(this.picked)};ColourPicker.prototype.setHue=function(a){parseInt(document.getElementById(\"H\").offsetLeft);a=a.clienty-parseInt(document.getElementById(\"H\").offsetTop);this.picked.H=scale(a,this.size,0,this.max.H);this.update(this.picked)};\n",
       "ColourPicker.prototype.incHue=function(a){this.picked.H+=a;this.picked.H=clamp(this.picked.H,0,this.max.H);this.update(this.picked)};ColourPicker.prototype.setOpacity=function(a){parseInt(document.getElementById(\"O\").offsetLeft);a=a.clienty-parseInt(document.getElementById(\"O\").offsetTop);this.picked.A=1-clamp(a/this.size,0,1);this.update(this.picked)};ColourPicker.prototype.updateString=function(a){a||(a=prompt(\"Edit colour:\",this.colour.html()));a&&(this.colour=new Colour(a),this.update(this.colour.HSV()))};\n",
       "ColourPicker.prototype.update=function(a){this.picked=a;this.colour=new Colour(a);rgba=this.colour.rgbaObj();rgbaStr=this.colour.html();bgcol=new Colour({H:a.H,S:100,V:100,A:255});document.getElementById(\"pickRGB\").innerHTML=this.colour.printString();document.getElementById(\"pickCUR\").style.background=rgbaStr;document.getElementById(\"pickCUR\").style.backgroundColour=rgbaStr;document.getElementById(\"SV\").style.backgroundColor=bgcol.htmlHex();document.getElementById(\"Hslide\").style.top=a.H/360*this.size-\n",
       "this.oh+\"px\";document.getElementById(\"SVslide\").style.top=Math.round(this.size-a.V/100*this.size-this.sv)+\"px\";document.getElementById(\"SVslide\").style.left=Math.round(a.S/100*this.size-this.sv)+\"px\";document.getElementById(\"Oslide\").style.top=this.size*(1-a.A)-this.oh-1+\"px\"};\n",
       "function GradientEditor(a,b,c,d,e){this.canvas=a;this.callback=b;this.premultiply=c;this.changed=!0;this.inserting=!1;this.element=this.editing=null;this.spin=0;this.scrollable=e;d||(this.picker=new ColourPicker(this.save.bind(this),this.cancel.bind(this)));this.palette=new Palette(null,c);this.canvas.mouse=new Mouse(this.canvas,this);this.canvas.oncontextmenu=\"return false;\";this.canvas.oncontextmenu=function(){return!1}}\n",
       "GradientEditor.prototype.read=function(a){this.palette=new Palette(a,this.premultiply);this.reset();this.update(!0)};GradientEditor.prototype.update=function(a){this.changed=!0;this.palette.draw(this.canvas,!0);!a&&this.callback&&this.callback(this)};GradientEditor.prototype.get=function(a,b){if(b&&!this.changed)return!1;this.changed=!1;this.palette.draw(a,!1);return!0};\n",
       "GradientEditor.prototype.insert=function(a,b,c){this.inserting=!0;var d=new Colour;this.editing=this.palette.newColour(a,d);this.picker=new ColourPicker(this.save.bind(this),this.cancel.bind(this));this.picker.pick(d,b,c)};GradientEditor.prototype.editBackground=function(a){this.editing=-1;var b=findElementPos(a);this.element=a;this.picker=new ColourPicker(this.save.bind(this),this.cancel.bind(this));this.picker.pick(this.palette.background,b[0]+32,b[1]+32)};\n",
       "GradientEditor.prototype.edit=function(a,b,c){this.picker=new ColourPicker(this.save.bind(this),this.cancel.bind(this));\"number\"==typeof a?(this.editing=a,this.picker.pick(this.palette.colours[a].colour,b,c)):\"object\"==typeof a&&(this.cancel(),this.element=a,b=new Colour(a.style.backgroundColor),a=findElementPos(a),this.picker.pick(b,a[0]+32,a[1]+32));this.update()};\n",
       "GradientEditor.prototype.save=function(a){null!=this.editing&&(0<=this.editing?this.palette.colours[this.editing].colour.setHSV(a):this.palette.background.setHSV(a));if(this.element){var b=new Colour(0);b.setHSV(a);this.element.style.backgroundColor=b.html();if(this.element.onchange)this.element.onchange()}this.reset();this.update()};GradientEditor.prototype.cancel=function(){0<=this.editing&&this.inserting&&this.palette.remove(this.editing);this.reset();this.update()};\n",
       "GradientEditor.prototype.reset=function(){this.inserting=!1;this.element=this.editing=null};\n",
       "GradientEditor.prototype.click=function(a,b){if(a.ctrlKey){for(var c=0;c<this.palette.colours.length;c++)this.palette.colours[c].position=1-this.palette.colours[c].position;this.update();return!1}this.scrollable||(b.x=b.clientx);if(null!=b.slider)return b.slider=null,this.palette.sort(),this.update(),!1;var d=this.canvas;if(d.getContext){this.cancel();d.getContext(\"2d\");var e=findElementPos(d)[1]+30,c=this.palette.inRange(b.x,this.palette.slider.width,d.width);0<=c?0==a.button?this.edit(c,a.clientX-\n",
       "128,e):2==a.button&&(this.palette.remove(c),this.update()):this.insert(b.x/d.width,a.clientX-128,e)}return!1};GradientEditor.prototype.down=function(a,b){return!1};\n",
       "GradientEditor.prototype.move=function(a,b){if(!b.isdown)return!0;this.scrollable||(b.x=b.clientx);if(null==b.slider){var c=this.palette.inDragRange(b.x,this.palette.slider.width,this.canvas.width);0<c&&(b.slider=c)}null==b.slider?b.isdown=!1:(1>b.x&&(b.x=1),b.x>this.canvas.width-1&&(b.x=this.canvas.width-1),this.palette.colours[b.slider].position=b.x/this.canvas.width,this.update(!0))};\n",
       "GradientEditor.prototype.wheel=function(a,b){this.timer?clearTimeout(this.timer):this.canvas.style.cursor=\"wait\";this.spin+=.01*a.spin;var c=this;this.timer=setTimeout(function(){c.cycle(c.spin);c.spin=0},150)};GradientEditor.prototype.leave=function(a,b){};\n",
       "GradientEditor.prototype.cycle=function(a){this.canvas.style.cursor=\"default\";this.timer=null;for(var b=1;b<this.palette.colours.length-1;b++){var c=this.palette.colours[b].position,c=c+a;0>=c&&(c+=1);1<=c&&--c;this.palette.colours[b].position=c}this.palette.sort();this.update()};\n",
       "//Bounding box only in WebGL\n",
       "//Urgent TODO: move shared code with draw.js into common file\n",
       "//maintaining two copies currently!\n",
       "function initBox(el, cmd_callback) {\n",
       "  //console.log(\"INITBOX: \" + el.id);\n",
       "  var canvas = document.createElement(\"canvas\");\n",
       "  canvas.id = \"canvas\";\n",
       "\n",
       "  if (!el) el = document.body.firstChild;\n",
       "  el.parentElement.appendChild(canvas);\n",
       "  //canvas.style.cssText = \"position: absolute; z-index: 0; margin: 0px; padding: 0px; border: none; display: block;\"\n",
       "  //canvas.style.cssText = \"position: absolute; width: 100%; height: 100%; margin: auto; top: 0; left: 0; bottom: 0; right: 0; z-index: 51; border: none;\"\n",
       "  canvas.style.cssText = \"position: absolute; width: 100%; height: 100%; margin: 0px; padding: 0px; top: 0; left: 0; bottom: 0; right: 0; z-index: 11; border: none;\"\n",
       "  viewer = new Viewer(canvas);\n",
       "\n",
       "  //Canvas event handling\n",
       "  canvas.mouse = new Mouse(canvas, new MouseEventHandler(canvasMouseClick, canvasMouseWheel, canvasMouseMove, canvasMouseDown, null, null, canvasMousePinch));\n",
       "  //Following two settings should probably be defaults?\n",
       "  canvas.mouse.moveUpdate = true; //Continual update of deltaX/Y\n",
       "  //canvas.mouse.setDefault();\n",
       "  canvas.mouse.wheelTimer = true; //Accumulate wheel scroll (prevents too many events backing up)\n",
       "  defaultMouse = document.mouse = canvas.mouse;\n",
       "\n",
       "  //Attach viewer object to canvas\n",
       "  canvas.viewer = viewer;\n",
       "\n",
       "  //Command callback function\n",
       "  viewer.command = cmd_callback;\n",
       "\n",
       "  return viewer;\n",
       "}\n",
       "\n",
       "function updateBox(viewer, loaderfn) {\n",
       "  //Loader callback\n",
       "  loaderfn(function(data) {viewer.loadFile(data);});\n",
       "}\n",
       "\n",
       "function canvasMouseClick(event, mouse) {\n",
       "  if (mouse.element.viewer.rotating)\n",
       "    mouse.element.viewer.command('' + mouse.element.viewer.getRotationString());\n",
       "  else\n",
       "    mouse.element.viewer.command('' + mouse.element.viewer.getTranslationString());\n",
       "\n",
       "  if (mouse.element.viewer.rotating) {\n",
       "    mouse.element.viewer.rotating = false;\n",
       "    //mouse.element.viewer.reload = true;\n",
       "  }\n",
       "\n",
       "  //Clear the webgl box\n",
       "  mouse.element.viewer.clear();\n",
       "\n",
       "  return false;\n",
       "}\n",
       "\n",
       "function canvasMouseDown(event, mouse) {\n",
       "  return false;\n",
       "}\n",
       "\n",
       "function canvasMouseMove(event, mouse) {\n",
       "  if (!mouse.isdown || !mouse.element.viewer) return true;\n",
       "  mouse.element.viewer.rotating = false;\n",
       "\n",
       "  //Switch buttons for translate/rotate\n",
       "  var button = mouse.button;\n",
       "\n",
       "  //console.log(mouse.deltaX + \",\" + mouse.deltaY);\n",
       "  switch (button)\n",
       "  {\n",
       "    case 0:\n",
       "      mouse.element.viewer.rotateY(mouse.deltaX/5);\n",
       "      mouse.element.viewer.rotateX(mouse.deltaY/5);\n",
       "      mouse.element.viewer.rotating = true;\n",
       "      break;\n",
       "    case 1:\n",
       "      mouse.element.viewer.rotateZ(Math.sqrt(mouse.deltaX*mouse.deltaX + mouse.deltaY*mouse.deltaY)/5);\n",
       "      mouse.element.viewer.rotating = true;\n",
       "      break;\n",
       "    case 2:\n",
       "      var adjust = mouse.element.viewer.modelsize / 1000;   //1/1000th of size\n",
       "      mouse.element.viewer.translate[0] += mouse.deltaX * adjust;\n",
       "      mouse.element.viewer.translate[1] -= mouse.deltaY * adjust;\n",
       "      break;\n",
       "  }\n",
       "\n",
       "  mouse.element.viewer.draw();\n",
       "\n",
       "  return false;\n",
       "}\n",
       "\n",
       "var zoomTimer;\n",
       "var zoomClipTimer;\n",
       "var zoomSpin = 0;\n",
       "\n",
       "function canvasMouseWheel(event, mouse) {\n",
       "  if (event.shiftKey) {\n",
       "    var factor = event.spin * 0.01;\n",
       "    if (window.navigator.platform.indexOf(\"Mac\") >= 0)\n",
       "      factor *= 0.1;\n",
       "    if (zoomClipTimer) clearTimeout(zoomClipTimer);\n",
       "    zoomClipTimer = setTimeout(function () {mouse.element.viewer.zoomClip(factor);}, 100 );\n",
       "  } else {\n",
       "    if (zoomTimer) \n",
       "      clearTimeout(zoomTimer);\n",
       "    zoomSpin += event.spin;\n",
       "    zoomTimer = setTimeout(function () {viewer.zoom(zoomSpin*0.01); zoomSpin = 0;}, 100 );\n",
       "  }\n",
       "  return false; //Prevent default\n",
       "}\n",
       "\n",
       "function canvasMousePinch(event, mouse) {\n",
       "  if (event.distance != 0) {\n",
       "    var factor = event.distance * 0.0001;\n",
       "    mouse.element.viewer.zoom(factor);\n",
       "    //Clear the box after a second\n",
       "    setTimeout(function() {mouse.element.viewer.clear();}, 1000);\n",
       "  }\n",
       "  return false; //Prevent default\n",
       "}\n",
       "\n",
       "//This object encapsulates a vertex buffer and shader set\n",
       "function Renderer(gl, colour) {\n",
       "  this.gl = gl;\n",
       "  if (colour) this.colour = colour;\n",
       "\n",
       "    //Line renderer\n",
       "    this.attributes = [\"aVertexPosition\"],\n",
       "    this.uniforms = [\"uColour\"]\n",
       "    this.attribSizes = [3 * Float32Array.BYTES_PER_ELEMENT];\n",
       "\n",
       "  this.elements = 0;\n",
       "  this.elementSize = 0;\n",
       "  for (var i=0; i<this.attribSizes.length; i++)\n",
       "    this.elementSize += this.attribSizes[i];\n",
       "}\n",
       "\n",
       "Renderer.prototype.init = function() {\n",
       "  //Compile the shaders\n",
       "  this.program = new WebGLProgram(this.gl, \"line-vs\", \"line-fs\");\n",
       "  if (this.program.errors) console.log(this.program.errors);\n",
       "  //Setup attribs/uniforms (flag set to skip enabling attribs)\n",
       "  this.program.setup(this.attributes, this.uniforms, true);\n",
       "\n",
       "  return true;\n",
       "}\n",
       "\n",
       "Renderer.prototype.updateBuffers = function(view) {\n",
       "  //Create buffer if not yet allocated\n",
       "  if (this.vertexBuffer == undefined) {\n",
       "    //Init shaders etc...\n",
       "    if (!this.init()) return;\n",
       "    this.vertexBuffer = this.gl.createBuffer();\n",
       "    this.indexBuffer = this.gl.createBuffer();\n",
       "  }\n",
       "\n",
       "  //Bind buffers\n",
       "  this.gl.bindBuffer(this.gl.ARRAY_BUFFER, this.vertexBuffer);\n",
       "  this.gl.bindBuffer(this.gl.ELEMENT_ARRAY_BUFFER, this.indexBuffer);\n",
       "\n",
       "  this.box(view.min, view.max);\n",
       "  this.elements = 24;\n",
       "}\n",
       "\n",
       "Renderer.prototype.box = function(min, max) {\n",
       "  var vertices = new Float32Array(\n",
       "        [\n",
       "          min[0], min[1], max[2],\n",
       "          min[0], max[1], max[2],\n",
       "          max[0], max[1], max[2],\n",
       "          max[0], min[1], max[2],\n",
       "          min[0], min[1], min[2],\n",
       "          min[0], max[1], min[2],\n",
       "          max[0], max[1], min[2],\n",
       "          max[0], min[1], min[2]\n",
       "        ]);\n",
       "\n",
       "  var indices = new Uint16Array(\n",
       "        [\n",
       "          0, 1, 1, 2, 2, 3, 3, 0,\n",
       "          4, 5, 5, 6, 6, 7, 7, 4,\n",
       "          0, 4, 3, 7, 1, 5, 2, 6\n",
       "        ]\n",
       "     );\n",
       "  this.gl.bufferData(this.gl.ARRAY_BUFFER, vertices, this.gl.STATIC_DRAW);\n",
       "  this.gl.bufferData(this.gl.ELEMENT_ARRAY_BUFFER, indices, this.gl.STATIC_DRAW);\n",
       "}\n",
       "\n",
       "Renderer.prototype.draw = function(webgl) {\n",
       "  if (!this.elements) return;\n",
       "\n",
       "  if (this.program.attributes[\"aVertexPosition\"] == undefined) return; //Require vertex buffer\n",
       "\n",
       "  webgl.use(this.program);\n",
       "  webgl.setMatrices();\n",
       "\n",
       "  //Bind buffers\n",
       "  this.gl.bindBuffer(this.gl.ARRAY_BUFFER, this.vertexBuffer);\n",
       "  this.gl.bindBuffer(this.gl.ELEMENT_ARRAY_BUFFER, this.indexBuffer);\n",
       "\n",
       "  //Enable attributes\n",
       "  for (var key in this.program.attributes)\n",
       "    this.gl.enableVertexAttribArray(this.program.attributes[key]);\n",
       "\n",
       "  if (this.colour)\n",
       "    this.gl.uniform4f(this.program.uniforms[\"uColour\"], this.colour[0], this.colour[1], this.colour[2], this.colour[3]);\n",
       " \n",
       "  //Line box render\n",
       "  this.gl.vertexAttribPointer(this.program.attributes[\"aVertexPosition\"], 3, this.gl.FLOAT, false, 0, 0);\n",
       "  this.gl.drawElements(this.gl.LINES, this.elements, this.gl.UNSIGNED_SHORT, 0);\n",
       "\n",
       "  //Disable attribs\n",
       "  for (var key in this.program.attributes)\n",
       "    this.gl.disableVertexAttribArray(this.program.attributes[key]);\n",
       "\n",
       "  this.gl.bindBuffer(this.gl.ARRAY_BUFFER, null);\n",
       "  this.gl.bindBuffer(this.gl.ELEMENT_ARRAY_BUFFER, null);\n",
       "  this.gl.useProgram(null);\n",
       "}\n",
       "\n",
       "//This object holds the viewer details and calls the renderers\n",
       "function Viewer(canvas) {\n",
       "  this.canvas = canvas;\n",
       "  if (!canvas) {alert(\"Invalid Canvas\"); return;}\n",
       "  try {\n",
       "    this.webgl = new WebGL(this.canvas, {antialias: true, premultipliedAlpha: false});\n",
       "    this.gl = this.webgl.gl;\n",
       "    canvas.addEventListener(\"webglcontextlost\", function(event) {\n",
       "      event.preventDefault();\n",
       "      console.log(\"CONTEXT LOSS DETECTED\");\n",
       "    }, false);\n",
       "\n",
       "  } catch(e) {\n",
       "    //No WebGL\n",
       "    console.log(\"No WebGL: \" + e);\n",
       "  }\n",
       "\n",
       "  this.vis = {};\n",
       "  this.view = 0; //Active view\n",
       "\n",
       "  this.rotating = false;\n",
       "  this.translate = [0,0,0];\n",
       "  this.rotate = quat4.create();\n",
       "  quat4.identity(this.rotate);\n",
       "  this.fov = 45;\n",
       "  this.focus = [0,0,0];\n",
       "  this.centre = [0,0,0];\n",
       "  this.near_clip = this.far_clip = 0.0;\n",
       "  this.modelsize = 1;\n",
       "  this.scale = [1, 1, 1];\n",
       "  this.orientation = 1.0; //1.0 for RH, -1.0 for LH\n",
       "\n",
       "  if (!this.gl) return;\n",
       "\n",
       "  //Create the renderers\n",
       "  this.border = new Renderer(this.gl, [0.5,0.5,0.5,1]);\n",
       "\n",
       "  this.gl.enable(this.gl.DEPTH_TEST);\n",
       "  this.gl.depthFunc(this.gl.LEQUAL);\n",
       "  //this.gl.depthMask(this.gl.FALSE);\n",
       "  this.gl.enable(this.gl.BLEND);\n",
       "  //this.gl.blendFunc(this.gl.SRC_ALPHA, this.gl.ONE_MINUS_SRC_ALPHA);\n",
       "  //this.gl.blendFuncSeparate(this.gl.SRC_ALPHA, this.gl.ONE_MINUS_SRC_ALPHA, this.gl.ZERO, this.gl.ONE);\n",
       "  this.gl.blendFuncSeparate(this.gl.SRC_ALPHA, this.gl.ONE_MINUS_SRC_ALPHA, this.gl.ONE, this.gl.ONE_MINUS_SRC_ALPHA);\n",
       "}\n",
       "\n",
       "Viewer.prototype.checkPointMinMax = function(coord) {\n",
       "  for (var i=0; i<3; i++) {\n",
       "    this.vis.views[this.view].min[i] = Math.min(coord[i], this.vis.views[this.view].min[i]);\n",
       "    this.vis.views[this.view].max[i] = Math.max(coord[i], this.vis.views[this.view].max[i]);\n",
       "  }\n",
       "  //console.log(JSON.stringify(this.vis.views[this.view].min) + \" -- \" + JSON.stringify(this.vis.views[this.view].max));\n",
       "}\n",
       "\n",
       "Viewer.prototype.loadFile = function(source) {\n",
       "  //Skip update to rotate/translate etc if in process of updating\n",
       "  //if (document.mouse.isdown) return;\n",
       "\n",
       "  //Replace data\n",
       "  this.vis = JSON.parse(source);\n",
       "\n",
       "  //Always set a bounding box, get from objects if not in view\n",
       "  var objbb = false;\n",
       "  if (!this.vis.views[this.view].min || !this.vis.views[this.view].max) {\n",
       "    this.vis.views[this.view].min = [Number.MAX_VALUE, Number.MAX_VALUE, Number.MAX_VALUE];\n",
       "    this.vis.views[this.view].max = [-Number.MAX_VALUE, -Number.MAX_VALUE, -Number.MAX_VALUE];\n",
       "    objbb = true;\n",
       "  }\n",
       "  //console.log(JSON.stringify(this.vis.views[this.view],null, 2));\n",
       "\n",
       "  if (this.vis.views[this.view]) {\n",
       "    this.near_clip = this.vis.views[this.view].near || 0;\n",
       "    this.far_clip = this.vis.views[this.view].far || 0;\n",
       "    this.orientation = this.vis.views[this.view].orientation || 1;\n",
       "    this.axes = this.vis.views[this.view].axis == undefined ? true : this.vis.views[this.view].axis;\n",
       "  }\n",
       "\n",
       "  if (this.vis.properties.resolution && this.vis.properties.resolution[0] && this.vis.properties.resolution[1]) {\n",
       "    this.width = this.vis.properties.resolution[0];\n",
       "    this.height = this.vis.properties.resolution[1];\n",
       "    //this.canvas.style.width = \"\";\n",
       "    //this.canvas.style.height = \"\";\n",
       "  } else {\n",
       "    //Autosize when img loaded\n",
       "    this.width = 0;\n",
       "    this.height = 0;\n",
       "  }\n",
       "\n",
       "  //Process object data\n",
       "  for (var id in this.vis.objects) {\n",
       "    //Apply object bounding box\n",
       "    if (objbb && this.vis.objects[id].min)\n",
       "      this.checkPointMinMax(this.vis.objects[id].min);\n",
       "    if (objbb && this.vis.objects[id].max)\n",
       "      this.checkPointMinMax(this.vis.objects[id].max);\n",
       "  }\n",
       "\n",
       "  this.updateDims(this.vis.views[this.view]);\n",
       "\n",
       "  //Update display\n",
       "  if (!this.gl) return;\n",
       "  this.draw();\n",
       "  this.clear();\n",
       "}\n",
       "\n",
       "\n",
       "Viewer.prototype.clear = function() {\n",
       "  if (!this.gl) return;\n",
       "  this.gl.clear(this.gl.COLOR_BUFFER_BIT | this.gl.DEPTH_BUFFER_BIT);\n",
       "}\n",
       "\n",
       "Viewer.prototype.draw = function() {\n",
       "  if (!this.canvas) return;\n",
       "\n",
       "  //Get the dimensions from the current canvas\n",
       "  if (this.width != this.canvas.offsetWidth || this.height != this.canvas.offsetHeight) {\n",
       "    this.width = this.canvas.offsetWidth;\n",
       "    this.height = this.canvas.offsetHeight;\n",
       "    //Need to set this too for some reason\n",
       "    this.canvas.width = this.width;\n",
       "    this.canvas.height = this.height;\n",
       "    if (this.gl) {\n",
       "      this.gl.viewportWidth = this.width;\n",
       "      this.gl.viewportHeight = this.height;\n",
       "      this.webgl.viewport = new Viewport(0, 0, this.width, this.height);\n",
       "    }\n",
       "  }\n",
       "  if (!this.gl) return;\n",
       "\n",
       "  this.gl.viewport(0, 0, this.gl.viewportWidth, this.gl.viewportHeight);\n",
       "  //this.gl.clearColor(1, 1, 1, 0);\n",
       "  this.gl.clearColor(0.5, 0.5, 0.5, 0);\n",
       "  this.gl.clear(this.gl.COLOR_BUFFER_BIT | this.gl.DEPTH_BUFFER_BIT);\n",
       "\n",
       "  this.webgl.view(this);\n",
       "\n",
       "  //Render objects\n",
       "  this.border.draw(this.webgl);\n",
       "\n",
       "}\n",
       "\n",
       "Viewer.prototype.rotateX = function(deg) {\n",
       "  this.rotation(deg, [1,0,0]);\n",
       "}\n",
       "\n",
       "Viewer.prototype.rotateY = function(deg) {\n",
       "  this.rotation(deg, [0,1,0]);\n",
       "}\n",
       "\n",
       "Viewer.prototype.rotateZ = function(deg) {\n",
       "  this.rotation(deg, [0,0,1]);\n",
       "}\n",
       "\n",
       "Viewer.prototype.rotation = function(deg, axis) {\n",
       "  //Quaterion rotate\n",
       "  var arad = deg * Math.PI / 180.0;\n",
       "  var rotation = quat4.fromAngleAxis(arad, axis);\n",
       "  rotation = quat4.normalize(rotation);\n",
       "  this.rotate = quat4.multiply(rotation, this.rotate);\n",
       "}\n",
       "\n",
       "Viewer.prototype.getRotation = function() {\n",
       "  return [this.rotate[0], this.rotate[1], this.rotate[2], this.rotate[3]];\n",
       "}\n",
       "\n",
       "Viewer.prototype.getRotationString = function() {\n",
       "  //Return current rotation quaternion as string\n",
       "  var q = this.getRotation();\n",
       "  return 'rotation ' + q[0] + ' ' + q[1] + ' ' + q[2] + ' ' + q[3];\n",
       "}\n",
       "\n",
       "Viewer.prototype.getTranslationString = function() {\n",
       "  return 'translation ' + this.translate[0] + ' ' + this.translate[1] + ' ' + this.translate[2];\n",
       "}\n",
       "\n",
       "Viewer.prototype.reset = function() {\n",
       "  if (this.gl) {\n",
       "    this.updateDims(this.vis.views[this.view]);\n",
       "    this.draw();\n",
       "  }\n",
       "\n",
       "  this.command('reset');\n",
       "}\n",
       "\n",
       "Viewer.prototype.zoom = function(factor) {\n",
       "  if (this.gl) {\n",
       "    this.translate[2] += factor * this.modelsize;\n",
       "    this.draw();\n",
       "  }\n",
       "\n",
       "  this.command('' + this.getTranslationString());\n",
       "  //this.command('zoom ' + factor);\n",
       "}\n",
       "\n",
       "Viewer.prototype.zoomClip = function(factor) {\n",
       "  if (this.gl) {\n",
       "     var near_clip = this.near_clip + factor * this.modelsize;\n",
       "     if (near_clip >= this.modelsize * 0.001)\n",
       "       this.near_clip = near_clip;\n",
       "    this.draw();\n",
       "  }\n",
       "\n",
       "  this.command('zoomclip ' + factor);\n",
       "}\n",
       "\n",
       "Viewer.prototype.updateDims = function(view) {\n",
       "  if (!view) return;\n",
       "  var oldsize = this.modelsize;\n",
       "  this.dims = [view.max[0] - view.min[0], view.max[1] - view.min[1], view.max[2] - view.min[2]];\n",
       "  this.modelsize = Math.sqrt(this.dims[0]*this.dims[0] + this.dims[1]*this.dims[1] + this.dims[2]*this.dims[2]);\n",
       "\n",
       "  this.focus = [view.min[0] + 0.5*this.dims[0], view.min[1] + 0.5*this.dims[1], view.min[2] + 0.5*this.dims[2]];\n",
       "  this.centre = [this.focus[0],this.focus[1],this.focus[2]];\n",
       "\n",
       "  this.translate = [0,0,0];\n",
       "  if (this.modelsize != oldsize) this.translate[2] = -this.modelsize*1.25;\n",
       "\n",
       "  if (this.near_clip == 0.0) this.near_clip = this.modelsize / 10.0;   \n",
       "  if (this.far_clip == 0.0) this.far_clip = this.modelsize * 10.0;\n",
       "\n",
       "  quat4.identity(this.rotate);\n",
       "\n",
       "  if (view) {\n",
       "    //Initial rotation\n",
       "    if (view.rotate) {\n",
       "      if (view.rotate.length == 3) {\n",
       "        this.rotateZ(-view.rotate[2]);\n",
       "        this.rotateY(-view.rotate[1]);\n",
       "        this.rotateX(-view.rotate[0]);\n",
       "      } else if (view.rotate.length == 4) {\n",
       "        this.rotate = quat4.create(view.rotate);\n",
       "      }\n",
       "    }\n",
       "\n",
       "    //Translate\n",
       "    if (view.translate) {\n",
       "      this.translate[0] = view.translate[0];\n",
       "      this.translate[1] = view.translate[1];\n",
       "      this.translate[2] = view.translate[2];\n",
       "    }\n",
       "\n",
       "    //Scale\n",
       "    if (view.scale) {\n",
       "      this.scale[0] = view.scale[0];\n",
       "      this.scale[1] = view.scale[1];\n",
       "      this.scale[2] = view.scale[2];\n",
       "    }\n",
       "\n",
       "    //Focal point\n",
       "    if (view.focus) {\n",
       "      this.focus[0] = this.centre[0] = view.focus[0];\n",
       "      this.focus[1] = this.centre[1] = view.focus[1];\n",
       "      this.focus[2] = this.centre[2] = view.focus[2];\n",
       "    }\n",
       "\n",
       "  }\n",
       "\n",
       "  //console.log(\"DIMS: \" + min[0] + \" to \" + max[0] + \",\" + min[1] + \" to \" + max[1] + \",\" + min[2] + \" to \" + max[2]);\n",
       "  //console.log(\"New model size: \" + this.modelsize + \", Focal point: \" + this.focus[0] + \",\" + this.focus[1] + \",\" + this.focus[2]);\n",
       "  //console.log(\"Translate: \" + this.translate[0] + \",\" + this.translate[1] + \",\" + this.translate[2]);\n",
       "\n",
       "  if (!this.gl) return;\n",
       " \n",
       "  //Create the bounding box vertex buffer\n",
       "  this.border.updateBuffers(this.vis.views[this.view]);\n",
       "}\n",
       "\n",
       "var kernel;\n",
       "//Maintain a list of interactor instances opened by id\n",
       "var instances = [];\n",
       "var delay = 0;\n",
       "\n",
       "function initLoad() {\n",
       "  if (parent.IPython) {\n",
       "    kernel = parent.IPython.notebook.kernel;\n",
       "    //This ensures we have imported the lavavu module in notebook scope\n",
       "    //required for control event pass-back\n",
       "    //(NOTE: this looks a bit hacky because it is, but it seems kernel.execute loops only work as a one liners)\n",
       "    kernel.execute('lavavu = None if not \"lavavu\" in globals() else lavavu');\n",
       "    kernel.execute('import sys');\n",
       "    kernel.execute('modules = dict(sys.modules)');\n",
       "    kernel.execute('for m in modules: lavavu = modules[m].lavavu if \"lavavu\" in dir(modules[m]) else lavavu');\n",
       "  }\n",
       "}\n",
       "\n",
       "function getUrl() {\n",
       "  var loc = window.location;\n",
       "  var baseUrl = loc.protocol + \"//\" + loc.hostname + (loc.port ? \":\" + loc.port : \"\");\n",
       "  return baseUrl;\n",
       "}\n",
       "\n",
       "function do_action(id, val, element) {\n",
       "  //Get instance id from parent element, next level up or just use 0\n",
       "  //console.log(\"DO_ACTION \" + id + \" ( \" + val + \" ) \" + element.id);\n",
       "  var inst_id = element.parentElement.dataset.id;\n",
       "  if (inst_id == undefined)\n",
       "    inst_id = element.parentElement.parentElement.dataset.id;\n",
       "  if (inst_id == undefined)\n",
       "    inst_id = 0;\n",
       "  else\n",
       "    inst_id = parseInt(inst_id);\n",
       "\n",
       "  //if (inst_id-1 > instances.length) return;\n",
       "  if (inst_id+1 > instances.length) {console.log(\"Invalid ID for action instance: \" + (inst_id+1) + ' > ' + instances.length); return;}\n",
       "  if (!instances[inst_id]) {console.log(\"NO ACTION INSTANCE: [\" + inst_id + \"] \" + instances[inst_id]); return;}\n",
       "\n",
       "  //Call action on interactor which owns control\n",
       "  instances[inst_id].do_action(id, val);\n",
       "}\n",
       "\n",
       "function redisplay(id) {\n",
       "  //If instance doesn't exist yet, try again in 0.5 seconds... only seems to be needed in chrome\n",
       "  if (id+1 > instances.length || !instances[id] || instances[id].box === undefined) {\n",
       "    if (delay > 5) return; //Give up after 5 tries\n",
       "    console.log(\"Delaying redisplay \" + id);\n",
       "    setTimeout(function() {redisplay(id);}, 500); \n",
       "    delay++;\n",
       "    if (id+1 > instances.length) console.log(\" -- Invalid ID \" + (id+1) + ' > ' + instances.length);\n",
       "    else if (!instances[id]) console.log(\" -- No instance: [\" + id + \"] \" + instances[id]);\n",
       "    else if (!instances[id].box) console.log(\" -- No box viewer: [\" + id + \"] \" + instances[id]);\n",
       "    return;\n",
       "  }\n",
       "\n",
       "  var that = instances[id];\n",
       "  if (that.img) {\n",
       "    //Call get_image\n",
       "    that.get_image();\n",
       "    //Update the box size by getting state\n",
       "    updateBox(that.box, function(onget) {that.get_state(onget);});\n",
       "  }\n",
       "}\n",
       "\n",
       "function WindowInteractor(id) {\n",
       "  //Always check for lavavu module loaded, in case IPython notebook has been reloaded\n",
       "  initLoad();\n",
       "  //Check instances list valid, wipe if adding initial interactor\n",
       "  if (!instances || id==0) instances = [];\n",
       "  //Store self in list and save id\n",
       "  this.id = id;\n",
       "  instances[id] = this;\n",
       "  \n",
       "  console.log(\"New interactor: \" + this.id + \" : \" + instances.length + \" :: \" + instances[this.id]);\n",
       "  //Interactor class, handles javascript side of window control\n",
       "  //Takes viewer id\n",
       "  // - set as active target for commands\n",
       "  // - init webgl bounding box\n",
       "  // - display initial image\n",
       "  this.img = document.getElementById(\"imgtarget_\" + this.id);\n",
       "\n",
       "  //No window\n",
       "  if (!this.img) return;\n",
       "\n",
       "  //Initial image\n",
       "  //(Init WebGL bounding box interaction on load)\n",
       "  var that = this;\n",
       "  this.get_image(function() {\n",
       "    //console.log('In image loaded callback 0');\n",
       "    //Init on image load with callback function to execute commands\n",
       "    that.box = initBox(that.img, function(cmd) {that.execute(cmd, true);});\n",
       "    //console.log(\"Box init on \" + that.id);\n",
       "    //Update the box size by getting state\n",
       "    updateBox(that.box, function(onget) {that.get_state(onget);});\n",
       "    //Clear onload\n",
       "    that.img.onload = null;\n",
       "  });\n",
       "}\n",
       "\n",
       "WindowInteractor.prototype.execute = function(cmd, instant) { //, viewer_id) {\n",
       "  //console.log(\"execute: \" + cmd);\n",
       "  if (kernel) {\n",
       "    kernel.execute('lavavu.control.windows[' + this.id + '].commands(\"' + cmd + '\")');\n",
       "    this.get_image();\n",
       "  } else {\n",
       "    //HTTP interface\n",
       "    //Replace newlines with semi-colon first\n",
       "    cmd = cmd.replace(/\\n/g,';');\n",
       "    if (instant && this.img) {\n",
       "      var url = getUrl() + \"/icommand=\" + cmd + \"?\" + new Date().getTime();\n",
       "      this.img.onload = null;\n",
       "      this.img.src = url;\n",
       "    } else {\n",
       "      var url = getUrl() + \"/command=\" + cmd + \"?\" + new Date().getTime()\n",
       "      x = new XMLHttpRequest();\n",
       "      x.open('GET', url, true);\n",
       "      x.send();\n",
       "      this.get_image();\n",
       "    }\n",
       "\n",
       "  }\n",
       "  return false;\n",
       "}\n",
       "\n",
       "WindowInteractor.prototype.set_prop = function(obj, prop, val) {\n",
       "  this.execute(\"select \" + obj + \"; \" + prop + \"=\" + val);\n",
       "  this.get_image();\n",
       "}\n",
       "\n",
       "WindowInteractor.prototype.do_action = function(id, val) {\n",
       "  if (kernel) {\n",
       "    //Non-numeric, add quotes and strip newlines\n",
       "    if (typeof(val) == 'string' && (\"\" + parseFloat(val) !== val)) {\n",
       "      val = val.replace(/\\n/g,';');\n",
       "      val = '\"' + val + '\"';\n",
       "    }\n",
       "    kernel.execute('cmds = lavavu.control.action(' + id + ',' + val + ')');\n",
       "    kernel.execute('if len(cmds): lavavu.control.windows[' + this.id + '].commands(cmds)');\n",
       "    /*/Debug callback\n",
       "    var callbacks = {'output' : function(out) {\n",
       "        if (!out.content.data) {alert(\"NO DATA\"); return;}\n",
       "        data = out.content.data['text/plain']\n",
       "          alert(\"RESULT: \" + data);\n",
       "      }\n",
       "    };*/\n",
       "    //kernel.execute('len(lavavu.control.actions)', {iopub: callbacks}, {silent: false});\n",
       "    //kernel.execute('str(lavavu.control.actions[' + id + '])', {iopub: callbacks}, {silent: false});\n",
       "    //kernel.execute('lavavu.control.action(' + id + ',' + val + ')', {iopub: callbacks}, {silent: false});\n",
       "    //kernel.execute('str(lavavu.control.actions[' + id + '][\"args\"][0])', {iopub: callbacks}, {silent: false});\n",
       "    //kernel.execute('cmds', {iopub: callbacks}, {silent: false});\n",
       "    this.get_image();\n",
       "  } else {\n",
       "    //HTML control actions via http\n",
       "    actions[id](val);\n",
       "  }\n",
       "\n",
       "  //Reload state\n",
       "  var that = this;\n",
       "  if (that.img)\n",
       "    updateBox(this.box, function(onget) {that.get_state(onget);});\n",
       "}\n",
       "\n",
       "WindowInteractor.prototype.get_image = function(onload) {\n",
       "  if (!this.img) return;\n",
       "  //console.log(\"get_img: \" + this.id);\n",
       "  if (kernel) {\n",
       "    var that = this;\n",
       "    var callbacks = {'output' : function(out) {\n",
       "        //Skip first message we get without data\n",
       "        //if (!out.content.data) {console.log(JSON.stringify(out)); return;}\n",
       "        if (!out.content.data) return;\n",
       "        data = out.content.data['text/plain']\n",
       "        data = data.substring(1, data.length-1)\n",
       "        //console.log(\"Got image: \" + data.length);\n",
       "        if (that.img) {\n",
       "          that.img.onload = onload;\n",
       "          that.img.src = data;\n",
       "        }\n",
       "      }\n",
       "    };\n",
       "\n",
       "    kernel.execute('lavavu.control.windows[' + this.id + '].frame()', {iopub: callbacks}, {silent: false});\n",
       "  } else {\n",
       "    //if (!this.img) this.img = document.getElementById('imgtarget_0');\n",
       "    var url = getUrl() + \"/image?\" + new Date().getTime();\n",
       "    if (this.img) {\n",
       "      this.img.onload = onload;\n",
       "      this.img.src = url;\n",
       "    }\n",
       "  }\n",
       "}\n",
       "\n",
       "WindowInteractor.prototype.get_state = function(onget) {\n",
       "  //console.log(\"get_state: \" + this.id);\n",
       "  if (kernel) {\n",
       "    var callbacks = {'output' : function(out) {\n",
       "        if (!out.content.data) return;\n",
       "        data = out.content.data['text/plain']\n",
       "        data = data.replace(/(?:\\\\n)+/g, \"\");\n",
       "        data = data.substring(1, data.length-1)\n",
       "        onget(data);\n",
       "      }\n",
       "    };\n",
       "    kernel.execute('lavavu.control.windows[' + this.id + '].getState()', {iopub: callbacks}, {silent: false});\n",
       "  }\n",
       "  else {\n",
       "    var url = getUrl() + \"/getstate\"\n",
       "    x = new XMLHttpRequest();\n",
       "    x.onload = function() { \n",
       "      if(x.status == 200)\n",
       "        onget(x.response);\n",
       "      else  \n",
       "        console.log(\"Ajax Request Error: \" + url + \", returned status code \" + x.status + \" \" + x.statusText);\n",
       "    } \n",
       "    x.open('GET', url, true);\n",
       "    x.send();\n",
       "  }\n",
       "}\n",
       "\n",
       "\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "step = 1\n",
      "step = 2\n",
      "step = 3\n",
      "step = 4\n",
      "step = 5\n",
      "step = 6\n",
      "step = 7\n",
      "step = 8\n",
      "step = 9\n",
      "step = 10\n",
      "step = 11\n",
      "step = 12"
     ]
    }
   ],
   "source": [
    "while step < md.maxSteps:\n",
    "    \n",
    "    solver.solve(nonLinearIterate=True, nonLinearTolerance=md.nltol)\n",
    "    \n",
    "    # main\n",
    "    time,step,files_this_step, next_image_step = main_update(next_image_step)\n",
    "    \n",
    "    #markers / markerLines\n",
    "    if step % md.swarmUpdate == 0:\n",
    "        markerLine_update()    \n",
    "    \n",
    "    #checkpoint\n",
    "    if step % md.checkpointEvery == 0:\n",
    "        cp.saveObjs(step, time)\n",
    "        cp.saveDicts(step, time)\n",
    "        \n",
    "    #particles\n",
    "    if step % md.swarmUpdate == 0:\n",
    "        swarm_update()\n",
    "        \n",
    "        \n",
    "    #Viz\n",
    "    if files_this_step:\n",
    "        viz_update()    \n",
    "        \n",
    "    #xdmfs\n",
    "    if files_this_step:\n",
    "        xdmfs_update()\n",
    "    \n",
    "        \n",
    "    print 'step =',step\n",
    "    \n",
    "\n",
    "    \n",
    "\n",
    "print 'step =',step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#Done!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  },
  "widgets": {
   "state": {},
   "version": "1.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}