fract_macrosystems.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "150 6.0\n",
      "150 6.016666666666667\n",
      "150 6.033333333333333\n",
      "150 6.05\n",
      "150 6.066666666666666\n",
      "150 6.083333333333333\n",
      "150 6.1\n",
      "150 6.116666666666666\n",
      "150 6.133333333333333\n",
      "150 6.1499999999999995\n",
      "150 6.166666666666666\n",
      "150 6.183333333333333\n",
      "150 6.199999999999999\n",
      "150 6.216666666666666\n",
      "150 6.2333333333333325\n",
      "150 6.249999999999999\n",
      "150 6.266666666666666\n",
      "150 6.283333333333332\n",
      "150 6.299999999999999\n",
      "150 6.3166666666666655\n",
      "150 6.333333333333332\n",
      "150 6.349999999999999\n",
      "150 6.366666666666665\n",
      "150 6.383333333333332\n",
      "150 6.399999999999999\n",
      "150 6.416666666666665\n",
      "150 6.433333333333332\n",
      "150 6.449999999999998\n",
      "150 6.466666666666665\n",
      "150 6.483333333333332\n",
      "150 6.499999999999998\n",
      "150 6.516666666666665\n",
      "150 6.533333333333331\n",
      "150 6.549999999999998\n",
      "150 6.566666666666665\n",
      "150 6.583333333333331\n",
      "150 6.599999999999998\n",
      "150 6.6166666666666645\n",
      "150 6.633333333333331\n",
      "150 6.649999999999998\n",
      "150 6.666666666666664\n",
      "150 6.683333333333331\n",
      "150 6.6999999999999975\n",
      "150 6.716666666666664\n",
      "150 6.733333333333331\n",
      "150 6.749999999999997\n",
      "150 6.766666666666664\n",
      "150 6.7833333333333306\n",
      "150 6.799999999999997\n",
      "150 6.816666666666664\n",
      "150 6.83333333333333\n",
      "150 6.849999999999997\n",
      "150 6.866666666666664\n",
      "150 6.88333333333333\n",
      "150 6.899999999999997\n",
      "150 6.916666666666663\n",
      "150 6.93333333333333\n",
      "150 6.949999999999997\n",
      "150 6.966666666666663\n",
      "150 6.98333333333333\n",
      "150 6.9999999999999964\n",
      "150 7.016666666666663\n",
      "150 7.03333333333333\n",
      "150 7.049999999999996\n",
      "150 7.066666666666663\n",
      "150 7.0833333333333295\n",
      "150 7.099999999999996\n",
      "150 7.116666666666663\n",
      "150 7.133333333333329\n",
      "150 7.149999999999996\n",
      "150 7.1666666666666625\n",
      "150 7.183333333333329\n",
      "150 7.199999999999996\n",
      "150 7.216666666666662\n",
      "150 7.233333333333329\n",
      "150 7.249999999999996\n",
      "150 7.266666666666662\n",
      "150 7.283333333333329\n",
      "150 7.299999999999995\n",
      "150 7.316666666666662\n",
      "150 7.333333333333329\n",
      "150 7.349999999999995\n",
      "150 7.366666666666662\n",
      "150 7.383333333333328\n",
      "150 7.399999999999995\n",
      "150 7.416666666666662\n",
      "150 7.433333333333328\n",
      "150 7.449999999999995\n",
      "150 7.4666666666666615\n",
      "150 7.483333333333328\n",
      "150 7.499999999999995\n",
      "150 7.516666666666661\n",
      "150 7.533333333333328\n",
      "150 7.5499999999999945\n",
      "150 7.566666666666661\n",
      "150 7.583333333333328\n",
      "150 7.599999999999994\n",
      "150 7.616666666666661\n",
      "150 7.6333333333333275\n",
      "150 7.649999999999994\n",
      "150 7.666666666666661\n",
      "150 7.683333333333327\n",
      "150 7.699999999999994\n",
      "150 7.716666666666661\n",
      "150 7.733333333333327\n",
      "150 7.749999999999994\n",
      "150 7.76666666666666\n",
      "150 7.783333333333327\n",
      "150 7.799999999999994\n",
      "150 7.81666666666666\n",
      "150 7.833333333333327\n",
      "150 7.849999999999993\n",
      "150 7.86666666666666\n",
      "150 7.883333333333327\n",
      "150 7.899999999999993\n",
      "150 7.91666666666666\n",
      "150 7.9333333333333265\n",
      "150 7.949999999999993\n",
      "150 7.96666666666666\n",
      "150 7.983333333333326\n",
      "150 7.999999999999993\n",
      "150 8.016666666666659\n",
      "150 8.033333333333326\n",
      "150 8.049999999999994\n",
      "150 8.06666666666666\n",
      "150 8.083333333333325\n",
      "150 8.099999999999993\n",
      "150 8.11666666666666\n",
      "150 8.133333333333326\n",
      "150 8.149999999999991\n",
      "150 8.166666666666659\n",
      "150 8.183333333333326\n",
      "150 8.199999999999992\n",
      "150 8.216666666666658\n",
      "150 8.233333333333325\n",
      "150 8.249999999999993\n",
      "150 8.266666666666659\n",
      "150 8.283333333333324\n",
      "150 8.299999999999992\n",
      "150 8.31666666666666\n",
      "150 8.333333333333325\n",
      "150 8.34999999999999\n",
      "150 8.366666666666658\n",
      "150 8.383333333333326\n",
      "150 8.399999999999991\n",
      "150 8.416666666666657\n",
      "150 8.433333333333325\n",
      "150 8.449999999999992\n",
      "150 8.466666666666658\n",
      "150 8.483333333333324\n",
      "150 8.499999999999991\n",
      "150 8.516666666666659\n",
      "150 8.533333333333324\n",
      "150 8.54999999999999\n",
      "150 8.566666666666658\n",
      "150 8.583333333333325\n",
      "150 8.59999999999999\n",
      "150 8.616666666666656\n",
      "150 8.633333333333324\n",
      "150 8.649999999999991\n",
      "150 8.666666666666657\n",
      "150 8.683333333333323\n",
      "150 8.69999999999999\n",
      "150 8.716666666666658\n",
      "150 8.733333333333324\n",
      "150 8.74999999999999\n",
      "150 8.766666666666657\n",
      "150 8.783333333333324\n",
      "150 8.79999999999999\n",
      "150 8.816666666666656\n",
      "150 8.833333333333323\n",
      "150 8.84999999999999\n",
      "150 8.866666666666656\n",
      "150 8.883333333333322\n",
      "150 8.89999999999999\n",
      "150 8.916666666666657\n",
      "150 8.933333333333323\n",
      "150 8.949999999999989\n",
      "150 8.966666666666656\n",
      "150 8.983333333333324\n",
      "150 8.99999999999999\n",
      "150 9.016666666666655\n",
      "150 9.033333333333323\n",
      "150 9.04999999999999\n",
      "150 9.066666666666656\n",
      "150 9.083333333333321\n",
      "150 9.099999999999989\n",
      "150 9.116666666666656\n",
      "150 9.133333333333322\n",
      "150 9.149999999999988\n",
      "150 9.166666666666655\n",
      "150 9.183333333333323\n",
      "150 9.199999999999989\n",
      "150 9.216666666666654\n",
      "150 9.233333333333322\n",
      "150 9.24999999999999\n",
      "150 9.266666666666655\n",
      "150 9.28333333333332\n",
      "150 9.299999999999988\n",
      "150 9.316666666666656\n",
      "150 9.333333333333321\n",
      "150 9.349999999999987\n",
      "150 9.366666666666655\n",
      "150 9.383333333333322\n",
      "150 9.399999999999988\n",
      "150 9.416666666666654\n",
      "150 9.433333333333321\n",
      "150 9.449999999999989\n",
      "150 9.466666666666654\n",
      "150 9.48333333333332\n",
      "150 9.499999999999988\n",
      "150 9.516666666666655\n",
      "150 9.53333333333332\n",
      "150 9.549999999999986\n",
      "150 9.566666666666654\n",
      "150 9.583333333333321\n",
      "150 9.599999999999987\n",
      "150 9.616666666666653\n",
      "150 9.63333333333332\n",
      "150 9.649999999999988\n",
      "150 9.666666666666654\n",
      "150 9.68333333333332\n",
      "150 9.699999999999987\n",
      "150 9.716666666666654\n",
      "150 9.73333333333332\n",
      "150 9.749999999999986\n",
      "150 9.766666666666653\n",
      "150 9.78333333333332\n",
      "150 9.799999999999986\n",
      "150 9.816666666666652\n",
      "150 9.83333333333332\n",
      "150 9.849999999999987\n",
      "150 9.866666666666653\n",
      "150 9.883333333333319\n",
      "150 9.899999999999986\n",
      "150 9.916666666666654\n",
      "150 9.93333333333332\n",
      "150 9.949999999999985\n",
      "150 9.966666666666653\n",
      "150 9.98333333333332\n",
      "150 9.999999999999986\n",
      "150 10.016666666666652\n",
      "150 10.033333333333319\n",
      "150 10.049999999999986\n",
      "150 10.066666666666652\n",
      "150 10.083333333333318\n",
      "150 10.099999999999985\n",
      "150 10.116666666666653\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "Cell \u001b[0;32mIn [1], line 1120\u001b[0m\n\u001b[1;32m   1092\u001b[0m \u001b[39m\"\"\"\u001b[39;00m\n\u001b[1;32m   1093\u001b[0m \u001b[39mjobs = []\u001b[39;00m\n\u001b[1;32m   1094\u001b[0m \u001b[39mfor x in range(1000):\u001b[39;00m\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m   1099\u001b[0m \u001b[39m\tjobs[x].join()\u001b[39;00m\n\u001b[1;32m   1100\u001b[0m \u001b[39m\"\"\"\u001b[39;00m\n\u001b[1;32m   1105\u001b[0m \u001b[39m\"\"\"\u001b[39;00m\n\u001b[1;32m   1106\u001b[0m \u001b[39mfrom multiprocessing import Pool\u001b[39;00m\n\u001b[1;32m   1107\u001b[0m \u001b[39mpool = Pool(8)\u001b[39;00m\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m   1116\u001b[0m \n\u001b[1;32m   1117\u001b[0m \u001b[39m\"\"\"\u001b[39;00m\n\u001b[0;32m-> 1120\u001b[0m do_sim() \u001b[39m#to run it one \u001b[39;00m\n\u001b[1;32m   1121\u001b[0m output\u001b[39m.\u001b[39mclose()\n\u001b[1;32m   1122\u001b[0m output2\u001b[39m.\u001b[39mclose()\n",
      "Cell \u001b[0;32mIn [1], line 1058\u001b[0m, in \u001b[0;36mdo_sim\u001b[0;34m()\u001b[0m\n\u001b[1;32m   1056\u001b[0m \u001b[39mfor\u001b[39;00m t \u001b[39min\u001b[39;00m times:\n\u001b[1;32m   1057\u001b[0m \t\u001b[39mprint\u001b[39m(J, t)\n\u001b[0;32m-> 1058\u001b[0m \t[sim_act(lizard,t,J) \u001b[39mfor\u001b[39;00m lizard \u001b[39min\u001b[39;00m population]\n\u001b[1;32m   1060\u001b[0m tot_act, tot_smr, tot_move, tot_ave_ate \u001b[39m=\u001b[39m\u001b[39m0.\u001b[39m,\u001b[39m0.\u001b[39m,\u001b[39m0.\u001b[39m, \u001b[39m0.\u001b[39m\n\u001b[1;32m   1061\u001b[0m \u001b[39mfor\u001b[39;00m lizard \u001b[39min\u001b[39;00m population:\n",
      "Cell \u001b[0;32mIn [1], line 1058\u001b[0m, in \u001b[0;36m<listcomp>\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m   1056\u001b[0m \u001b[39mfor\u001b[39;00m t \u001b[39min\u001b[39;00m times:\n\u001b[1;32m   1057\u001b[0m \t\u001b[39mprint\u001b[39m(J, t)\n\u001b[0;32m-> 1058\u001b[0m \t[sim_act(lizard,t,J) \u001b[39mfor\u001b[39;00m lizard \u001b[39min\u001b[39;00m population]\n\u001b[1;32m   1060\u001b[0m tot_act, tot_smr, tot_move, tot_ave_ate \u001b[39m=\u001b[39m\u001b[39m0.\u001b[39m,\u001b[39m0.\u001b[39m,\u001b[39m0.\u001b[39m, \u001b[39m0.\u001b[39m\n\u001b[1;32m   1061\u001b[0m \u001b[39mfor\u001b[39;00m lizard \u001b[39min\u001b[39;00m population:\n",
      "Cell \u001b[0;32mIn [1], line 1034\u001b[0m, in \u001b[0;36mdo_sim.<locals>.sim_act\u001b[0;34m(lizard, t, J)\u001b[0m\n\u001b[1;32m   1031\u001b[0m \t\u001b[39m#lizard.tot_dist_moved += lizard.moved\u001b[39;00m\n\u001b[1;32m   1033\u001b[0m \u001b[39melif\u001b[39;00m lizard\u001b[39m.\u001b[39mtlow \u001b[39m<\u001b[39m\u001b[39m=\u001b[39m lizard\u001b[39m.\u001b[39mtb \u001b[39m<\u001b[39m\u001b[39m=\u001b[39mlizard\u001b[39m.\u001b[39mctmax:\n\u001b[0;32m-> 1034\u001b[0m \tlizard\u001b[39m.\u001b[39;49mmove_reg(t,J,elev, slope, aspect, veget)\n\u001b[1;32m   1035\u001b[0m \tlizard\u001b[39m.\u001b[39mtb \u001b[39m=\u001b[39m lizard\u001b[39m.\u001b[39mupdate_tb()\n\u001b[1;32m   1036\u001b[0m \tlizard\u001b[39m.\u001b[39menergy_balance \u001b[39m+\u001b[39m\u001b[39m=\u001b[39m lizard\u001b[39m.\u001b[39msmr()\n",
      "Cell \u001b[0;32mIn [1], line 799\u001b[0m, in \u001b[0;36mdo_sim.<locals>.T_e.move_reg\u001b[0;34m(self, t, J, elev, slope, aspect, veget)\u001b[0m\n\u001b[1;32m    797\u001b[0m az\u001b[39m=\u001b[39mradians(\u001b[39mself\u001b[39m\u001b[39m.\u001b[39mazimuth(t, J))\n\u001b[1;32m    798\u001b[0m al\u001b[39m=\u001b[39mradians(\u001b[39mself\u001b[39m\u001b[39m.\u001b[39maltitude(t, J))\n\u001b[0;32m--> 799\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mte \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mt_e(t,J, sun\u001b[39m=\u001b[39mhorizon_master(az, elev[\u001b[39mint\u001b[39;49m(y1)][\u001b[39mint\u001b[39;49m(x1)], x\u001b[39m=\u001b[39;49mx1, y\u001b[39m=\u001b[39;49my1), a_s\u001b[39m=\u001b[39m\u001b[39m90.\u001b[39m\u001b[39m-\u001b[39m slope[\u001b[39mint\u001b[39m(y1)][\u001b[39mint\u001b[39m(x1)], alpha_s\u001b[39m=\u001b[39maspect[\u001b[39mint\u001b[39m(y1)][\u001b[39mint\u001b[39m(x1)], veg\u001b[39m=\u001b[39mveget[\u001b[39mint\u001b[39m(y1)][\u001b[39mint\u001b[39m(x1)])\n\u001b[1;32m    800\u001b[0m up_tb \u001b[39m=\u001b[39m  \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mupdate_tb()\n\u001b[1;32m    801\u001b[0m diff_tb_new \u001b[39m=\u001b[39m \u001b[39mabs\u001b[39m(\u001b[39mself\u001b[39m\u001b[39m.\u001b[39mt_target \u001b[39m-\u001b[39m up_tb)\n",
      "Cell \u001b[0;32mIn [1], line 968\u001b[0m, in \u001b[0;36mdo_sim.<locals>.horizon_master\u001b[0;34m(azimuth, elevation, x, y)\u001b[0m\n\u001b[1;32m    967\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mhorizon_master\u001b[39m(azimuth, elevation, x, y):\n\u001b[0;32m--> 968\u001b[0m \t\u001b[39mif\u001b[39;00m find_horizon(azimuth,x,y) \u001b[39m>\u001b[39m elevation:\n\u001b[1;32m    969\u001b[0m \t\t\u001b[39m#print True\u001b[39;00m\n\u001b[1;32m    970\u001b[0m \t\t\u001b[39mreturn\u001b[39;00m \u001b[39m0\u001b[39m \u001b[39m#True is shaded, so we want '0' for full sun\u001b[39;00m\n\u001b[1;32m    971\u001b[0m \t\u001b[39melse\u001b[39;00m:\n\u001b[1;32m    972\u001b[0m \t\t\u001b[39m#print False\u001b[39;00m\n",
      "Cell \u001b[0;32mIn [1], line 984\u001b[0m, in \u001b[0;36mdo_sim.<locals>.find_horizon\u001b[0;34m(azimuth, x, y)\u001b[0m\n\u001b[1;32m    982\u001b[0m \tx \u001b[39m=\u001b[39m x \u001b[39m+\u001b[39m hyp\u001b[39m*\u001b[39msin(azimuth)\u001b[39m*\u001b[39mcellsize \u001b[39m#added int # check whether this should be mulitplied by cell size\u001b[39;00m\n\u001b[1;32m    983\u001b[0m \ty \u001b[39m=\u001b[39m y \u001b[39m+\u001b[39m hyp\u001b[39m*\u001b[39mcos(azimuth)\u001b[39m*\u001b[39mcellsize\n\u001b[0;32m--> 984\u001b[0m \t\u001b[39mif\u001b[39;00m data[math\u001b[39m.\u001b[39;49mfloor(y)][math\u001b[39m.\u001b[39;49mfloor(x)] \u001b[39mnot\u001b[39;49;00m \u001b[39min\u001b[39;49;00m z_data:\n\u001b[1;32m    985\u001b[0m \t\tz_data\u001b[39m.\u001b[39mappend(data[math\u001b[39m.\u001b[39mfloor(y)][math\u001b[39m.\u001b[39mfloor(x)])\n\u001b[1;32m    986\u001b[0m h_max \u001b[39m=\u001b[39m \u001b[39m-\u001b[39mpi\u001b[39m/\u001b[39m\u001b[39m2.\u001b[39m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "from time import time as what_time\n",
    "from random import choice, gauss\n",
    "import numpy as np\n",
    "from numpy.random import *\n",
    "from scipy.stats.mstats import *\n",
    "from scipy.interpolate import interp1d\n",
    "from numpy.random import uniform, seed, normal, binomial\n",
    "from numpy import pi, min, max, percentile, radians, degrees\n",
    "from numpy.fft import ifft2\n",
    "import multiprocessing\n",
    "from time import time\n",
    "from scipy.interpolate import *\n",
    "#from matplotlib.mlab import prctile as percent\n",
    "from numpy import flipud\n",
    "from math import exp, log, sin, cos, sqrt, acos, asin, atan, atan2\n",
    "#from numpy import arctan as atan\n",
    "#from numpy import arctan2 as atan2\n",
    "#from numpy import arcsin as asin\n",
    "#from numpy import arccos as acos\n",
    "import pandas as pd ## added by LN\n",
    "import math\n",
    "\n",
    "seed() #99\n",
    "\n",
    "h = 0. #input(\"elevational range? \") #10.\n",
    "high = 1\n",
    "#D = 3 #input(\"fractal dimension (between 2 and 3)? \") #1 # D = 3 - H, for 0 < H < 1\n",
    "#cdef double H #= 3. - D\n",
    "N = 100 #input(\"spatial extent? \") #128\n",
    "MU= 0.\n",
    "SD = 1.\n",
    "J =  150 #input(\"day? \")\n",
    "from time import time\n",
    "\n",
    "## read in elev, slope, aspect, veget files\n",
    "## LR2\n",
    "#elev = pd.DataFrame(pd.read_csv('/Users/laurenneel/Desktop/Grass_RGB_DSM_Oct22/LR2/LR2_elevation'))\n",
    "#elev = elev.iloc[8:]\n",
    "#elev.reset_index(drop=True, inplace=True)\n",
    "#\n",
    "#slope = pd.DataFrame(pd.read_csv('/Users/laurenneel/Desktop/Grass_RGB_DSM_Oct22/LR2/LR2_slope'))\n",
    "#slope=slope.iloc[8:]\n",
    "#slope.reset_index(drop=True, inplace=True)\n",
    "#\n",
    "#aspect = pd.DataFrame(pd.read_csv('/Users/laurenneel/Desktop/Grass_RGB_DSM_Oct22/LR2/LR2_aspect'))\n",
    "#aspect = aspect.iloc[8:]\n",
    "#aspect.reset_index(drop=True, inplace=True)\n",
    "#\n",
    "#\n",
    "#veget = pd.DataFrame(pd.read_csv('/Users/laurenneel/Desktop/Grass_RGB_DSM_Oct22/LR2/LR2_VARI'))\n",
    "#veget=veget.iloc[8:]\n",
    "#veget.reset_index(drop=True, inplace=True)\n",
    "\n",
    "def read_it():\n",
    "\t\tx=[]\n",
    "\t\twhile True:\n",
    "\t\t\tline=finput.readline()\n",
    "\t\t\tif not line: break\n",
    "\t\t\tx.append(line.split())\n",
    "\t\tfinput.close()\n",
    "\t\treturn(x)\n",
    "\n",
    "\n",
    "#i=0\n",
    "def char2int(y):\n",
    "\t\ti=0\n",
    "\t\twhile i < len(y):\n",
    "\t\t\tj = 0\n",
    "\t\t\twhile j < len(y[i]):\n",
    "\t\t\t\ty[i][j] = float(y[i][j])\n",
    "\t\t\t\tj = j + 1\n",
    "\t\t\ti = i + 1\n",
    "\t\treturn y\n",
    "\n",
    "path_txt_input_folder='/Users/laurenneel/Desktop/Grass_RGB_DSM_Oct22/LR2/'\n",
    "finput = open(path_txt_input_folder + \"LR2_elevation\", \"r\")\n",
    "ELEV = []\n",
    "ELEV = read_it()\n",
    "HEADER1 = ELEV[:8]          #creates raster header\n",
    "ELEV = ELEV[8:]\n",
    "ELEV=char2int(ELEV)\n",
    "elev = np.array(ELEV, dtype=float)\n",
    "del ELEV\n",
    "\n",
    "finput = open(path_txt_input_folder + \"LR2_slope\", \"r\")\n",
    "SLOPE = []\n",
    "SLOPE= read_it()\n",
    "HEADER2 = SLOPE[:8]          #creates raster header\n",
    "SLOPE = SLOPE[8:]\n",
    "SLOPE = char2int(SLOPE)\n",
    "slope = np.array(SLOPE, dtype=float)\n",
    "del SLOPE\n",
    "\n",
    "finput = open(path_txt_input_folder + \"LR2_VARI\", \"r\")\n",
    "VEG = []\n",
    "VEG= read_it()\n",
    "HEADER4 = VEG[:8]          #creates raster header\n",
    "VEG = VEG[8:]\n",
    "VEG = char2int(VEG)\n",
    "veget = np.array(VEG, dtype=float)\n",
    "del VEG\n",
    "\n",
    "finput = open(path_txt_input_folder + \"LR2_aspect\", \"r\")\n",
    "ASPECT = []\n",
    "ASPECT= read_it()\n",
    "HEADER4 = ASPECT[:8]          #creates raster header\n",
    "ASPECT = ASPECT[8:]\n",
    "ASPECT = char2int(ASPECT)\n",
    "aspect = np.array(ASPECT, dtype=float)\n",
    "del ASPECT\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "def do_sim():\n",
    "\n",
    "\tSIGMA = 5.673 * 10**(-8) # Stephan-Boltmann constant\n",
    "\tSOLAR = 1360. # Solar constant\n",
    "\t#filename = raw_input(\"site? \")\n",
    "\t\n",
    "\t#t =  input(\"time of day? \")\n",
    "\n",
    "\tclimate_change = 0 #input(\"Climate change scenario? 0 = no, 1= yes: \")\n",
    "\tif climate_change==1:\n",
    "\t\tminTemps = np.array([-5., -5.2, -2.7, 1.3, 5.8, 11.5 , 16.2, 18.5, 17.3, 12.7, 6.6, -1.0, -5.0, -5.2]) +3.\n",
    "\t\tmaxTemps = np.array([14.4,15.2, 17.6, 21.4, 26.6, 31.2, 35.6, 36.0, 34.7, 31.2, 25.4, 18.8, 14.4, 15.2]) +3.\n",
    "\t\tg=3.\n",
    "\t\tclimname = \"change\"\n",
    "\telse:\n",
    "\t\tminTemps = np.array([-5.,-5.2, -2.7, 1.3, 5.8, 11.5, 16.2, 18.5, 17.3, 12.7, 6.6, -1.0, -5.0, -5.2]) \n",
    "\t\tmaxTemps = np.array([14.4,15.2, 17.6, 21.4, 26.6, 31.2, 35.6, 36.0, 34.7, 31.2, 25.4, 18.8, 14.4, 15.2])\n",
    "\t\tg=0.\n",
    "\t\tclimname = \"norm\"\n",
    "\n",
    "\t#write_to_file = raw_input(\"Write thermal map to: \")\n",
    "\n",
    "\t# DATA\n",
    "\n",
    "\t# Create air temperatures and soil temperatures\n",
    "\t# from splined climate records\n",
    "\n",
    "\t# Air temperature data from Zion NP, east entrance\n",
    "\n",
    "\tdays = np.array([-16, 15, 46, 75, 106, 136, 167, 197, 228, 259, 289, 320, 350, 381])\n",
    "\n",
    "\tmaxT = UnivariateSpline(days, maxTemps, k=3)\n",
    "\tminT = UnivariateSpline(days, minTemps, k=3)\n",
    "\n",
    "\t# Ground temperature data from splines of raw data from 24 days <0.5 C from average MIN & MAX for that day\n",
    "\thours = np.array([0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23.,24.])\n",
    "\tjanTemps = np.array([-7.01,-7.77,-7.89,-8.28,-8.72,-8.99,-9.32,-9.74,-8.6,-5.61,-1.41,2.74,6.3,8.37,8.65,6.85,3.44,0.28,-2.03,-3.38,-4.68,-5.22,-5.91,-6.82,-7.01]) + g\n",
    "\tfebTemps = np.array([-5.15,-6.06,-6.16,-6.59,-6.46,-7.39,-7.75,-7.71,-5.5,-1.36,3.6,7.96,11.43,13.52,13.71,12.09,8.77,4.29,1.23,-0.84,-1.93,-2.98,-3.64,-4.49,-5.11]) + g\n",
    "\tmarTemps = np.array([-1.24,-2.36,-2.52,-3.05,-2.81,-3.97,-4.31,-3.5,0.33,5.38,11.2,15.95,19.81,22.14,23.03,21.03,17.54,12.34,8,5.08,3.05,1.76,0.62,-0.76,-1.14]) + g\n",
    "\taprTemps = np.array([1.7,0.86,0.09,-0.33,-0.84,-1.12,-1.23,1.34,5.94,11.43,17.02,22.06,25.62,27.8,28.09,26.95,23.81,18.94,13.68,9.61,6.88,5.31,3.89,2.94,1.85]) + g\n",
    "\tmayTemps = np.array([7.78,6.76,6.03,5.11,4.91,4.31,5.1,8.89,13.94,19.55,25.38,30.42,34.1,36.79,37,35.81,32.65,28.4,22.89,17.18,14.25,12.02,10.67,8.97,7.84]) + g\n",
    "\tjunTemps = np.array([11.87,10.92,10.14,9.26,8.96,8.49,9.92,13.89,19.07,24.5,30.68,35.81,39.59,42.44,42.66,41.56,38.49,34.26,28.76,22.88,19.22,16.69,14.81,12.9,12.08]) + g\n",
    "\tjulTemps = np.array([17.73,16.67,15.93,14.97,14.78,14.27,15.26,19.37,24.45,30.41,35.81,40.79,44.37,46.65,47.28,46.63,43.27,38.9,33.85,28.06,24.45,21.91,20.33,18.91,17.92]) + g\n",
    "\taugTemps = np.array([16.17,15.02,14.3,14.42,13.15,12.65,13.28,16.21,21.04,26.26,32.55,37.76,41.68,43.98,44.6,43.14,40.12,35.4,29.63,24.99,21.99,19.85,18.38,17.11,15.59]) + g\n",
    "\tsepTemps = np.array([9.99,9.23,8.59,7.83,7.53,7.07,6.87,8.59,12.79,18.13,23.99,29.24,33.2,36.07,36.17,34.71,31.08,26.37,20.91,17.34,14.67,13.43,12.11,11.29,10.11]) + g\n",
    "\toctTemps = np.array( [1.92,1.03,0.67,0.16,-0.75,-0.73,-1.13,-0.97,2.13,6.87,12.26,17.56,21.76,24.31,25.28,23.18,19.49,14.15,10.33,7.82,6.02,4.76,3.68,2.3,1.99]) + g\n",
    "\tnovTemps = np.array( [-4.09,-3.99,-5.08,-5.51,-5.54,-6.28,-6.62,-7.01,-5.96,-1.85,2.78,7.13,11.2,13.47,13.88,12.1,8.58,4.35,1.83,0.2,-0.61,-1.98,-2.79,-3.79,-4.07]) + g\n",
    "\tdecTemps = np.array([-9.46,-9.49,-10.35,-10.75,-10.47,-11.46,-11.79,-12.02,-11.42,-8.67,-4.58,0.16,3.44,5.67,6.01,4.15,0.66,-2.1,-4.52,-5.83,-6.52,-7.66,-8.35,-9.23,-9.47]) + g\n",
    "\n",
    "\tjanSpline = UnivariateSpline(hours,janTemps,k=3)\n",
    "\tfebSpline = UnivariateSpline(hours,febTemps,k=3)\n",
    "\tmarSpline = UnivariateSpline(hours,marTemps,k=3)\n",
    "\taprSpline = UnivariateSpline(hours,aprTemps,k=3)\n",
    "\tmaySpline = UnivariateSpline(hours,mayTemps,k=3)\n",
    "\tjunSpline = UnivariateSpline(hours,junTemps,k=3)\n",
    "\tjulSpline = UnivariateSpline(hours,julTemps,k=3)\n",
    "\taugSpline = UnivariateSpline(hours,augTemps,k=3)\n",
    "\tsepSpline = UnivariateSpline(hours,sepTemps,k=3)\n",
    "\toctSpline = UnivariateSpline(hours,octTemps,k=3)\n",
    "\tnovSpline = UnivariateSpline(hours,novTemps,k=3)\n",
    "\tdecSpline = UnivariateSpline(hours,decTemps,k=3)\n",
    "\n",
    "\tdef make_elevation(h, high, D, N, MU, SD):\n",
    "\t\tseed(100) # just to get same fractal dimensions <-normally this changes\n",
    "\t\tArand = high\n",
    "\t\ti0, j0 = 0., 0.\n",
    "\t\tH = 3. - D\n",
    "\t\tA = np.zeros((N,N),dtype = complex) #complex number plane\n",
    "\t\tfor i in range(N//2):\n",
    "\t\t\tfor j in range(N//2):\n",
    "\t\t\t\tphase = 2. * pi  * uniform(0,high)/Arand\n",
    "\t\t\t\tif i != 0 or j != 0:\n",
    "\t\t\t\t\trad = pow(i*i + j*j,-(H + 1.)/2.) * normal(MU, SD)\n",
    "\t\t\t\telse:\n",
    "\t\t\t\t\trad = 0.\n",
    "\t\t\t\tA[i][j] = complex(rad * cos(phase), rad * sin(phase))\n",
    "\t\t\t\tif  i == 0.:\n",
    "\t\t\t\t\ti0 = 0.\n",
    "\t\t\t\telse:\n",
    "\t\t\t\t\ti0 = N-i #i added int(N-i)\n",
    "\t\t\t\tif j0 == 0.:\n",
    "\t\t\t\t\tj0 = 0.\n",
    "\t\t\t\telse:\n",
    "\t\t\t\t\tj0 = N-j #i added int(N-i)\n",
    "\t\t\t\tif A[int(i0)][int(j0)] == complex(rad * cos(phase), -rad*sin(phase)):\n",
    "\t\t\t\t\tbreak\n",
    "\t\t\t\n",
    "\t\tA[N//2][0]= complex(A[N//2][0].real, 0)\n",
    "\t\tA[0][N//2]= complex(A[N//2][0].real, 0)\n",
    "\t\tA[N//2][N//2]= complex(A[N//2][0].real, 0)\n",
    "\n",
    "\t\tfor i in range(1, int(N/2 -1)):\n",
    "\t\t\tfor j in range(1, int(N/2 -1)):\n",
    "\t\t\t\tphase = 2.* pi* uniform(0,high)/Arand\n",
    "\t\t\t\trad = pow(i*i + j*j, -(H+1)/2.) * normal(MU,SD)\n",
    "\t\t\t\tA[i][N-j] = complex(rad * cos(phase), rad * sin(phase))\n",
    "\t\t\t\tA[N-i][j] = complex(rad * cos(phase), -rad * sin(phase))\n",
    "\t\tX = ifft2(A).real\n",
    "\t\tX= h * (X - min(X.real)) / (max(X.real) - min(X.real))\n",
    "\t\tz = X.view(float)\n",
    "\t\treturn z\n",
    "\n",
    "\tdef make_temp_vegetation(h, high, D, N, MU, SD):\n",
    "\t\tprint(D)\n",
    "\t\tH = 3. - D\n",
    "\t\ttemp_v = make_elevation(h, high, D, N, MU, SD)\n",
    "\t\treturn temp_v\n",
    "\n",
    "\tdef make_vegetation(temp_v, pcnt):\n",
    "\t\treturn np.array(temp_v > percentile(temp_v, 100.-pcnt))#, int)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\t# e is the center cell to a moving grid.\n",
    "\t# Because we want to preserve the fractal nature of the surface,\n",
    "\t# we must move as if over a tesselation.\n",
    "\t#\n",
    "\t#\ta\tb\tc\n",
    "\t#\n",
    "\t#\td\te\tf\n",
    "\t#\n",
    "\t#\tg\th\ti\n",
    "\n",
    "\tdef calc_slope(elev_map, x, y, N):\n",
    "\t\tif (0<x<N-1) and (0 < y < N-1):\n",
    "\t\t\ta = elev_map[y-1][x-1]\n",
    "\t\t\tb = elev_map[y-1][x]\n",
    "\t\t\tc = elev_map[y-1][x+1]\n",
    "\t\t\td = elev_map[y][x-1]\n",
    "\t\t\te = elev_map[y][x]\n",
    "\t\t\tf = elev_map[y][x+1]\n",
    "\t\t\tg = elev_map[y+1][x-1]\n",
    "\t\t\th = elev_map[y+1][x]\n",
    "\t\t\ti = elev_map[y+1][x+1]\n",
    "\t\tif (x == 0) and (0 < y < N-1):\n",
    "\t\t\ta = elev_map[y-1][N-1]\n",
    "\t\t\tb = elev_map[y-1][x]\n",
    "\t\t\tc = elev_map[y-1][x+1]\n",
    "\t\t\td = elev_map[y][N-1]\n",
    "\t\t\te = elev_map[y][x]\n",
    "\t\t\tf = elev_map[y][x+1]\n",
    "\t\t\tg = elev_map[y+1][N-1]\n",
    "\t\t\th = elev_map[y+1][x]\n",
    "\t\t\ti = elev_map[y+1][x+1]\n",
    "\t\tif (x == 0) and (y ==  N-1):\n",
    "\t\t\ta = elev_map[y-1][N-1]\n",
    "\t\t\tb = elev_map[y-1][x]\n",
    "\t\t\tc = elev_map[y-1][x+1]\n",
    "\t\t\td = elev_map[y][N-1]\n",
    "\t\t\te = elev_map[y][x]\n",
    "\t\t\tf = elev_map[y][x+1]\n",
    "\t\t\tg = elev_map[0][N-1]\n",
    "\t\t\th = elev_map[0][x]\n",
    "\t\t\ti = elev_map[0][x+1]\n",
    "\t\tif (x == 0) and (y == 0):\n",
    "\t\t\ta = elev_map[N-1][N-1]\n",
    "\t\t\tb = elev_map[N-1][x]\n",
    "\t\t\tc = elev_map[N-1][x+1]\n",
    "\t\t\td = elev_map[y][N-1]\n",
    "\t\t\te = elev_map[y][x]\n",
    "\t\t\tf = elev_map[y][x+1]\n",
    "\t\t\tg = elev_map[y+1][N-1]\n",
    "\t\t\th = elev_map[y+1][x]\n",
    "\t\t\ti = elev_map[y+1][x+1]\n",
    "\t\tif (x == N-1) and (0 < y < N-1):\n",
    "\t\t\ta = elev_map[y-1][x-1]\n",
    "\t\t\tb = elev_map[y-1][x]\n",
    "\t\t\tc = elev_map[y-1][0]\n",
    "\t\t\td = elev_map[y][x-1]\n",
    "\t\t\te = elev_map[y][x]\n",
    "\t\t\tf = elev_map[y][0]\n",
    "\t\t\tg = elev_map[y+1][x-1]\n",
    "\t\t\th = elev_map[y+1][x]\n",
    "\t\t\ti = elev_map[y+1][0]\n",
    "\t\tif (x == N-1) and (y ==  N-1):\n",
    "\t\t\ta = elev_map[y-1][x-1]\n",
    "\t\t\tb = elev_map[y-1][x]\n",
    "\t\t\tc = elev_map[y-1][0]\n",
    "\t\t\td = elev_map[y][x-1]\n",
    "\t\t\te = elev_map[y][x]\n",
    "\t\t\tf = elev_map[y][0]\n",
    "\t\t\tg = elev_map[0][x-1]\n",
    "\t\t\th = elev_map[0][x]\n",
    "\t\t\ti = elev_map[0][0]\n",
    "\t\tif (x == N-1) and (y == 0):\n",
    "\t\t\ta = elev_map[N-1][x-1]\n",
    "\t\t\tb = elev_map[N-1][x]\n",
    "\t\t\tc = elev_map[N-1][0]\n",
    "\t\t\td = elev_map[y][x-1]\n",
    "\t\t\te = elev_map[y][x]\n",
    "\t\t\tf = elev_map[y][0]\n",
    "\t\t\tg = elev_map[y+1][x-1]\n",
    "\t\t\th = elev_map[y+1][x]\n",
    "\t\t\ti = elev_map[y+1][0]\n",
    "\t\tif (0 < x < N-1) and (y == 0):\n",
    "\t\t\ta = elev_map[N-1][x-1]\n",
    "\t\t\tb = elev_map[N-1][x]\n",
    "\t\t\tc = elev_map[N-1][x+1]\n",
    "\t\t\td = elev_map[y][x-1]\n",
    "\t\t\te = elev_map[y][x]\n",
    "\t\t\tf = elev_map[y][x+1]\n",
    "\t\t\tg = elev_map[y+1][x-1]\n",
    "\t\t\th = elev_map[y+1][x]\n",
    "\t\t\ti = elev_map[y+1][x+1]\n",
    "\t\tif (0 < x < N-1) and (y == N-1):\n",
    "\t\t\ta = elev_map[y-1][x-1]\n",
    "\t\t\tb = elev_map[y-1][x]\n",
    "\t\t\tc = elev_map[y-1][x+1]\n",
    "\t\t\td = elev_map[y][x-1]\n",
    "\t\t\te = elev_map[y][x]\n",
    "\t\t\tf = elev_map[y][x+1]\n",
    "\t\t\tg = elev_map[0][x-1]\n",
    "\t\t\th = elev_map[0][x]\n",
    "\t\t\ti = elev_map[0][x+1]\n",
    "\t\t\n",
    "\t\tx_cell_size, y_cell_size = 1,1\n",
    "\t\tdzdx = ((c + 2*f + i) - (a + 2*d + g)) / (8 * x_cell_size)\n",
    "\t\tdzdy = ((g + 2*h + i) - (a + 2*b + c)) / (8 * y_cell_size)\n",
    "\t\trise_run = sqrt(dzdx**2 +dzdy**2)\n",
    "\t\tslope_degrees = atan(rise_run) * 57.29578\n",
    "\t\treturn slope_degrees\n",
    "\n",
    "\tdef calc_aspect(elev_map, x, y, N):\n",
    "\t\tif (0<x<N-1) and (0 < y < N-1):\n",
    "\t\t\ta = elev_map[y-1][x-1]\n",
    "\t\t\tb = elev_map[y-1][x]\n",
    "\t\t\tc = elev_map[y-1][x+1]\n",
    "\t\t\td = elev_map[y][x-1]\n",
    "\t\t\te = elev_map[y][x]\n",
    "\t\t\tf = elev_map[y][x+1]\n",
    "\t\t\tg = elev_map[y+1][x-1]\n",
    "\t\t\th = elev_map[y+1][x]\n",
    "\t\t\ti = elev_map[y+1][x+1]\n",
    "\t\tif (x == 0) and (0 < y < N-1):\n",
    "\t\t\ta = elev_map[y-1][N-1]\n",
    "\t\t\tb = elev_map[y-1][x]\n",
    "\t\t\tc = elev_map[y-1][x+1]\n",
    "\t\t\td = elev_map[y][N-1]\n",
    "\t\t\te = elev_map[y][x]\n",
    "\t\t\tf = elev_map[y][x+1]\n",
    "\t\t\tg = elev_map[y+1][N-1]\n",
    "\t\t\th = elev_map[y+1][x]\n",
    "\t\t\ti = elev_map[y+1][x+1]\n",
    "\t\tif (x == 0) and (y ==  N-1):\n",
    "\t\t\ta = elev_map[y-1][N-1]\n",
    "\t\t\tb = elev_map[y-1][x]\n",
    "\t\t\tc = elev_map[y-1][x+1]\n",
    "\t\t\td = elev_map[y][N-1]\n",
    "\t\t\te = elev_map[y][x]\n",
    "\t\t\tf = elev_map[y][x+1]\n",
    "\t\t\tg = elev_map[0][N-1]\n",
    "\t\t\th = elev_map[0][x]\n",
    "\t\t\ti = elev_map[0][x+1]\n",
    "\t\tif (x == 0) and (y == 0):\n",
    "\t\t\ta = elev_map[N-1][N-1]\n",
    "\t\t\tb = elev_map[N-1][x]\n",
    "\t\t\tc = elev_map[N-1][x+1]\n",
    "\t\t\td = elev_map[y][N-1]\n",
    "\t\t\te = elev_map[y][x]\n",
    "\t\t\tf = elev_map[y][x+1]\n",
    "\t\t\tg = elev_map[y+1][N-1]\n",
    "\t\t\th = elev_map[y+1][x]\n",
    "\t\t\ti = elev_map[y+1][x+1]\n",
    "\t\tif (x == N-1) and (0 < y < N-1):\n",
    "\t\t\ta = elev_map[y-1][x-1]\n",
    "\t\t\tb = elev_map[y-1][x]\n",
    "\t\t\tc = elev_map[y-1][0]\n",
    "\t\t\td = elev_map[y][x-1]\n",
    "\t\t\te = elev_map[y][x]\n",
    "\t\t\tf = elev_map[y][0]\n",
    "\t\t\tg = elev_map[y+1][x-1]\n",
    "\t\t\th = elev_map[y+1][x]\n",
    "\t\t\ti = elev_map[y+1][0]\n",
    "\t\tif (x == N-1) and (y ==  N-1):\n",
    "\t\t\ta = elev_map[y-1][x-1]\n",
    "\t\t\tb = elev_map[y-1][x]\n",
    "\t\t\tc = elev_map[y-1][0]\n",
    "\t\t\td = elev_map[y][x-1]\n",
    "\t\t\te = elev_map[y][x]\n",
    "\t\t\tf = elev_map[y][0]\n",
    "\t\t\tg = elev_map[0][x-1]\n",
    "\t\t\th = elev_map[0][x]\n",
    "\t\t\ti = elev_map[0][0]\n",
    "\t\tif (x == N-1) and (y == 0):\n",
    "\t\t\ta = elev_map[N-1][x-1]\n",
    "\t\t\tb = elev_map[N-1][x]\n",
    "\t\t\tc = elev_map[N-1][0]\n",
    "\t\t\td = elev_map[y][x-1]\n",
    "\t\t\te = elev_map[y][x]\n",
    "\t\t\tf = elev_map[y][0]\n",
    "\t\t\tg = elev_map[y+1][x-1]\n",
    "\t\t\th = elev_map[y+1][x]\n",
    "\t\t\ti = elev_map[y+1][0]\n",
    "\t\tif (0 < x < N-1) and (y == 0):\n",
    "\t\t\ta = elev_map[N-1][x-1]\n",
    "\t\t\tb = elev_map[N-1][x]\n",
    "\t\t\tc = elev_map[N-1][x+1]\n",
    "\t\t\td = elev_map[y][x-1]\n",
    "\t\t\te = elev_map[y][x]\n",
    "\t\t\tf = elev_map[y][x+1]\n",
    "\t\t\tg = elev_map[y+1][x-1]\n",
    "\t\t\th = elev_map[y+1][x]\n",
    "\t\t\ti = elev_map[y+1][x+1]\n",
    "\t\tif (0 < x < N-1) and (y == N-1):\n",
    "\t\t\ta = elev_map[y-1][x-1]\n",
    "\t\t\tb = elev_map[y-1][x]\n",
    "\t\t\tc = elev_map[y-1][x+1]\n",
    "\t\t\td = elev_map[y][x-1]\n",
    "\t\t\te = elev_map[y][x]\n",
    "\t\t\tf = elev_map[y][x+1]\n",
    "\t\t\tg = elev_map[0][x-1]\n",
    "\t\t\th = elev_map[0][x]\n",
    "\t\t\ti = elev_map[0][x+1]\n",
    "\t\t\t\t\n",
    "\t\tdzdx = ((c + 2*f + i) - (a + 2*d + g)) / 8\n",
    "\t\tdzdy = ((g + 2*h + i) - (a + 2*b + c))  / 8\n",
    "\t\taspect = 57.29578 * atan2(dzdy, -dzdx)\n",
    "\t\n",
    "\t\tif aspect < 0:\n",
    "\t\t\tcell = 90.0 - aspect\n",
    "\t\telif aspect > 90.0:\n",
    "\t\t\tcell = 360.0 - aspect + 90.0\n",
    "\t\telse:\n",
    "\t\t\tcell = 90.0 - aspect\n",
    "\t\treturn cell\n",
    "\n",
    "\tdef make_maps(h,  high,  N,  MU,  SD):\n",
    "\t\tseed()\n",
    "\t\telev_range = 4. #uniform(0, h)\n",
    "\t\telev_fract = 2.2 #input(\"elev fract? \") #uniform(2,3)\n",
    "\t\tveg_range = uniform(0,100) #float(input(\"veg %? \")) #uniform(0,100)\n",
    "\t\tveg_fract = uniform(2,3) #float(input(\"veg fract? \")) #uniform(2,3)\n",
    "\t\t#H_elev = 3. - elev_fract\n",
    "\t\t#H_veg = 3. - veg_fract\n",
    "\t\tprint(elev_range, elev_fract, veg_range, veg_fract)\n",
    "\t\t#X = make_elevation(h, high, D, H, N, MU, SD)\n",
    "\t\tX = make_elevation(elev_range, high, elev_fract, N, MU, SD)\n",
    "\t\t_slope = np.zeros((N,N),dtype = float) #[range(len(X)) for i in range(len(X[0]))]\n",
    "\t\t_aspect = np.zeros((N,N),dtype = float) #[range(len(X)) for i in range(len(X[0]))]\n",
    "\t\n",
    "\t\tfor y in range(len(X)):\n",
    "\t\t\tfor x in range(len(X[y])):\n",
    "\t\t\t\t_slope[y][x] = calc_slope(X, x, y, N)\n",
    "\t\t\t\t_aspect[y][x] = calc_aspect(X, x, y, N)\n",
    "\t\n",
    "\t\tS = _slope #possibly read-in slope here and \n",
    "\t\tA = _aspect\n",
    "\t\ttV = make_temp_vegetation(1., high, veg_fract, N, MU, SD) # h needs to be veg specific\n",
    "\t\tV = make_vegetation(tV, veg_range)\n",
    "\t\n",
    "\t\treturn X, S, A, V, elev_range, elev_fract, veg_range, veg_fract\n",
    "\n",
    "\t#get_topo_data(name,limit)\n",
    "\t#print elev_header[4][1]\n",
    "\n",
    "\n",
    "\tclass T_e():\n",
    "\t\tdef __init__(self, emissivity = 0.95, phi = 34., longitude = 74., elev = 0., h = 0.07, d = 0.2, emmiss_animal = 0.95, sky = 0.7, r_soil = 0.3, emmiss_soil = 0.95, abs_animal_long = 0.97, abs_animal_short = 0.7, Ta_min0 = 20., Ta_min1 = 20., Ta_min2 = 20., Ta_max0 = 35., Ta_max1 = 35., Ta_max2 = 35., ampl = 15., wind = 0.1):\n",
    "\t\t\tself.SIGMA =  5.673e-8\n",
    "\t\t\tself.SOLAR = 1360.\n",
    "\t\t\tself.emissivity = emissivity\n",
    "\t\t\tself.phi = phi\n",
    "\t\t\tself.longitude = longitude\n",
    "\t\t\tself.elev = elev\n",
    "\t\t\tself.h = h\n",
    "\t\t\tself.d = d\n",
    "\t\t\tself.emmiss_animal = emmiss_animal\n",
    "\t\t\tself.char_dim = d\n",
    "\t\t\tself.sky = sky\n",
    "\t\t\tself.r_soil = r_soil\n",
    "\t\t\tself.emmiss_soil = emmiss_soil\n",
    "\t\t\tself.abs_animal_long = abs_animal_long\n",
    "\t\t\tself.abs_animal_short = abs_animal_short\n",
    "\t\t\tself.Ta_min0 = Ta_min0\n",
    "\t\t\tself.Ta_min1 = Ta_min1\n",
    "\t\t\tself.Ta_min2 = Ta_min2\n",
    "\t\t\tself.Ta_max0 = Ta_max0\n",
    "\t\t\tself.Ta_max1 = Ta_max1\n",
    "\t\t\tself.Ta_max2 = Ta_max2\n",
    "\t\t\tself.ampl = ampl\n",
    "\t\t\tself.wind = wind\n",
    "\t\t\tself.size = 10.\n",
    "\t\t\tself.h_t = exp(log(self.size) * 0.36 + 0.72)\n",
    "\t\t\tself.c_t = exp(log(0.42 + 0.44 * self.size))\n",
    "\t\t\tself.Z=0.\n",
    "\t\t\tself.TA= 0.\n",
    "\t\t\n",
    "\t\t\tself.x_min, self.y_min = 0., 0.\n",
    "\t\t\tself.x_max, self.y_max = 99., 99. #x_max, y_max\n",
    "\t\t\tself.position = {'x':rand() * self.x_max, 'y': rand() * self.y_max}\n",
    "\t\t\n",
    "\t\t\tself.thigh =  35.\n",
    "\t\t\tself.tlow = 29.\n",
    "\t\t\tself.topt = 32.\n",
    "\t\t\tself.mu = 0.\n",
    "\t\t\tself.kappa = 0.25\n",
    "\t\t\tself.beta_alpha = 1.\n",
    "\t\t\tself.beta_beta = 3.\n",
    "\t\t\tself.d_max = 5.*(np.log10(self.size)+1.)\n",
    "\t\t\tself.p1 = 0.5\n",
    "\t\t\tself.p2 = 0.0\n",
    "\t\t\tself.p3 = 0.0\n",
    "\t\t\tself.decisions = 6\n",
    "\t\t\tself.tb = 20.\n",
    "\t\t\tself.te = 20.\n",
    "\t\t\tself.orientation = vonmises(self.mu, self.kappa)\n",
    "\t\t\tself.dist = 0.\n",
    "\t\t\tself.ctmax = 39.\n",
    "\t\t\tself.moved = 0\n",
    "\t\t\tself.cost = 0\n",
    "\t\t\tself.totenergy = 0.\n",
    "\t\t\tself.dist = 0.\n",
    "\t\t\tself.active = 0\n",
    "\t\t\tself.t_target = 32.\n",
    "\t\t\tself.db = 0.\n",
    "\t\t\tself.de = 0.\n",
    "\t\t\tself.energy_balance = 0.\n",
    "\t\t\tself.moved = 0.\n",
    "\t\t\tself.tot_dist_moved = 0.\n",
    "\t\t\tself.total_activity = 0.\n",
    "\t\t\tself.ate = 0.\n",
    "\t\n",
    "\t\tdef tprops(self): #thermal preoperties\n",
    "\t\t\tself.h_t = exp(log(self.size) * 0.36 + 0.72)\n",
    "\t\t\tself.c_t = exp(0.42 + 0.44 * log(self.size))\n",
    "\n",
    "\t\tdef greybody(self, K):\n",
    "\t\t\treturn self.emissivity * self.SIGMA * (K + 273)**4\n",
    "\t\n",
    "        ## original MS fxn\n",
    "\t\t#def zenith(self ,t, J):\t# phi = lattitde, delta = solar declination, t0 = solar noon\n",
    "\t\t#\tself.Z = degrees(acos(sin(radians(self.phi)) * sin(radians(self.declination(J))) + cos(radians(self.phi)) * cos(radians(self.declination(J))) * cos(radians(15. * (t - self.noon(J))))))\n",
    "\t\t#\tif self.Z <0.: #or self.Z is None:  ##LN modified \n",
    "\t\t#\t\tself.Z = 0.\n",
    "\t\t#\tif self.Z >90: # or self.Z is None:  ##LN modified \n",
    "\t\t#\t\tself.Z = 90.\n",
    "\t\t#\treturn\n",
    "\n",
    "\t\tdef zenith(self ,t, J):\t# phi = lattitde, delta = solar declination, t0 = solar noon\n",
    "\t\t\tself.Z = degrees(acos(sin(radians(self.phi)) * sin(radians(self.declination(J))) + cos(radians(self.phi)) * cos(radians(self.declination(J))) * cos(radians(15. * (t - self.noon(J))))))\n",
    "\t\t\tif self.Z <0. or self.Z is None:  ##LN modified \n",
    "\t\t\t\tself.Z = 0.\n",
    "\t\t\tif self.Z >90 or self.Z is None:  ##LN modified \n",
    "\t\t\t\tself.Z = 90.\n",
    "\t\t\treturn self.Z # adding return statement so it always returns a value preventing 'NoneType\" error\n",
    "\t\n",
    "\t\tdef altitude(self, t, J):\n",
    "\t\t\treturn degrees(asin(sin(radians(self.phi)) * sin(radians(self.declination(J))) + cos(radians(self.phi)) * cos(radians(self.declination(J))) * cos(radians(15. * (t - self.noon(J))))))\n",
    "\t\t\n",
    "\t\tdef azimuth(self, t, J):\n",
    "\t\t\t#return degrees(asin(-cos(radians(self.declination(J))) * sin(radians(15. * (t - self.noon(J))))/cos(radians(self.altitude(t,J)))))\n",
    "\t\t\tself.zenith(t,J)\n",
    "\t\t\t#print self.Z, self.phi\n",
    "\t\t\tif t<self.noon(J):\n",
    "\t\t\t\t#return 180.- degrees(acos(-(sin(radians(self.declination(J)))-cos(radians(self.zenith(t,J)))*sin(radians(self.phi)))/(cos(radians(self.phi))*sin(radians(self.zenith(t,J))))))\n",
    "\t\t\t\treturn 180.- degrees(acos(-(sin(radians(self.declination(J)))-cos(radians(self.Z))*sin(radians(self.phi)))/(cos(radians(self.phi))*sin(radians(self.Z)))))\n",
    "\t\t\telse:\n",
    "\t\t\t\t#return 180.+ degrees(acos(-(sin(radians(self.declination(J)))-cos(radians(self.zenith(t,J)))*sin(radians(self.phi)))/(cos(radians(self.phi))*sin(radians(self.zenith(t,J))))))\n",
    "\t\t\t\treturn 180.+ degrees(acos(-(sin(radians(self.declination(J)))-cos(radians(self.Z))*sin(radians(self.phi)))/(cos(radians(self.phi))*sin(radians(self.Z)))))\n",
    "\t\t\t\n",
    "\t\tdef t_adj(self, t, J):\n",
    "\t\t\treturn t -12. + self.noon(J)\n",
    "\t\n",
    "\t\tdef declination(self, J):\t# delta\n",
    "\t\t\treturn degrees(asin(0.39785* sin(radians(278.97 + 0.9856 * J + 1.9165 * sin(radians(356.6 + 0.9856 * J))))))\n",
    "\t\n",
    "\t\tdef LC(self):\n",
    "\t\t\treturn self.longitude % 15\n",
    "\t\n",
    "\t\tdef f(self, J):\n",
    "\t\t\treturn 279.575 + 0.9856 * J\n",
    "\t\n",
    "\t\tdef ET(self, J):\n",
    "\t\t\treturn (-104.7 * sin(radians(self.f(J))) + 596.2 * sin(radians(2 * self.f(J))) + 4.3 * sin(radians(3 * self.f(J))) - 12.7 * sin(radians(4 * self.f(J))) -429.3 * cos(radians(self.f(J)))- 2.0 * cos(radians(2 * self.f(J))) + 19.3 * cos(radians(3 * self. f(J))))/3600.\n",
    "\t\n",
    "\t\tdef noon(self ,J):\n",
    "\t\t\treturn 12 + self.LC() / 15. - self.ET(J)\n",
    "\t\n",
    "\t\tdef bardd2(self, J):\n",
    "\t\t\treturn 1 + 2 * 0.01675 * cos(2 * pi / 365 * J)\n",
    "\t\n",
    "\t\tdef sol_noatm(self, t, J, a_s, alpha_s):\n",
    "\t\t\tif self.zenith(t,J)<90.:\n",
    "\t\t\t\ts = self.SOLAR * self.bardd2(J) * cos(self.a_o_i(t, J, a_s, alpha_s))##*cos(self.a_o_i(t, J, a_s, alpha_s))#(radians(self.zenith(t, J)))#cos(radians(self.zenith(t, J)))#\n",
    "\t\t\t\tif s >= 0: return s\n",
    "\t\t\t\telse: return 0. ##LN added .\n",
    "\t\t\telse: return 0. ##LN added .\n",
    "\t\n",
    "\t\tdef p_a(self):\n",
    "\t\t\treturn 101.3 * exp(-self.elev / 8200.)\n",
    "\t\n",
    "\t\tdef t_m(self, t, J):\n",
    "\t\t\t#return self.sky**(self.p_a() / (101.3 * cos(radians(self.zenith(t, J)))))\n",
    "\t\t\treturn self.sky**(self.p_a() / (101.3 * cos(radians(self.Z))))\n",
    "\t\t\n",
    "\t\tdef sol_atm(self, t, J, a_s=0, alpha_s=0):\n",
    "\t\t\treturn self.sol_noatm(t, J, a_s, alpha_s) * self.t_m(t, J)\n",
    "\t\n",
    "\t\tdef sol_diffuse(self, t, J, a_s, alpha_s):\n",
    "\t\t\treturn self.sol_noatm(t, J, a_s, alpha_s) * 0.3 * (1 - self.t_m(t, J))\n",
    "\t\n",
    "\t\tdef h_g(self, t, J, a_s=0, alpha_s=0):\n",
    "\t\t\treturn self.sol_diffuse(t, J, a_s, alpha_s) + self.sol_atm(t, J, a_s=0, alpha_s=0)\n",
    "\t\n",
    "\t\tdef short_ground(self, t, J): \n",
    "\t\t\treturn self.h_g(t, J) * self.r_soil\n",
    "\t\t\n",
    "\t\tdef gamma_t(self, t):\n",
    "\t\t\treturn 0.44 - 0.46 * sin(t * pi / 12. + 0.9) + 0.11 * sin(2. * t * pi / 12. + 0.9)\n",
    "\n",
    "\t\tdef t_air(self, t, J):\n",
    "\t\t\tif (t >= 0.) and (t <= 5.):\n",
    "\t\t\t\tself.TA = maxT(J - 1) * self.gamma_t(self.t_adj(t, J)) + minT(J) * (1 - self.gamma_t(self.t_adj(t, J)))\n",
    "\t\t\t\treturn self.TA\n",
    "\t\t\tif (t > 5.) and (t <= 14.):\n",
    "\t\t\t\tself.TA = maxT(J) * self.gamma_t(self.t_adj(t, J)) + minT(J) * (1 - self.gamma_t(self.t_adj(t, J)))\n",
    "\t\t\t\treturn self.TA\n",
    "\t\t\tif (t > 14.) and (t <= 24.):\n",
    "\t\t\t\tself.TA = maxT(J) * self.gamma_t(self.t_adj(t, J)) + minT(J + 1) * (1 - self.gamma_t(self.t_adj(t, J)))\n",
    "\t\t\t\treturn self.TA\n",
    "\n",
    "\t\tdef t_ground(self, t, J, sun = 1.):\n",
    "\t\t\tif 15 > J >= 0: # for days between jan 1 and jan 15\n",
    "\t\t\t\treturn sun * ((J - 0.)/(14. - -17.) * (janSpline(t)[0] - decSpline(t)[0]) + decSpline(t)[0]) + (1.-sun)*self.t_air(t,J)\n",
    "\t\t\tif 46 > J >= 15: # for days between jan 15 and feb 15\n",
    "\t\t\t\treturn sun *((J - 15.)/(45. - 15.) * (febSpline(t)[0] - janSpline(t)[0]) + janSpline(t)[0])+ (1.-sun)*self.t_air(t,J)\n",
    "\t\t\tif 75 > J >= 46:\n",
    "\t\t\t\treturn sun *((J - 46.)/(74. - 46.) * (marSpline(t)[0] - febSpline(t)[0]) + febSpline(t)[0]) + (1.-sun)*self.t_air(t,J)\n",
    "\t\t\tif 106 > J >= 75:\n",
    "\t\t\t\treturn sun *((J - 75.)/(105. - 75.) * (aprSpline(t) - marSpline(t)) + marSpline(t)) + (1.-sun)*self.t_air(t,J)\n",
    "\t\t\tif 136 > J >= 106:\n",
    "\t\t\t\treturn sun *((J - 106.)/(135. - 106.) * (maySpline(t) - aprSpline(t)) + aprSpline(t)) + (1.-sun)*self.t_air(t,J)\n",
    "\t\t\tif 167 > J >= 136:\n",
    "\t\t\t\treturn sun *((J - 136.)/(166. - 136.) * (junSpline(t) - maySpline(t)) + maySpline(t)) + (1.-sun)*self.t_air(t,J)\n",
    "\t\t\tif 197 > J >= 167:\n",
    "\t\t\t\treturn sun *((J - 167.)/(196. - 167.) * (julSpline(t) - junSpline(t)) + junSpline(t)) + (1.-sun)*self.t_air(t,J)\n",
    "\t\t\tif 228 > J >= 197:\n",
    "\t\t\t\treturn sun *((J - 197.)/(227. - 197.) * (augSpline(t) - julSpline(t)) + julSpline(t)) + (1.-sun)*self.t_air(t,J)\n",
    "\t\t\tif 259 > J >= 228:\n",
    "\t\t\t\treturn sun *((J - 228.)/(258. - 228.) * (sepSpline(t)[0] - augSpline(t)[0]) + augSpline(t)[0]) + (1.-sun)*self.t_air(t,J)\n",
    "\t\t\tif 289 > J >= 259:\n",
    "\t\t\t\treturn sun *((J - 259.)/(288. - 259.) * (octSpline(t)[0] - sepSpline(t)[0]) + sepSpline(t)[0]) + (1.-sun)*self.t_air(t,J)\n",
    "\t\t\tif 320 > J >= 289:\n",
    "\t\t\t\treturn sun *((J - 289.)/(319. - 289.) * (novSpline(t)[0] - octSpline(t)[0]) + octSpline(t)[0]) + (1.-sun)*self.t_air(t,J)\n",
    "\t\t\tif 350 > J >= 320:\n",
    "\t\t\t\treturn sun *((J - 320.)/(349. - 320.) * (decSpline(t)[0] - novSpline(t)[0]) + novSpline(t)[0])+ (1.-sun)*self.t_air(t,J)\n",
    "\t\t\tif 367 > J >= 350:\n",
    "\t\t\t\treturn sun *((J - 350.)/(380. - 350.) * (janSpline(t)[0] - decSpline(t)[0]) + decSpline(t)[0]) + (1.-sun)*self.t_air(t,J)\n",
    "\n",
    "\t\tdef long_atmos(self, t, J): #### UPDATE FOR MICROCLIM INPUT DATA...SUN Ta at 2m input!!\n",
    "\t\t\treturn 53.1e-14 * (self.t_air(t, J) + 273.)**6\n",
    "\t\n",
    "\t\tdef t_ave(self):\n",
    "\t\t\ttemp=[]\n",
    "\t\t\tfor i in range(24):\n",
    "\t\t\t\ttemp.append(self.t_air(i, J))\n",
    "\t\t\treturn mean(temp)\n",
    "\t\n",
    "\t\t#def t_ground(self, t, sun):\n",
    "\t\t#\treturn sun*(self.t_ave() + self.ampl * sin((pi / 12.) * (t - 8.))) + (1.- sun) * self.t_air(t)\n",
    "\t\n",
    "\t\tdef long_ground(self, t, J, sun): #### UPDATE FOR MICROCLIM INPUT DATA... sun/shade GROUND input\n",
    "\t\t\treturn self.emmiss_soil * self.SIGMA * (self.t_ground(t, J, sun) + 273)**4\n",
    "\t\n",
    "\t\tdef Ap_over_A(self, t,J,a_s, alpha_s):\n",
    "\t\t\tAp = 1.+ 4. * self.h * sin(self.a_o_i(t,J,a_s,alpha_s)) / (pi * self.d)\n",
    "\t\t\tA= 4.+ 4. * self.h / self.d\n",
    "\t\t\treturn Ap / A\n",
    "\t\n",
    "\t\tdef R_abs(self, t, J, sun, a_s, alpha_s): #### UPDATE FOR MICROCLIM INPUT DATA..... needs to call the right ground and air temp\n",
    "\t\t\treturn sun * self.abs_animal_short * (self.Ap_over_A(t,J,a_s,alpha_s) * self.sol_atm(t, J, a_s, alpha_s) + 0.5 *  self.sol_diffuse(t, J, a_s, alpha_s) + 0.5 * self.short_ground(t, J)) + 0.5 * self.abs_animal_long * (self.long_atmos(t, J) + self.long_ground(t, J, sun))\n",
    "\n",
    "\t\tdef Q_rad(self, t, J):\n",
    "\t\t\t#return self.emmiss_animal * self.SIGMA * (self.t_air(t,J) + 273.)**4\n",
    "\t\t\treturn self.emmiss_animal * self.SIGMA * (self.TA + 273.)**4\n",
    "\t\n",
    "\t\tdef g_Ha(self, wind):\n",
    "\t\t\t#return 1.4 * 0.135 * sqrt(self.wind / self.char_dim)\n",
    "\t\t\treturn 1.4 * 0.135 * sqrt(wind / self.char_dim)\n",
    "\t\n",
    "\t\tdef g_r(self, t, J, sun): #### input ground data in sun/shade here!!!\n",
    "\t\t\t#return 4.* self.SIGMA * ((self.t_air(t,J) + 273.)**3) / 29.3\n",
    "\t\t\treturn 4.* self.SIGMA * ((self.TA + 273.)**3) / 29.3\n",
    "\t\n",
    "\t\tdef a_o_i(self, t, J, a_s, alpha_s): #angle of incidence...NEED INPUTS for a_s = slope, alpha_s = aspect  with N=0 E=90\n",
    "\t\t\tself.zenith(t,J)\n",
    "\t\t\tif 0.<self.Z<=90.:\n",
    "\t\t\t\ttmp=acos(cos(radians(a_s)) * cos(radians(self.altitude(t,J))) * cos(radians(self.azimuth(t, J) - alpha_s)) + sin(radians(a_s)) * sin(radians(self.altitude(t,J))))\n",
    "\t\t\t\t#print tmp\n",
    "\t\t\t\tif tmp<=pi/2:\n",
    "\t\t\t\t\treturn tmp\n",
    "\t\t\t\tif tmp > pi/2.:\n",
    "\t\t\t\t\treturn pi/2.\n",
    "\t\t\t\telse: return pi/2\n",
    "\t\t\telse: return pi/2\n",
    "\t\t\n",
    "\t\tdef t_e(self, t, J, sun, a_s, alpha_s,veg): #UPDATE FOR MICROCLIM INPUT DATA .... t_ground should be replaced with Ta at 2cm... need to set it up for different sun/shade 2cm input\n",
    "\t\t\t#_wind = 0. #note, this is to mimic convection shadow effect\n",
    "\t\t\tif sun ==1 and veg == 0: #was or not and\n",
    "\t\t\t\tsun = 1\n",
    "\t\t\telse:\n",
    "\t\t\t\tsun = 0.375\n",
    "\t\t\tif veg == 1:\n",
    "\t\t\t\twind = 0.1\n",
    "\t\t\telse:\n",
    "\t\t\t\twind = 2.0\n",
    "\t\t\treturn self.t_ground(t, J, sun) + (self.R_abs(t, J, sun, a_s, alpha_s) - self.Q_rad(t, J)) / (29.3 * (self.g_Ha(wind) + self.g_r(t, J, sun)))\n",
    "\t\t\n",
    "\t\tdef mass(self):\n",
    "\t\t\treturn ((self.h * 100.) * pi * (0.5 * d * 100.)**2.) / 1000.\n",
    "\n",
    "\t\n",
    "\t\tdef update_tb(self):\n",
    "\t\t\tif self.tb <= self.te: \n",
    "\t\t\t\treturn self.tb + (1. - exp(-1. / self.h_t)) * (self.te - self.tb) #heating\n",
    "\t\t\tif self.tb  > self.te: \n",
    "\t\t\t\treturn self.tb + (1. - exp(-1. / self.c_t)) * (self.te - self.tb) #cooling\n",
    "\n",
    "\t\tdef mei(self, x):\n",
    "\t\t\t#mass = 10.\n",
    "\t\t\tmax_con=-328.7+22.82*x-0.32*x**2\n",
    "\t\t\tact_con= max_con * self.size /24.\n",
    "\t\t\tasin_de= 85.34-0.5*x\n",
    "\t\t\tde = (sin((3.14159/180)*asin_de))**2\n",
    "\t\t\treturn (act_con * 0.25*de/1000.*25600.)/60. #per minute\n",
    "\t\n",
    "\t\tdef smr(self):\n",
    "\t\t\tact_scope = 1. + 4. * self.moved/self.d_max #* 0.5\n",
    "\t\t\treturn (act_scope*20.2*10**(0.038*self.tb-1.771)*self.size**0.82)/60.\n",
    "\t\n",
    "\t\tdef net(self,x, dist=10.):\n",
    "\t\t\treturn self.mei(x) - self.smr(x,dist)\n",
    "\t\t\t\n",
    "\t\t\t\n",
    "\t\tdef mei(self):\n",
    "\t\t\tmax_con=-328.7+22.82*self.tb-0.32*self.tb**2\n",
    "\t\t\tact_con= max_con * self.size /24.\n",
    "\t\t\tasin_de= 85.34-0.5*self.tb\n",
    "\t\t\tde = (sin((3.14159/180)*asin_de))**2\n",
    "\t\t\treturn (act_con * 0.25*de/1000.*25600.)/60.\n",
    "\t\t\t\n",
    "\t\tdef  net(self):\n",
    "\t\t\t#E=mei(self.tb)- (dist/10.)*smr(self.tb,1)-(1.-dist/10.)*smr(self.tb,0)\n",
    "\t\t\treturn (self.moved/self.d_max)*self.smr()+(1.-self.dist/self.d_max)*self.smr()\n",
    "\t\t\t#return self.mei() - self.smr()\n",
    "\t\t\t#return E\t\n",
    "\n",
    "\t\tdef getDisp(self):\n",
    "\t\t\tincre = vonmises(self.mu, self.kappa)\n",
    "\t\t\tc = beta(self.beta_alpha, self.beta_beta) * self.d_max\n",
    "\t\t\ta, b = c * sin(self.orientation + incre), c * cos(self.orientation + incre)\n",
    "\t\t\treturn a, b, c, incre\n",
    "\t\t\n",
    "\t\tdef move_reg(self, t, J, elev, slope, aspect, veget):\n",
    "\t\t\tincre, x1, y1 = 0., self.position['x'], self.position['y'] #, 0.\n",
    "\t\t\t#self.active = 1\n",
    "\t\t\tlenx = len(elev[0])\n",
    "\t\t\tleny = len(elev)\n",
    "\t\t\taz=radians(self.azimuth(t, J))\n",
    "\t\t\tal=radians(self.altitude(t, J))\n",
    "\t\t\t#self.te = self.t_e(t,J, sun=horizon_master(az, elev[int(self.position['y'])][int(self.position['x'])], x=self.position['x'], y=self.position['y']), a_s=90.- slope[int(self.position['y'])][int(self.position['x'])], alpha_s=aspect[int(self.position['y'])][int(self.position['x'])], veg=veg) # at the moment this is just a holder from ast code\n",
    "\t\t\tdiff_tb = 50.\n",
    "\t\t\ttemp_c = 32.\n",
    "\t\t\twhile elev[int(y1)][int(x1)] == -9999 or slope[int(y1)][int(x1)] == -9999 or aspect[int(y1)][int(x1)] == -9999 or veget[int(y1)][int(x1)] == -9999:\n",
    "\t\t\t\ta, b, c, incre = self.getDisp()\n",
    "\t\t\t\tif self.thigh >= self.tb >= self.tlow:\n",
    "\t\t\t\t\tif binomial(1, self.p1): # remain still\n",
    "\t\t\t\t\t\t#self.active = 0 # x, y = x, y ... \n",
    "\t\t\t\t\t\tself.moved = 0.\n",
    "\t\t\t\t\telse: # move\n",
    "\t\t\t\t\t\tfor l in range(self.decisions): # needs to be one per second \n",
    "\t\t\t\t\t\t\ta, b, c, incre = self.getDisp()\n",
    "\t\t\t\t\t\t\tif (not (0.<=self.position['x']+a<=lenx-1)) or  (not (0.<=self.position['y']+b<=leny-1)):\n",
    "\t\t\t\t\t\t\t\twhile (not (0.<=self.position['x']+a<=lenx-1)) or  (not (0.<=self.position['y']+b<=leny-1)):\n",
    "\t\t\t\t\t\t\t\t\ta, b, c, incre = self.getDisp()\n",
    "\t\t\t\t\t\t\tx1 = self.position['x'] + a\n",
    "\t\t\t\t\t\t\ty1 = self.position['y'] + b\n",
    "\t\t\t\t\t\t\taz=radians(self.azimuth(t, J))\n",
    "\t\t\t\t\t\t\tal=radians(self.altitude(t, J))\n",
    "\t\t\t\t\t\t\tself.te = self.t_e(t,J, sun=horizon_master(az, elev[int(y1)][int(x1)], x=x1, y=y1), a_s=90.- slope[int(y1)][int(x1)], alpha_s=aspect[int(y1)][int(x1)], veg=veget[int(y1)][int(x1)])\n",
    "\t\t\t\t\t\t\tup_tb =  self.update_tb()\n",
    "\t\t\t\t\t\t\tdiff_tb_new = abs(self.t_target - up_tb)\n",
    "\t\t\t\t\t\t\tif diff_tb_new <= diff_tb: # and c <= temp_c:\n",
    "\t\t\t\t\t\t\t\tdiff_tb = diff_tb_new\n",
    "\t\t\t\t\t\t\t\ttemp_c = c\n",
    "\t\t\t\t\t\t\t\ttemp_x1, temp_y1 = x1, y1\n",
    "\t\t\t\t\t\t\telse:\n",
    "\t\t\t\t\t\t\t\tpass\t\n",
    "\t\t\t\t\t\t\tif l == (self.decisions -1):\n",
    "\t\t\t\t\t\t\t\tself.position['x'], self.position['y'], c = temp_x1, temp_y1, temp_c\n",
    "\t\t\t\t\t\t\t\tself.dist += temp_c\n",
    "\t\t\t\t\t\t\t\tself.orientation += incre\n",
    "\t\t\t\t\t\t\t\tself.moved = c\n",
    "\t\t\t\t\t\t\t\n",
    "\t\t\t\telif self.tb > self.thigh:\n",
    "\t\t\t\t\tif self.tb > self.thigh > self.t_e(t,J, sun=horizon_master(az, elev[int(y1)][int(x1)], x=x1, y=y1), a_s=90.- slope[int(y1)][int(x1)], alpha_s=aspect[int(y1)][int(x1)], veg=veget[int(y1)][int(x1)]):\n",
    "\t\t\t\t\t\tif binomial(1, self.p2):\n",
    "\t\t\t\t\t\t\tself.active = 0 #remain still\n",
    "\t\t\t\t\t\t\tself.moved = 0.\n",
    "\t\t\t\t\t\telse: # move\n",
    "\t\t\t\t\t\t\tfor l in range(self.decisions): # needs to be one per second \n",
    "\t\t\t\t\t\t\t\ta, b, c, incre = self.getDisp()\n",
    "\t\t\t\t\t\t\t\tif (not (0.<=self.position['x']+a<=lenx-1)) or  (not (0.<=self.position['y']+b<=leny-1)):\n",
    "\t\t\t\t\t\t\t\t\twhile (not (0.<=self.position['x']+a<=lenx-1)) or  (not (0.<=self.position['y']+b<=leny-1)):\n",
    "\t\t\t\t\t\t\t\t\t\ta, b, c, incre = self.getDisp()\n",
    "\t\t\t\t\t\t\t\tx1 = self.position['x'] + a\n",
    "\t\t\t\t\t\t\t\ty1 = self.position['y'] + b\n",
    "\t\t\t\t\t\t\t\taz=radians(self.azimuth(t, J))\n",
    "\t\t\t\t\t\t\t\tal=radians(self.altitude(t, J))\n",
    "\t\t\t\t\t\t\t\tself.te = self.t_e(t,J, sun=horizon_master(az, elev[int(y1)][int(x1)], x=x1, y=y1), a_s=90.- slope[int(y1)][int(x1)], alpha_s=aspect[int(y1)][int(x1)], veg=veget[int(y1)][int(x1)])\n",
    "\t\t\t\t\t\t\t\tup_tb =  self.update_tb()\n",
    "\t\t\t\t\t\t\t\tdiff_tb_new = abs(self.t_target - up_tb)\n",
    "\t\t\t\t\t\t\t\tif diff_tb_new <= diff_tb: # and c <= temp_c:\n",
    "\t\t\t\t\t\t\t\t\tdiff_tb = diff_tb_new\n",
    "\t\t\t\t\t\t\t\t\ttemp_c = c\n",
    "\t\t\t\t\t\t\t\t\ttemp_x1, temp_y1 = x1, y1\n",
    "\t\t\t\t\t\t\t\telse:\n",
    "\t\t\t\t\t\t\t\t\tpass\t\t\t\n",
    "\t\t\t\t\t\t\t\tif l == (self.decisions -1):\n",
    "\t\t\t\t\t\t\t\t\tself.position['x'], self.position['y'], c = temp_x1, temp_y1, temp_c\n",
    "\t\t\t\t\t\t\t\t\tself.dist += temp_c\n",
    "\t\t\t\t\t\t\t\t\tself.orientation += incre\n",
    "\t\t\t\t\t\t\t\t\tself.moved = c\n",
    "\t\t\t\t\telse: # move\n",
    "\t\t\t\t\t\tfor l in range(self.decisions): # needs to be one per second \n",
    "\t\t\t\t\t\t\ta, b, c, incre = self.getDisp()\n",
    "\t\t\t\t\t\t\tif (not (0.<=self.position['x']+a<=lenx-1)) or  (not (0.<=self.position['y']+b<=leny-1)):\n",
    "\t\t\t\t\t\t\t\twhile (not (0.<=self.position['x']+a<=lenx-1)) or  (not (0.<=self.position['y']+b<=leny-1)):\n",
    "\t\t\t\t\t\t\t\t\ta, b, c, incre = self.getDisp()\n",
    "\t\t\t\t\t\t\tx1 = self.position['x'] + a\n",
    "\t\t\t\t\t\t\ty1 = self.position['y'] + b\n",
    "\t\t\t\t\t\t\taz=radians(self.azimuth(t, J))\n",
    "\t\t\t\t\t\t\tal=radians(self.altitude(t, J))\n",
    "\t\t\t\t\t\t\tself.te = self.t_e(t,J, sun=horizon_master(az, elev[int(y1)][int(x1)], x=x1, y=y1), a_s=90.- slope[int(y1)][int(x1)], alpha_s=aspect[int(y1)][int(x1)], veg=veget[int(y1)][int(x1)])\n",
    "\t\t\t\t\t\t\tup_tb =  self.update_tb()\n",
    "\t\t\t\t\t\t\tdiff_tb_new = abs(self.t_target - up_tb)\n",
    "\t\t\t\t\t\t\tif diff_tb_new <= diff_tb:# and c <= temp_c:\n",
    "\t\t\t\t\t\t\t\tdiff_tb = diff_tb_new\n",
    "\t\t\t\t\t\t\t\ttemp_c = c\n",
    "\t\t\t\t\t\t\t\ttemp_x1, temp_y1 = x1, y1\n",
    "\t\t\t\t\t\t\telse:\n",
    "\t\t\t\t\t\t\t\tpass\t\n",
    "\t\t\t\t\t\t\tif l == (self.decisions -1):\n",
    "\t\t\t\t\t\t\t\tself.position['x'], self.position['y'], c = temp_x1, temp_y1, temp_c\n",
    "\t\t\t\t\t\t\t\tself.dist += temp_c\n",
    "\t\t\t\t\t\t\t\tself.orientation += incre\n",
    "\t\t\t\t\t\t\t\tself.moved = c\n",
    "\t\n",
    "\t\t\t\telse: #self.tb < self.tlow:\n",
    "\t\t\t\t\tif self.tb < self.tlow < self.t_e(t,J, sun=horizon_master(az, elev[int(y1)][int(x1)], x=x1, y=y1), a_s=90.- slope[y1][x1], alpha_s=aspect[int(y1)][int(x1)], veg=veget[int(y1)][int(x1)]):\n",
    "\t\t\t\t\t\tif binomial(1, self.p3):\n",
    "\t\t\t\t\t\t\tself.active = 0 #remain still\n",
    "\t\t\t\t\t\t\tself.moved = 0.\n",
    "\t\t\t\t\t\telse: # move\n",
    "\t\t\t\t\t\t\tfor l in range(self.decisions): # needs to be one per second \n",
    "\t\t\t\t\t\t\t\ta, b, c, incre = self.getDisp()\n",
    "\t\t\t\t\t\t\t\tif (not (0.<=self.position['x']+a<=lenx-1)) or  (not (0.<=self.position['y']+b<=leny-1)):\n",
    "\t\t\t\t\t\t\t\t\twhile (not (0.<=self.position['x']+a<=lenx-1)) or  (not (0.<=self.position['y']+b<=leny-1)):\n",
    "\t\t\t\t\t\t\t\t\t\ta, b, c, incre = self.getDisp()\n",
    "\t\t\t\t\t\t\t\tx1 = self.position['x'] + a\n",
    "\t\t\t\t\t\t\t\ty1 = self.position['y'] + b\n",
    "\t\t\t\t\t\t\t\taz=radians(self.azimuth(t, J))\n",
    "\t\t\t\t\t\t\t\tal=radians(self.altitude(t, J))\n",
    "\t\t\t\t\t\t\t\tself.te = self.t_e(t,J, sun=horizon_master(az, elev[int(y1)][int(x1)], x=x1, y=y1), a_s=90.- slope[int(y1)][int(x1)], alpha_s=aspect[int(y1)][int(x1)], veg=veget[int(y1)][int(x1)])\n",
    "\t\t\t\t\t\t\t\tup_tb =  self.update_tb()\t\n",
    "\t\t\t\t\t\t\t\tdiff_tb_new = abs(self.t_target - up_tb)\n",
    "\t\t\t\t\t\t\t\tif diff_tb_new <= diff_tb: # and c <= temp_c:\n",
    "\t\t\t\t\t\t\t\t\tdiff_tb = diff_tb_new\n",
    "\t\t\t\t\t\t\t\t\ttemp_c = c\n",
    "\t\t\t\t\t\t\t\t\ttemp_x1, temp_y1 = x1, y1\n",
    "\t\t\t\t\t\t\t\telse:\n",
    "\t\t\t\t\t\t\t\t\tpass\t\n",
    "\t\t\t\t\t\t\t\tif l == (self.decisions - 1):\n",
    "\t\t\t\t\t\t\t\t\tself.position['x'], self.position['y'], c = temp_x1, temp_y1, temp_c\n",
    "\t\t\t\t\t\t\t\t\tself.dist += temp_c\n",
    "\t\t\t\t\t\t\t\t\tself.orientation += incre\n",
    "\t\t\t\t\t\t\t\t\tself.moved = c\n",
    "\t\t\t\t\telse:\n",
    "\t\t\t\t\t\tfor l in range(self.decisions): # needs to be one per second \n",
    "\t\t\t\t\t\t\ta, b, c, incre = self.getDisp()\n",
    "\t\t\t\t\t\t\tif (not (0.<=self.position['x']+a<=lenx-1)) or  (not (0.<=self.position['y']+b<=leny-1)):\n",
    "\t\t\t\t\t\t\t\twhile (not (0.<=self.position['x']+a<=lenx-1)) or  (not (0.<=self.position['y']+b<=leny-1)):\n",
    "\t\t\t\t\t\t\t\t\ta, b, c, incre = self.getDisp()\n",
    "\t\t\t\t\t\t\tx1 = self.position['x'] + a\n",
    "\t\t\t\t\t\t\ty1 = self.position['y'] + b\n",
    "\t\t\t\t\t\t\taz=radians(self.azimuth(t, J))\n",
    "\t\t\t\t\t\t\tal=radians(self.altitude(t, J))\n",
    "\t\t\t\t\t\t\tself.te = self.t_e(t,J, sun=horizon_master(az, elev[int(y1)][int(x1)], x=x1, y=y1), a_s=90.- slope[int(y1)][int(x1)], alpha_s=aspect[int(y1)][int(x1)], veg=veget[int(y1)][int(x1)])\n",
    "\t\t\t\t\t\t\tup_tb =  self.update_tb()\n",
    "\t\t\t\t\t\t\tdiff_tb_new = abs(self.t_target - up_tb)\n",
    "\t\t\t\t\t\t\tif diff_tb_new <= diff_tb: # and c <= temp_c:\n",
    "\t\t\t\t\t\t\t\tdiff_tb = diff_tb_new\n",
    "\t\t\t\t\t\t\t\ttemp_c = c\n",
    "\t\t\t\t\t\t\t\ttemp_x1, temp_y1 = x1, y1\n",
    "\t\t\t\t\t\t\telse:\n",
    "\t\t\t\t\t\t\t\tpass\t\t\t\n",
    "\t\t\t\t\t\t\tif l == (self.decisions - 1):\n",
    "\t\t\t\t\t\t\t\tself.position['x'], self.position['y'], c = temp_x1, temp_y1, temp_c\n",
    "\t\t\t\t\t\t\t\tself.dist += temp_c\n",
    "\t\t\t\t\t\t\t\tself.orientation += incre\n",
    "\t\t\t\t\t\t\t\tself.moved = c\t\n",
    "\n",
    "\n",
    "    ### MAKE MAPS\n",
    "\t#elev, slope, aspect, veget, e_range, e_fract, v_range, v_fract = make_maps(h,  high,  N,  MU,  SD)  ###only need these first 4 variables -- assign \n",
    "\n",
    "\n",
    "\t\n",
    "\t# Save Maps\n",
    "\t#write_file = \"maps/elev\" +  \"_\" + '{:.4}'.format(e_fract) + \"_\" + '{:.3}'.format(e_range)+ \"_\" + '{:.3}'.format(v_fract+3.) + \"_\" + '{:.3}'.format(v_range)+ \".asc\"\n",
    "\t#output=open(write_file, 'w')\n",
    "\t#for i in range(N):\n",
    "\t#\tfor j in range(N):\n",
    "\t#\t\toutput.write(str(elev[i][j])), output.write(' ')\n",
    "\t#\toutput.write(\"\\n\")\n",
    "\t#output.close()\n",
    "\t#\n",
    "\t#write_file = \"maps/slope\" +  \"_\" +'{:.4}'.format(e_fract) + \"_\" + '{:.3}'.format(e_range)+ \"_\" + '{:.3}'.format(v_fract+3.) + \"_\" + '{:.3}'.format(v_range)+ \".asc\"\n",
    "\t#output=open(write_file, 'w')\n",
    "\t#for i in range(N):\n",
    "\t#\tfor j in range(N):\n",
    "\t#\t\toutput.write(str(slope[i][j])), output.write(' ')\n",
    "\t#\toutput.write(\"\\n\")\n",
    "\t#output.close()\n",
    "\t#\n",
    "\t#write_file = \"maps/aspect\" + \"_\" + '{:.4}'.format(e_fract) + \"_\" + '{:.3}'.format(e_range)+ \"_\" + '{:.3}'.format(v_fract+3.) + \"_\" + '{:.3}'.format(v_range)+ \".asc\"\n",
    "\t#output=open(write_file, 'w')\n",
    "\t#for i in range(N):\n",
    "\t#\tfor j in range(N):\n",
    "\t#\t\toutput.write(str(aspect[i][j])), output.write(' ')\n",
    "\t#\toutput.write(\"\\n\")\n",
    "\t#output.close()\n",
    "\t#\n",
    "\t#write_file = \"maps/veg\" + \"_\" + '{:.4}'.format(e_fract) + \"_\" + '{:.3}'.format(e_range)+ \"_\" + '{:.3}'.format(v_fract+3.) + \"_\" + '{:.3}'.format(v_range)+ \".asc\"\n",
    "\t#output=open(write_file, 'w')\n",
    "\t#for i in range(N):\n",
    "\t#\tfor j in range(N):\n",
    "\t#\t\tif veget[i][j]:\n",
    "\t#\t\t\toutput.write(str(1)), output.write(' ')\n",
    "\t#\t\telse:\n",
    "\t#\t\t\toutput.write(str(0)), output.write(' ')\t\t\t\t\n",
    "\t#\toutput.write(\"\\n\")\n",
    "\t#output.close()\n",
    "\n",
    "\n",
    "\tlimit = N\n",
    "\tdata = np.array(elev[:])\n",
    "\tdef horizon_master(azimuth, elevation, x, y):\n",
    "\t\tif find_horizon(azimuth,x,y) > elevation:\n",
    "\t\t\t#print True\n",
    "\t\t\treturn 0 #True is shaded, so we want '0' for full sun\n",
    "\t\telse:\n",
    "\t\t\t#print False\n",
    "\t\t\treturn 1\n",
    "\n",
    "\tcellsize = 1.\t\n",
    "\tdef find_horizon(azimuth, x, y):\n",
    "\t\tx, y=math.floor(x), math.floor(y)\n",
    "\t\tz_data = [data[y][x]]\n",
    "\t\thyp = 0.25# float(elev_header[4][1])/2. # 0.25#0.25# 0.25 \n",
    "\t\twhile 0 < x+ hyp*sin(azimuth) < len(data[0]) and 0 < y + hyp*cos(azimuth) < len(data):\n",
    "\t\t\t#print x+ hyp*sin(azimuth), y + hyp*cos(azimuth)\n",
    "\t\t\tx = x + hyp*sin(azimuth)*cellsize #added int # check whether this should be mulitplied by cell size\n",
    "\t\t\ty = y + hyp*cos(azimuth)*cellsize\n",
    "\t\t\tif data[math.floor(y)][math.floor(x)] not in z_data:\n",
    "\t\t\t\tz_data.append(data[math.floor(y)][math.floor(x)])\n",
    "\t\th_max = -pi/2.\n",
    "\t\tfor i in range(1, len(z_data)):\n",
    "\t\t\th_angle = atan((z_data[i]-z_data[0])/(i*cellsize))\n",
    "\t\t\tif h_angle > h_max:\n",
    "\t\t\t\th_max = h_angle\n",
    "\t\treturn h_max\n",
    "\t\n",
    "\t\t\n",
    "\t# Run activity sims\n",
    "\tshadows=np.copy(elev[:])\n",
    "\ttemperatures = np.copy(elev[:])\n",
    "\tactivity_temps = np.copy(elev[:])\n",
    "\tcritical_temps = np.copy(elev[:])\n",
    "\td_temps = np.copy(elev[:])\n",
    "\tperformance = np.copy(elev[:])\n",
    "\tact_sum=[]\n",
    "\tcrit_sum=[]\n",
    "\ttemp_temp = []\n",
    "\td_temp = []\n",
    "\tact_temp =[]\n",
    "\tcrit_temp = []\n",
    "\tp_temp = []\n",
    "\ttemp_percents, d_percents, p_percents=[],[],[]\n",
    "\n",
    "\ttimes = np.arange(6.,20.,1./60.).tolist()\n",
    "\tpop_size = 20\n",
    "\tpopulation = [T_e() for i in range(pop_size)] #population is a list, where each element is an instance of the \"T_e\" class that's assigned to variable \"lizard\" each iteration\n",
    "    #ppulation is defined and instantiated with a list of T_e class objects before sim_act function is called\n",
    "    # sim_act\n",
    "\n",
    "\toutput = open(name_of_sim+\".csv\", 'a')\n",
    "\toutput2 = open(name_of_sim2+\".csv\", 'a')\n",
    "\n",
    "\n",
    "\t\n",
    "\tdef sim_act(lizard, t,J=150.):\n",
    "\t\tlizard.zenith(t,J)\n",
    "\t\tif lizard.Z > 0.:\n",
    "\t\t\tlizard.te = lizard.t_e(t,J, sun=horizon_master(radians(lizard.azimuth(t, J)), elev[math.floor(lizard.position['y'])][math.floor(lizard.position['x'])], x=lizard.position['x'], y=lizard.position['y']), a_s=90.- slope[math.floor(lizard.position['y'])][math.floor(lizard.position['x'])], alpha_s=aspect[math.floor(lizard.position['y'])][math.floor(lizard.position['x'])], veg=veget[math.floor(lizard.position['y'])][math.floor(lizard.position['x'])])\n",
    "\t\t\t#print lizard.zenith(t,J)\n",
    "\t\t\tif lizard.tb < lizard.te < lizard.tlow:\n",
    "\t\t\t\tlizard.moved = 0.\n",
    "\t\t\t\tlizard.active = 0.\n",
    "\t\t\t\tlizard.tb = lizard.update_tb()\n",
    "\t\t\t\tlizard.energy_balance += lizard.smr()\n",
    "\t\t\t\t#lizard.tot_dist_moved += lizard.moved\n",
    "\t\t\n",
    "\t\t\telif lizard.tlow <= lizard.tb <=lizard.ctmax:\n",
    "\t\t\t\tlizard.move_reg(t,J,elev, slope, aspect, veget)\n",
    "\t\t\t\tlizard.tb = lizard.update_tb()\n",
    "\t\t\t\tlizard.energy_balance += lizard.smr()\n",
    "\t\t\t\tlizard.tot_dist_moved += lizard.moved\n",
    "\t\t\t\tlizard.ate += binomial(1, 0.05)\n",
    "\t\t\t\tlizard.total_activity += 1.\n",
    "\t\t\t\tlizard.active = 1.\n",
    "\t\t\telse:\n",
    "\t\t\t\tif lizard.te > lizard.ctmax:\n",
    "\t\t\t\t\tlizard.moved = 0.\n",
    "\t\t\t\t\tlizard.active = 0.\n",
    "\t\t\t\t\tlizard.energy_balance += lizard.smr()\n",
    "\t\t\t\t\t#lizard.tot_dist_moved += lizard.moved\n",
    "\t\t\t\telse: \n",
    "\t\t\t\t\tlizard.active = 0.\n",
    "\t\t\t\t\tlizard.tb = lizard.update_tb()\n",
    "\t\t\t\t\tlizard.energy_balance += lizard.smr()\n",
    "\t\t\t\t\t#lizard.tot_dist_moved += lizard.moved\n",
    "\t\t\t\t\t#lizard.ate += binomial(1, 0.05)\n",
    "\t\t\t\t\t#lizard.total_activity += 1.\n",
    "\t\t\toutput.write(str(t)+ \"\\t\" + str(lizard.position['x']) + \"\\t\" +str(lizard.position['y']) + \"\\t\" + str(lizard.te) + \"\\t\" +str(lizard.tb) + \"\\t\" + str(lizard.active) +\"\\t\" + str(lizard.mei())+\"\\t\" + str(lizard.smr()) + \"\\t\" + str(lizard.net())  + \"\\t\"+ str(lizard.moved)+ \"\\t\"+ str(lizard.ate)+\"\\n\")\n",
    "            #output.write(str(e_range) + \"\\t\" + str(e_fract) + \"\\t\" + str(v_range) + \"\\t\" + str(v_fract) +\"\\t\" + str(t)+ \"\\t\" + str(lizard.position['x']) + \"\\t\" +str(lizard.position['y']) + \"\\t\" + str(lizard.te) + \"\\t\" +str(lizard.tb) + \"\\t\" + str(lizard.active) +\"\\t\" + str(lizard.mei())+\"\\t\" + str(lizard.smr()) + \"\\t\" + str(lizard.net())  + \"\\t\"+ str(lizard.moved)+ \"\\t\"+ str(lizard.ate)+\"\\n\")\n",
    "\tfor t in times:\n",
    "\t\tprint(J, t)\n",
    "\t\t[sim_act(lizard,t,J) for lizard in population]\n",
    "\n",
    "\ttot_act, tot_smr, tot_move, tot_ave_ate =0.,0.,0., 0.\n",
    "\tfor lizard in population:\n",
    "\t\ttot_act += float(lizard.total_activity)\n",
    "\t\ttot_smr += lizard.energy_balance\n",
    "\t\ttot_move += lizard.tot_dist_moved\n",
    "\t\ttot_ave_ate += float(lizard.ate)\n",
    "\toutput2.write(str(t)+ \"\\t\" + str(lizard.position['x']) + \"\\t\" +str(lizard.position['y']) + \"\\t\" + str(lizard.te) + \"\\t\" +str(lizard.tb) + \"\\t\" + str(lizard.active) +\"\\t\" + str(lizard.mei())+\"\\t\" + str(lizard.smr()) + \"\\t\" + str(lizard.net())  + \"\\t\"+ str(lizard.moved)+ \"\\t\"+ str(lizard.ate)+\"\\n\")\n",
    "\t\t\t\n",
    "st1= time()\n",
    "name_of_sim=\"details\"#raw_input(\"name of file? \")\n",
    "output = open(name_of_sim+\".csv\", 'w')\n",
    "name_of_sim2=\"summary\"#raw_input(\"name of file? \")\n",
    "output2 = open(name_of_sim2+\".csv\", 'w')\n",
    "#output2.write(\"elev\\tefract\\tveg\\tvfract\\tactive\\tsmr\\tmoved\\tate\\n\") #MS original\n",
    "output2.write(\"elev\\tactive\\tsmr\\tmoved\\tate\\n\")\n",
    "output2.close()\n",
    "#output.write(\"elev\\tefract\\tveg\\tvfract\\ttime\\tX\\tY\\tte\\ttb\\tactive\\tmei\\tsmr\\tnet\\tmoved\\tate\\n\") #MS original\n",
    "output.write(\"elev\\ttime\\tX\\tY\\tte\\ttb\\tactive\\tmei\\tsmr\\tnet\\tmoved\\tate\\n\")\n",
    "output.close()\n",
    "\n",
    "\t#for t in times:\n",
    "    #    print(J, t)\n",
    "    #    for lizard in population:\n",
    "    #        print(\"lizard x position:\", lizard.position['x'])\n",
    "    #        print(\"lizard y position:\", lizard.position['y'])\n",
    "    #        print(\"elev array size:\", len(elev))\n",
    "    #        print(\"slope array size:\", len(slope))\n",
    "    #        print(\"aspect array size:\", len(aspect))\n",
    "    #        print(\"veget array size:\", len(veget))\n",
    "    #        sim_act(lizard,t,J)\n",
    "    \n",
    "\n",
    "\"\"\"\n",
    "jobs = []\n",
    "for x in range(1000):\n",
    "\tp=multiprocessing.Process(target=do_sim, args=())\n",
    "\tjobs.append(p)\n",
    "\tp.start()\n",
    "for x in range(10):\n",
    "\tjobs[x].join()\n",
    "\"\"\"\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "from multiprocessing import Pool\n",
    "pool = Pool(8)\n",
    "for _ in range(256): #was 1000 == number of maps , or the simulation number\n",
    "\tpool.apply_async(do_sim, ())\n",
    "pool.close()\n",
    "\n",
    "pool.join()\n",
    "\n",
    "st2 = time()\n",
    "print(st2-st1)\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "\n",
    "do_sim() #to run it one \n",
    "output.close()\n",
    "output2.close()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\t\t### original fxn \n",
    "\t\t\n",
    "\t\t#def move_reg(self, t, J, elev, slope, aspect, veget):\n",
    "\t\t#\tincre, x1, y1 = 0., self.position['x'], self.position['y'] #, 0.\n",
    "\t\t#\t#self.active = 1\n",
    "\t\t#\tlenx = len(elev[0])\n",
    "\t\t#\tleny = len(elev)\n",
    "\t\t#\taz=radians(self.azimuth(t, J))\n",
    "\t\t#\tal=radians(self.altitude(t, J))\n",
    "\t\t#\t#self.te = self.t_e(t,J, sun=horizon_master(az, elev[int(self.position['y'])][int(self.position['x'])], x=self.position['x'], y=self.position['y']), a_s=90.- slope[int(self.position['y'])][int(self.position['x'])], alpha_s=aspect[int(self.position['y'])][int(self.position['x'])], veg=veg) # at the moment this is just a holder from ast code\n",
    "\t\t#\tdiff_tb = 50.\n",
    "\t\t#\ttemp_c = 32.\n",
    "\t\t#\tif self.thigh >= self.tb >= self.tlow:\n",
    "\t\t#\t\tif binomial(1, self.p1): # remain still\n",
    "\t\t#\t\t\t#self.active = 0 # x, y = x, y ... \n",
    "\t\t#\t\t\tself.moved = 0.\n",
    "\t\t#\t\telse: # move\n",
    "\t\t#\t\t\tfor l in range(self.decisions): # needs to be one per second \n",
    "\t\t#\t\t\t\ta, b, c, incre = self.getDisp()\n",
    "\t\t#\t\t\t\tif (not (0.<=self.position['x']+a<=lenx-1)) or  (not (0.<=self.position['y']+b<=leny-1)):\n",
    "\t\t#\t\t\t\t\twhile (not (0.<=self.position['x']+a<=lenx-1)) or  (not (0.<=self.position['y']+b<=leny-1)):\n",
    "\t\t#\t\t\t\t\t\ta, b, c, incre = self.getDisp()\n",
    "\t\t#\t\t\t\tx1 = self.position['x'] + a\n",
    "\t\t#\t\t\t\ty1 = self.position['y'] + b\n",
    "\t\t#\t\t\t\taz=radians(self.azimuth(t, J))\n",
    "\t\t#\t\t\t\tal=radians(self.altitude(t, J))\n",
    "\t\t#\t\t\t\tself.te = self.t_e(t,J, sun=horizon_master(az, elev[int(y1)][int(x1)], x=x1, y=y1), a_s=90.- slope[int(y1)][int(x1)], alpha_s=aspect[int(y1)][int(x1)], veg=veget[int(y1)][int(x1)])\n",
    "\t\t#\t\t\t\tup_tb =  self.update_tb()\n",
    "\t\t#\t\t\t\tdiff_tb_new = abs(self.t_target - up_tb)\n",
    "\t\t#\t\t\t\tif diff_tb_new <= diff_tb: # and c <= temp_c:\n",
    "\t\t#\t\t\t\t\tdiff_tb = diff_tb_new\n",
    "\t\t#\t\t\t\t\ttemp_c = c\n",
    "\t\t#\t\t\t\t\ttemp_x1, temp_y1 = x1, y1\n",
    "\t\t#\t\t\t\telse:\n",
    "\t\t#\t\t\t\t\tpass\t\n",
    "\t\t#\t\t\t\tif l == (self.decisions -1):\n",
    "\t\t#\t\t\t\t\tself.position['x'], self.position['y'], c = temp_x1, temp_y1, temp_c\n",
    "\t\t#\t\t\t\t\tself.dist += temp_c\n",
    "\t\t#\t\t\t\t\tself.orientation += incre\n",
    "\t\t#\t\t\t\t\tself.moved = c\n",
    "\t\t#\t\t\t\t\n",
    "\t\t#\telif self.tb > self.thigh:\n",
    "\t\t#\t\tif self.tb > self.thigh > self.t_e(t,J, sun=horizon_master(az, elev[int(y1)][int(x1)], x=x1, y=y1), a_s=90.- slope[int(y1)][int(x1)], alpha_s=aspect[int(y1)][int(x1)], veg=veget[int(y1)][int(x1)]):\n",
    "\t\t#\t\t\tif binomial(1, self.p2):\n",
    "\t\t#\t\t\t\tself.active = 0 #remain still\n",
    "\t\t#\t\t\t\tself.moved = 0.\n",
    "\t\t#\t\t\telse: # move\n",
    "\t\t#\t\t\t\tfor l in range(self.decisions): # needs to be one per second \n",
    "\t\t#\t\t\t\t\ta, b, c, incre = self.getDisp()\n",
    "\t\t#\t\t\t\t\tif (not (0.<=self.position['x']+a<=lenx-1)) or  (not (0.<=self.position['y']+b<=leny-1)):\n",
    "\t\t#\t\t\t\t\t\twhile (not (0.<=self.position['x']+a<=lenx-1)) or  (not (0.<=self.position['y']+b<=leny-1)):\n",
    "\t\t#\t\t\t\t\t\t\ta, b, c, incre = self.getDisp()\n",
    "\t\t#\t\t\t\t\tx1 = self.position['x'] + a\n",
    "\t\t#\t\t\t\t\ty1 = self.position['y'] + b\n",
    "\t\t#\t\t\t\t\taz=radians(self.azimuth(t, J))\n",
    "\t\t#\t\t\t\t\tal=radians(self.altitude(t, J))\n",
    "\t\t#\t\t\t\t\tself.te = self.t_e(t,J, sun=horizon_master(az, elev[int(y1)][int(x1)], x=x1, y=y1), a_s=90.- slope[int(y1)][int(x1)], alpha_s=aspect[int(y1)][int(x1)], veg=veget[int(y1)][int(x1)])\n",
    "\t\t#\t\t\t\t\tup_tb =  self.update_tb()\n",
    "\t\t#\t\t\t\t\tdiff_tb_new = abs(self.t_target - up_tb)\n",
    "\t\t#\t\t\t\t\tif diff_tb_new <= diff_tb: # and c <= temp_c:\n",
    "\t\t#\t\t\t\t\t\tdiff_tb = diff_tb_new\n",
    "\t\t#\t\t\t\t\t\ttemp_c = c\n",
    "\t\t#\t\t\t\t\t\ttemp_x1, temp_y1 = x1, y1\n",
    "\t\t#\t\t\t\t\telse:\n",
    "\t\t#\t\t\t\t\t\tpass\t\t\t\n",
    "\t\t#\t\t\t\t\tif l == (self.decisions -1):\n",
    "\t\t#\t\t\t\t\t\tself.position['x'], self.position['y'], c = temp_x1, temp_y1, temp_c\n",
    "\t\t#\t\t\t\t\t\tself.dist += temp_c\n",
    "\t\t#\t\t\t\t\t\tself.orientation += incre\n",
    "\t\t#\t\t\t\t\t\tself.moved = c\n",
    "\t\t#\t\telse: # move\n",
    "\t\t#\t\t\tfor l in range(self.decisions): # needs to be one per second \n",
    "\t\t#\t\t\t\ta, b, c, incre = self.getDisp()\n",
    "\t\t#\t\t\t\tif (not (0.<=self.position['x']+a<=lenx-1)) or  (not (0.<=self.position['y']+b<=leny-1)):\n",
    "\t\t#\t\t\t\t\twhile (not (0.<=self.position['x']+a<=lenx-1)) or  (not (0.<=self.position['y']+b<=leny-1)):\n",
    "\t\t#\t\t\t\t\t\ta, b, c, incre = self.getDisp()\n",
    "\t\t#\t\t\t\tx1 = self.position['x'] + a\n",
    "\t\t#\t\t\t\ty1 = self.position['y'] + b\n",
    "\t\t#\t\t\t\taz=radians(self.azimuth(t, J))\n",
    "\t\t#\t\t\t\tal=radians(self.altitude(t, J))\n",
    "\t\t#\t\t\t\tself.te = self.t_e(t,J, sun=horizon_master(az, elev[int(y1)][int(x1)], x=x1, y=y1), a_s=90.- slope[int(y1)][int(x1)], alpha_s=aspect[int(y1)][int(x1)], veg=veget[int(y1)][int(x1)])\n",
    "\t\t#\t\t\t\tup_tb =  self.update_tb()\n",
    "\t\t#\t\t\t\tdiff_tb_new = abs(self.t_target - up_tb)\n",
    "\t\t#\t\t\t\tif diff_tb_new <= diff_tb:# and c <= temp_c:\n",
    "\t\t#\t\t\t\t\tdiff_tb = diff_tb_new\n",
    "\t\t#\t\t\t\t\ttemp_c = c\n",
    "\t\t#\t\t\t\t\ttemp_x1, temp_y1 = x1, y1\n",
    "\t\t#\t\t\t\telse:\n",
    "\t\t#\t\t\t\t\tpass\t\n",
    "\t\t#\t\t\t\tif l == (self.decisions -1):\n",
    "\t\t#\t\t\t\t\tself.position['x'], self.position['y'], c = temp_x1, temp_y1, temp_c\n",
    "\t\t#\t\t\t\t\tself.dist += temp_c\n",
    "\t\t#\t\t\t\t\tself.orientation += incre\n",
    "\t\t#\t\t\t\t\tself.moved = c\n",
    "#\n",
    "\t\t#\telse: #self.tb < self.tlow:\n",
    "\t\t#\t\tif self.tb < self.tlow < self.t_e(t,J, sun=horizon_master(az, elev[int(y1)][int(x1)], x=x1, y=y1), a_s=90.- slope[y1][x1], alpha_s=aspect[int(y1)][int(x1)], veg=veget[int(y1)][int(x1)]):\n",
    "\t\t#\t\t\tif binomial(1, self.p3):\n",
    "\t\t#\t\t\t\tself.active = 0 #remain still\n",
    "\t\t#\t\t\t\tself.moved = 0.\n",
    "\t\t#\t\t\telse: # move\n",
    "\t\t#\t\t\t\tfor l in range(self.decisions): # needs to be one per second \n",
    "\t\t#\t\t\t\t\ta, b, c, incre = self.getDisp()\n",
    "\t\t#\t\t\t\t\tif (not (0.<=self.position['x']+a<=lenx-1)) or  (not (0.<=self.position['y']+b<=leny-1)):\n",
    "\t\t#\t\t\t\t\t\twhile (not (0.<=self.position['x']+a<=lenx-1)) or  (not (0.<=self.position['y']+b<=leny-1)):\n",
    "\t\t#\t\t\t\t\t\t\ta, b, c, incre = self.getDisp()\n",
    "\t\t#\t\t\t\t\tx1 = self.position['x'] + a\n",
    "\t\t#\t\t\t\t\ty1 = self.position['y'] + b\n",
    "\t\t#\t\t\t\t\taz=radians(self.azimuth(t, J))\n",
    "\t\t#\t\t\t\t\tal=radians(self.altitude(t, J))\n",
    "\t\t#\t\t\t\t\tself.te = self.t_e(t,J, sun=horizon_master(az, elev[int(y1)][int(x1)], x=x1, y=y1), a_s=90.- slope[int(y1)][int(x1)], alpha_s=aspect[int(y1)][int(x1)], veg=veget[int(y1)][int(x1)])\n",
    "\t\t#\t\t\t\t\tup_tb =  self.update_tb()\t\n",
    "\t\t#\t\t\t\t\tdiff_tb_new = abs(self.t_target - up_tb)\n",
    "\t\t#\t\t\t\t\tif diff_tb_new <= diff_tb: # and c <= temp_c:\n",
    "\t\t#\t\t\t\t\t\tdiff_tb = diff_tb_new\n",
    "\t\t#\t\t\t\t\t\ttemp_c = c\n",
    "\t\t#\t\t\t\t\t\ttemp_x1, temp_y1 = x1, y1\n",
    "\t\t#\t\t\t\t\telse:\n",
    "\t\t#\t\t\t\t\t\tpass\t\n",
    "\t\t#\t\t\t\t\tif l == (self.decisions - 1):\n",
    "\t\t#\t\t\t\t\t\tself.position['x'], self.position['y'], c = temp_x1, temp_y1, temp_c\n",
    "\t\t#\t\t\t\t\t\tself.dist += temp_c\n",
    "\t\t#\t\t\t\t\t\tself.orientation += incre\n",
    "\t\t#\t\t\t\t\t\tself.moved = c\n",
    "\t\t#\t\telse:\n",
    "\t\t#\t\t\tfor l in range(self.decisions): # needs to be one per second \n",
    "\t\t#\t\t\t\ta, b, c, incre = self.getDisp()\n",
    "\t\t#\t\t\t\tif (not (0.<=self.position['x']+a<=lenx-1)) or  (not (0.<=self.position['y']+b<=leny-1)):\n",
    "\t\t#\t\t\t\t\twhile (not (0.<=self.position['x']+a<=lenx-1)) or  (not (0.<=self.position['y']+b<=leny-1)):\n",
    "\t\t#\t\t\t\t\t\ta, b, c, incre = self.getDisp()\n",
    "\t\t#\t\t\t\tx1 = self.position['x'] + a\n",
    "\t\t#\t\t\t\ty1 = self.position['y'] + b\n",
    "\t\t#\t\t\t\taz=radians(self.azimuth(t, J))\n",
    "\t\t#\t\t\t\tal=radians(self.altitude(t, J))\n",
    "\t\t#\t\t\t\tself.te = self.t_e(t,J, sun=horizon_master(az, elev[int(y1)][int(x1)], x=x1, y=y1), a_s=90.- slope[int(y1)][int(x1)], alpha_s=aspect[int(y1)][int(x1)], veg=veget[int(y1)][int(x1)])\n",
    "\t\t#\t\t\t\tup_tb =  self.update_tb()\n",
    "\t\t#\t\t\t\tdiff_tb_new = abs(self.t_target - up_tb)\n",
    "\t\t#\t\t\t\tif diff_tb_new <= diff_tb: # and c <= temp_c:\n",
    "\t\t#\t\t\t\t\tdiff_tb = diff_tb_new\n",
    "\t\t#\t\t\t\t\ttemp_c = c\n",
    "\t\t#\t\t\t\t\ttemp_x1, temp_y1 = x1, y1\n",
    "\t\t#\t\t\t\telse:\n",
    "\t\t#\t\t\t\t\tpass\t\t\t\n",
    "\t\t#\t\t\t\tif l == (self.decisions - 1):\n",
    "\t\t#\t\t\t\t\tself.position['x'], self.position['y'], c = temp_x1, temp_y1, temp_c\n",
    "\t\t#\t\t\t\t\tself.dist += temp_c\n",
    "\t\t#\t\t\t\t\tself.orientation += incre\n",
    "\t\t#\t\t\t\t\tself.moved = c\t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "45eb5877883b2b57d84becf5be6e27beb4fde762018cb87fa112a453076e36a8"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}