# -*- coding: utf-8 -*-
"""
Created on Mon Mar  5 13:17:47 2018

@author: JHodges
"""

import numpy as np
import matplotlib
matplotlib.rcParams['ps.useafm'] = True
matplotlib.rcParams['pdf.use14corefonts'] = True
#matplotlib.rcParams['text.usetex'] = True
import matplotlib.pyplot as plt
from PIL import Image, ImageDraw
from scipy.ndimage.interpolation import zoom
import pandas as pd
import util_common as uc
import os

class FuelModel(object):
    ''' This class contains a fuel model for wildfire spread.
    
        Fields:
                
        Functions:

    '''
    __slots__ = ['id','code','name',
                 'fuelBedDepth','moistureOfExtinctionDeadFuel','heatOfCombustionDeadFuel','heatOfCombustionLiveFuel',
                 'fuelLoad1Hour','fuelLoad10Hour','fuelLoad100Hour','fuelLoadLiveHerb','fuelLoadLiveWood',
                 'savr1HourFuel','savrLiveHerb','savrLiveWood',
                 'isDynamic','isReserved']
    
    def __init__(self,idNumber,idCode,idName,
                 fuelBedDepth,moistureOfExtinctionDeadFuel,heatOfCombustionDeadFuel,heatOfCombustionLiveFuel,
                 fuelLoad1Hour,fuelLoad10Hour,fuelLoad100Hour,fuelLoadLiveHerb,fuelLoadLiveWood,
                 savr1HourFuel,savrLiveHerb,savrLiveWood,
                 isDynamic, isReserved):
        self.id = idNumber
        self.code = idCode
        self.name = idName
        self.fuelBedDepth = fuelBedDepth
        self.moistureOfExtinctionDeadFuel = moistureOfExtinctionDeadFuel
        self.heatOfCombustionDeadFuel = heatOfCombustionDeadFuel
        self.heatOfCombustionLiveFuel = heatOfCombustionLiveFuel
        self.fuelLoad1Hour = fuelLoad1Hour
        self.fuelLoad10Hour = fuelLoad10Hour
        self.fuelLoad100Hour = fuelLoad100Hour
        self.fuelLoadLiveHerb = fuelLoadLiveHerb
        self.fuelLoadLiveWood = fuelLoadLiveWood
        self.savr1HourFuel = savr1HourFuel
        self.savrLiveHerb = savrLiveHerb
        self.savrLiveWood = savrLiveWood
        self.isDynamic = isDynamic
        self.isReserved = isReserved
        
    def __str__(self):
        ''' This function prints summary information of the object when a
        string is requested.
        '''
        string = "Fuel Model\n"
        string = string + "\tID:\t\t%s\n"%(str(self.id))
        string = string + "\tCode:\t%s\n"%(str(self.code))
        string = string + "\tName:\t%s\n"%(str(self.name))
        return string
    
    def __repr__(self):
        ''' This function prints summary information of the object when a
        string is requested.
        '''
        return self.__str__()

class lifeStateIntermediate(object):
    __slots__ = ['dead','live']
    
    def __init__(self,fuelModel,moistureDead,moistureLive):
        
        savrDead, savrLive = getSAV(fuelModel)
        deadFraction, liveFraction, deadFractionTotal, liveFractionTotal = getDLFraction(fuelModel,moistureLive)
        loadDead, loadLive = getFuelLoad(fuelModel,moistureLive)
        heatCDead, heatCLive = getHeatOfCombustion(fuelModel)
        heatDead = np.zeros((len(savrDead),))+heatCDead
        heatLive = np.zeros((len(savrLive),))+heatCLive
        heatLive[liveFraction == 0] = 0
        moistureLive.extend([0,0])
        self.dead = self.calculateIntermediates(fuelModel,savrDead,loadDead,deadFraction,heatDead,moistureDead)
        self.live = self.calculateIntermediates(fuelModel,savrLive,loadLive,liveFraction,heatLive,moistureLive)
        
    def calculateIntermediates(self,fuelModel,savr,load,fraction,heat,moisture):
        totalSilicaContent = 0.0555 # Rothermel 1972
        silicaEffective = np.zeros((len(savr),))+0.01 # From behavePlus source, should be 0 if no fuel
        wn = np.zeros((len(savr),))
        weightedHeat = 0.0
        weightedSilica = 0.0
        weightedMoisture = 0.0
        weightedSavr = 0.0
        totalLoadForLifeState = 0.0
        """
        for i in range(0,len(moisture)):
            wn[i] = load[i]*(1.0-totalSilicaContent)
            weightedHeat = weightedHeat + fraction[i] * heat[i]
            weightedSilica = weightedSilica + fraction[i] * silicaEffective
            weightedMoisture = weightedMoisture + fraction[i]*moisture[i]
            weightedSavr = weightedSavr + fraction[i] * savr[i]
            totalLoadForLifeState = totalLoadForLifeState + load[i]
        """
        #print(fraction,moisture)
        
        wn = [x*(1.0-totalSilicaContent) for x in load] #wn[i] = load[i]*(1.0-totalSilicaContent)
        weightedHeat = np.dot(fraction,heat) #weightedHeat = weightedHeat + fraction[i] * heat[i]
        weightedSilica = np.dot(fraction,silicaEffective) #weightedSilica = weightedSilica + fraction[i] * silicaEffective
        weightedMoisture = np.dot(fraction,moisture) #weightedMoisture = weightedMoisture + fraction[i]*moisture[i]
        weightedSavr = np.dot(fraction,savr) #weightedSavr = weightedSavr + fraction[i] * savr[i]
        totalLoadForLifeState = np.sum(load) #totalLoadForLifeState = totalLoadForLifeState + load[i]
        if fuelModel.isDynamic and False:
            weightedFuelLoad = np.sum(wn)   # This gives better agreement with
                                            # behavePlus for dynamic fuel models;
                                            # however, the source code for
                                            # BehavePlus shows this should be 
                                            # weightedFuelLoad=np.dot(fraction,wn)
        else:
            weightedFuelLoad = np.dot(wn,fraction)
        return [weightedHeat,weightedSilica,weightedMoisture,weightedSavr,totalLoadForLifeState,weightedFuelLoad]

def buildFuelModels(allowDynamicModels=True,allowNonBurningModels=False):
    """
        fuelModelNumber, code, name
        fuelBedDepth, moistureOfExtinctionDeadFuel, heatOfCombustionDeadFuel, heatOfCombustionLiveFuel,
        fuelLoad1Hour, fuelLoad10Hour, fuelLoad100Hour, fuelLoadLiveHerb, fuelLoadLiveWood,
        savr1HourFuel, savrLiveHerb, savrLiveWood,
        isDynamic, isReserved
        - WMC 10/2015
    """
    fuelModels = dict()
    
    # Code FMx: Original 13 Fuel Models
    fuelModels["FM1"] = FuelModel(
            1, "FM1", "Short grass [1]",
            1.0, 0.12, 8000, 8000,
            0.034, 0, 0, 0, 0,
            3500, 1500, 1500,
            False, True)
    fuelModels["FM2"] = FuelModel(
            2, "FM2", "Timber grass and understory [2]",
            1.0, 0.15, 8000, 8000,
            0.092, 0.046, 0.023, 0.023,
            0,3000, 1500, 1500,
            False, True)
    fuelModels["FM3"] = FuelModel(
            3, "FM3", "Tall grass [3]",
            2.5, 0.25, 8000, 8000,
            0.138, 0, 0, 0, 0,
            1500, 1500, 1500,
            False, True)
    fuelModels["FM4"] = FuelModel(
            4, "FM4", "Chaparral [4]",
            6.0, 0.2, 8000, 8000,
            0.230, 0.184, 0.092, 0, 0.230,
            2000, 1500, 1500,
            False, True)
    fuelModels["FM5"] = FuelModel(
            5, "FM5", "Brush [5]",
            2.0, 0.20, 8000, 8000,
            0.046, 0.023, 0, 0, 0.092,
            2000, 1500, 1500,
            False, True)
    fuelModels["FM6"] = FuelModel(
            6, "FM6", "Dormant brush, hardwood slash [6]",
            2.5, 0.25, 8000, 8000,
            0.069, 0.115, 0.092, 0, 0,
            1750, 1500, 1500,
            False, True)
    fuelModels["FM7"] = FuelModel(
            7, "FM7", "Southern rough [7]",
            2.5, 0.40, 8000, 8000,
            0.052, 0.086, 0.069, 0, 0.017,
            1750, 1500, 1500,
            False, True)
    fuelModels["FM8"] = FuelModel(
            8, "FM8", "Short needle litter [8]",
            0.2, 0.3, 8000, 8000,
            0.069, 0.046, 0.115, 0, 0,
            2000, 1500, 1500,
            False, True)
    fuelModels["FM9"] = FuelModel(
            9, "FM9", "Long needle or hardwood litter [9]",
            0.2, 0.25, 8000, 8000,
            0.134, 0.019, 0.007, 0, 0,
            2500, 1500, 1500,
            False, True)
    fuelModels["FM10"] = FuelModel(
            10, "FM10", "Timber litter & understory [10]",
            1.0, 0.25, 8000, 8000,
            0.138, 0.092, 0.230, 0, 0.092,
            2000, 1500, 1500,
            False, True)
    fuelModels["FM11"] = FuelModel(
            11, "FM11", "Light logging slash [11]",
            1.0, 0.15, 8000, 8000,
            0.069, 0.207, 0.253, 0, 0,
            1500, 1500, 1500,
            False, True)
    fuelModels["FM12"] = FuelModel(
            12, "FM12", "Medium logging slash [12]",
            2.3, 0.20, 8000, 8000,
            0.184, 0.644, 0.759, 0, 0,
            1500, 1500, 1500,
            False, True)
    fuelModels["FM13"] = FuelModel(
            13, "FM13", "Heavy logging slash [13]",
            3.0, 0.25, 8000, 8000,
            0.322, 1.058, 1.288, 0, 0,
            1500, 1500, 1500,
            False, True)
    
    if not allowDynamicModels:
        return fuelModels
    else:
        pass
    # 14-89 Available for custom models

    if allowNonBurningModels:
        # Code NBx: Non-burnable
        # 90 Available for custom NB model  
        fuelModels["NB1"] = FuelModel(
                91, "NB1", "Urban, developed [91]",
                1.0, 0.10, 8000, 8000,
                0, 0, 0, 0, 0,
                1500, 1500, 1500,
                False, True)
        fuelModels["NB2"] = FuelModel(
                92, "NB2", "Snow, ice [92]",
                1.0, 0.10, 8000, 8000,
                0, 0, 0, 0, 0,
                1500, 1500, 1500,
                False, True)
        fuelModels["NB3"] = FuelModel(
                93, "NB3", "Agricultural [93]",
                1.0, 0.10, 8000, 8000,
                0, 0, 0, 0, 0,
                1500, 1500, 1500,
                False, True)
    
        # Indices 94-95 Reserved for future standard non-burnable models
    
        fuelModels["NB4"] = FuelModel(
                94, "NB4", "Future standard non-burnable [94]",
                1.0, 0.10, 8000, 8000,
                0, 0, 0, 0, 0,
                1500, 1500, 1500,
                False, True)
        fuelModels["NB5"] = FuelModel(
                95, "NB5", "Future standard non-burnable [95]",
                1.0, 0.10, 8000, 8000,
                0, 0, 0, 0, 0,
                1500, 1500, 1500,
                False, True)
    
        # Indices 96-97 Available for custom NB model
    
        fuelModels["NB8"] = FuelModel(
                98, "NB8", "Open water [98]",
                1.0, 0.10, 8000, 8000,
                0, 0, 0, 0, 0,
                1500, 1500, 1500,
                False, True)
        fuelModels["NB9"] = FuelModel(
                99, "NB9", "Bare ground [99]",
                1.0, 0.10, 8000, 8000,
                0, 0, 0, 0, 0,
                1500, 1500, 1500,
                False, True)

    # Code GRx: Grass
    # Index 100 Available for custom GR model
    f = 2000.0 / 43560.0
    fuelModels["GR1"] = FuelModel(
            101, "GR1", "Short, sparse, dry climate grass (D)",
            0.4, 0.15, 8000, 8000,
            0.10*f, 0, 0, 0.30*f, 0,
            2200, 2000, 1500,
            True, True)
    fuelModels["GR2"] = FuelModel(
            102, "GR2", "Low load, dry climate grass (D)",
            1.0, 0.15, 8000, 8000,
            0.10*f, 0, 0, 1.0*f, 0,
            2000, 1800, 1500,
            True, True)
    fuelModels["GR3"] = FuelModel(
            103, "GR3", "Low load, very coarse, humid climate grass (D)",
            2.0, 0.30, 8000, 8000,
            0.10*f, 0.40*f, 0, 1.50*f, 0,
            1500, 1300, 1500,
            True, True)
    fuelModels["GR4"] = FuelModel(
            104, "GR4", "Moderate load, dry climate grass (D)",
            2.0, 0.15, 8000, 8000,
            0.25*f, 0, 0, 1.9*f, 0,
            2000, 1800, 1500,
            True, True)
    fuelModels["GR5"] = FuelModel(
            105, "GR5", "Low load, humid climate grass (D)",
            1.5, 0.40, 8000, 8000,
            0.40*f, 0.0, 0.0, 2.50*f, 0.0,
            1800, 1600, 1500,
            True, True)
    fuelModels["GR6"] = FuelModel(
            106, "GR6", "Moderate load, humid climate grass (D)",
            1.5, 0.40, 9000, 9000,
            0.10*f, 0, 0, 3.4*f, 0,
            2200, 2000, 1500,
            True, True)
    fuelModels["GR7"] = FuelModel(
            107, "GR7", "High load, dry climate grass (D)",
            3.0, 0.15, 8000, 8000,
            1.0*f, 0, 0, 5.4*f, 0,
            2000, 1800, 1500,
            True, True)
    fuelModels["GR8"] = FuelModel(
            108, "GR8", "High load, very coarse, humid climate grass (D)",
            4.0, 0.30, 8000, 8000,
            0.5*f, 1.0*f, 0, 7.3*f, 0,
            1500, 1300, 1500,
            True, True)
    fuelModels["GR9"] = FuelModel(
            109, "GR9", "Very high load, humid climate grass (D)",
            5.0, 0.40, 8000, 8000,
            1.0*f, 1.0*f, 0, 9.0*f, 0,
            1800, 1600, 1500,
            True, True)
    # 110-112 are reserved for future standard grass models
    # 113-119 are available for custom grass models

    # Code GSx: Grass and shrub
    # 120 available for custom grass and shrub model
    fuelModels["GS1"] = FuelModel(
            121, "GS1", "Low load, dry climate grass-shrub (D)",
            0.9, 0.15, 8000, 8000,
            0.2*f, 0, 0, 0.5*f, 0.65*f,
            2000, 1800, 1800,
            True, True)
    fuelModels["GS2"] = FuelModel(
            122, "GS2", "Moderate load, dry climate grass-shrub (D)",
            1.5, 0.15, 8000, 8000,
            0.5*f, 0.5*f, 0, 0.6*f, 1.0*f,
            2000, 1800, 1800,
            True, True)
    fuelModels["GS3"] = FuelModel(
            123, "GS3", "Moderate load, humid climate grass-shrub (D)",
            1.8, 0.40, 8000, 8000,
            0.3*f, 0.25*f, 0, 1.45*f, 1.25*f,
            1800, 1600, 1600,
            True, True)
    fuelModels["GS4"] = FuelModel(
            124, "GS4", "High load, humid climate grass-shrub (D)",
            2.1, 0.40, 8000, 8000,
            1.9*f, 0.3*f, 0.1*f, 3.4*f, 7.1*f,
            1800, 1600, 1600,
            True, True)
    # 125-130 reserved for future standard grass and shrub models
    # 131-139 available for custom grass and shrub models

    # Shrub
    # 140 available for custom shrub model
    fuelModels["SH1"] = FuelModel(
            141, "SH1", "Low load, dry climate shrub (D)",
            1.0, 0.15, 8000, 8000,
            0.25*f, 0.25*f, 0, 0.15*f, 1.3*f,
            2000, 1800, 1600,
            True, True)
    fuelModels["SH2"] = FuelModel(
            142, "SH2", "Moderate load, dry climate shrub (S)",
            1.0, 0.15, 8000, 8000,
            1.35*f, 2.4*f, 0.75*f, 0, 3.85*f,
            2000, 1800, 1600,
            True, True)
    fuelModels["SH3"] = FuelModel(
            143, "SH3", "Moderate load, humid climate shrub (S)",
            2.4, 0.40, 8000., 8000.,
            0.45*f, 3.0*f, 0, 0, 6.2*f,
            1600, 1800, 1400,
            True, True)
    fuelModels["SH4"] = FuelModel(
            144, "SH4", "Low load, humid climate timber-shrub (S)",
            3.0, 0.30, 8000, 8000,
            0.85*f, 1.15*f, 0.2*f, 0, 2.55*f,
            2000, 1800, 1600,
            True, True)
    fuelModels["SH5"] = FuelModel(
            145, "SH5", "High load, dry climate shrub (S)",
            6.0, 0.15, 8000, 8000,
            3.6*f, 2.1*f, 0, 0, 2.9*f,
            750, 1800, 1600,
            True, True)
    fuelModels["SH6"] = FuelModel(
            146, "SH6", "Low load, humid climate shrub (S)",
            2.0, 0.30, 8000, 8000,
            2.9*f, 1.45*f, 0, 0, 1.4*f,
            750, 1800, 1600,
            True, True)
    fuelModels["SH7"] = FuelModel(
            147, "SH7", "Very high load, dry climate shrub (S)",
            6.0, 0.15, 8000, 8000,
            3.5*f, 5.3*f, 2.2*f, 0, 3.4*f,
            750, 1800, 1600,
            True, True)
    fuelModels["SH8"] = FuelModel(
            148, "SH8", "High load, humid climate shrub (S)",
            3.0, 0.40, 8000, 8000,
            2.05*f, 3.4*f, 0.85*f, 0, 4.35*f,
            750, 1800, 1600,
            True, True)
    fuelModels["SH9"] = FuelModel(
            149, "SH9", "Very high load, humid climate shrub (D)",
            4.4, 0.40, 8000, 8000,
            4.5*f, 2.45*f, 0, 1.55*f, 7.0*f,
            750, 1800, 1500,
            True, True)
    # 150-152 reserved for future standard shrub models
    # 153-159 available for custom shrub models

    # Timber and understory
    # 160 available for custom timber and understory model
    fuelModels["TU1"] = FuelModel(
            161, "TU1", "Light load, dry climate timber-grass-shrub (D)",
            0.6, 0.20, 8000, 8000,
            0.2*f, 0.9*f, 1.5*f, 0.2*f, 0.9*f,
            2000, 1800, 1600,
            True, True)
    fuelModels["TU2"] = FuelModel(
            162, "TU2", "Moderate load, humid climate timber-shrub (S)",
            1.0, 0.30, 8000, 8000,
            0.95*f, 1.8*f, 1.25*f, 0, 0.2*f,
            2000, 1800, 1600,
            True, True)
    fuelModels["TU3"] = FuelModel(
            163, "TU3", "Moderate load, humid climate timber-grass-shrub (D)",
            1.3, 0.30, 8000, 8000,
            1.1*f, 0.15*f, 0.25*f, 0.65*f, 1.1*f,
            1800, 1600, 1400,
            True, True)
    fuelModels["TU4"] = FuelModel(
            164, "TU4", "Dwarf conifer understory (S)",
            0.5, 0.12, 8000, 8000,
            4.5*f, 0, 0, 0, 2.0*f,
            2300, 1800, 2000,
            True, True)
    fuelModels["TU5"] = FuelModel(
            165, "TU5", "Very high load, dry climate timber-shrub (S)",
            1.0, 0.25, 8000, 8000,
            4.0*f, 4.0*f, 3.0*f, 0, 3.0*f,
            1500, 1800, 750,
            True, True)
    # 166-170 reserved for future standard timber and understory models
    # 171-179 available for custom timber and understory models
    # Timber and litter
    # 180 available for custom timber and litter models
    fuelModels["TL1"] = FuelModel(
            181, "TL1", "Low load, compact conifer litter (S)",
            0.2, 0.30, 8000, 8000,
            1.0*f, 2.2*f, 3.6*f, 0, 0,
            2000, 1800, 1600,
            True, True)
    fuelModels["TL2"] = FuelModel(
            182, "TL2", "Low load broadleaf litter (S)",
            0.2, 0.25, 8000, 8000,
            1.4*f, 2.3*f, 2.2*f, 0, 0,
            2000, 1800, 1600,
            True, True)
    fuelModels["TL3"] = FuelModel(
            183, "TL3", "Moderate load conifer litter (S)",
            0.3, 0.20, 8000, 8000,
            0.5*f, 2.2*f, 2.8*f, 0, 0,
            2000, 1800, 1600,
            True, True)
    fuelModels["TL4"] = FuelModel(
            184, "TL4", "Small downed logs (S)",
            0.4, 0.25, 8000, 8000,
            0.5*f, 1.5*f, 4.2*f, 0, 0,
            2000, 1800, 1600,
            True, True)
    fuelModels["TL5"] = FuelModel(
            185, "TL5", "High load conifer litter (S)",
            0.6, 0.25, 8000, 8000,
            1.15*f, 2.5*f, 4.4*f, 0, 0,
            2000, 1800, 160,
            True, True)
    fuelModels["TL6"] = FuelModel(
            186, "TL6", "High load broadleaf litter (S)",
            0.3, 0.25, 8000, 8000,
            2.4*f, 1.2*f, 1.2*f, 0, 0,
            2000, 1800, 1600,
            True, True)
    fuelModels["TL7"] = FuelModel(
            187, "TL7", "Large downed logs (S)",
            0.4, 0.25, 8000, 8000,
            0.3*f, 1.4*f, 8.1*f, 0, 0,
            2000, 1800, 1600,
            True, True)
    fuelModels["TL8"] = FuelModel(
            188, "TL8", "Long-needle litter (S)",
            0.3, 0.35, 8000, 8000,
            5.8*f, 1.4*f, 1.1*f, 0, 0,
            1800, 1800, 1600,
            True, True)
    fuelModels["TL9"] = FuelModel(
            189, "TL9", "Very high load broadleaf litter (S)",
            0.6, 0.35, 8000, 8000,
            6.65*f, 3.30*f, 4.15*f, 0, 0,
            1800, 1800, 1600,
            True, True)
    # 190-192 reserved for future standard timber and litter models
    # 193-199 available for custom timber and litter models
    # Slash and blowdown
    # 200 available for custom slash and blowdown model
    fuelModels["SB1"] = FuelModel(
            201, "SB1", "Low load activity fuel (S)",
            1.0, 0.25, 8000, 8000,
            1.5*f, 3.0*f, 11.0*f, 0, 0,
            2000, 1800, 1600,
            True, True)
    fuelModels["SB2"] = FuelModel(
            202, "SB2", "Moderate load activity or low load blowdown (S)",
            1.0, 0.25, 8000, 8000,
            4.5*f, 4.25*f, 4.0*f, 0, 0,
            2000, 1800, 1600,
            True, True)
    fuelModels["SB3"] = FuelModel(
            203, "SB3", "High load activity fuel or moderate load blowdown (S)",
            1.2, 0.25, 8000, 8000,
            5.5*f, 2.75*f, 3.0*f, 0, 0,
            2000, 1800, 1600,
            True, True)
    fuelModels["SB4"] = FuelModel(
            204, "SB4", "High load blowdown (S)",
            2.7, 0.25, 8000, 8000,
            5.25*f, 3.5*f, 5.25*f, 0, 0,
            2000, 1800, 1600,
            True, True)
    
    return fuelModels

def buildFuelModelsIdx():
    fuelModels = np.empty((256,),dtype=object)
    fuelModels[1] = 'FM1'
    fuelModels[2] = 'FM2'
    fuelModels[3] = 'FM3'
    fuelModels[4] = 'FM4'
    fuelModels[5] = 'FM5'
    fuelModels[6] = 'FM6'
    fuelModels[7] = 'FM7'
    fuelModels[8] = 'FM8'
    fuelModels[9] = 'FM9'
    fuelModels[10] = 'FM10'
    fuelModels[11] = 'FM11'
    fuelModels[12] = 'FM12'
    fuelModels[13] = 'FM13'
    # 14-89 Available for custom models
    fuelModels[91] = 'NB1'
    fuelModels[92] = 'NB2'
    fuelModels[93] = 'NB3'
    # Indices 94-95 Reserved for future standard non-burnable models
    # Indices 96-97 Available for custom NB model
    fuelModels[98] = 'NB8'
    fuelModels[99] = 'NB9'
    # Index 100 Available for custom GR model
    fuelModels[101] = 'GR1'
    fuelModels[102] = 'GR2'
    fuelModels[103] = 'GR3'
    fuelModels[104] = 'GR4'
    fuelModels[105] = 'GR5'
    fuelModels[106] = 'GR6'
    fuelModels[107] = 'GR7'
    fuelModels[108] = 'GR8'
    fuelModels[109] = 'GR9'
    # 110-112 are reserved for future standard grass models
    # 113-119 are available for custom grass models
    # 120 available for custom grass and shrub model
    fuelModels[121] = 'GS1'
    fuelModels[122] = 'GS2'
    fuelModels[123] = 'GS3'
    fuelModels[124] = 'GS4'
    # 125-130 reserved for future standard grass and shrub models
    # 131-139 available for custom grass and shrub models
    # 140 available for custom shrub model
    fuelModels[141] = 'SH1'
    fuelModels[142] = 'SH2'
    fuelModels[143] = 'SH3'
    fuelModels[144] = 'SH4'
    fuelModels[145] = 'SH5'
    fuelModels[146] = 'SH6'
    fuelModels[147] = 'SH7'
    fuelModels[148] = 'SH8'
    fuelModels[149] = 'SH9'
    # 150-152 reserved for future standard shrub models
    # 153-159 available for custom shrub models
    # 160 available for custom timber and understory model
    fuelModels[161] = 'TU1'
    fuelModels[162] = 'TU2'
    fuelModels[163] = 'TU3'
    fuelModels[164] = 'TU4'
    fuelModels[165] = 'TU5'
    # 166-170 reserved for future standard timber and understory models
    # 171-179 available for custom timber and understory models
    # 180 available for custom timber and litter models
    fuelModels[181] = 'TL1'
    fuelModels[182] = 'TL2'
    fuelModels[183] = 'TL3'
    fuelModels[184] = 'TL4'
    fuelModels[185] = 'TL5'
    fuelModels[186] = 'TL6'
    fuelModels[187] = 'TL7'
    fuelModels[188] = 'TL8'
    fuelModels[189] = 'TL9'
    # 190-192 reserved for future standard timber and litter models
    # 193-199 available for custom timber and litter models
    # 200 available for custom slash and blowdown model
    fuelModels[201] = 'SB1'
    fuelModels[202] = 'SB2'
    fuelModels[203] = 'SB3'
    fuelModels[204] = 'SB4'
    return fuelModels

def getFuelModel(fuelModel):
    fuelModels = buildFuelModels(allowDynamicModels=True,allowNonBurningModels=True)
    return fuelModels[fuelModel]

def getMoistureContent(m1h,m10h,m100h,lhm,lwm):
    moistureDead = [m1h,m10h,m100h,m1h]
    moistureLive = [lhm,lwm]
    
    moistureDead = [x/100 for x in moistureDead]
    moistureLive = [x/100 for x in moistureLive]
    
    return moistureDead, moistureLive
    
def getSAV(fuelModel):
    # In behavePlus, there is a conversion to surfaceAreaToVolumeUnits
    savrDead = [fuelModel.savr1HourFuel, 109.0, 30.0, fuelModel.savrLiveHerb]
    savrLive = [fuelModel.savrLiveHerb, fuelModel.savrLiveWood, 0.0, 0.0]
    return savrDead, savrLive

def getFuelLoad(fuelModel,moistureLive):
    loadDead = [fuelModel.fuelLoad1Hour,
                fuelModel.fuelLoad10Hour,
                fuelModel.fuelLoad100Hour,
                0.0]
    loadLive = [fuelModel.fuelLoadLiveHerb,
                fuelModel.fuelLoadLiveWood,
                0.0,
                0.0]
    #print(loadDead)
    #print(loadLive)
    if fuelModel.isDynamic:
        if moistureLive[0] < 0.30:
            loadDead[3] = loadLive[0]
            loadLive[0] = 0.0
        elif moistureLive[0] <= 1.20:
            #print(loadLive[0] * (1.333 - 1.11 * moistureLive[0]))
            loadDead[3] = loadLive[0] * (1.333 - 1.11 * moistureLive[0])
            #loadDead[3] = loadLive[0] * (1.20 - moistureLive[0])/0.9
            loadLive[0] = loadLive[0] - loadDead[3]
        #print(loadLive)
        #print(loadDead)
    #print(loadDead)
    #print(loadLive)
    
    return loadDead, loadLive

def getHeatOfCombustion(fuelModel):
    heatOfCombustionDead = fuelModel.heatOfCombustionDeadFuel
    heatOfCombustionLive = fuelModel.heatOfCombustionLiveFuel
    return heatOfCombustionDead, heatOfCombustionLive

def getDLFraction(fuelModel,moistureLive):
    fuelDensity = 32.0 # Rothermel 1972
    savrDead, savrLive = getSAV(fuelModel)
    loadDead, loadLive = getFuelLoad(fuelModel,moistureLive)
    #print(loadDead)
    #print(savrDead)
    surfaceAreaDead = [x*y/fuelDensity for x,y in zip(loadDead,savrDead)]
    surfaceAreaLive = [x*y/fuelDensity for x,y in zip(loadLive,savrLive)]
    #print(surfaceAreaDead)
        
    totalSurfaceAreaDead = np.sum(surfaceAreaDead)
    totalSurfaceAreaLive = np.sum(surfaceAreaLive)
    
    fractionOfTotalSurfaceAreaDead = totalSurfaceAreaDead/(totalSurfaceAreaDead+totalSurfaceAreaLive)
    fractionOfTotalSurfaceAreaLive = 1.0 - fractionOfTotalSurfaceAreaDead
    
    if totalSurfaceAreaDead > 1.0e-7:
        deadFraction = [x/totalSurfaceAreaDead for x in surfaceAreaDead]
    else:
        deadFraction= [0 for x in surfaceAreaDead]
    
    if totalSurfaceAreaLive > 1.0e-7:
        liveFraction = [x/totalSurfaceAreaLive for x in surfaceAreaLive]
    else:
        liveFraction= [0 for x in surfaceAreaLive]
        
    return deadFraction, liveFraction, fractionOfTotalSurfaceAreaDead, fractionOfTotalSurfaceAreaLive

def getMoistOfExt(fuelModel,moistureDead,moistureLive):
    loadDead, loadLive = getFuelLoad(fuelModel,moistureLive)
    savrDead, savrLive = getSAV(fuelModel)
    
    moistOfExtDead = fuelModel.moistureOfExtinctionDeadFuel
    
    fineDead = 0.0
    fineLive = 0.0
    fineFuelsWeightingFactor = 0.0
    weightedMoistureFineDead = 0.0
    fineDeadMoisture = 0.0
    fineDeadOverFineLive = 0.0
    for i in range(0,len(loadDead)):
        if savrDead[i] > 1.0e-7:
            fineFuelsWeightingFactor = loadDead[i] * np.exp(-138.0/savrDead[i])
        fineDead = fineDead + fineFuelsWeightingFactor
        weightedMoistureFineDead = weightedMoistureFineDead + fineFuelsWeightingFactor * moistureDead[i]
    if fineDead > 1.0e-7:
        fineDeadMoisture = weightedMoistureFineDead / fineDead
    
    for i in range(0,len(loadLive)):
        if savrLive[i] > 1.0e-7:
            fineLive = fineLive + loadLive[i]*np.exp(-500.0/savrLive[i])
    
    if fineLive > 1.0e-7:
        fineDeadOverFineLive = fineDead / fineLive
    
    moistOfExtLive = (2.9 * fineDeadOverFineLive * (1.0 - (fineDeadMoisture) / moistOfExtDead)) - 0.226
    #print("MoEL:",moistOfExtLive)
    
    if moistOfExtLive < moistOfExtDead:
        moistOfExtLive = moistOfExtDead
    
    return moistOfExtDead, moistOfExtLive

def getCharacteristicSAVR(fuelModel,intermediates,moistureLive):
    deadFraction, liveFraction, deadFractionTotal, liveFractionTotal = getDLFraction(fuelModel,moistureLive)
    
    weightedSavrLive = intermediates.live[3]
    weightedSavrDead = intermediates.dead[3]
    
    sigma = deadFractionTotal * weightedSavrDead + liveFractionTotal * weightedSavrLive
    
    return sigma

def getPackingRatios(fuelModel,intermediates,moistureLive):
    fuelDensity = 32.0 # Rothermel 1972
    
    sigma = getCharacteristicSAVR(fuelModel,intermediates,moistureLive)
    
    totalLoadForLifeStateLive = intermediates.live[4]
    totalLoadForLifeStateDead = intermediates.dead[4]
    
    totalLoad = totalLoadForLifeStateLive + totalLoadForLifeStateDead
    
    depth = fuelModel.fuelBedDepth
    bulkDensity = totalLoad / depth
    
    packingRatio = totalLoad / (depth * fuelDensity)
    sigma = round(sigma,0)
    optimumPackingRatio = 3.348 / (sigma**0.8189)
    #packingRatio = round(packingRatio,4)
    relativePackingRatio = packingRatio / optimumPackingRatio
    
    return packingRatio, relativePackingRatio, bulkDensity

def getWeightedFuelLoads(fuelModel,intermediates):
    weightedFuelLoadDead = intermediates.dead[5]
    weightedFuelLoadLive = intermediates.live[5]
    return weightedFuelLoadDead, weightedFuelLoadLive

def getWeightedHeats(fuelModel,intermediates):
    weightedHeatDead = intermediates.dead[0]
    weightedHeatLive = intermediates.live[0]
    return weightedHeatDead, weightedHeatLive
    
def getWeightedSilicas(fuelModel,intermediates):
    weightedSilicaDead = intermediates.dead[1]
    weightedSilicaLive = intermediates.live[1]
    return weightedSilicaDead, weightedSilicaLive

def getHeatSink(fuelModel,moistureDead,moistureLive,bulkDensity):
    savrDead, savrLive = getSAV(fuelModel)
    qigDead = np.zeros((len(savrDead),))
    qigLive = np.zeros((len(savrLive),))
    deadFraction, liveFraction, deadFractionTotal, liveFractionTotal = getDLFraction(fuelModel,moistureLive)
    heatSink = 0
    
    for i in range(0,len(savrDead)):
        if savrDead[i] > 1.0e-7:
            qigDead[i] = 250 + 1116.0 * (moistureDead[i])
            heatSink = heatSink + deadFractionTotal*deadFraction[i]*qigDead[i]*np.exp(-138.0/savrDead[i])
        if savrLive[i] > 1.0e-7:
            qigLive[i] = 250 + 1116.0 * (moistureLive[i])
            heatSink = heatSink + liveFractionTotal*liveFraction[i]*qigLive[i]*np.exp(-138.0/savrLive[i])
    heatSink = heatSink * bulkDensity
    return heatSink

def getHeatFlux(fuelModel,moistureDead,moistureLive,sigma,packingRatio):
    if sigma < 1.0e-7:
        heatFlux = 0.0
    else:
        heatFlux = np.exp((0.792 + 0.681 * sigma**0.5)*(packingRatio + 0.1)) / (192 + 0.2595 * sigma)
    return heatFlux

def getWeightedMoistures(fuelModel,intermediates):
    weightedMoistureDead = intermediates.dead[2]
    weightedMoistureLive = intermediates.live[2]
    return weightedMoistureDead, weightedMoistureLive

def getEtaM(fuelModel,intermediates,MoED,MoEL):
    weightedMoistureDead, weightedMoistureLive = getWeightedMoistures(fuelModel,intermediates)
    
    def calculateEtaM(weightedMoisture,MoE):
        relativeMoisture = 0.0
        if MoE > 0.0:
            relativeMoisture = weightedMoisture / MoE
        if weightedMoisture > MoE or relativeMoisture > 1.0:
            etaM = 0
        else:
            etaM = 1.0 - (2.59 * relativeMoisture) + (5.11 * (relativeMoisture**2))-(3.52*(relativeMoisture**3))
        return etaM
    etaMDead = calculateEtaM(weightedMoistureDead,MoED)
    etaMLive = calculateEtaM(weightedMoistureLive,MoEL)
    
    return etaMDead, etaMLive

def getEtaS(fuelModel,intermediates):
    weightedSilicaDead, weightedSilicaLive = getWeightedSilicas(fuelModel,intermediates)
    
    def calculateEtaS(weightedSilica):
        etaSDen = weightedSilica ** 0.19
        if etaSDen < 1e-6:
            etaS = 0.0
        else:
            etaS = 0.174 / etaSDen
        return min([etaS,1.0])
    
    etaSDead = calculateEtaS(weightedSilicaDead)
    etaSLive = calculateEtaS(weightedSilicaLive)
    
    return etaSDead, etaSLive
    
def getSurfaceFireReactionIntensity(fuelModel,sigma,relativePackingRatio,MoED,MoEL,intermediates):
    aa = 133.0 / (sigma ** 0.7913) # Albini 1976
    gammaMax = (sigma ** 1.5) / (495.0+(0.0594*(sigma**1.5)))
    gamma = gammaMax * (relativePackingRatio**aa) * np.exp(aa * (1.0-relativePackingRatio))
    weightedFuelLoadDead, weightedFuelLoadLive = getWeightedFuelLoads(fuelModel,intermediates)
    weightedHeatDead, weightedHeatLive = getWeightedHeats(fuelModel,intermediates)
    #MoEL = 1.99
    etaMDead, etaMLive = getEtaM(fuelModel,intermediates,MoED,MoEL)
    etaSDead, etaSLive = getEtaS(fuelModel,intermediates)
    
    #print("gamma:",gamma)
    #print("weightedFuelLoadDead/Live:",weightedFuelLoadDead,weightedFuelLoadLive)
    #print("weightedHeatDead/Live:",weightedHeatDead,weightedHeatLive)
    #print("etaMDead/Live",etaMDead,etaMLive)
    #print("etaSDead/Live",etaSDead,etaSLive)
    """
    print("gamma",gamma)
    print("weightedFuelLoadDead",weightedFuelLoadDead)
    print("weightedHeatDead",weightedHeatDead,weightedHeatLive)
    print("etaMDead",etaMDead,etaMLive)
    print("etaSDead",etaSDead,etaSLive)
    """
    reactionIntensityDead = gamma * weightedFuelLoadDead * weightedHeatDead * etaMDead * etaSDead
    reactionIntensityLive = gamma * weightedFuelLoadLive * weightedHeatLive * etaMLive * etaSLive
    
    #reactionIntensityDead = 7505
    
    reactionIntensity = reactionIntensityDead+reactionIntensityLive
    
    #print("Dead Fuel Reaction Intensity: %.0f"%(reactionIntensityDead))
    #print("Live Fuel Reaction Intensity: %.0f"%(reactionIntensityLive))
    #print("Reaction Intensity: %.0f"%(reactionIntensity))
    
    return reactionIntensity, reactionIntensityDead, reactionIntensityLive

def getNoWindNoSlopeSpreadRate(reactionIntensity,heatFlux,heatSink):
    if heatSink < 1.0e-7:
        Rstar = 0.0
    else:
        Rstar = reactionIntensity*heatFlux/heatSink
    #print("HeatSource:",reactionIntensity*heatFlux)
    #print("NoWindNoSlopeSpredRate:",Rstar)
    
    return Rstar

def convertFtMinToChHr(R):
    R = R/1.100
    return R

def calculateMidflameWindSpeed(fuelModel,windSpeed,canopyCover,canopyHeight,crownRatio,
                               windHeightInputMode='TwentyFoot'):
    if windHeightInputMode == 'TenMeter':
        windSpeed = windSpeed/ 1.15
    depth = fuelModel.fuelBedDepth
    canopyCrownFraction = crownRatio * canopyCover / 3.0
    
    if canopyCover < 1.0e-7 or canopyCrownFraction < 0.05 or canopyHeight < 6.0:
        sheltered = False
    else:
        sheltered = True
    
    if sheltered:
        waf = 0.555 / (((canopyCrownFraction * canopyHeight)**0.5)*np.log((20.0+0.36*canopyHeight) / (0.13 * canopyHeight)))
    elif depth > 1.0e-7:
        waf = 1.83 / np.log((20.0+0.36 * depth) / (0.13 * depth))
    else:
        waf = 1.0
    midflameWindSpeed = waf * windSpeed
    return midflameWindSpeed, waf
        
def calculateWindFactor(sigma,relativePackingRatio,mfWindSpeed):
    windC, windB, windE = getWindIntermediates(sigma)
    mfWindSpeed = mfWindSpeed*88 # Convert mph to ft/min
    if mfWindSpeed < 1.0e-7:
        phiW = 0.0
    else:
        phiW = (mfWindSpeed**windB) * windC * (relativePackingRatio**(-windE))
    return phiW

def getWindIntermediates(sigma):
    windC = 7.47 * np.exp(-0.133 * (sigma**0.55))
    windB = 0.02526 * (sigma ** 0.54)
    windE = 0.715 * np.exp(-0.000359*sigma)
    return windC, windB, windE

def calculateSlopeFactor(slope,packingRatio,isAngle=False,isDegree=True):
    if isAngle:
        if isDegree:
            slope = slope/180.0*3.1415926535
        slopex = np.tan(slope)
    else:
        slopex = slope
    phiS = 5.275 * (packingRatio**(-0.3)) * (slopex**2)
    return phiS

def calculateROS(Rstar,phiW,phiS):
    R = Rstar * (1+phiW+phiS)
    return R

def calculateDirectionOfMaxSpread(windDir,aspect,Rstar,phiS,phiW):
    correctedWindDir = windDir-aspect
    windDirRadians = correctedWindDir * 3.1415926535 / 180.0
    slopeRate = Rstar*phiS
    windRate = Rstar*phiW
    
    x = slopeRate + (windRate * np.cos(windDirRadians))
    y = windRate * np.sin(windDirRadians)
    
    rateVector = ((x**2)+(y**2))**0.5
    
    forwardSpreadRate = Rstar + rateVector
    
    azimuth = np.arctan2(y,x) * 180.0 / 3.1415926535
    if azimuth < -1.0e-20:
        azimuth = azimuth + 360
        
    azimuth = azimuth + aspect + 180.0
    
    if azimuth >= 360.0:
        azimuth = azimuth - 360.0
    
    return azimuth, forwardSpreadRate

def calculateWindSpeedLimit(reactionIntensity,phiS):
    windSpeedLimit = 0.9 * reactionIntensity
    if phiS > 0.0:
        if phiS > windSpeedLimit:
            phiS = windSpeedLimit
    return windSpeedLimit

def calculateEffectiveWindSpeed(forwardSpreadRate,Rstar,relativePackingRatio,sigma,windSpeedLimit=9001):
    windC, windB, windE = getWindIntermediates(sigma)
    phiEffectiveWind = forwardSpreadRate/Rstar - 1.0
    effectiveWindSpeed = ((phiEffectiveWind*(relativePackingRatio**windE)) / windC)**(1/windB)
    effectiveWindSpeed = effectiveWindSpeed / 88
    if effectiveWindSpeed > windSpeedLimit/88:
        effectiveWindSpeed = windSpeedLimit/88
    return effectiveWindSpeed
    
def getResidenceTime(sigma):
    if sigma < 1.0e-7:
        residenceTime = 0.0
    else:
        residenceTime = 384.0/sigma
    return residenceTime

def calculateFireBasicDimensions(effectiveWindSpeed,forwardSpreadRate):
    #print("***EFF,",effectiveWindSpeed)#*88/60)
    if effectiveWindSpeed > 1.0e-7:
        fireLengthToWidthRatio = 1.0 + (0.25 * effectiveWindSpeed)#*88/60)
    else:
        fireLengthToWidthRatio = 1.0
    #print("default fl2wr:",fireLengthToWidthRatio)
    #print("default ecc:",(1-(1/fireLengthToWidthRatio)**2)**0.5)
    #fireLengthToWidthRatio = 1.174
    #fireLengthToWidthRatio = 2.25
    #fireLengthToWidthRatio = 1.174 # with effective wind speed 15.7 mi/h
    #fireLengthToWidthRatio = 1.161 # with effective wind speed 8.5 mi/h
    #fireLengthToWidthRatio = 1.145 # with effective wind speed 5.0 mi/h
    x = (fireLengthToWidthRatio**2) - 1.0
    if x > 0.0:
        eccentricity = (x**0.5) / fireLengthToWidthRatio
        #eccentricity = (1-(1/fireLengthToWidthRatio)**2)**0.5
    else:
        eccentricity = 0.0
    #eccentricity = 0.9045
    #print("modded fl2wr:",fireLengthToWidthRatio)
    #print("modded ecc:",eccentricity)
    backingSpreadRate = forwardSpreadRate * (1.0-eccentricity) / (1.0+eccentricity)
    
    ellipticalB = (forwardSpreadRate + backingSpreadRate) / 2.0
    ellipticalC = ellipticalB - backingSpreadRate
    if fireLengthToWidthRatio > 1e-7:
        ellipticalA = ellipticalB / fireLengthToWidthRatio
    else:
        ellipticalA = 0.0
    
    return fireLengthToWidthRatio, eccentricity, backingSpreadRate, ellipticalA, ellipticalB, ellipticalC

def calculateFireFirelineIntensity(forwardSpreadRate,reactionIntensity,residenceTime):
    firelineIntensity = forwardSpreadRate * reactionIntensity * residenceTime / 60.0
    return firelineIntensity

def calculateFlameLength(firelineIntensity):
    flameLength = max([0.0,0.45*(firelineIntensity**0.46)])
    return flameLength

def calculateSpreadRateAtVector(forwardSpreadRate,eccentricity,dirRmax,dirOfInterest):
    if forwardSpreadRate > 0.0:
        beta = abs(dirRmax - dirOfInterest)
        #print("%.1f,%.1f,%.1f"%(dirRmax,dirOfInterest,beta))
        if beta > 180.0:
            beta = (360-beta)
        betaRad = beta * 3.1415926535/180.0
        dirFactor = ((np.cos(betaRad)+1)/2)
        # This is the equation according to the BehavePlus source code:
        rosVector = forwardSpreadRate * (1.0-eccentricity) / (1.0-eccentricity* np.cos(betaRad))
        
        # This is the equaiton I have found to match BehavePlus results:
        rosVector = ((1-abs(betaRad)/3.1415926535)*forwardSpreadRate * dirFactor)
        
        # Combining the two smooths out the peak
        rosVector = ((1-abs(betaRad)/3.1415926535)*forwardSpreadRate * dirFactor + (abs(betaRad)/3.1415926535)*rosVector)
        
        #eccentricity = 0.9
        #rosVector = forwardSpreadRate * (1.0-eccentricity) / (1.0-eccentricity*dirFactor)
        #if beta < 30:
        
        #rosVector = ((1-abs(betaRad)/3.1415926535)*forwardSpreadRate * dirFactor + (betaRad/3.1415926535)*rosVector)
        #print(dirOfInterest,betaRad,rosVector)
    else:
        rosVector = 0.0
    return rosVector

def calculateSpreadRateAtVector2(forwardSpreadRate,backSpreadRate,eccentricity,dirRmax,dirOfInterest):
    if forwardSpreadRate > 0.0:
        beta = abs(dirRmax - dirOfInterest)
        #print("%.1f,%.1f,%.1f"%(dirRmax,dirOfInterest,beta))
        if beta > 180.0:
            beta = (360-beta)
        if abs(beta) > 0.1:
            betaRad = beta * 3.1415926535/180.0
            dirFactor = ((np.cos(betaRad)+1)/2)
            
            # This is the equation according to the BehavePlus source code:
            rosVector = forwardSpreadRate * (1.0-eccentricity) / (1.0-eccentricity* np.cos(betaRad))
            
            # This is the equaiton I have found to match BehavePlus results:
            #rosVector = ((1-abs(betaRad)/3.1415926535)*forwardSpreadRate * dirFactor)
            #rosVector = ((1-abs(betaRad)/3.1415926535)*(forwardSpreadRate-backSpreadRate) * dirFactor)+backSpreadRate
            
            # Combining the two smooths out the peak
            #rosVector = ((1-abs(betaRad)/3.1415926535)*forwardSpreadRate * dirFactor + (abs(betaRad)/3.1415926535)*rosVector)
            
            #eccentricity = 0.9
            #rosVector = forwardSpreadRate * (1.0-eccentricity) / (1.0-eccentricity*dirFactor)
            #if beta < 30:
            
            #rosVector = ((1-abs(betaRad)/3.1415926535)*forwardSpreadRate * dirFactor + (betaRad/3.1415926535)*rosVector)
            #print(dirOfInterest,betaRad,rosVector)
        else:
            rosVector = forwardSpreadRate
        if rosVector < backSpreadRate:
            rosVector = backSpreadRate
    else:
        rosVector = 0.0
    return rosVector
    
    
def scaleRandomValue(mn,mx):
    value = np.random.random()*(mx-mn)+mn
    return value

def getRandomConditions(params,allowDynamicModels=True):
    paramsRand = dict()
    for key in params.keys():
        if params[key][0] == None:
            minValue = params[key][1]
            maxValue = params[key][2]
            if key == 'model':
                fuelModels = list(buildFuelModels(allowDynamicModels=minValue,allowNonBurningModels=maxValue).keys())
                value = fuelModels[np.random.randint(0,len(fuelModels))]
            else:
                value = scaleRandomValue(minValue,maxValue)
        else:
            value = params[key][0]
        paramsRand[key] = value
    return paramsRand

def orderParams(params,toPrint=False):
    model = params['model']
    canopyCover = params['canopyCover']*100
    canopyHeight = params['canopyHeight']
    crownRatio = params['crownRatio']
    m1h = params['m1h']
    m10h = params['m10h']
    m100h = params['m100h']
    lhm = params['lhm']
    lwm = params['lwm']
    windSpeed = params['windSpeed']
    windDir = params['windDir']
    slope = params['slope']*100
    aspect = params['aspect']
    
    orderedParams = [model,canopyCover,canopyHeight,crownRatio,m1h,m10h,m100h,lhm,lwm,windSpeed,windDir,slope,aspect]
    
    if toPrint:
        print("************************************************************")
        print("Starting simulation")
        print("model:\t\t\t%s"%(model))
        print("canopyCover:\t\t%.2f"%(canopyCover))
        print("canopyHeight:\t\t%.2f"%(canopyHeight))
        print("crownRatio:\t\t%.2f"%(crownRatio))
        print("m1h:\t\t\t%.2f"%(m1h))
        print("m10h:\t\t\t%.2f"%(m10h))
        print("m100h:\t\t\t%.2f"%(m100h))
        print("lhm:\t\t\t%.2f"%(lhm))
        print("lwm:\t\t\t%.2f"%(lwm))
        print("windSpeed:\t\t%.2f"%(windSpeed))
        print("windDir:\t\t%.2f"%(windDir))
        print("slope:\t\t\t%.2f"%(slope))
        print("aspect:\t\t\t%.2f"%(aspect))
        
    return orderedParams

def getROSfromParams(params,toPrint=False,maxOnly=False):
    model = params['model']
    canopyCover = params['canopyCover']
    canopyHeight = params['canopyHeight']
    crownRatio = params['crownRatio']
    m1h = params['m1h']
    m10h = params['m10h']
    m100h = params['m100h']
    lhm = params['lhm']
    lwm = params['lwm']
    windSpeed = params['windSpeed']
    windDir = params['windDir']
    slope = params['slope']
    aspect = params['aspect']
    
    orderParams(params,toPrint=toPrint)
    
    directions = np.linspace(0,360,361)

    fuelModel = getFuelModel(model)
    
    moistureDead, moistureLive = getMoistureContent(m1h,m10h,m100h,lhm,lwm)
    loadDead, loadLive = getFuelLoad(fuelModel,moistureLive)
    savrDead, savrLive = getSAV(fuelModel)
    deadFraction, liveFraction, deadFractionTotal, liveFractionTotal = getDLFraction(fuelModel,moistureLive)
    if toPrint:
        print(deadFraction)
        print(liveFraction)
    moistOfExtDead, moistOfExtLive = getMoistOfExt(fuelModel,moistureDead,moistureLive)
    
    heatDead, heatLive = getHeatOfCombustion(fuelModel)
    
    intermediates = lifeStateIntermediate(fuelModel,moistureDead,moistureLive)
    
    sigma = getCharacteristicSAVR(fuelModel,intermediates,moistureLive)
    packingRatio, relativePackingRatio, bulkDensity = getPackingRatios(fuelModel,intermediates,moistureLive)
    heatSink = getHeatSink(fuelModel,moistureDead,moistureLive,bulkDensity)
    heatFlux = getHeatFlux(fuelModel,moistureDead,moistureLive,sigma,packingRatio)
    reactionIntensity, reactionIntensityDead, reactionIntensityLive  = getSurfaceFireReactionIntensity(fuelModel,sigma,relativePackingRatio,moistOfExtDead,moistOfExtLive,intermediates)
    Rstar = getNoWindNoSlopeSpreadRate(reactionIntensity,heatFlux,heatSink)
    
    mfWindSpeed, waf = calculateMidflameWindSpeed(fuelModel,windSpeed,canopyCover,canopyHeight,crownRatio)
    phiW = calculateWindFactor(sigma,relativePackingRatio,mfWindSpeed)
    phiS = calculateSlopeFactor(slope,packingRatio)
    dirRmax, forwardSpreadRate = calculateDirectionOfMaxSpread(windDir,aspect,Rstar,phiS,phiW)
    windSpeedLimit = calculateWindSpeedLimit(reactionIntensity,phiS)
    
    effectiveWindSpeed = calculateEffectiveWindSpeed(forwardSpreadRate,Rstar,relativePackingRatio,sigma,windSpeedLimit=windSpeedLimit)
    residenceTime = getResidenceTime(sigma)
    #effectiveWindSpeed = 3.9
    
    fireLengthToWidthRatio, eccentricity, backingSpreadRate, eA, eB, eC = calculateFireBasicDimensions(effectiveWindSpeed,forwardSpreadRate)
    firelineIntensity = calculateFireFirelineIntensity(forwardSpreadRate,reactionIntensity,residenceTime)
    flameLength = calculateFlameLength(firelineIntensity)
    rosVectors = []
    R = calculateSpreadRateAtVector(forwardSpreadRate,eccentricity,dirRmax,dirRmax)
    R = convertFtMinToChHr(R)

    if toPrint:
        print("************************************************************")
        print("Rate of Spread:\t\t\t\t\t%.1f\tch/h"%(R))
        print("Reaction Intensity:\t\t\t\t%.0f\tBtu/ft2/min"%(reactionIntensity))
        print("Surface Fire Dir of Max Spread (from north):\t%.0f\tdeg"%(dirRmax))
        print("Midflame Wind Speed:\t\t\t\t%.1f\tmi/h"%(mfWindSpeed))
        print("Wind Adjustment Factor:\t\t\t\t%.2f"%(waf))
        print("Effective Wind Speed:\t\t\t\t%.1f\tmi/h"%(effectiveWindSpeed))
        print("Live Fuel Moisture of Extinction:\t\t%.0f"%(moistOfExtLive*100))
        print("Characteristic SA/V:\t\t\t\t%s\tft2/ft3"%(int(sigma)))
        print("Bulk Density:\t\t\t\t\t%.4f\tlbs/ft3"%(bulkDensity))
        print("Packing Ratio:\t\t\t\t\t%.4f"%(packingRatio))
        print("Relative Packing Ratio:\t\t\t\t%.4f"%(relativePackingRatio))
        print("Dead Fuel Reaction Intensity:\t\t\t%.0f\tBtu/ft2/min"%(reactionIntensityDead))
        print("Live Fuel Reaction Intensity:\t\t\t%.0f\tBtu/ft2/min"%(reactionIntensityLive))
        print("Surface Fire Wind Factor:\t\t\t%.1f"%(phiW))
        print("Slope Factor:\t\t\t\t\t%.1f"%(phiS))
        print("Heat Source:\t\t\t\t\t%.0f\tBtu/ft2/min"%(heatFlux*reactionIntensity*(1+phiS+phiW)))
        print("Heat Sink:\t\t\t\t\t%.1f\tBtu/ft3"%(heatSink))
        print("Dead Herbaceous Fuel Load:\t\t\t%.2f\tton/ac"%(loadDead[3]*21.78))
        print("Live Fuel Load Remainder:\t\t\t%.2f\tton/ac"%(loadLive[0]*21.78))
        print("Total Dead Fuel Load:\t\t\t\t%.2f\tton/ac"%(np.sum(loadDead)*21.78))
        print("Total Live Fuel Load:\t\t\t\t%.2f\tton/ac"%(np.sum(loadLive)*21.78))
        print("Dead Fuel Load Portion:\t\t\t\t%.2f"%(np.sum(loadDead)/(np.sum(loadDead)+np.sum(loadLive))*100))
        print("Live Fuel Load Portion:\t\t\t\t%.2f"%(np.sum(loadLive)/(np.sum(loadDead)+np.sum(loadLive))*100))
        print("************************************************************")

    if maxOnly:
        return dirRmax, R
    
    for dirOfInterest in directions:
        rosVector = calculateSpreadRateAtVector2(forwardSpreadRate,backingSpreadRate,eccentricity,dirRmax,dirOfInterest)
        rosVector = convertFtMinToChHr(rosVector)
        rosVectors.append(rosVector)
    
    
    rosVectors = np.array(rosVectors)
    #R = calculateROS(Rstar, phiW, phiS)

    return directions, rosVectors

def cartCoords(thetaRad,rosVectorsKmHr):
    coords = np.zeros((len(rosVectorsKmHr),2))
    x = -1*np.array(rosVectorsKmHr)*np.sin(thetaRad)
    y = -1*np.array(rosVectorsKmHr)*np.cos(thetaRad)
    coords[:,0] = x
    coords[:,1] = y

    return coords

def rothermelOuputToImg(theta,R,resX=50,resY=50):
    coords = cartCoords(theta,R.copy())
       
    (coords[:,0],coords[:,1]) = (coords[:,0]+resX/2,coords[:,1]+resY/2)
    coords = np.array(coords,dtype=np.int32)
    
    coordsTuple = []
    for c in coords:
        coordsTuple.append((c[0],c[1]))
    
    img = Image.new('LA',(resX,resY))
    draw = ImageDraw.Draw(img)
    draw.polygon(coordsTuple,fill='black',outline=None)
    img = np.copy(np.asarray(img)[:,:,1])
    #img[int(resX/2),int(resY/2)] = 125
    
    return img

def rothermelOuputToImgMulti(theta,Rbase,times,resX=50,resY=50):
    img = Image.new('LA',(resX,resY))
    draw = ImageDraw.Draw(img)
    for t in times:
        coords = cartCoords(theta,Rbase.copy()*t)
           
        (coords[:,0],coords[:,1]) = (coords[:,0]+resX/2,coords[:,1]+resY/2)
        coords = np.array(coords,dtype=np.int32)
        
        coordsTuple = []
        for c in coords:
            coordsTuple.append((c[0],c[1]))
        draw.polygon(coordsTuple,fill=(t,t),outline=(t,t))
    img = np.copy(np.asarray(img)[:,:,1])
    #img[int(resX/2),int(resY/2)] = 125
    
    return img
    

def convertChHrToKmHour(R):
    if R is list:
        for r in R:
            r = r*(1.1)*(60.0/5280.0)*(1.60934)
    else:
        R = R*1.1 # Convert to ft/min
        R = R*60.0/5280.0 # Convert to mi/hour
        R = R*1.60934 # Convert to km/hour
    return R

def convertDegToRad(theta):
    if theta is list:
        for r in theta:
            r = r*3.1415926535/180.0
    else:
        theta = theta*3.1415926535/180.0
    return theta

def slopeToElevImg(phi,phiDir,resX=50,resY=50):
    phiDirRad = phiDir*3.1415926535/180.0
    slopeX = phi*np.sin(phiDirRad)
    slopeY = -phi*np.cos(phiDirRad)
    img = np.zeros((2,2))
    #img[img == 0] = np.nan
    img[0,0] = -resX/2*slopeX+resY/2*slopeY
    img[0,-1]= resX/2*slopeX+resY/2*slopeY
    img[-1,0] = -resX/2*slopeX-resY/2*slopeY
    img[-1,-1] = resX/2*slopeX-resY/2*slopeY
    img = zoom(img,resX/2,order=1)
    return img

def visualizeInputImgs(directions,rosVectors,params,resX=50,resY=50,toPlot=True):
    rosVectorsKmHr = convertChHrToKmHour(rosVectors)
    directionsRad = convertDegToRad(directions)
    
    x = -1*np.array(rosVectorsKmHr)*np.sin(directionsRad)
    y = np.array(rosVectorsKmHr)*np.cos(directionsRad)
    
    img6 = rothermelOuputToImg(directionsRad,rosVectorsKmHr*6.0,resX=resX,resY=resY)
    img12 = rothermelOuputToImg(directionsRad,rosVectorsKmHr*12.0,resX=resX,resY=resY)
    img18 = rothermelOuputToImg(directionsRad,rosVectorsKmHr*18.0,resX=resX,resY=resY)
    img24 = rothermelOuputToImg(directionsRad,rosVectorsKmHr*24.0,resX=resX,resY=resY)
    
    elevImg = slopeToElevImg(params['slope'],params['aspect'],resX=resX,resY=resY)
    windDirRad = params['windDir']*3.1415926536/180.0
    windX = np.zeros((resX,resY))+params['windSpeed']*np.sin(windDirRad)
    windY = np.zeros((resX,resY))-params['windSpeed']*np.cos(windDirRad)
    lhmImg = np.zeros((resX,resY))+params['lhm']
    lwmImg = np.zeros((resX,resY))+params['lwm']
    m1hImg = np.zeros((resX,resY))+params['m1h']
    m10hImg = np.zeros((resX,resY))+params['m10h']
    m100hImg = np.zeros((resX,resY))+params['m100h']
    
    canopyCoverImg = np.zeros((resX,resY))+params['canopyCover']
    canopyHeightImg = np.zeros((resX,resY))+params['canopyHeight']
    crownRatioImg = np.zeros((resX,resY))+params['crownRatio']
    
    modelImg = np.zeros((resX,resY))+params['modelInd']
    
    fireImages = [img6,img12,img18,img24]
    modelInputs = [elevImg,windX,windY,lhmImg,lwmImg,m1hImg,m10hImg,m100hImg,canopyCoverImg,canopyHeightImg,crownRatioImg,modelImg]
    
    if toPlot:
        plt.figure(figsize=(12,12))
        plt.suptitle('Fuel Model:%s'%(params['model']))
        plt.subplot(4,4,1)
        plt.imshow(img12,cmap='jet')
        plt.colorbar()
        plt.title('Fire at 12 hours')
        
        plt.subplot(4,4,2)
        plt.imshow(img24,cmap='jet')
        plt.colorbar()
        plt.title('Fire at 24 hours')
        
        plt.subplot(4,4,3)
        plt.imshow(elevImg,cmap='jet')
        plt.colorbar()
        plt.title('Elevation')
        
        plt.subplot(4,4,4)
        plt.imshow(windX,cmap='jet',vmin=-20,vmax=20)
        plt.colorbar()
        plt.title('WindX')
    
        plt.subplot(4,4,5)
        plt.imshow(windY,cmap='jet',vmin=-20,vmax=20)
        plt.colorbar()
        plt.title('WindY')
        
        plt.subplot(4,4,6)
        plt.imshow(lhmImg,cmap='jet',vmin=30,vmax=150)
        plt.colorbar()
        plt.title('Live Herbaceous Moisture')
        
        plt.subplot(4,4,7)
        plt.imshow(lwmImg,cmap='jet',vmin=30,vmax=150)
        plt.colorbar()
        plt.title('Live Woody Moisture')
        
        plt.subplot(4,4,8)
        plt.imshow(m1hImg,cmap='jet',vmin=0,vmax=40)
        plt.colorbar()
        plt.title('1-hour Moisture')
        
        plt.subplot(4,4,9)
        plt.imshow(canopyCoverImg,cmap='jet',vmin=0,vmax=1)
        plt.colorbar()
        plt.title('Canopy Cover')
        
        plt.subplot(4,4,10)
        plt.imshow(canopyHeightImg,cmap='jet',vmin=1,vmax=20)
        plt.colorbar()
        plt.title('Canopy Height')
        
        plt.subplot(4,4,11)
        plt.imshow(crownRatioImg,cmap='jet',vmin=0,vmax=1)
        plt.colorbar()
        plt.title('Crown Ratio')
        
        plt.subplot(4,4,12)
        plt.imshow(modelImg,cmap='jet',vmin=0,vmax=52)
        plt.colorbar()
        plt.title('Model')
        #plt.plot(x,y)
        #plt.plot(0,0,'ok')
        #xRange = x.max()-x.min()
        #yRange = y.max()-y.min()
        #plt.xlim([x.min()-xRange/2,x.max()+xRange/2])
        #plt.ylim([y.min()-yRange/2,y.max()+yRange/2])
        #plt.title('Rate of Spread')

    return fireImages, modelInputs
    


def visualizeInputValues(directions,rosVectors,params,resX=50,resY=50):
    rosVectorsKmHr = convertChHrToKmHour(rosVectors)
    directionsRad = convertDegToRad(directions)
    
    x = -1*np.array(rosVectorsKmHr)*np.sin(directionsRad)
    y = np.array(rosVectorsKmHr)*np.cos(directionsRad)
    
    imgFire = rothermelOuputToImgMulti(directionsRad,rosVectorsKmHr,[48,42,36,30,24,18,12,6],resX=resX,resY=resY)
    imgFire[25,25] = 0
    elevImg = slopeToElevImg(params['slope'],params['aspect'],resX=resX,resY=resY)
    
    windDirRad = params['windDir']*3.1415926536/180.0
    windSpeed = params['windSpeed']
    windX = 1.0*windSpeed*np.sin(windDirRad)
    windY = -1.0*windSpeed*np.cos(windDirRad)
    windYs = [windX,windY]
    windXs = np.arange(len(windYs))
    windNames = ('E+','N+')
    windLimits = [-20,20]
    
    moistYs = [params['m1h'],params['m10h'],params['m100h'],params['lhm']/5,params['lwm']/5]
    moistXs = np.arange(len(moistYs))
    moistNames = ('m1h','m10h','m100h','lhm/5','lwm/5')
    moistLimits = [0,60]
    
    canopyYs = [params['canopyCover'],params['canopyHeight']/20,params['crownRatio']]
    canopyXs = np.arange(len(canopyYs))
    canopyNames = ('Cover (%)','Height (ft/20)','Ratio (%)')
    canopyLimits = [0,1]
    
    modelYs = [params['modelInd'],0]
    modelXs = np.arange(len(modelYs))
    modelNames = (str(params['model']),'')
    modelLimits = [0,52]
    
    plt.figure(figsize=(10,14))
    plt.suptitle('Fuel Model:%s'%(params['model']))
    plt.subplot(3,2,1)
    plt.imshow(imgFire,cmap='gray_r')
    c = plt.colorbar(ticks=[48,36,24,12,0])
    plt.title('Fire spread')
    plt.xlabel('km')
    plt.ylabel('km')
    c.ax.set_label('Hours')
    
    #plt.subplot(3,3,2)
    #plt.imshow(img24,cmap='jet')
    #plt.colorbar()
    #plt.title('Fire at 24 hours')
    
    plt.subplot(3,2,3)
    plt.imshow(elevImg,cmap='jet')
    plt.colorbar()
    plt.title('Elevation Difference [km]')
    
    plt.subplot(3,2,4)
    plt.bar(windXs,windYs,align='center');
    plt.xticks(windXs,windNames);
    plt.ylabel('WindSpeed (mph)');
    plt.ylim(windLimits)

    plt.subplot(3,2,5)
    plt.bar(moistXs,moistYs,align='center');
    plt.xticks(moistXs,moistNames);
    plt.ylabel('Moisture (%)');
    plt.ylim(moistLimits)
    
    plt.subplot(3,2,6)
    plt.bar(canopyXs,canopyYs,align='center');
    plt.xticks(canopyXs,canopyNames);
    plt.ylabel('Canopy (%)');
    plt.ylim(canopyLimits)
    
    plt.subplot(3,2,2)
    plt.bar(modelXs,modelYs,align='center');
    plt.xticks(modelXs,modelNames);
    plt.ylabel('Model Rank');
    plt.ylim(modelLimits)
    plt.xlim([-0.95,0.95])
    #plt.plot(x,y)
    #plt.plot(0,0,'ok')
    #xRange = x.max()-x.min()
    #yRange = y.max()-y.min()
    #plt.xlim([x.min()-xRange/2,x.max()+xRange/2])
    #plt.ylim([y.min()-yRange/2,y.max()+yRange/2])
    #plt.title('Rate of Spread')
    return imgFire

def paramListTodict(paramsRaw):
    params = dict()
    params['model'] = paramsRaw[0]
    params['canopyCover'] = float(paramsRaw[1])/100
    params['canopyHeight'] = float(paramsRaw[2])
    params['crownRatio'] = float(paramsRaw[3])
    params['m1h'] = float(paramsRaw[4])
    params['m10h'] = float(paramsRaw[5])
    params['m100h'] = float(paramsRaw[6])
    params['lhm'] = float(paramsRaw[7])
    params['lwm'] = float(paramsRaw[8])
    params['windSpeed'] = float(paramsRaw[9])
    params['windDir'] = float(paramsRaw[10])
    params['slope'] = float(paramsRaw[11])/100
    params['aspect'] = float(paramsRaw[12])
    return params

def getStandardParams():
    # model,canopyCover/100,Height,Ratio,m1h,m10h,m100h,lhm,l2m,windSpeed,windDir,slope,aspect
    #paramList = ['FM1',0,0,0.5,8,6,4,60,60,10,0,0.5,0]
    paramList = ['FM1',0,0,0.5,8,9,10,60,60,10,0,0.5,0]
    params = paramListTodict(paramList)
    return params

def determineFastestModel(params=None,toPrint=False):
    if params is None:
        params = getStandardParams()
    fuelModels = buildFuelModels(allowDynamicModels=True,allowNonBurningModels=True)
    
    updatedModels = []
    Rs = []
    
    for fuelModel in list(fuelModels.keys()):
        params['model'] = fuelModel
        direction, R = getROSfromParams(params,maxOnly=True)
        updatedModels.append(fuelModel)
        Rs.append(R)
    
    numZero = len(np.where(np.array(Rs) <= 0.01)[0])
    
    inds = np.argsort(Rs)
    updatedModelsSort = np.array(updatedModels)[inds]
    RsSort = np.sort(Rs)
    
    modelIndexDict = dict()
    
    for i in range(0,len(inds)):
        value = max(0,i-numZero+1)
        modelIndexDict[updatedModelsSort[i]] = value
        if toPrint:
            print("Model = %s,\tR = %.2f"%(updatedModelsSort[i],RsSort[i]))
    
    return modelIndexDict

def rearrangeDatas(datas):
    sz = datas[0].shape
    szrs = sz[0]*sz[1]
    datasNew = np.zeros((szrs*len(datas),))
    
    for i in range(0,len(datas)):
        datasNew[i*szrs:(i+1)*szrs] = np.reshape(datas[i],(szrs,))
    return datasNew
            
def getStandardParamsInput():
    paramsInput = dict()
    paramsInput['model']        = [None,True,False] # string
    paramsInput['canopyCover']  = [None,0.0,1.0]    # percent (0-1)
    paramsInput['canopyHeight'] = [None,1.0,20.0]   # ft (1-20)
    paramsInput['crownRatio']   = [None,0.1,1.0]    # fraction (0.1-1)
    paramsInput['m1h']          = [None,1.0,40.0]  # percent (1-60)
    paramsInput['m10h']         = [None,1.0,40.0]     # percent (1-60)
    paramsInput['m100h']        = [None,1.0,40.0]     # percent (1-60)
    paramsInput['lhm']          = [None,30.0,100.0] # percent (30-300)
    paramsInput['lwm']          = [None,30.0,100.0] # percent (30-300)
    paramsInput['windSpeed']    = [None,0.0,30.0]   # mph (0-30)
    paramsInput['windDir']      = [None,0.0,360.0]  # degrees (0-360)
    paramsInput['slope']        = [None,0.0,1.0]    # fraction (0-1)
    paramsInput['aspect']       = [None,0.0,360.0]  # degrees (0-360)
    paramsInput['Mth']          = [None,5,9] # integer
    paramsInput['Day']          = [None,0,31] # integer
    paramsInput['Pcp']          = [None,0.3,10.9] # mm
    paramsInput['mTH']          = [None,400,600] # 24-Hour
    paramsInput['xTH']          = [None,1200,1500] # 24-Hour
    paramsInput['mT']           = [None,2.0,16.6] # degrees C
    paramsInput['xT']           = [None,28.9,37.2] # degrees C
    paramsInput['mH']           = [None,39.2,50.0] # Percent
    paramsInput['xH']           = [None,39.2,50.0] # Percent
    paramsInput['PST']          = [None,0,2400] # Precipitation Start Time
    paramsInput['PET']          = [None,0,2400] # Precipitation End Time
    paramsInput['startTime']    = [None,0,24] # Fire start hour
    
    return paramsInput

def manyFiresInputFigure(modelInputs):

    fig, ax = plt.subplots(figsize=(8,8))
    a = []
    lims = [[-30,30],[-20,20],[-20,20],[30,150],[30,150],[0,30],[0,30],[0,30],[0,1],[0,20],[0,1],[0,53]]
    names = ['Elevation','East Wind','North Wind','Live Herbaceous Moisture','Live Woody Moisture','1-Hour Moisture','10-Hour Moisture','100-Hour Moisture','Canopy Cover','Canopy Height','Crown Ratio','Fuel Model']
    textOffset = [0]
    #modelInputs = [elevImg,windX,windY,lhmImg,lwmImg,m1hImg,m10hImg,m100hImg,canopyCoverImg,canopyHeightImg,crownRatioImg,modelImg]
    for i in range(len(modelInputs)-1,-1,-1):
        img = modelInputs[i].copy()
        img[-1,-1] = lims[i][0]
        img[-1,-2] = lims[i][1]
        oi = OffsetImage(img, zoom = 2.0, cmap='jet')
        box = AnnotationBbox(oi, (-0.5*i,1*i), frameon=True)
        a.append(ax.add_artist(box))
        ax.annotate(names[i],xy=(-0.5*i-1.1,1*i-0.9),xycoords='data',textcoords='data',xytext=(-0.5*i-4-(len(names[i])-10)*0.1,1*i-0.85),arrowprops=dict(facecolor='black',shrink=0.05))
    i = -1
    oi = OffsetImage(imgFire, zoom = 2.0, cmap='jet')
    box = AnnotationBbox(oi, (-0.5*i,1*i), frameon=True)
    ax.annotate('Fire Perimiter',xy=(-0.5*i-1.1,1*i-0.9),xycoords='data',textcoords='data',xytext=(-0.5*i-4,1*i-0.85),arrowprops=dict(facecolor='black',shrink=0.05))
    a.append(ax.add_artist(box))
    plt.xlim(-2,6.15)
    plt.ylim(-1.9,12.2)
    plt.xlim(-9.0,1.4)
    #plt.ylim(-50,50)
    plt.axis('off')
    plt.tight_layout()
    plt.savefig('..%soutputs%sinputsExampleManyFires.png'%(os.sep, os.sep),dpi=300)

def makeFirePerimetersFigure(imgFire):
    import skimage.transform as sktf
    import skimage.filters as skfi
    from mpl_toolkits.axes_grid1 import make_axes_locatable
    oi = skfi.gaussian(imgFire,sigma=1.0,preserve_range=True)
    imgFire = visualizeInputValues(directions,rosVectors,params,resX=250,resY=250)
    imgFire[125:126,125:126] = 0
    imgFire = imgFire[25:175,100:]
    imgFire = imgFire[::-1,:]
    #oi = OffsetImage(imgFire, zoom = 2.0, cmap='jet')
    plt.figure(figsize=(12,12))
    ax = plt.gca()
    fs=32
    im = ax.imshow(imgFire,cmap='hot_r')
    plt.gca().invert_yaxis()
    plt.xlabel('km',fontsize=fs)
    plt.ylabel('km',fontsize=fs)
    plt.tick_params(labelsize=fs)
    plt.xticks([0,20,40,60,80,100,120,140])
    plt.yticks([0,20,40,60,80,100,120,140])
    divider = make_axes_locatable(ax)
    cax = divider.append_axes("right",size="5%", pad=0.05)
    c = plt.colorbar(im,ticks=[48,36,24,12,0],cax=cax)
    #plt.title('Fire spread')
    plt.tick_params(labelsize=fs)
    plt.ylabel('Hours',fontsize=fs)
    #c.ax.set_label(fontsize=fs)
    plt.tight_layout()
    
    plt.savefig('..%soutputs%sexampleFirePerimiter.eps'%(os.sep, os.sep))

if __name__ == "__main__":
    ''' case0: Generate 1 set of random inputs and visualize the results.
        case1: Generate set of 100 random inputs and save inputs for validation
               with behavePlus.
        case2: Re-generate prediction with same random inputs.
        case3: Re-generate single validation output.
        case4: Generate validation plots.
        case5: Generate neural network dataset
    '''
    
    case = 5
    
    paramsInput = dict()
    paramsInput['model']        = [None,True,False] # string
    paramsInput['canopyCover']  = [None,0.0,1.0]    # percent (0-1)
    paramsInput['canopyHeight'] = [None,1.0,20.0]   # ft (1-20)
    paramsInput['crownRatio']   = [None,0.1,1.0]    # fraction (0.1-1)
    paramsInput['m1h']          = [None,1.0,40.0]  # percent (1-60)
    paramsInput['m10h']         = [None,1.0,40.0]     # percent (1-60)
    paramsInput['m100h']        = [None,1.0,40.0]     # percent (1-60)
    paramsInput['lhm']          = [None,30.0,100.0] # percent (30-300)
    paramsInput['lwm']          = [None,30.0,100.0] # percent (30-300)
    paramsInput['windSpeed']    = [None,0.0,30.0]   # mph (0-30)
    paramsInput['windDir']      = [None,0.0,360.0]  # degrees (0-360)
    paramsInput['slope']        = [None,0.0,1.0]    # fraction (0-1)
    paramsInput['aspect']       = [None,0.0,360.0]  # degrees (0-360)
    paramsInput['Tmin']         = [None,2.0,16.6] # degrees C
    paramsInput['Tmax']         = [None,28.9,37.2] # degrees C
    
    resX = 50
    resY = 50
    
    """

    """
    #params['m1h'] = 40.0
    #params['windSpeed'] = 0.0 # mph (0-30)
    #params['windDir'] = -135.0 # degrees (0-360)
    #params['slope'] = 0.0 # fraction (0-1)
    #params['aspect'] = 135
    
    if case == 0:
        params = getRandomConditions(paramsInput,allowDynamicModels=True)
        params['model'] = 'TU2'
        params['canopyCover'] = 0.0607 # percent (0-1)
        params['canopyHeight'] = 17.46 # ft (1-20)
        params['crownRatio'] = 0.99 # fraction (0-1)
        params['m1h'] = 8.4 # percent (1-100)
        params['m10h'] = 6 # percent (1-100)
        params['m100h'] = 4 # percent (1-100)
        params['lhm'] = 82.75 # percent (30-300)
        params['lwm'] = 75.98 # percent (30-300)
        params['windSpeed'] = 12.08 # mph (0-30)
        params['windDir'] = 223.57 # degrees (0-360)
        params['slope'] = 0.9942 # fraction (0-1)
        params['aspect'] = 248.29 # degrees (0-360)
        directions, rosVectors = getROSfromParams(params,toPrint=True)
        visualizeInputImgs(directions,rosVectors,params,resX=resX,resY=resY)
        visualizeInputValues(directions,rosVectors,params,resX=resX,resY=resY)
    elif case == 1:
        allParams = []
        allDirections = []
        allRosVectors = []
        for i in range(0,1000):
            params = getRandomConditions(paramsInput,allowDynamicModels=True)
            directions, rosVectors = getROSfromParams(params)
            allParams.append(orderParams(params))
            allDirections.append(directions)
            allRosVectors.append(rosVectors)
        allParams = np.array(allParams).T
        #pd.DataFrame(allParams[1:,:],columns=allParams[0,:]).astype(float).round(2).to_csv('../rothermelData/validationInputs.csv')
        #pd.DataFrame(allDirections).T.to_csv('../rothermelData/validationDirections.csv')
        #pd.DataFrame(allRosVectors).T.to_csv('../rothermelData/validationRosVectors.csv')
    elif case == 2:
        allParams = pd.read_csv('../rothermelData/validationInputs.csv')
        allDirections = []
        allRosVectors = []
        for i in range(1,allParams.values.shape[1]):
            paramsRaw = allParams.values[:,i]
            params = paramListTodict(paramsRaw)
            directions, rosVectors = getROSfromParams(params)
            allParams.append(orderParams(params))
            allDirections.append(directions)
            allRosVectors.append(rosVectors)
        allParams = np.array(allParams).T
        pd.DataFrame(allDirections).T.to_csv('../rothermelData/validationDirections.csv')
        pd.DataFrame(allRosVectors).T.to_csv('../rothermelData/validationRosVectors.csv')
    elif case == 3:
        numToRepeat = 5
        allParams = pd.read_csv('../rothermelData/validationInputs.csv')
        paramsRaw = allParams.values[:,numToRepeat]
        params = paramListTodict(paramsRaw)
        directions, rosVectors = getROSfromParams(params,toPrint=True)
        
        behaveResults = pd.read_csv('../rothermelData/validationBehaveOutputs.csv')
        behaveDirections = behaveResults.values[:,0]
        behaveRos = behaveResults.values[:,numToRepeat]
        
        rosVectorsResample = np.interp(behaveDirections,directions,rosVectors)
        rmse = np.mean((rosVectorsResample-behaveRos)**2)**0.5
        
        plt.figure(figsize=(4,4))
        plt.plot(directions,rosVectors,label='prediction')
        plt.plot(behaveDirections,behaveRos,label='behavePlus')
        plt.legend()
        ylim = [min([np.min(rosVectors),np.min(behaveRos),0]),max([np.max(rosVectors),np.max(behaveRos),1.0])]
        plt.ylim(ylim)
        plt.title(str(i)+': '+str(np.round(rmse,2)))
        
    elif case == 4:
        behaveResults = pd.read_csv('../rothermelData/validationBehaveOutputs.csv')
        behaveDirections = behaveResults.values[:,0]
        allDirections = pd.read_csv('../rothermelData/validationDirections.csv')
        allRosVectors = pd.read_csv('../rothermelData/validationRosVectors.csv')
        
        rmses = []
        toPlot = True
        
        for i in range(1,51):
            behaveRos = behaveResults.values[:,i]
            directions = allDirections.values[:,i]
            rosVectors = allRosVectors.values[:,i]
            rosVectorsResample = np.interp(behaveDirections,directions,rosVectors)
            rmse = np.mean((rosVectorsResample-behaveRos)**2)**0.5
            rmses.append(rmse)
            if toPlot:
                plt.figure(figsize=(4,4))
                plt.plot(directions,rosVectors,label='prediction')
                plt.plot(behaveDirections,behaveRos,label='behavePlus')
                plt.legend()
                ylim = [min([np.min(rosVectors),np.min(behaveRos),0]),max([np.max(rosVectors),np.max(behaveRos),1.0])]
                plt.ylim(ylim)
                plt.title(str(i)+': '+str(np.round(rmse,2)))
    elif case == 5:
        outdir = '../rothermelData/'
        nsbase = outdir+'data'
        modelIndexDict = determineFastestModel()
        datasIn = []
        datasOut = []
        i = 0
        k = 0
        t1 = uc.tic()
        while i <= 0:
        #for i in range(0,4000):
            params = getRandomConditions(paramsInput,allowDynamicModels=True)
            if i == 0:
                params['aspect'] = 160
                params['model'] = 'FM1'
                params['slope'] = 0.805
                params['m1h'] = 5.26
                params['m10h'] = 6.26
                params['m100h'] = 7.26
                params['lhm'] = 69
                params['lwm'] = 49
                params['canopyCover'] = 0.7
                params['canopyHeight'] = 14
                params['crownRatio'] = 0.2
                params['windDir'] = 34
                params['windSpeed'] = 13.5

            directions, rosVectors = getROSfromParams(params,toPrint=False)
            params['modelInd'] = modelIndexDict[params['model']]
            if True:
                fireImages, modelInputs = visualizeInputImgs(directions,rosVectors,params,resX=resX,resY=resY,toPlot=False)
                for j in range(0,len(fireImages)-1):
                    data = [fireImages[j]]
                    data.extend(modelInputs)
                    datasIn.append(rearrangeDatas(data))
                    datasOut.append(rearrangeDatas([fireImages[j+1]]))
            
            if i % 1000 == 0 and False:
                datasIn = np.squeeze(datasIn)
                datasOut = np.squeeze(datasOut)
                uc.dumpPickle([datasIn,datasOut],outdir+'dataRemakeTest'+str(len(datasOut))+'_'+str(k)+'.pkl')
                datasIn = []
                datasOut = []
                k = k + 1
            i = i + 1
        print(uc.toc(t1))
        #assert False, "Stopped"
        #uc.dumpPickle([datasIn,datasOut],outdir+'dataBehaveMoist'+str(len(datasOut))+'_'+str(k)+'.pkl')
        imgFire = visualizeInputValues(directions,rosVectors,params,resX=resX,resY=resY)
        


        
    elif case == 6:
        import glob
        outdir = '../rothermelData/'
        dataFile = outdir+'dataBehaveMoist3000'
        files = glob.glob(dataFile+'*.pkl')
        ns = outdir+'behaveMoistData'
        allIn = []
        allOut = []
        for i in range(0,len(files)):
            [inData,outData] = uc.readPickle(files[i])
            allIn.extend(inData)
            allOut.extend(outData)
        datas = (inData,outData)
    
    elif case == 7:
        from matplotlib.offsetbox import OffsetImage, AnnotationBbox
        
        networkRaw = np.loadtxt('exampleNetworkRaw0.csv',delimiter=',')
        networkProcessed = np.loadtxt('exampleNetworkProcessed0.csv',delimiter=',')
        
        fs = 16
        
        fig, ax = plt.subplots(figsize=(16,8))
        a = []
        lims = [[-30,30],[-20,20],[-20,20],[30,150],[30,150],[0,30],[0,30],[0,30],[0,1],[0,20],[0,1],[0,53]]
        names = ['Elevation','East Wind','North Wind','Live Herbaceous Moisture','Live Woody Moisture','1-Hour Moisture','10-Hour Moisture','100-Hour Moisture','Canopy Cover','Canopy Height','Crown Ratio','Fuel Model']
        textOffset = [0]
        #modelInputs = [elevImg,windX,windY,lhmImg,lwmImg,m1hImg,m10hImg,m100hImg,canopyCoverImg,canopyHeightImg,crownRatioImg,modelImg]
        for i in range(len(modelInputs)-1,-1,-1):
            img = modelInputs[i].copy()
            img[-1,-1] = lims[i][0]
            img[-1,-2] = lims[i][1]
            oi = OffsetImage(img, zoom = 2.0, cmap='hot_r')
            box = AnnotationBbox(oi, (-0.5*i,1*i), frameon=True)
            a.append(ax.add_artist(box))
            ax.annotate(names[i],xy=(-0.5*i-1.5,1*i-1.6),xycoords='data',textcoords='data',xytext=(-0.5*i-5-(len(names[i])-10)*0.1,1*i-1.525),arrowprops=dict(facecolor='black',shrink=0.05),fontsize=fs)
        i = -1
        fireImages[1][fireImages[1] == 255] = 1
        #fireImages[1][-1,-1] = 2
        norm = matplotlib.colors.Normalize(vmin=0, vmax=2)
        oi = OffsetImage(fireImages[1], zoom = 2.0, cmap='hot_r',norm=norm)
        box = AnnotationBbox(oi, (-0.5*i,1*i), frameon=True)
        ax.annotate('Initial Burn Map',xy=(-0.5*i-1.5,1*i-1.6),xycoords='data',textcoords='data',xytext=(-0.5*i-5,1*i-1.525),arrowprops=dict(facecolor='black',shrink=0.05),fontsize=fs)
        a.append(ax.add_artist(box))
        
        '''
        oi = OffsetImage([[0,0],[0,0]], zoom = 2.0, cmap='hot_r')
        box = AnnotationBbox(oi, (-0.5*i,1*i), frameon=True)
        ax.annotate('Initial Burn Map',xy=(-0.5*i-1.5,1*i-1.6),xycoords='data',textcoords='data',xytext=(-0.5*i-5,1*i-1.525),arrowprops=dict(facecolor='black',shrink=0.05),fontsize=fs)
        a.append(ax.add_artist(box))
        '''
        
        
        
        i = 6# 6
        imgX = 3.50
        ax.annotate('',xy=(imgX,1*i),xycoords='data',textcoords='data',xytext=(imgX-2,1*i),arrowprops=dict(facecolor='black',shrink=0.01,width=50,headwidth=100,headlength=30))
        ax.annotate('Convolutional\nNeural Network',xy=(imgX-0.9,1*i-0.9),xycoords='data',textcoords='data',xytext=(imgX-2.5,1*i-2.85),fontsize=fs)
        imgX = 5.5
        i = 3
        oi = OffsetImage(networkRaw, zoom = 2.0, cmap='hot_r')
        box = AnnotationBbox(oi, (imgX,1*i), frameon=True)
        ax.annotate('Probability of\nFire',xy=(imgX-0.8,1*i-0.9),xycoords='data',textcoords='data',xytext=(imgX-1.1,1*i-2.85),fontsize=fs)
        a.append(ax.add_artist(box))
        imgX = 5.5
        i = 10
        oi = OffsetImage(1-networkRaw, zoom = 2.0, cmap='hot_r')
        box = AnnotationBbox(oi, (imgX,1*i), frameon=True)
        ax.annotate('Probability of\nNot Fire',xy=(imgX-0.8,1*i-0.9),xycoords='data',textcoords='data',xytext=(imgX-1.1,1*i-2.85),fontsize=fs)
        a.append(ax.add_artist(box))
        
        i = 6# 6
        imgX = 9.5
        ax.annotate('',xy=(imgX,1*i),xycoords='data',textcoords='data',xytext=(imgX-2,1*i),arrowprops=dict(facecolor='black',shrink=0.01,width=50,headwidth=100,headlength=30))
        ax.annotate('Post Processing',xy=(imgX-1.1,1*i-1.4),xycoords='data',textcoords='data',xytext=(imgX-2.25,1*i-2.85),fontsize=fs)
        
        imgX = 11.5
        i = 6
        oi = OffsetImage(networkProcessed, zoom = 2.0, cmap='hot_r')
        box = AnnotationBbox(oi, (imgX,1*i), frameon=True)
        ax.annotate('Burn Map\nAfter 6 Hours',xy=(imgX-1.5,1*i-0.9),xycoords='data',textcoords='data',xytext=(imgX-1.15,1*i-2.85),fontsize=fs)
        a.append(ax.add_artist(box))
        
        
        plt.ylim(-2.6,12.5)
        plt.xlim(-10,12.5)
        #plt.ylim(-50,50)
        plt.axis('off')
        plt.tight_layout()
        plt.savefig('..%soutputs%sinputsExampleSingleFire.eps'%(os.sep, os.sep))
        
        makeFirePerimetersFigure(imgFire)
    