TROUBLESHOOT.ipynb

{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [
        {
          "ename": "ValueError",
          "evalue": "zero-size array to reduction operation minimum which has no identity",
          "output_type": "error",
          "traceback": [
            "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
            "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
            "Cell \u001b[0;32mIn [1], line 1158\u001b[0m\n\u001b[1;32m   1154\u001b[0m output\u001b[39m.\u001b[39mwrite(\u001b[39m\"\u001b[39m\u001b[39melev\u001b[39m\u001b[39m\\t\u001b[39;00m\u001b[39mtime\u001b[39m\u001b[39m\\t\u001b[39;00m\u001b[39mX\u001b[39m\u001b[39m\\t\u001b[39;00m\u001b[39mY\u001b[39m\u001b[39m\\t\u001b[39;00m\u001b[39mte\u001b[39m\u001b[39m\\t\u001b[39;00m\u001b[39mtb\u001b[39m\u001b[39m\\t\u001b[39;00m\u001b[39mactive\u001b[39m\u001b[39m\\t\u001b[39;00m\u001b[39mmei\u001b[39m\u001b[39m\\t\u001b[39;00m\u001b[39msmr\u001b[39m\u001b[39m\\t\u001b[39;00m\u001b[39mnet\u001b[39m\u001b[39m\\t\u001b[39;00m\u001b[39mmoved\u001b[39m\u001b[39m\\t\u001b[39;00m\u001b[39mate\u001b[39m\u001b[39m\\n\u001b[39;00m\u001b[39m\"\u001b[39m)\n\u001b[1;32m   1155\u001b[0m output\u001b[39m.\u001b[39mclose()\n\u001b[0;32m-> 1158\u001b[0m do_sim() \u001b[39m#to run it one \u001b[39;00m\n\u001b[1;32m   1159\u001b[0m output\u001b[39m.\u001b[39mclose()\n\u001b[1;32m   1160\u001b[0m output2\u001b[39m.\u001b[39mclose()\n",
            "Cell \u001b[0;32mIn [1], line 1063\u001b[0m, in \u001b[0;36mdo_sim\u001b[0;34m()\u001b[0m\n\u001b[1;32m   1061\u001b[0m \t\t\t\u001b[39mif\u001b[39;00m t \u001b[39m>\u001b[39m\u001b[39m=\u001b[39m \u001b[39m6\u001b[39m \u001b[39mand\u001b[39;00m t \u001b[39m<\u001b[39m \u001b[39m20\u001b[39m:\n\u001b[1;32m   1062\u001b[0m \t\t\t\t\u001b[39mfor\u001b[39;00m lizard \u001b[39min\u001b[39;00m population:\n\u001b[0;32m-> 1063\u001b[0m \t\t\t\t\tsim_act(lizard, t, J)\n\u001b[1;32m   1064\u001b[0m \t\t\t\t\toutput_writer\u001b[39m.\u001b[39mwriterow({\n\u001b[1;32m   1065\u001b[0m                     \u001b[39m'\u001b[39m\u001b[39mtime\u001b[39m\u001b[39m'\u001b[39m: t,\n\u001b[1;32m   1066\u001b[0m                     \u001b[39m'\u001b[39m\u001b[39mJ\u001b[39m\u001b[39m'\u001b[39m: J,\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m   1078\u001b[0m                     \u001b[39m'\u001b[39m\u001b[39menergy_balance\u001b[39m\u001b[39m'\u001b[39m: lizard\u001b[39m.\u001b[39menergy_balance\n\u001b[1;32m   1079\u001b[0m                 })\n\u001b[1;32m   1082\u001b[0m \t\u001b[39mwith\u001b[39;00m \u001b[39mopen\u001b[39m(\u001b[39m'\u001b[39m\u001b[39moverall_averages.csv\u001b[39m\u001b[39m'\u001b[39m, mode\u001b[39m=\u001b[39m\u001b[39m'\u001b[39m\u001b[39mw\u001b[39m\u001b[39m'\u001b[39m, newline\u001b[39m=\u001b[39m\u001b[39m'\u001b[39m\u001b[39m'\u001b[39m) \u001b[39mas\u001b[39;00m output_file: \n",
            "Cell \u001b[0;32mIn [1], line 1032\u001b[0m, in \u001b[0;36mdo_sim.<locals>.sim_act\u001b[0;34m(lizard, t, J)\u001b[0m\n\u001b[1;32m   1028\u001b[0m \t\u001b[39m#lizard.tot_dist_moved += lizard.moved\u001b[39;00m\n\u001b[1;32m   1029\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[1;32m   1030\u001b[0m \t\u001b[39m#print(\"case 2: lizard is active and moves\")\u001b[39;00m\n\u001b[1;32m   1031\u001b[0m \t\u001b[39m#te_list.append(lizard.tb)\u001b[39;00m\n\u001b[0;32m-> 1032\u001b[0m \tlizard\u001b[39m.\u001b[39;49mmove_reg(t,J,elev, slope, aspect, veget)\n\u001b[1;32m   1033\u001b[0m \tlizard\u001b[39m.\u001b[39mtb \u001b[39m=\u001b[39m lizard\u001b[39m.\u001b[39mupdate_tb()\n\u001b[1;32m   1034\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 840\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    837\u001b[0m \t\t\t\t\t\u001b[39mself\u001b[39m\u001b[39m.\u001b[39mmoved \u001b[39m=\u001b[39m c\n\u001b[1;32m    838\u001b[0m \u001b[39m#print(self.te_list)\u001b[39;00m\n\u001b[1;32m    839\u001b[0m \u001b[39m#print(self)\u001b[39;00m\n\u001b[0;32m--> 840\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mte_list_min\u001b[39m=\u001b[39m\u001b[39mmin\u001b[39;49m(\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mte_list)\n\u001b[1;32m    841\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mte_list_max\u001b[39m=\u001b[39m\u001b[39mmax\u001b[39m(\u001b[39mself\u001b[39m\u001b[39m.\u001b[39mte_list)\n\u001b[1;32m    842\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mte_list\u001b[39m=\u001b[39m[]\n",
            "File \u001b[0;32m<__array_function__ internals>:180\u001b[0m, in \u001b[0;36mamin\u001b[0;34m(*args, **kwargs)\u001b[0m\n",
            "File \u001b[0;32m~/Documents/repo_energetics/energetics-model/.venv/lib/python3.8/site-packages/numpy/core/fromnumeric.py:2918\u001b[0m, in \u001b[0;36mamin\u001b[0;34m(a, axis, out, keepdims, initial, where)\u001b[0m\n\u001b[1;32m   2802\u001b[0m \u001b[39m@array_function_dispatch\u001b[39m(_amin_dispatcher)\n\u001b[1;32m   2803\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mamin\u001b[39m(a, axis\u001b[39m=\u001b[39m\u001b[39mNone\u001b[39;00m, out\u001b[39m=\u001b[39m\u001b[39mNone\u001b[39;00m, keepdims\u001b[39m=\u001b[39mnp\u001b[39m.\u001b[39m_NoValue, initial\u001b[39m=\u001b[39mnp\u001b[39m.\u001b[39m_NoValue,\n\u001b[1;32m   2804\u001b[0m          where\u001b[39m=\u001b[39mnp\u001b[39m.\u001b[39m_NoValue):\n\u001b[1;32m   2805\u001b[0m     \u001b[39m\"\"\"\u001b[39;00m\n\u001b[1;32m   2806\u001b[0m \u001b[39m    Return the minimum of an array or minimum along an axis.\u001b[39;00m\n\u001b[1;32m   2807\u001b[0m \n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m   2916\u001b[0m \u001b[39m    6\u001b[39;00m\n\u001b[1;32m   2917\u001b[0m \u001b[39m    \"\"\"\u001b[39;00m\n\u001b[0;32m-> 2918\u001b[0m     \u001b[39mreturn\u001b[39;00m _wrapreduction(a, np\u001b[39m.\u001b[39;49mminimum, \u001b[39m'\u001b[39;49m\u001b[39mmin\u001b[39;49m\u001b[39m'\u001b[39;49m, axis, \u001b[39mNone\u001b[39;49;00m, out,\n\u001b[1;32m   2919\u001b[0m                           keepdims\u001b[39m=\u001b[39;49mkeepdims, initial\u001b[39m=\u001b[39;49minitial, where\u001b[39m=\u001b[39;49mwhere)\n",
            "File \u001b[0;32m~/Documents/repo_energetics/energetics-model/.venv/lib/python3.8/site-packages/numpy/core/fromnumeric.py:86\u001b[0m, in \u001b[0;36m_wrapreduction\u001b[0;34m(obj, ufunc, method, axis, dtype, out, **kwargs)\u001b[0m\n\u001b[1;32m     83\u001b[0m         \u001b[39melse\u001b[39;00m:\n\u001b[1;32m     84\u001b[0m             \u001b[39mreturn\u001b[39;00m reduction(axis\u001b[39m=\u001b[39maxis, out\u001b[39m=\u001b[39mout, \u001b[39m*\u001b[39m\u001b[39m*\u001b[39mpasskwargs)\n\u001b[0;32m---> 86\u001b[0m \u001b[39mreturn\u001b[39;00m ufunc\u001b[39m.\u001b[39;49mreduce(obj, axis, dtype, out, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mpasskwargs)\n",
            "\u001b[0;31mValueError\u001b[0m: zero-size array to reduction operation minimum which has no identity"
          ]
        }
      ],
      "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",
        "import csv\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",
        "\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_elev_1m\", \"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_1m\", \"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_1m\", \"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_1m\", \"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",
        "#### modifying the shape of input drone files\n",
        "elev = elev[elev != -9999]\n",
        "shape = elev.shape[0]\n",
        "\n",
        "if shape > (132 * 132):\n",
        "  shape = (132, 132)\n",
        "  elev = np.reshape(elev[:shape[0]*shape[1]], shape)\n",
        "else:\n",
        "  elev = np.reshape(elev, (int(np.sqrt(shape)), int(np.sqrt(shape))))\n",
        "  \n",
        "#### aspect  \n",
        "aspect = aspect[aspect != -9999]\n",
        "shape1 = aspect.shape[0]\n",
        "\n",
        "if shape1 > (132 * 132):\n",
        "  shape1 = (132, 132)\n",
        "  aspect = np.reshape(aspect[:shape1[0]*shape1[1]], shape1)\n",
        "else:\n",
        "  aspect = np.reshape(aspect, (int(np.sqrt(shape1)), int(np.sqrt(shape1))))\n",
        "\n",
        "#### slope  \n",
        "slope = slope[slope != -9999]\n",
        "shape2 = slope.shape[0]\n",
        "\n",
        "if shape2 > (132 * 132):\n",
        "  shape2 = (132, 132)\n",
        "  slope = np.reshape(slope[:shape2[0]*shape2[1]], shape2)\n",
        "else:\n",
        "  slope = np.reshape(slope, (int(np.sqrt(shape2)), int(np.sqrt(shape2))))\n",
        "  \n",
        "  \n",
        "  #### veget  \n",
        "veget = veget[veget != -9999]\n",
        "shape3 = veget.shape[0]\n",
        "\n",
        "if shape3 > (132 * 132):\n",
        "  shape3 = (132, 132)\n",
        "  veget = np.reshape(veget[:shape3[0]*shape3[1]], shape3)\n",
        "else:\n",
        "  veget = np.reshape(veget, (int(np.sqrt(shape3)), int(np.sqrt(shape3))))\n",
        "  \n",
        "  \n",
        "  \n",
        "##### slope  \n",
        "#slope = slope[slope != -9999]\n",
        "#shape = slope.shape[0]\n",
        "#\n",
        "#if shape > (132 * 132):\n",
        "#  shape = (132, 132)\n",
        "#  slope = np.reshape(elev[:shape[0]*shape[1]], shape)\n",
        "#else:\n",
        "#  slope = np.reshape(slope, (int(np.sqrt(shape)), int(np.sqrt(shape))))\n",
        "#  \n",
        "#  \n",
        "#  #### veget  \n",
        "#veget = veget[veget != -9999]\n",
        "#shape = veget.shape[0]\n",
        "#\n",
        "#if shape > (132 * 132):\n",
        "#  shape = (132, 132)\n",
        "#  veget = np.reshape(elev[:shape[0]*shape[1]], shape)\n",
        "#else:\n",
        "#  veget = np.reshape(veget, (int(np.sqrt(shape)), int(np.sqrt(shape))))\n",
        "#    \n",
        "  \n",
        "  \n",
        "  \n",
        "  \n",
        "  \n",
        "  \n",
        "#### BRING in MICROCLIMATE DATA\n",
        "micro_df = pd.read_csv(\"/Users/laurenneel/Desktop/Grass_RGB_DSM_Oct22/microclim_input_era5/hourly_df_LR_2021.csv\")\n",
        "micro_df['datetime']= pd.to_datetime(micro_df['datetime'])\n",
        "\n",
        "#micro_df_minute_res = micro_df.set_index('datetime')\n",
        "#micro_df_minute_res = micro_df_minute_res.asfreq('T')\n",
        "#micro_df_minute_res = micro_df_minute_res.interpolate(method='time')\n",
        "#micro_df_minute_res['new_hour'] = micro_df_minute_res.index.hour\n",
        "#micro_df_minute_res['ordinal_day'] = micro_df_minute_res.index.dayofyear\n",
        "\n",
        "\n",
        "micro_df_minute_res = micro_df.set_index('datetime')\n",
        "micro_df_minute_res = micro_df_minute_res.asfreq('T')\n",
        "micro_df_minute_res = micro_df_minute_res.interpolate(method='time')\n",
        "micro_df_minute_res['new_hour'] = micro_df_minute_res.index.hour\n",
        "micro_df_minute_res['day_of_year'] = micro_df_minute_res.index.dayofyear       \n",
        "# create an array that goes from 0 to 24 in increments of 0.0166666667 and repeat it 525601 times\n",
        "time_array = np.linspace(0, 24, int(24 / 0.0166666667), endpoint=False)\n",
        "micro_df_minute_res['hour_float'] = np.tile(time_array, len(micro_df_minute_res) // len(time_array) + 1)[:len(micro_df_minute_res)]\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#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",
        "\t#def 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",
        "\t#def 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",
        "\t#def 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",
        "\n",
        "\n",
        "\tclass T_e():\n",
        "\t\tdef __init__(self, emissivity = 0.95, phi = 34., longitude = 109., elev = 0., h = 0.07, d = 0.2, emmiss_animal = 0.95, sky = 0.7, r_soil = 0.3, emmiss_soil = 0.94, 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 = 132., 132. #x_max = ncols, y_max = nrows #### CHANGED FROM SEARS 99. for each 8feb!\n",
        "\t\t\tself.position = {'x':rand() * self.x_max, 'y': rand() * self.y_max}\n",
        "\t\t\n",
        "\t\t\tself.thigh =  39. #sears had 35.#Upper voluntary max\n",
        "\t\t\tself.tlow = 29. #sears had 29.\n",
        "\t\t\tself.topt = 33. #sears had 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.min_te = 20.\n",
        "\t\t\tself.max_te = 20.\n",
        "\t\t\tself.te_list =[]\n",
        "\t\t\tself.te_list_min= 20. #LN added\n",
        "\t\t\tself.te_list_max = 20.\n",
        "\t\t\t#self.tb_range_min = 20. #LN added\n",
        "\t\t\t#self.tb_range_max = 20. #LN added\n",
        "\t\t\t#self.te_max_list = [] #LN added\n",
        "\t\t\tself.orientation = vonmises(self.mu, self.kappa)\n",
        "\t\t\tself.dist = 0.\n",
        "\t\t\tself.ctmax = 40.\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 = 35. #previously 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\treturn self.sky**(self.p_a() / (101.3 * cos(radians(self.zenith(t, J)))))\n",
        "\t\t\t#return self.sky**(self.p_a() / (101.3 * cos(radians(self.Z))))\n",
        "\t\t\n",
        "\t\tdef sol_atm(self, t, J, a_s, alpha_s): #def 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, alpha_s): ##LN changed args from def 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, alpha_s) #and self.sol_atm(t, J, a_s=0, alpha_s=0)\n",
        "\t\n",
        "\t\tdef short_ground(self, t, J, a_s, alpha_s): \n",
        "\t\t\treturn self.h_g(t, J, a_s, alpha_s) * 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\t#def 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, sun, Tg_shade, Tg_sun):\n",
        "\t\t\tif Tg_shade is None:\n",
        "\t\t\t\tTg_shade = micro_df_minute_res.loc[J, \"Tg_shade\"].loc[t]\n",
        "\t\t\tif Tg_sun is None:\n",
        "\t\t\t\tTg_sun= micro_df_minute_res.loc[J, \"Tg_sun\"].loc[t]\n",
        "\t\t\tif sun <= 0.5: # in shade\n",
        "\t\t\t\treturn Tg_shade\n",
        "\t\t\tif sun > 0.5: # in sun\n",
        "\t\t\t\treturn Tg_sun\n",
        "\n",
        "\t\tdef long_atmos(self, Ta_2m):## UPDATE FOR MICROCLIM INPUT DATA...SUN Ta at 2m input!!\n",
        "\t\t\tif Ta_2m is None:\n",
        "\t\t\t\tTa_2m = micro_df_minute_res.loc[J, \"Ta_2m\"].loc[t]\n",
        "\t\t\t#Ta_2m= micro_df.Ta_2m\n",
        "\t\t\treturn 53.1e-14 * (Ta_2m + 273.)**6\n",
        "        \n",
        "\t\t#def 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 np.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, sun, Tg_shade, Tg_sun): #### UPDATE FOR MICROCLIM INPUT DATA... sun/shade GROUND input\n",
        "\t\t\tif Tg_shade is None:\n",
        "\t\t\t\tTg_shade = micro_df_minute_res.loc[J, \"Tg_shade\"].loc[t]\n",
        "\t\t\tif Tg_sun is None:\n",
        "\t\t\t\tTg_sun= micro_df_minute_res.loc[J, \"Tg_sun\"].loc[t]\n",
        "\t\t\tif sun <= 0.5:\n",
        "\t\t\t\treturn self.emmiss_soil * self.SIGMA * (Tg_shade + 273)**4\n",
        "\t\t\tif sun > 0.5:\n",
        "\t\t\t\treturn self.emmiss_soil * self.SIGMA * (Tg_sun + 273)**4\n",
        "\t\t\t#return self.emmiss_soil * self.SIGMA * (self.t_ground(sun, micro_df) + 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\t#def R_abs_simple(self, t, J, sun)\n",
        "\t\t\n",
        "\t\tdef R_abs(self, t, J, sun, a_s, alpha_s, Ta_2m, Tg_shade, Tg_sun): #### UPDATE FOR MICROCLIM INPUT DATA..... needs to call the right ground and air temp\n",
        "\t\t\tif Ta_2m is None:\n",
        "\t\t\t\tTa_2m = micro_df_minute_res.loc[J, \"Ta_2m\"].loc[t]\n",
        "\t\t\tif Tg_shade is None:\n",
        "\t\t\t\tTg_shade = micro_df_minute_res.loc[J, \"Tg_shade\"].loc[t]\n",
        "\t\t\tif Tg_sun is None:\n",
        "\t\t\t\tTg_sun = micro_df_minute_res.loc[J, \"Tg_sun\"].loc[t]\n",
        "\t\t\tabsorbed_radiation = 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, a_s, alpha_s)) + 0.5 * self.abs_animal_long * (self.long_atmos(Ta_2m) + self.long_ground(sun, Tg_shade, Tg_sun))\n",
        "\t\t\treturn absorbed_radiation\n",
        "\t\t\t#return 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, a_s, alpha_s)) + 0.5 * self.abs_animal_long * (self.long_atmos(Ta_2m) + self.long_ground(sun, Tg_shade, Tg_sun))\n",
        "\n",
        "\t\tdef Q_rad(self, Ta_shade, Ta_sun, sun):\n",
        "\t\t\tif Ta_shade is None: \n",
        "\t\t\t\tTa_shade = micro_df_minute_res.loc[J, \"Ta_shade\"].loc[t]\n",
        "\t\t\tif Ta_sun is None:\n",
        "\t\t\t\tTa_sun = micro_df_minute_res.loc[J, \"Ta_sun\"].loc[t]\n",
        "\t\t\tif sun <= 0.5:\n",
        "\t\t\t\treturn self.emmiss_animal * self.SIGMA * (Ta_shade + 273.)**4\n",
        "\t\t\tif sun > 0.5:\n",
        "\t\t\t\treturn self.emmiss_animal * self.SIGMA * (Ta_sun + 273.)**4\n",
        "\n",
        "\t\t\t#return self.emmiss_animal * self.SIGMA * (self.t_air(t,J) + 273.)**4\n",
        "\t\t\t#return 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, sun, Ta_shade, Ta_sun): #### input ground data in sun/shade here!!!\n",
        "\t\t\tif Ta_shade is None:\n",
        "\t\t\t\tTa_shade = micro_df_minute_res.loc[J, \"Ta_shade\"].loc[t]\n",
        "\t\t\tif Ta_sun is None:\n",
        "\t\t\t\tTa_sun = micro_df_minute_res.loc[J, \"Ta_sun\"].loc[t]\n",
        "\t\t\tif sun <= 0.5:\n",
        "\t\t\t\treturn 4.* self.SIGMA * ((Ta_shade + 273.)**3) / 29.3\n",
        "\t\t\tif sun > 0.5:\n",
        "\t\t\t\treturn 4.* self.SIGMA * ((Ta_sun + 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, Tg_shade, Tg_sun, Ta_2m, Ta_shade, Ta_sun, 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 == 1: #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\toperative_temp = self.t_ground(sun, Tg_shade, Tg_sun) + (self.R_abs(t, J, sun, a_s, alpha_s, Ta_2m, Tg_shade, Tg_sun) - self.Q_rad(Ta_shade, Ta_sun, sun)) / (29.3 * (self.g_Ha(wind) + self.g_r(sun, Ta_shade, Ta_sun)))\n",
        "\n",
        "\t\t\treturn operative_temp\n",
        "\t\t\t#return self.t_ground(sun, Tg_shade, Tg_sun) + (self.R_abs(t, J, sun, a_s, alpha_s, Ta_2m, Tg_shade, Tg_sun) - self.Q_rad(Ta_shade, Ta_sun, sun)) / (29.3 * (self.g_Ha(wind) + self.g_r(sun, Ta_shade, Ta_sun)))\n",
        "\t\t\n",
        "\t\t#def get_min_te(self, t, J, Tg_shade, Tg_sun, Ta_2m, Ta_shade, Ta_sun, sun, a_s, alpha_s,veg):\n",
        "\t\t#\t#sun=0\n",
        "\t\t#\t#veg=0\n",
        "\t\t#\tmin_te = Ta_shade + ((self.radiation_abs(julian, hour, s, A_S, tau, Ta, Tg, A_L) - self.E_S*STEFAN_BOLTZMANN*Ta**4) / (cp*(self.gr(Ta, cp)+self.gha(u))))\n",
        "\t\t#\tmin_te = self.t_ground(sun, Tg_shade, Tg_sun) + (self.R_abs(t, J, sun, a_s, alpha_s, Ta_2m, Tg_shade, Tg_sun) - self.Q_rad(Ta_shade, Ta_sun, sun)) / (29.3 * (self.g_Ha(wind) + self.g_r(sun, Ta_shade, Ta_sun)))\n",
        "   #Ta +# ((self.radiation_abs(julian, hour, s, A_S, tau, Ta, Tg, A_L) - self.E_S*STEFAN_BOLTZMANN*Ta**4) / (cp*(self.gr(Ta, cp)+self.gha(u))))\n",
        "\t\t#\treturn min_te\n",
        "\n",
        "\t\tdef get_max_te(self, t, J, Tg_shade, Tg_sun, Ta_2m, Ta_shade, Ta_sun, sun, a_s, alpha_s,veg):\n",
        "\t\t\t#sun=1\n",
        "\t\t\t#veg=0\n",
        "\t\t\tmax_te = max(self.t_e(t, J, Tg_shade, Tg_sun, Ta_2m, Ta_shade, Ta_sun, sun, a_s, alpha_s,veg))\n",
        "\t\t\treturn max_te\n",
        "\n",
        "\t\tdef mass(self):\n",
        "\t\t\treturn ((self.h * 100.) * pi * (0.5 * self.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\tnew_tb = self.tb + (1. - exp(-1. / self.h_t)) * (self.te - self.tb) #heating\n",
        "\t\t\t\treturn new_tb\n",
        "\t\t\t\t#return 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\tnew_tb = self.tb + (1. - exp(-1. / self.c_t)) * (self.te - self.tb) #cooling\n",
        "\t\t\t\treturn new_tb\n",
        "\t\t\t\t#return 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_list=[]\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, Tg_shade, Tg_sun, Ta_2m, Ta_shade, Ta_sun, 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#if self.te > -100.0:\n",
        "\t\t\t\t\t\t\t#print(self)\n",
        "\t\t\t\t\t\t\tself.te_list.append(self.te) #append to te_list every line where self.t_e\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, Tg_shade, Tg_sun, Ta_2m, Ta_shade, Ta_sun, 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\tself.te_list.append(self.t_e(t,J, Tg_shade, Tg_sun, Ta_2m, Ta_shade, Ta_sun, 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)])) #te_list.append(self.t_e(t,J, Tg_shade, Tg_sun, Ta_2m, Ta_shade, Ta_sun, 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, Tg_shade, Tg_sun, Ta_2m, Ta_shade, Ta_sun, 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\tself.te_list.append(self.te)\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, Tg_shade, Tg_sun, Ta_2m, Ta_shade, Ta_sun, 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\tself.te_list.append(self.te)\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, Tg_shade, Tg_sun, Ta_2m, Ta_shade, Ta_sun, 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\tself.te_list.append(self.t_e(t,J, Tg_shade, Tg_sun, Ta_2m, Ta_shade, Ta_sun, 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)])) #te_list.append(self.t_e(t,J, Tg_shade, Tg_sun, Ta_2m, Ta_shade, Ta_sun, 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, Tg_shade, Tg_sun, Ta_2m, Ta_shade, Ta_sun, 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\tprint(1.)\n",
        "\t\t\t\t\t\t\t\tself.te_list.append(self.te)\n",
        "\t\t\t\t\t\t\t\tprint(self.te_list)\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, Tg_shade, Tg_sun, Ta_2m, Ta_shade, Ta_sun, 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\tself.te_list.append(self.te)\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#print(self.te_list)\n",
        "\t\t\t#print(self)\n",
        "\t\t\tself.te_list_min=min(self.te_list)\n",
        "\t\t\tself.te_list_max=max(self.te_list)\n",
        "\t\t\tself.te_list=[]\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",
        "\n",
        "# hmax is the max horizon angle in the direction of theta. If hmax is greater than elevation (where the lizard is), then something is blocking the sun from hitting liz (liz is shaded)\n",
        "# if hmax > elevation, lizard is in SHADE --- see Dozier, Bruno, and Downey, 1981. A faster solution to the horizon problem. Computers & Geosciences.\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",
        "\t#shadows=np.copy(elev[:])\n",
        "\t#temperatures = np.copy(elev[:])\n",
        "\t#activity_temps = np.copy(elev[:])\n",
        "\t#critical_temps = np.copy(elev[:])\n",
        "\t#d_temps = np.copy(elev[:])\n",
        "\t#performance = np.copy(elev[:])\n",
        "\t#act_sum=[]\n",
        "\t#crit_sum=[]\n",
        "\t#temp_temp = []b\n",
        "\t#d_temp = []\n",
        "\t#act_temp =[]\n",
        "\t#crit_temp = []\n",
        "\t#p_temp = []\n",
        "\t#temp_percents, d_percents, p_percents=[],[],[]\n",
        "\n",
        "\t#times = np.arange(6.,20.,1./60.).tolist() \n",
        "#\ttimes = \n",
        "\t#pop_size = 20\n",
        "\tpop_size = 1\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",
        "\t#output = open(name_of_sim+\".csv\", 'a')\n",
        "\t#output2 = open(name_of_sim2+\".csv\", 'a')\n",
        "\n",
        "\n",
        "\t#times = np.arange(6.,20.,1./60.).tolist() \n",
        "\n",
        " \n",
        "#\tdef sim_act(lizard, t, J):\n",
        "#\t\tlizard.zenith(t,J=J) #J=J\n",
        "#\t\tprint(\"zenith angle: \", lizard.Z)\n",
        "#\t\tif lizard.Z > 0.:\n",
        "#\t\t\tTa_sun = micro_df_minute_res.loc[index, \"Ta_sun\"]\n",
        "#\t\t\tTa_shade = micro_df_minute_res.loc[index, \"Ta_shade\"]\n",
        "#\t\t\tTg_sun = micro_df_minute_res.loc[index, \"Tg_sun\"]\n",
        "#\t\t\tTg_shade = micro_df_minute_res.loc[index, \"Tg_shade\"]\n",
        "#\t\t\tTa_2m = micro_df_minute_res.loc[index, \"Ta_2m\"]\n",
        "#\t\t\tprint(\"Before floor y:\", lizard.position['y'])\n",
        "#\t\t\tprint(\"Before floor x:\", lizard.position['x'])\n",
        "#\t\t\tlizard.te = lizard.t_e(t,J, Tg_shade, Tg_sun, Ta_2m, Ta_shade, Ta_sun, 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\tprint(\"Lizard operative temp: \", lizard.te)\n",
        "#\t\t\tprint(lizard.zenith(t,J))\n",
        "#\t\t\tif lizard.tb < lizard.te < lizard.tlow:\n",
        "#\t\t\t\tprint(\"case 1: lizard is not active and has not moved\")\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\telif lizard.tlow <= lizard.tb <=lizard.ctmax:\n",
        "#\t\t\t\tprint(\"case 2: lizard is active and moves\")\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\tprint(\"case 3: lizard is not active and has not moved\")\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\tprint(\"case 4: lizard is not active and not moving\")\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\tprint(str(lizard.position['x']), str(lizard.position['y']))\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",
        "\n",
        "\n",
        "## sim_act funciton only tells you info about the temps the lizard chooses.. nothing about other temps in environment\n",
        "\n",
        "\tdef sim_act(lizard, t, J):\n",
        "\t\tlizard.zenith(t,J=J) #J=J\n",
        "\t\t#print(\"zenith angle: \", lizard.Z)\n",
        "\t\t#print(\"Julian day is: \", J)\n",
        "\t\t#print(\"time is:\", t)\n",
        "\t\tif lizard.Z > 0.:\n",
        "\t\t\tTa_sun = micro_df_minute_res.loc[index, \"Ta_sun\"]\n",
        "\t\t\tTa_shade = micro_df_minute_res.loc[index, \"Ta_shade\"]\n",
        "\t\t\tTg_sun = micro_df_minute_res.loc[index, \"Tg_sun\"]\n",
        "\t\t\tTg_shade = micro_df_minute_res.loc[index, \"Tg_shade\"]\n",
        "\t\t\tTa_2m = micro_df_minute_res.loc[index, \"Ta_2m\"]\n",
        "\t\t\t#print(\"Before floor y:\", lizard.position['y'])\n",
        "\t\t\t#print(\"Before floor x:\", lizard.position['x'])\n",
        "\t\t\t#lizard.min_te=lizard.get_min_te(t,J, Tg_shade, Tg_sun, Ta_2m, Ta_shade, Ta_sun, 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=0.)\n",
        "\t\t\t#lizard.max_te=lizard.get_max_te(t,J, Tg_shade, Tg_sun, Ta_2m, Ta_shade, Ta_sun, 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\tlizard.te = lizard.t_e(t,J, Tg_shade, Tg_sun, Ta_2m, Ta_shade, Ta_sun, 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 operative temp: \", lizard.te)\n",
        "\t\t\t#print(lizard.zenith(t,J))\n",
        "\t\t\t#print(\"Sun value is:\", 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']))\n",
        "\t\t\t#print(\"Elev is:\", elev[math.floor(lizard.position['y'])][math.floor(lizard.position['x'])])\n",
        "\t\t\tif lizard.tb < lizard.te < lizard.tlow:\n",
        "\t\t\t\t#print(\"case 1: lizard is not active and has not moved\")\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\telif lizard.tlow <= lizard.tb <=lizard.ctmax:\n",
        "\t\t\t\t#print(\"case 2: lizard is active and moves\")\n",
        "\t\t\t\t#te_list.append(lizard.tb)\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\t#print(\"case 3: lizard is not active and has not moved\")\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\t#print(\"case 4: lizard is not active and not moving\")\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",
        "\n",
        "     \n",
        "\twith open('output.csv', mode='w', newline='') as output_file: ## ALL DATA\n",
        "\t\tfieldnames = ['time', 'J', 'x', 'y', 'te', 'tb', 'total_activity', 'mei', 'net', 'moved', 'smr', 'energy_balance']\n",
        "\t\toutput_writer = csv.DictWriter(output_file, fieldnames=fieldnames)\n",
        "\t\toutput_writer.writeheader()\n",
        "\t\tfor index, row in micro_df_minute_res.iterrows(): \n",
        "\t\t\tt = row['hour_float']\n",
        "\t\t\t#J = 200.\n",
        "\t\t\tJ = row['day_of_year'] \n",
        "\t\t\tif t >= 6 and t < 20:\n",
        "\t\t\t\tfor lizard in population:\n",
        "\t\t\t\t\tsim_act(lizard, t, J)\n",
        "\t\t\t\t\toutput_writer.writerow({\n",
        "                    'time': t,\n",
        "                    'J': J,\n",
        "                    'x': lizard.position['x'],\n",
        "                    'y': lizard.position['y'],\n",
        "                    'te': lizard.te,\n",
        "\t\t\t\t\t#'min_te': lizard.te_list_min,\n",
        "                    #'max_te': lizard.te_list_max,\n",
        "                    'tb': lizard.tb,\n",
        "                    'total_activity': lizard.total_activity,\n",
        "                    'mei': lizard.mei(),\n",
        "                    'net': lizard.net(),\n",
        "                    'moved': lizard.moved,\n",
        "                    'smr': lizard.smr(),\n",
        "                    'energy_balance': lizard.energy_balance\n",
        "                })\n",
        "      \n",
        "     \n",
        "\twith open('overall_averages.csv', mode='w', newline='') as output_file: \n",
        "\t\tfieldnames = ['te', 'tb', 'total_activity', 'mei', 'net', 'moved', 'smr', 'energy_balance']\n",
        "\t\toutput_writer = csv.DictWriter(output_file, fieldnames=fieldnames)\n",
        "\t\toutput_writer.writeheader()\n",
        "\n",
        "    # calculate overall averages\n",
        "\t\tte_avg = sum([lizard.te for lizard in population]) / len(population)\n",
        "\t\ttb_avg = sum([lizard.tb for lizard in population]) / len(population)\n",
        "\t\ttotal_activity_avg = sum([lizard.total_activity for lizard in population]) / len(population)\n",
        "\t\tmei_avg = sum([lizard.mei() for lizard in population]) / len(population)\n",
        "\t\tnet_avg = sum([lizard.net() for lizard in population]) / len(population)\n",
        "\t\tmoved_avg = sum([lizard.moved for lizard in population]) / len(population)\n",
        "\t\tsmr_avg = sum([lizard.smr() for lizard in population]) / len(population)\n",
        "\t\tenergy_balance_avg = sum([lizard.energy_balance for lizard in population]) / len(population)\n",
        "  \n",
        "\t\toutput_writer.writerow({\n",
        "        'te': te_avg,\n",
        "        'tb': tb_avg,\n",
        "        'total_activity': total_activity_avg,\n",
        "        'mei': mei_avg,\n",
        "        'net': net_avg,\n",
        "        'moved': moved_avg,\n",
        "        'smr': smr_avg,\n",
        "        'energy_balance': energy_balance_avg\n",
        "    })\n",
        "\n",
        "\t\twith open('daily_averages.csv', mode='a', newline='') as output_file:\n",
        "\t\t\tfieldnames = ['day_of_year', 'te_avg', 'tb_avg', 'total_activity_avg', 'mei_avg', 'net_avg', 'moved_avg', 'smr_avg', 'energy_balance_avg']\n",
        "\t\t\toutput_writer = csv.DictWriter(output_file, fieldnames=fieldnames)\n",
        "\t\t\toutput_writer.writeheader()\n",
        "\n",
        "\t    # iterate over days and calculate daily averages\n",
        "\t\t\tfor day in range(1, 367):\n",
        "\t\t\t\tlizards_on_day = [lizard for lizard in population if lizard.day_of_year == day]\n",
        "\t\t\t\t#print(f\"Day {day}: {len(lizards_on_day)} lizards\")\n",
        "\t\t\t\tif lizards_on_day:\n",
        "\t\t\t\t\tte_avg = sum([lizard.te for lizard in lizards_on_day]) / len(lizards_on_day)\n",
        "\t\t\t\t\ttb_avg = sum([lizard.tb for lizard in lizards_on_day]) / len(lizards_on_day)\n",
        "\t\t\t\t\ttotal_activity_avg = sum([lizard.total_activity for lizard in lizards_on_day]) / len(lizards_on_day)\n",
        "\t\t\t\t\tmei_avg = sum([lizard.mei() for lizard in lizards_on_day]) / len(lizards_on_day)\n",
        "\t\t\t\t\tnet_avg = sum([lizard.net() for lizard in lizards_on_day]) / len(lizards_on_day)\n",
        "\t\t\t\t\tmoved_avg = sum([lizard.moved for lizard in lizards_on_day]) / len(lizards_on_day)\n",
        "\t\t\t\t\tsmr_avg = sum([lizard.smr() for lizard in lizards_on_day]) / len(lizards_on_day)\n",
        "\t\t\t\t\tenergy_balance_avg = sum([lizard.energy_balance() for lizard in lizards_on_day]) / len(lizards_on_day)\n",
        "\t\t\t\t\toutput_writer.writerow({'day_of_year': day, 'te_avg': te_avg, 'tb_avg': tb_avg, 'total_activity_avg': total_activity_avg, 'mei_avg': mei_avg, 'net_avg': net_avg, 'moved_avg': moved_avg, 'smr_avg': smr_avg, 'energy_balance_avg': energy_balance_avg})\n",
        "\n",
        "\n",
        "#\tfor index, row in micro_df_minute_res.iterrows(): \n",
        "#\t\tt = row['hour_float']\n",
        "#\t\tJ = row['day_of_year']\n",
        "#\t\tif t >= 6 and t < 20:\n",
        "#\t\t\tprint(J, t)\n",
        "#\t\t\t[sim_act(lizard, t, J) for lizard in population]\n",
        "\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",
        "\n",
        "do_sim() #to run it one \n",
        "output.close()\n",
        "output2.close()"
      ]
    }
  ],
  "metadata": {
    "anaconda-cloud": {},
    "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"
    },
    "vscode": {
      "interpreter": {
        "hash": "45eb5877883b2b57d84becf5be6e27beb4fde762018cb87fa112a453076e36a8"
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 4
}